International Journal of Partial Differential EquationsVolume 2013 (2013), Article ID 476873, 13 pageshttp://dx.doi.org/10.1155/2013/476873

Research Article

A Numerical Method for Solving 3D Elasticity Equations with Sharp-Edged Interfaces

Liqun Wang,¹ Songming Hou,² and Liwei Shi²

¹Department of Mathematics, College of Science, China University of Petroleum, Beijing 102249, China²Department of Mathematics and Statistics, Louisiana Tech University, Ruston, LA 71272, USA

Received 11 April 2013; Revised 16 June 2013; Accepted 25 June 2013

Academic Editor: Yuri N. Skiba

Copyright © 2013 Liqun Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Interface problems occur frequently when two or more materials meet. Solving elasticity equations with sharp-edged interfaces in three dimensions is a very complicated and challenging problem for most existing methods. There are several difficulties: the coupled elliptic system, the matrix coefficients, the sharp-edged interface, and three dimensions. An accurate and efficient method is desired. In this paper, an efficient nontraditional finite element method with nonbody-fitting grids is proposed to solve elasticity equations with sharp-edged interfaces in three dimensions. The main idea is to choose the test function basis to be the standard finite element basis independent of the interface and to choose the solution basis to be piecewise linear satisfying the jump conditions across the interface. The resulting linear system of equations is shown to be positive definite under certain assumptions. Numerical experiments show that this method is second order accurate in the norm for piecewise smooth solutions. More than 1.5th order accuracy is observed for solution with singularity (second derivative blows up).

1. Introduction

The macroscopic mechanical response of most elastic materials, such as wood, metal, rubber, rock, and bone, can be described by the elasticity theory. However, the numerical solution for an elasticity problem with interface is highly nontrivial, because the coefficients in the partial differential equation for elasticity system are discontinuous, and the solution to the PDE system also needs to satisfy the jump conditions across each material interface. Examples of elasticity interface problem include the minimum-compliance-design problems [1], the microstructural evolution [2], and the atomic interactions [3]. Designing an accurate and efficient method for these problems is a difficult job, especially when the interface is not smooth.

In this paper, we consider a three-dimensional elasticity problem with sharp-edged interfaces and matrix coefficients. The model of interest is as follows.

Consider an open bounded domain . Let be an interface of codimension , which divides into disjoint open subdomains, and , hence . Assume that the boundary and the boundary of each subdomain are Lipschitz continuous as submanifolds. Since are Lipschitz continuous, so is . A unit normal vector of can be defined as a.e. on ; see Section 1.5 in [4].

We seek solutions of the variable coefficient elliptic equation away from the interface given by in which denotes the spatial variables and is the gradient operator. The coefficient is assumed to be a matrix that is uniformly elliptic on each disjoint subdomain, and , and its components are continuously differentiable on each disjoint subdomain, but they may be discontinuous across the interface . The right-hand side , is assumed to lie in .

Given functions and along the interface , we prescribe the jump conditions:

The “” superscripts refer to limits taken from within the subdomains .

Finally, we prescribe boundary conditions for a given function on the boundary .

The setup of the problem is illustrated in Figure 1.

Figure 1: Setup of the problem.

We consider the three-dimensional elasticity interface problem () in this paper. Although the real problem would have a three-by-three coupled PDE system and six jump conditions, for simplicity of coding, we use two-by-two PDE system to demonstrate the effectiveness of our method. Note that the mathematical challenge is in the 3D geometry when generalized from the corresponding 2D problem. Changing the system from two-by-two to three-by-three does not bring in new mathematical challenge.

Interface problems have been studied by many researchers and several techniques have been developed. Finite element methods with body-fitted grid use a mesh triangulation to capture the boundary [5–8]. However, in many situations, such as when the boundary is moving, the mesh generation may be both computationally expensive and challenging. A uniform Cartesian grid is preferred.

In [9, 10], to simulate the flow patten of blood in the heart, Peskin proposed the “immersed boundary” method, which used numerical approximation of -function for singular sources on the interface. In [11], in order to compute two-phase flow, a level-set method was combined with the “immersed boundary” method. These methods are simple to be used but difficult to achieve high-order accuracy.

In addition, a large class of finite difference methods have been proposed. The main idea is to use difference scheme and stencils carefully near the interface to incorporate jump conditions and achieve high-order local truncation error using Taylor expansion. Using finite difference scheme typically requires taking high-order derivatives of jump conditions and interface in Taylor expansion. Also property of the discretized linear system is hard to analyze for interface problem with general jump condition. The “immersed interface” method was proposed in [12]. This method incorporates the interface conditions into the finite difference scheme near the interface to achieve second-order accuracy based on a Taylor expansion in a local coordinate system. Second-order differentiation of the interface is needed. The Matched Interface and Boundary (MIB) method was developed in [13] and improved to handle sharp-edged interfaces in two dimensions in [14] and in three dimensions in [15]. With an elegant treatment, second-order accuracy was achieved in the norm.

The existing finite element schemes on unfitted meshes are usually designed by modifying the finite element basis near the interface. Examples are immersed finite element method [16, 17], adaptive immersed interface method [18], and extended finite element method [19–21].

An elasticity system can be solved by both finite difference and finite element methods. In [22, 23], a finite difference method called immersed interface method is proposed to solve elasticity problems with nonhomogeneous jump conditions. While second-order accuracy was achieved, the condition number of the discrete system is quite large especially in the nearly incompressible case ( is large) compared with that obtained from finite element formulations.

In [17], an immersed-interface finite element method was developed for scalar elliptic interface problems with nonhomogeneous jump conditions. In [24], the immersed-interface finite element method was developed for the elasticity system with homogeneous and nonhomogeneous jump conditions. The Soblev space theory provides strong theoretical foundations for convergence analysis for finite element methods.

There are some other approaches to solve the elliptic interface problems. In particular, some recent work [14, 25, 26] can handle sharp-edged interfaces. However, these have not been developed to solve the elasticity interface problems.

In [27], a nontraditional finite element method is proposed to solve the 2D elasticity interface problems with sharp-edged interfaces.

In this paper, based on the method in [27–30], we propose a numerical method for solving the 3D elasticity problem with sharp-edged interfaces. Compared with [27], this is not a trivial extension, as the geometry in three dimensions is more complicated. In two dimensions, the interface could only cut a triangle in one way (if not counting the case the interface hits a grid point as a separate case). In three dimension, There are two types of plane segments see Figure 3. We proved that the resulting linear system is nonsymmetric but positive definite under certain assumption. The method is simpler compared with that developed in [24] and can be applied for more general problems since we allow to be matrices. The condition number grows with , the same as the case without interface. The condition number also grows only linearly with the ratio between and .

The rest of the paper is organized as follows: Section 2 presents the weak form of the elasticity system. In Section 3, we discuss our new numerical method and some theoretical analysis. In Section 4, some examples are presented to demonstrate the performance of our method. We conclude in Section 5.

2. The Weak Formulation

We modify the weak formulation in [28, 30, 31]. We are going to use the usual Sobolev spaces . For , instead of the usual inner product, we choose one which is better suited to our problem:

Let be any closed Lipschitz continuous hypersurface of dimension in , where the overline denotes the closure of a set. Let denote the restriction operator from to . This restriction operator is well defined and bounded, because it is closed Lipschitz continuous (see Theorem  2.4.2 in [32]) and is dense in . Throughout this section, we shall always assume that our interface data and are the restrictions of functions and on and then limited on , respectively. That is on , We shall always assume that our boundary data can be obtained: assume that there exists a function so that is given as, on ,

Equation (6) could be thought as a compatibility condition between and . To simplify the notation, from now on, we will drop the tildes.

We will construct a unique solution of the problem in the class

If , then and . Note that can be identified with .

Definition 1. A function is a weak solution of (1)–(3), if satisfies for all , where or equivalently.

Definition 2. A function is a weak solution of (1)–(3), if satisfies, for all , A classical solution of (1)–(3), , is necessarily a weak solution. Because all the subdomains’ boundaries are Lipschitz continuous, the integration by parts is legal in each subdomain, ; see Theorem  1.5.3.1 in [33].

We have the following theorem.

Theorem 3. If , and , , and , then there exists a unique weak solution of (1)–(3) in .

Proof. See Theorem 2.1 in [31].

3. Numerical Method

We define and choose test function and redefine the gradient and divergence operator

Then (1) can be written as The jump condition equation (2) can be reformed as and boundary condition is

For simplicity, we will discuss the following properties under the form of (14), (15), and (16).

For ease of discussion, Section 3, and for accuracy testing in Section 4, we assume that , , and are smooth on the closure of , and are smooth on and , but they may be discontinuous across the interface . However, , , and are kept to be Lipschitz continuous. We assume that there is a Lipschitz continuous and piecewise smooth level-set function on , where , , and . A unit vector can be obtained on , which is a unit normal vector of at pointing from to .

In this paper, we restrict ourselves to a cube cell domain in 3D space, and is a matrix that is uniformly elliptic in each subdomain. Given positive integers , , and , set , , and . We define a uniform Cartesian grid for , and . Each is called a grid point. For the case , , or , a grid point is called a boundary point, otherwise it is called an interior point. The grid size is defined as .

Two sets of grid functions are needed, and they are denoted by

We cut every cube cell region into six tetrahedron regions. Collecting all those tetrahedron regions, we obtain a uniform tetrahedralization ; see Figure 2.

Figure 2: Cubic mesh and the tetrahedralization.

Figure 3: Interface tetrahedron, left: Case 1, right: Case 2.

If , we count the grid point as in ; otherwise we count it as in . We call an edge (an edge of a tetrahedron in the tetrahedralization) an interface edge if two of its ends (vertices of tetrahedrons in the tetrahedralization) belong to different subdomains; otherwise we call it a regular edge.

We call a cell an interface cell if its vertices belong to different subdomains. In the interface cell, we write . and are separated by a plane segment, denoted by . There are two kinds of plane segment; see Figure 3. The vertices of the plane segment are located on the interface , and their locations can be calculated from the linear interpolations of the discrete level-set functions . The vertices of are located in and the vertices of are located in . and are approximations of the regions of and , respectively. We call a cell a regular cell if all its vertices belong to the same subdomain, either or . For a regular cell, we also write , where (empty set) if all vertices of are in and if all vertices of are in . Clearly in a regular cell, and and are approximations of the regions and , respectively. We use and to represent the volumes of and , respectively.

Two extension operators are needed. The first one is . For any , is a standard continuous piecewise linear function, which is a linear function in every tetrahedron cell and matches on grid points. Clearly such a function set, denoted by , is a finite dimensional subspace of . The second extension operator is constructed as follows. For any with at boundary points, is a piecewise linear function and matches on grid points. It is a linear function in each regular cell, just like the first extension operator in a regular cell. In each interface cell, it consists of two pieces of linear functions, one is on and the other is on . The location of its discontinuity in the interface cell is the plane segment . Note that the vertices of the plane segment are located on the interface , and hence the interface condition is enforced exactly at the vertices of . In each interface cell, the interface condition is enforced with the value at the barycenter of for Case 1 and at the points labeled “9” and “10” in Figure 3 for Case 2. Clearly such a function is not continuous in general, and neither is the set of such functions a linear space. We denote the set of such functions as , which should be thought as an approximation of the solution class (see (7)) plus the restriction of . Similar versions of such extension can be found in the literature [16, 34]. In Case 1, the number of unknowns matches with the number of the equations and the extension is unique. In Case 2, the local construction problem is solved in a least square sense. The values of the points on the interface can be represented as a linear combination of a constant and the 8 unknown value of the vertices on the tetrahedron cell.

We propose the following method.

Method 4. Find a discrete function such that on boundary points and so that for all , we have Note that on the boundary is the same as on the boundary.

Since our solution bases and test function bases are different, the matrix for the linear system generated by Method 4 is not symmetric in the presence of an interface. However, we can prove that it is positive definite under certain assumptions. In all our numerical examples, is positive definite numerically.

Theorem 5. If is continuous and positive definite, then the matrix for the linear system generated by Method 4 is positive definite.

Proof. For vector we have Here we can choose specific jump conditions since changing , will not affect the matrix . In this case we have since So if is positive definite, then , for all . Therefore , which implies . Thus, is positive definite.

In some applications [24], the matrix is only semipositive definite with zero determinant. The previous Theorem does not apply. We prove later that when the matrix is of certain form that frequently appeared in applications and semipositive definite, the matrix generated by our method is still positive definite.

Theorem 6. If , and , , , and , then the matrix for the linear system generated by Method 4 is positive definite.

Proof. Suppose for a contradiction that is not positive definite. Then there is a vector and such that . Let Let Then Since for all , So if and only if , , and . Then , for all .However, implies that .Since , without loss of generality, we assume that . If we choose point such that and , , then , a contradiction.Therefore for all ; thus, is positive definite.

Remark 7. A positive definite matrix has positive determinant and is therefore invertible. It also has an factorization where and , and are lower triangular. The linear system can be solved efficiently.

4. Numerical Experiments

Consider the problem on the paralle lepiped domain . is an interface and prescribed by the zero level-set of a level-set function . The unit normal vector of is pointing from to .

In all numerical experiments later, the level-set function , the coefficients , , and , and the solutions, are given. Hence is obtained as a proper Dirichlet boundary condition, since the solutions are given.

All errors in solutions are measured in the norm in the whole domain .

We present three numerical examples to demonstrate the effectiveness of our method.

Example 1. The level-set function , the coefficients , , , and , and the solution , in Example 1 are given as follows: Table 1 shows the error on different grids. The numerical result shows second-order accuracy in the norm for the solution. Table 2 includes the condition number, which grows with order , the same as the case without interface.

Table 1: Example 1.

Err in
Order

0.00664

0.00174
1.93

0.00046
1.90

0.00012
1.90

Table 2: Example 1 with condition number.

Err in
Order
Condition number

0.00944

80.49853

0.00302
1.64

0.00080
1.92

Example 2. The level-set function , the coefficients , , , and , and the solution , in Example 2 are given as follows: Table 3 shows the error on different grids. The numerical result shows more than 1.5th order accuracy in the norm for the solution. The degeneracy is due to the poor regularity of the solution.

Table 3: Example 2.

Err in
Order

0.00641

0.00210
1.61

0.00069
1.61

0.00022
1.63

Example 3. The level-set function , the coefficients , , , and , and the solution , in Example 3 are given as follows: Table 4 shows the error on different grids. The numerical result shows second-order accuracy in the norm for the solution. This type of examples is also presented in [22, 27] in two dimensions. The coefficient matrices for the PDE meet the assumption in Theorem 6.

Table 4: Example 3.

Err in
Order

0.11372

0.03207
1.83

0.00839
1.93

0.00215
1.96

Example 4. The level-set function , the coefficients , , , and , and the solution , in Example 4 are given as follows:

Table 5 shows the error on different grids. The numerical result shows second-order accuracy in the norm for the solution. Compared with Example 3, the difference is that we have variable coefficient matrices of the same structure as those constant coefficient matrices in Example 3.

Table 5: Example 4.

Err in
Order

0.13994

0.04011
1.80

0.01036
1.95

0.00263
1.98

Example 5. The level-set function , the coefficients , , , and , and the solution , in Example 5 are given as follows: Table 6 shows the error on different grids. The numerical result shows second-order accuracy in the norm for the solution. Compared with Example 4, the difference is that we used general variable coefficient matrices without the special structure in Example 4 to demonstrate that our method can handle more general problems.

Table 6: Example 5.

Err in
Order

0.09240

0.02498
1.89

0.00635
1.96

0.00160
1.99

Example 6. The level-set function , the coefficients , , , and , and the solution , in this example are given as follows:

When it is grids, we get the result of Table 7. When it is grids, we get the result of Table 8. We plot the condition number versus the ratio between and for 11-by-11-by-11 grid in Figure 4 and 22-by-22-by-22 grid in Figure 5. The correlation coefficients for 11-by-11-by-11, 22-by-22-by-22 grids are 0.999990 and 0.999987, which are clearly linear relations.

Table 7: Example 6. grids.

Err in
Condition number

0.07981

0.20388

0.39149

0.54457

0.67088

0.77679

0.86828

Table 8: Example 6. grids.

Err in
Condition number

0.02173

0.08743

0.15283

0.21488

0.27731

0.33511

0.38877

Figure 4: Condition number versus the ratio between and for 11-by-11-by-11 grid.

Figure 5: Condition number versus the ratio between and for 22-by-22-by-22 grid.

5. Conclusion

In this paper, we modified the method in [27] to solve the 3D elasticity problem with sharp-edged interfaces. Due to the more complicated geometry in 3D, this is not a trivial extension. We proved that the matrix for the linear system generated by our method is positive definite (but not symmetric) under certain assumption. Through numerical experiments, our method achieved close to second-order accuracy in the norm for piecewise smooth solutions, and we can handle the difficulties of sharp-edged interfaces and singular solutions. The condition number of the large sparse linear system grows with order , the same as the case without interface. The condition number also grows only linearly with the ratio between and .

Acknowledgment

L. Wang’s research is supported by the Science Foundation of China University of Petroleum, Beijing (no. YJRC-2013-48). S. Hou’s research is supported by NSF grant DMS-1317994.

References

T. Lin, D. Sheen, and X. Zhang, “A locking-free immersed finite element method for planar elasticity interface problems,” Journal of Computational Physics, vol. 247, pp. 228–247, 2013.

H.-J. Jou, P. H. Leo, and J. S. Lowengrub, “Microstructural evolution in inhomogeneous elastic media,” Journal of Computational Physics, vol. 131, no. 1, pp. 109–148, 1997.

H. Gao, Y. Huang, and F. F. Abraham, “Continuum and atomistic studies of intersonic crack propagation,” Journal of the Mechanics and Physics of Solids, vol. 49, no. 9, pp. 2113–2132, 2001.

P. Grisvard, Elliptic Problems in Nonsmooth Domains—Monographs and Studies in Mathematics, Pitman Advanced Publishing Program, 1985.

A. Quarteroni, Numerical Models for Differential Problems, Springer, Berlin, Germany, 2009.

A. Quarteroni and R. Sacco, Numerical Approximation of Partial Differential Equations, Springer, Berlin, Germany, 1994.

R. H. Nochetto, M. Paolini, and C. Verdi, “An adaptive finite element method for two-phase stefan problem in two space dimensions, part II: implementation and numerical experiments,” SIAM Journal on Scientific and Statistical Computing, vol. 12, pp. 1207–1244, 1991.

A. Schmidt, “Computation of three dimensional dendrites with finite elements,” Journal of Computational Physics, vol. 125, no. 2, pp. 293–312, 1996.

C. S. Peskin, “Numerical analysis of blood flow in the heart,” Journal of Computational Physics, vol. 25, no. 3, pp. 220–252, 1977.

C. S. Peskin, “Improved volume conservation in the computation of flows with immersed elastic boundaries,” Journal of Computational Physics, vol. 105, no. 1, pp. 33–46, 1993.

M. Sussman, “A level set approach for computing solutions to incompressible two-phase flow,” Journal of Computational Physics, vol. 114, no. 1, pp. 146–159, 1994.

R. J. Leveque and Z. Li, “Immersed interface method for elliptic equations with discontinuous coefficients and singular sources,” SIAM Journal on Numerical Analysis, vol. 31, no. 4, pp. 1019–1044, 1994.

Y. C. Zhou, S. Zhao, M. Feig, and G. W. Wei, “High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources,” Journal of Computational Physics, vol. 213, no. 1, pp. 1–30, 2006.

S. Yu, Y. Zhou, and G. W. Wei, “Matched interface and boundary (MIB) method for elliptic problems with sharp-edged interfaces,” Journal of Computational Physics, vol. 224, no. 2, pp. 729–756, 2007.

S. N. Yu and G. W. Wei, “Three-dimensional matched interface and boundary (MIB) method for treating geometric singularities,” Journal of Computational Physics, vol. 227, pp. 602–632, 2007.

Z. Li, T. Lin, and X. Wu, “New Cartesian grid methods for interface problems using the finite element formulation,” Numerische Mathematik, vol. 96, no. 1, pp. 61–98, 2003.

Y. Gong, B. Li, and Z. Li, “Immersed-interface finite-element methods for elliptic interface problems with nonhomogeneous jump conditions,” SIAM Journal on Numerical Analysis, vol. 46, no. 1, pp. 472–495, 2007.

Z. Chen, Y. Xiao, and L. Zhang, “The adaptive immersed interface finite element method for elliptic and Maxwell interface problems,” Journal of Computational Physics, vol. 228, no. 14, pp. 5000–5019, 2009.

N. Moës, J. Dolbow, and T. Belytschko, “A finite element method for crack growth without remeshing,” International Journal for Numerical Methods in Engineering, vol. 46, no. 1, pp. 131–150, 1999.

G. J. Wagner, N. Moës, W. K. Liu, and T. Belytschko, “The extended finite element method for rigid particles in Stokes flow,” International Journal for Numerical Methods in Engineering, vol. 51, no. 3, pp. 293–313, 2001.

J. Chessa and T. Belytschko, “An extended finite element method for two-phase fluids,” Journal of Applied Mechanics, vol. 70, no. 1, pp. 10–17, 2003.

X. Yang, B. Li, and Z. Li, “The immersed interface method for elasticity problems with interfaces,” Dynamics of Continuous, Discrete and Impulsive Systems A, vol. 10, no. 5, pp. 783–808, 2003.

X. Yang, Immersed interface method for elasticity problems with interfaces [Ph.D. thesis], North Carolina State University.

Y. Gong and Z. Li, “Immersed interface finite element methods for elasticity interface problems with non-homogeneous jump conditions,” Numerical Mathematics, vol. 3, no. 1, pp. 23–39, 2010.

M. Oevermann, C. Scharfenberg, and R. Klein, “A sharp interface finite volume method for elliptic equations on Cartesian grids,” Journal of Computational Physics, vol. 228, no. 14, pp. 5184–5206, 2009.

J. Bedrossian, J. H. von Brecht, S. Zhu, E. Sifakis, and J. M. Teran, “A second order virtual node method for elliptic problems with interfaces and irregular domains,” Journal of Computational Physics, vol. 229, no. 18, pp. 6405–6426, 2010.

S. Hou, Z. Li, L. Wang, and W. Wang, “A numerical method for solving elasticity equations with interfaces,” Communications in Computational Physics, vol. 12, no. 2, pp. 595–612, 2012.

S. Hou, W. Wang, and L. Wang, “Numerical method for solving matrix coefficient elliptic equation with sharp-edged interfaces,” Journal of Computational Physics, vol. 229, no. 19, pp. 7162–7179, 2010.

S. Hou, L. Wang, and W. Wang, “A numerical method for solving the elliptic interface problems with multi-domains and triple junction points,” Journal of Computational Mathematics, vol. 30, no. 5, pp. 504–516, 2012.

S. Hou, P. Song, L. Wang, and H. Zhao, “A weak formulation for solving elliptic interface problems without body fitted grid,” Journal of Computational Physics, vol. 249, pp. 80–95, 2013.

S. Hou and X.-D. Liu, “A numerical method for solving variable coefficient elliptic equation with interfaces,” Journal of Computational Physics, vol. 202, no. 2, pp. 411–445, 2005.

J. Necas, Introduction to the Theory of Nonlinear Elliptic Equations, Teubner-Texte zur Mathematik, Band 52, 1983.

X. U.-D. Liu and T. C. Sideris, “Convergence of the ghost fluid method for elliptic equations with interfaces,” Mathematics of Computation, vol. 72, no. 244, pp. 1731–1746, 2003.

X.-D. Liu, R. P. Fedkiw, and M. Kang, “A boundary condition capturing method for Poisson’s equation on irregular domains,” Journal of Computational Physics, vol. 160, no. 1, pp. 151–178, 2000.