Coupling regularization and adaptive hp-BEM for the solution of a delamination problem

Coupling regularization and adaptive hp-BEM for the solution of a delamination problem

Nina Ovcharova N. Ovcharova Universität der Bundeswehr München, D-85577 Neubiberg/Munich, Germany
1L. Banz Department of Mathematics, University of Salzburg, Hellbrunner Straße 34, 5020 Salzburg, Austria
2
Lothar Banz N. Ovcharova Universität der Bundeswehr München, D-85577 Neubiberg/Munich, Germany
1L. Banz Department of Mathematics, University of Salzburg, Hellbrunner Straße 34, 5020 Salzburg, Austria
2
2email: nina.ovcharova@unibw.de
4email: lothar.banz@sbg.ac.at
Abstract

In this paper, we couple regularization techniques with the adaptive -version of the boundary element method (-BEM) for the efficient numerical solution of linear elastic problems with nonmonotone contact boundary conditions. As a model example we treat the delamination of composite structures with a contaminated interface layer. This problem has a weak formulation in terms of a nonsmooth variational inequality. The resulting hemivariational inequality (HVI) is first regularized and then, discretized by an adaptive -BEM. We give conditions for the uniqueness of the solution and provide an a-priori error estimate. Furthermore, we derive an a-posteriori error estimate for the nonsmooth variational problem based on a novel regularized mixed formulation, thus enabling -adaptivity. Various numerical experiments illustrate the behavior, strengths and weaknesses of the proposed high-order approximation scheme.

Keywords:
Hemivariational inequality Regularization techniques Delamination problem a-priori and a-posteriori error estimates -BEM
65N38 74M15

1 Introduction

The interest in robust and efficient numerical methods for the solution of nonsmooth problems in nonmonoton contact like adhesion and delamination problems increases constantly. The nonsmoothness comes from the nonsmooth data of the problems itself, in particular from nonmonotone multivalued physical laws involved in the boundary conditions that lead to nonsmooth functionals in the variational formulation. There are several approaches to treat this non-differentiability. We can first discretize by finite elements, and then solve the nonconvex optimization problems by novel nonsmooth optimization methods Noll_Ovcharova , or first regularize the nonsmooth functional, then discretize by finite element methods and finally use standard optimization solvers Ovcharova2012 . Since the interesting nonmonotone behavior takes place only in a boundary part of the involved linear elastic bodies, the domain hemivariational inequaliy can also be formulated as a boundary hemivariational inequality via the Poincaré-Steklov operator. Therefore, the boundary element methods (BEM) are the methods of choice for discretization. There are number of works exploiting the - and the high-order -version of the BEM for monotone contact problems. Advanced -BEM discretizations, a-priori and a-posteriori error estimates for unilateral Signorini problems or contact problems with monotone friction are established in banz2013posteriori ; Banz2015 ; chernov2008hp ; MS-1 ; maischak2007adaptive ; stephan2009hp . However, successfully applied -BEM techniques for nonsmooth nonconvex variational inequalities are still missing. Pure -BEM with lowest polynomial degree has been investigated for the first time by Nesemann and Stephan Nesemann for unilateral contact with adhesion. They study the existence and uniqueness, and introduce a residual error indicator.

In this paper, we extend the approximation scheme from Ovcharova2015 (dealing with the - version of the BEM) to the -version of the BEM for the regularized problem. More precisely, the nonsmooth variational inequality is first regularized and then discretized by -adaptive BEM. For the Galerkin solution of the regularized problem we provide an a-priori error estimate and derive a reliable a-posteriori error estimate based on a novel regularized mixed formulation and thus, enabling -adaptivity. All our theoretical results are illustrated with various numerical experiments for a contact problem with adhesion. We also provide conditions for the uniqueness of the solution that sharpen the results due to Nesemann .

Notation: We denote with , and such alike generic constants, which can take different values at different positions. We use bold font to indicate vector-valued variables, e.g. . If there are too many indices, one of them is written as a superscript, e.g. .

2 A nonmonotone boundary value problem from delamination

Let be a bounded domain with Lipschitz boundary . We assume that the boundary is decomposed into three disjoint open parts , and such that and, moreover, the measures of and are positive. We fix an elastic body occupying . The body is subject to volume force and , , is a gap function associating every point with its distance to the rigid obstacle measured in the direction of the unit outer normal vector . The body is fixed along , surface tractions act on , and on the part a nonmonotone, generally multivalued boundary condition holds.

Further, denotes the linearized strain tensor and stands for the stress tensor, where is the Hooke tensor, assumed to be uniformly positive definite with coefficients. The boundary stress vector can be decomposed further into the normal, respectively, the tangential stress:

 σn=σ(u)n⋅n,σt=σ(u)n−σnn.

By assuming that the structure is symmetric and the forces applied to the upper and lower part are the same, it suffices to consider only the upper half of the specimen represented by , see Fig. 2 left for the 2D benchmark problem.

The delamination problem under consideration is the following: Find such that

 −divσ(u) =f inΩ (1a) u =0 onΓD (1b) σ(u)n =t onΓN (1c) un ≤g onΓC (1d) σt(u) =0 onΓC (1e) −σn(u) ∈∂f(un) onΓC. (1f)

The contact law (1f), written as a differential inclusion by means of the Clarke subdifferential Clarke of a locally Lipschitz function , describes the nonmonotone, multivalued behavior of the adhesive. More precisely, is the physical law between the normal component of the stress boundary vector and the normal component of the displacement on . A typical zig-zagged nonmonotone adhesion law is shown in Figure 1(b).
To give a variational formulation of the above boundary value problem we define

 H1D(Ω)={v∈H1(Ω):v|ΓD=0},K={v∈H1D(Ω):vn≤g a.e.~{}on ΓC}

and introduce the -coercive and continuous bilinear form of linear elasticity

 a(u,v)=∫Ωσ(u):ε(v)dx.

Multiplying the equilibrium equation (1a) by , integrating over and applying the divergence theorem yields

 ∫Ωσ(u):ε(v−u)dx=∫Ωf⋅(v−u)dx+∫Γσ(u)n⋅(v−u)ds.

From the definition of the Clarke subdifferential, the nonmonotone boundary condition (1f) is equivalent to

 −σn(u)(vn−un)≤f0(un;vn−un).

Here, the notation stands for the generalized directional derivative of at in direction . Substituting by on , using on the decomposition

 σ(u)n⋅(v−u)=σt(u)⋅(vt−ut)+σn(u)(vn−un)

and taking into account that on no tangential stresses are assumed, c.f. (1e), we obtain the hemivariational inequality (HVI): Find such that

 a(u,v−u)+∫ΓCf0(un(s);vn(s)−un(s))ds≥∫Ωf⋅(v−u)dx +∫ΓNt⋅(v−u)ds∀v∈K. (2)

3 Boundary integral operator formulation

Since the main difficulties of the boundary value problem, namely the nonmonotone adhesion law (1f) and the unilateral contact condition (1d) appear on the boundary, it might be viewed advantageous to formulate (2) as a boundary integral problem. To this end, we introduce the free boundary part and recall the Sobolev spaces HsWendl :

 H1/2(Γ)={v∈L2(Γ):∃v′∈H1(Ω),trv′=v},H1/2(Γ0)={v=v′|Γ0:∃v′∈H1/2(Γ)},~H1/2(Γ0)={v=v′|Γ0:∃v′∈H1/2(Γ),suppv′⊂Γ0}

with the standard norms

 ∥u∥H1/2(Γ0)=infv∈H1/2(Γ),v|Γ0=u∥v∥H1/2(Γ)and∥u∥~H1/2(Γ0)=∥u0∥H1/2(Γ),

where is the extension of onto by zero. The Sobolev space of negative order on are defined by duality as

 H−1/2(Γ0)=(~H1/2(Γ0))∗and~H−1/2(Γ0)=(H1/2(Γ0))∗.

Moreover, from (HsWendl, , Lemma 4.3.1) we have the inclusions

 ~H1/2(Γ0)⊂H1/2(Γ0)⊂L2(Γ0)⊂~H−1/2(Γ0)⊂H−1/2(Γ0).

For the solution of (1a) with we have the following representation formula, also known as Somigliana’s identity, see e.g. Kleiber :

 (3)

where is the fundamental solution of the Navier-Lamé equation defined by

 G(x,y)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩λ+3μ4πμ(λ+2μ)(log|x−y|I+λ+μλ+3μ(x−y)(x−y)⊤|x−y|2), if d=2λ+3μ8πμ(λ+2μ)(|x−y|−1I+λ+μλ+3μ(x−y)(x−y)⊤|x−y|3), if d=3

with the Lamé constants depending on the material parameters, i.e. the modulus of elasticity  and the Poisson’s ratio :

 λ=Eν1−ν2,μ=E1+ν.

Here, stands for the traction operator with respect to defined by , where is the unit outer normal vector at . Letting in (3), we obtain the well-known Calderón operator

 (uTxu)=(12I−KVW12I+K′)(uTxu)+(N0fN1f),

with the single layer potential , the double layer potential , its formal adjoint , and the hypersingular integral operator defined for as follows:

 (Vϕ)(x) :=∫ΓG(x,y)ϕ(y)dsy, (Kϕ)(x) :=∫ΓTyG⊤(x,y)ϕ(y)dsy (K′ϕ)(x) :=Tx∫ΓG(x,y)ϕ(y)dsy, (Wϕ)(x) :=−Tx(Kϕ)(x).

The Newton potentials are given for by

From Co-1 it is known that the linear operators

 V: H−1/2+σ(Γ)→H1/2+σ(Γ), K: H1/2+σ(Γ)→H1/2+σ(Γ) K′: H−1/2+σ(Γ)→H−1/2+σ(Γ), W: H1/2+σ(Γ)→H−1/2+σ(Γ)

are well-defined and continuous for . Moreover, is symmetric and positive definite (elliptic on ) in and, if the capacity of is smaller than 1, also in . That can always be achieved by scaling, since the capacity (or conformal radius or transfinite diameter) of is smaller than 1, if is contained in a disc with radius (see e.g. SlSp ; Steinbach ). is symmetric and positive semidefinite with kernel (elliptic on ). Hence, since is invertible, we obtain by taking the Schur complement of the Calderon projector that

 Txu=Pu−Nf (4)

with the symmetric Poincaré-Steklov operator and its Newton potential given by

 Pu=Wu+(K′+12I)V−1(K+12I)u, Nf=(K′+12I)V−1N0f−N1f.

Note that if , maps to its traction and, therefore, the Poincaré-Steklov operator is sometimes called the Dirichlet-to-Neumann mapping. Moreover, induces a symmetric bilinear form on , and is continuous and -elliptic, see e.g. cf , i.e. there exist constants , such that

 ∥Pv∥H−1/2(Γ)≤CP∥v∥H1/2(Γ)∀v∈H1/2(Γ),⟨Pv,v⟩Γ0≥cP∥v∥~H1/2(Γ0)∀v∈~H1/2(Γ0).

Here, is the -duality pairing between the involved spaces. Using the boundary functional sets

 V=~H1/2(Γ0)andKΓ={v∈V:vn≤g a.e.~{}on ΓC}

we obtain as in the domain based cases the boundary integral hemivariational inequality (BIHVI), Problem (: Find such that

 ⟨Pu,v−u⟩Γ0+∫ΓCf0(un(s);vn(s)−un(s))ds≥⟨t,v−u⟩ΓN +⟨Nf,v−u⟩Γ0∀v∈KΓ. (5)

To shorten the right hand side we introduce the linear functional

 ⟨F,v⟩Γ0=∫ΓNt⋅vds+⟨Nf,v⟩Γ0.

4 Regularization of the nonsmooth functional

In this section, we review from Ovcharova2012 ; Ov_Gw a class of smoothing approximations for nonsmooth functions that can be re-expressed by means of the plus function . The approximation is based on smoothing of the plus function via convolution.

Let be a small regularization parameter. The smoothing approximation of a locally Lipschitz function is defined via convolution by

 ^f(x,ε)=∫Rf(x−εt)ρ(t)dt,

where is a probability density function such that

 κ=∫R|t|ρ(t)dt<∞

and

We consider a class of nonsmooth functions that can be reformulated by means of the plus function. To this class of functions belong maximum, minimum or nested max-min function. If is a maximum function, for example, , then can be reformulated by using the plus function as

 f(x)=g1(x)+p(g2(x)−g1(x)), (6)

The smoothing function of is then defined by

 S(x,ε)=g1(x)+^p(g2(x)−g1(x),ε), (7)

where is the smoothing function of via convolution.

Choose, for example, the Zang probability density function

 ρ(t)={1if−12≤t≤120otherwise.

Then

 ^p(t,ε)=∫Rp(t−εs)ρ(s)ds=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩0ift<−ε212ε(t+ε2)2if−ε2≤t≤ε2tift>ε2

and hence,

 S(x,ε)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩g1(x)if g2(x)−g1(x)≤−ε212ε[g2(x)−g1(x)]2+12(g2(x)+g1(x))+ε8if |g2(x)−g1(x)|≤ε2g2(x)if g2(x)−g1(x)≥ε2. (8)

The relation (7) can be extended to the maximum function of continuous functions , i.e.

 f(x)=max{g1(x),g2(x),…,gm(x)} (9)

by

 S(x,ε)=g1(x)+^p(g2(x)−g1(x)+…+^p(gm(x)−gm−1(x),ε),ε). (10)

The major properties of the function in (10) are listed in the following lemma:

Lemma 1

(Ovcharova2012, , lemma 10) Let and be defined as in (9), (10), respectively, with . Then there holds:

(i) For any and for all ,

 |S(x,ε)−f(x)|≤(m−1)κε.

(ii) The function is continuously differentiable on and for any and there exist such that and

 ∂S(x,ε)∂x=Sx(x,ε)=m∑i=1Λig′i(x). (11)

Moreover,

 {limz→x,ε→0+Sx(z,ε)}⊆∂f(x). (12)
Remark 1

Since these type of functions can be handled as above.

Assumption 1

Assume that there exists positive constants such that for all

 |g′i(x)| ≤ci(1+|x|) (13a) g′i(x)⋅(−x) ≤di|x|. (13b)

By using (11) - (12) and (13a) - (13b), the following auxiliary result can be easily deduced.

Lemma 2

For any it holds that

 |Sx(x,ε)⋅z| ≤c(1+|x|)|z| (14a) Sx(x,ε)⋅(−x) ≤d|x| (14b) limsupz→x,ε→0+Sx(ε,x)ξ ≤f0(x;ξ). (14c)

Further, we introduce the functional defined by

 Jε(u)=∫ΓCS(un(s),ε)ds.

The regularized problem of () is given by: Find such that

 ⟨Puε,v−uε⟩Γ0+⟨DJε(uε),v−uε⟩ΓC≥⟨F,v−uε⟩Γ0∀v∈KΓ, (15)

where is the Gâteaux derivative of the functional and is given by

 ⟨DJε(u),v⟩ΓC=∫ΓCSx(un(s),ε)vn(s)ds.

To simplify the notations we introduce defined by

 φ(u,v)=∫ΓCf0(un(s);vn(s)−un(s))ds∀u,v∈H1/2(Γ). (16)

The regularized version of is therewith ,

The existence result for (5), resp. (15), is due to Gwinner_PhD and relies in both cases on the pseudomonotonicity of and , respectively. For details we refer the reader to Ovcharova2012 ; Ov_Gw .

Finally, we recall that the functional , where is a real reflexive Banach space, is pseudomonotone if (weakly) in and imply that, for all , we have

5 Uniqueness results

In this section, we discuss the uniqueness of solution of the boundary hemivariational inequality and of the corresponding regularization problem. The main results are presented in Theorem 5.2 which is based on the abstract uniqueness Theorem 5.1 from Ovcharova2015 that gives a sufficient condition for uniqueness. Similar uniqueness result but for the regularized problem is derived in Theorem 5.3. We are also dealing with the question which classes of smoothing functions preserve the property of unique solvability of the original problem as .

5.1 Uniqueness of the BIHVI

Assumption 2

Assume that there exists an such that for any it holds

 φ(u,v)+φ(v,u)≤α0∥u−v∥2V. (17)

To make the results of uniqueness self-consistent, we introduce from Ovcharova2015 the following theorem.

Theorem 5.1

(Ovcharova2015, , Theorem 5.1) Under the assumption (17) with , there exists a unique solution to the BIHVI problem , which depends Lipschitz continuously on the right hand side .

Proof  Assume that , are two solutions of . Then the inequalities below hold:

 ⟨Pu−F,v−u⟩Γ0+φ(u,v) ≥0∀v∈KΓ ⟨P~u−F,v−~u⟩Γ0+φ(~u,v) ≥0∀v∈KΓ.

Setting in the first inequality and in the second one, and summing up the resulting inequalities, we get

 ⟨Pu−P~u,~u−u⟩Γ0+φ(u,~u)+φ(~u,u)≥0. (18)

From the coercivity of the operator and the assumption (17) we obtain

 cP∥u−~u∥2V≤φ(u,~u)+φ(~u,u)≤α0∥u−~u∥2V.

Now let and denote Analogously to (18), we find that

 ⟨Pu1−F1−Pu2+F2,u2−u1⟩Γ0+φ(u1,u2)+φ(u2,u1)≥0.

Hence,

 cP∥u1−u2∥2V≤φ(u1,u2)+φ(u2,u1)+⟨F1−F2,u2−u1⟩Γ0

and by (17),

 (cP−α0)∥u1−u2∥2V≤⟨F1−F2,u2−u1⟩Γ0≤∥F1−F2∥V∗∥u1−u2∥V.

Also, since we deduce that

 ∥u1−u2∥V≤1cP−α0∥F1−F2∥V∗,

which concludes the proof of the theorem. ∎

Further, we present a class of locally Lipschitz functions for which the crucial assumption (17) is satisfied. Let be a function such that

 (ξ∗−η∗)(ξ−η)≥−α0|ξ−η|2∀ξ∗∈∂f(ξ),∀η∗∈∂f(η) (19)

for any and some . From the definition of the Clarke generalized derivative Clarke we get

 f0(ξ;η−ξ)=maxξ∗∈∂f(ξ)ξ∗(η−ξ).

Rewriting (19) as

 ξ∗(η−ξ)+η∗(ξ−η)≤α0|ξ−η|2

we find

 f0(ξ;η−ξ)+f0(η;ξ−η)≤α0|ξ−η|2. (20)

Hence,

 φ(u,v)+φ(v,u) =∫ΓCf0(un;vn−un)ds+∫ΓCf0(vn;un−vn)ds ≤α0∥un−vn∥2L2(ΓC)≤α0∥u−v∥2V (21)

and consequently the assumption (17) is satisfied provided that is sufficiently small ().

Next, we show that if includes only non-negative jumps, then the condition (19) is globally satisfied. Whereas for negative jumps, (19) holds only locally.

Example 1

Let be a continuous, piecewise function such that its first derivative has finite non-negative jumps at the points of discontinuity (see Fig. 3). This means that , where

 b–(xJi):=f′(xJi−0)=limh→0−f(xJi+h)−f(xJi)h

and

 ¯¯b(xJi):=f′(xJi+0)=limh→0+f(xJi+h)−f(xJi)h.

Hence,

 ∂f(xJi)=[b–(xJi),¯¯b(xJi)].

Setting the finite set of open intervals , where the function is smooth with the Lipschitz constant , we define

 α0=maxicLi. (22)

Let be the set of jags between and with , where . Then, for any we have the following

 b(x1)−b(x2)x1−x2 = b(x1)−b(x2)+k∑i=1¯¯b(xJi)−k∑i=1¯¯b(xJi)x1−x2 ≥ b(x1)−b(x2)+k∑i=1b–(xJi)−k∑i=1¯¯b(xJi)x1−x2 = b(x1)−¯¯b(xJk)−b(x2)+b–(xJ1)+k∑i=2(b–(xJi)−¯¯b(xJi−1))x1−x2 ≥ −α0(x1−xJk)−α0(xJ1−x2)−α0k∑i=2(xJi−xJi−1)x1−x2 = −α0(x1−x2)x1−x2=−α0,

from which the assumption (19) follows immediately.

Remark 2

If the graph of consists of several decreasing straight line segments and non-negative jumps, then the value in (19) is the steepest decreasing slope.

The next example treats the case of negative jumps.

Example 2

Let be a continuous, piecewise function such that its first derivative has finite negative jumps at the points of discontinuity (see Fig. 4). In this case,

 ∂f(xJi)=[¯¯b(xJi),b–(xJi)].

Then, for any , we have

 b(x1)−b(x2)x1−x2≥−α0+k∑i=1cJix1−x2, (23)

with as defined in (22) and standing for the value of the jump in the point .

Further, for the sake of simplicity, we consider a function with one negative jump and estimate as follows:

 cJ(x1−x2)=cJ(x1−xJ)+cJ(xJ−x2)=cJ(x1−xJ)++cJ(xJ−x2)+. (24)

Using for any , we get

 (x1−x2)2≥(x1−x2−14)+≥(x1−x2−(xJ−x2))+=(x1−xJ)+∀x2≤xJ−14

as well as

 (x1−x2)2≥(x1−x2−14)+≥(x1−x2−(x1−xJ))+=(xJ−x2)+∀x1≥xJ+14.

Hence, the condition (19) is fulfilled with , but only for those satisfying and