A coupled domain–boundary type meshless method for phase-field modelling of dendritic solidification with the fluid flow

Abstract

Purpose

This study aims to simulate the dendritic growth in Stokes flow by iteratively coupling a domain and boundary type meshless method.

Design/methodology/approach

A preconditioned phase-field model for dendritic solidification of a pure supercooled melt is solved by the strong-form space-time adaptive approach based on dynamic quadtree domain decomposition. The domain-type space discretisation relies on monomial augmented polyharmonic splines interpolation. The forward Euler scheme is used for time evolution. The boundary-type meshless method solves the Stokes flow around the dendrite based on the collocation of the moving and fixed flow boundaries with the regularised Stokes flow fundamental solution. Both approaches are iteratively coupled at the moving solid–liquid interface. The solution procedure ensures computationally efficient and accurate calculations. The novel approach is numerically implemented for a 2D case.

Findings

The solution procedure reflects the advantages of both meshless methods. Domain one is not sensitive to the dendrite orientation and boundary one reduces the dimensionality of the flow field solution. The procedure results agree well with the reference results obtained by the classical numerical methods. Directions for selecting the appropriate free parameters which yield the highest accuracy and computational efficiency are presented.

Originality/value

A combination of boundary- and domain-type meshless methods is used to simulate dendritic solidification with the influence of fluid flow efficiently.

Citation

Dobravec, T., Mavrič, B., Zahoor, R. and Šarler, B. (2023), "A coupled domain–boundary type meshless method for phase-field modelling of dendritic solidification with the fluid flow", International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 33 No. 8, pp. 2963-2981. https://doi.org/10.1108/HFF-03-2023-0131

1. Introduction

The modelling of dendritic growth is of great importance for predicting the microstructure of many metallic materials (Kurz et al., 2019, 2021). Microstructure evolution is closely linked to material properties (Campbell, 2003); hence, one can use numerical modelling to design and optimise high-quality castings. Different approaches are used for modelling dendritic solidification, for example, the cellular automaton model (Reuther and Rettenmayr, 2014; Dobravec et al., 2017), level-set method (Gibou et al., 2003; Tan and Zabaras, 2006) and phase-field (PF) method (Chen, 2002; Boettinger et al., 2002; Dong et al., 2017; Karma and Tourret, 2016). This study uses the PF method, a powerful approach for solving various free boundary problems in materials science (Provatas and Elder, 2010; Steinbach, 2009). The present study tackles the dendritic solidification of pure melts with Stokes flow around the dendrite. We solve a PF model similar to the well-established PF model by Beckermann et al. (1999), which consists of energy and mass conservation equations and PF and Navier–Stokes liquid momentum equations. Our work slightly differs from the work by Beckermann et al. (1999); we solve Stokes instead of the Navier–Stokes liquid momentum equation. Additionally, we use non-linear preconditioning of the PF (Glasner, 2001; Boukellal et al., 2021) to ensure numerical stability when using larger node spacings.

The main aim of the present study is to develop a novel meshless approach to solve the considered PF model. Meshless methods (Atluri, 2004; Liu, 2009; Liu and Gu, 2005) represent an alternative to the mesh-based finite difference, finite volume and finite element methods. Contrary to mesh-based methods, a pre-defined mesh is not a prerequisite for solving the governing equations when using meshless methods. We differentiate between the domain- and boundary-type meshless methods (Liu and Gu, 2005). In the case of domain-type methods, the whole computational domain is discretised by the computational nodes. In the case of the boundary-type methods, the computational nodes are distributed on the boundary of the computational domain only. Examples of domain-type weak-form meshless methods are the element-free Galerkin method (Belytschko et al., 1994) and the radial point-interpolation method (Liu and Gu, 2001). In the group of the domain-type meshless strong-form methods, also known as the meshless collocation methods, we find, for example, the diffuse approximate method (Sadat and Prax, 1996; Hatić et al., 2019; Talat et al., 2018) and the radial basis function generated finite difference (RBF-FD) method (Flyer et al., 2016; Bayona et al., 2017), also known as the local radial basis function collocation method (Šarler and Vertnik, 2006; Kosec and Šarler, 2011; Vertnik et al., 2019; Mramor et al., 2014; Hanoglu and Šarler, 2018; Mavrič and Šarler, 2015). Examples of boundary-type meshless methods are the local boundary integral equation method (Zhu et al., 1998), the boundary-point interpolation method (Gu and Liu, 2002), the boundary radial point interpolation method (Gu and Liu, 2003), the non-singular method of fundamental solutions (Liu and Šarler, 2018), method of regularised sources (MRS) (Wang et al., 2016), and modified MRS (MRSM) (Rek et al., 2021).

In the present study, we develop a novel approach combining domain- and boundary-type meshless methods. The inspiration for the development of such an approach is twofold. First, previous research (Dobravec et al., 2020, 2022) has demonstrated that using the domain-type meshless RBF-FD method in combination with space-time adaptive approach ensures high accuracy and computational efficiency for solving PF and energy conservation equations. Second, the solution of the Stokes flow around an obstacle using the boundary-type meshless method MRSM (Rek et al., 2021) is computationally much less demanding than the traditional approaches for solving momentum and mass conservation equations (Beckermann et al., 1999; Jeong et al., 2001). In the present numerical model, the domain-type approach solves the PF and energy conservation equations and calculates the position of the solid–liquid interface. The boundary-type approach is set over the fluid domain’s moving solid–liquid interface and exterior boundaries for solving the Stokes flow around the dendrite.

2. Governing equations

We consider the solidification of pure supercooled melt in the 2D computational domain Ω with the boundary Γ. We study a simplified case with constant density ρ, specific heat at constant pressure c_p, and thermal conductivity k. The latent heat of melting and the melting temperature are denoted as L_m and T_m, respectively. We use the dimensionless PF model (Karma and Rappel, 1998), where the spatial and temporal coordinates are measured in units of the PF interface thickness and the PF characteristic attachment time, respectively. The PF interface thickness is defined as:

(1) W0=d01α1λ,

where d₀ is the thermal capillary length while α₁ and λ stand for a constant and the free parameter of the PF model, respectively. The PF characteristic attachment time is given as:

(2) τ0=d02DTα2α12λ3,

where α₂ stands for a constant of the PF model and D_T = k/(ρc_p) for the thermal diffusivity. The PF constants are equal to α₁ = 0.8839 and α₂ = 0.6267 (Karma and Rappel, 1998). The selection of free parameter λ has to yield W₀ much smaller than the diffusion length of solidification to ensure a valid PF model.

The PF model (Beckermann et al., 1999) constrains PF values in the interval −1 ≤ ϕ ≤ 1, where ϕ = 1 and ϕ = −1 denote solid and liquid phases, respectively. We use the preconditioned PF (Glasner, 2001):

(3) ψ=2tanh⁡−1(ϕ),

to increase numerical stability for larger node spacings. The energy conservation equation in terms of dimensionless temperature θ = (T − T_m)/(L_m/c_p) reads as:

(4) ∂θ∂t=D∇2θ+1−ϕ2(1+ϕ2∂ψ∂t−v·∇θ),

where D and v = (v_x, v_y) stand for the dimensionless thermal diffusivity and velocity, respectively. Dimensionless D and v are measured in units of W02/τ0 and W₀/τ₀, respectively. The PF equation reads as (Boukellal et al., 2021):

(5) a2(n)∂ψ∂t=2(ϕ−λ(1−ϕ)θ)+2a(n)∇a(n)·∇ψ−2ϕ∇ψ·a(n)+∇·a(n)+a2(n)(∇2ψ−2ϕ|∇ψ|2),

where a(n) and a(n) represent the anisotropy functions. They depend on the normal to the solid–liquid interface:

(6) n=(nx,ny)=−∇ψ|∇ψ|.

We consider the cubic anisotropy of the surface energy. In this case, anisotropy functions read as:

(7) a(n)=1−3ϵ4+4ϵ4(nx4+ny4),

and

(8) a(n)=16ϵ4|∇ψ|a(n)(nx(nx4+ny4−nx2),ny(nx4+ny4−ny2)),

where ϵ₄ stands for the anisotropy strength of the interface energy.

In the melt (ϕ < 0), we consider incompressible Newtonian Stokes flow. The mass and momentum conservation equations read as

(9) ∇·v=0,

and

(10) −∇p+ν∇2v+f=0,

where p, ν and f stand for the dimensionless rescaled pressure, the kinematic viscosity and the body force, respectively. Rescaled pressure is defined as p = P/ρ, where P stands for the pressure. Dimensionless p, ν and f are measured in units of W02/τ02, W02/τ0 and W0/τ02, respectively.

3. Numerical method

3.1 Solution of phase-field and energy conservation equations

The PF and energy conservation equations are solved by space-time adaptive approach (Dobravec et al., 2022) based on dynamic quadtree domain decomposition. Node distribution with fixed node spacing is generated in each quadtree sub-domain. The constant ratio m_Ω between the characteristic size of the quadtree domain and node spacing ensures space adaptivity, as seen on the left in Figure 1. The free parameters of the space-time adaptive approach are the minimum spacing h, the ratio m_Ω, the maximum number of different node spacings n_h, the maximum number of different time steps n_Δt, the overlapping parameter n_overlap and the type of node distribution (Dobravec et al., 2022). The possible types of node distribution are regular and scattered.

The forward Euler scheme and the RBF-FD method are applied to discretise the PF and energy conservation equations in the computational nodes from a quadtree sub-domain. We calculate the minimum stable time step in the forward Euler scheme as:

(11) Δt=αΔtmin⁡(hmax⁡(|v|),14h2max⁡(D,1/(1−ϵ4))),

where α_Δt stands for the time step stability parameter. We use adaptive time-stepping to increase computational efficiency. The stable time step depends on the node spacing; hence, finer time steps are used in quadtree sub-domains with finer node spacing, as seen in Figure 1.

The core of the RBF-FD method is the RBF interpolation of the field values in local support domains. We use polyharmonic spline (PHS) interpolation, i.e. we apply PHSs as RBFs when constructing the RBF-FD method. A PHS Φ is defined as:

(12) Φ(r)=rn, n=1,3,5, …,

where n is the (odd) PHS degree. The PHSs have to be augmented by N_aug monomials to ensure a well-defined interpolation problem (Fasshauer, 2007) and good convergence properties of the RBF-FD method (Flyer et al., 2016; Bayona et al., 2017). A monomial p is defined as:

(13) p1(r˜)=1, p2(r˜)=x˜, p3(r˜)=y˜, p4(r˜)=x˜2 ,…,

where r˜=(x˜,y˜). The PHS interpolation, augmented by monomials, has gained popularity as a choice for the RBF-FD method in the last years because of good performance (Flyer et al., 2016; Bayona et al., 2017; Dobravec et al., 2020).

In the construction of the RBF-FD method, we have to find a local support domain _lΩ for each computational node _lr from a quadtree sub-domain, as seen in Figure 1. Domain _lΩ is defined as a set of nodes {_lr_i} consisting of a computational node _lr and its N – 1 nearest neighbours. Suppose _lr is closest to r among the computational nodes from the quadtree sub-domain; we approximate arbitrary scalar field η at r as:

(14) η(r)≈∑i=1Nlαi Φ(|r−lri|lh)+∑i=1NauglαN+i pi(r−lrlh),

where _lα_i stands for an interpolation coefficient and _lh for the characteristic size of a local support domain. Applying equation (14) at N nodes from a local support domain yields an underdetermined system of equations. Hence, we add additional relations (Dobravec et al., 2022) to ensure a well-determined system (Iserles, 2000). The interpolation from equation (14) is used for calculating finite-difference-like coefficients _lw_k of any linear differential operator D applied on η at _lr:

(15) Dη(lr)≈∑k=1Nlwk η(lrk).

The details of the RBF-FD method and space-time adaptive approach are given in Dobravec et al. (2022).

3.2 Solution of stokes flow

The Stokes flow around the evolving dendrite is solved by the meshless boundary-type method MRSM (Rek et al., 2021). The MRSM has a basis in the method of fundamental solutions (MFS) (Cheng and Hong, 2020; Liu and Šarler, 2018; Šarler, 2006). In the MFS for 2D Stokes flow, the velocity and pressure are given as a sum of M trial functions for velocity vj* and pressure pj* (Rek and Šarler, 2021):

(16) v(r)=∑j=1Mvj*(r),

and

(17) p(r)=∑j=1Mpj*(r).

A trial function is a linear combination of Stokeslets, i.e. fundamental solutions for the Stokes flow. For example, a trial function for pressure is:

(18) pj*(r)=βjxpx*(r,sj)+βjypy*(r,sj),

where px*(r,s) and py*(r,s) are Stokeslets for pressure. Coefficients βjx and βjy are unknown coefficients determined by collocating the solution of the problem at the M boundary nodes {r_i}. The Stokeslets are singular at the source; therefore, the source points {s_i} must be positioned outside the computational domain. Each collocation point r_i on the boundary has a corresponding source point s_i some distance away from r_i in the outward-normal direction:

(19) si=ri+hifsnΓ,

where h_i stands for the spacing between boundary nodes and n_Γ for the outward-facing normal to the boundary of the computational domain at r_i while f_s is a free parameter controlling the distance between r_i and s_i. The MRSM uses regularised Stokeslets (Wen et al., 2017) for trial functions vj*(r) and pj*(r) from equations (16) and (17). The so-called blobs, i.e. bell-shaped functions with shape parameter ϵ, are used in the derivation of the regularised Stokeslets, as seen in Figure 2. A blob is reduced to delta function for ϵ → 0. We set the shape parameter as ϵ = f_ϵh_i, where f_ϵ stands for the free parameter controlling the shape of a blob. The MRSM with f_s = 0 is known as MRS (Wen et al., 2017). The MRS is suitable for solving problems with Dirichlet boundary conditions. The reader can find the details of the MRS and MRSM in Rek et al. (2021) and Wen et al. (2017).

3.3 Coupling domain- and boundary-type methods

The solution procedure consists of initialisation and iteration parts. In the initialisation part, we set the initial conditions for ψ and θ in the computational domain. The iteration part consists of two coupling steps. First, the MRSM calculates the Stokes flow in the computational domain using the nodes at the boundary of the computational domain and the nodes at the solid–liquid interface. Second, the RBF-FD-based adaptive approach solves the PF and energy conservation equations using the Stokes velocity.

The nodes on the solid–liquid interface ψ = 0 are calculated with the following algorithm. In each quadtree sub-domain with the minimum spacing h, a regular node distribution with spacing h is created, and the values of ψ are interpolated to the regular nodes. For each regular node r_i, we check whether the sign of ψ changes when we move one node to the east or to the north. If the change of sign is detected, the following position becomes a boundary node on the solid–liquid interface in the MRSM:

(20) r=ri−(ψ∇ψ|∇ψ|)|ri.

The spacing between the boundary nodes on the solid–liquid interface in the MRSM is approximately equal to the minimum spacing h in the RBF-FD method. Both methods require fine enough spacing h to properly describe the features of the solid–liquid interface. In subsection 4.3, we investigate the influence of the minimum spacing h on the accuracy in the case of diffusion-controlled growth and choose the optimal spacing h. In subsection 4.4, we use the optimal h to analyse the MRSM in the case of convection-diffusion-controlled growth.

3.4 Selection of free numerical parameters

Previous research (Dobravec et al., 2022, 2020) analyses the influence of the many free numerical parameters of the forward Euler scheme, the RBF-FD method and the space-time adaptive approach on the accuracy and computational efficiency in solving PF and energy conservation equations. However, the preconditioned PF model was not considered previously. Hence, we thoroughly repeat the assessment of the RBF-FD method for the case of preconditioning.

As mentioned in subsection 3.1, the space-time adaptive approach has the following free parameters: h, n_h, n_Δt, m_Ω, n_overlap and the type of node distribution. We set n_overlap = 1 and n_Δt = 2; such configuration yields good accuracy and computational efficiency (Dobravec et al., 2022). We test the minimum spacings in quadtree sub-domains h = 0.4, h = 0.8 and h = 1.2. The following sets of free parameters are used (n_h = 6, m_Ω = 9), (n_h = 5, m_Ω = 9) and (n_h = 4, m_Ω = 12) for h = 0.4, h = 0.8 and h = 1.2, respectively. We test the performance using regular and scattered node distributions. The forward Euler scheme has a single free parameter α_Δt. Value α_Δt = 0.3 yields sufficiently small time steps, i.e. further reduction of α_Δt does not increase the method’s accuracy. The RBF-FD method has the following free parameters: n, N and N_aug. We use fifth-degree polyharmonic splines (n = 5) and second-order augmentation with monomials (N_aug = 6). As Dobravec et al. (2020, 2022), we test the performance for N = 9, N = 13 and N = 21 nodes in local support domains.

A numerical model can use larger node spacings when using preconditioning compared to the non-preconditioned PF model. However, when we use space adaptivity, the preconditioning yields stability issues far from the solid–liquid interface, where large node spacings cannot resolve the model. Solving a non-preconditioned PF model by space adaptive algorithm does not experience this problem because the PF is a constant far away from the solid–liquid interface. We tackle this problem by applying the following restriction (Gong et al., 2018):

(21) ψ={ψcψ>ψc−ψcψ<−ψcψ−ψc≤ψ≤ψc,

where ψ_c is the numerical cut-off parameter. Numerical experiments show that ψ_c = 12 yields stable and accurate results. The adaptive algorithm refines the areas in the computational domain where |ψ| < 11 and de-refines the areas where |ψ| > 11.5. The algorithm ensures the minimum node spacing h in quadtree sub-domains where the solid–liquid interface lies. Care is taken to keep the quadtree balanced in the refinement/de-refinement procedure (Dobravec et al., 2022).

As mentioned in subsection 3.2, the MRSM has two free parameters: f_s and f_ϵ. We test the performance of the method for f_s ∈ [0.01,5.12] and f_ϵ ∈ [0.01,5.12]. The spacing between the solid–liquid boundary nodes, set according to equation (20), is approximately equal to the minimum spacing between the computational nodes h. We set the spacing between the nodes at the boundary of the computational domain as h_Γ = f_bh, where f_b is a free parameter. We test the performance for f_b ∈ [1,64]. To save computational time, we execute the MRSM every f_exe-th iteration of the PF and energy conservation equations, where f_exe is a free parameter. We test the performance for f_exe ∈ [1,128].

3.5 Numerical implementation

The novel numerical approach is implemented in the programming language Fortran 2008 and compiled with the Intel Visual Studio Compiler 19.0. The OpenMP (Chapman et al., 2008) application programming interface accelerates the calculations. The DGSEV routine from the LAPACK library (Anderson et al., 1987) solves the system of linear equations in the MRSM. Programming language Python with the libraries Matplotlib and Numpy is used for the post-processing and graphical presentation of the numerical results.

4. Results

4.1 Problem definition

We solve the test case by Beckermann et al. (1999) to test our newly developed numerical model. The test case considers the growth of dendrite from a supercooled melt in a square computational domain Ω = [−L/2, L/2] × [−L/2, L/2], where L stands for the size of the computational domain. The initial condition for PF is a circular nucleus with the origin r₀ and the radius r₀. We set the initial conditions for the PF and energy conservation as:

(22) ψ(r,t=0)=r0−|r−r0|,  θ(r,t=0)=−Δ,

where Δ stands for the initial supercooling. Zero flux Neumann boundary conditions are applied for ψ and θ:

(23) ∇ψ|Γ·nΓ=0, ∇θ|Γ·nΓ=0.

For the velocity, the test case prescribes the inlet Dirichlet boundary conditions on the north part of Γ, the mixed symmetry boundary conditions on the east and west part of Γ, the outlet Neumann boundary conditions on the south part of Γ and the no-slip Dirichlet boundary condition on the solid–liquid interface:

(24) v|Γnorth=(0,−vin), vx|Γeast,west=0, ∂vy∂x|Γeast,west=0, ∂v∂y|Γsouth=0, v|solid-liquid=0,

where v_in represents the absolute value of the inlet velocity. The defined initial and boundary conditions correspond to the dimensionless variables, defined in Section 2.

The performance of the numerical model is tested for the diffusion (v_in = 0) and for the convection-diffusion (v_in = 1) controlled growth. Figure 3 shows the results of the simulations at t = 130. The PF in Ω for v_in = 0 and v_in = 1 is shown on the top-left and top-right sub-figures, respectively. The refinement at the solid–liquid interface and de-refinement in the bulk of the solid phase is seen in the top-left figure. The melt velocity vectors are plotted when v_in = 1. Table 1 contains the simulation parameters used.

4.2 Assessment of the results

A dendrite grows equally fast in all four directions in diffusion-controlled growth, as seen in Figure 3. In diffusion-convection-controlled growth, the dendrite grows faster in the upstream direction and slower in the downstream direction. In contrast, the growth velocity in the direction perpendicular to the fluid flow appears similar to the diffusion-controlled case. The dendrite’s trunk is thicker and thinner in the upstream and downstream directions, respectively. One can also see how the west and east trunks are no longer symmetric. It is evident that a dendrite grows quicker in the directions that provide a faster release of latent heat. The fluid flow increases the temperature gradient in the melt in the upstream direction and decreases it in the downstream direction, as seen in the bottom-left in Figure 3. The absolute value of the velocity field is shown on the bottom-right of Figure 3. The melt slows down near the dendrite surface due to the no-slip boundary condition. It is largely accelerated near the east and west part of Γ as the dendrite occupies an increasingly larger portion of the computational domain, previously filled by the fluid.

Figure 4 shows the rescaled growth velocity vtip/vtip0 of the south, west and north dendrite tip as a function of time for v_in = 0 and v_in = 1. Velocity v_tip stands for the dimensionless velocity at a dendrite tip. Velocity vtip0 is the dimensionless Green’s function analytical velocity for v_in = 0, tabulated in Karma and Rappel (1998) as vtip0W0d0/(τ0DT)=0.017. The figure also shows the reference solutions for v_in = 1 reported by Beckermann et al. (1999). The steady-state growth velocity is very close to the analytical growth velocity in all three directions for v_in = 0. For v_in = 1, our results agree well with the reference results for the north and south tip and the analytical growth velocity for the west tip. The highest difference between our and reference results is observed for the north tip velocity between t = 25 and t = 75. Remember that the reference solution was obtained by solving the Navier–Stokes equations while we consider Stokes flow. Therefore, we cannot expect our results to converge to the exact same numerical values as those in Beckermann et al. (1999).

The model uses the following numerical parameters to obtain results from Figures 3 and 4: h = 0.8, N = 13, f_s = 5.12, f_ϵ = 0.16 and f_b = f_exe = 8. A scattered node distribution is generated in each quadtree sub-domain. The MRSM with f_s > 0 is used on the Γ, where the positioning of the source nodes is trivial. In the boundary nodes at the solid–liquid interface, the MRSM with f_s = 0, i.e. the MRS (Wen et al., 2017), is applied to avoid complications with positioning source points. Previous research shows that the MRS is suitable for solving Stokes flow with Dirichlet boundary conditions (Wen et al., 2017). The method’s free parameters’ impact on accuracy is thoroughly analysed in the following two sub-sections. First, in sub-section 4.3, the influence of h and N on the accuracy is investigated for the RBF-FD method when using either regular or scattered node distribution for the case of diffusion-controlled growth. Sub-section 4.4 investigates the influence of f_s, f_ϵ, f_exe and f_b on the accuracy in the MRSM for the case of convection-diffusion-controlled growth.

4.3 Diffusion-controlled growth (v_in = 0)

This subsection formally performs the same tests as in Dobravec et al. (2022, 2020), where the influence of node distribution, size of local support and node spacing on the accuracy for solving the non-preconditioned PF model is investigated. Here, we analyse how the RBF-FD method performs for solving the preconditioned PF model and select the appropriate free parameters, which we will use in the following sub-section in the analysis of the MRSM. Figure 5 shows the steady-state growth velocity as a function of h for three different values of N using regular (left) and scattered (right) node distribution in the case of non-rotated (top) and rotated (bottom) dendrite. The rotated dendrite is rotated for π/4 with respect to the coordinate system to analyse the influence of the discretisation-induced anisotropy.

We can see how the velocity converges towards the analytical solution when reducing h. For h = 0.4, the velocity agrees very well with the reference analytical solution using both node distributions for rotated and non-rotated dendrites. In the case of the non-rotated dendrite for h > 0.4, increasing N increases the accuracy using both node distributions. We observe the same behaviour in the case of scattered node distribution for the rotated dendrite. However, the behaviour when using regular node distribution in the case of the rotated dendrite is quite different. While the increase of N from N = 9 to N = 13 increases the accuracy, the increase of N from N = 13 to N = 21 decreases it.

In the case of the rotated dendrite, the results are much more sensitive to N when using regular node distribution. Similar results are also observed in the previous research (Dobravec et al., 2022, 2020); regular node distribution is much more prone to discretisation-induced anisotropy when considering growth in the arbitrary preferential growth direction. However, using the preconditioned PF model is more robust than the non-preconditioned PF model for both node distributions, especially for regular node distribution. For instance, for a similar test case with Δ = 0.65 and D = 1, which is analysed in Dobravec et al. (2022), the dendrite velocity of π/4-rotated dendrite is more than 40% higher compared to the non-rotated dendrite using N = 9, h = 0.8 and regular node distribution. In the present study, the worst-case deviation from the analytical velocity is only around 10% at h = 1.2.

The computational complexity of the numerical model increases with N and decreases with h. The configuration with h = 0.8, N = 13 and scattered node distribution is chosen to analyse the MRSM in the following sub-section. This configuration yields a good compromise between accuracy and computational efficiency. It takes around 30 s to finish the simulation with such configuration on an HP ZBook laptop with the hexacore Intel Core i7-9750H 2.6-4.5 GHz processor.

4.4 Convection-diffusion-controlled growth (v_in = 1)

In this sub-section, the influence of the MRSM’s free parameters on the accuracy is analysed; the RBF-FD method’s method parameters are fixed in this study (h = 0.8, N = 13 and scattered node distribution). The execution of the MRSM is a computationally expensive task. Each execution consists of constructing and solving the system of linear equations and calculating velocity at each computational node. A reduction of the number of executions and boundary nodes on Γ is needed to speed up the calculations.

Figure 6 on the left shows the relative difference between the tip velocity at f_exe > 1 compared to the velocity at f_exe = 1 for three growth directions. Naturally, increasing f_exe increases the difference. At f_exe = 8, the difference for the north tip is ≈10⁻³ while other directions experience lower errors. Figure 6 on the right shows the relative difference between the tip velocity at f_b > 1 compared to the velocity at f_b = 1 for three growth directions. As for f_exe, increasing f_b decreases the accuracy. Value f_b = 8 yields a difference below ≈10⁻³ for all three directions and is used in further calculations. The results suggest that values f_exe = f_b = 8 represent a good compromise between accuracy and computational efficiency. It takes around 6 min to finish the simulation with such parameters on an HP ZBook laptop with the hexacore Intel Core i7-9750H 2.6-4.5 GHz processor.

Figure 7 shows the velocity at the tip of a dendrite at t = 100 for three growth directions as a function of f_s at different values of f_ϵ. We can see that value f_ϵ = 5.12 yields too large an error in all three directions. Value f_ϵ = 2.56 yields good results for south and west directions at f_s ≥ 2.56; however, the error in the north direction is too high for f_s ≥ 2.56. For f_ϵ < 2.56, the growth velocity is no longer changing for f_s ≥ 2.56 in all three directions. This stagnation occurs for even lower values of f_s in the west and north directions. Our results are closest to the reference solutions from the literature in all three directions for f_s > 2.56 and f_ϵ < 1.28.

5. Conclusions

A novel numerical approach combining domain- and boundary-type meshless methods for the PF modelling of dendritic solidification with fluid flow is presented. This original approach uses the domain-type RBF-FD method for the spatial discretisation of PF and energy conservation equations. The boundary-type MRSM calculates the Stokes flow around evolving dendrite. The forward Euler scheme is used for the time-stepping of PF and energy conservation equations. The approach uses the space-time adaptive algorithm to accelerate the calculations. Non-linear preconditioning is applied to ensure stability when using larger node spacings. We first test the RBF-FD method in the case of diffusion-controlled growth. We next analyse the MRSM method in the case of convection-diffusion-controlled growth.

In the case of diffusion-controlled growth, we perform a similar analysis as Dobravec et al. (2022) and investigate the influence of free numerical parameters of the RBF-FD method on the accuracy. We observe similar behaviour as Dobravec et al. (2022); the accuracy increases with reduced node spacing h and the increased size of a local support domain N. We repeat the same analysis for a dendrite rotated for π/4 concerning the axes of the coordinate system to investigate the discretisation-induced anisotropy. As also seen in Dobravec et al. (2022), the method is more sensitive to N when using regular node distribution. The best results are for the rotated and non-rotated dendrite observed when using the minimum tested spacing h = 0.4. At that spacing, the results are closest to the reference analytical solution and almost independent of N for both node distributions. Spacing h = 0.4, however, yields long computational times. With increased h, the results depend more on N and the type of node distribution used. Configuration with h = 0.8 and N = 13 yields the same results for rotated and non-rotated dendrites using both node distributions and, therefore, represents a good compromise between accuracy and computational efficiency.

In the case of convection-diffusion controlled growth, we investigate the influence of the free parameters of the MRSM on accuracy. We first check how the accuracy is affected by executing the MRSM every f_exe-th iteration of the PF and energy conservation equation. The error introduced by this optimisation is below ≈10⁻³ for f_exe = 8. We next check how the accuracy depends on the boundary spacing parameter f_b, i.e. the ratio between the node spacing at the boundary of the computational domain and the solid–liquid interface. Value f_b = 8 yields error below ≈10⁻³. Using f_exe = f_b = 8 hugely reduces the computational time of a simulation while sustaining good accuracy. Finally, we analyse the influence of the free parameter controlling the distance between a boundary and a source point f_s and the free parameter controlling the shape of the blob function f_ϵ on the accuracy. The MRSM with f_s = 0 is applied at the solid–liquid interface to avoid problems with source node positioning. The method returns the best results for f_s > 2.56 and f_ϵ < 1.28.

The main originality and novelty of the present approach is the successful coupling between domain- and boundary-type meshless methods for modelling dendritic growth with Stokes flow. Our results agree well with the published reference results. The use of the MRSM for solving Stokes flow, space-time adaptive algorithm, and non-linear preconditioning of the PF provide a computationally efficient numerical tool. The approach has many free parameters, influencing the accuracy and computational efficiency. This paper has proposed a suitable selection of these parameters based on the performed numerical experiments. The numerical model can be straightforwardly extended to 3D using 3D regularised Stokeslets. Octree instead of quadtree algorithm has to be applied in 3D space-time adaptive approach. The RBF-FD method is, on the contrary, dimension-independent.

Figures

Figure 1.

A scheme of space-time adaptive meshless PF modelling of dendritic solidification

Figure 2.

A scheme of the computational domain, boundary nodes and source nodes in the MRSM

Figure 3.

The PF for v_in = 0 (top-left) and v_in = 1 (top-right); red and blue represent solid and liquid phases. The temperature (bottom-left) and the absolute value of velocity (bottom-right) for v_in = 1; red and blue represent high and low values of the fields

Figure 4.

Growth velocity at dendrite tips for v_in = 0 (left) and v_in = 1 (right)

Figure 5.

Growth velocity at t = 100 as a function of h for three values of N using regular (left) and scattered (right) node distribution in the case of non-rotated (top) and rotated (bottom) dendrite

Figure 6.

The relative difference between the tip velocity at f_exe = 1 and f_exe > 1 (left) and between f_b = 1 and f_b > 1 (right) for three growth directions at t = 100

Figure 7.

Growth velocity at t = 100 as a function of f_s for different values of f_ϵ in the north (top-left), south (top-right) and west (bottom) directions

Table 1.

Simulation parameters

Computational domain
Size of domain (L)	230.4
Physical problem
Strength of anisotropy (ϵ4)	0.05
Initial supercooling (Δ)	0.55
Center of nucleus (r0)	(0, 0)
Radius of nucleus (r0)	3
Prandtl number (Pr = µ/D)	23.1
Inlet velocity (vin)	1
PF model
Constant (α1)	0.8839
Constant (α2)	0.6267
Coupling parameter (λ)	4/α2
Dimensionless diffusivity (D)	4

Source: Beckermann et al. (1999)