Matched Interface and Boundary Method for Elasticity Interface Problems

Matched Interface and Boundary Method for Elasticity Interface Problems

Bao Wang, Kelin Xia and Guo-Wei Wei 111 Corresponding author. Email:

Department of Mathematics, Michigan State University, East Lansing, MI 48824, USA
Department of Electrical and Computer Engineering,
Michigan State University, East Lansing, MI 48824, USA
July 4, 2019

Elasticity theory is an important component of continuum mechanics and has had widely spread applications in science and engineering. Material interfaces are ubiquity in nature and man-made devices, and often give rise to discontinuous coefficients in the governing elasticity equations. In this work, the matched interface and boundary (MIB) method is developed to address elasticity interface problems. Linear elasticity theory for both isotropic homogeneous and inhomogeneous media is employed. In our approach, Lam’s parameters can have jumps across the interface and are allowed to be position dependent in modeling isotropic inhomogeneous material. Both strong discontinuity, i.e., discontinuous solution, and weak discontinuity, namely, discontinuous derivatives of the solution, are considered in the present study. In the proposed method, fictitious values are utilized so that the standard central finite different schemes can be employed regardless of the interface. Interface jump conditions are enforced on the interface, which in turn, accurately determines fictitious values. We design new MIB schemes to account for complex interface geometries. In particular, the cross derivatives in the elasticity equations are difficult to handle for complex interface geometries. We propose secondary fictitious values and construct geometry based interpolation schemes to overcome this difficulty. Numerous analytical examples are used to validate the accuracy, convergence and robustness of the present MIB method for elasticity interface problems with both small and large curvatures, strong and weak discontinuities, and constant and variable coefficients. Numerical tests indicate second order accuracy in both and norms.

Keywords:  Elasticity equations; Elasticity interface problems; Spatial-dependent shear modulus; Matched interface and boundary.

1 Introduction

Elasticity interface problems play significant roles in continuum mechanics in which elasticity theory and related governing partial differential equations (PDEs) are commonly employed to describe various material behaviors. For this class of problems, an interface description in the elasticity theory is indispensable whenever there are voids, pores, inclusions, dislocations, cracks or composite structures in materials [5, 9, 24, 23]. Elasticity interface problems are particularly important in tissue engineering, biomedical science and biophysics [28, 29, 31]. In many situations, the interface is not static such as fluid-structure interfacial boundaries [26]. Discontinuities in material properties often occur over the interface. Mathematically, there are two types of discontinuities, namely, strong discontinuities and weak discontinuities. Strong discontinuities are referred to situations where the displacement has jumps across the interface. In contrast, weak discontinuities are concerned with jumps in the gradient of the displacement, whereas the displacement is still continuous. In the linear elasticity, the stress-strain relation is governed by the constitutive equation. For isotropic homogeneous material, constitutive equations can be determined with any two terms of bulk modulus, Young’s modulus, Lam’s first parameter, shear modulus, Poisson’s ratio, and P-wave modulus [1]. If these moduli are position dependent functions, related constitutive equations can be used to describe elasticity property of isotropic inhomogeneous media. In seismic wave equations, inhomogeneity is accounted by assuming Lam’s parameters to be a position dependent function [22]. This model is also used in the elasticity analysis of biomolecules [28, 29, 31].

The study of the analytical solution for elasticity interface problems dated back to Eshelby in 1950s [6, 7]. Working on inclusion and inhomogeneity problems, Eshelby found that for an infinite and elastically isotropic system with an ellipsoidal inhomogeneity, the eigenstrain distribution is uniform inside the inhomogeneity when it is subjected to a uniformly applied stress [6, 7]. Much progress has been made on this area in the past few decades. Recently, semianalytic approaches for finding stress tensors have been proposed for arbitrarily shaped inhomogeneity [19].

Computationally, elasticity interface problems are more difficult than the corresponding Poisson interface problems because of the vector equation and cross derivatives. However, many numerical methods have been designed for elasticity interface problems. Based on meshes used, these methods can be classified into two types, i.e., algorithms relied on body-fitting meshes and algorithms based on special interface schemes. For the first type, meshes are generated to fit to the geometry of the interface without cutting through the interface. Therefore, adaptive meshes with local refinement techniques are frequently employed [34]. In the second type of algorithms, meshes are allowed to cut through the interface and particular schemes are designed to incorporate the interface information into the element shape function or discretization scheme. Immersed interface method (IIM) [15] has been used to solve elasticity interface problems for isotropic homogenous media [35, 12]. In this finite difference based algorithm, a local optimization scheme is designed for irregular grid points and the finial linear equation with a non-symmetric matrix is solved by special solvers like BICG and GRMES. Second order accuracy is obtained [35]. A second-order sharp numerical method has been developed for linear elasticity equations [25]. Finite element based methods are also proposed for elasticity interface problems. Among them, the partition of unity method (PUM), the generalized finite element method (GFEM) and extended finite element method (XFEM) are developed to capture the non-smooth property of the solution over the interface by adding enrichment functions to the approximation [24, 23, 9]. Through the weak enforcement of the continuity, discontinuous Galerkin based methods have been employed to simulate strong and weak discontinuities [13, 2, 20]. Recently, immerse finite element (IFM) method has been proposed to solve elasticity problems with inhomogeneous jump conditions [16, 33, 3]. In this approach, finite element basis functions are adjusted locally to satisfy the jump conditions across the interface. Sharp-edged interface is considered for a special elasticity interface problem [14]. Lin, Sheen and Zhang have proposed a bilinear IFM and further modified it to a locking-free version [18, 17]. For both compressible and nearly incompressible media, this method works well and offers second order accuracy. Recently, immersed meshfree Galerkin method has also been proposed for composite solids [30]. Most recently, a Nitsche type method has been proposed for elasticity interface problems [21]. Given the importance of elasticity interface problems in science and engineering, it is expected that more efficient numerical methods will be developed for this class of problems in the near future.

The matched interface and boundary (MIB) method was originally developed for solving Maxwell’s equations [41] and elliptic interface problems [36, 37, 47, 46, 11]. A unique feature of the MIB method is that it provides a systematic procedure to achieve arbitrarily high order convergence for simple interfaces [41, 47] and second order accuracy for arbitrarily complex interface geometry [36, 37]. The essential idea is to introduce fictitious values at irregular mesh points which form fictitious domains [27] so that standard finite difference schemes can still be used across the interface. The lowest order interface jump conditions are iteratively enforced at the interface which determines fictitious values on fictitious domains. Typically, whenever possible, a high-dimensional interface problem is split into simple one-dimensional (1D) interface problems, similar to our earlier discrete singular convolution algorithm [27]. Due to the great flexibility in the construction of fictitious approximations, the MIB method has been shown to deliver up to 16th order accuracy for simple interfaces [41, 47] and robust second order accuracy for arbitrarily complex interface geometry with geometric singularities (i.e., non-smooth interfaces with Lipschitz continuity) [36, 37] and singular sources [11]. In the past decades, MIB method has been applied to a variety of problems. In computational biophysics, an MIB based Poisson-Boltzmann solver, MIBPB [4], has been constructed for the analysis of the electrostatic potential of biomolecules [36, 11, 44], molecular dynamics [10] and charge transport phenomenon [42, 43]. Zhao has developed robust MIB schemes for the Helmholtz problems [40, 39]. A second order accurate MIB method is constructed by Zhou and coworkers to solve the Navier-Stokes equations with discontinuous viscosity and density [45]. Recently, the MIB method has been used to solve elliptic equations with multi-material interfaces [32].

The objective of the present paper is to introduce the MIB method for solving elasticity interface problems. We consider both strong and weak discontinuities for isotropic homogeneous and inhomogeneous media. Computationally, the cross derivative terms in the elasticity model give rise to a new challenge for the MIB method when the interface geometry is complex. To overcome this difficulty, we modify the tradition fictitious definition and redefine fictitious values. With the MIB dimension splitting technique, a new fictitious representation is generated for each irregular mesh point based on elastic jump conditions and local geometry. Secondary fictitious values are constructed by the interpolation of these fictitious values and function values. We have designed schemes to deal with both small curvature and large curvature for complex interface geometries. To validate our method, analytical tests for different types of discontinuities and interface geometries are constructed. We demonstrate the second order accuracy of our MIB schemes for elasticity interface problems.

The rest of this paper is organized as follows. The basic setting of elasticity interface problems is presented in Section 2. The linear elasticity equations, interface jump conditions and constitutive laws are discussed in detail to facilitate further consideration. Section 3 is devoted to the construction of MIB algorithms. General fictitious schemes are proposed for elasticity interface problems. Secondary fictitious values are introduced for cross derivative terms. Our method is extensively validated by analytical tests with complex interface geometries in Section 4. This paper ends with a conclusion.

2 Formulation of the elasticity interface problem

In this section, the elasticity problem with material interfaces is formulated. First, the governing equations of the linear elasticity interface problem are derived. Then the weak solution to the governing equation is defined. Based on the weak solution, the interface conditions are derived for the linear elasticity interface problem. Finally, the Dirichlet boundary condition is employed to make the linear elasticity interface problem computationally well-posed.

2.1 Linear elasticity equations

When solid objects are subjected to external or internal loads, they deformed and lead to stress. If the deformation of the solid is relatively small, linear relationships between the components of stress and strain are maintained. Consequently, linear elasticity theory is valid. In practice, linear elasticity theory is applicable to a wide range of natural and engineering materials, and thus extensively used in structural analysis and engineering design.

Figure 1 illustrates the displacement in a two-dimensional (2D) elastic motion. The displacement under an infinitesimal perturbation in position can be approximated by the linear term


where is the position of a point of the un-deformed elastic body, is the th component of the displacement vector and is the relative displacement

Figure 1: Deformation in the continuum where two points and are deformed to and respectively.

The Cauchy’s infinitesimal strain tensor, or the strain tensor for simplicity, is defined as


where is the th element of . In a compact notation, the linear strain tensor is given by


where is the position of the point of the deformed elastic body. Obviously, the strain tensor describes the total displacement. Physically, the Hooke’s law states that the strain must lead to stress. Mathematically, the constitutive equation between strain and stress tensors is given by


where is the th element of stress tensor and is the th element of elastic moduli or stiffness tensor , which is a fourth order tensor describing properties of the material. In the constitutive equation, the stress tensor is expressed through the contraction between the strain tensor and the stiffness tensor.

For isotropic homogeneous media, there is no preferred direction in the stiffness tensor. Therefore, the stress-strain relation can be dramatically simplified. We utilize Lamé’s parameter and shear modulus to simplify the constitutive equation as


where is the kronecker function, is the trace of the strain tensor. In a compact notation, the stress tensor is


where is the identity tensor.

In practical applications, one is often interested in the description of elasticity motion. By the Newton’s second law, the motion of elasticity body is governed by


where is the external force on the elastic body. This equation can be more rigorously derived from the variation principle [28]. The static state of the elastic motion is then governed by:


In many applications, Eq. (9) is solved to obtain the deformation under a given force.

2.2 Interface jump conditions

Consider the static state of two-phase elastic body motions in domain . Let us suppose that the two-phase elastic motions are separated by interface , which separates the whole domain into two sub-domains and , i.e., , as shown in Figure 2.

Figure 2: Illustration of a two-phase elasticity interface problem. The interface separates the whole computational domain into two parts and .

2.2.1 Weak solution of homogeneous elasticity equations

In this section, we define the weak solution to homogeneous linear elasticity equations without the body force


For simplicity, we temporally consider the first equation in homogeneous elasticity equations (10), i.e.,

where denotes the first column of the matrix of the stress tensor , and let , thus the first component of the homogeneous elasticity equation can be written as


To define the weak solution to the homogeneous elasticity equations, we choose the test function space to be . Multiplying Eq. (11) by and integrating the obtained equation over the whole domain yield

Applying the integration by parts to the above equation, and note that is of compact support, one has


Similarly, a form similar to Eq. (12) can be derived for the other equation in Eq. (10).

Definition 2.1.

Weak Solution is said to be the weak solution of the equation provided equation (12) holds for all .

2.2.2 Interface jump conditions

In this part, interface jump conditions are formulated for linear elasticity equation (9).

Lemma 2.1.

For 2D inhomogeneous linear elasticity equation (9), if the force term has a potential function , i.e., , then it can be expressed in a conservative form, i.e., homogeneous form.


Without loss of generality, we only prove the statement for the first equation in equation (9), which is:


where is the first column of the tonsorial matrix .

For convenient, we denote , then equation (13) becomes:

thus, there exists another 2 by 2 tensor , where only depends on potential function , such that


Furthermore, Eq. (9) is equivalent to Eq. (14). ∎

Theorem 2.1.

For 2D linear elasticity equations (9), if the source term has a potential function , i.e., , then across the interface, the weak solution should satisfy following interface conditions


where is a vector-valued function, is the difference of quantity “*” across the interface and is the normal direction of the interface.


By the above lemma, there exists another second order tensor , such that . Without loss of generality, we only prove the interface condition for the first equation in equations and denote it as

Since , and the interface is of measure zero, hence , the following equation holds


Integrating by part applied to the first term in the right hand side of Eq. (16) yields

In region , note , therefore:


where means that the value was evaluated by taking the limit from the region .

Similarly, we have

where the minus sign occurs because the outward normal for is the inward normal for . Therefore, we have


Since test function is arbitrary, we have . In the same manner we can show that , where . Therefore we have . According to the above lemma, we also have:

Similarly, we can prove the following result.

Theorem 2.2.

Let be a nd order tensor in , for elasticity equations


where is an -dimensional vector-valued function, and .

If the force term has a potential function , i.e., , then across the interface, the weak solution should satisfy the following interface conditions


where is an -dimensional vector-valued function and is the normal direction of the interface.

Remark 1.

In interface condition , is a given vector-valued function on the interface , which measures the jump of the traction across the interface .

Moreover, physically, we usually enforce a jump conditions to ensure that the material has no fracture in the weak discontinuity setting, i.e., it is continuous across the interface

However, fractures commonly occurs for many materials which is known as the strong discontinuity in the elasticity mechanic literature. Therefore, a known fracture is often applied

where and are fracture components for and , respectively.

2.3 Elasticity interface problem

Based on the above discussion, two dimensional static interface problems can be formulated as:


where with is the displacement field. Generally, if vector does not equal 0, it is called the strong discontinuity. Otherwise, we arrived at the weak discontinuity. Vector is the unit outer normal vector of the interface , and is a given vector valued function on the interface which measures the jump of the traction across the interface . Here is a 2D vector-valued function which denotes the body force on the elastic object, and is also a 2D vector-valued function to determine the Dirichlet boundary condition.

For isotropic elasticity problems, the stress-strain relation is given by



is the linear strain, and and are two Lamé’s parameters, which can be constants or spatially dependent functions. These two situations are discussed below.

2.3.1 General formulation for homogeneous media

A special case is that, all the material parameters are constants or piecewise constants, i.e., Lamé’s parameters and , Young modulus , and Poisson’s ratio are all constants or piecewise constants. In this case, the above parameters satisfies the following relationships

Particularly, we assume that shear modulus and Poisson’s ratio are given, respectively, by:


Usually, Poisson’s ratios satisfy constraints and . When they are close to , the material is near incompressible. Otherwise, the material is compressible.

The governing equations for the linear elasticity motion in the homogeneous media is given by:


with the interface jump condition defined on


2.3.2 General formulation for inhomogeneous media

Spatially dependent Lamé’s parameters frequently occur in many practical applications. The spatial dependence can be described as


Linear elasticity motion in inhomogeneous media is governed by:


Here these two equations are defined in the domain .

Remark 2.

As discussed earlier, the non-fracture conditions can be relaxed in both theoretical modeling and numerical analysis




3 Methods and algorithms

In this section we develop the second order MIB method for the elasticity interface problem with irregular interface. We consider a rectangular domain . Let be the grid size in griding the rectangular domain and suppose the grid points to be:

Here and are the total numbers of grid points in the - and -directions, respectively. In the standard central finite difference(CFD) discretizations scheme, due to the existence of interface and possibly discontinuous of the solution, the direct use of the CFD scheme may decrease the numerical accuracy, which cannot guarantee the second order convergence of the numerical solution. The main idea of the MIB method, is to replace the referred function values that from the different side of the interface in the (CFD) discretizations schemes by the fictitious values, which are the combinations of function values and interface conditions for capturing the interface conditions and discontinuity of the solutions.

In order to handle the interface in the elasticity interface problem, the regular and irregular grid points in the finite difference scheme of the elasticity equation should be distinguished.

Definition 3.1.

A grid point is said to be a regular grid point provided all the grid points referred in the discretization of the elasticity equations are at the same sub-domain as grid point ; otherwise the grid point is called an irregular grid point.

For the irregular grid point , if the function values at the different side of the interface applied to the finite difference scheme, the accuracy of the numerical solution will be reduced. In the MIB finite difference scheme, we replace function values which are referred in the different side of the interface with fictitious values.

In the following part of this section we discuss how to extend the solution across the interface such that the second order MIB discretization works. To clarify, we call the discretizations of ordinary derivatives and as a five-point stencil, and the discretizations of the cross derivatives as a nine-point stencil. We discuss how to determine fictitious values in these two stencils separately.

3.1 General algorithms for fictitious value

In this section, we develop the MIB method for the discretization of the ordinary derivatives in elasticity equations via fictitious values. Here ordinary derivatives refer to first or second order derivatives along a single direction. In contrast, a cross derivative involves at least two directions and requires special treatments near the interface. For simplicity, we only construct the scheme for the constant material parameter case, while the case of spatially dependent parameters can be treated in the same manner.

3.1.1 Fictitious scheme for regular interface

The interface is said to be regular at the grid point if the interface is locally parallel or vertical to the mesh lines that passing through the grid point . For the regular interface, we can determine fictitious values via an iteratively approach, which ensures that the MIB method can be made into arbitrarily high order. A detail description of the MIB procedure has been given for iteratively finding fictitious values for regular interface [47].

In this section, we focus on the method for estimating fictitious values for a five-point stencil in the second order finite difference scheme of elasticity equations.

Figure 3: Illustration of a regular interface. The interface (red line) is vertical to the horizontal mesh line. In the standard second order MIB finite difference scheme, fictitious values are required at grid points and . To determine the fictitious value at and , one needs to compute and . To this end, and are interpolated by the function values at and fictitious values at ; and are interpolated by the function values at and fictitious values at .

As is shown in Figure 3, the interface is locally vertical to the -mesh at point . Here we only discuss the case of a vertical regular interface, while the case of a parallel regular interface can be handled in the same manner.

We employ the second order accurate interpolation scheme to interpolate one-sided values at point . Four fictitious values are involved, i.e., and , here represents the fictitious value for at the grid point , similar for others. In the following, we present the method to represent fictitious values by interface conditions and function values.

First, the outer normal direction , when the interface is locally vertical to the -mesh, the interface conditions (28)-(31) become


The second order interpolation approximation of , , and can be done by using the function and fictitious values at the grid points , and


where are the coefficients in the Lagrangian interpolation and can be obtained by an appropriate algorithm in the literature [8]. The interpolation approximation of , , and can be handled in the same manner.

Plugging the above finite difference approximation into interface conditions (37)-(40) and solving the resulting approximation of interface conditions, one can represent fictitious values at grid points and by function values at grid points , , and , and interface conditions at .

3.1.2 Fictitious scheme for interface with small curvature

We describe the MIB method for the extension of the function values across the interface which is irregular but with relatively small curvatures in this subsection. Interface conditions are given by Eqs. (28)-(31).

Figure 4 illustrates the situation where fictitious values along the -direction are to be found. The case for determining fictitious values along the -direction is similar.

Figure 4: Illustration of determining fictitious values along the -direction. Fictitious values at two yellow grid point are to be found. To this end, derivatives at need to be calculated. Two red points ( and ) are the points referred along the main direction. Six cyan points are auxiliary points used to approximate values at and . Values at can be computed from the fictitious value at and function values at and . Then interpolated values at and are employed to estimate derivatives at .

When the interface is locally not aligned to an -mesh line or -mesh line, a two-step tactic in the MIB approach can be used. First, one takes the derivative along the tangential direction of the interface to introduce one more set of interface conditions. Then one can eliminate some derivatives at the interface that are difficult to compute. This flexibility directly leads to the efficiency in handling versatile and difficult interface geometries. Finally, one obtains the representation of fictitious values according to interface conditions after the aforementioned two steps.

Lemma 3.1.

For a given function defined on the domain , i.e., the domain is separated by the interface into two parts. If the function is piecewise continuous along the interface , then its tangential derivative is continuous along the interface.

Let the tangential direction be . By taking a derivative along the tangential direction, the obtained new set of interface jump conditions along the tangential direction are:




As shown in Figure 4, let the angle between the normal direction and -mesh be , in this case, the normal direction is , and the tangential direction is , where . For the irregular interface, one has and .

If we define a vector as,


The interface condition along the tangential direction (45)-(46) can be represented as,


Furthermore, additional interface conditions along the normal direction can be written as


where is the wave modulus.

In the standard second-order MIB finite difference scheme, two fictitious values occur for each of elasticity equations near the interface. A total of four fictitious values are to be estimated. Since there are six interface conditions (28)-(31) and (45)-(46), one can use them to eliminate two of eight derivatives. The selection of the derivatives should be based on the local geometry of the interface. A main principle is to eliminate the derivatives that are most difficult to compute due to the local geometric constraint.

Lemma 3.2.

The following matrix is of full rank,

provided that and , i.e., the local interface does not parallel or vertical to the local mesh directions.

According to the above lemma, when the local interface is irregular, the third and fourth rows of the above matrix should be used to minus the linear combinations of the first and second rows so as to replace two interface conditions and leave two interface conditions to compute the fictitious values.

There are four ways to eliminate derivatives at the interface as discussed below.

  • The elimination of the derivatives and generates the following combined interface conditions:

  • The elimination of the derivatives and generates the following combined interface conditions:

  • The elimination of the derivatives and generates the following combined interface conditions:

  • The elimination of the derivatives and generates the following combined interface conditions:


For a given irregular grid point, according to the local interface geometry, one of the replaced interface conditions from above should be chosen to represent fictitious values. Here, we discuss the interface shown in Figure 4. In this case, one-sided derivatives and are to be eliminated. Therefore, the first set of interface conditions is employed. Note that the interpolation approximation of one-sided function values and derivatives referred in the first set of interface conditions for is given as: