Deformation problem for glued elastic bodies and an alternative iteration method

TL;DR

This paper models the deformation of glued elastic bodies using linear elasticity with adhesive forces, proving solution existence, analyzing an iterative solution method, and numerically validating convergence.

Contribution

It introduces a variational framework for the deformation problem and demonstrates the convergence of an alternating energy minimization method.

Findings

Unique existence of weak solutions established.

Alternating iteration method shown to be an energy minimization process.

Numerical results confirm convergence of the proposed methods.

Abstract

We study a mathematical model for deformation of glued elastic bodies in 2D or 3D, which is a linear elasticity system with adhesive force on the glued surface. We reveal a variational structure of the model and prove the unique existence of a weak solution based on it. Furthermore, we also consider an alternating iteration method and show that it is nothing but an alternating minimizing method of the total energy. The convergence of a monolithic formulation and the alternating iteration method are numerically studied with the finite element method.

Tables

Table 1. Table 1: Elapsed time of monolithic and Gauss-Seidel methods by FreeFem++ and Dissection solver with 4 cores. CPU time for factorization are also shown within parentheses.

mesh	$h_{\max }$	solver	# unknowns	# iteration	time (sec.)
					factorization		iteration	total
$20\times 20\times 20$	0.36551	monolithic	30,870	—	1.364	( 3.147 )	—	1.364
		Gauss-Seidel	14,742+16,128	12	1.339	( 2.752 )	0.245	1.585
$30\times 30\times 30$	0.24729	monolithic	97,743	—	6.901	( 18.588 )	—	6.901
		Gauss-Seidel	46,965+50,778	15	5.082	( 12.425 )	1.286	6.368
$60\times 60\times 60$	0.12783	monolithic	748,470	—	174.123	( 603.567 )	—	174.123
		Gauss-Seidel	368,928+379,542	20	118.973	( 316.482 )	25.559	144.532

Equations

$$\Omega \setminus \Gamma ={\Omega }_1\cup {\Omega }_2,{\Omega }_1\cap {\Omega }_2=\varnothing ,$$

$$\partial \Omega ={\partial }_D\Omega \cup \overline{{\partial }_N\Omega },{\partial }_D\Omega \cap \overline{{\partial }_N\Omega }=\varnothing ,{\mathcal{H}}^{d-1}(\overline{{\partial }_N\Omega }\setminus {\partial }_N\Omega )=0,$$

$$\left\{\begin{array}[]{ll}-\text{div}\,\sigma(u)=f(x)&(x\in\Omega\setminus\Gamma),\\ u=g(x)&(x\in\partial_{D}\Omega),\\ \sigma(u)\nu=q(x)&(x\in\partial_{N}\Omega),\\ \sigma(u_{1})\nu=\zeta(x)[u]=\sigma(u_{2})\nu&(x\in\Gamma).\end{array}\right.$$

$$e(u):=\frac{1}{2}\left(\nabla u^{\mathrm{T}}+{\left(\nabla u^{\mathrm{T}}\right)}^{\mathrm{T}}\right)\in {\mathbb{R}}_{\mathrm{sym}}^{d\times d}.$$

$$\sigma (u):=\mathcal{C}e(u)={\left(c_{klmn}{e}_{mn}(u)\right)}_{k,l}\in {\mathbb{R}}_{\mathrm{sym}}^{d\times d},$$

$$c_{klmn}(x){\xi }_{kl}{\xi }_{mn}\ge c_{\ast }|\xi {|}^2(\text{a.e.}x\in \Omega ,\xi \in {\mathbb{R}}_{\mathrm{sym}}^{d\times d}).$$

$$X_i:=H^1({\Omega }_i;{\mathbb{R}}^d),V_i:=\{u_i\in X_i;u_i{|}_{{\partial }_D{\Omega }_i}=0\}(i=1,2),$$

$$X:=H^1(\Omega \setminus \Gamma ;{\mathbb{R}}^d)\cong X_1\times X_2,V:=\{u\in X;u{|}_{{\partial }_D\Omega }=0\}\cong V_1\times V_2.$$

$$V(g):=V+g=\{u\in X;u-g\in V\},V_i(g_i):=V_i+g_i(i=1,2).$$

$$u\in V(g),a_0(u,v)=l_0(v)\text{mbox}forallv\in V,$$

$$a_0(u,v):={\int }_{\Omega \setminus \Gamma }\sigma (u):e(v)dx+{\int }_{\Gamma }\zeta [u]\cdot [v]ds(u,v\in X),$$

$$l_0(v):={\int }_{\Omega \setminus \Gamma }f\cdot vdx+{\int }_{{\partial }_N\Omega }q\cdot vds,(v\in X).$$

$$a_0(u,v)=\sum _{i=1}^2{a}_i(u_i,v_i)+a_{\Gamma }([u],[v])(u,v\in X),$$

$$a_i(u_i,v_i):={\int }_{{\Omega }_i}\sigma (u_i):e(v_i)dx,(i=1,2),a_{\Gamma }(u,v):={\int }_{\Gamma }\zeta u\cdot vds.$$

$$l_0(v)=\sum _{i=1}^2{l}_i(v_i),l_i(v_i):={\int }_{{\Omega }_i}f\cdot v_{i}dx+{\int }_{{\partial }_N{\Omega }_i}q\cdot v_{i}ds,(v_i\in X_i,i=1,2)$$

$$\tilde{u}\in V,a_0(\tilde{u},v)=l_0(v)+a_0(g,v)(v\in V).$$

$$u^{\ast }=\underset{u\in V(g)}{\arg \min }E(u),$$

$$E(u):=\frac{1}{2}{a}_0(u,u)-l_0(u).$$

$$a^{\ast }:=\underset{v,w\in V}{\sup }\frac{a_0(v,w)}{\Vert v{\Vert }_X\Vert w{\Vert }_X}<\infty .$$

$$u_h\in V_h(g):=V_h+g,a_0(u_h,v_h)=l_0(v_h)(v_h\in V_h).$$

$$\Vert u-u_h{\Vert }_X\le \frac{a^{\ast }}{a_{\ast }}\underset{v_h\in V_h(g)}{\inf }\Vert u-v_h{\Vert }_X.$$

$$a_1(u_1^{\ast },v_1)+a_{\Gamma }(u_1^{\ast },v_1)=a_{\Gamma }(u_2^{\ast },v_1)+l_1(v_1){(}^{\forall }{v}_1\in V_1),$$

$$a_2(u_2^{\ast },v_2)+a_{\Gamma }(u_2^{\ast },v_2)=a_{\Gamma }(u_1^{\ast },v_2)+l_2(v_2){(}^{\forall }{v}_2\in V_2).$$

$$a_1(u_1^m,v_1)+a_{\Gamma }(u_1^m,v_1)=a_{\Gamma }(u_2^{m-1},v_1)+l_1(v_1){(}^{\forall }{v}_1\in V_1,m=1,2,\cdots ),$$

$$a_2(u_2^m,v_2)+a_{\Gamma }(u_2^m,v_2)=a_{\Gamma }(u_1^m,v_2)+l_2(v_2){(}^{\forall }{v}_2\in V_2,m=1,2,\cdots ).$$

$$E((u_1^{m-1},u_2^{m-1}))\ge E((u_1^m,u_2^{m-1}))\ge E((u_1^m,u_2^m))(m=1,2,\cdots ).$$

$$E(u^{m,m-1})-E(u^{m,m})$$

$$=$$

$$\ge$$

$$=$$

$$[A_{\mu \nu }^{(i)}{]}_{kl}$$

$$[M^{(ij)}{]}_{kl}$$

$$\left[\begin{array}{cccc}{A}_{II}^{(1)}& A_{IB}^{(1)}& & \\ A_{BI}^{(1)}& \text{mbox}{A}_{BB}^{(1)}+M^{(11)}\text{mbox}& -M^{(12)}& \\ & -M^{(21)}& A_{BB}^{(2)}+M^{(22)}\text{mbox}& A_{BI}^{(2)}\\ & & A_{IB}^{(2)}& A_{II}^{(2)}\end{array}\right]\left[\begin{array}{c}{u}_{1,I}\\ u_{1,B}\\ u_{2,B}\\ u_{2,I}\end{array}\right]=\left[\begin{array}{c}{f}_{1,I}\\ f_{1,B}\\ f_{2,B}\\ f_{2,I}\end{array}\right],$$

$$M=M^{(11)}=M^{(12)}=M^{(21)}=M^{(22)}.$$

Full text

1 Introduction

We consider a mathematical model which describes deformation of two elastic bodies glued to each other on a surface. The understanding of such glued structure or adhesive bonding process is important in industrial and scientific applications, especially in the case that the mechanical bonding technique exhibits its disadvantages comparing with the adhesive one, e.g. bonding between soft materials, or one between very small scaled materials. The importance of mathematical modeling and numerical simulation is increasing in the design of desirable mechanical properties of composite materials with glued layer structure.

In mathematics, M. Frémond Fremond2002 and T. Roubíček et al. RSZ2009 proposed mathematical models of such glued structure and its delamination process. R. Scala Scala2017 studied more extended delamination model with kinetic and viscoelastic terms and proved existence of a solution. For further mathematical studies on the delamination process, we refer the above works and references therein.

In this paper, we concentrate on the stationary deformation problem of the glued structure, which is a linear elasticity system with adhesive force on the glued surface. A similar stationary problem also appears in the implicit time discretization of the above delamination models Yoneda2018 .

The outline of this paper is as follows. In Section 2, we describe the deformation model and give a definition of a weak solution. Section 3 is devoted to review several known consequences from the coercivity of a bilinear form; the existence and the uniqueness of a weak solution, a variational principle, and an error estimate of a finite element approximation.

For the purpose of efficient numerical computation of the obtained weak form of our model in 2D and 3D, we propose an alternating iteration method in Section 4. The alternating iteration method was proposed in Yoneda2018 and was used to simulate the vibration-delamination model proposed in Scala2017 . We will show that it is nothing but an alternating energy minimization procedure. In particular, it generates a sequence of displacements which monotonically decreases the total energy (Theorem 4.1).

In Section 5, we consider finite element discretization. We give discrete forms of the monolithic method and the alternating iteration method, and prove that the sequence generated by the alternating iteration method converges to the discrete solution by the monolithic method. Those theoretical results is also supported by numerical experiments in three dimensional setting.

2 Deformation of glued elastic bodies

We consider a bounded Lipschitz domain $\Omega\subset{\mathbb{R}}^{d},~{}(d=2,3)$. We suppose that

[TABLE]

where $\Gamma$ is a Lipschitz surface which is the common boundary of two disjoint Lipschitz domains $\Omega_{1}$ and $\Omega_{2}$ as shown in Fig. 2 and Fig. 2. We denote by $\nu$ the unit normal vector on $\Gamma$ pointing from $\Omega_{1}$ into $\Omega_{2}$, and the one on $\partial\Omega$ pointing outward. We suppose that the boundary $\partial\Omega$ is decomposed to the following disjoint portions:

[TABLE]

where $\partial_{D}\Omega$ and $\partial_{N}\Omega$ are relatively open subsets of $\partial\Omega$ and ${\mathcal{H}}^{d-1}$ denotes the $d-1$ dimensional Hausdorff measure. We also define $\partial_{D}\Omega_{i}:=\partial_{D}\Omega\cap\partial\Omega_{i}$ and $\partial_{N}\Omega_{i}:=\partial_{N}\Omega\cap\partial\Omega_{i}$ for $i=1,2$. We assume that $\partial_{D}\Omega$ is a nonempty relatively open subset of $\partial\Omega$, and ${\mathcal{H}}^{d-1}(\partial_{D}\Omega_{1})$ or ${\mathcal{H}}^{d-1}(\partial_{D}\Omega_{2})$ is positive. Without loss of generality, we always suppose ${\mathcal{H}}^{d-1}(\partial_{D}\Omega_{1})>0$ (Fig. 2, Fig. 2).

In this paper, we consider the following stationary deformation model of two elastic bodies $\Omega_{1}$ and $\Omega_{2}$ which are glued by an adhesive on $\Gamma$. We consider an adhesion force but ignore a friction force on the interface.

The problem is to find $u:\Omega\setminus\Gamma\to{\mathbb{R}}^{d}$ such that:

[TABLE]

The meanings of the above symbols are as follows. We use the Einstein summation convention in this section. For matrices $\xi=(\xi_{kl}),~{}\eta=(\eta_{kl})\in{\mathbb{R}}^{d\times d}$, we denote their component-wise inner product by $\xi:\eta:=\xi_{kl}\eta_{kl}$ and the norm by $|\xi|:=\sqrt{\xi:\xi}$.

The solution $u$ is a displacement field on $\Omega\setminus\Gamma=\Omega_{1}\cup\Omega_{2}$. We denote $u|_{\Omega_{i}}$ by $u_{i}$ for $i=1,2$, and often write $u=(u_{1},u_{2})$. The symmetric gradient of $u$ is defined by

[TABLE]

The stress tensor $\sigma(u)\in{\mathbb{R}}^{d\times d}_{\mathrm{sym}}$ satisfies the constitutive relation

[TABLE]

where ${\mathcal{C}}:={\mathcal{C}}(x)=(c_{klmn}(x))\in{\mathbb{R}}^{d\times d\times d\times d}$ is the elasticity tensor with the symmetries $c_{klmn}=c_{mnkl}=c_{lkmn},(1\leq k,l,m,n\leq d)$.

We assume that ${\mathcal{C}}\in L^{\infty}(\Omega;{\mathbb{R}}^{d\times d\times d\times d})$ and there exists $c_{*}>0$ such that

[TABLE]

The first equation of (1) is the force balance equation in each subdomain $\Omega_{i}$, where $f$ is a given body force. The second and third equations are Dirichlet and Neumann boundary conditions, where $g$ is a given displacement on $\partial_{D}\Omega$ and $q$ is a given surface traction on $\partial_{N}\Omega$.

In the fourth equation, $\zeta(x)\geq 0$ is a given adhesive parameter on the glued surface $\Gamma$ which represents the strength of adhesive bonding. The adhesive force at $x\in\Gamma$ is assumed to be $\zeta(x)[u]$, where $[u(x)]:=(u_{2}(x)-u_{1}(x))$ is the gap of the displacement $u=(u_{1},u_{2})$ on $\Gamma$. It should be balanced with the surface traction force $\sigma(u_{1}(x))\nu$ on $\Omega_{1}$ side and also with $\sigma(u_{2}(x))\nu$ on $\Omega_{2}$ side.

To consider a weak formulation of (1), we introduce the following spaces.

[TABLE]

For $g=(g_{1},g_{2})\in X$, we also define affine spaces:

[TABLE]

Definition 2.1** (Weak solution)**

We suppose that a Dirichlet boundary data $g\in X$, a body force $f\in L^{2}(\Omega;{\mathbb{R}}^{d})$, a surface traction $q\in L^{2}(\partial_{N}\Omega;{\mathbb{R}}^{d})$, and an adhesive coefficient $\zeta\in L^{\infty}(\Gamma)$, $\zeta(x)\geq 0$ are given. Then, we call $u$ a weak solution of (1) if

[TABLE]

where

[TABLE]

The bilinear form $a_{0}$ and the linear form $l_{0}$ are decomposed into sum of the following subforms:

[TABLE]

We remark that, if $c_{klmn}\in C^{1}(\overline{\Omega_{i}})$ for $i=1,2$, then a strong solution of (1), i.e., $u\in H^{2}(\Omega\setminus\Gamma;{\mathbb{R}}^{d})$ which satisfies (1) almost everywhere on $\Omega\setminus\Gamma$ or on $\partial\Omega\cup\Gamma$, is a weak solution. On the other hand, if a weak solution belongs to $H^{2}(\Omega\setminus\Gamma;{\mathbb{R}}^{d})$, then it is a strong solution.

3 Unique existence of a weak solution

For establishing the unique existence of a weak solution, the coercivity of the bilinear form $a_{0}(u,v)$ defined in (2) is essential.

Theorem 3.1** (Coercivity of $a_{0}$)**

We suppose that $\zeta(x)\geq 0$ and $\|\zeta\|_{L^{\infty}(\Gamma)}>0$. Then there exists $a_{*}>0$ such that $a_{0}(v,v)\geq a_{*}\|v\|_{X}^{2}$ holds for all $v\in V$.

A slightly long proof of this theorem using an argument by contradiction was given in Yoneda2018 and another simpler proof will be given in our forthcoming paper. We postpone the proof to it and here we just remark that the case of ${\mathcal{H}}^{d-1}(\partial_{D}\Omega_{i})>0$ for both $i=1,2$ is relatively easily shown by using Körn’s second inequality D-L1976 .

As a consequence of Theorem 3.1, we immediately have the unique existence of a weak solution.

Theorem 3.2** (Unique existence)**

We suppose that Dirichlet boundary data $g\in X$, a body force $f\in L^{2}(\Omega;{\mathbb{R}}^{d})$, a surface traction $q\in L^{2}(\partial_{N}\Omega;{\mathbb{R}}^{d})$ are given. Then, under the conditions of Theorem 3.1, there exists a unique weak solution to (1).

Proof

∎Under the conditions, from the definitions of $a_{0}$ and $l_{0}$, we can show that $a_{0}$ is a continuous symmetric bilinear form on $X\times X$ and $l_{0}$ is a continuous linear form on $X$.

We set $\tilde{u}:=u-g$. Then $u$ is a weak solution to (1) if and only if

[TABLE]

From the Lax-Milgram lemma Ciarlet1978 and Theorem 3.1, there exists a unique $\tilde{u}$ which satisfies (3). Hence the unique existence of the weak solution has been proved. ∎

A variational principle for the above symmetric Lax-Milgram type problem is also well known. The weak solution $u^{*}\in V(g)$ is a unique minimizer of the following energy:

[TABLE]

where

[TABLE]

In Section 5, we consider a finite element approximation for our model (1). So called Céa’s lemma Ciarlet1978 implies the following error estimate. We define

[TABLE]

Proposition 3.3

Under the assumptions of Theorem 3.2, we suppose that $V_{h}$ is a closed subspace of $V$. Then there uniquely exists $u_{h}$ such that

[TABLE]

Furthermore, it satisfies

[TABLE]

The problem (6) corresponds to the finite element scheme. If $V_{h}$ is a space of piece-wise linear element (P1 element) on a regular triangular mesh, it is known that $\inf_{v_{h}\in V_{h}(g)}\|u-v_{h}\|_{X}=O(h)$ as the mesh size $h$ tends to [math] under suitable regularity for $u$ and the triangular mesh Ciarlet1978 .

4 Alternating iteration method

We remark that $u^{*}=(u_{1}^{*},u_{2}^{*})\in V(g)$ is a weak solution to (1) if and only if

[TABLE]

We consider the following alternating method.

**Gauss-Seidel type scheme

**For given $u_{2}^{0}\in V_{2}(g)$, and for $m=0,1,2,\cdots$, find $u^{m}=(u_{1}^{m},u_{2}^{m})\in V(g)$ such that

[TABLE]

We call the above alternating iteration method “Gauss-Seidel type” by analogy with an iterative solver for linear systems. The unique solvability of each $u_{i}^{m}$ is clear from Körn’s second inequality.

The following theorem tells us that the Gauss-Seidel type scheme defines a sequence $\{u^{m}\}$ which monotonically decreases the total energy $E(u)$ defined in (5).

Theorem 4.1

The obtained sequence $\{u^{m}=(u_{1}^{m},u_{2}^{m})\}_{m}$ by the Gauss-Seidel type scheme satisfies the following energy decay property:

[TABLE]

Proof

For simplicity, we denote $(u_{1}^{m},u_{2}^{n})\in V(g)$ by $u^{m,n}$. The second inequality is shown as follows. Setting $v=(v_{1},v_{2}):=u^{m,m-1}-u^{m,m}=(0,u_{2}^{m-1}-u_{2}^{m})$, we have

[TABLE]

where we have used $v_{1}=0$ and (10) for the last equality. The first inequality in (11) is shown in the same way. ∎

Since the weak solution $u^{*}$ is the minimizer of the total energy $E$ as written in (4), the sequence generated by the Gauss-Seidel type scheme is expected to approximate $u^{*}$. We will study it numerically in the next section.

5 Numerical solution

First we recall the assumption ${\mathcal{H}}^{d-1}(\partial_{D}\Omega)>0$. In this section we only deal with the case $\partial_{D}\Omega\subset\partial\Omega_{1}$ and $\zeta(x)>0$ on $\Gamma$.

5.1 A matrix representation of the monolithic formulation

We briefly review a matrix representation of the monolithic formulation. Let $\Lambda^{(i)}$ be an index for finite element basis function and be decomposed as $\Lambda^{(i)}_{I}\cup\Lambda^{(i)}_{B}$, corresponding to nodes in the subdomain $\overline{\Omega_{i}}\setminus\Gamma$ and on the common boundary $\Gamma$.

We define the following stiffness matrices in each subdomain $\Omega_{i}$ using the bilinear forms $\{a_{i}(\cdot,\cdot)\}_{i=1,2}$ defined in Section 2.

[TABLE]

Here combination of four mass matrices $\{M^{(i\,j)}\}$ provides a matrix representation of the bilinear form $a_{\Gamma}([\cdot],[\cdot])$. A matrix representation of the monolithic formulation reads

[TABLE]

where the right hand side consists of the body force $f$ and inhomogenous Dirichlet and Neumann data $g$ and $q$.

Remark 5.1

The monolithic method can be computed by $LDU$-factorization with any symmetric permutation because the stiffness matrix is positive definite thanks to the coercivity of the weak form with ${\mathcal{H}}^{d-1}(\partial_{D}\Omega)>0$ (Theorem 3.1, Proposition 3.3).

5.2 Alternating iterative

method in discrete form

In the following, we suppose that the nodal points and the surface meshes of the mesh decomposition of domains $\Omega_{1}$ and $\Omega_{2}$ coincide on the interface $\Gamma$, which leads to

[TABLE]

Let us define an inner product of vector $u_{i,\,\mu}$ and $v_{i,\,\mu}$ as $(u_{i,\,\mu},v_{i,\,\mu}):=\sum_{k\in\Lambda_{\mu}^{(i)}}(u_{\mu})_{k}(v_{\mu})_{k}$ for $i=1,2$ and $\mu\in\{I,B\}$, and denote the standard $\ell^{2}$-norm by $\|u_{i,\,\mu}\|=(u_{i,\,\mu},u_{i,\,\mu})^{1/2}$. Since $\zeta(x)>0$ on $\Gamma$, the mass matrix $M$ is positive definite and $(Mu_{i,B},v_{i,B})$ becomes an inner product with weight $M$. Hence, we denote a norm with the weight $M$ by $\|u_{i,B}\|_{M}=(M\,u_{i,B},u_{i,B})^{\frac{1}{2}}$. Then there exist $\beta_{1}>0$ and $\beta_{2}>0$ such that $\beta_{1}\|u_{i,B}\|_{M}^{2}\leq\|u_{i,B}\|^{2}\leq\beta_{2}\|u_{i,B}\|_{M}^{2}$ holds for any $u_{i,B}$.

We prepare two matrices in each subdomain to describe the linear system in a simpler way,

[TABLE]

Lemma 5.2

There exists $\alpha_{1}>0$ such that

[TABLE]

holds for any vector $[v_{1,I}\ \ v_{1,B}]^{T}$.

Proof

∎Since ${\cal A}^{(1)}$ in $\Omega_{1}$ is positive definite due to the Dirichlet boundary $\partial_{D}\Omega\subset\partial\Omega_{1}$, there exists $\alpha_{1}>0$ and the first inequality holds. The second one is clear from the definition of $\beta_{1}$. ∎

The Gauss-Seidel type iteration defined in $(\ref{uGS1})$ and $(\ref{uGS2})$ is written in the following matrix presentation.

Algorithm 5.3** (Gauss-Seidel type iteration)**

Let $u_{2,B}^{0}$ be an initial guess of $u_{2}$ on $\Gamma$. From obtained data $u_{2,B}^{m-1}$ of $m-1$-step, approximate solution $[u_{1,I}^{m}\ u_{1,B}^{m}]^{T}$ and $[u_{2,I}^{m}\ u_{2,B}^{m}]^{T}$ are generated by solving following two problems successively,

[TABLE]

Let $[u_{1,I}^{m}\ u_{1,B}^{m}]^{T}$ and $[u_{2,I}^{m}\ u_{2,B}^{m}]^{T}$ ($m=1,2,\cdots$) be the solution of the Gauss-Seidel type iteration and let $[u_{1,I}^{*}\ u_{1,B}^{*}]^{T}$ and $[u_{2,I}^{*}\ u_{2,B}^{*}]^{T}$ be the one of the monolithic system (12). We define the error between them by

[TABLE]

We have the following convergence estimates for the Gauss-Seidel type iteration.

Lemma 5.4

The error on the boundary admits the following estimates:

[TABLE]

where $r:=1/(1+\alpha_{1}\beta_{1})<1$.

Proof

From the definition of the errors (13), they satisfy the following linear systems:

[TABLE]

Taking inner product of $[e_{1,I}^{m}\ \ e_{1,B}^{m}]^{T}$ with the left equation of (14), and of $[e_{2,I}^{m}\ \ e_{2,B}^{m}]^{T}$ with the right, we obtain

[TABLE]

From Lemma 5.2 and (15), we have

[TABLE]

This gives the first inequality. The second inequality is obtained from (16) with positive semi-definiteness of ${\cal A}^{(2)}$. ∎∎

Theorem 5.5

There exists $C>0$ such that the following error estimate holds:

[TABLE]

where $r:=1/(1+\alpha_{1}\beta_{1})<1$.

Proof

From Lemma 5.4, we have

[TABLE]

These inequalities imply

[TABLE]

Since the matrices $\widetilde{\cal A}^{(i)}$ $(i=1,2)$ are invertible, from (14), we obtain

[TABLE]

Together with the estimate (17) and with the fact that $M$ is positive definite, there exists a $C>0$ such that the assertion of the theorem holds. ∎∎

Remark 5.6

The Gauss-Seidel type iteration is straightforwardly extended to SOR type iteration by introducing a relaxation parameter.

5.3 Numerical results

Now we show numerical verification on convergence of Gauss-Seidel type iteration using a manufactured solution,

[TABLE]

in $\Omega_{1}=(0,4)\times(0,2)\times(0,4)$, $\Omega_{2}=(0,4)\times(2,4)\times(0,4)$, and $\partial_{D}\Omega=\{(x,y,z)\,;\,0<x<4,y=0,0<z<4\}$ with corresponding inhomogeneous Dirichlet data $g(x)$, the load $f(x)$ on $\Omega$, and the surface traction $q(x)$ on $\partial_{N}\Omega$. For simplicity, we put $\zeta(x)\equiv 1$. Figure 3 shows relative errors of finite element solution discretized with P1 element solved by the monolithic formulation, which ensures the 1st order approximation error of the solution, $O(h)$ with mesh size $h$.

Convergence of Gauss-Seidel type iteration to the solution by the monolithic formulation with the relative error measured by $\|\cdot\|_{H^{1}(\Omega)}$ is shown in the left of Figure 4, and relative error to the manufactured solution in the right of Figure 4. Here mesh subdivisions with $h_{\max}=0.36551$ in $20\times 20\times 20$, and $h_{\max}=0.12878$ in $60\times 60\times 60$ are used. We can see the convergence does not depend on the mesh size, and the same relative error as one by monolithic formulation to the manufactured solution is obtained after certain iterations, though Gauss-Seidel type iteration continues to converge.

5.4 Computational efficiency

We used FreeFem++ software package and Dissection sparse direct solver SuzukiRoux2014 to obtain finite element solution. Table 1 shows computational time of the direct solver for the monolithic formulation and the Gauss-Seidel type iteration, using Intel Core processor i7-6770HQ with 4 cores running at 2.60GHz. $LDU$-factorization is performed before starting iteration and forward/backward substitutions are repeated because the stiffness matrix in each subdomain does not change during the iteration. Since computational complexity of $LDU$-factorization of sparse matrix with P1 finite element is more than $O(N^{2})$ with number of unknowns $N$, factorization cost for sub-matrices in Gauss-Seidel type is less than half, $2\times(N/2)^{2}/N^{2}=1/2$ of the one for monolithic formulation. We can observe shorter CPU time for factorization of Gauss-Seidel type, which is shown as number in parentheses, when the problem size is enough large. In Dissection solver, numerical factorization is fully parallelized but there exist some sequential processing parts, e.g. fill-in analysis of the sparse matrix, which masks speed-up of elapsed time for the factorization in Gauss-Seidel type. We also observe that elapsed time of iteration in Gauss-Seidel type with selected iterating number to obtain appropriately approximate solution is less than time of the factorization. Hence if we can find a reasonable criteria to stop iteration, the Gauss-Seidel type iteration becomes more efficient than monolithic formulation.

6 Conclusion

We considered a stationary deformation problem for glued elastic bodies and have established solvability of its weak formulation. We proposed a kind of alternating iteration scheme to approximate the problem and showed that the scheme has a nature of alternating minimizing algorithm with respect to the total energy.

We proved the convergence for the Gauss-Seidel type iteration in a rate of $O(r^{m})$ with $r\in(0,1)$ in discrete setting. Computational efficiency of the monolithic and alternating iterative algorithms have been verified with a three dimensional problem. The alternating iterative method requires smaller computational resource than the monolithic method and also it has an advantage in computation time when the degree of freedom is sufficiently large.

Acknowledgements.

The authors are grateful to Prof. Frederic Hecht for his useful comments on numerical computation with FreeFem++. This work is supported by JSPS KAKENHI Grant Number JP16H03946 and JP17H02857.