Solving 2D Linear Isotropic Elastodynamics by Means of Scalar Potentials: A New Challenge for Finite Elements

HAL Id: hal-01803536 https://inria.hal.science/hal-01803536 Submitted on 30 May 2018 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Solving 2D linear isotropic elastodynamics by means of scalar potentials: a new challenge for finite elements Jorge Albella Martínez, Sébastien Imperiale, Patrick Joly, Jerónimo Rodríguez To cite this version: Jorge Albella Martínez, Sébastien Imperiale, Patrick Joly, Jerónimo Rodríguez. Solving 2D linear isotropic elastodynamics by means of scalar potentials: a new challenge for finite elements. Journal of Scientific Computing, 2018, �10.1007/s10915-018-0768-9�. �hal-01803536� Solving 2D linear isotropic elastodynamics by means of scalar potentials: a new challenge for finite elements Jorge Albella Mart´ınez1, S´ebastien Imperiale2,3, Patrick Joly2,4, and Jer´onimo Rodr´ıguez1 1Departamento de Matem´atica Aplicada, Universidade de Santiago de Compostela, 15706 Santiago de Compostela, Spain 2Inria, Universit´e Paris-Saclay, France 3LMS, Ecole Polytechnique, CNRS, Universit´e Paris-Saclay, France 4UMA, Ensta, CNRS, Universit´e Paris-Saclay, France 5IMAT, Universidade de Santiago de Compostela, 15706 Santiago de Compostela, Spain 6ITMATI, Campus Sur, 15706 Santiago de Compostela, Spain May 7, 2018 Abstract In this work we present a method for the computation of numerical solutions of 2D homogeneous isotropic elastodynamics equations by solving scalar wave equations. These equations act on the potentials of a Helmholtz decomposition of the displacement field and are decoupled inside the propagation domain. We detail how these equations are coupled at the boundary depending on the nature of the boundary condition satisfied by the displacement field. After presenting the case of rigid boundary conditions, that presents no specific difficulty, we tackle the challenging case of free surface boundary conditions that presents severe stability issues if a straightforward approach is used. We introduce an adequate functional framework as well as a time domain mixed formulation to circumvent these issues. Numerical results confirm the stability of the proposed approach. 1 Introduction In this paper our goal is to revisit a very classical question, namely, the numerical solution of elastodynamics equations in isotropic media, which govern the 1 propagation of elastic waves in solids, in the time domain. As a matter of fact there exist already many numerical methods for solving these equations. To begin with, for instance, in the framework of finite elements in space and finite differences in time, standard conforming finite elements (possibly high order) methods for the pure displacement formulation of the elastodynamics system (a second order hyperbolic system) as well as mixed finite element methods for the equivalent velocity-stress formulation of the same system (first order differential system). This space discretization is then coupled to explicit finite difference time stepping that is subject to a CFL stability condition. On the other hand, in many classical physics text books authors used the well-known Helmholtz decomposition of vector fields (write a vector field as the sum of a gradient and a curl) to compute analytical solutions in homogeneous isotropic media. Such a decomposition relates elastodynamic equations to two wave equations and enlightens the decomposition of the wave field as the sum of pressure waves (P-waves, that are gradients of a pressure potential ϕP) and shear waves (S-waves, that are curls of a shear potential ϕS) that propagate independently with different velocities, the velocity VPof the P-waves being larger than VSthe velocity of the S-waves. In the 2D case, to which we will restrict ourselves for simplicity, the simplification is that both pressure and shear potentials are scalar. However the extension to 3D does not pose a priori any additional conceptual difficulty and will be the object of further developments. In a piecewise homogeneous media, such a decomposition is valid locally and the different types of waves recouple at boundaries and interfaces. This is the main source of complexity of the propagation process. Looking at the literature, it seems that very few works have been devoted on the exploitation of this idea for finite element computations (however, one can find a few references concerning finite differences computations, see [1] in which a finite difference scheme is constructed with approximation properties independent of the ratio VP/VS) although it has been used in other domains of physics, in particular in fluid mechanics (current-vorticity formulations [2], chapter 2, [3] and [4]). The first motivation of the present work is an intellectual curiosity: could we use potentials to solve isotropic elastodynamics with finite elements? There is also a more relevant motivation concerning applications. This would concern the propagation of elastic waves in nearly incompressible media, such as soft tissues, in which P-waves propagate much faster than S-waves. In such a case, it is well known that displacement-based methods, which do not distinguish both waves along the calculation process, are greatly penalized by large values of the ratio VP/VSdue to the CFL condition (assuming that explicit time integrators are used). Let us explain this by a simple computation. Let us consider a ddimensional isotropic homogeneous medium of characteristic length Lon each direction, subject to a source term involving a minimal time scale T?. This source generates two different minimal wavelenghts, the S-wavelength λS= VS ?, which is much smaller than the P-wavelength λP=VP ? if VP/VSis large. Assuming that we consider P1or Q1finite elements on a quasi-regular mesh of step size h, for accuracy reasons, hshould be chosen proportional to λS, that is, h∝λS. On the other hand, considering a leap-frog time discretization for 2 instance, the time step is constrained by the stability condition that involves the fastest velocity VP, that is, ∆t∝h/VP. Considering a time interval integration [0, T], the number of time steps is T/∆tand since, one uses an explicit scheme, the cost of each iteration is proportional to number of degrees of freedom, namely Ld/hd. As a consequence, we can roughly estimate the overall computational cost as Cost ∝LdT hdT?∝Ld λd S T T? VP VS ,(1) where Ld/λd Srepresents the size of the problem in space and T/T?the size of the simulation in time. Clearly the penalizing factor is VP/VSwhich would not appear when solving a standard scalar wave equation. Potential formulations will a priori authorize the use of different meshes for both the pressure and shear potentials in view of adapting the mesh size to each wave length which is smaller for Sthan for Pwaves which would result in an important saving of the computation cost for large ratio VP/VS. This will be explained with more details in Section 2.3. A by product of this approach is that one could benefit of well-known techniques for the numerical treatment of the standard scalar wave equation such as for example the use of perfectly match layers (PMLs) for the treatment of unbounded domains. Indeed long time stable implementations of PMLs for isotropic elastodynamic equations raise some difficulties especially if the ratio VP/VSis large (even though the so-called C-PML solve the long time stability issue as shown in [5]). As the reader can expect, the main source of difficulties is the treatment of boundaries and interfaces because, contrary to the interior equations, boundary and transmission conditions are not easily expressed in terms of these potentials. The main goal and the main challenge we wish to address in this paper is the treatment of the couplings induced by these various surface conditions, that we need to handle in a guaranteed stable way and possibly with hopefully no influence on the CFL condition after time discretization. In [6, 7] we first addressed the case of a homogeneous medium with a clamped boundary, that is to say the Dirichlet boundary condition, for which we succeeded in achieving this goal. The approach and results of this paper will be recalled in Section 2. This second paper follows the philosophy of the previous work and aims at treating the free surface boundary condition, or Neumann boundary condition, that appears as much more challenging. This work is also preparatory to the treatment of interface conditions: there are two of them, one is of Dirichlet type, the second one of Neumann type and we anticipate that their treatment would rely on both treatments of Dirichlet and Neumann boundary conditions. The outline of the rest of the article is as follows. In Section 2, we first shall recap how to reduce the solution of 2D isotropic elastodynamics equations to two scalar wave equations and more importantly, explain how to treat the Dirichlet boundary condition as it has been done in [6]. The main section of this paper is Section 3 where we treat the free surface boundary condition. In Sect. 3.2 we show that the most naive approach directly inspired from the treatment of the Dirichlet condition gives rise to serious numerical stability problems af- 3 ter discretization. This is linked to the proposed variational formulation of the continuous model, the apparently natural functional space being too large and authorizing the development of unstable surface modes after space discretization. More precisely it appears that the new mass bilinear form which contains an additional boundary term (this is the main difference between Dirichlet and Neumann problems) fails to be positive contrary to the stiffness matrix (which remains the same as for the Dirichlet condition). So the key idea for the circumventing the problem is to find a smaller (but still sufficiently large) variational space in which we recover the positivity of the mass bilinear form. This is precisely the object of Section 3.3. For the construction of this space, we are guided by the comparison of the energy naturally associated to the new potentials formulation with the classical elastic energy associated to the displacement formulation. The definition of this new space is quite implicit and involves the solution of some elastostatic problem. Fortunately, such a space can be characterized as the orthogonal (with respect to the new mass bilinear form) of another subspace which is itself isomorphic to a space of scalar functions defined on the boundary. We can exploit this characterization by proposing a mixed variational formulation in which the above-mentioned orthogonality relation is treated as a constraint leading to the introduction of Lagrange multipliers as functions defined along the boundary. The resulting formulation is proven to be stable: this is the major achievement of this paper. Finally we show some numerical experiments that confirm the theoretical results previously obtained. 2 Decomposition into potentials: the case of a Dirichlet boundary condition This section has been added for pedagogical purpose, for making the paper self contained and for preparing Section 3. Section 2.1 recaps very standard material while Sections 2.2 and 2.3 are a summary to what has been done in [6]. 2.1 Decomposition into potentials in homogeneous media Preliminary notation. Throughout the paper we will work in 2D and x= (x1, x2) will denote the space variable. We shall use bold letters for representing vector fields such as u= (u1, u2) for the displacement field in a elastic body or v= (v1, v2) for the velocity field (v=∂tu). Ordinary letters will be used for scalar fields such as the components of the vector fields or the forthcoming potentials to be introduced. Finally, underlined bold letters will be used for 2×2 tensor fields such as the deformation or strain tensor ε(u) = εij(u) where 1 ≤i, j ≤2 or the stress tensor σ=σijthat represents the internal efforts inside the body. 4 Mathematical model. Let us briefly recap the 2D elastodynamics equations. First, the time variation of the displacement field uis governed by the fundamental law in mechanics ρ ∂2 tu−div σ=f,(2) where div σis the vector field defined by div σi=∂jσij(u), (with Einstein’s convention for summation over the repeated indices) ρ=ρ(x)≥ ρ0>0 is the density of the body that might depend on the xvariable for heterogeneous media and the source term f∈L1 loc(R+,(L2(Ω))2). Equation (2) must be completed by constitutive laws that relates the displacement field to the stress tensor. In an isotropic medium this is given by Hooke’s law which involves the (non negative) Lam´e parameters λ(x) and µ(x) σ=σ(u) := λdiv u I + 2 µε(u),(3) where Iis the 2 ×2 identity matrix and (we use again Einstein’s convention) div u=∂juj, εij(u) = 1 2∂iuj+∂jui,1≤i, j ≤2. One can eliminate the unknown σ(u) by substituting (3) in (2) and obtain a second order system in u. In the homogeneous case, i.e. when λ,µand ρare constant, we easily compute that div (σ(u)) = (λ+ 2µ)∇div u−µcurl curl u,(4) so that the equations can be written as follows (see [8, 9] for instance) ρ ∂2 tu−(λ+ 2µ)∇div u+µcurl curl u=f,(5) where we have introduced the two curl operators in 2D defined by curl u:=∂1u2−∂2u1,for the scalar curl of a vector field u, curl ϕ:=∂2ϕ, −∂1ϕ,for the vector curl of a scalar field ϕ. (6) Equation (5) is completed, in the presence of boundaries, with boundary conditions (see later) and, for the sake of simplicity, vanishing initial conditions u(t= 0) = 0, ∂tu(t= 0) = 0.(7) Decomposition into potentials. We are now going to introduce two scalar potentials ϕPand ϕSthat realize a Helmholtz decomposition, when there is no source term, of the velocity field v=∂tu, ρ ∂tϕP= (λ+ 2µ) div u, ρ ∂tϕS=−µcurl u.(8) 5 Substituting (8) into (5) leads, to ∂tv−∇ϕP−curl ϕS−g= 0 where g(t) = 1 ρZt 0 f(s)ds. (9) Imposing vanishing initial conditions for the potentials ϕP(t= 0) = 0, ϕS(t= 0) = 0,(10) we get (since v(t= 0) = 0) v=∇ϕP+curl ϕS+g,(11) which provides a Helmholtz decomposition [10] of the vector field vwhen gvanishes. To obtain the equations satisfied by the potentials we simply substitute (11) into the two equations in (8) differentiated in time to get two scalar wave equations for ϕPand ϕS 1 V2 P ∂2 tϕP−∆ϕP= div g,1 V2 S ∂2 tϕS−∆ϕS=−curl g,(12) with VP(resp. VS) the velocity of the P-waves (resp. S-waves) defined by VP=sλ+ 2µ ρ,P waves velocity, VS=rµ ρ,S waves velocity.(13) From (8) we obtain the initial conditions for the time derivative of the potentials ∂tϕP(t= 0) = 0, ∂tϕS(t= 0) = 0.(14) Note that in the free space (in absence of any boundary), the two wave equations in (12) are fully decoupled. 2.2 Decomposition into potentials for a clamped domain We now consider a 2D homogeneous isotropic propagation domain Ω R2, for instance, Ω bounded, with boundary Γ = ∂Ω that we assume to be clamped which means that equations (5) are completed with the boundary condition v=0,in Γ.(15) Proceeding as in the previous section, we introduce ϕPand ϕSvia equations (8) so that, inside Ω, the Helmholtz decomposition (11) holds and the potentials satisfy the scalar wave equations in (12). These equations must be completed by boundary conditions traducing (15). We assume in the following that Γ is a finite union of piecewise C1closed curves and thus admits almost everywhere a unit normal outward vector nand unit tangent vector τin such a way that the frame (τ,n) is a direct frame so that, if n= (n1, n2), then τ= (n2,−n1). For 6 any sufficiently smooth scalar field ϕwe have the following identities for traces in Γ curl ϕ·n=−∂τϕ, curl ϕ·τ=∂nϕ, (16) where as usual ∂nϕ=∇ϕ·nand ∂τϕ=∇ϕ·τ. Thus, writing that v=0is equivalent to writing v·n= 0 and v·τ= 0, which leads, according to (11), to the following boundary conditions for ϕPand ϕS ∂nϕP=∂τϕS−g·n, ∂nϕS=−∂τϕP−g·τ.(17) Note that the two essential conditions (15) for the displacement formulation become two natural conditions that couple the two potentials ϕPand ϕS. 2.3 A numerical approach for the Dirichlet problem: A recap Variational formulation. We first recall how to establish a weak formulation for the boundary value problem (12, 17). Assuming that the solution is sufficiently smooth, we can multiply the equations (12) by test functions ψP and ψSin H1(Ω), integrate by parts and use (17) to replace the normal derivatives of the potentials by tangential derivatives. After summation of the two resulting equations we can propose a first abstract variational formulation of the problem. To do so we introduce ϕ= (ϕP, ϕS) and ψ= (ψP, ψS) to get    Find ϕ(t) : R+−→ H1(Ω)2such that (ϕ, ∂tϕ)(t= 0) = (0,0) and d2 dt2mΩ(ϕ(t),ψ) + a(ϕ(t),ψ) = l(t, ψ),∀ψ∈H1(Ω)2, (18) where the linear form l(t, ·) is given by l(t, ψ) = −ZΩ g·∇ψP+curl ψSdx,(19) and the mass bilinear form mΩ(·,·) decouples ϕPand ϕS      mΩ(ϕ,ψ) = mP(ϕP, ψP) + mS(ϕS, ψS) mQ(ϕQ, ψQ) = 1 V2 QZΩ ϕQψQdx, Q ∈ {P, S}.(20) The stiffness bilinear form a(·,·) is given by a(ϕ,ψ) = aΩ(ϕ,ψ) + aΓ(ϕ,ψ),(21) where the volumic bilinear form aΩ(·,·) decouples in the same way as mΩ(·,·)      aΩ(ϕ,ψ) = aP(ϕP, ψP) + aS(ϕS, ψS) aQ(ϕQ, ψQ) = ZΩ∇ϕQ·∇ψQdx, Q ∈ {P, S},(22) 7 and the coupling surface bilinear form aΓ(·,·) is defined by aΓ(ϕ,ψ) = ZΓ∂τϕPψS−∂τϕSψPdγ,(23) where the integrals in the boundary should be interpreted as duality products between elements in H1 2(Γ) and its dual H−1 2(Γ). All the above bilinear forms are symmetric (for aΓ(·,·) use integration by parts along the boundary), however in order that (18) fits the classical theory of second order partial differential equations [11], some adequate positivity / coercivity for the forms mΩ(·,·) and a(·,·) need to be checked. The positivity of mΩ(·,·) is clear but the positivity of the a(·,·) is not obvious from (21) but relies on the following lemma Lemma 2.1. One has the identity: for any ϕ,ψ∈H1(Ω)2×H1(Ω)2 a(ϕ,ψ) = ZΩ∇ϕP+curl ϕS·∇ψP+curl ψSdx.(24) Proof. Let us denote ˜a(ϕ,ψ) the right hand side of (24). We obtain after expansion, using curl ϕS·curl ψS=∇ϕS·∇ψSand (22) ˜a(ϕ,ψ) = aΩ(ϕ,ψ) + ZΩ∇ϕP·curl ψSdx+ZΩ∇ψP·curl ϕSdx. Next we observe that (Green’s formula and div curl = 0) ZΩ∇ψP·curl ϕSdx=ZΓ ψPcurl ϕS·ndγ=−ZΓ ψP∂τϕSdγ.(25) In the same way ZΩ∇ϕP·curl ψSdx=−ZΓ ϕP∂τψSdγ=ZΓ ∂τϕPψSdγ,(26) after integration by parts along the boundary. To conclude we add (25) and (26) to infer that by definition of aΓ(·,·) (see (23)) ˜a(ϕ,ψ) = aΩ(ϕ,ψ) + aΓ(ϕ,ψ) = a(ϕ,ψ), according to the definition (21).  This lemma proves that the bilinear form a(·,·) is positive, but also suggests that H1(Ω)2is not the appropriate variational space for the weak formulation because of the coercivity requirement. That is why we introduce the space V:= ϕ= (ϕP, ϕS)∈L2(Ω)2such that ∇ϕP+curl ϕS∈L2(Ω)2.(27) Interpreting ϕas a vector field whose first component is ϕPand its second component is ϕS, one notices that ∇ϕP+curl ϕS= div ϕ −curl ϕ!(28) 8 3 The case of a free surface boundary condition In this section, we consider the case of a free surface boundary condition (the corresponding problem will be referred to as the Neumann problem in the sequel), that is, written in terms of the unknow v σ(v)n≡λdiv v n + 2 µε(v)n=0,on Γ.(51) In Section 3.1 we will recap some preliminary results on the Neumann problem. As it will be shown in Section 3.2, a naive extension of the technique explained in the previous section leads to an unstable variational formulation. The functional space in which it is set appears to be too large. To overcome this problem, in Section 3.3 this space is constrained in such a way that the new formulation is stabilized and suitable for finite element approximations. For the sake of simplicity (the reader will easily convince himself that this is not restrictive), we shall assume that Ω is bounded and simply connected, thus that Γ is a closed curve (see Figure 1). We shall also assume that Γ is parameterized by x(s)∈W1,∞(0, L) where sis the curvilinear abscissa along Γ and Lis the total length of Γ. Finally, we shall use the following notation for denoting a particular primitive of a function defined on Γ (where we arbitrarily particularize the point associated to s= 0 but this choice has no influence) ∀η∈L2(Γ),Iη(s) := Zs 0 η(σ) dσ∈H1(Γ).(52) It is clear that Ican be extended as a linear continuous operator I ∈ L(H−1/2(Γ), H1/2(Γ)).(53) 3.1 Preliminary recaps on the Neumann problem In this section we are interested in the problem (ρ ∂2 tu−div σ(u) = f,in Ω, σ(u)n=0,on Γ, (54) completed with the initial conditions (7). In this problem, a particular role is played by the 3 dimensional space of the so-called rigid displacements R(Ω) = nwR∈L2(Ω)2/ε(wR) = 0o =na(x2,−x1)t+ (b1, b2)t,(a, (b1, b2)) ∈R×R2o. (55) We introduce the spaces L2 R(Ω) = nw∈L2(Ω)2/ZΩ w·wRdx= 0,∀wR∈R(Ω)o, H1 R(Ω) = H1(Ω)2∩L2 R(Ω). (56) 15 so that we have the direct sums L2(Ω)2=L2 R(Ω) ⊕R(Ω), H1(Ω)2=H1 R(Ω) ⊕R(Ω),(57) which are orthogonal in L2(Ω)2. A classical but important property of problem (54) is provided in the following lemma. Lemma 3.1. If f(·, t)∈L2 R(Ω),∀t≥0, then ∀t≥0,u(·, t)∈H1 R(Ω).(58) Proof. Multiply (54) by wRand integrate over Ω. Using Green’s formula, d2 dt2ρZΩ u(·, t)·wRdx= 0 ∀wR∈R(Ω).(59) One concludes using the initial conditions.  In the sequel of this section we restrict ourselves to source terms satisfying f(·, t)∈L2 R(Ω),∀t≥0.(60) This is not restrictive due to the following remark. Remark 3.2. For a general source term f=fR+f⊥ R, with fR(·, t)∈L2 R(Ω) and f⊥ R(·, t)∈R(Ω), it is easy to see that the solution uof (54) can be decomposed as u=uR+u⊥ R, where uRis the solution of (54) with source term given by fRand u⊥ Ris given by u⊥ R=1 ρZt 0 (t−s)f⊥ R(·, s) ds. Another important property is Korn’s inequality in H1 R(Ω) (see [18]): Proposition 3.3. There exists a constant CΩ>0such that ∀w∈H1 R(Ω),kwk2 H1(Ω) ≤ CΩZΩ|ε(w)|2dx.(61) 3.2 The naive approach. Stability issues 3.2.1 The free boundary condition with potentials Of course, the first step consists in rewriting the free boundary condition (51) in terms of the potentials defined by (8) and (10) (as we did in section 2.2 by transforming the Dirichlet condition (15) into (17)). However, from (11), we see that, a priori, the condition (51) leads to an equation involving the second order space derivatives of the potentials, thus not well adapted for a finite element formulation. To overcome this, let us assume for a while that we know the value vΓof the velocity field on the boundary, or, in other words, that we want 16 to treat a non homogeneous Dirichlet boundary condition with v=vΓon the boundary. Then, proceeding as in section 2.2 for obtaining (17), we would use the non homogeneous boundary condition to obtain ∂nϕP=∂τϕS−g·n+vΓ·n, ∂nϕS=−∂τϕP−g·τ+vΓ·τ.(62) Then, since we do not know vΓ, we would like to compute it as a function of ϕ using the free boundary condition (51). To do so, we first remark that ε(v) = 2 div v I −curl v J +H(v),(63) where J= 0 1 −1 0 !,H(v) = −∂2v2∂1v2 ∂2v1−∂1v1!,(64) so that the free boundary condition (51) can be rewritten λ+ 2µdiv v n −µcurl v τ + 2 µH(v)n= 0,on Γ.(65) We deduce from (8), differentiating in time, that ρ ∂2 tϕP= (λ+ 2µ) div v, ρ ∂2 tϕS=−µcurl v,in Ω,(66) so that, assuming sufficient smoothness, we have on the boundary λ+ 2µdiv v n −µcurl v τ =ρ ∂2 tϕPn+ρ ∂2 tϕSτ,on Γ. After projection on coordinate axes, using τ= (n2,−n1), this can be rewritten λ+ 2µdiv v n −µcurl v n =ρ ∂2 tϕPn1+∂2 tϕSn2 ∂2 tϕPτ1+∂2 tϕSτ2 ≡ρ ∂2 tϕ·n ∂2 tϕ·τ!. In the same way, one computes that H(v)n|Γ= ∂τv2|Γ −∂τv1|Γ!≡ ∂τvΓ,2 −∂τvΓ,1 on Γ. Thus, recalling that µ=ρ V 2 S, the boundary condition (65) can be rewritten as ∂τvΓ,1=1 2V2 S ∂2 tϕ·τ, ∂τvΓ,2=−1 2V2 S ∂2 tϕ·n,on Γ.(67) This allows us, as desired, to compute vΓin terms of (ϕP, ϕS), up to an additive constant, using the operator I(see (52)). More precisely, P0(Γ) denoting the space of constant functions on Γ, using the fact that Γ is a closed curve and the initial conditions (14), equations (67) are easily seen to be equivalent to          vΓ,2+1 2V2 SI∂2 tϕ·n∈P0(Γ),ZΓ ϕ·n= 0, vΓ,1−1 2V2 SI∂2 tϕ·τ∈P0(Γ),ZΓ ϕ·τ= 0. (68) 17 In the sequel, the last two columns of (68) will be referred to as gauge conditions. To summarize this section, we have shown that vsatisfies (51) ⇔there exists (vΓ,ϕ) such that (62,68) are satisfied.(69) The form (62, 68) of the free boundary condition is the one that is useful for establishing the variational formulation of the problem (see section 3.2.2). Remark 3.4. In more general situations where the boundary has Nc>1con- nected components (each of them being smooth enough and closed) we would introduce one operator such as the one in (52) per component. In consequence, the potentials should satisfy 2Ncgauge conditions similar to those of the last column of (68). 3.2.2 Naive variational formulation According to section 3.2.1 (and also section 2), the problem we want to solve is                                                                Find ϕ= (ϕP, ϕS):Ω×R+→R2,vΓ: Γ ×R+→R2/ 1 V2 P ∂2 tϕP−∆ϕP= div g,in Ω ×R+,(i) 1 V2 S ∂2 tϕS−∆ϕS=−curl g,in Ω ×R+,(ii) ∂nϕP=∂τϕS−g·n+vΓ·n,on Γ ×R+,(iii) ∂nϕS=−∂τϕP−g·τ+vΓ·τ,on Γ ×R+,(iv) vΓ,2+1 2V2 SI∂2 tϕ·n∈P0(Γ),on Γ ×R+,(v) vΓ,1−1 2V2 SI∂2 tϕ·τ∈P0(Γ),on Γ ×R+,(vi) ZΓ ϕ·n=ZΓ ϕ·τ= 0,(vii) (70) completed with the initial conditions ϕ(·,0) = 0, ∂tϕ(·,0) = 0.(71) We are going to provide a variational formulation of (70, 71) which naturally eliminates vΓand provides a problem in ϕonly. We first take into account the last equation (70)(vii) by seeking ϕ(·, t) in V0, where V0is defined by V0:=nϕ∈Vs.t. ZΓ ϕ·n=ZΓ ϕ·τ= 0o.(72) 18 Next, we consider test functions ψ= (ψP, ψS) in V0and proceed as in Section 2.3 to obtain a variational formulation of the problem. The main difference comes from the additional boundary terms in the integration by parts due to the presence of vΓin the equations (70)(iii) and (iv). We obtain d2 dt2mΩ(ϕ(t),ψ)−ZΓvΓ·nψP+vΓ·τψSdσ+a(ϕ(t),ψ) = l(t, ψ).(73) Next, we apply the following identity (tricky but straightforward, its verification is left to the reader) vΓ·nψP+vΓ·τψS=vΓ,1ψ·n+vΓ,2ψ·τ,(74) so that, by using (70)(v) and (vi), we can eliminate vΓthanks to the fact that ψ∈V0. More precisely −ZΓvΓ·nψP+vΓ·τψSdσ=mΓ∂2 tϕ,ψ,(75) where we define the bilinear form mΓ(ϕ,ψ) := 1 2V2 SZΓI(ϕ·n)ψ·τ−I(ϕ·τ)ψ·ndγ.(76) Finally substituting (75) into (73), we see that ϕis solution of the following variational problem    Find ϕ(t) : R+−→ V0, satisfying (71) and such that d2 dt2m(ϕ(t),ψ) + a(ϕ(t),ψ) = l(t, ψ),∀ψ∈V0, (77) where the new mass bilinear form m(·,·) is defined by m(ϕ,ψ) = mΩ(ϕ,ψ) + mΓ(ϕ,ψ).(78) 3.2.3 Well-posedness issues At a first glance, the variational problem (77) looks like a nice hyperbolic variational problem in the sense of theory of Lions-Magenes [11]. We already saw that the bilinear form a(·,·) is continuous and coercive in V. Another good point is that the bilinear form mΓ(·,·) is symmetric (so m(·,·) is too) due to the observation that, using integration by parts along Γ, we can write, for all (ϕ,ψ)∈V0×V0 mΓ(ϕ,ψ) = 1 2V2 SZΓI(ϕ·n)ψ·τ+I(ψ·n)ϕ·τdγ.(79) In addition one observes that mΓ(·,·) is continuous in V0because the operator Imaps continuously H−1/2(Γ) into H1/2(Γ). Another nice property of m(·,·) is given by the following lemma. 19 Lemma 3.5. We have the injectivity result (i)m(ϕ,ψ)=0 ∀ψ∈V0=⇒(ii)ϕ= 0.(80) Proof. Let ϕsatisfying (i), since D(Ω)2⊂V0, in particular ∀ψ∈ D(Ω)2, m(ϕ,ψ) = mΩ(ϕ,ψ) = 0 thus ϕ= 0 by density of D(Ω)2in L2(Ω)2. However, all these properties are not sufficient to fit Lions-Magenes theory which also requires the positivity of m(·,·). Unfortunately, this fails to be true: Theorem 3.6. Assume that there is a part of the boundary Γthat is of class C2. Then, therere exists ψ∈V0such that m(ψ,ψ)<0.(81) Proof. Assume that the function x(s) that parametrizes Γ (see the begining of section 3) satisfies x(s)∈C2(a, b),for some [a, b]⊂[0, L]. Let n(s) be the unit normal vector to Γ at point x(s), outgoing with respect to Ω, and c(s) be the curvature of Γ at this point. Let us define ν+such that ν+=1 2sup s∈(a,b) 1 |c(s)|. By elementary differential geometry it is well known that the map (s, ν)∈(a, b)×(0, ν+)→x(s)−νn(s)∈R2 is injective. Moreover, there exists 0 < ν∗≤ν+such that Ω∗ a,b := {x(s)−νn(s), s ∈(a, b), ν ∈(0, ν∗)} ⊂ Ω and that (s, ν)→x(s) + νn(s) defines a change of variable from (a, b)×(0, ν∗) into Ω∗ a,b with jacobian J(s, ν) := 1 −ν c(s) that is uniformly bounded by J∗on (a, b)×(0, ν∗). Let θ∈ D(a, b) and θ6= 0 such that Zb a θ(s)ds = 0 and χ∈C∞(R+) such that supp χ⊂[0,1] and χ(0) = 1. Let 0 < δ < ν∗be a small parameter devoted to tend to 0, we define ψδ∈C1(Ω)2as ψδ(x) = (θ(s)τ(s)−θ0(s)n(s)χ(ν/δ),if x=x(s)−νn(s)∈Ω∗ a,b, 0,else, 20 Thanks to the assumption on θ,ψδbelongs to V0. On the one hand mΩ(ψδ,ψδ)≤δ(J∗/ V 2 S)kχk2 L2(R+)kθk2 H1(a,b).(82) On the other hand, along Γ, ψδ·τ=θ(s) for s∈(a, b) and 0 otherwise. In the same way ψδ·n=−θ0(s) for s∈(a, b) and 0 otherwise, so that I(ψδ·n) = −θ(s) for s∈(a, b) and 0 otherwise. Therefore mΓ(ψδ,ψδ) = −kθk2 L2(a,b)/V 2 S.(83) Comparing (82) and (83) it is clear that m(ψδ,ψδ) is strictly negative for δ small enough.  Remark 3.7. The technical assumption on Theorem 3.6 about the local regularity of Γis most likely unnecessary (and at the same time not very restrictive in practice) but needed for the proof above. Figure 1: Left: Definition of the parametrization of the boundary. Right: Definition of the curvilinear coordinates and notations for the construction of the function ψδin the proof of Theorem 3.6. As one can expect, the property of Theorem 3.6 is the cause of severe instabilities for any finite element approximation of the variational formulation (77). This will be put in evidence in the next section. 3.2.4 Numerical instabilities of Galerkin discretizations Introduction. We introduce a finite dimensional approximation of the space V0that will be denoted by V0hthat is constructed as follows V0h=Vh∩V0(84) where Vh=VP,h ×VS,h is a Galerkin approximation of the space H1(Ω)2. The problem to solve is    Find ϕh(t) : R+−→ V0hsuch that (ϕh, ∂tϕh)(t= 0) = (0,0) and d2 dt2mh(ϕh(t),ψh) + a(ϕh(t),ψh) = l(t, ψh),∀ψh∈V0h, (85) 21 where mh(·,·) is an approximation of the bilinear form m(·,·) (for instance it can be computed using quadrature formulae). At the algebraic level the formulation takes the following form Mh d2Φh dt2+AhΦh=Fh,Mh:= MΩ h+MΓ h,Ah:= AΩ h+AΓ h(86) where as in (43), Φt h= (ΦP,h,ΦS,h)tis the vectors of degrees of freedom of ϕP,h and ϕS,h. If one thinks of continuous Lagrange finite elements for instance, the matrix MΓ hhas the following structure MΓ h= MΓ,P hMΓ,P S h (MΓ,P S h)tMΓ,S h!, where the only non-zero values entries of the matrices (MΓ,P h,MΓ,S h,MΓ,P S h) correspond to the degrees of freedom located on the boundary Γ (note that due to the double integral on the boundary, each degree of freedom on the boundary is coupled with all the other degrees of freedom on the boundary). Despite of the injectivity property (80), it is not clear that the matrix Mhis invertible (which is a necessary property if one wants to use an explicit scheme in time). Indeed the proof was done at the continuous level using density properties of smooth compactly supported functions. At the discrete level, i.e. in finite dimensional space, such an argument can not be used any more. Moreover, even in the case where Mhis invertible, it is most likely that, because of the result of Theorem 3.6, the solution of the semi-discrete evolution problem, which is given by Φh(t) = Zt 0 (M−1 hAh)−1 2sin (M−1 hAh)1 2(t−s)M−1 hFh(s)ds, will blow up exponentially in time (i.e. the semi-discrete scheme (86) is unstable). Such a phenomenen is linked to the existence of strictly negative eigenvalues for the following symmetric eigenvalue problem (note that, as soon as Mh is invertible, such eigenvalues are necessarily real) Find Ψh6= 0 and λ∈Rsuch that AhΨh=λMhΨh.(87) Of course, the rate of exponential blow up will be, at least, given by the most negative eigenvalue λof M−1 hAh. More precisely, if σhdenotes the spectrum of M−1 hAhwith minimum value λ−(h)<0, then k(M−1 hAh)−1 2sin (M−1 hAh)1 2tk ≥ C et√|C(h)|,(88) where C(h)∈Ris such that 0 <C(h)<|λ−(h)|. In the following, we come back in more details on these invertibility and stability issues. In particular, in 22 the next paragraph we consider a simplified toy problem associated to a simple geometry. In this case, the two issues (invertibility and lack of positivity) can be studied analytically. Then, in the last paragraph of this section we consider a more general finite element discretization and illustrate the instability of the scheme through numerical computations. Remark 3.8. Note that because of the gauge condition (70)(vii) the two spaces are coupled at the boundary. In practice the gauge condition is imposed using a Lagrange multiplier as usually done when imposing zero average condition [19]. Study of a particular toy problem. We want to solve elastodynamic equations (2) on the cylinder Ω = (0,+∞)×(−π, π) (89) obtained when identifying the upper and lower boundaries, where x1= 0 is the free boundary. In consequence, we impose periodic boundary conditions u(t, x1, π) = u(t, x1,−π), ∂2u(t, x1, π) = ∂2u(t, x1,−π),∀x1∈R+,(90) in such a way that the boundary of Ω can be identified to a circle (this is a closed curve). Note that this particular example does not completely fit the assumptions made on the domain’s geometry since Ω is unbounded, however the reader will easily convince himself that this is not an essential issue. Next we discretize the space H1(Ω) as follows. We denote Vh⊂H1(R+) the uniform discretization of R+by P1finite elements of length h, i.e. Vh:= ψh∈H1(R+,C) s.t. ψh(x1) = +∞ X j=0 ψjwj(x1)(91) where the set {wj}are the piecewise affine functions that satisfy wj(kh) = δkj (where δis the Kronecker symbol) so that ψj=ψh(jh). Note that this space is isomorphic to the space `2(N). In the direction x2, we use a spectral method consisting in truncating the natural Fouirier series expansion of a function in L2(−π, π) at order L > 0, where 2π/L is thus the minimal oscillation length allowed in the approximate space. This corresponds to the following Galerkin approximation space Vh=VL h×VL hwith VL h:= nϕh∈H1(Ω,R) s.t. ϕh= L X `=−L ϕ` h(x1)e−i`x2with ϕ` h∈Vho.(92) Note that the approximation parameter is the couple (h, L) where his the space step in x1, devoted to tend to 0, and Lthe truncation parameter in frequency, devoted to tend to +∞(see also remark 3.9). It is then immediate to see that the space V0h(see (84)) is nothing but V0h=nϕh= (ϕP,h, ϕS,h)∈Vhs.t. ϕ0,0 P=ϕ0,0 S= 0o⊂V0. 23 where ϕ0,0 P:= ϕ0 P,h(0) and ϕ0,0 S:= ϕ0 S,h(0), according to the notation introduced above. Remark 3.9. If one makes an analogy with Q1finite elements for instance, this would corresponds to discretise Ωuniformly with rectangular finite elements of length h1=hin the direction x1and h2=π/2Lin the direction x2. Then we look for the solution of (85) with the following expression of the bilinear form a(·,·) and m(·,·) that take into account the fact that we deal with complex valued functions,                    mh(ϕh,ψh) = 1 V2 PZπ −πIR+ ϕP,hψP,h dx1dx2+1 V2 SZπ −πIR+ ϕS,hψS,h dx1dx2 −1 2V2 SZπ −πZx2 −π ϕP,h dsψS,h +Zx2 −π ψP,h dsϕS,h dx2, a(ϕh,ψh) = ZΩ∇ϕP,h +curl ϕS,h·∇ψP,h +curl ψS,h dx, where the symbol His used to account for the use of the following quadrature formula IR+ ϕhψhdx1:= h 2ϕ0ψ0+h +∞ X j=1 ϕjψj,kϕhk2 h:= IR+|ϕh|2dx1,(93) that allows us to get mass lumping. Since we have used a spectral approximation in x2, the previous problem decouples as a family of 2L+1 problems in 1Dwith respective unknowns {ϕ` h:= (ϕ` P,h, ϕ` S,h)}. More precisely, choosing ψh=e−i`x2ψ` h(x1) with ψ` h:= (ψ` P,h, ψ` S,h) as test functions in (85) and using the orthogonality properties of trigonometric functions, we get that, for each `∈ {−L, . . . , L},ϕ` h(t) : R+7→ Vh×Vhsatisfies the following 1D variational problem: for all ψ` h∈Vh×Vh, d2 dt2m`ϕ` h,ψ` h+a`ϕ` h,ψ` h=l`(t, ψ` h),(94) 24 The item (ii) will guide the construction of VNand is expected to guarantee the well-posedness of (108) and the stability of its Galerkin approximation. 3.3.1 Construction of the new variational space The energy naturally associated with (77), which is conserved as soon as the right hand side vanishes, is EN(t):=ED(t) + 1 2mΓ(∂tϕ(t), ∂tϕ(t)) =1 2[m(∂tϕ(t), ∂tϕ(t)) + a(ϕ(t),ϕ(t)) ]. (110) The positivity of this energy is obviously related to the positivity of the bilinear form m(·,·) on a space to which ϕbelongs. For the Dirichlet problem, the (obvious) positivity of the energy ED(t) was confirmed by the identities (i)ρ 2mΩ(∂tϕ, ∂tϕ) = Ep(t) and (ii)ρ 2a(ϕ,ϕ) = Ec(t),(111) where the potential energy Ep(·) (resp. kinetic energy Ec(·)) are defined by (37) (resp. (36)). We expect similar identities for the solution of the free boundary problem. In fact, it is clear that the identity (111)(ii) still holds for the solution of the Neumann problem because its proof does not refer to the Dirichlet condition: it only uses (11) which is valid independently of the boundary condition. We just have to obtain an equivalent of (111)(i) with m(·,·) instead of mΩ(·,·), i.e., ρ 2m(∂tϕ, ∂tϕ) = Ep(t).(112) This is related to the identity (38) in Lemma 2.3 which yields Ep(t) = ρ 2mΩ(∂tϕ, ∂tϕ)−2µZΓ u2∂τu1dγ,(113) so that we have to check that, since 2µ= 2V2 Sρ, mΓ(∂tϕ, ∂tϕ) = −4V2 SZΓ u2∂τu1dγ.(114) Here is where we are going to use the fact that usatisfies the free boundary condition. More precisely, integrating in time (67), we have ∂τu1=1 2V2 S ∂tϕ·τ, u2+1 2V2 SI∂tϕ·n∈P0(Γ).(115) Then, since ϕ(t) belongs to V0, (114) follows. The reader should notice that which was important for obtaining (112) is that we could write (see (8)) ∂tϕP=V2 Pdiv u, ∂tϕS=−V2 Scurl u, 31 with ua smooth enough function satisfying the free boundary condition (51). The above observation gives the idea of the construction of the space VN. Let us first introduce the space D:={w∈H1(Ω)2such that div σ(w)∈L2(Ω)2} ≡ {w∈H1(Ω)2/−V2 P∇(div w) + V2 Scurl (curl w)∈L2(Ω)2}, (116) where the second line comes from using (4). This space is a Hilbert space for the norm: kwk2 D=kwk2 H1(Ω) +kdiv σ(w)k2 L2(Ω).(117) Note that, by construction of Dand definition of V, ∀u∈D,Fu:= V2 Pdiv u,−V2 Scurl u∈V.(118) Next, we consider the closed subspace of Dof vector fields that are orthogonal to the rigid displacements and satisfy the free boundary condition DN:={w∈D∩H1 R(Ω) such that σ(w)n= 0 on Γ},(119) and finally the space VN:= F(DN)≡nV2 Pdiv w,−V2 Scurl w,w∈DNo.(120) A first remarkable property of this space is given by the lemma Lemma 3.12. The space VNis a subspace of V0. Proof. By (118) we already know that VN⊂V. So we simply have to check the gauge conditions in (72). Let ψ∈VN, i.e., ψ= (V2 Pdiv w,−V2 Scurl w),w∈DN.(121) The proof is essentially a matter of reproducing the computation in section 3.2 for proving (67), with ψinstead of ∂2 tϕand winstead of v. Simply note that σ(w)n=0on Γ follows from w∈DNwhile the equivalent of (66) is nothing but (121). Then we obtain the equivalent of (67), namely, ∂τw1=1 2V2 S ψ·τ, ∂τw2=−1 2V2 S ψ·n,on Γ.(122) Finally, the gauge conditions are simply obtained by integrating the above equalities on Γ.  We shall use later another nice property of the space VN: Lemma 3.13. For any ψ= (ψP, ψS)∈VN,∇ψP+curl ψS∈L2 R(Ω). 32 Proof. Let ψ∈VN. Then, ψ= (V2 Pdiv w,−V2 Scurl w) with w∈DN. By (4), div σ(w) = ρ∇ψP+curl ψS, that we multiply by any wR∈R(Ω) and integrate over Ω to obtain, using Green’s formula ZΩ σ(w) : ε(wR) dx+ZΓ σ(w)n·wRdγ=−ρZΩ∇ψP+curl ψS·wRdx= 0, since ε(wR) = 0 and σ(w)n= 0 for w∈DN. Now, we remark that, thanks to the property (58) (see Lemma 3.1) and Lemma 3.12, the vector of potentials ϕ(·, t) belongs to VNfor all t≥0. This implies that the space VNsatisfies the requirement (108). Next we prove that the space VNsatisfies the requirement (109). Theorem 3.14. The bilinear form m(·,·)is positive definite in the space VN. Furthermore, there exists C>0, only depending on Ω,λ,µand ρsuch that m(ψ,ψ)≥CZΩ|ψ|2dx,∀ψ∈VN. Proof. Let ψ= (ψP, ψS)∈VN, i.e., ψ= (V2 Pdiv w,−V2 Scurl w),w∈DN. Applying the identity (38) ZΩ|ε(w)|2dx=ZΩ|div w|2dx+1 2ZΩ|curl w|2dx−2ZΓ w2∂τw1dγ so that we deduce from (3) that ZΩ σ(w) : ε(w) dx= (λ+ 2µ)ZΩ|div w|2dx+µZΩ|curl w|2dx −4µZΓ w2∂τw1dγ. Since λ+ 2µ=ρ V 2 Pand µ=ρ V 2 S, we compute from the definition of ψthat (λ+ 2µ)ZΩ|div w|2dx+µZΩ|curl w|2dx=ρ mΩ(ψ,ψ). On the other hand, using (122) (see the proof of Lemma 3.12) we obtain −4µZΓ w2∂τw1dγ=ρ mΓ(ψ,ψ). which, combined to the two previous equalities, gives m(ψ,ψ) = 1 ρZΩ ε(w) : σ(w) dx. Therefore, since, by (3) again, σ(w) : ε(w) = λdivw2+ 2 µ|ε(w)|2≥2µ|ε(w)|2, 33 we get m(ψ,ψ)≥2V2 SZΩ|ε(w)|2dx. In the space DN⊂H1 R(Ω), we can use the Korn’s inequality (61) to obtain ZΩ|ε(w)|2dx≥2CΩV2 SZΩ|∇w|2dx. This allows us to conclude since, obviously, from the definition of ψ, ZΩ|∇w|2dx≥e CZΩ|ψ|2dx, with e Conly depending on VPand VS. At this level, we have identified a good space VNsatisfying the requirements at the beginning of this section. However, the definition of this space (120) is quite theoretical and implicit and thus rather hard to use it numerically: in particular, it refers to displacement fields, which we want precisely to avoid. The goal of the next section is to give a more suitable characterization of VN. 3.3.2 A characterization of the new space Let ΠRbe the L2(Ω)2orthogonal projection onto the space L2 R(Ω) ≡R(Ω)⊥. Given ψ∈V, let us denote w?:= SNψthe solution of the following elastostatic problem (whose well-posedness follows from Fredholm alternative)   −div σ(w?) = −ρΠR∇ψP+curl ψS,in Ω, w?∈H1 R(Ω),σ(w?)n=0,on Γ. (123) The reader will easily verify that, by construction, SN∈ L(V;D) and Im(SN)⊂DN. Even more we have the following result: Lemma 3.15. The operator SNis a left inverse of Frestricted to the space DN. More precisely, ∀w∈DN,w=SNFw. As a consequence, the image of SNcoincides with the space DN. Proof. Let w∈DN. As −div σ(w) = −ρ∇ψP+curl ψSwith ψ=Fw (cf. the proof of Lemma 3.13), we deduce from Lemma 3.13 that div σ(w)∈L2 R(Ω). Then div σ(w) coincides with its own projection onto L2 R(Ω) and we can write div σ(w) = ρΠR∇ψP+curl ψS.(124) 34 Moreover, since w∈DNwe have w∈H1 R(Ω),σ(w)n=0on Γ.(125) Finally, (125) and (124) prove that w=SNψ, i.e., w=SNFw. Of course, writing w=SNFwfor any w∈DN, proves that DN⊂ ImSN. Since F ∈ L(D,V) we can define T:= F ◦SN∈ L(V),and ImT=VN,(126) where the equality derives from the definition of VNand Lemma 3.15. We summarize in Figure 6 the images and preimages of the operator introduced so far, i.e. F,SNand T. Figure 6: Representation of the images and pre-images of the operator F,SN and T. Note that, by definition of SN(see (123)) and F(118) Tψ= ( V2 Pdiv w?,−V2 Scurl w?) where w?is the solution of (123).(127) A straighforward, but important, consequence of Lemma 3.15 is that Tis a projector into VN. Indeed, ∀ψ∈V,T2ψ=F SNF SNψ=FSNFSNψ.(128) Since SNψ∈DN, by Lemma 3.15, SNFSNψ=SNψand thus ∀ψ∈V,T2ψ=FSNψ=Tψ⇐⇒ T =T2.(129) 35 As for any projector, we can write the direct sum V= ker T ⊕ImT= ker T ⊕VN.(130) Contrary to Vor VN, it is possible to give an explicit description of the space ker Twhich is completely independent of the spaces Dor DN. More precisely, Lemma 3.16. The kernel of the operator Tis characterized by ker T={ψo= (ψo P, ψo S)∈V/∇ψo P+curl ψo S∈R(Ω)}.(131) Proof. Let ϕ∈ker T, and w?=SNϕ, solution of (123). Then Tϕ= 0 means that div w?= curl w?= 0, thus div σ(w?) = 0by (4). In consequence of (123), ΠR∇ϕP+ curl ϕS=0. The reverse implication is trivial.  The decomposition in (130) is not orthogonal in the classical sense but orthogonal with respect to m(·,·) in the sense of the following theorem. Theorem 3.17. We have VN=ker T⊥,m where ker T⊥,m := {ψ∈V/∀ψo∈ker T, m(ψ,ψo)=0}. Proof. Step 1:VN⊂ker T⊥,m. Let ϕ= (ϕP, ϕS)∈VNand ψo= (ψo P, ψo S)∈ker T. We know that, by definition of VN,ϕ=V2 Pdiv w,−V2 Scurl wwith w∈DN. Let us compute m(ϕ,ψ) = mΩ(ϕ,ψo) + mΓ(ϕ,ψo). We have, by Green’s formulas, 1 V2 PZΩ ϕPψo Pdx=ZΩ div wψo Pdx=−ZΩ w·∇ψo Pdx+ZΓ (w·n)ψo Pdγ, 1 V2 SZΩ ϕSψo Sdx=−ZΩ curl wψo Sdx=−ZΩ w·curl ψo Sdx+ZΓ (w·τ)ψo Sdγ. Let us add the two equalities. Since by Lemma 3.16, ∇ψo P+curl ψo S∈R(Ω), while w∈DN⊂L2 R(Ω), the volume integral vanishes and mΩ(ϕ,ψo) = ZΓ(w·n)ψo P+ (w·τ)ψo Sdγ. (132) According to the proof of Lemma 3.12, we can use (122), with ϕinstead of ψ, which gives, for come constants C1and C2, 1 2V2 SI(ϕ·τ) = w1+C1,−1 2V2 SI(ϕ·n) = w2+C2,on Γ.(133) Substituting (133) into the expression (76) of mΓ(·,·), we get, using the gauge conditions for ψo, mΓ(ϕ,ψo) = −ZΓw1(ψo·τ) + w2(ψo·n)dγ. (134) 36 Finally, adding (132) and (134) gives m(ϕ,ψo) = 0 thanks to the identity (74). Step 2:(ker T)⊥,m ⊂VN. Let ϕ∈(ker T)⊥,m, i.e. such that m(ϕ,ψo) = 0 for all ψo∈ker T. Let ψ∈V. Using (130), we decompose ϕ,ψas: ϕ=ϕo+ϕN,ψ=ψo+ψN,(ϕo,ψo)∈(ker T)2,(ϕN,ψN)∈V2 N. Our goal is to prove that ϕo= 0. Note that by step 1, m(ϕo,ψN) = m(ϕN,ψo)=0.(135) Therefore, we compute m(ϕo,ψ) = m(ϕo,ψo) + m(ϕo,ψN) = m(ϕo,ψo). Next, since ϕo=ϕ−ϕNand ϕ∈(ker T)⊥,m, m(ϕo,ψ) = m(ϕ,ψo)−m(ϕN,ψo) = (135) m(ϕ,ψo) = 0. Thus m(ϕo,ψ)=0,∀ψ∈V. In particular, for ψ∈(D(Ω))2⊂V, mΩ(ϕo,ψ) = m(ϕo,ψ) = 0. We conclude by density of D(Ω)2in L2(Ω)2that ϕo=0, thus ϕ∈VN. 3.3.3 A first stabilized mixed formulation In this section, we are going to exploit Theorem 3.17, namely ϕ∈VN⇐⇒ ϕ∈Vand m(ϕ,ψ)=0,∀ψ∈ker T, by reinterpreting the last condition as an equality constraint on ϕ. As it is usual for treating equality constraints (see [19, 20, 21]), we are going to introduce a Lagrange multiplier in the space ker T. This leads us to write the following mixed problem            Find (ϕ(t),ϕo(t)) : R+−→ V×ker Tsatisfying (71) and d2 dt2m(ϕ(t),ψ) + a(ϕ(t),ψ) + m(ϕo(t),ψ) = l(t, ψ),∀ψ∈V, m(ϕ(t),ψo) = 0,∀ψo∈ker T. (136) Then we have the following equivalence theorem between the above mixed problem and the variational problem (108) posed in the space VN. Theorem 3.18. The problem (136) admits a unique solution given by (ϕ(t),0) where ϕ(t)is the solution of the problem (108). 37 Proof. Let ϕ(t) be the solution of (108). First, it is clear that it satisfies the second equation in (136) since ϕ(t)∈VNand VN= (ker T)⊥,m (Theorem 3.17). This also implies that d2 dt2m(ϕ(t),ψo) = 0,∀ψo∈ker T. On the other hand, by definition (31) of a(·,·), we have a(ϕ(t),ψo) = ZΩ∇ϕP(t) + curl ϕS(t)·∇ψo P+curl ψo Sdx. By Lemma 3.16, we know that ∇ψo P+curl ψo S∈R(Ω). Then, since ϕ(t)∈VN we deduce a(ϕ(t),ψo) = 0 from Lemma 3.13. Finally, thanks to the assumption (60) on the source term, defined by (9, 19), we have l(t, ψo) = 0. Thus, we conclude that d2 dt2m(ϕ(t),ψo) + a(ϕ(t),ψo) = l(t, ψo),∀ψo∈ker T. Since ϕis the solution (108), we also have d2 dt2m(ϕ(t),ψ) + a(ϕ(t),ψ) = l(t, ψ),∀ψ∈VN, and then, by linearity and using the decomposition (130) of V, we obtain the above equality also for all ψ∈Vwhich is nothing but the first equation in (136) with ϕo=0. We thus have proven that (ϕ(t),0) is solution of (136). It remains to prove the uniqueness of solutions of (136). Let (ϕ,ϕo) be a solution of (136) with l(·,·) = 0 . Then, by restricting the test function ψin the second line of (136) to ψ∈VN, using again the m−orthogonality of VNand ker T, we deduce that ϕis the solution of (108) with l(·,·) = 0. Thus ϕ=0. We then deduce from the second line in (136) again that m(ϕo(t),ψ)=0,∀ψ∈V, which leads to ϕo=0due to the injectivity property (80) (see Lemma 3.5).  The reader which is not familiar with mixed variational formulation could be surprised that we introduce an additionnal unknown that we know is 0. All the interest of this new formulation is when Galerkin discretization is concerned: after discretization we get a stable problem and the discrete approximation of the unknown ϕois no longer 0(see Remark 3.23). 3.3.4 Characterization of the multipliers space ker T Even though the mixed variational problem (136) is nicer than (108) in the sense that any reference to displacement fields is removed, it is still not completely 38 satisfactory for finite element approximation since we need a priori to construct a Galerkin approximation space for ker T. In order to work around this problem we are going to characterize the space ker Tup to an explicit three dimensional subspace (see Lemma 3.19) as the image by an explicit mapping Eof a space Mof functions along the boundary Γ (see (143)). This is then satisfactory from the numerical point of view because it is easy to approximate the space M with finite elements (defined on the boundary Γ) while the operator Eis easy to approximate numerically. This will result in an alternative reformulation of the mixed problem which will be given in section 3.3.5. We first begin with a lemma. Lemma 3.19. The space ker Tcan be decomposed as the direct sum ker T=KR⊕K0,(137) where K0is the closed space of Vof the so-called harmonic fields defined by K0=nψ= (ψP, ψS)∈V/∇ψP+curl ψS=0o ≡nψ∈V/div ψ= curl ψ= 0o, (138) and KRis the 3dimensional space KR= spanϕ1, ϕ2, ϕ3where ϕ1=1 2(x1, x2)t,ϕ2=1 2(x2,−x1)t,ϕ3=1 2(0, x2 1+x2 2)t.(139) Proof. First, we easily check that div ϕ1 −curl ϕ1=1 0,div ϕ2 −curl ϕ2=0 1,div ϕ3 −curl ϕ3=x2 −x1.(140) By (55), if ψ∈ker T, there exist (a1, a2, b)∈Rsuch that div ψ −curl ψ=a11 0+a20 1+bx2 −x1. Then ψ−a1ϕ1+a2ϕ2+bϕ3∈K0. To conclude is suffices to remark, using the formulae (140), that KR∩K0=∅. Next, we show that the space K0can be identified with a space of functions defined on the boundary Γ. To do so, we introduce the classical normal trace map γ∈ LH(div,Ω), H−1 2(Γ), γ :ϕ−→ ϕ·n|Γ.(141) Theorem 3.20. The map γis an isomorphism from K0onto M:= nν∈H−1 2(Γ) /ZΓ νdγ= 0o. 39 Proof. The fact that γϕ∈Mfor ϕ∈K0follows from Green’s formula ZΓ ϕ·ndγ=ZΩ div ϕdγ= 0,since div ϕ= 0 (second line of (138)). To conclude, it suffices to show that for any ν∈M, there exists a unique ϕ∈K0such that γϕ=ν. For the existence, let pbe the unique solution of the Neumann problem (note that ν∈Myields the compatibility condition required for the unique solvability of this problem)      Find p∈H1(Ω)/Rsuch that −∆p= 0,in Ω ∂np=ν, in Γ, (142) Then setting ϕ=∇pwe have div ϕ= curl ϕ= 0 and ϕ·n=∂np=νon Γ, hence ϕ∈K0. For the uniqueness, we have simply to remark that curl ψ= 0 in Ω implies that ψ=∇q, with q∈H1(Ω)/R(see Theorem 2.9 in [2]). Then if ψ·n= 0 on Γ and div ψ= 0 we have ∆q= 0 as well as ∂nq= 0 therefore q= 0 and ψ=0.  The above proof shows that the inverse of the map γis the lifting operator E ∈ L(M, K0), where Eµ:= ∇p, with pthe unique solution of (142). In other words, Theorem 3.20 can be rephrased as K0=Eν / ν ∈M.(143) This characterization brings a new light on the absence of positivity of the bilinear form m(·,·). Corollary 3.21. For all ψ∈K0,m(ψ,ψ)≤0. Proof. From (143) we know that there exists νsuch that ψ=Eν=∇pwith pthe unique solution of (142), then m(∇p, ∇p) = 1 V2 PZΩ|∂1p|2dx+1 V2 SZΩ|∂2p|2dx−1 V2 SZΓI(∇p·τ)∇p·ndγ. Since I(∇p·τ) = p+cfor a constant cand since ∇p·n=νhas zero average along Γ ZΓI(∇p·τ)∇p·n=ZΩ|∇p|2dx+ZΩ p∆pdx=ZΩ|∇p|2dx. Therefore combining the two previous equation we obtain m(∇p, ∇p) = 1 V2 P−1 V2 SZΩ|∂1p|2dx≤0.  40 quasi-regular mesh of approximately 16000 triangles, the time step is ∆t= 0.01 and the time scheme used is the explicit scheme (148). We plot Figure 8 and 9 snapshots of the obtained solution. For comparison we also plot a snapshot of the velocity field e vhobtained by the standard P1finite element discretization of the elastodynamics equations (2). The results obtained are stable in time, even for long time of simulations, moreover the reconstructed velocity field defined by vh=∇ϕP,h +curl ϕS,h −gshow good agreements with the direct computations of the velocity field. −0.06 −0.04 −0.02 0 0.02 0.04 0.06 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 P,h(1,·) S,h(1,·) S,h(1.5,·) P,h(1.5,·) Figure 8: Snapshot of (ϕP,h, ϕS,h), solution of problem (148) for different time of simulation. 4 Conclusions and perspectives We have presented a method for the computation of solutions of an isotropic elastodynamics problem by solving scalar decoupled wave equations acting on the potentials of a Helmholtz decomposition of the displacement field. We detailed how these equations are coupled at the boundary and how this coupling 47 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 |P,h |(1.5,·) |curl S,h|(1.5,·) | vh|(1.5,·) |vh|(1.5,·) Figure 9: Snapshot of the reconstructed velocity field from the solution of problem (148) as well as the solution e vhcomputed by solving (2). takes into account the nature of the boundary conditions satisfied by the displacement field. Although the case of rigid boundary conditions presented no specific difficulty, the challenge appeared to deal with the free surface boundary conditions since severe stability issues was revealed. We solve these issues by introducing an adequate functional framework in which one has to look for the solutions to avoid these instabilities. A mixed formulation was constructed to ensure, using a Lagrange multiplier, that the sought solutions indeed belongs to the adequate functional space. Numerical results confirm the stability of the proposed approach. Among the perspectives of this work we can mention, first, the construction and analysis of an efficient numerical method, including the analysis of the discrete lifting operator introduced in the last section. Then, the analysis of the transmission problem between two isotropic media could be addressed, however we do not expect theoretical difficulties since one can see transmission conditions as both heterogeneous Dirichlet and Neumann bound- 48 ary conditions. Finally, the 3D case should be addressed, the difficulty being that the potential corresponding to the shear waves in the Helmholtz decomposition is no longer scalar and is associated with a gauge condition that should be taken into account. Acknowledgments The research of the first and fourth authors was partially funded by FEDER and the Spanish Ministry of Science and Innovation through grants MTM2013-43745-R and MTM2017-86459-R and by Xunta de Galicia through grant ED431C 2017/60. References [1] J. Virieux. P-SV wave propagation in heterogeneous media: Velocitystress finite-diference method. Geophysics, 51(4):889–901, 1986. [2] Girault V. and Raviart P.-A. Finite element methods for Navier-Stokes equations: theory and algorithms, volume 5. Springer Science & Business Media, 2012. [3] Glowinski R. and Pironneau O. Numerical methods for the first biharmonic equation and the two-dimensional Stokes problem. SIAM Rev., 21(2):167– 212, 1979. [4] Babuska I., Osborn J., and Pitk¨aranta J. Analysis of mixed methods using mesh dependent norms. Mathematics of Computation, 35(152):1039–1062, 1980. [5] Komatitsch D. and Martin R. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation. Geophysics, 72(5):SM155–SM167, 2007. [6] Burel A., Imperiale S., and Joly P. Solving the homogeneous isotropic linear elastodynamics equations using potentials and finite elements. The case of the rigid boundary condition. Numerical Analysis and Applications, 5(2):136–143, 2012. [7] Burel A. Contributions `a la simulation num´erique en ´elastodynamique: d´ecouplage des ondes P et S, mod`eles asymptotiques pour la travers´ee de couches minces. PhD thesis, Universit´e Paris Sud-Paris XI, 2014. [8] Gurtin M. E. An introduction to continuum mechanics, volume 158. Academic press, 1982. [9] Ciarlet P. G. Elasticit´e tridimensionnelle, volume 1. Masson, 1986. [10] Alonso Rodr´ıguez A. and Valli A. Eddy Current Approximation of Maxwell Equations: Theory, Algorithms and Applications, volume 4. Springer Science & Business Media, 2010. 49 [11] Lions J.-L. and E. Magenes. Non-Homogeneous Boundary Value Problems and Applications, volume 1. Springer Berlin Heidelberg, 1972. [12] Monk P. Finite Element Methods for Maxwell’s Equations. Oxford University Press, 2003. [13] Cherif A., Bernardi C., Dauge M., and Girault V. Vector potentials in threedimensional non-smooth domains. Mathematical Methods in the Applied Sciences, 21(9):823–864, 1998. [14] Ciarlet P. G. The finite element method for elliptic problems. SIAM, 2002. [15] Cohen G. Higher-Order Numerical Methods for Transient Wave Equations. Springer, 2002. [16] Komatitsch D. and Tromp J. Introduction to the spectral element method for three-dimensional seismic wave propagation. Geophysical Journal International, 139(3):806–822, 1999. [17] Cohen G., Joly P., Roberts J. E., and Tordjman N. Higher order triangular finite elements with mass lumping for the wave equation. SIAM Journal on Numerical Analysis, 38(6):2047–2078, 2001. [18] Ciarlet P. G. On Korn’s inequality. Chinese Annals of Mathematics-Series B, 31(5):607–618, 2010. [19] Pavel Bochev and Richard B Lehoucq. On the finite element solution of the pure Neumann problem. SIAM review, 47(1):50–66, 2005. [20] Alfredo Berm´udez de Castro, Dolores G´omez, and Pilar Salgado. Mathematical models and numerical simulation in electromagnetism, volume 74. Springer, 2014. [21] Wei Jiang, Na Liu, Yifa Tang, and Qing Huo Liu. Mixed finite element method for 2d vector Maxwell’s eigenvalue problem in anisotropic media. Progress In Electromagnetics Research, 148:159–170, 2014. [22] Franco Brezzi and Michel Fortin. Mixed and Hybrid Finite Element Methods. Springer-Verlag New York, Inc., New York, NY, USA, 1991. 50