Extended Smoothed Boundary Method for Solving Partial Differential Equations with General Boundary Conditions on Complex Boundaries

# Extended Smoothed Boundary Method for Solving Partial Differential Equations with General Boundary Conditions on Complex Boundaries

Hui-Chia Yu Hsun-Yi Chen K. Thornton Department of Materials Science and Engineering, the University of Michigan, Ann Arbor, USA
###### Abstract

In this article, we describe an approach for solving partial differential equations with general boundary conditions imposed on arbitrarily shaped boundaries. A continuous function, the domain parameter, is used to modify the original differential equations such that the equations are solved in the region where a domain parameter takes a specified value while boundary conditions are imposed on the region where the value of the domain parameter varies smoothly across a short distance. The mathematical derivations are straightforward and generically applicable to a wide variety of partial differential equations. To demonstrate the general applicability of the approach, we provide four examples herein: (1) the diffusion equation with both Neumann and Dirichlet boundary conditions; (2) the diffusion equation with both surface diffusion and reaction; (3) the mechanical equilibrium equation; and (4) the equation for phase transformation with the presence of additional boundaries. The solutions for several of these cases are validated against corresponding analytical and semi-analytical solutions. The potential of the approach is demonstrated with five applications: surface-reaction-diffusion kinetics with a complex geometry, Kirkendall-effect-induced deformation, thermal stress in a complex geometry, phase transformations affected by substrate surfaces, and a self-propelled droplet.

###### keywords:
smoothed boundary method; diffuse interface method; complex microstructure; image-based simulation

, , and K. Thornton

## 1 Introduction

The smoothed boundary method [1, 2, 3] has recently been demonstrated as a powerful tool for solving diffusion equations with no-flux boundary conditions imposed at irregular boundaries within the computational domain. The method’s origin can be traced to the embedded boundary method and the immersed interface (boundary) method, both of which embed a more complicated domain in a computational box with simpler geometry. These methods are advantageous because they eliminate the need for a structural mesh when solving partial differential equations within the embedded geometries because the grid system is obtained by a discretization of the regular computational box. (For an overview, see Refs. [4, 5, 6, 7, 8, 9].) To impose boundary conditions at the immersed interfaces, a discretized Dirac delta function is employed to distribute a singular source over nearby grid points. Various studies have examined optimal discretization of the Dirac delta function [7, 10, 11]. Similarly, the level set method can also be considered an immersed-interface-type method because the boundary defined by the contour of the zero level set is embedded within a regular computational box. Although the level set method was developed mainly for tracking moving boundaries [12], it is also applicable for solving partial differential equations with boundary conditions imposed at the zero-level-set contour using a technique similar to the immersed interface method [13, 14]. In addition to the methods above, the phase field approach possesses certain similarities to embedding interfaces within the computational box and also has the significant advantage of avoiding the need to explicitly track the interfaces. However, phase field methods are not widely employed in simulations that involve explicit boundary conditions along interfaces. While the Gibbs-Thomson boundary condition is automatically imposed in the standard phase field model, there have been only a few studies in which phase field models were used with explicit boundary conditions at interfaces. For example, in solidification problems, equilibrium conditions, such as equilibrium temperature or concentration [15, 16], are imposed at solid-liquid interfaces in which the order parameter field and the temperature field are coupled via a latent heat term. Except for this type of phase field model, the direct application of boundary conditions at interfaces is rarely used because the construction of boundary conditions requires the tedious process of formally including an additional energy term in the energy functional, as suggested by Cahn [17]. Examples of imposed boundary conditions in the phase field model using modified energy functionals can be found in the recent works of Warren et al. [18, 19].

In contrast to the techniques for distributing a singular source of boundary conditions to grid points near the interfaces in the immersed interface method, the smoothed boundary method spreads the zero-thickness boundary into a finite-thickness diffuse interface using a phase-field-like, continuously transitioning domain indicator function (hereinafter termed “the domain parameter”). Mathematically, this method approximates a Heaviside step function as a hyperbolic tangent function having one specified uniform value in a domain and continuously changing its value across the interface to another value specifying the other domain. Therefore, boundary conditions are straightforwardly distributed among the grid points residing within the interfacial regions in which the domain parameter varies smoothly across a short distance. This method has been successfully employed in simulating diffusion processes [20, 21] and wave propagation [1, 2, 3, 22, 23] constrained within geometries described by a domain parameter with a no-flux boundary condition imposed on the diffuse interfaces. Similar approaches have also been proposed to solve differential equations constrained in domains defined by order parameters in the phase field model [24, 25, 26]. These works demonstrated the potential for this type of numerical method that circumvents the difficulties associated with constructing a finite element mesh (e.g., meshing the surface and then building a volumetric mesh based on the surface mesh or combining regular subdomains that can be easily meshed). Such an approach is particularly useful when complex structures are involved. However, the method was only applicable to no-flux boundary conditions, and no further extensions to other types of equations or boundary conditions have been reported. Recently, Lowengrub and coworkers [27, 28, 29, 30, 31, 32, 33] developed an alternative formulation for solving partial differential equations with various boundary conditions, based on asymptotic analyses commonly conducted in phase field modeling, which is different from the general derivation of the smoothed boundary method presented in this paper. Although such an implementation for imposing boundary conditions differs from the ‘formal’ practice suggested by Cahn [17], it dramatically simplifies the formulation, provides a justification of the method, and increases the applicability of the approach.

In this study, we provide a mathematically consistent smoothed boundary method and a precise derivation for the equations, such that the method is generalized from its limited original application to a wide range of differential equations and boundary conditions. We consider the following specific equations: (1) the diffusion equation with Neumann and/or Dirichlet boundary conditions; (2) the bulk diffusion equation coupled with surface diffusion and reaction; (3) the mechanical equilibrium equation for linear elasticity; and (4) the Allen-Cahn or Cahn-Hilliard equations with contact angles as boundary conditions. The method is especially useful for three-dimensional image-based simulations because of its efficiency and flexibility in handling complex geometries without structural-mesh techniques.

## 2 Background

The method is based on a diffuse interface description of different phases, similar to the continuously transitioning order parameters in the phase field method [34, 35, 36, 37, 38, 39] often employed in simulating phase transformations and microstructural evolutions in materials. Phase field models are based on thermodynamics and kinetics of multiphase system, in which phases (e.g., liquid, solid, vapor, or two solids or liquids with different compositions) are described by one or more order parameters with prescribed bulk values for each phase. At the interface, the order parameter changes in a controlled manner. Asymptotic analyses [39] can be used to show that the phase-field governing equations approach the corresponding sharp interface equations in the sharp interface limit.

Despite the advantages of phase-field-type diffuse interface methods for front tracking problems, we focus on another important advantage for efficiently solving differential equations within diffuse-interface-defined domains. Here, we adopt the concept to describe internal domain boundaries with an order-parameter-like domain parameter, which may or may not be stationary and takes a value of 1 inside the domain of interest and 0 outside. The equations are solved where the domain parameter is 1, with boundary conditions imposed where the domain parameter is at an intermediate value (approximately 0.5). Figure 1 schematically illustrates the sharp and diffuse interfaces. In the conventional sharp interface description, the domain of interest is and is bounded by a zero-thickness boundary, denoted ; see Fig. 1(a). Within , the partial differential equations are solved according to the boundary conditions imposed at . Conversely, in the diffuse interface description, we employ a continuous domain parameter, which is uniformly 1 within the domain of interest and uniformly 0 outside. In this case, the originally sharp domain boundary is smoothed to yield a diffuse interface with a finite thickness given by . The system thus determines the boundary by variation of the domain parameter. In addition, the gradient of the domain parameter automatically determines the inward normal vector of the contour level sets of ; see Fig. 1(c). Our goal is to solve partial differential equations within the region where while imposing boundary conditions at the narrow transitioning interfacial region where . However, the convention can be reversed such that the domain is defined by , in which case the following derivation could be modified by replacing with accordingly. This could be used to solve a problem where multiple equations govern different regions within the computational domain. Furthermore, these equations can be coupled through the shared boundary conditions, making the method highly versatile.

## 3 Formulation

### 3.1 General Approach

The general approach is as follows. The domain parameter describes the domain of interest ( inside the domain, and outside). The transition between the two values described is must be smooth so that the gradient is well defined. In this work, we have assumed the domain parameter to take the form of a hyperbolic tangent function for three reasons. First, it can be numerically implemented with ease. Second, it is consistent with the solution to the phase field equations, and thus this choice allows coupling of the two approaches; that is, one could for example simulate microstructural evolution using the phase field model, but impose interfacial flux boundary condition using the smoothed boundary method. Third, when given a non-smooth microstructural data, one could use the phase field equations to obtain the domain parameter from the discontinuous data. Other forms of domain parameters are possible, as long as the transition is monotonic and the gradient of the domain parameter has a narrow peak in the interfacial region. A better convergence could possibly be obtained using a function that has a more confined gradient; however, examining other forms of the domain parameter is beyond the scope of this paper.Ó As an example, we consider the Laplacian of the function, . As the first step in deriving formulation for the Neumann boundary condition, is multiplied by the domain parameter, . Using identities of the product rule of differentiation such as:

 ψ∇2H=∇⋅(ψ∇H)−∇ψ⋅∇H, (1)

we obtain terms proportional to . Because the inward unit normal of the boundary (pointing to the regions where ), , is given by , such terms can be written in terms of , and thus reformulated to be the Neumann boundary condition imposed on the diffuse interface.

Similarly, to derive the smoothed boundary formulation for the Dirichlet boundary condition, the differential equation is multiplied by the square of the domain parameter. Again using mathematical identities, , where , we obtain:

 ψ2∇2H=ψ∇⋅(ψ∇H)−[∇ψ⋅∇(ψH)−H|∇ψ|2]. (2)

Note that , associated with appearing in the last term, is the boundary value, , imposed on the diffuse interface. Specific details of the derivation depend on the equation to which the approach is applied, and we therefore provide four examples below.

### 3.2 Diffusion Equation

The first example is the diffusion equation with Neumann and/or Dirichlet boundary conditions. The Neumann boundary condition is appropriate, for example, as the no-flux boundary condition, whereas the Dirichlet boundary condition is necessary when the diffusion equation is solved with a fixed concentration on the boundaries. For Fick’s Second Law of diffusion, the original governing equation is expressed as:

 ∂C∂t=−∇⋅→j+S=∇⋅(D∇C)+S, (3)

where is the flux vector, is the diffusion coefficient, is the concentration, is the source term, and is time. Instead of directly solving the diffusion equation, we multiply Eq. (3) by , the domain parameter that describes the domain in which diffusion occurs, and use the identity to obtain the smoothed boundary formulated diffusion equation:

 ψ∂C∂t=∇⋅(ψD∇C)−∇ψ⋅(D∇C)+ψS. (4)

Next, we consider the boundary condition in this formulation. The Neumann boundary condition is the inward flux across the domain boundary, mathematically the normal gradient of at the diffuse interface, and is treated as:

 DBN≡D∂C∂n≡→n⋅→j=−∇ψ⋅(D∇C)|∇ψ|, (5)

where is the unit inward normal vector at the boundary defined in the diffuse interface description. Note that the flux at the interface is equal to . Equation (5) is rearranged to become and substituted back into Eq. (4); thus, we obtain:

 ∂C∂t=1ψ∇⋅(ψD∇C)+|∇ψ|ψDBN+S, (6)

with the Neumann boundary condition appearing in the second term. When a no-flux boundary is imposed, the second term vanishes and the resulting equation is the same as that proposed in Refs. [1, 2, 3, 20, 22].

To impose the Dirichlet boundary condition, we manipulate the original governing equation in a procedure similar to the derivation of Eq. (6). Multiplying Eq. (4) by and using the identity to replace the second term, we obtain:

 ∂C∂t=1ψ∇⋅(ψD∇C)−1ψ2D[∇ψ⋅∇(ψC)−BD|∇ψ|2]+S, (7)

where is the Dirichlet boundary condition imposed at the diffuse interface to replace , associated with in the third term. The convergence to the imposed Neumann and Dirichlet boundary conditions is shown in Appendix A.

In this method, the boundary gradient, , and the boundary value, , are not specified to be constant values; they can vary spatially and/or temporally or be functions of or other parameters. In addition, it is convenient to use weighting factors to combine Eqs. (6) and (7) to impose Neumann and Dirichlet boundary conditions simultaneously to yield mixed (or Robin) boundary conditions. The equation then becomes:

 ∂C∂t=1ψ∇⋅(ψD∇C)+|∇ψ|ψDBNWN−Dψ2[∇ψ⋅∇(ψC)−BD|∇ψ|2]WD+S, (8)

where and are the spatially dependent weighting factors for the Neumann and Dirichlet boundary conditions, respectively (). These factors can be a linear combination when imposing Robin boundary conditions or be employed to impose Neumann and Dirichlet boundary conditions at different regions of the interface. Moreover, a small nonzero value () should be added to the domain parameter appearing in the denominators to avoid singularities resulting from the terms and in regions where .

### 3.3 Coupled Surface-Bulk Diffusion

The second example demonstrates that surface diffusion can be incorporated into the smoothed boundary equation derived above. For this case, we take the set of equations that includes the surface reaction, surface diffusion and bulk diffusion to describe an oxygen reduction model in a solid oxide fuel cell (SOFC) cathode [40]. The oxygen-vacancy concentration, , on the cathode surface is governed by Fick’s Second Law for surface diffusion:

 −Db∂C∂n=κC−Ds∇2% sC+L∂C∂t, (9)

where is the coordinate along the inward unit normal vector of the surface and is the surface Laplacian. The parameters , , , and are the bulk diffusivity, reaction rate, surface diffusivity, and accumulation coefficient, respectively. Thus, the term on the left-hand side represents the flux from the bulk, and the terms on the right-hand side represent the surface reaction, surface diffusion and surface concentration accumulation, respectively [40]. Here, these parameters are all assumed to be constant for simplicity. In the bulk of the cathode material, the oxygen-vacancy diffusion is governed by Fick’s Second Law for bulk diffusion:

 ∂C∂t=Db∇2C. (10)

To simulate the oxygen-vacancy concentration evolution in the cathode, the two diffusion equations, Eqs. (9) and (10), are coupled and must be solved simultaneously, with the flux normal to the cathode surface as a common boundary condition. Recently, a similar set of equations was formulated using anothor diffuse interface approach combined with an asymptotic analysis [30], which led to two differential equations coupled by a common boundary condition. We show below that we can eliminate the need for solving two separate equations by applying the smoothed boundary formulation described herein to obtain a single equation that governs both surface and bulk effects.

Similar to the derivation of Eq. (6), we first multiply Eq. (10) by and apply the product rule of differentiation to obtain the bulk diffusion equation containing a boundary term, , similar to Eq. (4). As in Eq. (5), the normal gradient at the diffuse interface is defined by . Substituting this relationship back into Eq. (9) and rearranging terms gives:

 ∇ψ⋅∇C=|∇ψ|Db[κC−Ds∇2sC+L∂C∂t]. (11)

Using this relation, we replace the boundary term to obtain the smoothed boundary formulated equation:

 ∂C∂t=1ψDb∇⋅(ψ∇C)−|∇ψ|ψ[κC−Ds∇2sC+L∂C∂t], (12)

which combines the bulk diffusion and surface diffusion terms into a single equation used in the examples presented in Sections 4.2 and 5.1. Again, a small nonzero value should be added to the domain parameter appearing in the denominators. In the bulk ( and ), Eq. (12) reduces to Eq. (10). When the interfacial thickness approaches zero, it reduces to Eq. (9) at the interface (), as proven in Appendix A.

The surface Laplacian () is calculated according to the surface gradient given by:

 ∇s=(I−n⊗n)∇=(I−∇ψ|∇ψ|⊗∇ψ|∇ψ|)∇, (13)

where is the unity tensor, ‘’ is the Dyadic product, and is the inward unit normal vector of the diffuse interface, as used in Ref. [30]. (In indicial notation, the surface gradient is expressed as: , where is the Kronecker delta ( if , and if ). The repeated indices indicate summation over the index. See Appendix B for details.) To simulate only surface diffusion on a diffuse-interface-described geometry, one can simply eliminate the bulk-related and reaction terms in Eq. (9) to obtain , such that the concentration evolves only over the interfacial region.

### 3.4 Mechanical Equilibrium

The smoothed boundary method can also be applied to the mechanical equilibrium equation. When a solid body is in mechanical equilibrium, all forces acting on the body are balanced in all directions, as represented by , where is a stress tensor component, which is the force per unit area along th axis on the surface whose normal vector is along the th axis. Repeated indices indicate summation over the index. For linear elasticity, the stress tensor is given by the generalized form of Hooke’s Law: , where is the elastic constant tensor and is a scalar body force, such as thermal expansion () or misfit eigenstrain (, where and are the lattice constants of the precipitate and matrix phases, respectively), depending on the governing physics. The total strain tensor is defined by the gradients of displacements as , where is the infinitesimal displacement in the th direction. Substituting Hooke’s Law and the total strain back into the mechanical equilibrium equation gives:

 ∂∂xjCijkl12(∂uk∂xl+∂ul∂xk)=∂∂xj(ρCijklδkl). (14)

Multiplying Eq. (14) by the domain parameter that distinguishes the elastic solid region () from the environment () and using the product rule of differentiation yields the smoothed boundary formulation:

 (15)

see Appendix C for details of the derivation.

The traction exerted on the solid surface is defined by , where is the inward unit normal of the solid surface. (In indicial notation, and .) Therefore, the surface traction force is given by:

 (16)

Substituting Eq. (16) into Eq. (15) yields the smoothed boundary formulated mechanical equilibrium equation with a traction boundary condition on the solid surface:

 ∂∂xj[ψCijkl12(∂uk∂xl+∂ul∂xk)]+|∇ψ|Ni=∂∂xj(ψρCijklδkl), (17)

where can be treated as an effective body force in the th direction.

For linear elasticity problems with prescribed displacements at the solid surface, one can perform the smoothed boundary formulation, as in the derivation of the Dirichlet boundary condition in Section 3.2, by multiplying Eq. (14) by and using the product rule to obtain:

 ψ∂∂xj[ψCijkl12(∂uk∂xl+∂ul∂xk)]−{(∂ψ∂xj)[Cijkl12(∂ψuk∂xl+∂ψul∂xk)]−(∂ψ∂xj)Cijkl12(uk∂ψ∂xl+ul∂ψ∂xk)}=ψ2∂∂xj(ρCijklδkl), (18)

where the displacements and appearing in the third term on the left-hand side should be the boundary values of the displacements at the solid surface; see Appendix C for the derivation. A similar formulation for solving the mechanical equilibrium equation within a domain defined by a phase-field-like order parameter can also be obtained by the asymptotic approach by matching terms of different orders [41].

### 3.5 Equations for Phase Transformations with Additional Boundaries

Phase transformations affected by a mobile or immobile surface or other boundaries are of importance in many materials processes, including heterogeneous nucleation that occurs at material interfaces [18, 19]. Maintaining a proper contact angle at the three-phase boundary (where the interface between the two phases meets the surface) is necessary to capture the dynamics accurately because the contact angle represents the difference between the surface energies (tensions) of the different phase boundaries. Although researchers have previously developed methods for imposing contact-angle boundary conditions on sharp domain walls [18, 19], here we show that a similar model with diffuse domain walls can be obtained simply by applying the approach described above. Below, we assume that the boundary is immobile, but this assumption can be easily removed by describing the evolution of the domain parameter as dictated by the physics of the system.

In the Allen-Cahn and Cahn-Hilliard equations of the phase field model, the total free energy has the following form [34, 35]:

 F=∫Ω[f(ϕ)+ϵ22|∇ϕ|2]dΩ, (19)

where is the phase field order parameter commonly used to define different phases, is a double-well free energy functional (in terms of ), is the gradient energy coefficient, and is the domain of interest. At the extremum of the functional , the variational derivative of the total free energy vanishes: . This requirement provides the following conditions: , which can be reformulated as , by multiplying both sides by . We thus find a useful equality for deriving the contact angle boundary condition: ; see Appendix D for details.

In the smoothed boundary method, we introduce a domain parameter to incorporate boundary conditions into the original governing equation. As mentioned previously, the level sets of this domain parameter describe the diffuse boundaries and should satisfy . On , we impose a contact angle, , such that , where is the unit normal vector of the phase boundary (pointing to regions where ). We can thus derive the following equation for the boundary condition:

 ∇ψ⋅∇ϕ=−|∇ψ|cosθ√2fϵ. (20)

This contact-angle boundary condition is similar to that suggested by Warren et al. [19] for contacting a sharp interface, in which a Dirac delta function replaces .

The chemical potential that drives the morphological evolution is defined by the variational derivative of the total free energy of the system: . We can apply the smoothed boundary formulation to the chemical potential by multiplying it by the domain parameter and applying the product rule to obtain:

 μ=∂f∂ϕ−ϵ2ψ∇⋅(ψ∇ϕ)+ϵ2ψ∇ψ⋅∇ϕ=∂f∂ϕ−ϵ2ψ∇⋅(ψ∇ϕ)−ϵ|∇ψ|ψ√2fcosθ, (21)

where Eq. (20) was used in the third term.

For a nonconserved order parameter in the phase field models, the evolution is governed by the Allen-Cahn equation [36] (also known as the time-dependent Ginzburg-Landau equation [37]), in which the order parameter evolves according to the local chemical potential variation:

 ∂ϕ∂t=−Mμ=−M(∂f∂ϕ−ϵ2ψ∇⋅(ψ∇ϕ)−ϵ|∇ψ|ψ√2fcosθ), (22)

where is the mobility coefficient.

For a conserved order parameter, the evolution of the order parameter is governed by the Cahn-Hilliard equation, where the rate of the order parameter change is equal to the divergence of the its flux, which is proportional to the gradient of the chemical potential [34, 35], . The smoothed boundary formulation (derived in a similar manner as in the derivation in Section 3.2) is given by . Note that is the flux of the conserved order parameter; thus, the second term represents the fluxes normal to the domain boundary (equivalent to Eq. (5)). The smoothed boundary formulation of the Cahn-Hilliard equation is thus written as:

 ∂ϕ∂t=1ψ∇⋅[ψM∇(∂f∂ϕ−ϵ2ψ∇⋅(ψ∇ϕ)−ϵ|∇ψ|ψ√2fcosθ)]+|∇ψ|ψJn, (23)

where . For a closed system, is zero. Note that a small nonzero value should to be added to the domain parameter in the denominators in Eqs. (22) and (23) to avoid division by zero.

## 4 Validation of the Presented Approach

We herein demonstrate the validity and accuracy of the approach introduced in Sections 3.2 and 3.3 and the phase transformation with the presence of additional surface in Section 3.5.

### 4.1 1D Diffusion Equation

First, we performed a 1D simulation to demonstrate that the Neumann and Dirichlet boundary conditions were satisfied on two different sides of the domain. Fick’s second diffusion equation with given source and sink terms was solved within the domain defined by . The diffusion coefficient was set at 1, and the source and sink strengths were 0.02 and 0.01, respectively. On the right boundary of the diffusion domain, the gradient of was set at -0.1, and on the left boundary, the value of was set at 0.4. We selected the 1D computational box for and used a hyperbolic tangent function for the continuous domain parameter :

 (24)

where is the coefficient for adjusting the interfacial thickness. The interfacial thickness is given approximately by where the interfacial region is defined by the range, . The left and right interfaces are located at and , respectively. We applied the smoothed boundary formulation, as in the derivation of Eq. (8), to reformulate the original diffusion equation, , to:

 (25)

where is the Heaviside function used to specify the choice of the boundary condition and is the midpoint of the diffusion domain. Therefore, the second and third terms only apply to the right and left interfaces, respectively. The initial concentration was everywhere in the computational box. A standard central finite difference scheme in space and an Euler explicit time scheme were employed in the simulations.

Figure 2 shows the concentration profiles recorded at four different times (in solid blue lines). The domain parameter is plotted in the red line (the red circular markers indicate the position of grid points). The computational box was discretized to 1,200 grid points (), and was taken to be , such that the interfacial thickness is approximately . The parameters are given as Case 1b in Table 1. On the right interface, it can be clearly observed that at all times (except for a rapid change from to in the very early transient period). In the early period, the concentration even took negative values to satisfy the gradient boundary condition imposed at the right interface. In contrast, the concentration remained at 0.4 at the left interface during the entire diffusion process (except in the very early transient period, during which changed from 0 to 0.4). The analytical solution for the original sharp interface equation is also plotted for comparison, showing excellent agreement between the two methods. This result clearly demonstrates that both Neumann and Dirichlet boundary conditions are satisfied on the diffuse interfaces, and the smoothed boundary formulated equation reproduces the same result to the corresponding sharp interface version.

To further analyze the effects of interfacial thickness and discretization resolution on the smoothed boundary method, various simulations were conducted. In the first case, various interfacial thicknesses were selected, as in Case 1 in Table 1, while the grid size was kept at . Figure 3(a) shows the concentration distributions at (nearly equilibrium) for (Case 1b) and (Case 1f), for which the interfaces approximately span and , respectively. It is clear that the calculated concentration deviates farther from the analytical solution when the interfacial thickness is greater (as shown in the derivation in Appendix A). Figure 3(b) illustrates the relative errors during concentration evolution for various values of . Here, the relative error is defined by the root-mean-square deviation (between the smoothed boundary result and the analytical solution) divided by the average analytical concentration. The results clearly show that the error increases as the interfacial thickness increases. As the concentration evolution approaches equilibrium, the errors also converge to their equilibrium values, as listed in Table 1. A scaling of the error to the interfacial thickness is observed as is varied from to . In addition, the deviation between the smoothed boundary results and the analytical solution is much larger near the left boundary than near the right boundary, indicating the error associated with a Dirichlet boundary condition is larger than that with a Neumann boundary condition; see Figs. 3(c)–(d).

In the second case, we examined the effect of varying without changing the number of grid points across the interface. This was accomplished by selecting various grid sizes while maintaining the ratio of interfacial thickness to the grid size at . Results similar to Case 1 were obtained (see Case 2 in Table 1 and Figs. 4(a)–(b)). Error increases with . Since the resolution of the interface is unchanged (i.e., the number of points across the interface is fixed), it implies that the increased interfacial thickness is the dominant source of error. However, the errors are in general smaller than those in Case 1, which may be due to the fact that the interfacial thickness is effectively reduced when the resolution is decreased (the parts of the interfacial regions where is near 0 and 1 are not resolved by large grid spacing). The same reason may explain the steep drop in the error for in Case 1a in Table 1, for which the rapid transition of is not properly resolved by the discretization.

In the third case, we selected various grid sizes to examine the effect of the resolution across the interface while maintaining the interfacial thickness (specifically, fixing at ). The results show that error decreases when a larger grid size is selected; see Case 3 in Table 1 and Figs. 4(c)–(d). This can be understood as follows. As we observed earlier, the smoothed boundary formulation reduces to the bulk partial differential equation far from the interface, where the gradient of the domain parameter vanishes. In the interfacial region, the bulk term and the boundary term together set the boundary condition, as shown in Appendix A. In between, there is a region where the bulk equation is affected by the boundary term, which is small because the gradient is small, but not negligible. When the resolution is sufficiently low, the domain parameter in these regions take the bulk values and vanishing the boundary term and thus increasing the accuracy. As can be observed in Fig. 4(d), the error behavior in such case is very different from other cases. In the specific example presented here, the discretized interface at the low resolution of is nearly a Heaviside step function, which yields smaller error than the high-resolution cases. (When the resolution is high enough, the error is not affected by the resolution.) Therefore, in the 1D case, an interface does not need to be fully resolved, and in fact the accuracy can be increased by not doing so. However, we found that numerical instability ensues when the resolution is further reduced. It has been determined that at least one point with an intermediate value between the two bulk values is required in order to achieve numerical stability.

Above argument applies to only 1D case or when interfaces are very flat in multiple dimensions. In one dimension, the curvature of the interface is zero (i.e., the interface is flat), and therefore the effect of curvature can be neglected and a good resolution across the interface (which provides the smoothness of the curved interface) is not required. This is not the case when sufficiently large curvature is present, and thus smoothed interface with about three grid points are required to obtain accurate results for 2D or 3D calculations.

In the fourth case, we varied the value of (the small value added to the denominators to avoid division by zero) from to , while the grid size and interfacial thickness were maintained at and , respectively. In practice, a smaller would lead to a less stable numerical implementation because the values of or become much larger, which requires a much smaller time step size. The results show that the error quickly converges to a small value when is smaller than ; see Case 4 in Table 1. This suggests that once is small enough to yield converged results, further reduction is unnecessary and should be avoided so that a larger time step can be employed.

In summarizing the above 1D test simulations, we found that the interfacial thickness is the dominant source of error. The errors are less sensitive to the resolution of the finite-differencing discretization (selection of ) and the parameter for singularity control (selection of ). When the diffuse interface is properly resolved, the error scales with the interfacial thickness. Moreover, in general, the error that results when a Dirichlet boundary condition is imposed is larger and more sensitive to the interfacial thickness than when a Neumann boundary condition is imposed. This behavior can be understood from the results of analysis in Appendix A, where the scaling of the errors can be found in Eqs. (37) and (39).

### 4.2 Surface Diffusion and Bulk Diffusion in a Cylinder

To further demonstrate the validity of the smoothed boundary method, we applied the method to simulate oxygen-vacancy diffusion in a cylinder, for which a cylindrical coordinate grid system was used. We solved the coupled surface-bulk diffusion problem using both the smoothed boundary and the original sharp interface formulations in the same grid system for comparison. For the smoothed boundary method, we used a hyperbolic tangent function, , of the continuous domain parameter to define a cylinder, where is the radial position, is the axial position, and is the cylinder radius. Therefore, the cylinder surface () where surface reaction and surface diffusion occur, is located at , the solid region () for bulk diffusion is defined at , and the environment () is defined at . We selected the cylinder radius and the cylinder axial length to be 1 and 12, respectively. The grid sizes were selected to be and , such that the cylinder contains 57 and 300 grid points in the radial and axial directions, respectively. (The computational box is larger than the cylinder in the radial direction, and contains 75 and 300 grid points in the radial and axial directions, respectively.) The interfacial thickness was selected to be by setting . Equation (12) was solved using a standard central finite difference scheme in cylindrical coordinate system and an Euler explicit time scheme; see Appendix B for the discretization scheme. The parameters for the diffusion equation were selected to be , , and . The boundary conditions at the computation box boundary were set at at =0 and at , with no gradient on the remaining two sides. For comparison, the original sharp interface equations, Eqs. (9) and (10), were solved using the same discretization scheme with the same grid system and resolution. For this case, the surface-reaction-diffusion boundary condition, Eq. (9), was explicitly imposed at the 57th grid points in the radial direction. The concentration evolution is implemented as follows. First, the surface concentration is updated by Eq. (9) according to the normal flux at the cylinder surface calculated from the normal gradient of the surface concentration obtained from Eq. (10). Next, the normal surface flux is calculated using Eq. (9) with the updated surface concentration. The cylinder concentration is then evolved according to Eq. (10) with the normal flux boundary condition. This procedure is repeated within the Euler explicit time scheme for the concentration evolution. For the smoothed boundary formulation, we simply solve a single equation that automatically includes coupled bulk and surface diffusion, Eq. (12).

Figures 5(a) and (b) show the steady-state concentration profiles of the sharp interface version for and , respectively. The concentration decays along the axial direction according to boundary values prescribed at the box boundaries. The diffusion front bends because of the surface reaction, such that the concentration is lower near the cylinder surface. Shown in Figs. 5(c) and (d) are the corresponding smoothed boundary results. For clarity, only the concentration in the region of is presented. The results from the two methods are in excellent agreement, clearly demonstrating the utility and validity of the smoothed boundary method for incorporating two sharp-interface equations into one smoothed boundary equation.

To further examine the effect of interfacial thickness, we included two other radial grid sizes in the simulations, i.e., and , such that the cylinders contain 29 and 15 radial grid points, respectively. By selecting the radial grid sizes in this way, each radial grid point in a lower-level resolution (thicker interface) overlaps with every other grid point in a higher-level resolution (thiner interface). The diffuse interface is maintained to span , and the axial grid size is kept at . Hereafter, we refer the three interfacial thicknesses to as the “thin” (), “medium” () and “thick” () interface cases. Shown in Figs. 5(e) and (f) are the steady-state concentration profiles for the thick-interface results, corresponding to the cases in Figs. 5(a) and (b). The results are still in reasonably good agreement with the original sharp interface results, even though the interfacial thickness is approximately 29.2% of the cylinder radius.

The relative errors of the thin, medium, and thick interface smoothed boundary results are plotted in Fig. 6. The relative errors are calculated by dividing the differences between the smoothed boundary and sharp interface results by the average concentration of the sharp interface results. The average concentration is calculated for the active region between the plane at and the plane on which the maximum concentration is 0.01. Note that only the errors within in the cylinder defined by are considered. For the thick-interface smoothed boundary result, the maximum local errors for and are approximately and (see Figs. 6(e) and (f)), whereas the average relative errors are and , respectively; see Table 2. The average relative errors, denoted by in Table 2, are calculated by dividing the root-mean-square deviation between the smoothed boundary and the sharp interface results by the average sharp interface concentration in the cylinder. The root-mean-square deviation and average concentration are calculated in the cylindrical coordinate system. As expected, the error increases as the interface becomes thicker (i.e., as increases). However, in contrast to the 1D simple diffusion test in Section 4.1, the behavior of the error is inconsistent across the parameter sets; see Table 2. This relatively complicated error behavior may originate from the coupling of the bulk and surface diffusion equations. In addition to the effect of the interfacial thickness, the error also increases with a larger reaction coefficient , which may be explained by the increase of the scaling coefficient for the error ( defined above Eq. (36) in Appendix A) when the given boundary value is larger. In addition, the gradient of the concentration near the boundary increases in magnitude with increasing boundary condition value, which can lead to a larger error.

One interesting phenomenon is observed in the high results ( and ). Although the errors in the bulk greatly increase with the interfacial thickness, the errors at the surface remain small; see Fig. 6(f) and Table 2. This indicates that the error originates from the boundary condition affecting the bulk solution, rather than from an increased error in the boundary condition value. The thicker interface thus leads to a larger bulk region that is affected by the boundary condition. Therefore, we compare the error associated with the bulk region and with the boundary condition. Here, the bulk errors, denoted by in Table 2, are calculated by the same method as the average relative error but exclude the grid points on the nominal cylinder surface (). The surface errors, denoted by in Table 2, are calculated by the same method but with only the grid points at the nominal cylinder surface. In the case where , the error at the surface even decreases with interfacial thickness.

### 4.3 Contact-Angle Boundary Condition

We performed simple 2D simulations to validate the smoothed boundary formulation for the contact-angle boundary condition at the three-phase boundary. Equations (22) and (23) were tested for nonconserved and conserved order parameters, respectively. The equations were solved using the central finite difference scheme and the Euler explicit time scheme. The computational box sizes are and , and the parameters used are and . A simple common double-well function was selected for the bulk free energy functional, , such that the steady-state phase field order parameter profile is determined by , where is the coordinate variable indicating the distance to the phase boundary, and the characteristic thickness of the diffuse interface is determined by . By setting , the characteristic thickness is controlled by , and the phase field interfacial energy is maintained at a constant value, . A horizontal diffuse-interface flat substrate surface is defined by the hyperbolic tangent function , such that is at , and gradually transitions from 0 to 1 from below to above the substrate surface. Here, we have two diffuse interfaces: one for the phase field order parameter and the other for the smoothed boundary domain parameter. Both thicknesses can affect the accuracy when imposing contact-angle boundary conditions. To verify the contact angle boundary conditions, various combinations of substrate surface thickness and phase boundary thickness were selected by adjusting the values of and . The initial phase boundary was placed vertically in the middle of the domain (), with Phase 1 () and Phase 0 () on the left and right halves, respectively. On the computational box boundaries, the normal gradients of the phase field order parameter were set at zero: at and , and at and , which can be interpreted as the no-flux boundary conditions.

In the first set of simulations, we evolved Eq. (22) for a nonconserved order parameter with a 60 contact angle. The result clearly shows a 60 contact angle at the three-phase boundary, as specified; see Fig. 7(a). The angle can be measured at the intersection between the two contours of and , as shown in Figs. 7(b) and (c). The 60 angle is maintained during the entire evolution, except for the very early transient period, when the contact angle changes from the initial 90 angle to the prescribed 60 angle. Because of the contact-angle boundary condition, the initially flat phase boundary bends and creates a negative curvature in Phase 1. As a result, the phase boundary moves toward Phase 0. Once the phase boundary evolves to a circular arc with a uniform curvature everywhere (other than regions in contact with the substrate), it moves at a uniform constant speed in a steady-state motion, and eventually only Phase 1 remains in the system.

In the second set of simulations, we evolved Eq. (23) for a conserved order parameter in a closed system with and a 120 contact angle. As expected, the phase boundary intersects the substrate surface at a 120 contact angle; see Fig. 7(d)–(f). In contrast to the Allen-Cahn-type dynamics, because of the conservation of the order parameter, the phase boundary near the substrate moves toward the left, whereas the phase boundary away from the substrate moves in the opposite direction. As a result, the phase boundary deforms into a curved shape. When the system reaches its equilibrium state, the phase boundary forms a circular arc with a uniform curvature everywhere (except where the phase is in contact with the substrate), such that the total surface energy is minimized; see Fig. 7(d) for .

Table 3 lists the average values of calculated by at the grid points within the three-phase boundary region defined by and in the steady state (i.e., when the phase boundary becomes a circular arc). These results again clearly show that the error (for a given phase boundary thickness) increases as the interfacial thickness increases. This can be understood based on the analysis of the 1D test results in Section 4.1 since the substrate surface is assumed to be flat in this test. If the substrate interface is curved, the resolution of the interface will have more influence on the error.

In contrast to the effect of domain boundary thicknesses, the error is relatively insensitive to the phase boundary thickness once the phase boundaries are properly resolved; see cases with in Table 3. However, when the phase boundary is too thin, the error tends to increase because of the loss of resolution in the phase-field-order-parameter gradient; see cases with in Table 3. In general, the results demonstrate that the contact-angle boundary condition is well imposed using the presented method. Even when the domain boundary thickness is as large as 16.74 grid spacings with , the contact angle only deviates less than 2 from the imposed values as long as the phase boundary is properly resolved. Here, the error is sensitive to the resolution of the interface because the phase boundary is curved unlike the substrate surface.

In addition to the contact angle, a no-flux boundary condition for a conserved order parameter is implicitly imposed at the substrate surface. The error associated with such a boundary condition was evaluated by examining the overall change in the value of the total order parameter. The conservation of the order parameter was met within a numerical error (well below 1% in most cases) in these validation simulations; see Table 3.

## 5 Applications

Although the details of the scientific calculations performed applying these methods to problems in materials science will be published elsewhere, it is worth presenting some of the results herein to demonstrate the potential of the method.

### 5.1 Oxygen-Vacancy Diffusion in SOFC Cathode

The first example is ionic transport through a complex microstructure. Here, ion diffusion is driven by a sinusoidal voltage perturbation. For the steady-state solution, the time dependence of the form , where is the angular frequency and , can be removed as in the equation derived by Lu et al. [40]. For this case, the smoothed boundary formulated equation is obtained from Eq. (12) to:

 ∇⋅(ψDb∇~C)−|∇ψ|(κ~C−Ds∇2s~C)=iψω~C, (26)

where is the complex concentration amplitude, which consists of real and imaginary parts and includes the amplitude of the concentration wave and the phase shift. Note that the surface accumulation term is ignored () here because its magnitude is usually very small in comparison with the bulk concentration [40]. This equation can be solved by an alternating direction line relaxation (ADLR) method in a second-order central-difference scheme in space; see Appendix B for the numerical implementation.

In this work, we adopted an experimentally reconstructed complex microstructure, the porous ceramic cathode and nonporous ceramic electrolyte of an SOFC, as the input geometry. The microstructure data is stored as a 3D array consisting of voxels that indicate the electrolyte (gadolinia-doped ceria: GDC), cathode (lanthanum strontium chromite: LSC) and pore phases by different values. To emphasize the convenience of image-based smoothed boundary simulations, we treat the center of each voxel as the location of the grid points in the calculation without further enhancement of the resolution from our initial reconstructed microstructure. For very high-accuracy scientific calculations, one can easily enhance the resolution by refining the grid sizes. To smooth the voxelated, discontinuous data, we first employed a level set distance function method [42] to determine the distances between grid points and the solid-pore interface, and then computed the hyperbolic tangent of the distance function to obtain the domain parameter profile; see Appendix E for details.

For simulations of the concentration distribution in the porous cathode, the regions containing nonporous electrolyte are excluded, such that the computational box only consists of grid points. The grid spacing is set at . The boundary conditions along the main diffusion direction (the -axis) on the computational box are and at , and and at . The boundary conditions on the remaining four sides are zero-gradient for both the real and imaginary parts. As a demonstration of the method, the length scale and physical material properties are nondimensionalized. Figure 8 shows the steady-state concentrations for the cases in which surface diffusion is excluded ( and ) and included ( and ) with and direct current (DC) loading (). In these cases, the imaginary part vanishes, and the solution of the real part is equivalent to that of a homogeneous Helmholtz-like equation with the right-hand side of Eq. (26) equal to zero. As shown in Fig. 8, the concentration decays from 1 to 0 along the -axis over the complex cathode microstructure to satisfy the boundary conditions imposed on the box boundaries and at the cathode-pore interfaces. The utilization lengths (i.e., the length over which the cathode material is active) of the two cases are similar, as predicted by Lu et al. [40] for a cathode with simplified cylindrical geometry, in which the effective diffusivities under DC loading with and without surface diffusion are found to be similar for the parameters given above. However, a slight difference in the concentration distributions of the two cases can be observed. Because of the faster transport path along the surface, the diffusion front with surface diffusion (Fig. 8(b)) is more planar compared with that without surface diffusion (Fig. 8(a)).

Figure 9 shows the real and imaginary parts of the steady-state concentration amplitude for the cases in which , and , with alternating current (AC) loading of the angular frequencies of and . The boundary conditions on the computational box are the same as in the DC loading case above. In the low frequency case (Figs. 9(a) and (d)), the real part of the concentration, which represents the amplitude of the concentration wave decays and forms a planar diffusion front within the utilization length, where the material is active (). Additionally, a negative value of the imaginary part occurs in the regions where the real part decays because of the phase shift resulting from the delayed response. The magnitude of the imaginary part then decays back to zero toward the inactive region. In the high frequency case, the enhancement of concentration along the surface is observed due to surface diffusion; this is evident in Figs. 9(b) and (c), which show larger values of the real part of the concentration amplitude within the utilization length (). The real part of the concentration amplitude quickly decays from the surface into the bulk. In contrast to the enhanced real part at the cathode surface, the magnitude of the imaginary part is small near the surface because surface diffusion reduces the response time and the phase change is thus decreased. In an analogy to the low-frequency response, a negative imaginary part occurs in the region where the real part decays. The magnitude of the imaginary part decays toward the inactive region. This behavior can be more clearly discerned in the magnified views in Figs. 9(c) and (f).

The smoothed boundary method can also be used to impose Dirichlet boundary conditions on irregular surfaces. For example, if the ionic diffusivity in the electrolyte is assumed to be much larger than that in the cathode, the concentration in the electrolyte will be nearly uniform. To simulate this scenario, we impose a fixed concentration at the electrolyte-cathode contacting surface as the boundary condition. We used the experimentally reconstructed array that contains a porous cathode and a nonporous electrolyte as our input geometry. The voxelated data were smoothed to the hyperbolic tangent domain parameter profile by the level set distance function method mentioned in Appendix E. Here, three domain parameters are employed to define the three regions: electrolyte (), cathode (), and pore (). The smoothed boundary formulated governing equation is obtained by modifying Eq. (12) to:

 ∂C∂t=∇⋅(ψ2Db∇C)ψ2−|∇ψ2|ψ2[κC−Ds∇2sC]WN−Dbψ22[∇ψ2⋅∇(ψ2C)−|∇ψ2|2BD]WD, (27)

where the weighting factors are given by and , such that the Neumann boundary condition (surface reaction and surface diffusion) is imposed only at the cathode-pore interface (), and the Dirichlet boundary condition (a prescribed concentration value) is imposed only at the electrolyte-cathode interface (). The exponent determines the transition profiles from the Neumann to the Dirichlet boundary conditions in the regions of three-phase boundaries. We selected for this numerical simulation. On the computational box boundaries, we set at and the zero-gradient boundary condition for the remaining five sides.

The same material parameters used in the cases of Fig. 8 were selected. Figures 10 (a) and (b) illustrate the reconstructed SOFC complex microstructure and irregular surfaces defined by the values and gradients of the domain parameters, respectively. Figures 10(c) and (d) show the simulation results of the steady-state oxygen-vacancy concentration distributions with a fixed value of imposed at the cathode (LSC)-electrolyte (GDC) interfaces. The concentration distribution is very different from the ones shown in Fig. 8 because a larger portion of lateral diffusion occurs in the and directions, which results from the smaller contacting areas (compared to the cross-sectional area of LSC on the - plane in Fig. 8, where diffusion is mainly in the direction). As a result, the concentration drops rapidly within a short distance from the contacting areas, making the utilization length of the cathode material shorter and uneven.

### 5.2 Kirkendall-Effect Diffusion with a Moving Boundary

#### 5.2.1 Kirkendall-Effect-Induced Deformation Modeled by Navier-Stokes-Cahn-Hilliard Equations

The third application demonstrates the smoothed boundary method’s broad applicability by applying it to the coupled Navier-Stokes-Cahn-Hilliard equations [43, 44, 45, 46, 47, 48]. This particular formulation aims to solve diffusion problems with the Kirkendall effect with efficient and abundant vacancy sources and sinks in the bulk of a solid [49, 50, 51, 52, 53]. In this case, the solid experiences deformation because of vacancy generation and elimination. The Navier-Stokes-Cahn-Hilliard equations are coupled with the smoothed boundary formulation of the diffusion equation derived in Section 3.2 as a model of plastic deformation because of volume expansion and contraction resulting from vacancy flow.

When the diffusing species of a binary substitutional alloy have different mobilities, the diffusion fluxes of the two species are unbalanced, creating a net vacancy flux toward the side containing the fast diffuser. Here, we denote the quantities associated with the slow diffuser, fast diffuser and vacancy by the subscripts , , and , respectively. Because of the accommodation/supply of excess/depleted vacancies, the solid locally expands/shrinks [54, 55, 56, 57, 58] when maintaining the vacancy mole fraction at its thermal-equilibrium value. We treat the solid as a very viscous fluid [59, 60, 61, 62, 63] with a much larger viscosity than that of the surrounding environment. In this case, we solve the Navier-Stokes-Cahn-Hilliard equations to update the shape of the material as follows [64]:

 (28a) ∇⋅v=gV, (28b) ∂ψ∂t−v⋅∇ψ=M∇2(∂f∂ψ−ϵ2∇2ψ), (28c)

where is the effective pressure, is the viscosity, is the velocity vector, is the number of dimensions, the superscript denotes the transpose, is the Cahn number reflecting the capillary force compared to the pressure gradient, is the vacancy generation rate per unit volume, and is the domain parameter indicating the solid phase for diffusion. One great advantage in employing a phase-field type equation is that it automatically maintains the profile of the domain parameter, , in the form of a hyperbolic tangent function because it is the equilibrium solution for the phase field equation (Eq. (28c)). Note that here we ignore the inertial force in the Navier-Stokes equation to obtain Eq. (28a) because the deformation is assumed to be a quasi-steady-state process. The vacancy generation rate that results in the local volume change (dilatational strain) is given by , where is the mole fraction of the fast diffuser, is the thermal-equilibrium vacancy mole fraction (which is assumed to be maintained throughout the solid in this model), is the diffusivity for vacancy flux associated with , and is the lattice site density of the solid. Here, is taken to be . The evolution of the fast diffuser mole fraction is governed by the advective Fick’s diffusion equation, written as:

 ∂XB∂t−v⋅∇XB=∇⋅(DVBB∇XB)−XBgV, (29)

where is the diffusivity for the fast diffuser flux associated with , and the advective term accounts for the lattice shift because of volume change. Because diffusing atoms cannot depart from the solid region, a no-flux boundary condition is imposed at the solid surface. Thus, the smoothed boundary formulation of Eq. (29) is written as:

 ∂XB∂t−v⋅∇XB=∇⋅(ψDVBB∇XB)ψ−XBgV. (30)

As the concentration evolves, the shape of the solid is also updated by Eq. (28c) and by iteratively solving Eqs. (28a) and (28b) through the application of a projection method [65, 66]; see Appendix F for the numerical implementation.

The slow and fast diffusers are initially placed in the left and right halves of the solid, respectively. We use their theoretically calculated diffusivities for this simulation [67, 68, 69, 70]. Here, we calculate the slow diffuser atomic hop frequency based on the material parameters of aluminum at 600 K, and set the fast diffuser atomic hop frequency four times larger than that of the slow diffuser. Figure 11 shows snapshots of the mole fraction profiles (left column) and velocity fields (right column) from a 2D simulation. As the fast diffuser diffuses from the right to the left side, the vacancy elimination and generation cause contraction and expansion on the right and left sides, respectively. As a result, the initially rectangular slab deforms into a bottle-shaped object.

#### 5.2.2 Kirkendall Void Growth with Localized Vacancy Sources

In another scenario in which the vacancy diffusion length is comparable to or smaller than the distance between vacancy sources and sinks, the explicit vacancy diffusion process must be considered [58, 70, 71, 72]. In this case, vacancies diffuse in the same manner as the atomic species. In the bulk of a solid devoid of vacancy sources/sinks, the concentration evolutions are governed by:

 ∂XV∂t=∇⋅(DVV∇XV+DVB∇XB), (31a) ∂XB∂t=∇⋅(DBV∇XV+DVBB∇XB). (31b)

Because the solid surfaces are very efficient vacancy sources/sinks [64, 72], we impose the thermal-equilibrium vacancy mole fraction at the solid surfaces as the Dirichlet boundary condition for solving Eq. (31). In this case, the smoothed boundary formulation of Eq. (31) is given by:

 ∂XV∂t=1ψ∇⋅[ψ(DVV∇XV+DVB∇XB)]−Kψ2, (32a) ∂XB∂t=1ψ∇⋅[ψ(DBV∇XV+DVBB∇XB)]+XB1−XeqVKψ2, (32b)

where . Because the vacancy generation and elimination in this scenario only occurs on the solid surfaces, internal volume change in the bulk is not considered. Therefore, instead of using a plastic deformation model as in the previous case, we adopt a typical Cahn-Hilliard type dynamics to model the shape change:

 ∂ψ∂t=M∇2(∂f∂ψ−ϵ2∇2ψ)+∇ψ|∇ψ|⋅→JV1−XeqV, (33)

where is the vacancy flux, and the last term represents the normal velocity of the solid surfaces because of vacancy injection into or ejection from the solid.

An example of the results obtained using this approach is the growth of a void in a rod [72, 73, 74, 75]. The above equations were solved using a central difference scheme in space and an implicit time scheme (see Appendix G). The fast diffuser was initially placed in the central region while the slow diffuser filled the outer region. A void was initially placed off-centered in the fast-diffuser region, where a 1D study found to be the likely nucleation site [72]. The vacancy mole fractions were fixed at the void and cylinder surfaces. Figure 12 shows snapshots of the fast diffuser mole fraction profile (normalized to the lattice density) and the vacancy mole fraction profile (normalized to its equilibrium value). As the fast diffuser diffuses outward, vacancies diffuse inward from the rod surface to the void surface, causing vacancy concentration enhancement and depletion in the central and outer regions, respectively. To maintain the equilibrium vacancy mole fraction at the rod and void surfaces, vacancies are injected and ejected at those surfaces. As a result, the rod radius increases, and the void grows. Similar dynamics were examined using a sharp interface approach [72], but this new method provides the flexibility in geometry to examine cases where a void initially forms off-centered.

### 5.3 Thermal Stress

Solid oxide fuel cells (SOFCs) usually operate at temperatures near C. Evaluating the thermal stress resulting from the differences in thermal expansion and elastic moduli is important for analyzing mechanical failure. We expand the generalized mechanical equilibrium equation, Eq. (17), for a linear, elastic and isotropic solid. (Note that the derivation for the mechanical equilibrium equation is general and is not limited to isotropic solids. We selected an isotropic model because of the lack of available crystallographic information among the experimental data.) The equation is discretized in a central finite difference scheme and numerically solved by an ADLR solver; see Appendix H for details.

The thermal expansion rates of the ceramic electrolyte (GDC) and cathode (LSC) are taken to be [76] and [77], such that the thermal expansions at operation temperature are and , respectively. (Here, we have assumed arbitrarily that the composite material is relaxed at a reference temperature, and assumed an operation temperature of above the reference temperature.) We chose the elastic constants of GDC to be isotropic (), and the values are GPa, GPa and GPa, calculated from a Young’s modulus of 250 GPa and a Poisson’s ratio of 0.334 [78, 79]; see Appendix H. The LSC phase is softer than the GDC phase, and its elastic constant is also assumed to be isotropic. The values are selected to be GPa, GPa and GPa, based on a Young’s modulus of 200 GPa and a Poisson’s ratio of 0.3 [77]. As in Section 5.1, we again use domain parameters to indicate the GDC phase ( inside the GDC and outside the GDC) and the LSC phase ( inside the LSC and outside the LSC). The entire solid phase is then represented by the sum of the two phases, . The body force term and elastic constant tensor are replaced by an interpolated, spatially dependent thermal expansion and elastic constant tensor according to the domain parameters; see Appendix H. The solid surface is assumed to be traction-free, .

In this simulation, we selected the same computational box as in the case of Fig. 10(a), containing grid points in the , , and directions, respectively. Each grid point represents a voxel in the experimentally obtained microstructure. The yellow color indicates the LSC phase, and the semitransparent cyan color indicates the GDC phase; see Fig. 10(a). The grid spacing is nm, such that the computational box spans m. We assumed a rigid computational box with frictionless boundaries on the six sides, which means that on the two - planes, on the two - planes, and on the two - planes of the computational box boundaries, where , and are the displacements along the , and axes, respectively. While this set of boundary conditions is not realistic for SOFC material environment, we chose it for the demonstration purpose in order to avoid overlaps with a future publication of physically based SOFC simulations.

Shown in Fig. 13(a) are the calculated mean stress distributions resulting from thermal expansion in a confined sample. The mean stress is defined by: , where the stress components are calculated according to the method provided in Appendix H. Here, we choose mean stress to illustrate the effective pressure in the solid. A negative mean stress indicates that the region is under compression. Despite a complicated stress distribution observed because of the complex geometry, the overall magnitude of the mean stress is roughly between 2 and 4 GPa, which can be roughly estimated by the product of Young’s modulus and the thermal expansion with an enhancement resulting from the porosity of the solid. Additionally, an overall larger stress in the GDC phase is observed, reflecting the GDC phase is harder than the LSC phase. Figure 13(b) shows the mean stress in the LSC phase and Fig. 13(c) shows the mean stress on the GDC surface after rotating the volume 180 around the -axis. Three types of stress enhancements can be observed in the simulation result. At the cathode-electrolyte contacting surfaces, stress is enhanced because of the mismatch of thermal expansion and elastic constants between the two materials; see the red arrows in Figs. 13(b) and (c). The second is the concentrated stress observed at the grooves on the electrolyte surface (not contacting the cathode), as shown by the white arrows in Fig. 13(c). The third type is the stress concentration effect at the bottlenecks in the cathode phase, where the stresses are roughly larger by a factor of three to four compared to the overall value, as shown by the green arrows in Fig. 13(b). The simulation results demonstrate that the smoothed boundary method can properly capture the linear elasticity and the geometric effects of the system based on a diffuse-interface defined geometry.

### 5.4 Phase Transformations in the Presence of a Foreign Surface

The Allen-Cahn equation describes the dynamics of a nonconserved order parameter, which can be taken as a model for the ordering of magnetic moments [38] and diffusionless phase transformations that involve only changes in crystalline order [38]. This equation can also be used as a model for evaporation-condensation dynamics [38, 39]. Here, we use the Allen-Cahn equation to examine the evaporation of a droplet on a rough surface. The domain parameter was given a ripple-like feature, as shown in Fig. 14, having a hyperbolic-tangent-like profile continuously transitioning through the substrate surface ( above the surface, and below the surface). The droplet phase was placed on top of the boundary, and its shape was evolved by the smoothed boundary formulation of the Allen-Cahn equation, Eq. (22), using the standard central difference scheme in space and an Euler explicit scheme in time. The simulation was performed in two dimensions, using the parameters , and , with a domain size of and . The contact angle was set at 135, and a zero-gradient boundary condition of is set at the computational box boundaries.

The evolution of the droplet surface as it evaporates is illustrated in Fig. 14(a) as a contour () plotted at equal intervals of 270 dimensionless time units. The color change from blue to red indicates various times from the initial to the final stages, respectively. As the surface evolves, it is clear that the contact angle is maintained, as shown in Fig. 14(b). The dynamics of the motion of the three-phase boundary are interesting in that the velocity changes depending on the angle of the surface (with respect to the horizontal axis), which can be inferred from the change in the density of the contours. Because the interfacial energy is assumed to be constant, the droplet would prefer to have a circular cap shape. However, the contact angle imposes another constraint at the three-phase boundary. When the orientation of the surface is such that both of these conditions are nearly met, the motion of the three-phase boundary is slow as the droplet evaporates. When the orientation becomes such that the shape of the droplet near the three-phase boundary must be deformed (compared to the circular cap), the three-phase boundary moves very quickly, which leads to an unsteady motion of the three-phase boundary. In contrast, at the top of the droplet far from the substrate, the curvature is barely affected by the angle of the substrate surface; thus, the phase boundary there moves at a speed inversely proportional to the radius.

### 5.5 Motion of a Droplet Due to Unbalanced Surface Tensions

In another example application, we modeled a self-propelled droplet. Here, two different contact-angle boundary conditions are imposed on the right and left sides of the droplet placed on a flat surface. The smoothed boundary formulation of the Cahn-Hilliard equation, Eq. (23), is used with in this simulation. The domain sizes are and . The parameters and computational box boundary condition are the same as in Section 5.4. The contact angle on the right side of the droplet is set to 45 and that on the left side to 60 by imposing position-dependent boundary conditions. Note that this setup is equivalent to the situation in which the substrate-environment, droplet-substrate and droplet-environment surface energies satisfy the conditions of Young’s equation:

 γse−γsd=γdecos60∘  for the left side, (34a) γse−γsd=γdecos45∘  for the right side, (34b)

where , and are the interfacial energies of the substrate-environment, droplet-substrate and droplet-environment interfaces, respectively. Therefore, this model can be used to simulate a case where the surface energies are spatially and/or temporally dependent on other fields, such as surface temperature or surface composition. This specific case applies when the wetted substrate behind the droplet have a higher interfacial energy than the pristine substrate, as in Ref. [80].

The evolution of the droplet surface is illustrated in Fig. 15. The droplet initially has the shape of a hemisphere, with a 90 contact angle with the substrate surface. The early evolution is marked by the evolution of the droplet shape as it relaxes to satisfy the contact-angle boundary condition, as seen in Fig. 15(a). The droplet then begins to accelerate. Once the contact angle reaches the prescribed value, it is maintained as the droplet moves toward the right; see Fig. 15(b). In the steady state, the droplet moves at constant speed without other effects present. Such motions of droplets have been observed and explained as a result of an unbalanced surface tension between the head portion (with a nonwetting surface) and tail portion (with a wetting surface) because of the resulting spatially varying composition and composition-dependent surface energy [80].

Figure 16 shows the relaxation of an initially hemispherical droplet on an irregular substrate surface in a 3D simulation. The contact-angle boundary condition imposed at the three-phase boundary is 135. The computational box sizes are and . As shown here, the droplet changes its shape to satisfy the imposed contact angle, and the droplet evolves into a shape for which the total surface energy is minimized. The behavior favoring dewetting imposed by the contact angle () is properly reflected in the lifting of the droplet, as shown in Figs. 16(a)–(c) and (d)–(f). During this relaxation process, the three-phase boundary moves toward the center as the droplet-substrate contacting area decreases, as shown in Fig. 16(a)–(c). This model and numerical method has been applied to simulate a nickel particle coarsening process in the complex channel within supporting porous ceramic microstructure (consisting of yttria-stabilized zirconia) in SOFC anodes, and to estimate the degradation of the anode material during SOFC operation [81].

## 6 Discussion and Conclusions

In this work, we demonstrated a generalized formulation of the smoothed boundary method. This method allows Neumann, Dirichlet, or mixed boundary conditions to be imposed on a diffuse interface to solve partial differential equations within the region where the domain parameter uniformly equals 1. The derivation of the method, as well as its implementation, is straightforward. The method can be used to solve differential equations numerically without complicated and time-consuming structural meshing of the domain of interest, as the domain boundary is specified by a spatially varying function. Instead, any grid system, including a regular Cartesian grid system, can be used with this method.

This smoothed boundary approach is flexible in coupling multiple differential equations. In Section 3.3, we demonstrated how this method can be used to couple bulk diffusion with surface reaction-diffusion into a single equation while the two equations serve as complementary boundary conditions. In principle, this method can be used to couple multiple differential equations in different regions defined by different domain parameters. For example, if the physics within a domain defined by are governed by a differential equation , the overall phenomenon will be then represented by , where the subscript ‘’ denotes the th domain and represents the entire computational box. When sharing the diffuse interfaces between domains, the physical quantities can be interconnected as boundary conditions for each equation in each domain. Therefore, this method could be used to simulate coupled multiphysical and/or multiple-domain problems such as fluid-solid interaction phenomena or diffusion in multi-material polycrystalline solids.

We further demonstrated the capability of applying the smoothed boundary method to moving boundary problems in Section 5.2. When the locations of domain boundaries are updated by a phase-field-type dynamics such that the domain parameters remain uniformly at 1 and 0 on either side of the interface, the smoothed boundary method can be conveniently employed to solve partial differential equations with moving boundaries. In addition to the phase-field-type dynamics, the smoothed boundary method is also applicable to moving boundary problems implementing the level set method [32, 33, 82], with the domain parameter obtained simply by taking the hyperbolic tangent of the distance function.

In addition to Neumann and Dirichlet boundary conditions, we also showed the capability of the smoothed boundary method for specifying contact angles between phase boundaries and domain boundaries (Sections 4.3, 5.4 and 5.5). This type of boundary condition is difficult to impose using conventional sharp interface models.

Although the smoothed boundary method has many advantages, as shown in the results in Section 4 and in the derivations in Appendix A, the nature of the diffuse interface inevitably introduces an error proportional to the interfacial thickness because we expand an originally zero-thickness boundary into a finite thickness interface. This spread of interface also leads to another error source depending on the resolution of the rapid transition of the domain parameter across the interfacial region. When numerically solving the smoothed boundary formulated equations, properly capturing the gradient of the domain parameter across the interface becomes very important. Based on our experience, 3 – 5 grid points are necessary to properly resolve the diffuse interfaces while ensuring that the errors are well-controlled. Moreover, when solving time-dependent equations, one singularity occurs because of the terms and used to impose the Neumann and Dirichlet boundary conditions, respectively. In practice, a small value is necessarily added to to avoid singularity resulting from division by zero. In our simulations, the errors were quickly saturated when the value added to was selected to be smaller than when 3 – 5 grid spacings were used for the interfacial regions, which suggests it is unnecessary to select a smaller value for the singularity-control term. However, when solving time-independent equations, such as the mechanical equilibrium equation and the steady-state diffusion equation, there are no singular terms in the equations. The small additional term is then merely used to condition the matrix solver. In this case, it can be on the order of numerical precision, such as .

Based on the general nature of the derivation, the smoothed boundary method is applicable to generalized boundary conditions, including time-dependent boundary values important for simulating the evolution of many physical systems. Because the domain boundaries are not specifically defined in the smoothed boundary method, this method can be applied to almost any geometry as long as it can be defined by the domain parameter. The developed method is thus a very powerful and convenient technique for solving differential equations in complex geometries that are often difficult and time-consuming to structurally mesh. As three-dimensional image-based calculations are increasingly prevalent in scientific and engineering research fields [83, 84, 85], where voxelated data from serial scanning or sectioning are often utilized and are difficult to render as meshes, the smoothed boundary method is expected to be widely employed to simulate and study physics in complex geometries defined by 2D pixelated and 3D voxelated data with a simple process of smoothing the domain boundaries.

Acknowledgements: HCY and KT thank the National Science Foundation for financial support under Grant Nos. 0511232, 0502737, and 0854905. HYC and KT thank the National Science Foundation for financial support under Grant Nos. 0542619 and 0907030. KT also acknowledges the support of the National Science Foundation under Grant No. 0746424. The authors thank John Lowengrub, Axel Voigt, Xiaofan Li, Anton Van der Ven and James Warren for valuable discussions and comments. The authors also thank Scott Barnett and Stuart Adler for providing the experimental 3D microstructures used in the demonstration.

## Appendix A Proof of Convergence for Neumann and Dirichlet Boundary Conditions for the Diffusion Equation

To demonstrate that the smoothed boundary formulated diffusion equation satisfies the assigned Neumann boundary condition (specifying the boundary flux or normal gradient), we use the one-dimensional version of Eq. (6) without loss of generality for cases with higher dimensions. By reorganizing terms and integrating over the interfacial region, we obtain:

 ∫ai+ξ/2ai−ξ/2ψ(∂C∂t−S)dx=ψD∂C∂x∣∣∣ai+ξ/2ai−ξ/2−∫ai+ξ/2ai−ξ/2∣∣∣∂ψ∂x∣∣∣DBNdx, (35)

where is the region of the interface and is the thickness of the interface. Following Refs. [2, 3, 20, 23], we introduce the mean value theorem of integrals, which states that for a continuous function there exists a constant value, , such that , where . By eliminating the second terms on the right-hand sides of Eqs. (6) and (35), the no-flux boundary condition can be imposed (); the resulting equation is similar to those proposed in Refs. [1, 2, 3, 20, 23]. However, here we retain the term to maintain the generality of the method. Therefore, the analysis presented herein leads to an extension of the original method that greatly extends its applicability.

Because the function on the left-hand side of Eq. (35) is continuous and finite within the interfacial region, we can relate its value to the interfacial thickness by , according to the mean value theorem of integrals. Using the conditions that at and at , the first term on the right-hand side of Eq. (35) is written as . Because for or , the bounds of the integral can be extended to and . Therefore, we can rewrite Eq. (35) as

 h0ξ=D∂C∂x∣∣∣ai+ξ/2−∫+∞−∞∣∣∣∂ψ∂x∣∣∣DBNdx, (36)

and by taking the limit of this expression as , we obtain:

 D∂C∂x∣∣∣ai=∫+∞−∞δ(x−ai)DBNdx=DBN∣∣∣ai, (37)

where and when takes the form of a hyperbolic tangent function and is the Dirac delta function. The Dirac delta function has the property that , providing the second equality in Eq. (37). Therefore, Eq. (37) clearly shows that the smoothed boundary method recovers the Neumann boundary condition at the boundary when the thickness of the diffuse boundary approaches zero. This convergence has been observed for both stationary and moving boundaries [20, 32, 33].

To demonstrate the convergence of the solution at the boundaries to the specified boundary value, we again use a one-dimensional version of the smoothed boundary formulated equation. Integrating Eq. (7) over the interfacial region and reorganizing terms, we obtain:

 (38)

Similar to the derivation of Eq. (37), the left-hand side of Eq. (38) is proportional to the interfacial thickness and approaches zero in the limit of . On the right-hand side of Eq. (38), the gradient of approaches the Dirac delta function, , as the interface thickness approaches zero. Therefore, we can reduce Eq. (38) to in the limit . By integrating over the interfacial region again, we obtain:

 −limξ→0h0ξ2D=C∣∣∣ai+ξ/2−∫ai+ξ/2ai−ξ/2BD∂ψ∂xdx, (39)

which gives in the limit of because