- Email: [email protected]
A discrete kinetic scheme to model anisotropic liquid–solid phase transitions
Applied Mathematics Letters 103 (2020) 106222
Contents lists available at ScienceDirect
Applied Mathematics Letters www.elsevier.com/locate/aml
A discrete kinetic scheme to model anisotropic liquid–solid phase transitions Dongke Sun School of Mechanical Engineering, Southeast University, Nanjing 211189, China
article
info
Article history: Received 20 December 2019 Received in revised form 5 January 2020 Accepted 5 January 2020 Available online 13 January 2020 Keywords: Discrete kinetic scheme Lattice Boltzmann Phase field Phase transition
abstract An extended Boltzmann equation is constructed to model phase transitions with anisotropic effects. It is directly transformed into a discrete kinetic scheme in both time and phase space. The new scheme can be used to derive anisotropic lattice Boltzmann-phase field equation and convective diffusion equations. A lattice Boltzmann algorithm with the anisotropic streaming-relaxation operator is proposed, and numerical simulations of liquid–solid phase transition are performed as examples. The results agree well with previous numerical data and analytical solutions. It demonstrates that the scheme has acceptable numerical accuracy and computational efficiency, and the scheme can be used to study anisotropic liquid– solid phase transitions. This work provides an alternative numerical approach to model and simulate phase transitions with relatively fast computational efficiency. © 2020 Elsevier Ltd. All rights reserved.
1. Introduction Numerical methods based on the kinetic theory which bridges micro- and macro-physics have achieved great success both in academic and industrial areas [1–3]. Among these methods, the lattice Boltzmann (LB) method [4,5] and the discrete Boltzmann (DB) method [6,7] can be considered as a special finite difference scheme for discretization of the Boltzmann equation. For some distinguished advantages, these methods are being developed along two complementary directions. The LB scheme is widely used to solve nonlinear partial differential equations (PDEs), such as convective-diffusion equations [8], Allen–Cahn equation [9], phase field (PF) equations [10,11], and fractional PDEs [12,13]. The DB method is adopted to kinetically model non-equilibrium behaviors of fluids and reveal underlying physical features of various complex flows [14]. Phase transitions with anisotropic effects have been studied for decades by various phenomenological modeling and simulation techniques [15–17]. However, little research has been done on modeling of phase transitions by fitting the kinetic theoretical framework. Therefore, we attempt to construct a kinetic theory based equation to describe anisotropic phase transitions, and make a priori derivation from the kinetic equation to macroscopic phase field and convective diffusion equations. A discrete kinetic scheme is proposed E-mail address: [email protected]. https://doi.org/10.1016/j.aml.2020.106222 0893-9659/© 2020 Elsevier Ltd. All rights reserved.
D. Sun / Applied Mathematics Letters 103 (2020) 106222
2
to solve these equations, and the liquid–solid phase transition is selected as a benchmark to examine the present scheme. Finally, a brief conclusion closes the paper. 2. Extended kinetic equation and discrete kinetic scheme According to the kinetic theory, fluid flows can be described by the Boltzmann equation with the particle distribution function (PDF), f (x, ξ, t), ∂f + ξ · ∇f + a · ∇ξ f = Ω (f ), (1) ∂t where f is defined as the probability of finding a particle moving with velocity ξ at position x and time t, and a is the acceleration evoked by external forces. Ω (f ) denotes a collision operator, and it is usually approximated by some simplified models. One of the most popular choices is the Bhatnagar–Gross–Krook (BGK) approximation, which is given as ) 1 ( f − f (eq) . (2) Ω (f ) = − τc Here, τc represents the relaxation time, f (eq) is the equilibrium PDF. To extend the Boltzmann equation to model anisotropic phase transitions, the particle distribution function is redefined as the finding probability of a molecule moving with ξ at x and t. It is similar as the PDF defined in the Boltzmann equation for fluid flows, but the new PDF is used to describe phases rather than fluids. As phase transitions usually occur with anisotropic properties, an excess factor, an , is introduced to describe phase changes evoked by the anisotropic effect. Besides the collision of particles, the particle should satisfy mass conservation, no matter there exist phase changes or not. The continuous equation of the kinetic model with an is constructed. Its tensor format is written as ) ∂f 1 ( ∂f ∂f + ξβ + aβ =− f − f (eq) . an (3) ∂t ∂xβ ∂ξβ τc Here, an is a function of interface normal vector n, and we have an = an (n). Eq. (3) can be considered as an extension of Boltzmann equation which plays a central role in the kinetic theory. The particle velocity ξ can be split into two parts, the macroscopic velocity u and the relative velocity ξ, i.e. ξ = u + ζ. Here, u is independent of velocity space while ζ depends on velocity space. Multiplying equation (3) by ξα2 /2 and letting S = ξα2 f /2, S (eq) = ξα2 f (eq) /2 yield ) ∂S ∂S 1 ( ∂S + ξβ + aβ = S − S (eq) , (4) an ∂t ∂xβ ∂ξβ τc Integrating the equation over velocity space, result in ∫ ∫ ∫ ∫ ( ) ∂ ∂ ∂S 3 1 3 3 an Sd ξ + ξβ Sd ξ + aβ dξ = ξα2 S − S (eq) d3 ξ. ∂t ∂xβ ∂ξβ τc
(5)
As the collision conserves total energy, the term at the right hand side of the above equation is 0. Setting ∫ ψ = Sd3 ξ yields the macroscopic equation an
∂ψ ∂(uβ ψ) ∂qβ + =− + Sψ ∂t ∂xβ ∂xβ
with
Sψ =
∂(uα σαβ ) + ρaβ uβ . ∂xβ
(6)
Here, σαβ is the stress tensor given as σαβ and qβ is the flux calculated as,
∂ =− ∂xβ
1 qβ = 2
∫
∫
ζα ζβ f d3 ξ,
ζα ζα ζβ f d3 ξ.
(7)
(8)
D. Sun / Applied Mathematics Letters 103 (2020) 106222
3
When compressive and deformation work is neglected, we have ∂(uα σαβ )/∂xβ = 0 and then Sψ = ρaβ uβ . The flux qβ is usually related to the gradient of a scalar variable in a certain field. One choice is setting qβ = −ηp ∂ψ/∂xβ with ηp as the diffusion coefficient. Correspondingly, Eq. (6) can be written as ( ) ∂ψ ∂(uβ ψ) ∂ ∂ψ an + = ηp + Sψ . (9) ∂t ∂xβ ∂xβ ∂xβ Substitution of an = a2s (n), ψ = ϕ, uβ ψ = W02 Nβ /τ0 , ηp = W02 a2s (n)/τ0 , and Sψ = Sϕ /τ0 in Eq. (9) results in ( ) ∂ϕ ∂ϕ 2 ∂Nβ 2 ∂ 2 2 − W0 = W0 as (n) + Sϕ (10) τ0 as (n) ∂t ∂xβ ∂xβ ∂xβ where N and n are the anisotropic vector and the interfacial normal vector, respectively. Finally, a general formula of the phase field equation is derived from the extended kinetic equation, i.e. Eq. (3). The convective diffusion equation without anisotropic effect for heat transfer can also be derived in the same way. Setting an = 1, ψ = θ, ηd = κh , and Sψ = Sθ , we obtain ( ) ∂ ∂θ ∂(uβ θ) ∂θ + = κh + Sθ , (11) ∂t ∂xβ ∂xβ ∂xβ where θ denotes dimensionless temperature, and κh is diffusion coefficient for heat transfer. In the absence of flow, u = 0 and Eq. (11) is reduced to the pure diffusion equation with a source term. Now, we start to derive the anisotropic LB-PF equation [10,11] from the extended continuous kinetic equation. Firstly, a discrete distribution function is defined as Si (x, t) = wi S(x, ci , t), whose evolution is governed by the following equation ) ∂Si ∂Si Sϕ 1 ( (eq) an + ciβ Si − Si + wi , =− (12) ∂t ∂xβ τc τ0 where wi and ci are the weight coefficients and discrete velocities of particles moving along the ith direction ∑ ∑ (eq) in discrete physical space. It can be proved that ϕ = i Si , Sϕ = i wi Sϕ and hi = wi S (eq) . Integrating Eq. (12) in the time interval [t, t + δt] along the characteristic line yields an Si (x + ci δt, t + δt) = Si (x, t) + (an − 1)Si (x + ci δt, t) ( ) 1 Sϕ (eq) Si (x, t) − Si (x, t) + wi δt, − τp (x, t) τ0
(13)
where τp is the dimensionless relaxation time and τp = τc /δt. Obviously, the above equation fits into the framework of lattice BGK (LBGK) model which is widely applied in various areas. Among the available LBGK models, the DnQb (n-dimensional b-velocity) models [4] are used to determine wi and ci . In the D2Q9 model, the discrete velocities and weight coefficients are defined by [ ] 0, c, 0, −c, 0, c, −c, −c, c [ci , i = 0, . . . , 8] = , (14) 0, 0, c, 0, −c, c, c, −c, −c and w0 = 4/9, w1−4 = 1/9, w5−8 = 1/36. The equilibrium distribution function is given by ] [ vβ2 (ciβ vβ )2 ciβ vβ (eq) − 2 , Si (x, t) = wi ϕ + 2 + cs 2c4s 2cs
(15)
where v = −W02 N /τ0 and c = δx/δt. The equilibrium function satisfies the following conditions which can be verified ∑ (eq) ∑ ∑ W2 (eq) (eq) Si = ϕ, ciα Si = − 0 Nα , ciα ciβ Si = ϕc2s δαβ . (16) τ 0 i i i
D. Sun / Applied Mathematics Letters 103 (2020) 106222
4
where δαβ is the Kronecker delta function. The dimensionless macroscopic scalar ϕ is computed through ∑ ϕ = i Si . We set c2s = c2 /3 for the D2Q9 model, and identify (τp − 1/2)c2s δt = a2s (n)W02 /τ0 to recover the phase field equation. For the convective diffusion equation, Eq. (11), the evolution equation of discrete BGK scheme is Ti (x + ci δt, t + δt) = Ti (x, t) −
1 [Ti (x, t) − Tieq (x, t)] + wi Sθ δt, τθ
(17)
where Ti is the distribution function for θ, τθ is the dimensionless single relaxation time, and the last term represents the discretization of Sθ . Correspondingly, the equilibrium distribution function is given by ] [ u2β (ciβ uβ )2 ciβ uβ eq + − 2 . (18) Ti (x, t) = wi θ 1 + c2s 2c4s 2cs ∑ The dimensionless macroscopic temperature θ is calculated by θ = i Ti , and the relaxation time τθ can be computed via (τθ − 1/2)c2s δt = κh . It should be emphasized that the convective diffusion equation and discrete evolution equation for solute transfer can be derived from the extended kinetic equation in the same way as deriving equations for phase field and heat transport. Therefore, the discrete kinetic scheme can be used to model evolution of multi-component systems with and without anisotropic effects. 3. Numerical simulation Numerical simulations are performed to examine accuracy and efficiency of the present scheme. The exact formula of a2s (n) should be determined firstly for it enables the diversity of crystal patterns in liquid–solid phase transitions. In this work, dendritic growth from the cubic crystal system is used as an example. The anisotropy function for such system is given as as (n) = 1 − 3εs + 4εs n4α .
(19)
Here, the vector n is defined as directing from solid side to liquid side at the interface, and it is computed by n(x, t) = −∇ϕ/|∇ϕ|. Correspondingly, the vector N is computed by [ ]T ∂as (n) ∂as (n) ∂as (n) 2 N = |∇ϕ| as (n) , , . (20) ∂(∂x ϕ) ∂(∂y ϕ) ∂(∂z ϕ) To couple the phase field and heat transfer equations, we assign Sϕ and Sθ as Sϕ = [ϕ − λθ(1 − ϕ2 )](1 − ϕ2 ) and Sθ (x, t) = ∂ϕ/2∂t. Here, λ is the coupling coefficient. For simplicity, dendritic growth in pure diffusion condition is considered as a benchmark, and thus u = 0. In the simulations, a circular seed is initially placed in the domain center. The preferring growth orientation of dendrite is set to be parallel with the coordinate axes. Symmetrical conditions are applied on the side walls of the computational domain. The 1024 × 1024 meshes are used in all of the simulations. The spacial step, time interval and thermal diffusion coefficient are set as δx = 0.4W0 , δt = 0.008τ0 and κh = 4.0W02 /τ0 , respectively. By setting W0 = 1 and τ0 = 1, the coefficient λ is therefore equal to λ = a1 W0 /d0 = 6.382 where a1 = 0.8839 and d0 = 0.1385W0 . The anisotropic strength is chosen as ε = 0.05. Fig. 1 shows the time history of tip growth velocities at initial temperature θ0 = −0.30 and θ0 = −0.55. The simulated tip growth velocities agree exactly with the previous data in reference [10]. Fig. 2 depicts comparison of the simulated dimensionless tip velocity as a function of |θ0 | with analytical data by the Lipton–Glicksman–Kurz model [17,18]. It shows that the profiles of tip velocity versus |θ0 | simulated by the present model have the same trends with the analytical predictions. Furthermore, the simulations are performed on a small windows workstation (Intel Xeon CPU E3-1245v6 @ 3.70 GHz, 64GB memory) achieving the computational speed of 2.34 MLUPS (Million Lattice Updates Per Second). It demonstrates that the present model has super good algorithmic efficiency.
D. Sun / Applied Mathematics Letters 103 (2020) 106222
5
Fig. 1. Simulated time history of tip velocity vg for dendritic growth in the conditions of θ0 = −0.30 and −0.55.
Fig. 2. Comparison of the simulated dimensionless tip velocity against |θ0 | with analytical data by the Lipton–Glicksman–Kurz model.
4. Conclusion An extended kinetic equation is constructed for modeling of anisotropic liquid–solid phase transitions. On the fundamental principle of the kinetic theory, the equation can recover the macroscopic phase field and convective diffusion equations directly. A discrete kinetic scheme is proposed through discretization of the equation in both time and phase space. A model with the anisotropic streaming-relaxation operator is presented, and numerical simulations of anisotropic phase transitions are performed as examples. The results agree well with the analytical solutions, which demonstrates the present scheme is accurate and efficient. The present discrete kinetic scheme provides an alternative numerical approach to study phase transitions with anisotropic effects. CRediT authorship contribution statement Dongke Sun: Conceptualization, Investigation, Data curation, Validation.
D. Sun / Applied Mathematics Letters 103 (2020) 106222
6
Acknowledgments Thanks to my friends who inspired me to complete this work. References [1] Z. Guo, K. Xu, R. Wang, Discrete unified gas kinetic scheme for all Knudsen number flows: low-speed isothermal case, Phys. Rev. E 88 (2013) 033305. [2] X.-P. Luo, H.-L. Yi, A discrete unified gas kinetic scheme for phonon Boltzmann transport equation accounting for phonon dispersion and polarization, Int. J. Heat Mass Transfer 114 (2017) 970–980. [3] J. Li, C. Zhong, D. Pan, C. Zhuo, A gas-kinetic scheme coupled with SST model for turbulent flows, Comput. Math. Appl. 78 (4) (2019) 1227–1242. [4] Y.H. Qian, D. d’Humi` eres, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17 (1992) 479. [5] D. d’Humi` eres, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 360 (2002) 437–451. [6] A. Xu, G. Zhang, Y. Gan, F. Chen, X. Yu, Lattice Boltzmann modeling and simulation of compressible flows, Front. Phys. 7 (5) (2012) 582. [7] A. Xu, G. Zhang, Y. Zhang, Discrete Boltzmann modeling of compressible flows, in: G.Z. Kyzas, A.C. Mitropoulos (Eds.), Kinetic Theory, InTech, Rijeka, 2018, Ch. 02. [8] Z. Chai, T.S. Zhao, Lattice Boltzmann model for the convection–diffusion equation, Phys. Rev. E 87 (2013) 063309. [9] Z. Chai, D. Sun, H. Wang, B. Shi, A comparative study of local and nonlocal Allen-Cahn equations with mass conservation, Int. J. Heat Mass Transfer 122 (2018) 631–642. [10] A. Younsi, A. Cartalade, On anisotropy function in crystal growth simulations using lattice Boltzmann equation, J. Comput. Phys. 325 (2016) 1–21. [11] D. Sun, H. Xing, X. Dong, Y. Han, An anisotropic lattice Boltzmann - phase field scheme for numerical simulations of dendritic growth with melt convection, Int. J. Heat Mass Transfer 133 (2019) 1240–1250. [12] J.G. Zhou, P.M. Haygarth, P.J.A. Withers, C.J.A. Macleod, P.D. Falloon, K.J. Beven, M.C. Ockenden, K.J. Forber, M.J. Hollaway, R. Evans, A.L. Collins, K.M. Hiscock, C. Wearing, R. Kahana, M.L. Villamizar Velez, Lattice Boltzmann method for the fractional advection-diffusion equation, Phys. Rev. E 93 (2016) 043310. [13] R. Du, D. Sun, B. Shi, Z. Chai, Lattice Boltzmann model for time sub-diffusion equation in Caputo sense, Appl. Math. Comput. 358 (2019) 80–90. [14] Aiguo Xu, Guang-cai Zhang, Yudong Zhang, Pei Wang, Yangjun Ying, Discrete Boltzmann model for implosion- and explosion-related compressible flow with spherical symmetry, Front. Phys. 13 (5) (2018) 135102. [15] T. Takaki, R. Sato, R. Rojas, M. Ohno, Y. Shibuta, Phase-field lattice Boltzmann simulations of multiple dendrite growth with motion, collision, and coalescence and subsequent grain growth, Comput. Mater. Sci. 147 (2018) 124–131. [16] H. Xing, K. Ankit, X. Dong, H. Chen, K. Jin, Growth direction selection of tilted dendritic arrays in directional solidification over a wide range of pulling velocity: a phase-field study, Int. J. Heat Mass Transfer 117 (2018) 1107–1114. [17] D.K. Sun, M.F. Zhu, S.Y. Pan, D. Raabe, Lattice Boltzmann modeling of dendritic growth in a forced melt convection, Acta Mater. 57 (2009) 1755–1767. [18] J. Lipton, M. Glicksman, W. Kurz, Dendritic growth into undercooled alloy metals, Mater. Sci. Eng. 65 (1) (1984) 57–63.
Contents lists available at ScienceDirect
Applied Mathematics Letters www.elsevier.com/locate/aml
A discrete kinetic scheme to model anisotropic liquid–solid phase transitions Dongke Sun School of Mechanical Engineering, Southeast University, Nanjing 211189, China
article
info
Article history: Received 20 December 2019 Received in revised form 5 January 2020 Accepted 5 January 2020 Available online 13 January 2020 Keywords: Discrete kinetic scheme Lattice Boltzmann Phase field Phase transition
abstract An extended Boltzmann equation is constructed to model phase transitions with anisotropic effects. It is directly transformed into a discrete kinetic scheme in both time and phase space. The new scheme can be used to derive anisotropic lattice Boltzmann-phase field equation and convective diffusion equations. A lattice Boltzmann algorithm with the anisotropic streaming-relaxation operator is proposed, and numerical simulations of liquid–solid phase transition are performed as examples. The results agree well with previous numerical data and analytical solutions. It demonstrates that the scheme has acceptable numerical accuracy and computational efficiency, and the scheme can be used to study anisotropic liquid– solid phase transitions. This work provides an alternative numerical approach to model and simulate phase transitions with relatively fast computational efficiency. © 2020 Elsevier Ltd. All rights reserved.
1. Introduction Numerical methods based on the kinetic theory which bridges micro- and macro-physics have achieved great success both in academic and industrial areas [1–3]. Among these methods, the lattice Boltzmann (LB) method [4,5] and the discrete Boltzmann (DB) method [6,7] can be considered as a special finite difference scheme for discretization of the Boltzmann equation. For some distinguished advantages, these methods are being developed along two complementary directions. The LB scheme is widely used to solve nonlinear partial differential equations (PDEs), such as convective-diffusion equations [8], Allen–Cahn equation [9], phase field (PF) equations [10,11], and fractional PDEs [12,13]. The DB method is adopted to kinetically model non-equilibrium behaviors of fluids and reveal underlying physical features of various complex flows [14]. Phase transitions with anisotropic effects have been studied for decades by various phenomenological modeling and simulation techniques [15–17]. However, little research has been done on modeling of phase transitions by fitting the kinetic theoretical framework. Therefore, we attempt to construct a kinetic theory based equation to describe anisotropic phase transitions, and make a priori derivation from the kinetic equation to macroscopic phase field and convective diffusion equations. A discrete kinetic scheme is proposed E-mail address: [email protected]. https://doi.org/10.1016/j.aml.2020.106222 0893-9659/© 2020 Elsevier Ltd. All rights reserved.
D. Sun / Applied Mathematics Letters 103 (2020) 106222
2
to solve these equations, and the liquid–solid phase transition is selected as a benchmark to examine the present scheme. Finally, a brief conclusion closes the paper. 2. Extended kinetic equation and discrete kinetic scheme According to the kinetic theory, fluid flows can be described by the Boltzmann equation with the particle distribution function (PDF), f (x, ξ, t), ∂f + ξ · ∇f + a · ∇ξ f = Ω (f ), (1) ∂t where f is defined as the probability of finding a particle moving with velocity ξ at position x and time t, and a is the acceleration evoked by external forces. Ω (f ) denotes a collision operator, and it is usually approximated by some simplified models. One of the most popular choices is the Bhatnagar–Gross–Krook (BGK) approximation, which is given as ) 1 ( f − f (eq) . (2) Ω (f ) = − τc Here, τc represents the relaxation time, f (eq) is the equilibrium PDF. To extend the Boltzmann equation to model anisotropic phase transitions, the particle distribution function is redefined as the finding probability of a molecule moving with ξ at x and t. It is similar as the PDF defined in the Boltzmann equation for fluid flows, but the new PDF is used to describe phases rather than fluids. As phase transitions usually occur with anisotropic properties, an excess factor, an , is introduced to describe phase changes evoked by the anisotropic effect. Besides the collision of particles, the particle should satisfy mass conservation, no matter there exist phase changes or not. The continuous equation of the kinetic model with an is constructed. Its tensor format is written as ) ∂f 1 ( ∂f ∂f + ξβ + aβ =− f − f (eq) . an (3) ∂t ∂xβ ∂ξβ τc Here, an is a function of interface normal vector n, and we have an = an (n). Eq. (3) can be considered as an extension of Boltzmann equation which plays a central role in the kinetic theory. The particle velocity ξ can be split into two parts, the macroscopic velocity u and the relative velocity ξ, i.e. ξ = u + ζ. Here, u is independent of velocity space while ζ depends on velocity space. Multiplying equation (3) by ξα2 /2 and letting S = ξα2 f /2, S (eq) = ξα2 f (eq) /2 yield ) ∂S ∂S 1 ( ∂S + ξβ + aβ = S − S (eq) , (4) an ∂t ∂xβ ∂ξβ τc Integrating the equation over velocity space, result in ∫ ∫ ∫ ∫ ( ) ∂ ∂ ∂S 3 1 3 3 an Sd ξ + ξβ Sd ξ + aβ dξ = ξα2 S − S (eq) d3 ξ. ∂t ∂xβ ∂ξβ τc
(5)
As the collision conserves total energy, the term at the right hand side of the above equation is 0. Setting ∫ ψ = Sd3 ξ yields the macroscopic equation an
∂ψ ∂(uβ ψ) ∂qβ + =− + Sψ ∂t ∂xβ ∂xβ
with
Sψ =
∂(uα σαβ ) + ρaβ uβ . ∂xβ
(6)
Here, σαβ is the stress tensor given as σαβ and qβ is the flux calculated as,
∂ =− ∂xβ
1 qβ = 2
∫
∫
ζα ζβ f d3 ξ,
ζα ζα ζβ f d3 ξ.
(7)
(8)
D. Sun / Applied Mathematics Letters 103 (2020) 106222
3
When compressive and deformation work is neglected, we have ∂(uα σαβ )/∂xβ = 0 and then Sψ = ρaβ uβ . The flux qβ is usually related to the gradient of a scalar variable in a certain field. One choice is setting qβ = −ηp ∂ψ/∂xβ with ηp as the diffusion coefficient. Correspondingly, Eq. (6) can be written as ( ) ∂ψ ∂(uβ ψ) ∂ ∂ψ an + = ηp + Sψ . (9) ∂t ∂xβ ∂xβ ∂xβ Substitution of an = a2s (n), ψ = ϕ, uβ ψ = W02 Nβ /τ0 , ηp = W02 a2s (n)/τ0 , and Sψ = Sϕ /τ0 in Eq. (9) results in ( ) ∂ϕ ∂ϕ 2 ∂Nβ 2 ∂ 2 2 − W0 = W0 as (n) + Sϕ (10) τ0 as (n) ∂t ∂xβ ∂xβ ∂xβ where N and n are the anisotropic vector and the interfacial normal vector, respectively. Finally, a general formula of the phase field equation is derived from the extended kinetic equation, i.e. Eq. (3). The convective diffusion equation without anisotropic effect for heat transfer can also be derived in the same way. Setting an = 1, ψ = θ, ηd = κh , and Sψ = Sθ , we obtain ( ) ∂ ∂θ ∂(uβ θ) ∂θ + = κh + Sθ , (11) ∂t ∂xβ ∂xβ ∂xβ where θ denotes dimensionless temperature, and κh is diffusion coefficient for heat transfer. In the absence of flow, u = 0 and Eq. (11) is reduced to the pure diffusion equation with a source term. Now, we start to derive the anisotropic LB-PF equation [10,11] from the extended continuous kinetic equation. Firstly, a discrete distribution function is defined as Si (x, t) = wi S(x, ci , t), whose evolution is governed by the following equation ) ∂Si ∂Si Sϕ 1 ( (eq) an + ciβ Si − Si + wi , =− (12) ∂t ∂xβ τc τ0 where wi and ci are the weight coefficients and discrete velocities of particles moving along the ith direction ∑ ∑ (eq) in discrete physical space. It can be proved that ϕ = i Si , Sϕ = i wi Sϕ and hi = wi S (eq) . Integrating Eq. (12) in the time interval [t, t + δt] along the characteristic line yields an Si (x + ci δt, t + δt) = Si (x, t) + (an − 1)Si (x + ci δt, t) ( ) 1 Sϕ (eq) Si (x, t) − Si (x, t) + wi δt, − τp (x, t) τ0
(13)
where τp is the dimensionless relaxation time and τp = τc /δt. Obviously, the above equation fits into the framework of lattice BGK (LBGK) model which is widely applied in various areas. Among the available LBGK models, the DnQb (n-dimensional b-velocity) models [4] are used to determine wi and ci . In the D2Q9 model, the discrete velocities and weight coefficients are defined by [ ] 0, c, 0, −c, 0, c, −c, −c, c [ci , i = 0, . . . , 8] = , (14) 0, 0, c, 0, −c, c, c, −c, −c and w0 = 4/9, w1−4 = 1/9, w5−8 = 1/36. The equilibrium distribution function is given by ] [ vβ2 (ciβ vβ )2 ciβ vβ (eq) − 2 , Si (x, t) = wi ϕ + 2 + cs 2c4s 2cs
(15)
where v = −W02 N /τ0 and c = δx/δt. The equilibrium function satisfies the following conditions which can be verified ∑ (eq) ∑ ∑ W2 (eq) (eq) Si = ϕ, ciα Si = − 0 Nα , ciα ciβ Si = ϕc2s δαβ . (16) τ 0 i i i
D. Sun / Applied Mathematics Letters 103 (2020) 106222
4
where δαβ is the Kronecker delta function. The dimensionless macroscopic scalar ϕ is computed through ∑ ϕ = i Si . We set c2s = c2 /3 for the D2Q9 model, and identify (τp − 1/2)c2s δt = a2s (n)W02 /τ0 to recover the phase field equation. For the convective diffusion equation, Eq. (11), the evolution equation of discrete BGK scheme is Ti (x + ci δt, t + δt) = Ti (x, t) −
1 [Ti (x, t) − Tieq (x, t)] + wi Sθ δt, τθ
(17)
where Ti is the distribution function for θ, τθ is the dimensionless single relaxation time, and the last term represents the discretization of Sθ . Correspondingly, the equilibrium distribution function is given by ] [ u2β (ciβ uβ )2 ciβ uβ eq + − 2 . (18) Ti (x, t) = wi θ 1 + c2s 2c4s 2cs ∑ The dimensionless macroscopic temperature θ is calculated by θ = i Ti , and the relaxation time τθ can be computed via (τθ − 1/2)c2s δt = κh . It should be emphasized that the convective diffusion equation and discrete evolution equation for solute transfer can be derived from the extended kinetic equation in the same way as deriving equations for phase field and heat transport. Therefore, the discrete kinetic scheme can be used to model evolution of multi-component systems with and without anisotropic effects. 3. Numerical simulation Numerical simulations are performed to examine accuracy and efficiency of the present scheme. The exact formula of a2s (n) should be determined firstly for it enables the diversity of crystal patterns in liquid–solid phase transitions. In this work, dendritic growth from the cubic crystal system is used as an example. The anisotropy function for such system is given as as (n) = 1 − 3εs + 4εs n4α .
(19)
Here, the vector n is defined as directing from solid side to liquid side at the interface, and it is computed by n(x, t) = −∇ϕ/|∇ϕ|. Correspondingly, the vector N is computed by [ ]T ∂as (n) ∂as (n) ∂as (n) 2 N = |∇ϕ| as (n) , , . (20) ∂(∂x ϕ) ∂(∂y ϕ) ∂(∂z ϕ) To couple the phase field and heat transfer equations, we assign Sϕ and Sθ as Sϕ = [ϕ − λθ(1 − ϕ2 )](1 − ϕ2 ) and Sθ (x, t) = ∂ϕ/2∂t. Here, λ is the coupling coefficient. For simplicity, dendritic growth in pure diffusion condition is considered as a benchmark, and thus u = 0. In the simulations, a circular seed is initially placed in the domain center. The preferring growth orientation of dendrite is set to be parallel with the coordinate axes. Symmetrical conditions are applied on the side walls of the computational domain. The 1024 × 1024 meshes are used in all of the simulations. The spacial step, time interval and thermal diffusion coefficient are set as δx = 0.4W0 , δt = 0.008τ0 and κh = 4.0W02 /τ0 , respectively. By setting W0 = 1 and τ0 = 1, the coefficient λ is therefore equal to λ = a1 W0 /d0 = 6.382 where a1 = 0.8839 and d0 = 0.1385W0 . The anisotropic strength is chosen as ε = 0.05. Fig. 1 shows the time history of tip growth velocities at initial temperature θ0 = −0.30 and θ0 = −0.55. The simulated tip growth velocities agree exactly with the previous data in reference [10]. Fig. 2 depicts comparison of the simulated dimensionless tip velocity as a function of |θ0 | with analytical data by the Lipton–Glicksman–Kurz model [17,18]. It shows that the profiles of tip velocity versus |θ0 | simulated by the present model have the same trends with the analytical predictions. Furthermore, the simulations are performed on a small windows workstation (Intel Xeon CPU E3-1245v6 @ 3.70 GHz, 64GB memory) achieving the computational speed of 2.34 MLUPS (Million Lattice Updates Per Second). It demonstrates that the present model has super good algorithmic efficiency.
D. Sun / Applied Mathematics Letters 103 (2020) 106222
5
Fig. 1. Simulated time history of tip velocity vg for dendritic growth in the conditions of θ0 = −0.30 and −0.55.
Fig. 2. Comparison of the simulated dimensionless tip velocity against |θ0 | with analytical data by the Lipton–Glicksman–Kurz model.
4. Conclusion An extended kinetic equation is constructed for modeling of anisotropic liquid–solid phase transitions. On the fundamental principle of the kinetic theory, the equation can recover the macroscopic phase field and convective diffusion equations directly. A discrete kinetic scheme is proposed through discretization of the equation in both time and phase space. A model with the anisotropic streaming-relaxation operator is presented, and numerical simulations of anisotropic phase transitions are performed as examples. The results agree well with the analytical solutions, which demonstrates the present scheme is accurate and efficient. The present discrete kinetic scheme provides an alternative numerical approach to study phase transitions with anisotropic effects. CRediT authorship contribution statement Dongke Sun: Conceptualization, Investigation, Data curation, Validation.
D. Sun / Applied Mathematics Letters 103 (2020) 106222
6
Acknowledgments Thanks to my friends who inspired me to complete this work. References [1] Z. Guo, K. Xu, R. Wang, Discrete unified gas kinetic scheme for all Knudsen number flows: low-speed isothermal case, Phys. Rev. E 88 (2013) 033305. [2] X.-P. Luo, H.-L. Yi, A discrete unified gas kinetic scheme for phonon Boltzmann transport equation accounting for phonon dispersion and polarization, Int. J. Heat Mass Transfer 114 (2017) 970–980. [3] J. Li, C. Zhong, D. Pan, C. Zhuo, A gas-kinetic scheme coupled with SST model for turbulent flows, Comput. Math. Appl. 78 (4) (2019) 1227–1242. [4] Y.H. Qian, D. d’Humi` eres, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17 (1992) 479. [5] D. d’Humi` eres, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 360 (2002) 437–451. [6] A. Xu, G. Zhang, Y. Gan, F. Chen, X. Yu, Lattice Boltzmann modeling and simulation of compressible flows, Front. Phys. 7 (5) (2012) 582. [7] A. Xu, G. Zhang, Y. Zhang, Discrete Boltzmann modeling of compressible flows, in: G.Z. Kyzas, A.C. Mitropoulos (Eds.), Kinetic Theory, InTech, Rijeka, 2018, Ch. 02. [8] Z. Chai, T.S. Zhao, Lattice Boltzmann model for the convection–diffusion equation, Phys. Rev. E 87 (2013) 063309. [9] Z. Chai, D. Sun, H. Wang, B. Shi, A comparative study of local and nonlocal Allen-Cahn equations with mass conservation, Int. J. Heat Mass Transfer 122 (2018) 631–642. [10] A. Younsi, A. Cartalade, On anisotropy function in crystal growth simulations using lattice Boltzmann equation, J. Comput. Phys. 325 (2016) 1–21. [11] D. Sun, H. Xing, X. Dong, Y. Han, An anisotropic lattice Boltzmann - phase field scheme for numerical simulations of dendritic growth with melt convection, Int. J. Heat Mass Transfer 133 (2019) 1240–1250. [12] J.G. Zhou, P.M. Haygarth, P.J.A. Withers, C.J.A. Macleod, P.D. Falloon, K.J. Beven, M.C. Ockenden, K.J. Forber, M.J. Hollaway, R. Evans, A.L. Collins, K.M. Hiscock, C. Wearing, R. Kahana, M.L. Villamizar Velez, Lattice Boltzmann method for the fractional advection-diffusion equation, Phys. Rev. E 93 (2016) 043310. [13] R. Du, D. Sun, B. Shi, Z. Chai, Lattice Boltzmann model for time sub-diffusion equation in Caputo sense, Appl. Math. Comput. 358 (2019) 80–90. [14] Aiguo Xu, Guang-cai Zhang, Yudong Zhang, Pei Wang, Yangjun Ying, Discrete Boltzmann model for implosion- and explosion-related compressible flow with spherical symmetry, Front. Phys. 13 (5) (2018) 135102. [15] T. Takaki, R. Sato, R. Rojas, M. Ohno, Y. Shibuta, Phase-field lattice Boltzmann simulations of multiple dendrite growth with motion, collision, and coalescence and subsequent grain growth, Comput. Mater. Sci. 147 (2018) 124–131. [16] H. Xing, K. Ankit, X. Dong, H. Chen, K. Jin, Growth direction selection of tilted dendritic arrays in directional solidification over a wide range of pulling velocity: a phase-field study, Int. J. Heat Mass Transfer 117 (2018) 1107–1114. [17] D.K. Sun, M.F. Zhu, S.Y. Pan, D. Raabe, Lattice Boltzmann modeling of dendritic growth in a forced melt convection, Acta Mater. 57 (2009) 1755–1767. [18] J. Lipton, M. Glicksman, W. Kurz, Dendritic growth into undercooled alloy metals, Mater. Sci. Eng. 65 (1) (1984) 57–63.