NONCONFORMING SPECTRAL ELEMENT METHOD FOR ELASTICITY INTERFACE PROBLEMS

Abstract

An exponentially accurate nonconforming spectral element method for elasticity systems with discontinuities in the coefficients and the flux across the interface is proposed in this paper. The method is least-squares spectral element method. The jump in the flux across the interface is incorporated (in appropriate Sobolev norm) in the functional to be minimized. The interface is resolved exactly using blending elements. The solution is obtained by the preconditioned conjugate gradient method. The numerical solution for different examples with discontinuous coefficients and non-homogeneous jump in the flux across the interface are presented to show the efficiency of the proposed method.

1. Introduction

The elliptic interface problems arises in many engineering problems, for example, in heat conduction or elasticity problems where the domain of definition is composed of different materials. In this paper we study the nonconforming spectral element method for elasticity interface problems. These problems have wide applications in continuum mechanics, multi-phase elasticity problems, etc.

There exists several methods in the literature to solve elliptic interface problems (see [26]). There are two types of finite element methods for elliptic interface problems: fitted and unfitted finite element methods. Fitted finite element method is a common approach where the mesh is fitted to the interface, so the interface conditions are satisfied in the weak formulation. The interface is approximated by the sides of isoparametric elements in the discretization. The accuracy of the method depends on the approximation quality of the interface. In this case, the method converges with optimal rates in h. If the mesh is not fitted to the interface, suboptimal convergence behavior will occur or may be the method does not converge at all. To avoid this difficulty, in [2] Babuska have formulated an equivalent minimization problem with all boundary and jump conditions incorporated in the cost functions.

Unfitted finite element methods are based on a mesh which is independent of interface. In [3] unfitted finite element method based on penalized problem as in [2] has been proposed. With an appropriate choice of penalty term the approximation converges to the solution at optimal rate (in h) in H1 norm. Only lower order finite elements have been studied in the literature (for more details see [27]). In [27], a conforming higher order finite element method has been analyzed for elliptic interface problems. One can also look for different formats of finite element methods for interface problems in [24,33,35,41]. Many nonconforming approaches are available in the literature like discontinuous Galerkin methods, Mortar finite element methods, etc (see [26]). In [8] a priori and a posteriori error estimates have been derived for discontinuous Galerkin method. Quasi optimal a priori estimates for interface problems even with lower smoothness conditions on the solution were derived.

This problem is also studied in the framework of least-squares finite element method [4,5,6,11]. In these formulations, the given differential equation is converted into first order partial differential system and a suitable least-squares formulation is applied. Optimal convergence rates in h have been shown. The first order system least-squares method (FOSLS) for linear elasticity problems has been proposed in [9,10]. Least-squares spectral element method has been proposed in [19,20]. In [31,32] iterative substructuring methods for spectral element discretizations of elliptic systems have been proposed. The method provides an efficient preconditioner with an optimal condition number. The extended finite element method (XFEM) or generalized finite element method (GFEM) is a useful method for approximating the solutions with singularities and solutions of interface problems. This method extends the FEM approach by enriching the solution space. The approximation consists of standard finite element approximation and the enrichment through the partition of unity concepts [17,18].

Immersed interface method has been widely studied for elliptic interface problems [28]. Finite difference based explicit jump immersed interface method for elasticity systems was described in [25]. An immersed finite element method for elasticity equations with interfaces has been studied in [29,39,40]. In this method the mesh is independent of the interface and basis functions are chosen such that they satisfies the interface conditions. Optimal convergence rates in h have been derived. In [30], linear and bilinear immersed finite elements for planar elasticity interface problems have been discussed. 2D linear, bilinear immersed finite elements which satisfy the interface jump conditions were used. Optimal convergence rates in h were shown in L2 and semi H1 norms. Details and complete citation list on immersed finite element methods can be found in [38]. In [12] an adaptive immersed interface finite element method for elasticity interface problem was presented.

In [22] Nitsche’s method has been described. In [23] a finite element solution of elliptic interface problem using an approach due to Nitsche has been proposed. The method allows for discontinuities, internal to the elements, in the approximation across the interface and it was shown to be second order accurate (in h) in L2 norm. In this method the interface conditions are satisfied weakly by means of variant of Nitsche’s method. A hp Nitsche’s method for interface problems with nonconforming unstructured finite element meshes have been proposed and error estimates with optimal bound in h and suboptimal bound in p by degree p1/2 were obtained in [13].

In this paper, we propose a least-squares spectral element method for elasticity interface problems based on the method proposed in [26]. In the least-squares formulation of the method, a solution is sought which minimizes the sum of the squares of a squared norms of the residuals in the partial differential equation and the sum of the residuals in the boundary conditions in fractional Sobolev norms and the sum of the jumps in the displacement and the flux across the interface in appropriate fractional Sobolev norms and enforce the continuity along the inter element boundaries by adding a term which measures the sum of the squares of the jump in the function and its derivatives in fractional Sobolev norms. The proposed numerical formulation is based on the regularity estimate for the interface problems stated in [7] and the stability estimate proved in [26].

This method is nonconforming (in terms of approximation). This formulation is different from the standard techniques in LSFEM used to convert the second order elliptic equations into first order system. The interface is resolved completely using blending elements [21]. Higher order spectral elements are used to approximate the solution. The spectral elements are the sum of tensor products of the polynomials of degree W in each variable. The solution is obtained using preconditioned conjugate gradient method (PCGM) without storing the stiffness matrix and load vector. Even though we do not store the matrix, the added advantage of the proposed method is the resulting stiffness matrix is symmetric and positive definite. The integrals involved in the residual computations are obtained efficiently and inexpensively [36] (a brief description is given in the Appendix).

The rest of the paper is organized as follows: In Section 2 the elasticity interface problem is defined. The discretization of the domain is given in Section 3 and the numerical scheme is described. Finally in Section 4 numerical results are presented for various examples.

 

2. Elasticity interface problem

In this section we state the elasticity interface problem on a domain Ω ⊆ ℝ2. First we define the function spaces which we need in the latter sections.

Denoting Hk(Ω), the usual Sobolev space of integer order k with the norm ∥.∥k,Ω as given below,

Further, let

denote the fractional Sobolev norm of order s, where 0 < s < 1. Here J denotes an interval contained in ℝ.

We denote vectors and vector spaces by bold characters. For example, u = (u1, u2)T, Hk(Ω) = Hk(Ω) × Hk(Ω), etc. The norms are given by for etc.

2.1. Linear elasticity system. Let x = (x, y) be a point in space, u = (u1(x, y), u2(x, y))T be the displacement vector and ϵ = (ϵij) be the strain tensor. If u1, u2 are the two planar displacement components, then the straindisplacement relation is given by

The relation between stresses and strains (from the Hooke’s law) is given by,

where λ and μ are the Lame coefficients, and

Let σ = (σij) be the stress tensor, f (x) = (f1, f2)T be the applied body forces, then the stress tensor satisfies the following partial differential equations,

From the above equations, we can re-write the above system as the system of plane elasticity equations of the following,

The Lame coefficients λ and μ are given by

where E is the Young’s modulus and ν is the Poisson’s ratio. Let us assume that the constants λ and μ have finite jumps across the interface; so does the flux σn. Now the elasticity interface problem is defined below.

2.2. The interface problem. Let Ω and be open bounded domains with boundaries respectively. Assume that the boundary Γ0 is sufficiently smooth. Further, let Let Now the elasticity interface problem can be written as follows:

where Lu = (L1u,L2u)T , f = (f1, f2)T, q = (q1, q2)T, g = (g1, g2)T are known vector functions. n is the unit outward normal to the interface Γ0. The jump [.] is defined as the difference of the limiting values from the outside of the interface to the inside. The coefficients λ and μ are piecewise constant, i.e.,

 

3. Discretization and Numerical Scheme

Considered the circular domain Ω1 such that where Ω is square whose boundary is as shown in Figure 1, for brevity. Let and the interface is Γ0 which is smooth as shown in Figure 1. The results presented are applicable to arbitrary smooth interfaces also.

FIGURE 1.The domain Ω and discretization.

Now the domain Ω1 and Ω2 are partitioned into finite number of quadrilateral subdomains (elements) such that the subdomain divisions match on the interface. The interface is completely resolved using blending elements [21].

Each element is mapped to the master square S = (−1, 1)2. Define an analytic map from the master square by (see [1,21])

A brief description of the map is given in the Appendix A1. Here and in the rest of this section l = 1, ..., p for i = 1 and l = 1, ..., q for i = 2.

Define the spectral element functions1 as the tensor product of polynomials of degree W in each variable ξ and η as

Then2 are given by

Now

Here is the Jacobian of the mapping from

Define3 (the differential operator in the transformed coordinates in the domains Ω1 and Ω2 respectively). Then

Define Let for i = 1, 2. Define4

Let γs be a side common to the two adjacent elements (as shown in Fig. 2(a) for i = 2). Assume that γs is the image of η = −1 under the mapping which maps and also the image of η = 1 under the mapping which maps By chain rule

FIGURE 2.Elements with common edges

Then the jumps along the inter-element boundaries are defined as

Here and in what follows, I is an interval (−1, 1).

As the division of the domain into subdomains match along the interface, we define the jump across the interface by taking it (a part of interface) as the common edge. Consider the elements (as shown in Fig. 2(b)) which have the common edge γs ⊆ Γ0. Let γs be the image of ξ = 1 under the mapping which maps and also the image of ξ = −1 under the mapping which maps Define

where are the tangential derivatives of respectively.

Now along the boundary (for some j) be the image of ξ = 1 under the mapping which maps Then

As defined earlier u1 = u |Ω1 and u2 = u |Ω2, so the boundary condition u = g on Γ in the discrete form will be u2 = g on be the image of the mapping corresponding to the side ξ = 1 and

where −1 ≤ η ≤ 1.

On the interface Γ0 we have [u] = 0 and [σn] = q. Let γs ⊆ Γ0 be the image of ξ = 1 under the mapping which maps and also the image of ξ = −1 under the mapping which maps Let

Let the space of spectral element functions. Define the functional

The approximate solution is chosen as the unique which minimizes the functional over all

The minimization problem leads to a symmetric and positive definite lineear system AZ = b. Where Z be a vector assembled from the values of at the Gauss-Lobatto-Legendre points arranged in lexicographic order for 1 ≤ k ≤ p, 1 ≤ l ≤ q. The solution is obtained by preconditioned conjugate gradient method. The action of a matrix on a vector in each iteration is obtained efficiently and inexpensively without storing the matrix A (since PCGM requires the action of a matrix on a vector). The details are shown in Appendix A2.

We used a preconditioner which was proposed in [15]. The preconditioner5 is block diagonal matrix, where each diagonal block corresponds to the H2 norm of the spectral element function representation of each component of the vector on a particular element which is mapped onto the master square S. The obtained solution of the preconditioned system is nonconforming. A set of corrections are made to the solution so that the corrected solution is conforming (see A3).

Let u be the exact solution and z be the approximate solution which is conforming. Let e = u − z. Then for W large enough we have the following error estimate H1 norm (since the jump in displacement across the interface is zero, u ∈ H1(Ω))

holds, where C and b are constants. Proof is very similar to the one proven in [26].

 

4. Numerical Results

To prove the effectiveness of the method we present the numerical results for the problem defined in Section 3. The relative error ∥e∥ER is defined as In all the examples, degree of the approximation polynomial is denoted by 6 W, ’DOF’ means the number of degrees of freedom and ’Iters’ means the total number of iterations required to compute the solution using PCGM. We have used the relative residual norm as a stopping criteria in PCGM. That is, the iteration process is stopped when the relative residual norm (ri is the residual in ith iteration, ∥.∥2 is vector norm) is less than the tolerance ϵ.

Example 1. Interface problem with homogeneous jump conditions: Consider the linear elasticity interface problem (5) (plane strain, see (4)) stated in Section 2 on a square domain [−0.75,0.75]2 with a circle centered at the origin of radius s as the interface as shown in Fig. 1. The coefficients λ, μ are given by

Chosen the data such that the given interface problem has the exact solution u = (u1, u2)

where Here we choose the radius of the circle Note that the solution u satisfies homogeneous jump conditions across the interface.

This problem have been studied in [39]. We discretized the domain into 9 quadrilateral elements as shown in Fig. 1. The conforming numerical solution has been obtained for various values of b (see (8)) for different degree of the approximating polynomial W. The relative error ∥e∥ER in percent, the iterations are tabulated for different values of b (see (8)) in Table 1 and Table 2.

TABLE 1.The relative error in percent and iterations for different W

TABLE 2.The relative error in percent and iterations for different W

For smaller values of b (see (8)) the iteration count is less but the iteration count is large when bb is large. More efficient preconditioner is under investigation. The log of relative error against the degree of the approximating polynomial W is drawn in Fig. 3 for b = 50, 100. The relation is almost linear. This shows the exponential accuracy of the method. The error decays exponentially for all values of b.

FIGURE 3.Log of the relative error against W for b = 50, 100

Example 2. Interface problem with nonhomogeneous jump in the flux: Consider the linear elasticity problem as defined in (5) (plane strain problem) on the same domain as considered in Example 1 with a circle centered at the origin of radius s as the interface. The coefficients λ, μ are given by

Chosen the data such that the given interface problem has the exact solution u = (u1, u2)

where Here we choose the radius of the circle

For any b, the solution u is continuous across the interface and [σn] = q = (q1, q2),

The domain is discretized as in the above example and the conforming solution is obtained for b=10, 0.1. The relative error ∥e∥ER in percent, the iterations are tabulated in Table 3. The log of relative error against W is drawn in Fig. 4 for b = 10, 0.1. Here one can see that the iteration count is large compared to the count in previous example (look at Table 1 for b = 10).

TABLE 3.The relative error in percent and iterations for different W

FIGURE 4.Log of the relative error against W for b = 0.1, 10

 

5. Conclusions

The proposed method is nonconforming and exponentially accurate. The interface is resolved exactly using blending elements. A small data has to be interchanged in between the elements for each iteration of the PCGM and the residuals in the normal equations can be obtained efficiently and inexpensively. The proposed method is efficient even when the jump in the coefficient is large. The numerical results shows that large differences in the coefficients leads to increase in the number of iterations of the PCGM. A more efficient preconditioner is under investigation. This method is also efficient on parallel computers. The method is applicable to arbitrary smooth interfaces too and the method can be extended to the singular case which is ongoing work.