Mass transfer analyses of reactive boundary schemes for lattice Boltzmann method with staircase approximation

Lattice Boltzmann (LB) methods with reactive boundary conditions are widely used in pore‑scale simulations of dissolution and ablation processes. The staircase approxi‑ mation of curved boundary is often employed because of its simplicity in handling solid structure changes. In this work, the mass transfer of two typical LB reactive boundary schemes are analyzed for the staircase boundary. The Type I boundary scheme is based on relations of local distribution functions and a wet‑node boundary mesh. The Type II boundary scheme adopts the half‑way bounce‑back scheme. Bound‑ ary concentrations are determined by finite difference, and a link‑wise boundary mesh is used. The analyses demonstrate that for straight boundaries, both the boundary schemes have accurate mass transfer rates, which means the mass transfer calculated by exchanges of distribution functions is the same as that calculated by reaction rates. For curved boundaries with staircase approximation, including interfacial normal direc‑ tions in the Type I boundary scheme can provide accurate mass transfer for inclined straight boundaries. However, if the staircase boundary geometry is used directly without normal directions, the reaction rate will be overestimated. One‑dimensional and two‑dimensional reaction‑diffusion processes with dissolution are simulated to validate the analyses. Both the boundary schemes work well for one‑dimensional simulations. For two‑dimensional simulations, the Type II boundary scheme signifi‑ cantly overestimates the reaction rate, and stronger artificial anisotropic effects are observed. The Type I boundary scheme with normal directions has better performance, but error still exists.


Introduction
Pore-scale studies of chemical reactions and porous structure evolutions are important in many research fields, such as mineral precipitation and dissolution, geologic carbon storage, and C/C composite ablation [1][2][3].Because of its convenience in dealing with complex boundary conditions, lattice Boltzmann (LB) methods have been widely used in numerical simulations of pore-scale heat and mass transfer with chemical reactions [4].By combining with a volume of pixel (VOP) representation of solid structures, the LB methods can also be used to simulate morphology evolutions of porous structures [5][6][7][8][9].The basic idea is to define a volume fraction of fluid in each node, which indicates whether the node is solid, fluid, or interface.The mass transfer in the fluid region is simulated by the LB model with reactive boundary conditions on the fluid-solid interface.The volume fractions of the interface nodes are changed based on the reaction rate.Then, the evolutions of solid structures can be achieved.
Proper treatments of the reactive boundary conditions on the interface are important for the accuracy of the above simulations.The boundary schemes proposed and employed in some of the existing LB researches are summarized in Table 1.For the mass transfer and structure evolution, there are two parts of the treatment of reactive boundary conditions.For the fluid side, reactive boundary schemes are needed for LB mass transfer models.For the solid side, the reaction rates should be calculated to update the volume fraction.
For the reactive boundary schemes of LB models, they can be summarized into three categories although there are differences in detail [14]. 1) The first category of schemes is to take advantage of the moments of distribution functions in LB models [5][6][7][8][9][10][11][12][13][14][15][16].Because the first-order moments in LB mass transfer models are related to concentration gradients, the reactive boundary condition can be rearranged into algebraic equations of the distribution functions.The unknown concentration and distribution functions can be solved locally.2) The second category of schemes is to calculate the unknown concentration on the solid boundary by finite difference of the concentration gradient, which is related to the reaction rate [17][18][19].Then, the boundary becomes a Dirichlet boundary and the bounce-back scheme can be employed.3) In the third category, the reaction rate is treated as a mass source and is added to the evolution equation [20,21].The first and second schemes have been widely used and validated, while the third scheme is less frequently used.It should be mentioned that the first and second boundary condition schemes can be combined for the Robin boundary condition.For example, Chen et al. [6] used finite difference to calculate the boundary concentration, and used the local moments of distribution functions to calculate the unknown distribution function.
Typically, the Cartesian mesh is used for the LB models.For straight boundaries, there are two types of boundary locations: 1) wet-node boundary located on lattice nodes; 2) link-wise boundary located on lattice links [22].The link-wise boundaries are often placed on the midway of the links, so that the half-way bounce-back boundary scheme can be easily employed [23,24].For curved boundaries, if the intersections between the boundaries and lattice links are known, the interpolations or extrapolations can be used to calculate the unknown distribution functions on the boundaries [25][26][27].The reactive boundary schemes are also extended to the curved boundary condition [12,13,19].3) Mass source in evolution equation Reaction rate calculation [20,21] However, for the simulations with solid structure evolution, only solid volume fractions on the boundary are known and the exact locations of the boundary are not explicitly given.A practical approach is to use the staircase approximation of the curved boundary, and to assume the boundary locates on the grid nodes or the midway of the links.
As for the fluid-solid boundary condition, the second part is the calculation of the mass transfer between the solid and fluid nodes to update the solid volume fractions.Because the equation of reaction rate is specified on the interfaces as the boundary condition, the reaction rates can be calculated and the solid volume fractions can be updated accordingly.This procedure has been employed in many works [5-9, 14-16, 20, 21], as shown in Table 1.However, an important parameter in calculating the reaction rate is the effective surface area for an interface node, because the reaction rate is proportional to it.Although it exists in many papers, the method to determine the effective surface area is not often provided.For example, Ju et al. [14] calculated it by the magnitude of the volume fraction gradient.Recently, Kashani et al. [7] and Izadi et al. [8] used the interface reconstruction techniques from the volume of fluid (VOF) method to determine the specific surface area geometrically.They showed that the VOP method can overestimate the reaction rate by more than 20% if the boundary was not properly treated [7].When structure evolutions are studied, the reconstruction of surface geometry can be complex.Therefore, to use a staircase approximation of the boundary is still appealing.
Meanwhile, if the mass transfer is observed from the fluid side, the exchange of distribution functions between fluid and solid nodes represents the mass transfer in LB simulations.Thus, the mass transfer can also be calculated from distribution functions.This treatment is consistent with the nature of LB models, and the total mass in both the fluid and solid regions is conserved.It is also easy to be implemented with the staircase approximation of the boundary, because the exact boundary geometry and specific surface area are not needed.Although the exchange of distribution functions is widely used for the calculations of momentum transfer [24], its application in mass transfer still needs further study.
Therefore, the present research aims to enrich the understanding of the LB reactive boundary schemes with staircase approximation.The exchange of distribution functions is used to determine the mass transfer of the boundary schemes.Its relations with the geometry of the boundary and the reaction rates are also analyzed.In the rest of the paper, the LB model for mass diffusion and two types of boundary schemes are described in Section 2. The analyses of the mass transfer of the two boundary schemes are conducted in Section 3. Finally, some numerical examples are used to further validate the analyses in Section 4. Conclusions are drawn in Section 5.

Lattice Boltzmann model for mass diffusion
Before analyzing the reactive boundary schemes, the LB model for mass or concentration transfer is briefly introduced.For simplicity, only diffusion of concentration is considered in this work, and the flow is not included.The D2Q5 model is chosen to be analyzed because it is simple and sufficient to model convection and diffusion processes of passive scalars [4,22].Compared with the D2Q9 model, it consumes less computational resources and can be more accurate and robust when the convection is not very strong [28].In addition, it prevents the complexity in dealing with the unknown diagonal velocities e 5-8 on boundaries.
The five discrete velocities of the D2Q5 model are e 0 = (0, 0), e 1-2 = c(±1, 0) and e 3-4 = c(0, ±1).The c = δx/δt is the lattice speed and δx and δt are the spatial and temporal steps.The evolution equation for the distribution function g i (x, t) at position x and time t is: where C is the concentration and τ is the dimensionless relaxation time.The equilibrium distribution g eq i (C) has a simple form: where w is the weight for the model.When w = 1/4 is used, the model becomes D2Q4.The concentration C is given by: The evolution equation ( 1) can be divided into collision and streaming steps: where the g * represents the distribution functions after the collision step and before the streaming step.By the Chapman-Enskog expansion analysis, it can be shown that the above model solves the concentration diffusion equation: where the diffusion coefficient D is related to the relaxation time τ by: Meanwhile, according to the Chapman-Enskog expansion, the following relations that are frequently used in the boundary schemes can be derived [10,14]: (1) where the non-equilibrium part of the distribution function is defined as: and the index i represents the direction opposite to i.The approximately equal symbol is used in Eqs. ( 8) and ( 9) because of the high-order terms in the expansion of g i .The highorder terms are usually ignored and the equal symbol can be used [10,[12][13][14].

Reactive boundary condition
As mentioned in the Introduction section, because to obtain the exact boundary position is difficult, the staircase approximation of boundaries is often used in the LB simulations with reactions and structure evolutions.There are two typical boundary positions, which are the link-wise boundary and the wet-node boundary, as shown in Fig. 1 [22].For the link-wise boundary, the solid-fluid boundary locates on the midway between the solid and fluid nodes, and the solid node locates in the center of a solid volume.For the wet-node boundary, the solid node locates on the boundary of the solid volume.The distribution functions on the solid nodes do not need to be defined for the link-wise boundary, but are necessary for the wet-node boundary.Usually, the boundary condition with local distribution function calculation is based on the wetnode boundary, which is denoted as Type I boundary in this work.The boundary condition with finite difference and bounce-back is based on the link-wise boundary, and is named Type II boundary in this work.The link-wise boundary is more compatible with the VOP method, because the solid volume fraction can be defined on the solid node.In Fig. 1, the g 3 on x w for Type I boundary and the g 3 on x f for Type II boundary are the unknown distribution functions that should be determined.
As for the reactive boundary, the following Robin boundary condition for linear heterogeneous reaction is considered [13]: where the normal vector n is pointing into the fluid region.The nonlinear reaction can be treated by similar procedures [29].For the one-dimensional example in Fig. 1, the equation becomes: (10)  For the first type of boundary schemes, we can assume that the g * 1 (x f , t − δt) moves to g 1 (x w , t) in a time step, which is denoted by g 1 for convenience.Then, Eqs. ( 8) and ( 11) are applied on the wall node x w : Substituting Eqs. ( 14) and ( 15) into Eq.( 13), the unknown g 3 (x w , t) can be solved as: where δx = cδt is used.
It should be mentioned that the g 3 (x w , t) cannot move back to g 3 (x f , t + δt) directly, because a collision step on x w is necessary before streaming.The after collision g * 3 is given by: Using Eq. ( 15) to replace the C w , we finally have: where the g 3 is calculated by Eq. ( 16).
This collision procedure is important for the mass conservation of the boundary scheme, which means the mass transfer calculated from distribution functions is the same as that calculated by the reaction rate.This procedure is seldom explicitly mentioned in the existing literature.The reason may be that for the wet-node boundary as shown in Fig. 1a, the wall node is treated as a fluid node, so the collision step is naturally included in the computation procedure.Details will be given in the following section.

b) Type II: Dirichlet boundary by finite difference
In the second type of boundary schemes, the finite difference is used to calculate the concentration on the boundary [18].The first-order finite difference is: (13) Substituting Eq. ( 19) into Eq.( 13), the concentration C w on the boundary can be solved as [18]: The second-order finite difference can be also employed by using the concentration on x f and x ff : The C w is then given by: Then, the half-way bounce-back scheme is used: It can be seen that if the first-order finite difference is used, this boundary scheme is still localized because only information on x f is needed.

Mass transfer for a straight boundary
In this section, the mass transfer of the boundary scheme for a straight boundary is firstly analyzed.Only one distribution function on the boundary is unknown.The results can be directly extended to the two-dimensional or three-dimensional analyses.Still using Fig. 1 as an example, the a 1 in Eq. ( 13) is assumed to be the diffusion coefficient D, and the boundary condition becomes: As mentioned in the Introduction, there are two ways to calculate the mass transfer from the solid to fluid region in a time step δt.The mass transfer can be calculated by the mass flux or reaction rate on the boundary, which is given by: On the other hand, the streaming of distribution functions also represents a mass exchange between solid and fluid nodes.For the Type I boundary scheme in Fig. 1a, if we assume there is an interface between x f and x w , the fluid region loses a mass (19 g 1 (x w , t) in the interval t − δt to t, and will receive a mass g * 3 (x w , t) in the interval t to t + δt.Thus, the mass transfer in a time step can be calculated by: The cδt and δx are included here because the unit of g is the same as the concentration, so a multiplication of volume cδtδx transforms it into mass.
As for the mass transfer in a time step for the Type II boundary scheme, it can be calculated by: If the mass transfer is consistent with the reactive boundary condition, the mass transfer calculated by Eqs. ( 26) and ( 27) should be the same as Eq. ( 25).Then, if the gain or loss of solid mass is calculated by the reaction rate, the total mass in all the fluid and solid regions will be conserved.

a) Type I: Local distribution function calculation
For the Type I boundary scheme, by substituting Eq. ( 18) into Eq.( 26), it can be found that: According to Eqs. ( 14) and (7), it can be derived that Eqs. ( 29) and (25) demonstrate that the Type I boundary scheme guarantees the conservation of mass.This analysis also indicates the importance of the collision step in Eq. (17).If the g 3 in Eq. ( 16) is directly used without the collision, Eq. ( 14) will provide a mistaken result as: A ratio (1-1/2τ) is missing.This inconsistence is also mentioned by Kang et al. [10], so they chose to use a different definition of diffusion coefficient on the boundary.The analyses demonstrate that considering the collision step can solve this inconsistence.

b) Type II: Dirichlet boundary by finite difference
For the mass transfer of the Type II boundary scheme, the following approximation of distribution function is used: ( (31) g i = wC − wτ δte i • ∇C, which can also be derived from Eqs. ( 14) and (15).It should be mentioned that Eq. ( 31) is for distribution functions after streaming and before collision.Thus, the g * 1 (x f , t) can be expressed as: Using Eq. ( 23), the mass transfer Eq. ( 27) is: If the first-order finite difference in Eq. ( 19) is used, Eq. ( 33) can be changed into: If we further assume that the concentration gradients on x w and x f are the same, Eq. ( 34) becomes: This additional assumption indicates a linear concentration distribution near the boundary, and it is also compatible with the first-order finite difference used on the boundary.Then, the conservation of mass transfer is also proved.
For transient simulations, the gradients on x w and x f can be different.Considering a simple example, if the initial concentration distribution is uniform, we have ∂C/∂x| f = 0.The mass transfer in the initial steps can have errors and will not be equal to the reaction rate.However, it can be expected that the errors will become smaller when the gradients near the boundary are developed during simulations.
For the second-order finite difference by Eq. ( 21), Eq. ( 33) can be changed into: which is different from Eq. ( 35).If a linear concentration distribution between x ff and x w can be assumed, Eq. ( 36) can be reduced to Eq. (35).Eq. ( 35) proves that the Type II boundary scheme with finite difference of concentration can also have correct mass transfer that equals to the reaction rate, although more assumptions are needed.

Comparison of the boundary schemes
Both the Type I and Type II boundary schemes use the same relations between opposite distribution functions, as given by Eqs. ( 15) and (23).However, the collision step is only employed for the Type I boundary scheme, and is not necessary for the Type II boundary scheme.Since both the boundary schemes have the same mass transfer rate, the relation between the boundary schemes is briefly discussed. (32) Here, Eq. ( 32) is still used to express the incoming g * 1 (x f , t − δt) .Then, both Eqs. ( 14) and ( 15) can be used to calculate the unknown g 3 , and the results should be equal.Thus, we have: which means: It can be seen that the finite difference of C w and C f in the Type I boundary scheme is a weighted average of the gradients at x f and x w , which is different from the finite difference in Eq. ( 19).If Eq. ( 39) is substituted into Eq.( 33) for the mass transfer, Eq. ( 30) can be obtained, which is consistent with the analysis in Section 3.1.
Therefore, the different treatments of finite difference near the solid boundary can be regarded as the differences between the two types of boundary schemes.Still, both the schemes can have mass transfer equal to the reaction rate.

Mass transfer for staircase approximation of curved boundary
The above analyses are based on straight boundaries.However, in the simulations of dissolution or ablation processes, curved boundaries appear and are often approximated by staircase geometries in LB simulations.A typical stair of two-dimensional boundary configuration is shown in Fig. 2, where the g 2 and g 3 are unknown distribution functions that should be determined.There is also a unit normal vector n on the solid node, which can be calculated by the distribution of the solid volume fraction.For the Type II scheme, the wall nodes on δx and δy surfaces are denoted by x wy and x wx , and the corresponding fluid nodes are x fy and x fx .For the Type I scheme, the solid node locates on the corner of the solid region.Then, different methods can be used to treat the reaction and mass transfer on this staircase boundary. (37) Fig. 2 A typical "stair" of the curved boundary Firstly, a rough approximation is to assume that the staircase boundary is the exact boundary geometry and the normal direction is parallel to the x or y direction.For both the Type I and Type II boundary schemes in Section 2.2, it can be found that the distribution functions on different directions are decoupled.Only distribution functions along the direction perpendicular to the boundary (x or y direction) are needed in the boundary schemes.Therefore, taking the Type II boundary scheme as the example, the two surfaces δx and δy can be treated individually.The concentrations on x wy and x wx can be calculated, and the unknown g 2 and g 3 are determined accordingly.For the Type I boundary, we can also assume that there are two concentrations on x w corresponding to each direction.According to the analyses in Section 3.1, the mass transfer on this staircase boundary is: Obviously, this treatment overestimates the reaction rate because the surface area of reaction for this cell is 2δx.
Secondly, there is another treatment which considers the normal direction [13,14].For the Type I boundary scheme, Ju et al. [14] suggested to use the following set of equations: where the normal vector is denoted by n = (n x , n y ).Eqs. ( 41)-( 44) have four equations for four unknowns: C w , n ⋅ ∇C w , g 2 and g 3 .Therefore, the unknown distribution functions can be solved.This boundary scheme is reduced to the one-dimensional Type I boundary scheme in Section 2.2 if the normal direction is along the x or y direction.
The analyses in Section 3.1 can be still used to study the mass transfer for this treatment.Because Eqs. ( 15), (17), and ( 18) are still valid for the staircase boundary, the mass transfer by streaming of distribution functions can be also calculated by Eq. (28).Considering both the x and y directions, we have the following relation for Fig. 2: In order to understand the mass transfer rate given by Eq. (45), Eq. ( 43) is rewritten as: (41) (46)

Thus, we have
If the surface geometry is the same as Fig. 3a, the normal vector becomes n = (−1, 1)/ √ 2 and the specific surface area is δa = √ 2δx .Using Eqs. ( 45) and (47), it can be proved that: Thus, for this special interface geometry, the mass transfer by the staircase approximation is the same as that calculated by the reaction rate.A more general result for the inclined straight boundary will be proved in Section 3.4, and it can be found that Eq. ( 48) is a simple case of the inclined straight boundary.
Generally, the surface geometry will be different from Fig. 3a during the structure evolution, as shown in Fig. 3b.The surface can be reconstructed by the methods in VOF [7,8].Under this general condition, the mass transfer calculated by Eq. ( 45) can be different from the reaction rate: However, if Eq. ( 49) is used to update the solid mass and volume fraction, the loss/gain of solid mass will be different from the mass transferred to/from the fluid region in LB models, which is described by Eq. (45).The overall mass conservation will be violated.

Special case of an inclined straight line
Ju et al. [14] conducted LB simulations of convection-diffusion processes in channels with inclined straight boundaries, and the Type I boundary scheme is used.The inclined angles in their study are tan(θ) = 1, 1/4, 1/2, 3/4, 4/3 and 5/3.The numerical results coincide well with the analytical results.According to the analyses in Section 3.3, we may provide an explanation for this accuracy for inclined straight boundaries. (47) Fig. 3 A sketch of the reaction area fine enough.However, for more general curved boundaries, it is hard to use straight lines to represent the curve with a finite size of mesh, and the error will still exist.

Numerical examples
In this section, three numerical examples will be used to verify the previous analyses.The D2Q4 LB model is used, and the weight in Eq. ( 2) is w = 1/4.

One-dimensional diffusion
A simple one-dimensional diffusion problem is firstly simulated to validate the boundary schemes in Section 2.2.The length of the region is 20δx.The diffusion coefficient is D = 0.1 and τ = 0.7.The concentration on the left boundary is C L = 0, while the reactive boundary on the right is: where C 0 = 2.0.If h is specified, the analytical result can be easily obtained.
For the Type I boundary scheme, 21 nodes are used and the Type I boundary scheme is employed on both the boundaries.For the concentration boundary on the left, Eq. ( 15) can be directly used to calculate the unknown g 1 .For the Type II boundary scheme, because there are two half grids on the left and right boundaries, only 19 nodes are used and the concentration on the right boundary is calculated by Eq. (20).
The results are compared with the analytical results in Fig. 5 for different h.It can be seen that both the Type I and Type II boundary schemes can provide correct boundary treatments for the simulations, and the numerical results coincide well with the analytical results.

One-dimensional reaction-diffusion
In this section, one-dimensional dissolution problems are simulated to validate the boundary conditions.The sketch of the problem is shown in Fig. 6.The left boundary (54) where r is a coefficient for the reaction rate.
The VOP method is used to treat the evolution of the solid position.A mass M(x s , 0) = M 0 is initially given in each solid node.For the solid node on the interface, mass losses calculated by Eq. ( 26) or Eq. ( 27) are subtracted from the solid node in each time step.If the mass M(x s , t) is less than zero, the solid node is changed into a fluid node.The distribution functions should also be initialized for this new fluid node.
Here the concentration on the new fluid node is extrapolated from the neighboring fluid nodes, and the non-equilibrium parts of the distribution functions are assumed to be the same as those of the neighboring fluid node.
The total length of the region is 100δx and the initial fluid region is 20δx.The initial concentration in the fluid region is C 0 .For the reaction Eq.(55), C 0 = 1.0,D = 0.01, r = 0.01 and M 0 = 2.0 are used.The reaction rate is constant, so the mass loss is linear.
The evolutions of the loss of solid mass are given in Fig. 7.The simulation results with the two boundary schemes and the linear analytical result are compared.It can be found that the solid mass loss in the simulation with the Type I boundary scheme is exactly the same as the analytical result.The Type II boundary scheme slightly overpredicts the reaction rate and the mass loss, although the difference is not significant.The reason is explained in Section 3.1.Because a uniform concentration is used as the initial condition, the Type II boundary scheme will overpredict the mass transfer in the beginning steps.The concentration distributions after 10000δt for the simulations are also compared in Fig. 8.It can be seen that the results of the two boundary schemes coincide well.
For the reaction Eq.(56), C 0 = 1.0,D = 0.1, r = 0.0005 and M 0 = 2.0 are used.If the diffusion is much faster than the reaction, we can approximate that the concentration in the fluid region reaches a quasi-steady linear distribution [30].Then, an approximate solution can be obtained, which is expressed as: where l is the length of the fluid region and l 0 = 20 is its initial value.The simulation results of l for two boundary schemes are compared with Eq. (57) in Fig. 9.It can be seen that the results coincide well with the analytical relation.Thus, the boundary schemes are validated.

Two-dimensional reaction-diffusion
In this section, the dissolution process of a two-dimensional circular disk is simulated to demonstrate the discussions in Section 3.3.The computational domain is 400δx × 400δx, and a disk with a radius close to 50δx is placed in the center of the domain.The method given by Thies [31] is used to calculate the normal vectors, which is based on the solid volume fractions.This method has been successfully employed in the free-surface LB model [32].To be consistent with the normal vector calculation, a previously developed free-surface LB model is used to obtain a smooth solid volume fraction distribution as the initial condition [33], which is shown in Fig. 10 and the normal vectors are also shown.For this problem, the constant reaction Eq. ( 55) is used and the parameters are C 0 = 1.0,D = 0.1, r = 0.0001 and M 0 = 2.0.The C 0 is the concentration on the boundaries of the domain.The governing equation for the radius R(t) is given by: The solution of Eq. ( 58) demonstrates that R should decrease linearly with time: For the Type I boundary scheme, Eqs. ( 41)-( 44) are used for the stair nodes.For the Type II boundary scheme, it is assumed that the stair nodes have two δx surfaces, as described by Eq. ( 40).The effective radius R is calculated based on the total mass of solid M according to M = πR 2 M 0 .The changes of the radii during simulations are shown in Fig. 11.It can be seen that the Type II boundary scheme with staircase approximation of solid boundary significantly overestimates the reaction rate.The Type I boundary scheme gives a better result compared with the Type II scheme.However, an overestimation of the reaction can be still observed.
The shapes of the disk and the concentration distribution after 2 × 10 5 δt are shown in Fig. 12.In addition, the shapes of the disk when the R is 40 are compared in Fig. 13.It can be seen that the dissolution process is not isotropic and the shape of the disk becomes close to a square.This anisotropic effect is weak for the Type I boundary scheme with normal directions, but it is more obvious for the Type II boundary scheme.The half-way bounce-back scheme tends to overestimate the reactions in 45°, 135°, 225° and 315° directions, where the solid boundary is rough.This phenomenon is consistent with the observation of artificial nucleation by Yu et al. [34], where the surface heat transfer is stronger in the aforementioned directions and anisotropic nucleation occurs on the cylinder surface.(59) R = R 0 − r M 0 t.
Fig. 11 The evolutions of radii in simulations with two boundary conditions

Conclusions
In the above sections, the mass transfer of two boundary schemes of the LB method for the reactive boundary conditions are discussed.It is proved that for straight boundaries, both the two boundary schemes have the correct mass transfer rate, which means the mass transfer calculated by the exchange of distribution functions is the same as that calculated by the reaction rate.It should be mentioned that the collision procedure is essential for the wet-node boundary of the Type I scheme, and the half-way bounce-back scheme can be used directly for the Type II scheme.
For curved boundaries that are approximated by the staircase mesh, using the normal directions in the Type I boundary scheme can provide accurate mass transfer for inclined straight boundaries with rational value of tan(θ).Otherwise, if the staircase geometry is assumed as the real geometry, the mass transfer or reaction rate will be obviously overestimated.Artificial anisotropic effects can be also observed.
It should be mentioned that although the analyses and numerical examples are onedimensional and two-dimensional, the results of the present work are also valid for three-dimensional D3Q7 or D3Q6 models.The reason is that the analyses of mass transfer for one direction only include the two distribution functions along that direction, so Fig. 12 The shapes of the disk and the concentration distributions after 2 × 10 5 δt Fig. 13 The shapes of the disk when the effective radius is 40 the different directions are decoupled.A pair of distribution functions for the z-direction can be added and the analyses in this work can be repeated.
Finally, if the geometry of the interface can be reconstructed during simulations, more accurate boundary locations can be obtained.The interpolation or extrapolation curved boundary schemes of the LB model can be used.However, further studies are still needed to understand the mass transfer of these boundary schemes and its relations with the reaction rate based on the interface geometry.

Fig. 1
Fig. 1 Two types of straight boundaries.The solid and hollow points represent the solid and fluid nodes, respectively.The shadow regions are the solid regions

Fig. 5
Fig. 5 Simulations of one-dimensional diffusion with different boundary schemes

Fig. 7 Fig. 8
Fig.7The solid mass losses for simulations with two boundary schemes and the analytical result for reaction Eq.(55)

Fig. 9 Fig. 10
Fig.9 The length of fluid region for simulations with two boundary schemes and the analytical result for reaction Eq.(56) R 2 = −2π rR.

Table 1
Reactive boundary schemes used in some of the existing publications