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

Tadej Dobravec (Department of Fluid Dynamics and Thermodynamics, Faculty of Mechanical Engineering, University of Ljubljana, Ljubljana, Slovenia)
Boštjan Mavrič (Laboratory for Simulation of Materials and Processes, Institute of Metals and Technology, Ljubljana, Slovenia and Department of Fluid Dynamics and Thermodynamics, Faculty of Mechanical Engineering, University of Ljubljana, Ljubljana, Slovenia)
Rizwan Zahoor (Department of Fluid Dynamics and Thermodynamics, Faculty of Mechanical Engineering, University of Ljubljana, Ljubljana, Slovenia)
Božidar Šarler (Department of Fluid Dynamics and Thermodynamics, Faculty of Mechanical Engineering, University of Ljubljana, Ljubljana, Slovenia and Laboratory for Simulation of Materials and Processes, Institute of Metals and Technology, Ljubljana, Slovenia)

International Journal of Numerical Methods for Heat & Fluid Flow

ISSN: 0961-5539

Article publication date: 8 June 2023

Issue publication date: 22 June 2023

385

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.

Keywords

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

Publisher

:

Emerald Publishing Limited

Copyright © 2023, Tadej Dobravec, Boštjan Mavrič, Rizwan Zahoor and Božidar Šarler.

License

Published by Emerald Publishing Limited. This article is published under the Creative Commons Attribution (CC BY 4.0) licence. Anyone may reproduce, distribute, translate and create derivative works of this article (for both commercial & non-commercial purposes), subject to full attribution to the original publication and authors. The full terms of this licence may be seen at http://creativecommons.org/licences/by/4.0/legalcode


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 cp, and thermal conductivity k. The latent heat of melting and the melting temperature are denoted as Lm and Tm, 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 d0 is the thermal capillary length while α1 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 α2 stands for a constant of the PF model and DT = k/(ρcp) for the thermal diffusivity. The PF constants are equal to α1 = 0.8839 and α2 = 0.6267 (Karma and Rappel, 1998). The selection of free parameter λ has to yield W0 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) ψ=2tanh1(ϕ),
to increase numerical stability for larger node spacings. The energy conservation equation in terms of dimensionless temperature θ = (TTm)/(Lm/cp) reads as:
(4) θt=D2θ+1ϕ2(1+ϕ2ψtv·θ),
where D and v = (vx, vy) stand for the dimensionless thermal diffusivity and velocity, respectively. Dimensionless D and v are measured in units of W02/τ0 and W0/τ0, 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)=13ϵ4+4ϵ4(nx4+ny4),

and

(8) a(n)=16ϵ4|ψ|a(n)(nx(nx4+ny4nx2),ny(nx4+ny4ny2)),
where ϵ4 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 nh, the maximum number of different time steps nΔt, the overlapping parameter noverlap 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 Naug 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 {lri} 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Φ(|rlri|lh)+i=1NauglαN+ipi(rlrlh),

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 lwk 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 {ri}. The Stokeslets are singular at the source; therefore, the source points {si} must be positioned outside the computational domain. Each collocation point ri on the boundary has a corresponding source point si some distance away from ri in the outward-normal direction:
(19) si=ri+hifsnΓ,
where hi stands for the spacing between boundary nodes and nΓ for the outward-facing normal to the boundary of the computational domain at ri while fs is a free parameter controlling the distance between ri and si. 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ϵhi, where fϵ stands for the free parameter controlling the shape of a blob. The MRSM with fs = 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 ri, 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, nh, nΔt, mΩ, noverlap and the type of node distribution. We set noverlap = 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 (nh = 6, mΩ = 9), (nh = 5, mΩ = 9) and (nh = 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 Naug. We use fifth-degree polyharmonic splines (n = 5) and second-order augmentation with monomials (Naug = 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: fs and fϵ. We test the performance of the method for fs ∈ [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Γ = fbh, where fb is a free parameter. We test the performance for fb ∈ [1,64]. To save computational time, we execute the MRSM every fexe-th iteration of the PF and energy conservation equations, where fexe is a free parameter. We test the performance for fexe ∈ [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 r0 and the radius r0. We set the initial conditions for the PF and energy conservation as:

(22) ψ(r,t=0)=r0|rr0|,θ(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,vyx|Γeast,west=0,vy|Γsouth=0,v|solid-liquid=0,
where vin 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 (vin = 0) and for the convection-diffusion (vin = 1) controlled growth. Figure 3 shows the results of the simulations at t = 130. The PF in Ω for vin = 0 and vin = 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 vin = 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 vin = 0 and vin = 1. Velocity vtip stands for the dimensionless velocity at a dendrite tip. Velocity vtip0 is the dimensionless Green’s function analytical velocity for vin = 0, tabulated in Karma and Rappel (1998) as vtip0W0d0/(τ0DT)=0.017. The figure also shows the reference solutions for vin = 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 vin = 0. For vin = 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, fs = 5.12, fϵ = 0.16 and fb = fexe = 8. A scattered node distribution is generated in each quadtree sub-domain. The MRSM with fs > 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 fs = 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 fs, fϵ, fexe and fb on the accuracy in the MRSM for the case of convection-diffusion-controlled growth.

4.3 Diffusion-controlled growth (vin = 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 (vin = 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 fexe > 1 compared to the velocity at fexe = 1 for three growth directions. Naturally, increasing fexe increases the difference. At fexe = 8, the difference for the north tip is ≈10−3 while other directions experience lower errors. Figure 6 on the right shows the relative difference between the tip velocity at fb > 1 compared to the velocity at fb = 1 for three growth directions. As for fexe, increasing fb decreases the accuracy. Value fb = 8 yields a difference below ≈10−3 for all three directions and is used in further calculations. The results suggest that values fexe = fb = 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 fs 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 fs ≥ 2.56; however, the error in the north direction is too high for fs ≥ 2.56. For fϵ < 2.56, the growth velocity is no longer changing for fs ≥ 2.56 in all three directions. This stagnation occurs for even lower values of fs in the west and north directions. Our results are closest to the reference solutions from the literature in all three directions for fs > 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 fexe-th iteration of the PF and energy conservation equation. The error introduced by this optimisation is below ≈10−3 for fexe = 8. We next check how the accuracy depends on the boundary spacing parameter fb, i.e. the ratio between the node spacing at the boundary of the computational domain and the solid–liquid interface. Value fb = 8 yields error below ≈10−3. Using fexe = fb = 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 fs and the free parameter controlling the shape of the blob function fϵ on the accuracy. The MRSM with fs = 0 is applied at the solid–liquid interface to avoid problems with source node positioning. The method returns the best results for fs > 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

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

Figure 1.

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

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

Figure 2.

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

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

Figure 3.

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

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

Figure 4.

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

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 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

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

Figure 6.

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

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

Figure 7.

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

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)

References

Anderson, E., Bai, Z., Bischof, C., Blackford, S., Demmel, J., Dongarra, J., Croz, J.D., Greenbaum, A., Hammarling, S., McKenney, A. and Sorensen, D. (1987), LAPACK Users’ Guide, 3rd ed., Society for Industrial and Applied Mathematics, Philadelphia.

Atluri, S.N. (2004), The Meshless Method (MLPG) for Domain and BIE Discretizations, Tech Science Press, Forsyth, GA.

Bayona, V., Flyer, N., Fornberg, B. and Barnett, G.A. (2017), “On the role of polynomials in RBF-FD approximations: II. Numerical solution of elliptic PDEs”, Journal of Computational Physics, Vol. 332, pp. 257-273.

Beckermann, C., Diepers, H.J., Steinbach, I., Karma, A. and Tong, X. (1999), “Modeling melt convection in phase-field simulations of solidification”, Journal of Computational Physics, Vol. 154 No. 2, pp. 468-496.

Belytschko, T., Lu, Y.Y. and Gu, L. (1994), “Element-free Galerkin methods”, International Journal for Numerical Methods in Engineering, Vol. 37 No. 2, pp. 229-256.

Boettinger, W.J., Warren, J.A., Beckermann, C. and Karma, A. (2002), “Phase-field simulation of solidification”, Annual Review of Materials Research, Vol. 32 No. 1, pp. 163-194.

Boukellal, A.K., Rouby, M. and Debierre, J.-M. (2021), “Tip dynamics for equiaxed Al-Cu dendrites in thin samples: phase-field study of thermodynamic effects”, Computational Materials Science, Vol. 186, p. 110051.

Campbell, J. (2003), Castings, Butterworth-Heinemann, Oxford, UK.

Chapman, B., Jost, G. and Pas, R. V D (2008), Using OpenMP: Portable Shared Memory Parallel Programming, MIT Press, Cambridge, MA.

Chen, L.-Q. (2002), “Phase-field models for microstructure evolution”, Annual Review of Materials Research, Vol. 32 No. 1, pp. 113-140.

Cheng, A.H.D. and Hong, Y. (2020), “An overview of the method of fundamental solutions – solvability, uniqueness, convergence, and stability”, Engineering Analysis with Boundary Elements, Vol. 120, pp. 118-152.

Dobravec, T., Mavrič, B. and Šarler, B. (2017), “A cellular automaton — finite volume method for the simulation of dendritic and eutectic growth in binary alloys using an adaptive mesh refinement”, Journal of Computational Physics, Vol. 349, pp. 351-375.

Dobravec, T., Mavrič, B. and Šarler, B. (2020), “Reduction of discretisation-induced anisotropy in the phase-field modelling of dendritic growth by meshless approach”, Computational Materials Science, Vol. 172, p. 109166.

Dobravec, T., Mavrič, B. and Šarler, B. (2022), “Acceleration of RBF-FD meshless phase-field modelling of dendritic solidification by space-time adaptive approach”, Computers and Mathematics with Applications, Vol. 126, pp. 77-99.

Dong, X., Xing, H., Weng, K. and Zhao, H. (2017), “Current development in quantitative phase-field modeling of solidification”, Journal of Iron and Steel Research International, Vol. 24 No. 9, pp. 865-878.

Fasshauer, G.F. (2007), Meshfree Approximation Methods with MATLAB, World Scientific Publishing, River Edge, NJ.

Flyer, N., Fornberg, B., Bayona, V. and Barnett, G.A. (2016), “On the role of polynomials in RBF-FD approximations: I. Interpolation and accuracy”, Journal of Computational Physics, Vol. 321, pp. 21-38.

Gibou, F., Fedkiw, R., Caflisch, R. and Osher, S. (2003), “A level set approach for the numerical simulation of dendritic growth”, Journal of Scientific Computing, Vol. 19 Nos 1/3, pp. 183-199.

Glasner, K. (2001), “Nonlinear preconditioning for diffuse interfaces”, Journal of Computational Physics, Vol. 174 No. 2, pp. 695-711.

Gong, T.Z., Chen, Y., Cao, Y.F., Kang, X.H. and Li, D.Z. (2018), “Fast simulations of a large number of crystals growth in centimeter-scale during alloy solidification via nonlinearly preconditioned quantitative phase-field formula”, Computational Materials Science, Vol. 147, pp. 338-352.

Gu, Y. and Liu, G. (2003), “A boundary radial point interpolation method (BRPIM) for 2-D structural analyses”, Structural Engineering and Mechanics, Vol. 15 No. 5.

Gu, Y.T. and Liu, G.R. (2002), “A boundary point interpolation method for stress analysis of solids”, Computational Mechanics, Vol. 28 No. 1, pp. 47-54.

Hanoglu, U. and Šarler, B. (2018), “Multi-pass hot-rolling simulation using a meshless method”, Computers and Structures, Vol. 194, pp. 1-14.

Hatić, V., Cisternas Fernández, M., Mavrič, B., Založnik, M., Combeau, H. and Šarler, B. (2019), “Simulation of a macrosegregation benchmark in a cylindrical coordinate system with a meshless method”, International Journal of Thermal Sciences, Vol. 142, pp. 121-133.

Iserles, A. (2000), Acta Numerica, Cambridge University Press, Cambridge, UK.

Jeong, J.-H., Goldenfeld, N. and Dantzig, J.A. (2001), “Phase field model for three-dimensional dendritic growth with fluid flow”, Physical Review E, Vol. 64 No. 4, p. 41602.

Karma, A. and Rappel, W.-J. (1998), “Quantitative phase-field modeling of dendritic growth in two and three dimensions”, Physical Review E, Vol. 57 No. 4, p. 4323.

Karma, A. and Tourret, D. (2016), “Atomistic to continuum modeling of solidification microstructures”, Current Opinion in Solid State and Materials Science, Vol. 20 No. 1, pp. 25-36.

Kosec, G. and Šarler, B. (2011), “H-adaptive local radial basis function collocation meshless method”, Computers, Materials and Continua, Vol. 26 No. 3, pp. 227-253.

Kurz, W., Fisher, D.J. and Trivedi, R. (2019), “Progress in modelling solidification microstructures in metals and alloys: dendrites and cells from 1700 to 2000”, International Materials Reviews, Vol. 64 No. 6, pp. 311-354.

Kurz, W., Rappaz, M. and Trivedi, R. (2021), “Progress in modelling solidification microstructures in metals and alloys. Part II: dendrites from 2001 to 2018”, International Materials Reviews, Vol. 66 No. 1, pp. 30-76.

Liu, G.R. (2009), Meshfree Methods: Moving Beyond the Finite Element Method, 2nd ed., CRC Press, Boca Raton, FL.

Liu, G.R. and Gu, Y.T. (2001), “A local radial point interpolation method (LRPIM) for free vibration analyses of 2-D solids”, Journal of Sound and Vibration, Vol. 246 No. 1, pp. 29-46.

Liu, G.R. and Gu, Y.T. (2005), An Introduction to Meshfree Methods and Their Programming, Springer Netherlands, Dordrecht, Netherlands.

Liu, Q.G. and Šarler, B. (2018), “Non-singular method of fundamental solutions for elasticity problems in three-dimensions”, Engineering Analysis with Boundary Elements, Vol. 96, pp. 23-35.

Mavrič, B. and Šarler, B. (2015), “Local radial basis function collocation method for linear thermoelasticity in two dimensions”, International Journal of Numerical Methods for Heat and Fluid Flow, Vol. 25 No. 6, pp. 1488-1510.

Mramor, K., Vertnik, R. and Šarler, B. (2014), “Simulation of laminar backward facing step flow under magnetic field with explicit local radial basis function collocation method”, Engineering Analysis with Boundary Elements, Vol. 49, pp. 37-47.

Provatas, N. and Elder, K. (2010), Phase-Field Methods in Materials Science and Engineering, Wiley-VCH, Weinheim, Germany.

Rek, Z. and Šarler, B. (2021), “The method of fundamental solutions for the Stokes flow with the subdomain technique”, Engineering Analysis with Boundary Elements, Vol. 128, pp. 80-89.

Rek, Z., Zahoor, R. and Šarler, B. (2021), “Modified method of regularized sources for potential flow”, Computers and Mathematics with Applications, Vol. 88, pp. 110-119.

Reuther, K. and Rettenmayr, M. (2014), “Perspectives for cellular automata for the simulation of dendritic solidification — a review”, Computational Materials Science, Vol. 95, pp. 213-220.

Sadat, H. and Prax, C. (1996), “Application of the diffuse approximation for solving fluid flow and heat transfer problems”, International Journal of Heat and Mass Transfer, Vol. 39 No. 1, pp. 214-218.

Šarler, B. (2006), “Solution of a two-dimensional bubble shape in potential flow by the method of fundamental solutions”, Engineering Analysis with Boundary Elements, Vol. 30 No. 3, pp. 227-235.

Šarler, B. and Vertnik, R. (2006), “Meshfree explicit local radial basis function collocation method for diffusion problems”, Computers and Mathematics with Applications, Vol. 51 No. 8, pp. 1269-1282.

Steinbach, I. (2009), “Phase-field models in materials science”, Modelling and Simulation in Materials Science and Engineering, Vol. 17 No. 7, p. 73001.

Talat, N., Mavrič, B., Belšak, G., Hatić, V., Bajt, S. and Šarler, B. (2018), “Development of meshless phase field method for two-phase flow”, International Journal of Multiphase Flow, Vol. 108, pp. 169-180.

Tan, L. and Zabaras, N. (2006), “A level set simulation of dendritic solidification with combined features of front-tracking and fixed-domain methods”, Journal of Computational Physics, Vol. 211 No. 1, pp. 36-63.

Vertnik, R., Mramor, K. and Šarler, B. (2019), “Solution of three-dimensional temperature and turbulent velocity field in continuously cast steel billets with electromagnetic stirring by a meshless method”, Engineering Analysis with Boundary Elements, Vol. 104, pp. 347-363.

Wang, K., Wen, S., Zahoor, R., Li, M. and Šarler, B. (2016), “Method of regularized sources for axisymmetric stokes flow problems”, International Journal of Numerical Methods for Heat and Fluid Flow, Vol. 26 Nos 3/4, pp. 1226-1239.

Wen, S., Wang, K., Zahoor, R., Li, M. and Šarler, B. (2017), “Method of regularized sources for two-dimensional Stokes flow problems based on rational or exponential blobs”, Computer Assisted Methods in Engineering and Science, Vol. 22 No. 3, pp. 289-300.

Zhu, T., Zhang, J. and Atluri, S.N. (1998), “A meshless local boundary integral equation (LBIE) method for solving nonlinear problems”, Computational Mechanics, Vol. 22 No. 2, pp. 174-186.

Acknowledgements

The Slovenian Research Agency (ARRS) supported this work under projects Z2-4479 (TD), Z2-2640 (BM), P2-0162, J2-4477 (RZ), and L2-3173 supported also by Štore-Steel company (BŠ). We thank Dr Zlatko Rek for valuable discussions regarding the boundary-type meshless numerical method.

Corresponding author

Božidar Šarler can be contacted at: bozidar.sarler@fs.uni-lj.si

Related articles