- Email: [email protected]
An embedded boundary method for elliptic and parabolic problems with interfaces and application to multi-material systems with phase transitions
Acta Mathematica Scientia 2010,30B(2):499–521 http://actams.wipm.ac.cn
AN EMBEDDED BOUNDARY METHOD FOR ELLIPTIC AND PARABOLIC PROBLEMS WITH INTERFACES AND APPLICATION TO MULTI-MATERIAL SYSTEMS WITH PHASE TRANSITIONS∗ Dedicated to Professor James Glimm on the occasion of his 75th birthday Shuqiang Wang1
Roman Samulyak1,2
Tongfei Guo1
1.Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook, NY 11794, USA 2.Computational Science Center, Brookhaven National Laboratory, Upton, NY 11973, USA E-mail: [email protected]; [email protected]; [email protected]
Abstract The embedded boundary method for solving elliptic and parabolic problems in geometrically complex domains using Cartesian meshes by Johansen and Colella (1998, J. Comput. Phys. 147, 60) has been extended for elliptic and parabolic problems with interior boundaries or interfaces of discontinuities of material properties or solutions. Second order accuracy is achieved in space and time for both stationary and moving interface problems. The method is conservative for elliptic and parabolic problems with fixed interfaces. Based on this method, a front tracking algorithm for the Stefan problem has been developed. The accuracy of the method is measured through comparison with exact solution to a two-dimensional Stefan problem. The algorithm has been used for the study of melting and solidification problems. Key words embedded boundary method; elliptic interface problem; front tracking; Stefan problem 2000 MR Subject Classification
1
65N06
Introduction
The classical Stefan problem provides mathematical foundation for the description of first order phase transition problems such as melting and solidification, vaporization and condensation etc. The numerical simulation of processes described by the Stefan problem has been an active area of research for last decades. Most of methods can be divided into two main classes: phase field models and sharp interface models. ∗ Received
December 1, 2009. This research is supported by the U.S. Department of Energy under Contract No. DE-AC02-98CH10886 and by the State of New York.
500
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
Phase field models and their numerical implementation have been extensively studied by Langer [13], Fix [7], Caginalp [3, 4] etc. In phase field models, a smooth function φ is defined for the identification of particular phase in the computational domain: for instance ρ = 0 in the liquid phase and ρ = 1 in the solid phase. Therefore it implicitly defines the location of the interface. A phase field model formulates an evolution equation for this scalar function coupled to a modified heat transfer equation. Phase field models can easily deal with great geometrical complexity of interfaces. But they often have difficulties in dealing with properties of real fluids such as the surface tension and phase transition regimes such the growth of dendrite structures. The thermodynamical basis of phase field models is still not universally accepted [26]. In the other class of methods, Sullivan et al. [11] performed a finite element simulation of the melting and solidification using a finite element discretization on deforming mesh that dynamically conformed to the phase boundary. They used an accurate numerical algorithm accounting for the interface curvature and obtained a good agreement with experiments on the growth of dendrites. Sethian and Strain [24] used a level set method for the interface tracking and included effects of anisotropic surface tension and the interface kinetics. They used an interface tracking method with curvilinear coordinates and performed simulations of stable and unstable solidification problems. Our method is based on front tracking developed for hyperbolic problems [9] and coupled hyperbolic-elliptic problems such as MHD in the low magnetic Reynolds number approximation [22]. The front tracking based FronTier code [5] is capable of dealing with the propagation and topological changes of a large number of geometrically complex interfaces in two and three dimensional spaces. It has been extensively tested and validated using problems of fluid mixing by Rayleigh-Taylor and Richtmyer-Meshkov instabilities [8, 16], and applied to numerous problems involving multiphase flows. The phase transition algorithm for compressible fluids that accounts for the pressure wave induced mass transfer has already been developed within the framework of front tracking [17]. In this paper, we develop the phase transition algorithm for incompressible flows described by the Stefan problem. This method is close to the work of Juric and Tryggvason [12] who also applied front tracking to the melting and solidification problem but differs by the spatial and temporal discretizations. The quality of sharp interface methods applied to the phase transition problem critically depends on the quality of their numerical schemes for solving parabolic equations in the presence of geometrically complex interfaces. Among the popular methods developed for elliptic and/or parabolic problems with Cartesian meshes, we would like to mention the immersed boundary method by Peskin [21] which uses the discrete delta function on domain boundaries, the immersed interface method by LeVeque and Li [14], the ghost fluid methods by Liu at al. [15], and the integral equation method by Mayo [18], Mckenney, Greengard and Mayo [20]. The series of papers on the Embedded Boundary Method (EBM) by Colella and collaborators [10, 19, 23] proposed a conservative, finite volume discretization of elliptic and parabolic equations in geometrically complex domains using Cartesian meshes. The main feature of the EBM method is the treatment of the solution as cell-centered on a rectangular grid, even when cell centers are outside of the elliptic problem domain. By using linear (in 2D) and bi-linear (in 3D) interpolation of fluxes in cells cut by the interface, well conditioned linear systems and second order accuracy were obtained. In this paper, the embedded boundary method is ex-
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
501
tended to solve the elliptic problems with interior boundaries, parabolic problems that contain fixed or moving interior surfaces (material interfaces) and applied for the solution of the Stefan problem with the method of front tracking. The paper is organized as follows. Section 2 contains the governing system of equations for the Stefan problem. The underlying elliptic and parabolic interface problem and the corresponding embedded boundary algorithms for stationary and moving interface problems are presented in Section 3. Section 4 describes the extension of the numerical algorithms to enable multiple material components within single Cartesian cell and the methods for function interpolation. Section 5 contains examples that validate the algorithms for elliptic and parabolic problems with interfaces as well as the Stefan problem. Section 6 describes an example of the application of proposed algorithms for the simulations of melting processes within the safety study of nuclear fuel rods. We conclude the paper with a summary and conclusions and perspectives for the future work.
2
Governing Equations for the Stefan Problem
Let us consider the melting and solidification processes in multi-material systems. Let Ω be a bounded domain containing multiple materials in the solid and liquid states. We define the material interface as a contact surface between different materials and the phase boundary a contact surface between different phases of the same material. We denote as ∂Ω the collection of interfaces, phase boundaries, and the external boundary of the domain Ω. The equation of the heat transport in Ω \ ∂Ω is [1] ∂ρCP T = ∇ · k∇T + Q, ∂t
(1)
where T is the temperature, ρ is the density, CP is the heat capacity at constant pressure, k is the thermal conductivity, and Q is the external heat source. If equation (1) is applied to a single phase, ρ, CP , and k are typically constants. But since the material properties change with the change of temperature, we do not restrict the algorithm by assuming constant values for these variables. The continuity of temperature and heat flux is satisfied at material interfaces: [T ]material interface = 0, ∂T k = 0, ∂n material interface
(2) (3)
where [.] denotes the jump of quantity across the interface and n is the normal to the interface. At the phase transition boundary, the following conditions are satisfied: [T ]phase
boundary
k
∂T ∂n
= 0,
Tphase
boundary
= TM ,
(4)
= −ρLv,
(5)
phase boundary
where TM is the melting temperature, L is the latent heat, and v is the velocity of the phase boundary. In the last equation, we assumed that the density does not change during the phase transition. It is the standard approximation in the theory of melting and solidification [1].
502
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
In addition, the Gibbs-Thomson condition must be satisfied at the phase boundary. The Gibbs-Thomson effect is the dependence of TM on the interface curvature. Its thermodynamic foundation can be found in [1]. The equation of the Gibbs-Thomson effect is TM = TM,0 −
Kγ , ΔSf
(6)
where TM,0 is the melting temperature at the flat interface, K is the interface curvature, γ is the interfacial tension, and ΔSf is the fusion entropy. Since the change of TM is small at moderate surface curvatures, it is usually negligible during melting. The unsteady solidification which occurs in undercooled liquids (when the liquid temperature is lower then TM ) results in the growth of dendrites of very complex shape. The inclusion of the Gibbs-Thomson effect is necessary for the accurate resolution of the dynamics of unsteady solidification. Finally, various Dirichlet or Neumann or mixed type boundary conditions can be applied on the outside boundaries of the multi-phase domain. We call the heat conduction problem in multi-material system without phase transitions, described by equations (1)–(3), the parabolic interface problem. A numerical algorithm for the parabolic interface problem based on the embedded boundary method is proposed in this paper. Phase boundaries described by the jump condition (5) decouple the problem in a set of parabolic problems with geometrically complex and moving external boundaries (phase boundaries) with the Dirichlet boundary condition Tboundary = TM (see Figure 2).
Fig. 1
Schematic of a multi-material domain with phase transition and its decoupling
into the parabolic interface problem (top right) containing the material interface and the moving boundary (the inner phase boundary with the Dirichlet boundary condition Tphase
boundary
= TM ) and the elliptic problem (bottom right) with the moving outside
boundary (phase boundary with the same Dirichlet boundary condition)
3
Elliptic and Parabolic Interface Problem
In this section, we extend the embedded boundary method for irregular domains [10, 19, 23] to solve the elliptic and parabolic interface problem with interior boundaries.
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
503
The elliptic partial differential equation is: ∇ β∇T = f,
(7)
where β is a continuous function except across the interior boundary and f is also assumed to be a given function which is continuous inside each part of the domain. To close the problem, boundary condition are needed for the exterior and interior boundary. Either Dirichlet or Neumann boundary can be given on the exterior boundary. For the interior boundary, we have the following jump conditions: [T ] = J1 (x),
(8)
∂T = J2 (x), (9) β ∂n where J1 and J2 are given functions of space variables. We generalize the parabolic interface problem introduced in the previous section in the following way. The partial differential equation is ∂T = ∇ β∇T + f, ∂t
(10)
where β and f are continuous functions except across interfaces. Either Dirichlet or Neumann boundary can be given on the exterior boundary. For the interior boundary, we generalize the boundary condition (3) to include the following: [T ] = J1 (x, t),
(11)
∂T = J2 (x, t), (12) β ∂n where J1 and J2 are given functions of space and time. These two jump conditions across the interior interface are similar with those for the elliptic equation (8), (9). The embedded boundary method [10, 19, 23] is a finite volume method for an irregular domain embedded on a Cartesian grid. When the EBM is used to solve the elliptic boundary value problem (without interfaces), the unknowns are defined at the cell centers, even for partial cells which have intersections with the boundary. We retain this main feature of the EBM method in our algorithm for the elliptic/parabolic interface problem: unknowns are defined at centers of regular grid cells. For a full cell (cell containing no external boundaries or interfaces) or a boundary cell (cell cut by the exterior boundary), one unknown is needed at the cell center and a standard finite volume method can be used to setup one algebraic equation as in [10, 19, 23]. For a partial cell (cell cut by the interface), more unknowns are needed to make the discretization of the elliptic/parabolic equation consistent with the two interface jump conditions. Figure 2 depicts the placement of unknowns in a partial cell with interface. The whole cell contains two partial cells (for material components A and B) which are separated by the interface. For each partial cell, there is one unknown defined at the cell center (T1 for the A side, T2 for the B side). In order to satisfy the two jump conditions (8), (9), two more unknowns Ta and Tb are defined at the center of the cell interior boundary (portion of the interface contained within the cell, called cell interface) for the two sides A and B. Therefore, for the partial cell with the interface, four unknowns are defined in total. Four algebraic equations
504
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
can be constructed for this kind of cells using the two jump conditions on the cell interior boundary and the PDE for the two partial cells.
Fig. 2
Placement of unknowns in a cell containing the interface
3.1
Equations for the jump conditions In this section, we describe algorithms for the discretization of the two jump conditions across the interface for the cell (i, j). A schematic of the corresponding stencil and states used in the interpolation method are shown in Figure 3. Two unknowns Ta and Tb are defined at the center of the cell interface for sides A and B, respectively.
Fig. 3
Stencil for the interface unknowns for the jump condition
Let us specify the direction of the normal to the interface as pointing from A to B. For the first jump condition (8), the discretization is simply Tb − Ta = J1 .
(13)
To discretize the second jump condition (9), we need to calculate the normal derivatives of the unknowns. One approach is to use the method given in [10, 19] which fits the quadratic polynomial into a stencil normal to the interface with states in the stencil nodes obtained via interpolation from interior points. Our method of fitting a quadratic polynomial with two variables is explained in detail in Section 4.2 in the later part of the paper. The main idea is to use 1 unknown in cell (i, j) (either cell center unknown or cell interface unknown depending on usage), 2 unknowns in the first layer neighbors of cell (i, j) and 3 unknowns in the second layer neighbors of cell (i, j) which in total give 6 unknowns for the 6 coefficients of the quadratic polynomial. For example, to construct the quadratic polynomial to calculate the flux for side A at the cell interface center where Ta is defined, the 6 unknown candidates are Ta , 2 cell center
No.2
505
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
unknowns Ti+1,j , Ti+1,j−1 in the first layer and 3 cell center unknowns Ti+2,j , Ti+2,j−1 , Ti+2,j−2 in the second layer. Similarly we can fit a quadratic polynomial with two variables for side B. Taking the normal derivative of the fitted polynomials, we obtain the normal derivatives ∂T at the cell interface center for the side A, ∂T ∂n |a , and the side B, ∂n |b . Using the second jump condition (9), we have ∂T ∂T (14) β · − β · = J2 . b ∂n b a ∂n a 3.2 Equations for the cell center unknowns The setup of algebraic equations corresponding to the discretization of the parabolic equation depends on the choice of the time stepping method and the discretization of the elliptic operator. The temporal discretization will be discussed in the next section. Here we describe the discretization of the elliptic operator ∇ β∇T = f.
(15)
For the two cell center unknowns defined at cell (i, j), we use the EBM technique to setup two equations separately. See Figure 4 for the stencil to setup the unknown for side A using the partial cell cdef .
Fig. 4
Stencil for the cell center unknowns
Integrating the equation (15) around the partial cell cdef and using the divergence theorem, we obtain the following expressions: ∇ · β∇T dxdy = β∇T · ndS = f dxdy, ∂(cdef )
cdef
or
β∇T · ndS +
cd
β∇T · ndS +
de
cdef
β∇T · ndS +
ef
β∇T · ndS =
fc
f dxdy, cdef
which is lcd · Fluxcd + lde · Fluxde + lef · Fluxef + lf c · Fluxf c =
f dxdy
(16)
cdef
where lmn is the length between m and n. Therefore we only need to calculate the flux across the cell edges or cell interface. For Fluxcd , a 2nd order derivative is calculated by using a T −Ti,j T −Ti+1,j linear interpolation of i,j−1 and i+1,j−1 to the center of cd (see [10]). For Fluxde , x x
506 we simply use Ti,j+1 −Ti,j x
ACTA MATHEMATICA SCIENTIA Ti+1,j −Ti,j to y Ti+1,j+1 −Ti+1,j x
Vol.30 Ser.B
calculate the derivative. For Fluxef , a linear interpolation of
and to the center of ef is used. Fluxf c , is calculated using ∂T ∂n |a used in (14). Note that we need to multiply the derivatives calculated with β at the corresponding edge or cell interface center to get the flux. In the same way, we calculate fluxes for the other partial cells. More details can be found in [10]. It is also worth noting that the two unknowns defined at the cell interface can be represented by the cell centers unknowns by solving the two discretized jump conditions (13), (14). This can result in a system of equations with only unknowns defined at cell centers. 3.3 Implicit time stepping methods The embedded boundary method is used with the implicit Runge-Kutta method (TGA, [25]) for solving the time dependent parabolic initial boundary value problem in ([19, 23]). Here we apply this time stepping technique to solve the parabolic interface problem. As in [19], in order to solve the following time-dependent initial boundary value problem: dU = LhI (t)U (t) + F (t), dt
(17)
U (0) = U 0 ,
(18)
where LhI (t) is the time-dependent inhomogeneous operator, we use the two step implicit RungeKutta method (TGA, [25]): U n+1 = (I − μ1 Lh1 (tnew ))−1 (I − μ2 LhI (tint ))−1 ·[(I +
μ3 LhI (told ))U n
+ (I +
μ4 LhH )f (tavg )t],
(19) (20)
where LhH (t) is LhI (t) with homogeneous boundary conditions and told = nt, tnew = (n + 1)t = told + μ1 + μ2 + μ3 , tint = tnew − μ1 = told + μ2 + μ3 , tavg = (told + tnew )/2 = told + μ1 + μ2 + μ4 , √ a − a2 − 4a + 2 t, μ1 = √ 2 a + a2 − 4a + 2 μ2 = t, 2 μ3 = (1 − a)t, 1 μ4 = ( − a)t. 2
√ To get a L0 -stable method with real arithmetic, it is required that 1/2 < a < 2 − 2. In this case, μ1 > 0, μ2 > 0, μ3 > 0 > μ4 . Therefore, told < tavg < tint < tnew . When a = 1/2, we obtain the Crank-Nicolson method with μ1 = μ4 = 0 and tint = tnew . It should be noted that the homogeneous boundary conditions of the operator LhH corresponding to the interior boundary jump condition (8), (9) are to set the jumps to be zero [T ] = 0, ∂T = 0. β ∂n
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
507
For the boundary conditions corresponding to the exterior boundary, we set T = 0 for Dirichlet boundary or set ∂T ∂n = 0 for Neumann boundary. Equation (19) are solved in two step: (I − μ2 LhI (tint ))U tint = (I + μ3 LhI (told ))U n + (I + μ4 LhH )f (tavg )t
(21)
(I − μ1 Lh1 (tnew ))U n+1 = U tint
(22)
Each equation is discretized using the previously developed embedded boundary method. For example, equation (21) is the following: U tint − μ2 ∇ β(tint )∇U tint = U n + μ3 ∇ β(tint )∇U n + t[f (tavg ) + μ4 ∇ β∇f (tavg )]. 3.4
Moving boundary problem
In this Section, we generalize the EBM algorithm for problems with moving interfaces. Our motivation is to develop an accurate numerical methods for the sharp interface classical Stefan problem. For the moving boundary problem, the EBM method for elliptic problems with moving exterior boundaries ([19, 23]) extrapolates the boundary condition located at the position of the intermediate step tint to the position of the new boundary at time tnew using Taylor series. It is impractical to apply similar method for interior boundaries with the jump conditions (8), (9). Instead we solve the first step of the implicit Runge-Kutta (21) using the intermediate boundary at time tint and solve the second step (22) using the new boundary at time tnew . The uncovered states for the intermediate and new interface are obtained via extrapolation. Our extrapolation method is different from those in ([19, 23]): the extrapolation is performed through the center of the nearest cell interface by fitting a quadratic polynomial with two variables which is the same as the fitting used in previous section to calculate the cell interface flux to discretize the parabolic equation. Figure 5 shows the extrapolation method. Both cell center states and cell interface states may need to be extrapolated when the interface is moved. The extrapolation method is the same for both cases. Note that we do not do the extrapolation in the normal to the interface direction and it is not necessary. In Figure 5, to extrapolate state to the new cell interface at a for side A, we first find the nearest old cell interface d, then find the candidates for fitting in the direction along line ad.
Fig. 5
Extrapolation method for moving boundary
508
4 4.1
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
Multiple Components and Interpolation Method Algorithms and coding for multiple components
For numerous practical applications, the assumption of at most two material components in a single cell is not correct even for sufficiently refined meshes. Figure 6 shows a computational domain with three material components that cross at a single point. Unless the intersection point is precisely on the mesh line, there is one cell containing three different components. In order to deal with this set of problems, we would like to extend our original method but preserve the central idea of the method that unknowns are defined at cell centers (even for irregular cells) and on cell interface centers. We first need to clarify the possible cases of cells in problems with multiple components with assumption that there is at most one crossing by the exterior or interior boundary on each cell edge. We also assume that no crossing occur at the nodes of the Cartesian grid. We already have external cells which is located outside the computational domain, internal cells without exterior or interior boundary, boundary cells containing exterior boundary, partial cells containing interior boundary with at most 2 components. We need to deal with cells with 3 components as in Figure 7 and cells 4 components as in Figure 8. These two cases will still exist even under mesh refinement. There are also cases which can be dealt with using mesh refinement as in Figure 9 and Figure 10. However, we also deal with them to make the code to be more robust.
Fig. 6
Computational domain with
3
3
1
2
Fig. 7
Cell with 3 components
three components A, B, C 3
4
2
1
1
2
1
2
Fig. 8
Cell with 4 components
Fig. 9
Cell with 2 components
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
509
The algorithm to setup the unknown matrix for solving elliptic equation with multiple components is the following: For each Cartesian cell, do the following: 1. for each cell interface, setup the equations related to the unknowns on the cell interface and the two components. 2. for each cell edge, setup the equations. 3. for each cell interior, setup the interior equation corresponding to that component. To setup the unknown matrix, we need to enumerate the cell interface. Each cell interface is identified by a identification number (ID) consisting of three numbers (a, b, c): a and b are the two ordered components for the two side of the cell interface; c is an integer starting from 0 which is needed when two cell interface has the same two components for cases as in Figure 9. The information stored at the cell center for a particular component can be identified by the component number. Each cell contains many information including number of cell unknowns, cell components, unknowns index, cell areas for each components, pointer to states, number of cell interface unknowns, cell interface length, cell interface normal, cell interface center. It is impractical to use a single data structure to store simple cells and multiple component cells (which need to store much more information). To save computational memory, two different data structure are used: one for cells holding at most 2 components and one for cells holding more complex cells as in Figures 7–10. More memory can be saved to use structures for cells with 1 component and cells with 2 components. It is notable that cell containing multiple components are in the order of O(1) if the exterior and interior boundary is not very complex. 4.2
Interpolation method
To deal with multiple components in a single cell and improve the robustness of the code, the interpolation method need to be modified. As in Figure 11, our task is to get a bilinear polynomial using three unknowns or biquadratic polynomial using six unknowns located at different locations. Taking the derivative of these polynomial in the normal direction of the cell interface, we get the derivative which is needed to setup the matrix for the elliptic equations. If quadratic interpolation is not possible, linear interpolation is used which is always possible under our assumption that there is at most one crossing on each cell edge by the exterior boundary or interior interface.
Fig. 10
Cell with 3 components
Fig. 11
Interpolation method
510
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
When selecting unknowns to do interpolation for one of the component, at most one of the unknowns of a cell (located at cell center or cell interface) with the same component is selected. Otherwise, their positions might be too close, leading to the singularity of the matrix for the interpolation coefficients. It is also required that at most one cell interface unknown (say the current cell (i, j)) is selected for an interpolation stencil. We classify the cells into different group of neighbors. Zeroth layer neighbor: for the cell (i, j), the zeroth layer neighbor is simply the current cell (i, j). First layer neighbor: for the cell (i, j), its first layer neighbors are cells with the cell index (m, n) satisfying max(|m − i|, |n − j|) = 1. Thus the first layer cells include (i − 1, j − 1), (i, j − 1), (j + 1, j − 1), (i − 1, j), (i, j), (i + 1, j), (i − 1, j + 1), (i, j + 1), (i + 1, j + 1). Second layer neighbor: the second layer neighbors for cell (i, j) are similarly defined as cells with cell index (m, n) satisfying max(|m − i|, |n − j|) = 2. Consecutive cells: the consecutive cells in this paper denote the consecutive cells in a given layer. For example, consider the first layer of cell (i, j): cell (i + 1, j) is consecutive with cell (i + 1, j + 1) and (i + 1, j − 1) in the first layer; it is not consecutive with any other cells in the first layer. Similarly, consider the second layer of cell (i, j): cell (i + 2, j + 2) is only consecutive with cell (i + 1, j + 2) and (i + 2, j + 1); it is not consecutive with any other cells in the second layer. Now we can select unknowns to do interpolation. Considering cell (i, j) with its neighbors, the cells in Figure 11 consist of cell (i, j), its first and second layer neighbors. To do bilinear interpolation, we select from cell (i, j) and two consecutive cells from its first layer neighbors, using only one unknown from each cell resulting to 3 unknowns in total for the three coefficients of linear interpolation. To do biquadratic interpolation, we select another three consecutive cells from its second layer of neighbors besides the three unknowns for the first order interpolation resulting to 6 unknowns in total for the six coefficients of biquadratic interpolation. Only cell center unknowns are used in the interpolation except the current cell (i, j) for which either the cell center or the cell interface unknown is selected depending on what is needed. When selecting candidate unknowns for interpolation, it is preferred to first check along the cell interface normal direction. In Figure 11, a candidate stencil for linear interpolation of solution with component B includes cell (i, j), (i + 1, j), (i + 1, j + 1) for linear interpolation. Add cell (i + 2, j + 1), (i + 2, j + 2), (i + 1, j + 2) for quadratic interpolation. We have the following three observations for boundary or partial cells. Proposition 1 For a cell containing component a at one of its four corners, it is always possible to find three unknowns with the same components to do bilinear interpolation. Proof Denote the current cell index as (i, j). Since it contains component a, there exists one cell interface with one side of its component as a. Thus we can pick the cell unknown or cell interface unknown of cell (i, j). Suppose that the down-right corner of the cell (i, j) has
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
511
component a. Then it must be the case that cell (i, j − 1) contains component a at its up-right corner. Similarly, cell (i + 1, j − 1), (i + 1, j) both contain the same components. Therefore, two more unknowns can be selected from either {(i, j − 1), (i + 1, j − 1)}, or {(i + 1, j − 1), (i + 1, j)}. Proposition 2 The matrix for the 3 selected unknowns for linear interpolation are always solvable (non-singular). Proof For cell (i, j), suppose that an interpolation stencil is (i, j) , (i + 1, j), (i + 1, j + 1). The interpolation algorithm will pick cell center unknown from the last two cell and either cell center or cell interface unknown from the first cell. It is easy to know that any point in cell (i, j) can not be on the same line passing through the center of the last two cell. Thus the matrix for this particular stencil is non-singular. Proposition 3 The matrix for the 6 selected unknowns for quadratic interpolation are always solvable (non-singular). Proof This statement can be verified by the same reasoning as above.
5
Numerical Results
In this section, we validate the proposed numerical algorithms by solving elliptic and parabolic interface problems with interior boundaries and comparing numerical solutions with exact analytical ones. We also discuss the accuracy and convergence of numerical schemes. 5.1 Example 1 elliptic interface problem We solve the elliptic interface problem inside [−1, 1] × [−1, 1]. The interior interface in Figure 12 is described by (x − 0.1)2 + (y − 0.1)2 = 0.5 + 0.2 × cos(θ), where θ is the angle between (x − 0.1, y − 0.1) and x-coordinate.
Fig. 12
Interface for example 1 elliptic interface problem
Denote the part inside the interface as A and the part out of the interface as B. We use the following coefficient ⎧ ⎨ 1 (x, y) ∈ A, β(x, y) = ⎩ 1 (x, y) ∈ B, 10
512
ACTA MATHEMATICA SCIENTIA
and exact solution T (x, y) =
⎧ ⎨ (x2 + y 2 ) 32
(x, y) ∈ A,
⎩e
(x, y) ∈ B.
x2 +y 2
Vol.30 Ser.B
The approximate solution is given in Figure 13. Table 1 shows the max errors under mesh refinement. The method is 2nd order accurate. Table 1
Max errors for example 1 elliptic interface problem
Fig. 13
5.2
Mesh Size
Max Error
Ratio
40 × 40
1.249277e-02
80 × 80
3.290830e-03
3.80
160 × 160
8.442124e-04
3.90
320 × 320
2.137630e-04
3.95
640 × 640
5.378009e-05
3.97
Approximate solution for example 1 (160 × 160)
Example 2 Elliptic interface problem with multiple components The computational domain is [−1, 1] × [−1, 1]. It is partitioned in to three different part A = {(x, y)|x < 0.3, y < 0.3}, B = {(x, y)|x ≥ 0.3, y < 0.3}, C = {(x, y)|y ≥ 0.3}. Figure 14 shows the computational domain with interior interfaces.
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
Fig. 14
Interface for Example 2 elliptic interface problem with multiple components Table 2 Max errors with mesh refinement for Example 2 elliptic interface problem with multiple components
Fig. 15
Mesh Size
Max Error
Ratio
20 × 20
1.560965e-02
40 × 40
3.942480e-03
3.96
80 × 80
9.901261e-04
3.98
160 × 160
2.480690e-04
3.99
320 × 320
6.208303e-05
4.00
Approximate solution for Example 2 elliptic interface problem with multiple components (160 × 160)
We use the following coefficient
⎧ ⎪ 1 ⎪ ⎪ ⎨ 1 β(x, y) = ⎪ 10 ⎪ ⎪ ⎩3
(x, y) ∈ A, (x, y) ∈ B, (x, y) ∈ C,
513
514
ACTA MATHEMATICA SCIENTIA
and exact solution
⎧ 2 2 32 ⎪ ⎪ ⎨ (x + y ) T (x, y) =
⎪ ⎪ ⎩
e
(x, y) ∈ A,
x2 +y 2
2(x2 + y 2 )
Vol.30 Ser.B
(x, y) ∈ B, 3 2
(x, y) ∈ C.
Table 2 shows the max errors and Figure 15 shows the approximate solution. 5.3 Example 3 elliptic interface problem in 3D The computational domain is [−1, 1] × [−1, 1] × [−1, 1]. The interface is a circle with center at (0, 0, 0) with radius r = 0.5. Let A be the region inside the interface and B the region outside the interface. We use the following as the coefficient: ⎧ ⎨ 1 (x, y) ∈ A, β(x, y) = ⎩ 10 (x, y) ∈ B, and exact solution: T (x, y) =
⎧ 3 ⎪ ⎨ (x2 + y 2 ) 2 2
2
(x, y) ∈ A,
3 2
(x + y ) 1 1 ⎪ ⎩ + (1 − ) × 10 10 8
(x, y) ∈ B.
Currently, the function interpolation method used for the 3D code is similar with those in [23]. Table 3 shows the max errors with mesh refinement. Table 3
Max error for Example 3 elliptic interface problem for 3D Mesh Size
Max Error
Ratio
10 × 10
1.077035e-01
20 × 20
3.451512e-02
3.12
40 × 40
9.139040e-03
3.78
80 × 80
2.362821e-03
3.87
160 × 160
6.045999e-04
3.91
5.4
Example 4 parabolic interface problem with fixed interface For this example, we solve the parabolic interface problem on a fixed interface. The computational domain is [−1, 1] × [−1, 1]. The interface is a circle with center at [−0.1, −0.1] with radius r = 0.5. Let A be the region inside the interface and B the region outside the interface. We use the following coefficient: ⎧ ⎨ 1 (x, y) ∈ A, β(x, y) = ⎩ 10 (x, y) ∈ B, and the exact solution: ⎧ x2 +y 2 ⎪ 4 exp(− 5(t+1) ) ⎪ ⎪ ⎪ ⎨ 5π(t + 1) T (x, y) = x2 +y 2 ⎪ ⎪ exp(− 2(t+1) ) ⎪ ⎪ ⎩ 2π(t + 1)
(x, y) ∈ A, (x, y) ∈ B,
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
515
where the formula for T (x, y), x ∈ A is picked from [23]. The heat equation is solved until t = 0.5. Table 4 shows the max errors at time t = 0.5. Figure 16 shows the approximate solution at time t = 0.5 with mesh size 160 × 160. The test confirms the second order accuracy of the method. Table 4
Fig. 16
5.5
Maximum errors at time t = 0.5 with mesh refinement for Example 4 parabolic interface problem with fixed interface, t = 0.9 h Mesh Size
Max Error
Ratio
20 × 20
3.709544e-04
40 × 40
8.830904e-05
4.20
80 × 80
2.257112e-05
3.91
160 × 160
5.705129e-06
3.96
Approximate solution for Example 4 parabolic interface problem with fixed interface
Example 5 parabolic interface problem with prescribed boundary motion In the second numerical test, we solve the parabolic interface problem with interface that moves with a given velocity. The computational domain is [−1, 1] × [−1, 1]. The initial interface is a circle centered at [−0.25, −0.25] with radius r = 0.4. The motion of the circle is defined simply by dx = 0.5, dt dy = 0.5. dt Other conditions are the same as Example 4. Thus the motion is simply a rigid motion. Table 5 shows the max errors under the mesh refinement. Figure 17 shows the approximate solution at time t = 0.2 with mesh size 320 × 320. The time step in the simulation was restricted by t h = 0.1. The increase of the time step t to h = 0.9 reduces the order of accuracy to approximately 1.2. To retain the second order accuracy, the time step must be lower than the stability limit t t max v ,β < 1. h h
516
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
The last expressions assumes en explicit method for the parabolic equation. Since the velocity of the moving interface is v = 0.5, t h < 0.9 satisfies the stability condition v
t < 1. h
However the increase of the time step to t h = 0.9 results in the loss of accuracy caused by the extrapolation of states to the location of cells swept during the interface motion. Table 5 Maximum errors at time t = 0.2 with mesh refinement for Example 5 parabolic interface problem with prescribed boundary motion, t = 0.1 h
Fig. 17
Mesh Size
Max Error
Ratio
40 × 40
4.103197e-05
80 × 80
1.008246e-05
4.07
160 × 160
2.895398e-06
3.49
320 × 320
6.377383e-07
4.54
Approximate solution at time t = 0.2 for Example 5 parabolic interface problem with prescribed boundary motion.
t h
= 0.1
5.6
Example 6 Stefan problem In this numerical example, we validate the Stefan problem algorithm using equations that possess an exact analytical solution. Consider a square domain [−1, 1] × [−1, 1] containing two phases (solid and liquid) of the same substance. The phases are separated by a circular interface that moves as follows 1 R0 (t) = 3 . − 4 (t + 4)4 + 23 4 The liquid is inside the circle. The temperature of the liquid satisfies the following analytical formula T (r) = (R0 (t)3 − r3 )(t + 1)3 ,
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
517
The temperature is equal to 0 on the interface. The analytical formula for the solid phase temperature is simply 0. The PDE is Tt = T + f, where f can be obtained by substituting the exact solution. The starting time is t0 = 0.4 and the final time T = 0.45. Figure 18 shows the maximum interface position errors and Figure 19 shows the maximum solution errors. Both plots demonstrate the 2nd order accuracy of the algorithm. −4 log(max intfc error) ~ log(dx) y=2.0153*log(dx) + 0.3868 −5
log(max intfc error)
−6
−7
−8
−9
−10 −5.5
−5
−4.5
−4
−3.5
−3
−2.5
−2
log(dx)
Fig. 18
Maximum interface position error for Example 6 Stefan problem at time T = 0.45 −3 log(max solution error) ~ log(dx) y=1.9153*log(dx) + 0.5717 −4
log(max solution error)
−5
−6
−7
−8
−9
−10 −5.5
−5
−4.5
−4
−3.5
−3
−2.5
−2
log(dx)
Fig. 19
6
Maximum solution error for Example 6 Stefan problem at time T = 0.45
Application
In this section, we illustrate the application of our algorithm to the study of accident scenarios in nuclear fuel rods. This study is being performed within multi-institutional project
518
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
on multiscale simulations of Generation-IV Sodium Fast Reactors [2] focusing on reactor fuel and core transient response under beyond-design and accident conditions. A nuclear fuel rod contains 8 mm diameter pellets of metallic or oxide uranium fuel in a stainless still cladding of 0.5 mm thickness. The total length of the fuel rod is about 2.5 m. A narrow (0.1 mm), irregular gap between the fuel and the cladding is used to transport fission gasses from the fuel to the gas plenum on the top of the reactor core. Fuel rods are assembled in hexagonal structures and places into liquid sodium coolant that circulates in the reactor core and transports the thermal energy from fuel rods to the heat exchangers. The cross-section of the fuel rod is shown schematically in Figure 20.
Fig. 20
Schematic of the nuclear fuel rod cross-section: (a) is the fuel pellet,
(b) is the gas gap, (c) is the stainless steel cladding, and (d) is the liquid sodium.
The temperature distribution across the fuel rod at normal operating conditions is shown in Figure 22. The power production of about 60 kW/m in the fuel rod is balanced by the heat removal by the sodium coolant so that the temperature is in the steady state. The picture shows rapid changes of the temperature gradients across the fuel rod due to sharp changes of the heat conductivity in the fuel, gap, and cladding. After years of operation, the surface of the fuel erodes and comes to a point-wise contact with the clad as shown in Figure 21. The interior of the gas gap is not resolved on the mesh level and an empirical formula for the temperature jump across the gap [6] is applied for the heat transport calculations. For the related problem of the transport of fission gases in the gas gap, the flow in porous media equations are solved.
Fig. 21
Shape of the gas gap between the fuel and cladding
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
3500
Melting Point Overheating State Normal State
3000
Temperature (K)
519
2500 2000 1500 1000 500 0
Fig. 22
0.1
0.2
0.3
0.4 0.5 Radius (cm)
0.6
0.7
0.8
0.9
Temperature across the fuel rod at normal operating conditions and transient overheating
During either transient overheating or loss of coolant accidents, the heat transfer balance is changed and the temperature in the fuel rod increases. It can cause melting of the fuel rod and even stainless steel cladding. However because of the change of cladding properties with the temperature increase, the cladding usually fails before melting, causing the ejection of molten fuel and fission gases into the coolant reservoir. Figures 23 and 22 depict the increase of temperature and melting of the fuel in the transient overheating accident.
Fig. 23
Temperature distribution across the fuel rod during transient overheating and the surface of the molten fuel
7
Conclusions
In this paper, we generalized the embedded boundary algorithm for solving elliptic and parabolic problems in geometrically complex domains [10] to include irregular interior boundaries or interfaces. The algorithm was validated by solving problems that possess exact analytical solutions and the second order accuracy in both space and time was confirmed. Based on this method, we developed an algorithm for solving the classical Stefan problem that describes first order phase transitions such as melting and solidification. The Stefan problem algorithm was also validated using an exact analytical solution. The combination of elliptic
520
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
and parabolic algorithms makes it possible to simulate the heat transfer and phase transitions in multi-material systems containing geometrically complex interfaces between different materials and phase boundaries or surfaces between different phases of the same material. The explicit tracking of material interfaces removes all restrictions on the magnitude of discontinuities of physical properties and temperature gradients across interfaces and phase boundaries and is especially beneficial for the study of materials at extreme conditions. We illustrated the use of the method in the study of accident scenarios in nuclear fuel rods in Generation-IV Sodium Fast Reactors. The embedded boundary method was used to simulate the temperature increase in the fuel assembly and melting of uranium fuel during the transient overpower and loss of coolant accidents. The algorithm was implemented in FronTier, a multiphysics code for the simulation of multiphase flows using the method of front tracking. In the future, the embedded boundary method will also be implemented in the incompressible hydrodynamics and magnetohydrodynamics components of the FronTier code. Acknowledgments The authors would like to thank James Glimm and Xiaolin Li for stimulating discussions. This manuscript has been authored in part by Brookhaven Science Associates, LLC, under Contract No. DE-AC02-98CH10886 with the U.S. Department of Energy. The United States Government retains, and the publisher, by accepting the article for publication, acknowledges, a world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for the United States Government purpose. This research utilized resources at the New York Center for Computational Sciences at Stony Brook University/Brookhaven National Laboratory. References [1] Alexiades V, Solomon A D. Mathematical Modeling of Melting and Freezing Processes. McGraw-Hill, 1993 [2] Bolotnov I A, Behafarid F, Shaver D, Antal S P, Jansen K E, Samulyak R, Podowski M Z. Interaction of computational tools for multiscale multiphysics simulation of generation-iv reactors//Proceedings of International Congress on Advances in Nuclear Power Plants, 2010. To be published [3] Caginalp G, Fife P. Higher-order phase field models and detailed anisotropy. Phys Rev B, 1986, 34: 4940–4943 [4] Caginalp G. Stefan and Hele-Shaw type models as asymptotic limits of the phase field equations. Phys Rev A, 1989, 39: 5887–5896 [5] Du Jian, Fix Brian, Glimm James, Jia Xicheng, Li Xiaolin, Li Yunhua, Wu Lingling. A simple package for front tracking. J Comput Phys, 2006, 213: 613–628 [6] Duderstadt J, Hamilton L. Nuclear Reactor Analysis. 2nd ed. John Wiley & Sons, Inc, 1976 [7] Fix G.//Fasano A, Promiceiro M, ed. Free Boundary Problems. London: Pitman, 1983: 580 [8] George E, Glimm J, Li X L, Li Y H, Liu X F. The influence of scale-breaking phenomena on turbulent mixing rates. Phys Rev E, 2006, 73: 1–5 [9] Glimm J, Grove J W, Li X -L, Shyue K -M, Zhang Q, Zeng Y. Three dimensional front tracking. SIAM J Sci Comp, 1998, 19: 703–727 [10] Johansen H, Colella P. A cartesian grid embedding boundary method for poisson’s equation on irregular domains. J Comput Phys, 1998, 147: 60–85 [11] Sullivan J M, Lynch D R, O’Neill K. Finite-element simulation of planare instabilities during solidification of an undercooled melt. J Comput Phys, 1987, 69: 81–111
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
521
[12] Juric D, Tryggvason G. A front tracking method for dendritic solidi-fication. J Comput Phys, 1996, 123: 127–148 [13] Langer.//Grinstein G, Mazenko G, ed. Directions in Condensed Matter Physics. Singapore: Worls Scientific, 1986: 164 [14] LeVeque R J, Li Z L. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM J Numer Anal, 1994, 31: 1019–1044 [15] Liu X -D, Fedkiw R, Kang M. A boundary condition capturing method for poisson’s equation on irregular domains. J Comput Phys, 2000, 160: 151–178 [16] Liu X F, Li Y H, Glimm J, Li X L. A front tracking algorithm for limited mass diffusion. J Comput Phys, 2007, 222: 644–653 [17] Lu T, Xu Z L, Samulyak R, Glimm J, Ji X M. Dynamic phase boundaries for compressible fluids. SIAM J Sci Comp, 2008, 30: 895–815 [18] Mayo A. The fast solution of poisson and the biharmonic equations on irregular regions. SIAM J Numer Anal, 1984, 21: 285–299 [19] McCorquodale P, Colella P, Johansen H. A cartesian grid embedded boundary method for the heat equation on irregular domains. J Comput Phys, 2001, 173: 620–635 [20] Mckenney A, Mayo A, Greengard L. A fast poisson solver for complex geometries. J Comput Phys, 1995, 118: 348–355 [21] Perskin Charles S. The immersed boundary method. Acta Numerica, 2002, 11: 479–517 [22] Samulyak R, Du J, Glimm J, Xu Z. A numerical algorithm for MHD of free surface flows at low magnetic reynolds numbers. J Comput Phys, 2007, 226: 1532–1546 [23] Schwartz P, Barad M, Colella P, Ligocki T. A cartesian grid embedded boundary method for the heat equation and Poisson’s equation in three dimensions. J Comput Phys, 2006, 211: 531–550 [24] Strain J, Sethian J. Crystal Growth and dendritic solidification. J Comput Phys, 1992, 98: 231–253 [25] Twizell E H, Gumel A B, Arigu M A. Second-order, l0-stable methods for the heat equation with timedependent boundary conditions. Adv Comput Math, 1996, 6: 333 [26] Wheeler A A.//Hurle D T J, ed. Handbook of Crystal Growth, Volume 1B. Amsterdam: North-Holland, 1993: 783
AN EMBEDDED BOUNDARY METHOD FOR ELLIPTIC AND PARABOLIC PROBLEMS WITH INTERFACES AND APPLICATION TO MULTI-MATERIAL SYSTEMS WITH PHASE TRANSITIONS∗ Dedicated to Professor James Glimm on the occasion of his 75th birthday Shuqiang Wang1
Roman Samulyak1,2
Tongfei Guo1
1.Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook, NY 11794, USA 2.Computational Science Center, Brookhaven National Laboratory, Upton, NY 11973, USA E-mail: [email protected]; [email protected]; [email protected]
Abstract The embedded boundary method for solving elliptic and parabolic problems in geometrically complex domains using Cartesian meshes by Johansen and Colella (1998, J. Comput. Phys. 147, 60) has been extended for elliptic and parabolic problems with interior boundaries or interfaces of discontinuities of material properties or solutions. Second order accuracy is achieved in space and time for both stationary and moving interface problems. The method is conservative for elliptic and parabolic problems with fixed interfaces. Based on this method, a front tracking algorithm for the Stefan problem has been developed. The accuracy of the method is measured through comparison with exact solution to a two-dimensional Stefan problem. The algorithm has been used for the study of melting and solidification problems. Key words embedded boundary method; elliptic interface problem; front tracking; Stefan problem 2000 MR Subject Classification
1
65N06
Introduction
The classical Stefan problem provides mathematical foundation for the description of first order phase transition problems such as melting and solidification, vaporization and condensation etc. The numerical simulation of processes described by the Stefan problem has been an active area of research for last decades. Most of methods can be divided into two main classes: phase field models and sharp interface models. ∗ Received
December 1, 2009. This research is supported by the U.S. Department of Energy under Contract No. DE-AC02-98CH10886 and by the State of New York.
500
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
Phase field models and their numerical implementation have been extensively studied by Langer [13], Fix [7], Caginalp [3, 4] etc. In phase field models, a smooth function φ is defined for the identification of particular phase in the computational domain: for instance ρ = 0 in the liquid phase and ρ = 1 in the solid phase. Therefore it implicitly defines the location of the interface. A phase field model formulates an evolution equation for this scalar function coupled to a modified heat transfer equation. Phase field models can easily deal with great geometrical complexity of interfaces. But they often have difficulties in dealing with properties of real fluids such as the surface tension and phase transition regimes such the growth of dendrite structures. The thermodynamical basis of phase field models is still not universally accepted [26]. In the other class of methods, Sullivan et al. [11] performed a finite element simulation of the melting and solidification using a finite element discretization on deforming mesh that dynamically conformed to the phase boundary. They used an accurate numerical algorithm accounting for the interface curvature and obtained a good agreement with experiments on the growth of dendrites. Sethian and Strain [24] used a level set method for the interface tracking and included effects of anisotropic surface tension and the interface kinetics. They used an interface tracking method with curvilinear coordinates and performed simulations of stable and unstable solidification problems. Our method is based on front tracking developed for hyperbolic problems [9] and coupled hyperbolic-elliptic problems such as MHD in the low magnetic Reynolds number approximation [22]. The front tracking based FronTier code [5] is capable of dealing with the propagation and topological changes of a large number of geometrically complex interfaces in two and three dimensional spaces. It has been extensively tested and validated using problems of fluid mixing by Rayleigh-Taylor and Richtmyer-Meshkov instabilities [8, 16], and applied to numerous problems involving multiphase flows. The phase transition algorithm for compressible fluids that accounts for the pressure wave induced mass transfer has already been developed within the framework of front tracking [17]. In this paper, we develop the phase transition algorithm for incompressible flows described by the Stefan problem. This method is close to the work of Juric and Tryggvason [12] who also applied front tracking to the melting and solidification problem but differs by the spatial and temporal discretizations. The quality of sharp interface methods applied to the phase transition problem critically depends on the quality of their numerical schemes for solving parabolic equations in the presence of geometrically complex interfaces. Among the popular methods developed for elliptic and/or parabolic problems with Cartesian meshes, we would like to mention the immersed boundary method by Peskin [21] which uses the discrete delta function on domain boundaries, the immersed interface method by LeVeque and Li [14], the ghost fluid methods by Liu at al. [15], and the integral equation method by Mayo [18], Mckenney, Greengard and Mayo [20]. The series of papers on the Embedded Boundary Method (EBM) by Colella and collaborators [10, 19, 23] proposed a conservative, finite volume discretization of elliptic and parabolic equations in geometrically complex domains using Cartesian meshes. The main feature of the EBM method is the treatment of the solution as cell-centered on a rectangular grid, even when cell centers are outside of the elliptic problem domain. By using linear (in 2D) and bi-linear (in 3D) interpolation of fluxes in cells cut by the interface, well conditioned linear systems and second order accuracy were obtained. In this paper, the embedded boundary method is ex-
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
501
tended to solve the elliptic problems with interior boundaries, parabolic problems that contain fixed or moving interior surfaces (material interfaces) and applied for the solution of the Stefan problem with the method of front tracking. The paper is organized as follows. Section 2 contains the governing system of equations for the Stefan problem. The underlying elliptic and parabolic interface problem and the corresponding embedded boundary algorithms for stationary and moving interface problems are presented in Section 3. Section 4 describes the extension of the numerical algorithms to enable multiple material components within single Cartesian cell and the methods for function interpolation. Section 5 contains examples that validate the algorithms for elliptic and parabolic problems with interfaces as well as the Stefan problem. Section 6 describes an example of the application of proposed algorithms for the simulations of melting processes within the safety study of nuclear fuel rods. We conclude the paper with a summary and conclusions and perspectives for the future work.
2
Governing Equations for the Stefan Problem
Let us consider the melting and solidification processes in multi-material systems. Let Ω be a bounded domain containing multiple materials in the solid and liquid states. We define the material interface as a contact surface between different materials and the phase boundary a contact surface between different phases of the same material. We denote as ∂Ω the collection of interfaces, phase boundaries, and the external boundary of the domain Ω. The equation of the heat transport in Ω \ ∂Ω is [1] ∂ρCP T = ∇ · k∇T + Q, ∂t
(1)
where T is the temperature, ρ is the density, CP is the heat capacity at constant pressure, k is the thermal conductivity, and Q is the external heat source. If equation (1) is applied to a single phase, ρ, CP , and k are typically constants. But since the material properties change with the change of temperature, we do not restrict the algorithm by assuming constant values for these variables. The continuity of temperature and heat flux is satisfied at material interfaces: [T ]material interface = 0, ∂T k = 0, ∂n material interface
(2) (3)
where [.] denotes the jump of quantity across the interface and n is the normal to the interface. At the phase transition boundary, the following conditions are satisfied: [T ]phase
boundary
k
∂T ∂n
= 0,
Tphase
boundary
= TM ,
(4)
= −ρLv,
(5)
phase boundary
where TM is the melting temperature, L is the latent heat, and v is the velocity of the phase boundary. In the last equation, we assumed that the density does not change during the phase transition. It is the standard approximation in the theory of melting and solidification [1].
502
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
In addition, the Gibbs-Thomson condition must be satisfied at the phase boundary. The Gibbs-Thomson effect is the dependence of TM on the interface curvature. Its thermodynamic foundation can be found in [1]. The equation of the Gibbs-Thomson effect is TM = TM,0 −
Kγ , ΔSf
(6)
where TM,0 is the melting temperature at the flat interface, K is the interface curvature, γ is the interfacial tension, and ΔSf is the fusion entropy. Since the change of TM is small at moderate surface curvatures, it is usually negligible during melting. The unsteady solidification which occurs in undercooled liquids (when the liquid temperature is lower then TM ) results in the growth of dendrites of very complex shape. The inclusion of the Gibbs-Thomson effect is necessary for the accurate resolution of the dynamics of unsteady solidification. Finally, various Dirichlet or Neumann or mixed type boundary conditions can be applied on the outside boundaries of the multi-phase domain. We call the heat conduction problem in multi-material system without phase transitions, described by equations (1)–(3), the parabolic interface problem. A numerical algorithm for the parabolic interface problem based on the embedded boundary method is proposed in this paper. Phase boundaries described by the jump condition (5) decouple the problem in a set of parabolic problems with geometrically complex and moving external boundaries (phase boundaries) with the Dirichlet boundary condition Tboundary = TM (see Figure 2).
Fig. 1
Schematic of a multi-material domain with phase transition and its decoupling
into the parabolic interface problem (top right) containing the material interface and the moving boundary (the inner phase boundary with the Dirichlet boundary condition Tphase
boundary
= TM ) and the elliptic problem (bottom right) with the moving outside
boundary (phase boundary with the same Dirichlet boundary condition)
3
Elliptic and Parabolic Interface Problem
In this section, we extend the embedded boundary method for irregular domains [10, 19, 23] to solve the elliptic and parabolic interface problem with interior boundaries.
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
503
The elliptic partial differential equation is: ∇ β∇T = f,
(7)
where β is a continuous function except across the interior boundary and f is also assumed to be a given function which is continuous inside each part of the domain. To close the problem, boundary condition are needed for the exterior and interior boundary. Either Dirichlet or Neumann boundary can be given on the exterior boundary. For the interior boundary, we have the following jump conditions: [T ] = J1 (x),
(8)
∂T = J2 (x), (9) β ∂n where J1 and J2 are given functions of space variables. We generalize the parabolic interface problem introduced in the previous section in the following way. The partial differential equation is ∂T = ∇ β∇T + f, ∂t
(10)
where β and f are continuous functions except across interfaces. Either Dirichlet or Neumann boundary can be given on the exterior boundary. For the interior boundary, we generalize the boundary condition (3) to include the following: [T ] = J1 (x, t),
(11)
∂T = J2 (x, t), (12) β ∂n where J1 and J2 are given functions of space and time. These two jump conditions across the interior interface are similar with those for the elliptic equation (8), (9). The embedded boundary method [10, 19, 23] is a finite volume method for an irregular domain embedded on a Cartesian grid. When the EBM is used to solve the elliptic boundary value problem (without interfaces), the unknowns are defined at the cell centers, even for partial cells which have intersections with the boundary. We retain this main feature of the EBM method in our algorithm for the elliptic/parabolic interface problem: unknowns are defined at centers of regular grid cells. For a full cell (cell containing no external boundaries or interfaces) or a boundary cell (cell cut by the exterior boundary), one unknown is needed at the cell center and a standard finite volume method can be used to setup one algebraic equation as in [10, 19, 23]. For a partial cell (cell cut by the interface), more unknowns are needed to make the discretization of the elliptic/parabolic equation consistent with the two interface jump conditions. Figure 2 depicts the placement of unknowns in a partial cell with interface. The whole cell contains two partial cells (for material components A and B) which are separated by the interface. For each partial cell, there is one unknown defined at the cell center (T1 for the A side, T2 for the B side). In order to satisfy the two jump conditions (8), (9), two more unknowns Ta and Tb are defined at the center of the cell interior boundary (portion of the interface contained within the cell, called cell interface) for the two sides A and B. Therefore, for the partial cell with the interface, four unknowns are defined in total. Four algebraic equations
504
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
can be constructed for this kind of cells using the two jump conditions on the cell interior boundary and the PDE for the two partial cells.
Fig. 2
Placement of unknowns in a cell containing the interface
3.1
Equations for the jump conditions In this section, we describe algorithms for the discretization of the two jump conditions across the interface for the cell (i, j). A schematic of the corresponding stencil and states used in the interpolation method are shown in Figure 3. Two unknowns Ta and Tb are defined at the center of the cell interface for sides A and B, respectively.
Fig. 3
Stencil for the interface unknowns for the jump condition
Let us specify the direction of the normal to the interface as pointing from A to B. For the first jump condition (8), the discretization is simply Tb − Ta = J1 .
(13)
To discretize the second jump condition (9), we need to calculate the normal derivatives of the unknowns. One approach is to use the method given in [10, 19] which fits the quadratic polynomial into a stencil normal to the interface with states in the stencil nodes obtained via interpolation from interior points. Our method of fitting a quadratic polynomial with two variables is explained in detail in Section 4.2 in the later part of the paper. The main idea is to use 1 unknown in cell (i, j) (either cell center unknown or cell interface unknown depending on usage), 2 unknowns in the first layer neighbors of cell (i, j) and 3 unknowns in the second layer neighbors of cell (i, j) which in total give 6 unknowns for the 6 coefficients of the quadratic polynomial. For example, to construct the quadratic polynomial to calculate the flux for side A at the cell interface center where Ta is defined, the 6 unknown candidates are Ta , 2 cell center
No.2
505
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
unknowns Ti+1,j , Ti+1,j−1 in the first layer and 3 cell center unknowns Ti+2,j , Ti+2,j−1 , Ti+2,j−2 in the second layer. Similarly we can fit a quadratic polynomial with two variables for side B. Taking the normal derivative of the fitted polynomials, we obtain the normal derivatives ∂T at the cell interface center for the side A, ∂T ∂n |a , and the side B, ∂n |b . Using the second jump condition (9), we have ∂T ∂T (14) β · − β · = J2 . b ∂n b a ∂n a 3.2 Equations for the cell center unknowns The setup of algebraic equations corresponding to the discretization of the parabolic equation depends on the choice of the time stepping method and the discretization of the elliptic operator. The temporal discretization will be discussed in the next section. Here we describe the discretization of the elliptic operator ∇ β∇T = f.
(15)
For the two cell center unknowns defined at cell (i, j), we use the EBM technique to setup two equations separately. See Figure 4 for the stencil to setup the unknown for side A using the partial cell cdef .
Fig. 4
Stencil for the cell center unknowns
Integrating the equation (15) around the partial cell cdef and using the divergence theorem, we obtain the following expressions: ∇ · β∇T dxdy = β∇T · ndS = f dxdy, ∂(cdef )
cdef
or
β∇T · ndS +
cd
β∇T · ndS +
de
cdef
β∇T · ndS +
ef
β∇T · ndS =
fc
f dxdy, cdef
which is lcd · Fluxcd + lde · Fluxde + lef · Fluxef + lf c · Fluxf c =
f dxdy
(16)
cdef
where lmn is the length between m and n. Therefore we only need to calculate the flux across the cell edges or cell interface. For Fluxcd , a 2nd order derivative is calculated by using a T −Ti,j T −Ti+1,j linear interpolation of i,j−1 and i+1,j−1 to the center of cd (see [10]). For Fluxde , x x
506 we simply use Ti,j+1 −Ti,j x
ACTA MATHEMATICA SCIENTIA Ti+1,j −Ti,j to y Ti+1,j+1 −Ti+1,j x
Vol.30 Ser.B
calculate the derivative. For Fluxef , a linear interpolation of
and to the center of ef is used. Fluxf c , is calculated using ∂T ∂n |a used in (14). Note that we need to multiply the derivatives calculated with β at the corresponding edge or cell interface center to get the flux. In the same way, we calculate fluxes for the other partial cells. More details can be found in [10]. It is also worth noting that the two unknowns defined at the cell interface can be represented by the cell centers unknowns by solving the two discretized jump conditions (13), (14). This can result in a system of equations with only unknowns defined at cell centers. 3.3 Implicit time stepping methods The embedded boundary method is used with the implicit Runge-Kutta method (TGA, [25]) for solving the time dependent parabolic initial boundary value problem in ([19, 23]). Here we apply this time stepping technique to solve the parabolic interface problem. As in [19], in order to solve the following time-dependent initial boundary value problem: dU = LhI (t)U (t) + F (t), dt
(17)
U (0) = U 0 ,
(18)
where LhI (t) is the time-dependent inhomogeneous operator, we use the two step implicit RungeKutta method (TGA, [25]): U n+1 = (I − μ1 Lh1 (tnew ))−1 (I − μ2 LhI (tint ))−1 ·[(I +
μ3 LhI (told ))U n
+ (I +
μ4 LhH )f (tavg )t],
(19) (20)
where LhH (t) is LhI (t) with homogeneous boundary conditions and told = nt, tnew = (n + 1)t = told + μ1 + μ2 + μ3 , tint = tnew − μ1 = told + μ2 + μ3 , tavg = (told + tnew )/2 = told + μ1 + μ2 + μ4 , √ a − a2 − 4a + 2 t, μ1 = √ 2 a + a2 − 4a + 2 μ2 = t, 2 μ3 = (1 − a)t, 1 μ4 = ( − a)t. 2
√ To get a L0 -stable method with real arithmetic, it is required that 1/2 < a < 2 − 2. In this case, μ1 > 0, μ2 > 0, μ3 > 0 > μ4 . Therefore, told < tavg < tint < tnew . When a = 1/2, we obtain the Crank-Nicolson method with μ1 = μ4 = 0 and tint = tnew . It should be noted that the homogeneous boundary conditions of the operator LhH corresponding to the interior boundary jump condition (8), (9) are to set the jumps to be zero [T ] = 0, ∂T = 0. β ∂n
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
507
For the boundary conditions corresponding to the exterior boundary, we set T = 0 for Dirichlet boundary or set ∂T ∂n = 0 for Neumann boundary. Equation (19) are solved in two step: (I − μ2 LhI (tint ))U tint = (I + μ3 LhI (told ))U n + (I + μ4 LhH )f (tavg )t
(21)
(I − μ1 Lh1 (tnew ))U n+1 = U tint
(22)
Each equation is discretized using the previously developed embedded boundary method. For example, equation (21) is the following: U tint − μ2 ∇ β(tint )∇U tint = U n + μ3 ∇ β(tint )∇U n + t[f (tavg ) + μ4 ∇ β∇f (tavg )]. 3.4
Moving boundary problem
In this Section, we generalize the EBM algorithm for problems with moving interfaces. Our motivation is to develop an accurate numerical methods for the sharp interface classical Stefan problem. For the moving boundary problem, the EBM method for elliptic problems with moving exterior boundaries ([19, 23]) extrapolates the boundary condition located at the position of the intermediate step tint to the position of the new boundary at time tnew using Taylor series. It is impractical to apply similar method for interior boundaries with the jump conditions (8), (9). Instead we solve the first step of the implicit Runge-Kutta (21) using the intermediate boundary at time tint and solve the second step (22) using the new boundary at time tnew . The uncovered states for the intermediate and new interface are obtained via extrapolation. Our extrapolation method is different from those in ([19, 23]): the extrapolation is performed through the center of the nearest cell interface by fitting a quadratic polynomial with two variables which is the same as the fitting used in previous section to calculate the cell interface flux to discretize the parabolic equation. Figure 5 shows the extrapolation method. Both cell center states and cell interface states may need to be extrapolated when the interface is moved. The extrapolation method is the same for both cases. Note that we do not do the extrapolation in the normal to the interface direction and it is not necessary. In Figure 5, to extrapolate state to the new cell interface at a for side A, we first find the nearest old cell interface d, then find the candidates for fitting in the direction along line ad.
Fig. 5
Extrapolation method for moving boundary
508
4 4.1
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
Multiple Components and Interpolation Method Algorithms and coding for multiple components
For numerous practical applications, the assumption of at most two material components in a single cell is not correct even for sufficiently refined meshes. Figure 6 shows a computational domain with three material components that cross at a single point. Unless the intersection point is precisely on the mesh line, there is one cell containing three different components. In order to deal with this set of problems, we would like to extend our original method but preserve the central idea of the method that unknowns are defined at cell centers (even for irregular cells) and on cell interface centers. We first need to clarify the possible cases of cells in problems with multiple components with assumption that there is at most one crossing by the exterior or interior boundary on each cell edge. We also assume that no crossing occur at the nodes of the Cartesian grid. We already have external cells which is located outside the computational domain, internal cells without exterior or interior boundary, boundary cells containing exterior boundary, partial cells containing interior boundary with at most 2 components. We need to deal with cells with 3 components as in Figure 7 and cells 4 components as in Figure 8. These two cases will still exist even under mesh refinement. There are also cases which can be dealt with using mesh refinement as in Figure 9 and Figure 10. However, we also deal with them to make the code to be more robust.
Fig. 6
Computational domain with
3
3
1
2
Fig. 7
Cell with 3 components
three components A, B, C 3
4
2
1
1
2
1
2
Fig. 8
Cell with 4 components
Fig. 9
Cell with 2 components
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
509
The algorithm to setup the unknown matrix for solving elliptic equation with multiple components is the following: For each Cartesian cell, do the following: 1. for each cell interface, setup the equations related to the unknowns on the cell interface and the two components. 2. for each cell edge, setup the equations. 3. for each cell interior, setup the interior equation corresponding to that component. To setup the unknown matrix, we need to enumerate the cell interface. Each cell interface is identified by a identification number (ID) consisting of three numbers (a, b, c): a and b are the two ordered components for the two side of the cell interface; c is an integer starting from 0 which is needed when two cell interface has the same two components for cases as in Figure 9. The information stored at the cell center for a particular component can be identified by the component number. Each cell contains many information including number of cell unknowns, cell components, unknowns index, cell areas for each components, pointer to states, number of cell interface unknowns, cell interface length, cell interface normal, cell interface center. It is impractical to use a single data structure to store simple cells and multiple component cells (which need to store much more information). To save computational memory, two different data structure are used: one for cells holding at most 2 components and one for cells holding more complex cells as in Figures 7–10. More memory can be saved to use structures for cells with 1 component and cells with 2 components. It is notable that cell containing multiple components are in the order of O(1) if the exterior and interior boundary is not very complex. 4.2
Interpolation method
To deal with multiple components in a single cell and improve the robustness of the code, the interpolation method need to be modified. As in Figure 11, our task is to get a bilinear polynomial using three unknowns or biquadratic polynomial using six unknowns located at different locations. Taking the derivative of these polynomial in the normal direction of the cell interface, we get the derivative which is needed to setup the matrix for the elliptic equations. If quadratic interpolation is not possible, linear interpolation is used which is always possible under our assumption that there is at most one crossing on each cell edge by the exterior boundary or interior interface.
Fig. 10
Cell with 3 components
Fig. 11
Interpolation method
510
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
When selecting unknowns to do interpolation for one of the component, at most one of the unknowns of a cell (located at cell center or cell interface) with the same component is selected. Otherwise, their positions might be too close, leading to the singularity of the matrix for the interpolation coefficients. It is also required that at most one cell interface unknown (say the current cell (i, j)) is selected for an interpolation stencil. We classify the cells into different group of neighbors. Zeroth layer neighbor: for the cell (i, j), the zeroth layer neighbor is simply the current cell (i, j). First layer neighbor: for the cell (i, j), its first layer neighbors are cells with the cell index (m, n) satisfying max(|m − i|, |n − j|) = 1. Thus the first layer cells include (i − 1, j − 1), (i, j − 1), (j + 1, j − 1), (i − 1, j), (i, j), (i + 1, j), (i − 1, j + 1), (i, j + 1), (i + 1, j + 1). Second layer neighbor: the second layer neighbors for cell (i, j) are similarly defined as cells with cell index (m, n) satisfying max(|m − i|, |n − j|) = 2. Consecutive cells: the consecutive cells in this paper denote the consecutive cells in a given layer. For example, consider the first layer of cell (i, j): cell (i + 1, j) is consecutive with cell (i + 1, j + 1) and (i + 1, j − 1) in the first layer; it is not consecutive with any other cells in the first layer. Similarly, consider the second layer of cell (i, j): cell (i + 2, j + 2) is only consecutive with cell (i + 1, j + 2) and (i + 2, j + 1); it is not consecutive with any other cells in the second layer. Now we can select unknowns to do interpolation. Considering cell (i, j) with its neighbors, the cells in Figure 11 consist of cell (i, j), its first and second layer neighbors. To do bilinear interpolation, we select from cell (i, j) and two consecutive cells from its first layer neighbors, using only one unknown from each cell resulting to 3 unknowns in total for the three coefficients of linear interpolation. To do biquadratic interpolation, we select another three consecutive cells from its second layer of neighbors besides the three unknowns for the first order interpolation resulting to 6 unknowns in total for the six coefficients of biquadratic interpolation. Only cell center unknowns are used in the interpolation except the current cell (i, j) for which either the cell center or the cell interface unknown is selected depending on what is needed. When selecting candidate unknowns for interpolation, it is preferred to first check along the cell interface normal direction. In Figure 11, a candidate stencil for linear interpolation of solution with component B includes cell (i, j), (i + 1, j), (i + 1, j + 1) for linear interpolation. Add cell (i + 2, j + 1), (i + 2, j + 2), (i + 1, j + 2) for quadratic interpolation. We have the following three observations for boundary or partial cells. Proposition 1 For a cell containing component a at one of its four corners, it is always possible to find three unknowns with the same components to do bilinear interpolation. Proof Denote the current cell index as (i, j). Since it contains component a, there exists one cell interface with one side of its component as a. Thus we can pick the cell unknown or cell interface unknown of cell (i, j). Suppose that the down-right corner of the cell (i, j) has
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
511
component a. Then it must be the case that cell (i, j − 1) contains component a at its up-right corner. Similarly, cell (i + 1, j − 1), (i + 1, j) both contain the same components. Therefore, two more unknowns can be selected from either {(i, j − 1), (i + 1, j − 1)}, or {(i + 1, j − 1), (i + 1, j)}. Proposition 2 The matrix for the 3 selected unknowns for linear interpolation are always solvable (non-singular). Proof For cell (i, j), suppose that an interpolation stencil is (i, j) , (i + 1, j), (i + 1, j + 1). The interpolation algorithm will pick cell center unknown from the last two cell and either cell center or cell interface unknown from the first cell. It is easy to know that any point in cell (i, j) can not be on the same line passing through the center of the last two cell. Thus the matrix for this particular stencil is non-singular. Proposition 3 The matrix for the 6 selected unknowns for quadratic interpolation are always solvable (non-singular). Proof This statement can be verified by the same reasoning as above.
5
Numerical Results
In this section, we validate the proposed numerical algorithms by solving elliptic and parabolic interface problems with interior boundaries and comparing numerical solutions with exact analytical ones. We also discuss the accuracy and convergence of numerical schemes. 5.1 Example 1 elliptic interface problem We solve the elliptic interface problem inside [−1, 1] × [−1, 1]. The interior interface in Figure 12 is described by (x − 0.1)2 + (y − 0.1)2 = 0.5 + 0.2 × cos(θ), where θ is the angle between (x − 0.1, y − 0.1) and x-coordinate.
Fig. 12
Interface for example 1 elliptic interface problem
Denote the part inside the interface as A and the part out of the interface as B. We use the following coefficient ⎧ ⎨ 1 (x, y) ∈ A, β(x, y) = ⎩ 1 (x, y) ∈ B, 10
512
ACTA MATHEMATICA SCIENTIA
and exact solution T (x, y) =
⎧ ⎨ (x2 + y 2 ) 32
(x, y) ∈ A,
⎩e
(x, y) ∈ B.
x2 +y 2
Vol.30 Ser.B
The approximate solution is given in Figure 13. Table 1 shows the max errors under mesh refinement. The method is 2nd order accurate. Table 1
Max errors for example 1 elliptic interface problem
Fig. 13
5.2
Mesh Size
Max Error
Ratio
40 × 40
1.249277e-02
80 × 80
3.290830e-03
3.80
160 × 160
8.442124e-04
3.90
320 × 320
2.137630e-04
3.95
640 × 640
5.378009e-05
3.97
Approximate solution for example 1 (160 × 160)
Example 2 Elliptic interface problem with multiple components The computational domain is [−1, 1] × [−1, 1]. It is partitioned in to three different part A = {(x, y)|x < 0.3, y < 0.3}, B = {(x, y)|x ≥ 0.3, y < 0.3}, C = {(x, y)|y ≥ 0.3}. Figure 14 shows the computational domain with interior interfaces.
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
Fig. 14
Interface for Example 2 elliptic interface problem with multiple components Table 2 Max errors with mesh refinement for Example 2 elliptic interface problem with multiple components
Fig. 15
Mesh Size
Max Error
Ratio
20 × 20
1.560965e-02
40 × 40
3.942480e-03
3.96
80 × 80
9.901261e-04
3.98
160 × 160
2.480690e-04
3.99
320 × 320
6.208303e-05
4.00
Approximate solution for Example 2 elliptic interface problem with multiple components (160 × 160)
We use the following coefficient
⎧ ⎪ 1 ⎪ ⎪ ⎨ 1 β(x, y) = ⎪ 10 ⎪ ⎪ ⎩3
(x, y) ∈ A, (x, y) ∈ B, (x, y) ∈ C,
513
514
ACTA MATHEMATICA SCIENTIA
and exact solution
⎧ 2 2 32 ⎪ ⎪ ⎨ (x + y ) T (x, y) =
⎪ ⎪ ⎩
e
(x, y) ∈ A,
x2 +y 2
2(x2 + y 2 )
Vol.30 Ser.B
(x, y) ∈ B, 3 2
(x, y) ∈ C.
Table 2 shows the max errors and Figure 15 shows the approximate solution. 5.3 Example 3 elliptic interface problem in 3D The computational domain is [−1, 1] × [−1, 1] × [−1, 1]. The interface is a circle with center at (0, 0, 0) with radius r = 0.5. Let A be the region inside the interface and B the region outside the interface. We use the following as the coefficient: ⎧ ⎨ 1 (x, y) ∈ A, β(x, y) = ⎩ 10 (x, y) ∈ B, and exact solution: T (x, y) =
⎧ 3 ⎪ ⎨ (x2 + y 2 ) 2 2
2
(x, y) ∈ A,
3 2
(x + y ) 1 1 ⎪ ⎩ + (1 − ) × 10 10 8
(x, y) ∈ B.
Currently, the function interpolation method used for the 3D code is similar with those in [23]. Table 3 shows the max errors with mesh refinement. Table 3
Max error for Example 3 elliptic interface problem for 3D Mesh Size
Max Error
Ratio
10 × 10
1.077035e-01
20 × 20
3.451512e-02
3.12
40 × 40
9.139040e-03
3.78
80 × 80
2.362821e-03
3.87
160 × 160
6.045999e-04
3.91
5.4
Example 4 parabolic interface problem with fixed interface For this example, we solve the parabolic interface problem on a fixed interface. The computational domain is [−1, 1] × [−1, 1]. The interface is a circle with center at [−0.1, −0.1] with radius r = 0.5. Let A be the region inside the interface and B the region outside the interface. We use the following coefficient: ⎧ ⎨ 1 (x, y) ∈ A, β(x, y) = ⎩ 10 (x, y) ∈ B, and the exact solution: ⎧ x2 +y 2 ⎪ 4 exp(− 5(t+1) ) ⎪ ⎪ ⎪ ⎨ 5π(t + 1) T (x, y) = x2 +y 2 ⎪ ⎪ exp(− 2(t+1) ) ⎪ ⎪ ⎩ 2π(t + 1)
(x, y) ∈ A, (x, y) ∈ B,
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
515
where the formula for T (x, y), x ∈ A is picked from [23]. The heat equation is solved until t = 0.5. Table 4 shows the max errors at time t = 0.5. Figure 16 shows the approximate solution at time t = 0.5 with mesh size 160 × 160. The test confirms the second order accuracy of the method. Table 4
Fig. 16
5.5
Maximum errors at time t = 0.5 with mesh refinement for Example 4 parabolic interface problem with fixed interface, t = 0.9 h Mesh Size
Max Error
Ratio
20 × 20
3.709544e-04
40 × 40
8.830904e-05
4.20
80 × 80
2.257112e-05
3.91
160 × 160
5.705129e-06
3.96
Approximate solution for Example 4 parabolic interface problem with fixed interface
Example 5 parabolic interface problem with prescribed boundary motion In the second numerical test, we solve the parabolic interface problem with interface that moves with a given velocity. The computational domain is [−1, 1] × [−1, 1]. The initial interface is a circle centered at [−0.25, −0.25] with radius r = 0.4. The motion of the circle is defined simply by dx = 0.5, dt dy = 0.5. dt Other conditions are the same as Example 4. Thus the motion is simply a rigid motion. Table 5 shows the max errors under the mesh refinement. Figure 17 shows the approximate solution at time t = 0.2 with mesh size 320 × 320. The time step in the simulation was restricted by t h = 0.1. The increase of the time step t to h = 0.9 reduces the order of accuracy to approximately 1.2. To retain the second order accuracy, the time step must be lower than the stability limit t t max v ,β < 1. h h
516
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
The last expressions assumes en explicit method for the parabolic equation. Since the velocity of the moving interface is v = 0.5, t h < 0.9 satisfies the stability condition v
t < 1. h
However the increase of the time step to t h = 0.9 results in the loss of accuracy caused by the extrapolation of states to the location of cells swept during the interface motion. Table 5 Maximum errors at time t = 0.2 with mesh refinement for Example 5 parabolic interface problem with prescribed boundary motion, t = 0.1 h
Fig. 17
Mesh Size
Max Error
Ratio
40 × 40
4.103197e-05
80 × 80
1.008246e-05
4.07
160 × 160
2.895398e-06
3.49
320 × 320
6.377383e-07
4.54
Approximate solution at time t = 0.2 for Example 5 parabolic interface problem with prescribed boundary motion.
t h
= 0.1
5.6
Example 6 Stefan problem In this numerical example, we validate the Stefan problem algorithm using equations that possess an exact analytical solution. Consider a square domain [−1, 1] × [−1, 1] containing two phases (solid and liquid) of the same substance. The phases are separated by a circular interface that moves as follows 1 R0 (t) = 3 . − 4 (t + 4)4 + 23 4 The liquid is inside the circle. The temperature of the liquid satisfies the following analytical formula T (r) = (R0 (t)3 − r3 )(t + 1)3 ,
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
517
The temperature is equal to 0 on the interface. The analytical formula for the solid phase temperature is simply 0. The PDE is Tt = T + f, where f can be obtained by substituting the exact solution. The starting time is t0 = 0.4 and the final time T = 0.45. Figure 18 shows the maximum interface position errors and Figure 19 shows the maximum solution errors. Both plots demonstrate the 2nd order accuracy of the algorithm. −4 log(max intfc error) ~ log(dx) y=2.0153*log(dx) + 0.3868 −5
log(max intfc error)
−6
−7
−8
−9
−10 −5.5
−5
−4.5
−4
−3.5
−3
−2.5
−2
log(dx)
Fig. 18
Maximum interface position error for Example 6 Stefan problem at time T = 0.45 −3 log(max solution error) ~ log(dx) y=1.9153*log(dx) + 0.5717 −4
log(max solution error)
−5
−6
−7
−8
−9
−10 −5.5
−5
−4.5
−4
−3.5
−3
−2.5
−2
log(dx)
Fig. 19
6
Maximum solution error for Example 6 Stefan problem at time T = 0.45
Application
In this section, we illustrate the application of our algorithm to the study of accident scenarios in nuclear fuel rods. This study is being performed within multi-institutional project
518
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
on multiscale simulations of Generation-IV Sodium Fast Reactors [2] focusing on reactor fuel and core transient response under beyond-design and accident conditions. A nuclear fuel rod contains 8 mm diameter pellets of metallic or oxide uranium fuel in a stainless still cladding of 0.5 mm thickness. The total length of the fuel rod is about 2.5 m. A narrow (0.1 mm), irregular gap between the fuel and the cladding is used to transport fission gasses from the fuel to the gas plenum on the top of the reactor core. Fuel rods are assembled in hexagonal structures and places into liquid sodium coolant that circulates in the reactor core and transports the thermal energy from fuel rods to the heat exchangers. The cross-section of the fuel rod is shown schematically in Figure 20.
Fig. 20
Schematic of the nuclear fuel rod cross-section: (a) is the fuel pellet,
(b) is the gas gap, (c) is the stainless steel cladding, and (d) is the liquid sodium.
The temperature distribution across the fuel rod at normal operating conditions is shown in Figure 22. The power production of about 60 kW/m in the fuel rod is balanced by the heat removal by the sodium coolant so that the temperature is in the steady state. The picture shows rapid changes of the temperature gradients across the fuel rod due to sharp changes of the heat conductivity in the fuel, gap, and cladding. After years of operation, the surface of the fuel erodes and comes to a point-wise contact with the clad as shown in Figure 21. The interior of the gas gap is not resolved on the mesh level and an empirical formula for the temperature jump across the gap [6] is applied for the heat transport calculations. For the related problem of the transport of fission gases in the gas gap, the flow in porous media equations are solved.
Fig. 21
Shape of the gas gap between the fuel and cladding
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
3500
Melting Point Overheating State Normal State
3000
Temperature (K)
519
2500 2000 1500 1000 500 0
Fig. 22
0.1
0.2
0.3
0.4 0.5 Radius (cm)
0.6
0.7
0.8
0.9
Temperature across the fuel rod at normal operating conditions and transient overheating
During either transient overheating or loss of coolant accidents, the heat transfer balance is changed and the temperature in the fuel rod increases. It can cause melting of the fuel rod and even stainless steel cladding. However because of the change of cladding properties with the temperature increase, the cladding usually fails before melting, causing the ejection of molten fuel and fission gases into the coolant reservoir. Figures 23 and 22 depict the increase of temperature and melting of the fuel in the transient overheating accident.
Fig. 23
Temperature distribution across the fuel rod during transient overheating and the surface of the molten fuel
7
Conclusions
In this paper, we generalized the embedded boundary algorithm for solving elliptic and parabolic problems in geometrically complex domains [10] to include irregular interior boundaries or interfaces. The algorithm was validated by solving problems that possess exact analytical solutions and the second order accuracy in both space and time was confirmed. Based on this method, we developed an algorithm for solving the classical Stefan problem that describes first order phase transitions such as melting and solidification. The Stefan problem algorithm was also validated using an exact analytical solution. The combination of elliptic
520
ACTA MATHEMATICA SCIENTIA
Vol.30 Ser.B
and parabolic algorithms makes it possible to simulate the heat transfer and phase transitions in multi-material systems containing geometrically complex interfaces between different materials and phase boundaries or surfaces between different phases of the same material. The explicit tracking of material interfaces removes all restrictions on the magnitude of discontinuities of physical properties and temperature gradients across interfaces and phase boundaries and is especially beneficial for the study of materials at extreme conditions. We illustrated the use of the method in the study of accident scenarios in nuclear fuel rods in Generation-IV Sodium Fast Reactors. The embedded boundary method was used to simulate the temperature increase in the fuel assembly and melting of uranium fuel during the transient overpower and loss of coolant accidents. The algorithm was implemented in FronTier, a multiphysics code for the simulation of multiphase flows using the method of front tracking. In the future, the embedded boundary method will also be implemented in the incompressible hydrodynamics and magnetohydrodynamics components of the FronTier code. Acknowledgments The authors would like to thank James Glimm and Xiaolin Li for stimulating discussions. This manuscript has been authored in part by Brookhaven Science Associates, LLC, under Contract No. DE-AC02-98CH10886 with the U.S. Department of Energy. The United States Government retains, and the publisher, by accepting the article for publication, acknowledges, a world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for the United States Government purpose. This research utilized resources at the New York Center for Computational Sciences at Stony Brook University/Brookhaven National Laboratory. References [1] Alexiades V, Solomon A D. Mathematical Modeling of Melting and Freezing Processes. McGraw-Hill, 1993 [2] Bolotnov I A, Behafarid F, Shaver D, Antal S P, Jansen K E, Samulyak R, Podowski M Z. Interaction of computational tools for multiscale multiphysics simulation of generation-iv reactors//Proceedings of International Congress on Advances in Nuclear Power Plants, 2010. To be published [3] Caginalp G, Fife P. Higher-order phase field models and detailed anisotropy. Phys Rev B, 1986, 34: 4940–4943 [4] Caginalp G. Stefan and Hele-Shaw type models as asymptotic limits of the phase field equations. Phys Rev A, 1989, 39: 5887–5896 [5] Du Jian, Fix Brian, Glimm James, Jia Xicheng, Li Xiaolin, Li Yunhua, Wu Lingling. A simple package for front tracking. J Comput Phys, 2006, 213: 613–628 [6] Duderstadt J, Hamilton L. Nuclear Reactor Analysis. 2nd ed. John Wiley & Sons, Inc, 1976 [7] Fix G.//Fasano A, Promiceiro M, ed. Free Boundary Problems. London: Pitman, 1983: 580 [8] George E, Glimm J, Li X L, Li Y H, Liu X F. The influence of scale-breaking phenomena on turbulent mixing rates. Phys Rev E, 2006, 73: 1–5 [9] Glimm J, Grove J W, Li X -L, Shyue K -M, Zhang Q, Zeng Y. Three dimensional front tracking. SIAM J Sci Comp, 1998, 19: 703–727 [10] Johansen H, Colella P. A cartesian grid embedding boundary method for poisson’s equation on irregular domains. J Comput Phys, 1998, 147: 60–85 [11] Sullivan J M, Lynch D R, O’Neill K. Finite-element simulation of planare instabilities during solidification of an undercooled melt. J Comput Phys, 1987, 69: 81–111
No.2
S.Q. Wang et al: AN EMBEDDED BOUNDARY METHOD
521
[12] Juric D, Tryggvason G. A front tracking method for dendritic solidi-fication. J Comput Phys, 1996, 123: 127–148 [13] Langer.//Grinstein G, Mazenko G, ed. Directions in Condensed Matter Physics. Singapore: Worls Scientific, 1986: 164 [14] LeVeque R J, Li Z L. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM J Numer Anal, 1994, 31: 1019–1044 [15] Liu X -D, Fedkiw R, Kang M. A boundary condition capturing method for poisson’s equation on irregular domains. J Comput Phys, 2000, 160: 151–178 [16] Liu X F, Li Y H, Glimm J, Li X L. A front tracking algorithm for limited mass diffusion. J Comput Phys, 2007, 222: 644–653 [17] Lu T, Xu Z L, Samulyak R, Glimm J, Ji X M. Dynamic phase boundaries for compressible fluids. SIAM J Sci Comp, 2008, 30: 895–815 [18] Mayo A. The fast solution of poisson and the biharmonic equations on irregular regions. SIAM J Numer Anal, 1984, 21: 285–299 [19] McCorquodale P, Colella P, Johansen H. A cartesian grid embedded boundary method for the heat equation on irregular domains. J Comput Phys, 2001, 173: 620–635 [20] Mckenney A, Mayo A, Greengard L. A fast poisson solver for complex geometries. J Comput Phys, 1995, 118: 348–355 [21] Perskin Charles S. The immersed boundary method. Acta Numerica, 2002, 11: 479–517 [22] Samulyak R, Du J, Glimm J, Xu Z. A numerical algorithm for MHD of free surface flows at low magnetic reynolds numbers. J Comput Phys, 2007, 226: 1532–1546 [23] Schwartz P, Barad M, Colella P, Ligocki T. A cartesian grid embedded boundary method for the heat equation and Poisson’s equation in three dimensions. J Comput Phys, 2006, 211: 531–550 [24] Strain J, Sethian J. Crystal Growth and dendritic solidification. J Comput Phys, 1992, 98: 231–253 [25] Twizell E H, Gumel A B, Arigu M A. Second-order, l0-stable methods for the heat equation with timedependent boundary conditions. Adv Comput Math, 1996, 6: 333 [26] Wheeler A A.//Hurle D T J, ed. Handbook of Crystal Growth, Volume 1B. Amsterdam: North-Holland, 1993: 783