A Coupling Approach for Linear Elasticity Problems with Spatially Noncoincident Interfaces

# A Coupling Approach for Linear Elasticity Problems with Spatially Noncoincident Interfaces

Pavel Bochev, James Cheung, Max Gunzburger, and Mauro Perego
###### Abstract

We present a new formulation based on the classical Dirichlet–Neumann formulation for interface coupling problems in linearized elasticity. By using Taylor series expansions, we derive a new set of interface conditions that allow our formulation to pass the linear consistency test. In addition, we propose an iterative method to determine the solution of our formulation. We demonstrate in our numerical results that we may achieve the desired piecewise linear finite element error bounds for both nonoverlapping domain decomposition problems as well as for interface coupling problems where the Lamé parameters of the structures differ.

## 1 Introduction

In many important physical applications, two systems are coupled together through a physical interface. It is on this interface that physical quantities, such as stresses, are transferred between the two systems. Mathematically, such systems are modeled by two sets of equations with a set of additional conditions posed at the interface. We refer to these formulations as interface coupling formulations. If the equations are the same on both sides of this interface, the interface coupling formulations fall under a class of problems called non–overlapping domain decomposition formulations.

In the continuous setting, the interface exists only as a single curve that lies between the two problem subdomains; however, in the finite element setting, this is generally not the case. Often, the two subdomains of the interface coupling problem are meshed separately, and the single interface in the continuous setting becomes two distinct interfaces. If the positions of the nodes in the two interfaces are matching, then there exists no problem in determining the solution of the discretized problem. Most often, classical methods for determining solutions to interface coupling problems in the continuous setting may be extended to the finite element formulation in this case. A thorough exposition of classical interface coupling methods may be found in  and .

If the two interfaces are spatially coincident with nonmatching nodes, one may introduce various operators that map values from one interface to the other without much difficulty. In fact, there exists an entire body of literature that discusses solution methods for this kind of discrete interface problem. Among the methods that deal with this kind of nonmatching interface problem, one of the most studied and utilized methods are mortar element methods [1, 12].

However, if the continuous interface is a curved surface, the interfaces generated by separate meshing will often not spatially coincide. The most prominent issue in this setting is that the interface conditions defined in the continuous setting has little meaning since there is no clear definition of an interface condition when the discrete interfaces are spatially noncoincident. When comparing the volume of literature available for solution methods of the case of spatially noncoincident interface coupling, one will find that the amount of available literature on this subject is quite sparse. Most often, the methods discussed in the previous paragraph are generalized to this case.

A desirable property of any discrete approximation to its continuous counterpart is that it is consistent. This means that the behavior of the discrete model should match the behavior of the continuous model. A common approach to test whether a discrete formulation is consistent is whether or not it can reproduce –degree polynomials, where is the polynomial degree used in the approximation space. As of current, no method in the literature is –degree consistent; however, there are a few that are linearly consistent. A novel approach presented in  uses an energy correction approach to remove excess energy that arises from the overlaps generated from the spatially noncoincident interfaces and adds additional energy to the deficit incurred by the gaps between the noncoincident interfaces. The complexity of the implementation of this method may be prohibitive since a mesh–like structure must be constructed between the discrete interfaces to link the two interfaces. The method proposed in  removes this complexity by perturbing the meshes at the discrete interfaces so that the area of the gaps and overlaps in the subdomain meshes are equal.  removes this perturbation requirement by requiring that the subdomain meshes are only overlapping. Because these methods are based on minimizing a variational principle over the entire problem domain, their applicability is limited to cases where the material constant on both sides of the domain are the same.

In this work, we present an interface coupling approach for spatially noncoincident interfaces for the case of coupled linearized elasticity equations. Through the use of Taylor series expansions, we introduce a new set of interface conditions based on the classical Dirichlet–Neumann formulation that allow us to pass the linear patch test in the domain decomposition problem. Subsequently, we introduce an iterative solution method that allows us to solve this formulation. The iterative method is particularly useful if a coupled solution must be determined from separate codes. To ensure that we achieve the expected accuracy from piecewise linear finite element methods, we utilize the Zhang–Naga gradient recovery operator  to recover a second order accurate stress tensor as Neumann data in the coupled formulation. An additional benefit of this new formulation is that, unlike the methods described in the previous paragraph, this method requires no intermediate meshing procedure, nor does it require any sort of mesh perturbation in the neighborhood of the discrete interfaces. In addition, this coupling approach is also applicable to interface coupling problems where the Lamé parameters differ between the coupled subdomains.

The remainder of this paper will be structured as follows. In Sect. 2 we establish the framework of our problem and define our monolithic coupling formulation. In Sect. 3, we introduce an iterative solution method for the coupling formulation. In Sect. 4, we will demonstrate that our coupling formulation achieves the desired finite element error bounds for both interface coupling and domain decomposition problems. And finally, in Sect. 5 we will provide some concluding statements and provide some insight in the future directions of this work. An appendix is provided at the end of this paper to give the reader insight on the gradient recovery method utilized in this work.

## 2 Statement of the Interface Problem

Assume that is a domain; i.e., a simply connected bounded subset of with a Lipschitz continuous boundary. We shall denote the boundary of as . In addition, let denote the virtual displacement of a material characterized by its Lamé parameters and . Then, using standard notation, we define

 ε(u)=12(∇u+(∇u)T)

as the strain tensor, and

 σ(u)=λ tr ε(u)I+2με(u)

as the stress tensor, where denotes the matrix trace operator. Then, the statement of linearized elasticity is given by

 −∇⋅σ(u) =f in Ω u =g on ∂Ω,

where is the external force applied to the structure, and is the initial displacement of the structure on its boundary.

### 2.1 The Continuous Interface Problem

Let us now assume that there exists two domains, , , such that the intersection of their boundaries is a simply connected curve. Let be the interface between , and be the complement of the interface on the boundary . In addition, let be the outward unit normal vector of . Let us denote as the virtual displacement over for a material characterized by its Lamé parameters and . Assume that there is no slip at ; then, virtual displacements is governed by the system

 −∇⋅σi(ui) =fi in Ωi (1a) ui =gi on Γi subject to us−um=0σs(us)nm−σm(um)nm=0on Γc, (1b)

where , is the external force applied to each structure, and is the initial displacement of each structure on . The interface conditions defined on in (1b) describes the transmission of stresses between the two structures. If and , then (1) may be referred to as a nonoverlapping domain decomposition problem because it models the same structural material over both subdomains.

Using the standard Sobolev space notation, let

 H1Γi(Ωi)={v∈H1(Ωi);v=0 on Γi},
 Vi=(H1Γi(Ωi))2,

and

 Vi,0=(H10(Ωi))2

where, again, . In addition, we shall define the interface trace space

 W=(H1/200(Γc))2.

In this work, we consider (1) in the weak form: seek such that

 ∫Ωsσs(us):ϵ(vs)dx =∫Ωsfs⋅vsdx ∀vs∈Vs,0 (2a) ∫Γc(us−um)⋅wds =0 ∀w∈W ∫Ωmσm(um):ϵ(vm)dx−∫Γc(σs(us)nm)⋅vmds=∫Ωmfm⋅vmdx∀vm∈Vm. (2b)

The set of equations (2) is the Dirichlet–Neumann coupling of the linearized elasticity equations, since the matching of the Dirichlet conditions on is weakly enforced in (2a) and the matching of the Neumann conditions on are enforced in (2b).

### 2.2 Finite Element Approximation

Let , , be the triangulation of , with mesh size and be their boundaries. Assume for now that is a simply connected curve; i.e., the interface between the triangulations and is spatially coincident. Let us denote this discretized interface and let . Also, let us define

 Vhi⊂H1Γi(Ωi,h)
 Vhi,0={v∈Vhi;v=0 on ∂Ωi,h},

and

 Wh⊂H1/200(Γch)

as standard piecewise linear approximation spaces. In addition, let , and be the finite dimensional product spaces. The statement of the finite element discretization of (2) becomes the following: seek such that

 ∫Ωs,hσs(uhs):ϵ(vhs)dx =∫Ωs,hfs⋅vsdx ∀vhs∈Vhs,0 (3a) ∫Γch(uhs−uhm)⋅whds =0 ∀wh∈Wh (3b)

where , , denotes the discretized virtual displacement.

#### 2.2.1 Spatially Noncoincident Interfaces

We shall now do away with the assumption that the interfaces between and are spatially coincident. Let , , be the discretization of that belongs to and . In addition, let and again be the the piecewise linear approximation spaces over redefined with the change of definition of . Also again, let and . Now let us define

 Whi⊂H1/200(Γci,h)

as the piecewise linear approximation space over and as the product space over the trace.

The difficulty in coupling problems with noncoincident discrete interfaces presents itself in the interface conditions. Because and are defined over separate interfaces, it becomes a troublesome ordeal to define a suitable set of interface conditions. In this work, we deal with this issue by extending the values related to to by using a Taylor series expansion.

##### Taylor Series Expansion Operators

Before defining our discrete interface conditions we must first define some preliminary material. Let and be points in and respectively, then defining

 h:=maxi=s,mhi,

we assume that there exists mappings

 ϕsm:xs∈Γcs,h→xm∈Γcm,handϕms:xm∈Γcm,h→xs∈Γcs,h

such that

 supxs∈Γcs,hd(xs,ϕsm(xs))≤Chandsupxm∈Γcm,hd(xm,ϕms(xm))≤Ch, (4)

where is an arbitrary constant. To simplify notation, we set and , where the dependence of on and on is implied.

We are now ready to define our Taylor expansion operators. Let be defined as

 Esmv(xs)=v(xs)+∇v(xs)(^xm−xs) (5)

and be defined as

 (6)

where is a rank two tensor . We emphasize that maps functions defined over to functions defined over , since depends on . On the contrary, maps functions defined over to functions defined over , because depends on .

##### A Noncoincident Interface Formulation

Central to our method’s ability to achieve second order accuracy in the –norm is the superconvergent gradient recovery operator. There are various methods in the literature that allow us to recover a superconvergent gradient; i.e., see [15, 14, 7]. In this work, we choose to use the Zhang–Naga gradient recovery operator [13, 8]. We shall denote as the recovered gradient.

Borrowing terminology from the mortar element method, we call the slave interface and the master interface. It is on the master interface, , that we wish to enforce some sort of continuity of the virtual displacements and the normal stresses. Let the extended strain tensor, , be defined as

 ¯¯¯ε(v):=12{Esm(Ghv)+[Esm(Ghv)]T},

where

 Gh[v1v2]=[(Ghv1)T(Ghv2)T]

is the superconvergent Jacobian constructed from the recovered gradient of the vector components of . Using the Taylor expansion operators defined in (5) and (6), we want to enforce

 Esmuhs=uhm(^xm) on Γcs,h, (7a) and σm(uhm)nm=¯¯¯σs(uhs)nm on Γcm,h, (7b)

where is the stress tensor constructed from the extended strain tensor. Again, note that (7a) is an equation posed over despite being defined over due to the dependence of on . Also note that (7b) is an equation posed over . We provide a 1D interpretation of these interface conditions in Fig. 3 to illustrate the meaning of these conditions. Figure 3: Illustrations of the interface conditions (7). Left: Gap at the interface. The dashed line represents the linear extrapolation of us. Right: Overlapping subdomains.
###### Remark.

In this work, we make the fundamental assumption that the master interface, , is an approximation of the continuous interface so that we can achieve optimal convergence rates in the following coupling approach. Otherwise, the error in the interface approximation will dominate the error of the polynomial approximation.

With the interface conditions (7a) and (7b) defined, we may pose the noncoincident interface coupling problem as the following: seek such that

 ∫Ωs,hσs(uhs):ε(vhs)dx =∫Ωs,hfs⋅vhsdx ∀vhs∈Vhs,0 (8a) ∫Γcs,m(Esmuhs−uhm(^xm))⋅whsds =0 ∀whs∈Whs ∫Ωm,hσm(uhm):ε(vhm)dx−∫Γcm,h(¯¯¯σ(uhs)nm)⋅vhmds=∫Ωm,hfm⋅vhmdx∀vhm∈Vhm. (8b)

In the remainder of the paper, we shall use a more compact notation to describe (8). Let the bilinear form , , be defined as

 ahi(v,w)=∫Ωi,hσi(v):ε(w)dx,

then (8) can be represented as: Seek such that

 ahs(uhs,vhs) =⟨fs,vs⟩Ωs,h ∀vhs∈Vhs,0 (9) ⟨Esmuhs−uhm(^xm),whs⟩Γcs,h =0 ∀whs∈Whs ahm(uhm,vhm)−⟨¯σ(uhs)nm,vhm⟩Γcm,h =⟨fm,vhm⟩Ωm,h ∀vhm∈Vhm,

where is the standard duality pairing with the integral taken over , where can be or . The “Dirichlet” subproblem in the above formulation is not a Dirichlet condition in the classical sense; it is perturbed by a oblique derivative term introduced by the Taylor series expansion. The interface condition used in the Neumann subproblem remains a classical natural boundary condition on the interface .

##### Linear Consistency

As stated in the introduction, linear consistency is an important property for a noncoincident coupling method to possess. In the following proposition, we show that (9) unconditionally possesses linear consistency.

###### Proposition 1.

Assume that and . If and are linear vector functions with the same analytic representation, then is a solution of (9).

###### Proof.

Let be any continuous lifting operator that maps into , then we may decompose any function into the following form

 v=v0+Rhv∣∣Γcm,h, (10)

where and . From this, we may represent (9) in the following form: Seek such that

 ⟨−∇⋅σ(qhs),vhs⟩Ωs,h=0∀vhs∈Vhs,0 (11a) ⟨−∇⋅σ(qhm),vh,0m⟩Ωm,h=0∀vhm∈Vhm,0 (11b) ⟨Esmqhs−uhm(^xm),whs⟩Γcs,h=0∀whs∈Whs (11c) ⟨−∇⋅σ(qhm),Rhwhm⟩Ωm,h+⟨σ(qhm)nm−¯¯¯σ(qhs)nm,whm⟩Γcm,h=0∀whm∈Whm, (11d)

where we have used the decomposition (10), and Green’s identity on the bilinear forms . The interface condition (11c) is satified due to the well–known linear polynomial preserving property of the Taylor expansion. It is clear that (11a) and (11b) are satisfied, since the second derivatives of all linear polynomials vanish. It remains to show that (11d) is satisfied. Seeing that from the preceding argument, (11d) becomes

 ⟨σ(qhm)nm−¯¯¯σ(qhs)nm,whm⟩Γcm,h=0∀whm∈Whm, (12)

The polynomial preserving property of the Zhang–Naga gradient recovery operator (cite paper and lemma) implies that the second derivatives in vanish; thus, (12) must be satisfied since and are linear vector functions with the same analytic representation. ∎

###### Remark (Second Order Accuracy).

The second order accuracy of (9) is attributed to the superconvergent stress tensor constructed using the Zhang–Naga gradient recovery operator in conjunction with the use of the second order Taylor series expansion to approximate the continuous interface condition on the slave and master interfaces.

## 3 An Iterative Solution Method

We shall now introduce a modified Dirichlet–Neumann iterative method to solve (9). Let be a relaxation parameter, then given an initial guess , for ,

 ahs(uhs,k,vhs) =⟨fs,vs⟩Ωs,h ∀vhs∈Vhs,0 (13a) ⟨uhs,k,whs⟩Γcs,h =⟨gk,whs⟩Γcs,h ∀whs∈Whs =⟨fm,vhm⟩Ωm,h ∀vhm∈Vhm gk+1=ω(uhm,k(^xm)−∇uhs,k(^xm−xs))+(1−ω)gk. (13b)

For now, we shall assume that (13) is convergent. We now show that, under this assumption, is the solution of (9).

For simplicity, let us denote as the converged solution of the iterative method (13). We see that in the converged limit, the update condition (13b) may be represented as

 limk→∞gk=uhm(^xm)−∇uhs(^xm−xs).

Taking and then inserting this limit into the second equation in (13) yields the interface condition

which is exactly the interface condition on found in the monolithic formulation (9). Hence, the limit of this iterative solution method solves the monolithic problem under the assumption that (13) is convergent.

We remark that the convergence of the (13) is dependent on the Lamé parameters over the slave and master domains along with the choice of the relaxation parameter . In particular, if is significantly larger than , the iterative method becomes nonconvergent regardless of the choice of .

###### Remark.

The convergence of the iterative solution method to the monolithic formulation implies that there exists a nonsingular solution branch of (9).

### 3.1 Implementation

In this subsection, we discuss some details pertaining to the computational implementation of (13). The key topics to be considered here are the mappings between the slave and master interfaces, treatment of the oblique derivative boundary condition, and the construction of the superconvergent gradient at the interface.

##### Defining ϕsm and ϕms

Recall that and are mappings that maps to . In the computational setting, and need only be defined over a finite number of points on each interface. One may choose to define these mappings for the nodes at the interface if it is desired to enforce interpolated boundary conditions; or alternatively, one may choose to define these mappings on the quadrature points of the interface if weak enforcement of the interface condition is desired. In this work, we choose to use the nearest–neighboring node approach to define and on the nodes of their respective ranges, i.e., and respectively. The definition of this approach is defined in the following paragraph.

Let and be arbitrary nodes of and respectively. Then for every node , its nearest–neighboring node is defined as

 ϕsm(zs)=argminzm∈Γcm,h d(zs,zm). (14)

Similarly, for every node , its nearest–neighboring node is defined as

 ϕms(zm)=argminzs∈Γcs,h d(zm,zs). (15)

If , , are approximations of the continuous interface , then it is clear that the distance condition (4) is satisfied using the nearest–neighboring node approach. Of course, this approach may yield greater than one nearest node. In such a case, it is sufficient to just chose a single node out of all possible node candidates that satisfies the minimum distance criterion.

###### Remark.

In the literature, there are methods that allows one to get a closer neighboring point on the neighboring interface. Some methods are detailed in .

##### The Dirichlet Subproblem

In (13), the Dirichlet interface condition on is enforced weakly. While this provides desirable mathematical properties (i.e., the Dirichlet data may be less regular), in practice, one may choose to instead enforce this condition in the strong form, i.e.,

 uhs,k∣∣Γcs,h=gk. (16)

In our numerical experiments, this is the condition that we enforce. Of course, we must now introduce some sort of operator that maps the term in the update condition of (13) into , which is again the piecewise linear approximation space over . In this work, we choose to interpolate the values of this term over the nodes of the . Subsequently, the update condition becomes

 gk+1=ωIs,h(uhm,k(^xm)−∇uhs,k⋅(^xm−xs))+(1−ω)gk∀whs∈Whs, (17)

where is the piecewise linear interpolation operator over . This interpolation is implied in the standard method of enforcing Dirichlet conditions in the finite element method.

##### The Neumann Subproblem

At first glance, it may appear that we must define the superconvergent Jacobian, , everywhere in ; however, in this application, is used only in the construction of . Therefore, we only require that is defined only over the triangles that contain the nodes in . We shall call this region . An illustration of is presented in Fig. 4.

After constructing over , we we may construct , and subsequently . While there is no issue approximating by using a quadrature rule, in this work, we choose to interpolate the values of at the nodes of and then use a quadrature rule to compute the natural boundary condition term of the variational form.

##### A Complete Description of the Algorithm

In the preceding paragraphs, we have given the mathematical statement of the modified Dirichlet–Neumann iterative approach and provided some key details in the implementation of this algorithm. We shall now provide a sketch of what a code implementation of (13) would look like.

###### Algorithm 1 (Modified Dirichlet–Neumann).

Define and , and and as the bases of and respectively. Given triangulations and , an error tolerance , and a relaxation parameter ,

1. Determine .

2. Obtain the element patch information for the nodes in for the superconvergent gradient algorithm.

3. Get the nearest–neighboring nodes of the interfaces and .

4. Set initial guess .

5. While :

1. Solve

 ahs(uhs,k,ψj)=⟨fs,ψj⟩Ωs,hfor j=1,…,M

with interface condition

 uhs,k∣∣Γcs,h=gk.
2. Compute over .

3. On :

1. Compute from formula (6).

2. Compute .

3. Compute .

4. Solve

 ahm(uhm,k,Ψj′)−⟨Im,h¯¯¯σ(uhs,k)nm,Ψj′⟩Γcm,h=⟨fm,Ψj′⟩Ωm,hfor j′=1,…,M′.
5. Compute

 gk+1=ωIs,h(uhm,k(^xm)−∇uhs,k⋅(^xm−xs))+(1−ω)gk.

## 4 Numerical Results

In this section, we will demonstrate that the iterative coupling procedure presented in the previous section allows us to achieve the error bounds expected from the piecewise linear finite element method.

For simplicity, let us define as the solution of the continuous interface coupling problem (2), and as the converged solution of the iterative formulation (13). In addition, let . We will report the error in the broken norms defined as

 ∥v∥L2(Ωh)=(∥vs∥2L2(Ωs,h)+∥vm∥2L2(Ωm,h))12

and

 ∥v∥H1(Ωh)=(∥vs∥2H1(Ωs,h)+∥vm∥2H1(Ωm,h))12

for all functions in . In addition, we will use the notation , to denote the number of elements on . Also, again, we define as the maximum element diameter in the of both triangulations and . Finally, let be the number of iterations required to reach convergence.

The convergence results in the following numerical experiments will be structured in the following manner. First, we shall discuss the domain decomposition problem, i.e., where the Lamé parameters of the subproblems are equal. In this test problem, we will first demonstrate that the linear patch test is passed for subdomains with noncoincident curved interfaces, and then we will demonstrate that we get the expected finite element convergence rates for piecewise linear approximation spaces. Second, we will present the convergence rates for the general interface coupling problem where the Lamé parameters of the solids are not necessarily equal. The error tolerance was used in all of the following numerical experiments. For convenience, we call the ratio between the number of elements on the slave interface and the number of elements on the master interface the “slave : master ratio”.

### 4.1 The Domain Decomposition Problem

For the sake of convenience, we set the Lamé parameters on both subdomains. The coupled problem domain in this test problem is a rectangle that occupies the region with discrete interfaces generated from the analytic representation , where at . is taken to be the triangulation of the subdomain on the left side of the continuous interface, and is taken to be the triangulation of the subdomain on the right side of the continuous interface.

#### 4.1.1 The Linear Patch Test

In this numerical experiment, we determine whether (13) can recover the linear vector function as an exact solution for the homogeneous problem. The results for some representative subdomain configurations are presented in Table 1. It is demonstrated that (13) for both spatially coincident and noncoincident interfaces. The relaxation parameter in the algorithm was set to . Figure 6: Left: The x–displacement solution of the linear patch test. Right: sin(πx)sin(2πy) as computed from the x–displacement of the solution of the iterative method.

#### 4.1.2 Convergence to a Manufactured Solution

Here, we determine the rate of convergence of the solution of the discrete problem to the manufactured solution . First, we present the convergence results for the case where , , coincide, and then we shall present the convergence history for the spatially noncoincident cases, where the ratio between the number of elements on the master interface and the slave interface vary. We have demonstrated in the Tables 2–4 that the iterative method presented in the previous section converges optimally for piecewise linear Lagrange elements. An illustration of the convergence rate can be found in Fig. 7.