Weak imposition of Signorini boundary conditions on the boundary element method
aa r X i v : . [ m a t h . NA ] M a y WEAK IMPOSITION OF SIGNORINI BOUNDARY CONDITIONSON THE BOUNDARY ELEMENT METHOD ∗ ERIK BURMAN † , STEFAN FREI ‡ , AND
MATTHEW W. SCROGGS § Abstract.
We derive and analyse a boundary element formulation for boundary conditionsinvolving inequalities. In particular, we focus on Signorini contact conditions. The Calder´on projec-tor is used for the system matrix and boundary conditions are weakly imposed using a particularvariational boundary operator designed using techniques from augmented Lagrangian methods. Wepresent a complete numerical a priori error analysis and present some numerical examples to illus-trate the theory.
Key words. boundary element methods, Nitsche’s method, Signorini problem, Calder´on pro-jector
AMS subject classifications.
1. Introduction.
The application of Nitsche techniques to deal with variationalinequalities has received increasing interest recently, starting from a series of worksby Chouly, Hild and Renard for elasticity problems with contact [7]. Their approachgoes back to an augmented Lagrangian formulation, that has first been introduced byAlart & Curnier [1].In a previous paper [2], we have shown how Nitsche techniques can be used toimpose Dirichlet, Neumann, mixed Dirichlet–Neumann or Robin conditions weaklywithin boundary element methods. By using the Calder´on projector, we were able toderive a unified framework that can be used for different boundary conditions.The purpose of this article is to extend these techniques to boundary conditionsinvolving inequalities, such as Signorini contact conditions. In particular, we considerthe Laplace equation with mixed Dirichlet and Signorini boundary conditions: Find u such that − ∆ u = 0 in Ω , (1.1a) u = g D on Γ D , (1.1b) u g C and ∂u∂ ν ψ C on Γ C , (1.1c) (cid:18) ∂u∂ ν − ψ C (cid:19)(cid:18) u − g C (cid:19) = 0 on Γ C . (1.1d)Here Ω ⊂ R denotes a polyhedral domain with outward pointing normal ν andboundary Γ := Γ D ∪ Γ C . We assume for simplicity that the boundary between Γ D and Γ C coincides with edges between the faces of Γ. Whenever it is ambiguous,we will write ν x for the outward pointing normal at the point x . We assume that g = ( g D in Γ D g C in Γ C ∈ L (Γ) and ψ C ∈ H / (Γ C ). ∗ Submitted to the editors 2019-08-15.
Funding:
Erik Burman was funded by the EPSRC grant EP/P01576X/1. Stefan Frei wasfunded by the DFG Research Scholarship 3935/1-1. † Department of Mathematics, University College London, UK ([email protected]). ‡ Department of Mathematics and Statistics, University of Konstanz, Germany ([email protected]). § E. BURMAN, S. FREI, M. W. SCROGGS
Observe that when Γ C = ∅ , there exists a unique solution to (1.1) by the Lax–Milgram lemma. In the case that meas(Γ C ) >
0, the theory of Lions and Stampacchia[12] for variational inequalities yields existence and uniqueness of solutions. We as-sume that u ∈ H / ǫ (Ω), for some ǫ > L (Γ)-error estimatebased on a duality argument. Maischak & Stephan [13] presented a posteriori errorestimates and an hp -adaptive algorithm for the Signorini problem. A priori errorestimates for a penalty-based hp algorithm were shown by Chernov, Maischak &Stephan [6]. Recently, an augmented Lagrangian approach has been presented incombination with a semi-smooth Newton method [22], and variational inequalitieshave been successfully used for time-dependent contact problems [9].We will consider an approach where the full Calder´on projector is used and theboundary conditions are included by adding properly scaled penalty terms to the twoequations. This results in formulations similar to the ones obtained for weak impo-sition of boundary conditions using Nitsche’s method [14]. The proposed frameworkis flexible and allows for the design of a range of different methods depending on thechoice of weights and residuals.An outline of the paper is as follows. In Section 2, we introduce the basic boundaryoperators that will be needed and review some of their properties. Then, in Section 3,we introduce the variational framework and review the results from [2] for the pureDirichlet problem. In Section 4, we show how the framework can be applied toSignorini boundary conditions and the mixed problem (1.1). The method is analysedin Section 5. We conclude by showing some numerical experiments in Section 6.
2. Boundary operators.
We define the Green’s function for the Laplace oper-ator in R by(2.1) G ( x , y ) = 14 π | x − y | . In this paper, we focus on the problem in R . Similar analysis can be used for problemsin R , in which case this definition should be replaced by G ( x , y ) = − log | x − y | / π .In the standard fashion (see e.g. [19, Chapter 6]), we define the single layerpotential operator, V : H − / (Γ) → H (Ω), and the double layer potential, K : H / (Γ) → H (Ω), for v ∈ H / (Γ), µ ∈ H − / (Γ), and x ∈ Ω \ Γ by( V µ )( x ) := Z Γ G ( x , y ) µ ( y ) d y , (2.2) ( K v )( x ) := Z Γ ∂G ( x , y ) ∂ ν y v ( y ) d y . (2.3)We define the space H (∆ , Ω) := { v ∈ H (Ω) : ∆ v ∈ L (Ω) } , and the Dirichlet andNeumann traces, γ D : H (Ω) → H / (Γ) and γ N : H (∆ , Ω) → H − / (Γ), by γ D f ( x ) := lim Ω ∋ y → x ∈ Γ f ( y ) , (2.4) γ N f ( x ) := lim Ω ∋ y → x ∈ Γ ν x · ∇ f ( y ) . (2.5) EAK IMPOSITION OF SIGNORINI BOUNDARY CONDITIONS ON BEM u = −K ( γ D u ) + V ( γ N u ) . It is also known [19, Lemma 6.6] that for all µ ∈ H − / (Γ), the function(2.7) u V µ := V µ satisfies − ∆ u V µ = 0 and(2.8) k u V µ k H (Ω) c k µ k H − / (Γ) . Similarly [19, Lemma 6.10], the function(2.9) u K v := K v satisfies − ∆ u K v = 0 for all v ∈ H / (Γ) and(2.10) k u K v k H (Ω) c k v k H / (Γ) . We define { γ D f } Γ and { γ N f } Γ to be the averages of the interior and exteriorDirichlet and Neumann traces of f . We define the single layer, double layer, ad-joint double layer, and hypersingular boundary integral operators, V : H − / (Γ) → H / (Γ), K : H / (Γ) → H / (Γ), K ′ : H − / (Γ) → H − / (Γ), and W : H / (Γ) → H − / (Γ), by ( K v )( x ) := { γ D K v } Γ ( x ) , ( V µ )( x ) := { γ D V µ } Γ ( x ) , (2.11a) ( W v )( x ) := − { γ N K v } Γ ( x ) , ( K ′ µ )( x ) := { γ N V µ } Γ ( x ) , (2.11b)where x ∈ Γ, v ∈ H / (Γ) and µ ∈ H − / (Γ) [19, Chapter 6].Next, we define the Calder´on projector by(2.12) C := (cid:18) (1 − σ ) Id − K VW σ Id + K ′ (cid:19) , where σ is defined for x ∈ Γ by [19, Equation 6.11](2.13) σ ( x ) = lim ǫ → π ǫ Z y ∈ Ω: | y − x | = ǫ d y . Recall that if u is a solution of (1.1) then it satisfies(2.14) C (cid:18) γ D uγ N u (cid:19) = (cid:18) γ D uγ N u (cid:19) . Taking the product of (2.14) with two test functions, and using the fact that σ = almost everywhere, we arrive at the following equations. h γ D u, µ i Γ = (cid:10) ( Id − K ) γ D u, µ (cid:11) Γ + h V γ N u, µ i Γ ∀ µ ∈ H − / (Γ) , (2.15) h γ N u, v i Γ = (cid:10) ( Id + K ′ ) γ N u, v (cid:11) Γ + h W γ D u, v i Γ ∀ v ∈ H / (Γ) . (2.16) E. BURMAN, S. FREI, M. W. SCROGGS
For a more compact notation, we introduce λ = γ N u and u = γ D u and theCalder´on form(2.17) C [( u, λ ) , ( v, µ )] := (cid:10) ( Id − K ) u, µ (cid:11) Γ + h V λ, µ i Γ + (cid:10) ( Id + K ′ ) λ, v (cid:11) Γ + h W u, v i Γ . We may then rewrite (2.15) and (2.16) as(2.18) C [( u, λ ) , ( v, µ )] = h u, µ i Γ + h λ, v i Γ . We will also frequently use the multitrace form, defined by(2.19) A [( u, λ ) , ( v, µ )] := − h K u, µ i Γ + h V λ, µ i Γ + h K ′ λ, v i Γ + h W u, v i Γ . Using this, we may rewrite (2.18) as(2.20) A [( u, λ ) , ( v, µ )] = h u, µ i Γ + h λ, v i Γ . To quantify the two traces we introduce the product space V := H / (Γ) × H − / (Γ)and the associated norm k ( v, µ ) k V := k v k H / (Γ) + k µ k H − / (Γ) . The continuity and coercivity of A are immediate consequences of the propertiesof the operators V , K , K ′ and W : Lemma
There exists
C > such that |A [( w, η ) , ( v, µ )] | C k ( w, η ) k V k ( v, µ ) k V ∀ ( w, η ) , ( v, µ ) ∈ V . There exists α > such that α (cid:16) | v | H / ∗ (Γ) + k µ k H − / (Γ) (cid:17) A [( v, µ ) , ( v, µ )] ∀ ( v, µ ) ∈ V . Proof.
See [2].
3. Discretisation and weak imposition of Dirichlet boundary condi-tions.
In this section, we introduce the discrete spaces and review briefly how (non-homogeneous) Dirichlet boundary conditions can be imposed weakly within the vari-ational formulations introduced above. For a detailed derivation, and for differentboundary conditions, we refer to [2].To reduce the number of constants that appear, we introduce the following nota-tion. • If ∃ C > a Cb , then we write a . b . • If a . b and b . a , then we write a h b .We assume that Ω is a polygonal domain with faces denoted by { Γ i } Mi =1 . Weintroduce a family of conforming, shape regular triangulations of Γ, {T h } h> , indexedby the largest element diameter of the mesh, h . We let T , .., T m ∈ T h be the trianglesof a triangulation. EAK IMPOSITION OF SIGNORINI BOUNDARY CONDITIONS ON BEM Figure 1 . A grid (left), the barycentric refinement of the grid (centre), and the dual grid(right). In a typical example, the initial grid will not be flat, and so the elements of the dual gridwill not necessarily be flat.
We consider the following finite element spacesP kh (Γ) := { v h ∈ C (Γ) : v h | T i ∈ P k ( T i ), for every T i ∈ T h } , DP lh (Γ) := { v h ∈ L (Γ) : v h | T i ∈ P l ( T i ), for every T i ∈ T h } , g DP lh (Γ) := { v h ∈ DP lh (Γ) : v h | Γ i ∈ C (Γ i ), for i = 1 , . . . , M } , where P k ( T i ) denotes the space of polynomials of order less than or equal to k on thetriangle T i .In addition, we consider the space DUAL h (Γ) of piecewise constant functions onthe barycentric dual grid, as shown in Figure 1. On non-smooth domains, these spaceshave lower order approximation properties than the standard space DP h (Γ), as givenin the following lemma. Lemma
Let µ ∈ H s (Γ) . If Γ consists of a finite number of smooth facesmeeting at edges, then inf η h ∈ DUAL h (Γ) k µ − η h k H − / (Γ) . h ξ +1 / k µ k H ξ (Γ) where ξ = min( , s ) . If Γ is smooth, then the same result holds with ξ = min(1 , s ) .Proof. See [16, Appendix 2].We observe that P kh (Γ) ⊂ H / (Γ), DP lh (Γ) ⊂ L (Γ), g DP lh (Γ) ⊂ L (Γ), andDUAL h (Γ) ⊂ L (Γ). We define the discrete product space V h := P kh (Γ) × Λ lh , where Λ lh can be any of the spaces DP lh (Γ) , g DP lh (Γ) or DUAL h (Γ). Let us for the moment assume thatΓ ≡ Γ D . Then, the basic idea is to add the following suitably weighted boundaryresidual to the weak formulation.(3.1) R Γ D ( u h , λ h ) := β / ( g D − u h ) . This is defined such that R Γ D ( u h , λ h ) = 0 is equivalent to the boundary condition(1.1b). We obtain an expression of the form(3.2) C [( u h , λ h ) , ( v h , µ h )] = h u h , µ h i Γ + h λ h , v h i Γ + h R Γ D ( u h , λ h ) , β v h + β µ h i Γ , E. BURMAN, S. FREI, M. W. SCROGGS or equivalently(3.3) A [( u h , λ h ) , ( v h , µ h )] = h u h , µ h i Γ + h λ h , v h i Γ + h R Γ D ( u h , λ h ) , β v h + β µ h i Γ , where β and β are problem dependent scaling operators that can be chosen as afunction of the physical parameters in order to obtain robustness of the method.For the Dirichlet problem, we choose β = β / , β = β − / , where differentchoices for β D in the range 0 β D . h − are possible. Inserting this into (3.3), weobtain the formulation:(3.4) A [( u, λ ) , ( v h , µ h )] − h λ h , v h i Γ D + h u h , µ h i Γ D + h β D u h , v h i Γ D = h g D , β D v h + µ h i Γ D . By formally identifying λ h with ∂ ν u h and µ h with ∂ ν v h , we obtain the classical (non-symmetric) Nitsche’s method (up to the multiplicative factor ).For a more compact notation, we introduce the boundary operator associatedwith the non-homogeneous Dirichlet condition(3.5) B D [( u h , λ h ) , ( v h , µ h )] := − h λ h , v h i Γ D + h u h , µ h i Γ D + h β D u h , v h i Γ D , the operator corresponding to the left-hand side(3.6) A D [( u h , λ h ) , ( v h , µ h )] := A [( u h , λ h ) , ( v h , µ h )] + B D [( u h , λ h ) , ( v h , µ h )]and the operator associated with the right-hand side(3.7) L D ( v h , µ h ) := h g D , β D v h + µ h i Γ D . Using these and (3.4), we arrive at the following boundary element formulation:Find ( u h , λ h ) ∈ V h such that A D [( u h , λ h ) , ( v h , µ h )] = L D ( v h , µ h ) ∀ ( v h , µ h ) ∈ V h . (3.8)We introduce the following B D -norm k ( v, µ ) k B D := k ( v, µ ) k V + β / k v k Γ D , and summarise the properties of the bilinear form A D in the following lemma. Lemma
Let W be a product Hilbert spacefor the primal and flux variables, such that V ⊂ W . The bilinear form has the followingproperties: Property 1 (Coercivity): If β D = 0 or if there exists β min > (independent of h ) such that β D > β min , then there exists α > such that ∀ ( v, µ ) ∈ W α k ( v, µ ) k B D A D [( v, µ ) , ( v, µ )] . Property 2 (Continuity):
There exists
M > such that ∀ ( w, η ) , ( v, µ ) ∈ W |A D [( v, µ ) , ( w, η )] | M k ( v, µ ) k B D k ( w, η ) k B D . Proof.
See [2, Section 4.1].
EAK IMPOSITION OF SIGNORINI BOUNDARY CONDITIONS ON BEM
4. Weak imposition of Signorini boundary conditions.
Recently Chouly,Hild and Renard [7, 8] showed how contact problems can be treated in the contextof Nitsche’s method. We will here show how we may use arguments similar to theirsin the present framework to integrate unilateral contact seamlessly. The result is anonlinear system to which one may apply Newton’s method or a fixed-point iterationin a straightforward manner. We prove existence and uniqueness of solutions to thenonlinear system and optimal order error estimates.For the derivation of the formulation on the contact boundary we will first omitthe Dirichlet part, letting Γ = Γ C . To impose the contact conditions, we recall thefollowing relations, introduced by Alart and Curnier [1], with [ x ] ± := ± max(0 , ± x ).( u − g C ) = (cid:2) ( u − g C ) − τ − ( λ − ψ C ) (cid:3) − on Γ C , (4.1) ( λ − ψ C ) = − [ τ ( u − g C ) − ( λ − ψ C )] + on Γ C , (4.2)for all τ >
0. It is straighforward [7] to show that each of these two conditions isequivalent to the contact boundary conditions (1.1c) and (1.1d).To simplify the notation, we introduce the operators P τ ( u h , λ h ) := τ ( u h − g C ) − ( λ h − ψ C ) and P τ ( u h , λ h ) := τ u h − λ h . Using (4.1), we arrive at the following boundary term for the contact conditions(4.3) R C ( u h , λ h ) = ( g C − u h ) + τ − [ P τ ( u h , λ h )] − . Alternatively, by using (4.2), we arrive at the following boundary term(4.4) R C ( u h , λ h ) = τ − (cid:0) ( ψ C − λ h ) − [ P τ ( u h , λ h )] + (cid:1) . By using the fact that x = [ x ] + + [ x ] − , it can be shown that (4.3) and (4.4) are equal.Substituting (4.3) into (3.3), and using the weights β = τ and β = 1, we obtain(4.5) A [( u h , λ h ) , ( v h , µ h )] + h µ h , u h i Γ C + (cid:10) τ u h − λ h , v h (cid:11) Γ C − (cid:10) [ P τ ( u h , λ h )] − , v h + τ − µ h (cid:11) Γ C = h g C , τ v h + µ h i Γ C . Using (4.4), we have(4.6) A [( u h , λ h ) , ( v h , µ h )] + h λ h , v h i Γ C + (cid:10) τ − λ h − u h , µ h (cid:11) Γ C + (cid:10) [ P τ ( u h , λ h )] + , v h + τ − µ h (cid:11) Γ C = (cid:10) ψ C , v h + τ − µ h (cid:11) Γ C . We see that (4.6) is similar to the non-symmetric version of the method pro-posed in [8] and (4.5) is similar to the non-symmetric Nitsche formulation for contactdiscussed in [5]. As pointed out in the latter reference, the two formulations are equiv-alent, with the same solutions. In what follows, we focus exclusively on the variant(4.6).Defining B C [( u h , λ h ) , ( v h , µ h )] := h λ h , v h i Γ C + (cid:10) τ − λ h − u h , µ h (cid:11) Γ C + (cid:10) [ P τ ( u h , λ h )] + , v h + τ − µ h (cid:11) Γ C , (4.7) L C ( v h , µ h ) := (cid:10) ψ C , v h + τ − µ h (cid:11) Γ C , (4.8) A C [( u h , λ h ) , ( v h , µ h )] := A [( u h , λ h ) , ( v h , µ h )] + B C [( u h , λ h ) , ( v h , µ h )] , (4.9)we arrive at the boundary element method formulation: Find ( u h , λ h ) ∈ V h such that A C [( u h , λ h ) , ( v h , µ h )] = L C ( v h , µ h ) ∀ ( v h , µ h ) ∈ V h . (4.10) E. BURMAN, S. FREI, M. W. SCROGGS
Combining theformulations for the Dirichlet and contact conditions, we arrive at the following bound-ary element method for the problem (1.1): Find ( u h , λ h ) ∈ V h such that(4.11) A D [( u h , λ h ) , ( v h , µ h )] + B C [( u h , λ h ) , ( v h , µ h )] = L D ( v h , µ h ) + L C ( v h , µ h ) ∀ ( v h , µ h ) ∈ V h , where A D , L D , B C and L C are defined in (3.6), (3.7), (4.7), and (4.8). For dis-cretisation, we use the assumptions and spaces introduced in Section 3. Note thatthe formulation (4.11) is consistent, i.e. the continuous solution ( u, λ ) to (1.1) fulfills(4.11) for all ( v h , µ h ) ∈ V h .
5. Analysis.
In this section, we prove the existence of unique solutions to thenonlinear system of equations (4.11) as well as optimal error estimates.We assume that the solution ( u, λ ) of (1.1) lies in W := H ǫ (Γ) × H ǫ (˜Γ) forsome ǫ ∈ (0 , / ∪ Mi =1 Γ i \ ∂ Γ i is the set of boundary points that lie in theinterior of the faces Γ i . As the normal vectors ν x are discontinuous between faces, wecan not expect a higher global regularity for λ .We define the distance function d C and norm k · k ∗ , for ( v, µ ) , ( w, η ) ∈ W , by d C (( v, µ ) , ( w, η )) := k ( v − w, µ − η ) k B D + k τ − (cid:0) µ − η + [ P τ ( v, µ )] + − [ P τ ( w, η )] + (cid:1) k Γ C , (5.1) k ( v, µ ) k ∗ := k ( v, µ ) k B D + k τ v k Γ C + k τ − µ k Γ C . (5.2)We note that due the appearance of [ · ] + in its second term, d C is not a norm. d C does provide a bound on the error however, as for all ( v, µ ) ∈ W , d C (( v, µ ) , (0 , > k ( v, µ ) k B D > k ( v, µ ) k V .When proving this section’s results, we will use properties of the [ · ] + functionthat are given in the following lemma. Lemma
For all a, b ∈ R , (cid:0) [ a ] + − [ b ] + (cid:1) (cid:0) [ a ] + − [ b ] + (cid:1) ( a − b ) , (5.3) | [ a ] + − [ b ] + | | a − b | . (5.4) Proof.
For a proof of these well-known properties see e.g. [7].We now prove a result analogous to the coercivity assumption in [2].
Lemma
If there is β min > , independent of h , such that β D > β min , thenthere is α > such that for all ( v, µ ) , ( w, η ) ∈ W , α ( d C (( v, µ ) , ( w, η ))) ( A + B D )[( v − w, µ − η ) , ( v − w, µ − η )]+ B C [( v, µ ) , ( v − w, µ − η )] − B C [( w, η ) , ( v − w, µ − η )] . Proof.
From the analysis of the Dirichlet problem (Lemma 3.2) we know thatwhen β D > β min > α k ( v − w, µ − η ) k B D ( A + B D )[( v − w, µ − η ) , ( v − w, µ − η )] . Introducing the notation δP := [ P τ ( v, µ )] + − [ P τ ( w, η )] + , we have(5.6) B C [( v, µ ) , ( v − w, µ − η )] − B C [( w, η ) , ( v − w, µ − η )]= τ − k µ − η k C + (cid:10) δP, v − w + τ − ( µ − η ) (cid:11) Γ C . EAK IMPOSITION OF SIGNORINI BOUNDARY CONDITIONS ON BEM τ − k µ − η + δP k C = τ − (cid:0) k µ − η k C + k δP k C + 2 h µ − η, δP i Γ C (cid:1) . Using (5.3), this implies the bound τ − k µ − η + δP k C τ − (cid:0) k µ − η k C + h δP, P τ ( v − w, µ − η ) i Γ C + 2 h µ − η, δP i Γ C (cid:1) . Observing that P τ ( v − w, µ − η ) + 2( µ − η ) = τ ( v − w ) + µ − η , we infer that(5.7) τ − k µ − η + δP k C B C [( v, µ ) , ( v − w, µ − η ) − B C [( w, η ) , ( v − w, µ − η )] . We conclude the proof by noting that( d C (( v, µ ) , ( w, η ))) . k ( v − w, µ − η ) k B D + τ − k µ − η + [ P τ ( v, µ )] + − [ P τ ( w, η )] + k C , and applying (5.5) and (5.7).Next, we prove a result analagous to the discrete coercivity assumption in [2]. Lemma
If there is β min > , independent of h , such that β D > β min , thenthere is α > such that for all ( v h , µ h ) ∈ V h , α (cid:16) k ( v h , µ h ) k B D + k τ − (cid:0) µ h + [ P τ ( v h , µ h )] + (cid:1) k Γ C (cid:17) ( A + B D + B C )[( v h , µ h ) , ( v h , µ h )] − (cid:10) [ P τ ( v h , µ h )] + , g C − τ − ψ C (cid:11) Γ C Proof.
The proof is similar to that of Lemma 5.2, but with µ h and v h instead of µ − η and v − w . The appearance of the data term in the right-hand side is due tothe relation τ − k [ P τ ( v h , µ h )] + k C + 2 τ − (cid:10) µ h , [ P τ ( v h , µ h )] + (cid:11) Γ C + τ − k µ h k C = τ − (cid:10) [ P τ ( v h , µ h )] + , P τ ( v h , µ h ) (cid:11) Γ C + τ − k µ h k C = (cid:10) [ P τ ( v h , µ h )] + , u h + τ − µ h (cid:11) Γ C − (cid:10) [ P τ ( v h , µ h )] + , g C − τ − ψ C (cid:11) Γ C + τ − k µ h k C = B C [( v h , µ h ) , ( v h , µ h )] − (cid:10) [ P τ ( v h , µ h )] + , g C − τ − ψ C (cid:11) Γ C . Using Lemmas 5.2 and 5.3, we may now prove that (4.11) is well-posed.
Theorem
The finite dimensional nonlinear system (4.11) admits a uniquesolution.Proof.
To prove the existence of a solution, we show the continuity and the posi-tivity of the nonlinear operator A + B D + B C . This allows us to apply Brouwer’s fixedpoint theorem, see eg [21, Chapter 2, Lemma 1.4].We define F : V h → V h , for ( v h , µ h ) ∈ V h , by h F ( v h , µ h ) , ( w h , η h ) i Γ = ( A + B D + B C )[( v h , µ h ) , ( w h , η h )] − L D ( w h , η h ) − L C ( w h , η h ) , E. BURMAN, S. FREI, M. W. SCROGGS for all ( w h , η h ) ∈ V h . We may write the non-linear system (4.11) as h F ( v h , µ h ) , ( w h , η h ) i Γ = 0 ∀ ( w h , η h ) ∈ V h . (5.8)For fixed h , by the equivalance of norms on discrete spaces, there exist c , c > v h , µ h ) ∈ V h , c k ( v h , µ h ) k Γ k ( v h , µ h ) k B D c k ( v h , µ h ) k Γ . To show positivity, we let ( v h , µ h ) ∈ V h . Using Lemma 5.3, we see that h F ( v h , µ h ) , ( v h , µ h ) i Γ > α k ( v h , µ h ) k B D + ατ − k µ h + [ P τ ( v h , µ h )] + k C + (cid:10) [ P τ ( v h , µ h )] + , g C − τ − ψ C (cid:11) Γ C − L D ( v h , µ h ) − L C ( v h , µ h ) . Using the Cauchy–Schwarz inequality and an arithmetic-geometric inequality, we seethat there exists C g C ψ C > (cid:10) [ P τ ( v h , µ h )] + , g C − τ − ψ C (cid:11) Γ C − L D ( v h , µ h ) − L C ( v h , µ h )= (cid:10) [ P τ ( v h , µ h )] + + µ h , g C − τ − ψ C (cid:11) Γ C − (cid:10) µ h , g C − τ − ψ C (cid:11) Γ C − h g D , β D v h + µ h i Γ D − (cid:10) ψ C , v h + τ − µ h (cid:11) Γ C > − C g C ψ C − α (cid:0) k ( v h , µ h ) k B D + τ − k µ h + [ P τ ( v h , µ h )] + k C (cid:1) . Using norm equivalence, we obtain h F ( v h , µ h ) , ( v h , µ h ) i Γ > α (cid:0) k ( v h , µ h ) k B D + τ − k µ h + [ P τ ( v h , µ h )] + k C (cid:1) − C g C ψ C > C ′ k ( v h , µ h ) k − C g C ψ C , for some C ′ >
0. We conclude that for all ( v h , µ h ) ∈ V h with k ( v h , µ h ) k > C g C ψ C C ′ + 1 , there holds h F ( v h , µ h ) , ( v h , µ h ) i Γ > v h , µ h ) , ( v h , µ h ) ∈ V h . We have for all ( w h , η h ) ∈ V h , (cid:10) F ( v h , µ h ) − F ( v h , µ h ) , ( w h , η h ) (cid:11) Γ = D(cid:2) P τ ( v h , µ h ) (cid:3) + − (cid:2) P τ ( v h , µ h ) (cid:3) + , w h + τ − η h E Γ C + (cid:10) µ h − µ h , w h + τ − η h (cid:11) Γ − (cid:10) v h − v h , µ h − µ h (cid:11) Γ C + ( A + B D )[( v h − v h , µ h − µ h ) , ( w h , η h )] (cid:0) τ k v h − v h k Γ C + k µ h − µ h k Γ C (cid:1) (cid:0) k w h k Γ C + τ − k η h k Γ C (cid:1) , where we have used (5.4). By norm equivalence, this means that (cid:10) F ( v h , µ h ) − F ( v h , µ h ) , ( w h , η h ) (cid:11) Γ k ( w h , η h ) k Γ C k ( v h − v h , µ h − µ h ) k Γ showing that F is continuous. EAK IMPOSITION OF SIGNORINI BOUNDARY CONDITIONS ON BEM u h , λ h )and ( u h , λ h ) are solutions to (4.11). We immediately see that α (cid:0) d C (cid:0) ( u h , λ h ) , ( u h , λ h ) (cid:1)(cid:1) = 0 , and we conclude that the solution is unique.We now proceed to prove the following best approximation result. Lemma
Let ( u, λ ) ∈ W be the solution of (1.1) and ( u h , λ h ) ∈ V h the solutionof (4.11) . Then there holds d C (( u, λ ) , ( u h , λ h )) C inf ( v h ,µ h ) ∈ V h k ( u − v h , λ − µ h ) k ∗ . Proof.
Using Lemma 5.2 and Galerkin orthogonality, we see that, for arbitrary( v h , µ h ) ∈ V h , α ( d C (( u, λ ) , ( u h , λ h ))) ( A + B D )[( u − u h , λ − λ h ) , ( u − u h , λ − λ h )]+ B C [( u, λ ) , ( u − u h , λ − λ h )] − B C [( u h , λ h ) , ( u − u h , λ − λ h )]= ( A + B D )[( u − u h , λ − λ h ) , ( u − v h , λ − µ h )]+ B C [( u, λ ) , ( u − v h , λ − µ h )] − B C [( u h , λ h ) , ( u − v h , λ − µ h )] . Next, we use B C [( u, λ ) , ( u − v h , λ − µ h )] − B C [( u h , λ h ) , ( u − v h , λ − µ h )]= (cid:10) λ − λ h + [ P τ ( u, λ )] + − [ P τ ( u h , λ h )] + , ( u − v h ) + τ − ( λ − µ h ) (cid:11) Γ C − h u − u h , λ − µ h i Γ C − h λ − λ h , u − v h i Γ C to show that( A + B D )[( u − u h , λ − λ h ) , ( u − u h , λ − λ h )]+ B C [( u, λ ) , ( u − u h , λ − λ h )] − B C [( u h , λ h ) , ( u − u h , λ − λ h )]= ( A + B D )[( u − u h , λ − λ h ) , ( u − v h , λ − µ h )] | {z } (I) − h u − u h , λ − µ h i Γ C − h λ − λ h , u − v h i Γ C | {z } (II) + (cid:10) λ − λ h + [ P τ ( u, λ )] + − [ P τ ( u h , λ h )] + , ( u − v h ) + τ − ( λ − µ h ) (cid:11) Γ C | {z } (III) . We estimate the three parts of the right-hand separately. For the first term, weuse the continuity of A + B D (Lemma 3.2) to obtain(I) M k ( u − u h , λ − λ h ) k B D k ( u − v h , λ − µ h ) k B D . For the second line, we use H / (Γ)– H − / (Γ) duality and the Cauchy–Schwarz in-equality to obtain(II) k ( u − u h , λ − λ h ) k B D k ( u − v h , λ − µ h ) k B D . E. BURMAN, S. FREI, M. W. SCROGGS
For the last term, we use the Cauchy–Schwarz inequality to get(III) k τ − / (cid:0) λ − λ h + [ P τ ( u, λ )] + − [ P τ ( u h , λ h )] + (cid:1) k Γ C · (cid:16) k τ / ( u − v h ) k Γ C + k τ − / ( λ − µ h ) k Γ C (cid:17) . Collecting these bounds, we see that d C (( u, λ ) , ( u h , λ h )) . d C (( u, λ ) , ( u h , λ h )) k ( u − v h , λ − µ h ) k ∗ . Dividing through by d C (( u, λ ) , ( u h , λ h )), and taking the infimum yields the desiredresult.We now prove the main result of this section, an a priori bound on the error ofthe solution of (4.11). Theorem
Let ( u, λ ) ∈ H s (Γ) × H r (˜Γ) for some s > , r > and ( u h , λ h ) ∈ P kh (Γ) × Λ lh be the solutions of (1.1) and the discrete problem (4.11) , respectively. Ifthere is β min > such that β min < β D . h − and τ h h − , then k ( u − u h , λ − λ h ) k V d C (( u, λ ) , ( u h , λ h )) . h ζ − / | u | H ζ (Γ) + h ξ +1 / | λ | H ξ (˜Γ) , where ζ = min( k + 1 , s ) and ξ = min( l + 1 , r ) for Λ lh ∈ { DP lh (Γ) , g DP lh (Γ) } and ζ = min(2 , s ) and ξ = min( , r ) for Λ lh = DUAL h (Γ) . Additionally, k ˜ u − ˜ u h k H (Ω) . h ζ − / | u | H ζ (Γ) + h ξ +1 / | λ | H ξ (˜Γ) , where ˜ u and ˜ u h are the solutions in Ω defined by (2.6) .Proof. First, we observe that for all ( v, µ ) and ( w, η ) in W k ( v − w, µ − η ) k V d C (( v, µ ) , ( w, η )) . Using standard approximation results for Λ lh ∈ { DP lh (Γ) , g DP lh (Γ) } (see eg [19, chapter10]) and Lemma 3.1 for Λ lh = DUAL h (Γ), we see thatinf ( v h ,µ h ) ∈ V h k ( u − v h , λ − µ h ) k V = inf v h ∈ P kh (Γ) k u − v h k H / (Γ) + inf µ h ∈ Λ lh (Γ) k λ − µ h k H − / (Γ) . h ζ − / | u | H ζ (Γ) + h ξ +1 / | λ | H ξ (˜Γ) , inf v h ∈ P kh (Γ) k u − v h k Γ . h ζ | u | H ζ (Γ) , inf µ h ∈ Λ lh k λ − µ h k Γ . h ξ | λ | H ξ (˜Γ) . Applying these to the definition of k · k ∗ givesinf ( v h ,µ h ) ∈ V h k ( u − v h , λ − µ h ) k ∗ . h ζ − / | u | H ζ (Γ) + h ξ +1 / | λ | H ξ (˜Γ) + β / h ζ | u | H ζ (Γ) + τ / h ζ | u | H ζ (Γ) + τ − / h ξ | λ | H ξ (˜Γ) . By means of Lemma 5.5 and the given choice of the parameters τ and β D this provesthe first assertion. The estimate in the domain Ω follows by using the relations (2.8)and (2.10).If λ is smooth enough and k = l , the bounds on τ can be replaced with h . τ . h − without reducing the order of convergence. EAK IMPOSITION OF SIGNORINI BOUNDARY CONDITIONS ON BEM
6. Numerical results.
We now demonstrate the theory with a series of nu-merical examples. In this section, we consider the following test problem. LetΩ = [0 , × [0 , × [0 ,
1] be the unit cube, Γ C := { ( x, y, z ) ∈ Γ : z = 1 } , andΓ D := Γ \ Γ C . Let g D = 0 , (6.1a) g C = ( sin( πx ) sin( πy ) sinh( √ π ) x sin( πy ) sinh( √ π ) x > , (6.1b) ψ C = ( √ π sin( πx ) sin( πy ) cosh( √ π ) x > √ π sin( πy ) cosh( √ π ) x < . (6.1c)It can be shown that u ( x, y, z ) = sin( πx ) sin( πy ) sinh( √ πz )is the solution to (1.1) with these boundary conditions.To solve the non-linear system (4.10), we will treat the nonlinear term explicitly.Therefore, we define B ′ C [( u, λ ) , ( v, µ )] := h λ, v i Γ C + (cid:10) τ − λ − u, µ (cid:11) Γ C (6.2)Note that B ′ C differs from B C only by the missing nonlinear term.We pick initial guesses ( u , λ ) ∈ V h and define ( u n +1 , λ n +1 ) ∈ V h , for n ∈ N , tobe the solution of(6.3) ( A + B D + B ′ C )[( u n +1 , λ n +1 ) , ( v h , µ h )]= L C ( v h , µ h ) − (cid:10) [ P τ ( u n , λ n )] + , v h + τ − µ h (cid:11) Γ C ∀ ( v h , µ h ) ∈ V h . This leads us to Algorithm 6.1, an iterative method for solving the contact problem.In all the computations in this section, we preconditioned the GMRES solverusing a mass matrix preconditioner applied blockwise from the left, as described in[3].
Algorithm 6.1
Iterative algorithm for solving the contact problemInput ( u , λ ), tol , maxiter for n ← maxiter do ( u n +1 , λ n +1 ) ← solution of (6.3), calculated using GMRES if k ( u n +1 , λ n +1 ) − ( u n , λ n ) k V < tol thenreturn ( u n +1 , λ n +1 ) end ifend for Inspired by the parameter choices in [2], we fix β D = 0 .
01 and look for suitablevalues of the parameter τ . Figure 2 shows how the error, number of outer iterations,and the average number of GMRES iterations inside each outer iteration changeas the parameter τ is varied, for both V h = P h (Γ) × DUAL h (Γ) (left, blue) and V h = P h (Γ) × DP h (Γ) (right, orange). Here, we see that the error and number ofouter iterations are lowest when τ is between around 1 and 10.4 E. BURMAN, S. FREI, M. W. SCROGGS10 − − τ E rr o r i n V n o r m − − τ E rr o r i n V n o r m − − τ № o f o u t e r i t e r a t i o n s − − τ № o f o u t e r i t e r a t i o n s − − τ A v e r ag e № o f G M R E S i t e r a t i o n s − − , , , τ A v e r ag e № o f G M R E S i t e r a t i o n s Figure 2 . The dependence of the error, number of outer iterations, and the average number ofGMRES iterations on τ , for the problem (1.1) with boundary conditions (6.1) on the unit cube with h = 2 − (triangles), h = 2 − . (diamonds), and h = 2 − (pentagons). Here we take u = λ = 0 , β D = 0 . , tol = 0 . , and maxiter = 50 . On the left (blue), we take ( u n , λ n ) , ( v h , µ h ) ∈ P h (Γ) × DUAL h (Γ) ; on the right (orange), we take ( u n , λ n ) , ( v h , µ h ) ∈ P h (Γ) × DP h (Γ) . Motivated by Figure 2 and the bounds in Theorem 5.6, we take τ = 0 . /h , andlook at the convergence as h is decreased. Figure 3 shows how the error and iterationcounts vary as h is decreased when V h = P h (Γ) × DUAL h (Γ) (left, blue circles) and V h = P h (Γ) × DP h (Γ) (right, orange squares).For V h = P h (Γ) × DUAL h (Γ), we observe slightly higher than the order 1 con- EAK IMPOSITION OF SIGNORINI BOUNDARY CONDITIONS ON BEM − − h E rr o r i n V n o r m − − h E rr o r i n V n o r m − − h № o f o u t e r i t e r a t i o n s − − h № o f o u t e r i t e r a t i o n s − − h A v e r ag e № o f G M R E S i t e r a t i o n s − − h A v e r ag e № o f G M R E S i t e r a t i o n s Figure 3 . The error, number of outer iterations and averge number of inner GMRES iterationfor the problem (1.1) with boundary conditions (6.1) on the unit cube as h is reduced. Here wetake u = λ = 0 , β D = 0 . , tol = 0 . , maxiter = 200 , and τ = 0 . /h . On the left (bluecircles), we take ( u n , λ n ) , ( v h , µ h ) ∈ P h (Γ) × DUAL h (Γ) ; on the right (orange squares), we take ( u n , λ n ) , ( v h , µ h ) ∈ P h (Γ) × DP h (Γ) . The dashed lines show order 1 convergence (left) and order1.5 convergence (right). vergence predicted by Theorem 5.6. In this case, the mass matrix preconditioner iseffective, as the number of GMRES iterations required inside each outer iteration isreasonably low, and only grows slowly as h is decreased. We believe that the effec-tiveness of the preconditioner for this choice of spaces is due to the spaces P h (Γ) andDUAL h (Γ) forming an inf-sup stable pair [18, Lemma 3.1].6 E. BURMAN, S. FREI, M. W. SCROGGS
When V h = P h (Γ) × DP h (Γ), Theorem 5.6 tells us to expect order 1.5 conver-gence. However, we observe a slightly lower order. This appears to be due to theill-conditioning of this system, and the mass matrix preconditioner being ineffective,leading to an inaccurate solution when using GMRES. In this case, the spaces P h (Γ)and DP h (Γ) do not form an inf-sup stable pair, and so the mass-matrix between themis not guaranteed to be invertible leading to a less effective preconditioner.In order to obtain order 1.5 convergence with a well-conditioned system, we couldlook for ( u h , λ h ) ∈ P h (Γ) × DP h (Γ) and test with ( v h , µ h ) ∈ DUAL h (Γ) × DUAL h (Γ),where DUAL h (Γ) is the space of piecewise linear functions on the dual grid that formsan inf-sup stable pair with the space DP h (Γ), as defined in [4]. With this choice ofspaces, we obtain the higher order convergence as in Theorem 5.6, while having stabledual pairings and hence more effective mass matrix preconditioning.For the problems discussed in [2], we have run numerical experiments using thisspace pairing and observe the full order convergence in a low number of iterations.A deeper investigation of this method using these dual spaces, and the adaption ofthe theory to this case, warrants future work.
7. Conclusions.
Based on our work in [2], we have analysed and demonstratedthe effectiveness of Nitsche type coupling methods for boundary element formulationsof contact problems.An open problem is preconditioning. While the iteration counts in the presentedexamples were already practically useful, for large and complex structures precondi-tioning is still essential. The hope is to use the properties of the Calder´on projectorto build effective operator preconditioning techniques for the presented Nitsche typeframeworks.Avenues of future research include looking at how this approach could be appliedto problems in linear elasticity, and an extension of this method to problems involvingfriction.
REFERENCES[1]
P. Alart and A. Curnier , A mixed formulation for frictional contact problems prone toNewton like solution methods , Computer Methods in Applied Mechanics and Engineering,92 (1991), pp. 353–375.[2]
T. Betcke, E. Burman, and M. W. Scroggs , Boundary element methods with weakly imposedboundary conditions , SIAM Journal on Scientific Computing, 41 (2019), pp. A1357–A1384.[3]
T. Betcke, M. W. Scroggs, and W. ´Smigaj , Product algebras for Galerkin discretizationsof boundary integral operators and their applications . submitted to ACM Transactions onMathematical Software, 2018.[4]
A. Buffa and S. H. Christiansen , A dual finite element complex on the barycentric refine-ment , Mathematics of Computation, 76 (2007), pp. 1743–1769.[5]
E. Burman, P. Hansbo, and M. G. Larson , The penalty-free Nitsche method and noncon-forming finite elements for the Signorini problem , SIAM Journal on Numerical Analysis,55 (2017), pp. 2523–2539.[6]
A. Chernov, M. Maischak, and E. Stephan , A priori error estimates for hp penalty BEM forcontact problems in elasticity , Computer Methods in Applied Mechanics and Engineering,196 (2007), pp. 3871–3880.[7]
F. Chouly and P. Hild , A Nitsche-based method for unilateral contact problems: numericalanalysis , SIAM Journal on Numerical Analysis, 51 (2013), pp. 1295–1307.[8]
F. Chouly, P. Hild, and Y. Renard , Symmetric and non-symmetric variants of Nitsche’smethod for contact problems in elasticity: theory and numerical experiments , Mathematicsof Computation, 84 (2015), pp. 1089–1112.[9]
H. Gimperlein, F. Meyer, C. ¨Ozdemird, and E. P. Stephan , Time domain boundary ele-ments for dynamic contact problems , Computer Methods in Applied Mechanics and Engi-neering, 333 (2018), pp. 147–175.EAK IMPOSITION OF SIGNORINI BOUNDARY CONDITIONS ON BEM [10] H. Han , A direct boundary element method for Signorini problems , Mathematics of Computa-tion, 55 (1990), pp. 115–128.[11]
H.-d. Han , The boundary finite element methods for Signorini problems , in Numerical Methodsfor Partial Differential Equations, Y.-I. Zhu and B.-Y. Guo, eds., Springer, 1987, pp. 38–49.[12]
J. L. Lions and G. Stampacchia , Variational inequalities , Communications on Pure andApplied Mathematics, 20 (1967), pp. 493–519.[13]
M. Maischak and E. P. Stephan , Adaptive hp -versions of BEM for Signorini problems ,Applied Numerical Mathematics, 54 (2005), pp. 425 – 449.[14] J. Nitsche , ¨Uber ein Variationsprinzip zur L¨osung von Dirichlet-Problemen bei Verwendungvon Teilr¨aumen, die keinen Randbedingungen unterworfen sind , Abhandlungen aus demMathematischen Seminar der Universit¨at Hamburg, 36 (1971), pp. 9–15.[15] H. Schmit and G. Schneider , Boundary element solution of the Dirichlet-Signorini problemby a penalty method , Applicable Analysis, 51 (1993), pp. 175–186.[16]
M. W. Scroggs , Efficient computation and applications of the Calder´on projector , PhD thesis,University College London, 2019.[17]
W. Spann , On the boundary element method for the Signorini problem of the Laplacian , Nu-merische Mathematik, 65 (1993), pp. 337–356.[18]
O. Steinbach , On a generalized L projection and some related stability estimates in Sobolevspaces , Numer Math, 90 (2002), pp. 775–786.[19] O. Steinbach , Numerical approximation methods for elliptic boundary value problems ,Springer, 2008. Finite and boundary elements.[20]
O. Steinbach , Boundary element methods for variational inequalities , Numerische Mathe-matik, 126 (2014), pp. 173–197.[21]
R. Temam , Navier-Stokes equations: Theory and numerical analysis , vol. 2 of Studies in math-ematics and its applications, North-Holland Publishing, 1977.[22]
S. Zhang and X. Li , An augmented Lagrangian method for the Signorini boundary valueproblem with BEM , Boundary Value Problems, 2016 (2016), p. 62.[23]