Xiaofeng Yang

Department of Mathematics,University of South Carolina,Columbia,SC 29208,USA.

Abstract.The phase-field dendritic crystal growth model is a highly nonlinear system that couples the anisotropic Allen-Cahn type equation and the heat equation.By combining the recently developed SAV(Scalar Auxiliary Variable)method with the linear stabilization approach,as well as a special decoupling technique,we arrive at a totally decoupled,linear,and unconditionally energy stable scheme for solving the dendritic model.We prove its unconditional energy stability rigorously and present various numerical simulations to demonstrate the stability and accuracy.

Key words:Phase-field,dendritic,stabilized-SAV method,anisotropy,Allen-Cahn,decoupled.

1 Introduction

The use of the phase-field method for investigating the process of dendritic crystal growth can be attributed to the pioneering modeling work by Halperin,Kobayashi,and Collins et.al.in[1–3],and see also the subsequent modeling/simulations in[4–15].In a typical phase-field dendritic crystal system,an order parameter(called phase-field variable)is usually introduced to define the physical state(liquid or solid)at each point and the total free energy incorporates a specific form of the conformational entropy with anisotropic spatial gradients.The system usually consists of two coupled nonlinear,second-order equations:the Allen-Cahn type equation with a gradient-dependent anisotropic coefficient,and the heat transfer equation.

In this paper,we consider numerical approximations for a phase-field dendritic crystal growth model which was proposed by Karma and Rappel in[12].It is well known that the main objective of algorithm design for phase-field related models is to construct efficient and easy-to-implement numerical schemes that can verify a discrete energy law.For the particular dendritic model proposed in[12],the associated difficulties to this aim lie on how to discretize three nonlinear terms,including the anisotropic coefficient,the cubic polynomial term,as well as the heat transfer term.Simple explicit treatment for these nonlinear terms will induce large spatial oscillations that may cause the computations easily blow up or loss of accuracy(shown in Figure 2,Figure 4,and Figure 5(b)).

We recall there exist plenty of time discretization methods that had been proved to be effective for solving the phase field models,see[15–38].However,for this particular model considered in this paper,most of the available schemes are either nonlinear which need some efficient iterative solvers,and/or do not preserve energy stability at all(cf.[7,39–44]and the references therein).Therefore,in this paper,by combining the recently developed SAV(Scalar Auxiliary Variable)method with the linear stabilization approach,as well as a special decoupling technique,we arrive at a fully-decoupled,stabilized-SAV scheme.The novelty of this scheme is that two linear stabilization terms are added in the SAV scheme,where one is used to remove the oscillations caused by the anisotropic coefficient,and the other is added to the latent heat transfer term in order to realize the decoupling.At each time step,one can only solve an elliptic system for the phase function,and a linear elliptic equation for the temperature.We then prove that the unconditionally energy stability of the scheme and present numerous numerical examples to illustrate its accuracy and stability numerically.

The rest of the paper is organized as follows.In Section 2,we give a briefintroduction of the governing PDE system for the phase-field anisotropic dendritic crystal growth model.In Section 3,we develop the scheme for solving the model,and rigorously prove the unconditional energy stability.Various numerical experiments are given in Section 4 to demonstrate the accuracy and efficiency of the proposed numerical scheme.Finally,some concluding remarks are given in Section 5.

2 Model equations

We give a brief description of the anisotropic phase-field dendritic crystal growth model proposed in[12].LetΩbe a smooth,open,bounded,connected domain in Rdwithd=2,3.A scalar phase-field functionφ(x,t)is introduced to label the liquid and solid phase,whereφ=1 for the solid andφ=-1 for the fluid.These two regions are connected by a smooth transitional layer with the thickness∈.The total free energy is postulated as follows,

in which,T(x,t)is the temperature,∈,λandKare all positive parameters,F(φ)=(φ2-1)2is the Ginzburg-Landau double well potential,κ(?φ)is a function describing the anisotropic property that depends on the direction of the outer normal vectornwhich is the interface normal defined as.For the 2D system,the anisotropy coefficientκ(?φ)is usually given by

The last two terms in(2.4)and(2.5)are two postulated terms which are not derived by the variational derivative of the total free energy.The functionp(φ)accounts for the generation of latent heat and it is a phenomenological functional taking the form preserving the minima ofφat±1 independently of the local value ofT.Forp(φ),there are two common choices:andp′(φ)=(1-φ2)2(cf.[2,12]);orandp′(φ)=1-φ2(cf.[45]),that imply the heat only transfers through the interface.

Without the loss of generality,one can adopt the periodic boundary condition or the no-flux homogenous Neumann boundary conditions in order to remove all complexities associated with the boundary integrals in this study,i.e.,

where n is the outward normal of the computational domainΩ.

The model equations(2.4)-(2.5)follows the dissipative energy law.By taking theL2inner product of(2.4)withφt,and of(2.5)with,using the integration by parts and combining the obtained two equalities,we obtain

3 Numerical schemes

The aim of this section is to construct decoupled,energy stable schemes to solving the system(2.4)-(2.5).It has been shown that spurious solutions may occur if a numerical scheme does not satisfy the discrete energy dissipation law when the spatial grid and time step sizes are not carefully chosen.In addition,with unconditionally energy stable schemes,one can use relatively large time steps,the size of which is dictated only by accuracy considerations,or a suitable adaptive time stepping.The system(2.4)-(2.5)is a coupled nonlinear model.While it is relatively easy to design some fully implicit schemes with the energy stability,it is very difficult to maintain the unconditional energy stability together with the decoupled feature,that can be seen from our recent paper[46]where we have established a linear and stable,but coupled scheme.

To realize the linear and unconditionally energy stable features,we combine the SAV approach developed in[34]with the stabilization techinque to treat the nonlinear gradient potential and the double-well potential.To realize the decoupling feature,we introduce an intermediate temperature that is the combination of the temperature at previous time step and an extra semi-explicit stabilizing term that includes the heat transfer term.As a consequence,the nonlinear heat transfer term in the heat equation vanishes since it is already included in the intermediate temperature.Finally,a decoupled,well-posed linear scheme is obtained and it can be proved to be unconditionally energy stable theoretically and numerically.

3.1 First order scheme

We construct a first-order time marching scheme for solving the system(3.3)-(3.5)based on the first-order backward formula,shown as follows.

3.2 Second order scheme

A formally decoupled,second-order scheme can be constructed similar to the firstorder scheme,which is based on the backward differentiation formula(BDF2)that reads as follows.

4 Numerical simulations

In this section,we present various numerical examples to validate the proposed schemes and demonstrate their accuracy,energy stability,and efficiency numerically.

4.1 Accuracy test

We first implement a numerical example with fourfold anisotropy(2.3)in 2D space Ω=[0,h1]×[0,h2]to test the convergence rates of the proposed first-order scheme(3.8)-(3.11),denoted by SSAV;and the second-order scheme(3.30)-(3.33),denoted by SSAV-2nd.For comparisons,we also compute the convergence rates by using the non-stabilized version of scheme,i.e.,scheme(3.8)-(3.11)but withS1=S2=0,denoted by SAV.

We assume the following two functions

to be the exact solutions,and impose some suitable force fields such that the given solutions can satisfy the system(2.4)-(2.5).The model parameters are set as follows,

We discretize the space usingNx=Ny=129 Fourier modes forxandydirections so that the errors from the spatial discretization are negligible compared to the temporal discretization errors.

In Figures 1 and 2,by using the the first-order scheme SSAV and its non-stabilized version SAV scheme,we plot theL2errors of the phase and temperature variables between the numerical solution and the exact solution att=0.2 with different time step sizes by varying the parameterτ0decreasingly.Some observed features are listed as follows.

·In Figure 1,we setand plot the errors computed by using the two firstorder schemes,SSAV and SAV.We observe these schemes not only present the good convergence rate that almost perfectly matches the first-order accuracy for the time step but also good approximations to the exact solution,regardless of whether they are stabilized or not.

Figure 4:The profile ofγ(n0)with∈4=0.25 and the initial condition(4.3).The two subfigures are from the different view angle.

Figure 5:(a)The 2D dynamical evolution of the phase variableφby using the initial condition(4.3)and the time stepδt=1e-2 where snapshots of the numerical approximation are taken at t=0,2,4,6,and 8.(b)-(c)Time evolution of the free energy functional(2.1)when using three combinations of linear stabilizers where(I)S1=S2=0;(II)S1=4,S2=0;and(III)S1=4,S2=4.

·In Figure 2,we setand plot the errors computed by using the schemes SSAV and SAV.We observe that the non-stabilized scheme SAV totally loses the accuracy for all tested time steps.On the contrary,the scheme SSAV is stable for all tested time steps and perform good approximations and the first-order accuracy all along.

Therefore,through these numerical tests,we conclude that(i)if the mobility parameterτ0is large,all schemes can solve the model well;(ii)if the mobility parameterτ0is small,the stabilized scheme SSAV overwhelmingly defeats its non-stabilized version SAV from the stability and/or accuracy.

4.2 Evolution of a circle

In this subsection,we consider the evolution of a circle driven by the fourfold anisotropy.We set the mobility parameterτ0=1,the anisotropic strength parameter∈4=0.25,and still use other model parameters from the previous example.The initial conditions forφandTread as

Figure 6:The dynamical evolutions of the dendritic crystal growth process with fourfold anisotropy where K=0.5 and default parameters of(4.5),computed by using the scheme SSAV and the time stepδt=1e-2.Snapshots of the numerical approximation are taken at t=0,40,60,80,and 120.

where(x0,y0,r0,∈0,T0)=(π,π,1.5,0.072,-0.55).The time step isδt=0.01 and the space is discretized by using 1292Fourier modes.

In Figure 4,we present the 2D profile ofγ(n0),where one can observe high oscillations appear almost everywhere.In Figure 5,we show the dynamics of how a circular shape interface with full orientations evolves to an anisotropic diamond shape with missing orientations at four corners.Snapshots of the phase field variableφare taken att=0,2,4,6,and 8.To show how the two stabilizersS1andS2can suppress high-frequency oscillations efficiently,in Figure 5(b)and(c),we test the performance of three combinations of stabilizers:(I)S1=0,S2=0;(II)S1=4,S2=0;and(III)S1=4,S2=4 by using the same time stepδt=0.01.With the stabilizers in(I)and(II),the energies increase with time and the computations blow up very quickly.With the stabilizer(III),the high-frequency spatial oscillations are efficiently eliminated and the energy decays monotically.

4.3 2D dendrite crystal growth

In this subsection,we investigate the dynamics on how the anisotropic coefficient affects the shape of the dendritic crystal.We initially deposit a single small crystal nucleus in center of the computed domain and observe how it grows heterogeneously.This is a benchmark simulation that had been extensively studied in[2,12,14,15].

Figure 7:The dynamical evolutions of the dendritic crystal growth process with fourfold anisotropy where K=0.6 and default parameters of(4.5),computed by the using the scheme SSAV and the time stepδt=1e-2.Snapshots of the numerical approximation are taken at t=0,40,60,80,and 120.

Figure 8:The dynamical evolutions of the dendritic crystal growth process with fourfold anisotropy where K=0.7 and default parameters of(4.5),computed by using the scheme SSAV and the time stepδt=1e-2.Snapshots of the numerical approximation are taken at t=0,40,60,80,and 120.

The initial condition reads as

where(x0,y0,r0,∈0,T0)=(π,π,0.02,0.072,-0.55),and the model parameters read as follows,

Figure 9:The dynamical evolutions of the dendritic crystal growth process with fourfold anisotropy where K=0.8 and default parameters of(4.5),computed by the using scheme SSAV and the time stepδt=1e-2.Snapshots of the numerical approximation are taken at t=0,40,60,80,and 120.

We discretize the space usingNx=Ny=513 Fourier modes.

4.3.1 Fourfold anisotropy

We first investigate the fourfold anisotropy case by varying the parameterKand fixing all other parameters from(4.5).We use the first-order scheme SSAV,and the time stepδtis set to be 0.001 for better accuracy.In Figure 6(a),we show snapshots of the phase variable at various times withK=0.5.We observe that the tiny nucleus grows up with time and finally evolves to an anisotropic shape with missing orientations at four corners due to the anisotropic effect.In Figure 6(b),snapshots of the temperature variableTare presented.The interfacial contour of the temperature actually agrees well with the phase-field interface since the latent heat transfers only through the interface.

We then increase the parameterKto 0.6,0.7,and 0.8 in Figure 7,Figure 8,and Figure,9,respectively.WhenK=0.6,in Figure 7,we observe that the tiny nucleus finally grows into a star-shape with four fat branches.WhenK=0.7 and 0.8,star-shapes with much sharper tips and thinner branches are formed.To get more detailed evolution of the dendrites,in Fig 10(a),we summarize the contour of the interface{φ=0}every 20 time units from the initial moment for the above four cases.It can be seen that the shape of the dendrites is significantly influenced by the parameterK.All these numerical results demonstrate similar features to those obtained in[2,8,12,15].In Figure 10(b)and(c),for the above four cases,we summarize the evolutions of the logarithm of the total free energy which monotonically decays and the radius of the crystal that is measured by an equivalent radius of a circle with the same area.We observe that the decaying speed of the total energy and growing speed of the area become slower when the parameterKis larger.

Figure 10:(a)The summary of the contour of the interface{φ=0}every 20 time units from the initial moment for the four different parameter K,where,(a)K=0.5,(b)K=0.6,(c)K=0.7,and(d)K=0.8.(b)Time evolutions of the logarithm of the free energy functional.(c)The size of the dendritic crystals changing with time,where the crystal size is measured by an equivalent radius of a circle with the same area.

Figure 11:(a)The summary of the contour of the interface{φ=0}every 20 time units from the initial moment for the four anisotropic strengths∈4 and K=0.8,where,(a)∈4=0.01,(b)∈4=0.03,(c)∈4=0.05,and(d)∈4=0.07.(b)Time evolutions of the logarithm of the free energy functional.(c)The size of the dendritic crystals changing with time,where the crystal size is measured by an equivalent radius of a circle with the same area.

Figure 12:Time evolutions of the total free energy(2.1)computed by six different time step sizes until t=10 by using(a)the first-order SSAV scheme and(b)the second-order SSAV-2nd scheme.

Figure 13:The dynamical evolutions of the dendritic crystal growth process with sixfold anisotropy where K=0.6 and default parameters of(4.5),computed by using the scheme SSAV-2nd and the time stepδt=1e-2.Snapshots of the numerical approximation are taken at t=20,40,60,120,and 160.

We further investigate the effects of the anisotropic lengths∈4by fixingK=0.8 and all other parameters from(4.5).In Figure 11,we summarize the contour of the interface{φ=0}every 20 time units from the initial moment for four cases of∈4=0.01,∈4=0.03,∈4=0.05,and∈4=0.07,respectively.We observe that the thickness of formed branches are affected by the magnitude of∈4,and larger value of it can bring more slender pattern of branches.

Figure 14:The dynamical evolutions of the dendritic crystal growth process with sixfold anisotropy where K=0.65 and default parameters of(4.5),computed by the scheme SSAV-2nd and the time stepδt=1e-2.Snapshots of the numerical approximation are taken at t=20,40,60,120,and 160.

Figure 15:The dynamical evolutions of the dendritic crystal growth process with sixfold anisotropy where K=0.65 and default parameters of(4.5),computed by the scheme SSAV-2nd and the time stepδt=1e-2.Snapshots of the numerical approximation are taken at t=20,40,60,120,and 160.

We finally show the two developed schemes,the first-order SSAV scheme(3.8)-(3.11),and the second-order SSAV-2nd(3.30)-(3.33),are unconditionally energy stable numerically.In Figure 12,we plot the evolution curves of the total free energy(2.1)computed by six different time step sizes untilt=10 for∈4=0.05 andK=0.5.For all tested time steps,the obtained energy curves show the monotonic decays that confirm these two algorithms are unconditionally stable.We also observe that,for the first-order scheme SSAV,whenδt≤0.00125,the energy curves coincide very well.But whenδt＞0.00125,the energy curves deviate viewable away from others.This means the adopted time step size for the first-order scheme SSAV should not be larger than 0.00125,in order to get reasonably good accuracy(we setδt=0.001 in the computations).For second-order scheme SSAV-2nd,all energy curves coincide very well which implies the second-order scheme can provide much better accuracy than the first-order scheme.

Figure 16:The dynamical evolutions of the dendritic crystal growth process with sixfold anisotropy where K=0.75 and default parameters of(4.5),computed by using the scheme SSAV-2nd and the time stepδt=1e-2.Snapshots of the numerical approximation are taken at t=20,40,60,120,and 160.

Figure 17:(a)The summary of the contour of the interface{φ=0}every 20 time units from the initial moment for the sixfold anisotropy by using the four different parameter K,where,(a)K=0.6,(b)K=0.65,(c)K=0.7,and(d)K=0.75.(b)Time evolutions of the logarithm of the free energy functional.(c)The size of the dendritic crystals changing with time,where the crystal size is measured by an equivalent radius of a circle with the same area.

Figure 18:The dynamical evolutions of the dendritic crystal growth process with the sixfold anisotropy and K=0.75 computed by the scheme SSAV-2nd and the time stepδt=1e-2 where three randomly deposited nuclei are used as initial conditions.Snapshots of the numerical approximation are taken at t=0,40,60,100,and 140.

4.3.2 Sixfold anisotropy

In this subsection,we consider the the sixfold anisotropy by settingm=6 in(2.2).To compare with the fourfold case,we vary the parameterKwith the same initial condition(4.4)and order parameters from(4.5).We use the second-order scheme SSAV-2nd,and the time stepδtis set to be 0.01 for better accuracy.

We set the latent heat parameterKto be 0.6,0.7,0.75,and 0.8 in Figure 13–16,respectively.In Figure 13 withK=0.6,the circular nucleus initially grows into the regular hexagon shape,and then it becomes a snow-flake shape,but with few subtle microstructures.WhenK=0.65,in Figure 14,the obtained snow-flake pattern is filled with numerous subtle microstructures which are formed due to the anisotropy in the heat transfer process.When we further increase the parameterKin Figure 15 and Figure 16 withK=0.7 andK=0.75,respectively,we observe the less subtle microstructures and sub-branches.These simulations are qualitatively consistent with[2]where a slightly different model was used.To get more detailed evolution of the dendrites,in Figure 17(a),we summarize the contour of the interface{φ=0}every 20 time units from the initial moment for the above four cases,and present the evolutions of the logarithm of the total free energy and the radius of the crystal in Figure 17(b)and(c),respectively.

Finally,in Figure 18,three tiny nuclei are randomly deposited initially in the computed domain with the sixfold anisotropy andK=0.75.We observe that three dendrites are finally formed but with plenty of squeezed branches during the formation process.In the last subfigure of each figure,we show the corresponding temperature fieldT.

Figure 19:2D dynamical evolutions ofφfor the spinodal decomposition example,where the initial condition is(4.6)and m=0.Snapshots are taken at t=100,200,300,400,1000.

Figure 20:2D dynamical evolutions ofφfor the spinodal decomposition example,where the initial condition is(4.6)and m=3.Snapshots are taken at t=100,200,300,400,1000.

4.4 Spinodal decompostion

In this example,we study the phase separation dynamics that is called spinodal decomposition using the developed second-order scheme SSAV-2nd.By considering a homogeneous binary mixture,the spontaneous growth of the concentration fluctuations can lead the system from the homogeneous to the two-phase state.

We set the initial condition as the randomly perturbed concentration field as follows,

where the rand(x,y)is the random number in[-1,1]that follows the normal distribution.We set the time stepδt=0.001 for better accuracy,and discretize the space usingNx=Ny=257 Fourier modes.The model parameters are set as follows,

In Figure 19,Figure 20,Figure 21,and Figure 22,we setm=0,3,4,and 5,respectively,and present the profiles of the phase field variableφ.From the shapes of each case,we conclude that the number of sides of the obtained polygons is determined bymexclusively.In Figure 23,we present the evolution of the total free energy functional(2.1)for these cases together.

5 Concluding remarks

Figure 21:2D dynamical evolutions ofφfor the spinodal decomposition example,where the initial condition is(4.6)and m=4.Snapshots are taken at t=100,200,300,400,1000.

Figure 22:2D dynamical evolutions ofφfor the spinodal decomposition example,where the initial condition is(4.6)and m=5.Snapshots are taken at t=100,200,300,400,1000.

Figure 23:(a)Time evolution of the total free energy functional(2.1)for the spinodal decomposition example with four different parameters of m=0,3,4,and 5.(b)A close-up view for t∈[0,30].

In this paper,we have developed two totally decoupled,linear,and unconditionally energy stable schemes for solving the anisotropic dendritic phase-field model.The schemes are developed by combining the stabilized-SAV approach and a novel decoupling technique.We add two linear stabilization terms that can suppress high-frequency oscillations caused by the anisotropic coefficient efficiently,and also introduce an intermediate temperature that help to realize the total decoupling.Compared to the existed schemes for the anisotropic model,the proposed decoupled schemes(i)conquer the inconvenience from nonlinearities by linearizing the nonlinear terms in a new way,(ii)possess advantages of easy implementations and lower computational cost,(iii)are provably and/or numerically unconditionally energy stable,and thus(iv)allow for large time steps in computations.We further numerically verify the accuracy in time and present numerous numerical results for some benchmark numerical simulations.