A discontinuous Galerkin method for non-linear

0 downloads 0 Views 5MB Size Report
Lina Homsi · Ludovic Noels. Received: date / Accepted: date. Abstract A coupled .... for convection-dominated thermo-poro-mechanics. Furthermore, Zheng et al.
Noname manuscript No. (will be inserted by the editor)

A discontinuous Galerkin method for non-linear electro-thermo-mechanical problems Application to shape memory composite materials Lina Homsi · Ludovic Noels

Received: date / Accepted: date

Abstract A coupled Electro-Thermo-Mechanical Discontinuous Galerkin (DG) method is developed considering the non-linear interactions of electrical, thermal, and mechanical fields. In order to develop a stable discontinuous Galerkin formulation the governing equations are expressed in terms of energetically conjugated fields gradients and fluxes. Moreover, the DG method is formulated in finite deformations and finite fields variations. The multi-physics DG formulation is shown to satisfy the consistency condition, and the uniqueness and optimal convergence rate properties are derived under the assumption of small deformation. First the numerical properties are verified on a simple numerical example, and then the framework is applied to simulate the response of smart composite materials in which the shape memory effect of the matrix is triggered by the Joule effect. Keywords Discontinuous Galerkin Method · Electro-thermo-mechanics · non-linear elliptic problems · smart composites

1 Introduction The emergence of complex multi-physics materials, such as multi-functional and shape-memory materials, see the reviews [37, 36, 18, 51] as a non-exhaustive list, has motivated the development of reliable, accurate, and efficient multi-physics numerical models. As an example, Shape Memory Polymer (SMP) [9, 44] can fix a temporary deformed configuration and recover their initial shape upon application of a stimulus such as temperature [34], light [35], electric field [17], magnetic field [56], water [29], or solvent [41]. Temperature triggered SMP take advantage of a property change at the glass transition temperature Tg : Below Tg , the movement of the polymer segments is frozen and the polymers are considered to be in a glassy L. Homsi, L. Noels University of Li` ege, Department of Aerospace and Mechanical Engineering, Computational & Multiscale Mechanics of Materials (CM3), Quartier Polytech 1, All´ ee de la D´ ecouverte 9, B-4000 Li` ege, Belgium Tel.: +32 4 366 48 26 Fax: +32 4 366 95 05 E-mail: {Lina.Homsi,L.Noels}@ulg.ac.be

2

Lina Homsi, Ludovic Noels

state. Once they are heated above Tg the chains become weak and the polymers are considered to be in a rubbery state, such that the materials can be deformed with minimal force. However, SMP have the drawback of low strength and stiffness when they are used for structural applications. This drawback can be overcome by dispersing continuous or discontinuous reinforcements throughout a polymer matrix, leading to Shape Memory Polymer Composites (SMPC) [45, 14]. Besides improving the mechanical properties, the fillers can improve the shape memory recovery stress and, in addition, act as triggering mechanisms. For example, carbon fibers exhibit conductivity which can be exploited as a shape memory triggering mechanism since the increase of temperature required to trigger the Shape Memory effect of the matrix is obtained through Joule effect by applying an electric current. Henceforth SMPC are seeing a growing interest in the area of deployable structures, sensors, actuators etc However, several difficulties arise when modeling the response of this kind of materials. On the one hand, multi-physics numerical models involving strong coupling are required. On the other hand, as SMP are capable of large deformations (high recovery strain) and as the field variations are important, the numerical methods ought to be formulated in a finite deformation and fields variation setting. As a result the governing equations involve strong non-linear coupling. Muliana et al. [47] have studied the time dependent response of active piezoelectric fiber and polymer composite materials in a multi-scale approach. Rothe et al. [55] have considered Electro-Thermo-Elasticity in a small strain setting, where they have focused on the numerical treatment of the monolithic approach, with the development of a one dimensional analytical solution in the purpose of code verification. In particular Zhupanska et al. [74] have discussed the governing equations describing electromagnetic, thermal, and mechanical field interactions. Nevertheless these contributions are still limited to small deformation settings. In this work, a multi-field coupling resolution strategy is used for the resolution of electrical, energy, and momentum conservation equations by means of the Discontinuous Galerkin (DG) Finite Element Method (FEM) to solve ElectroThermo-Mechanical (ETM) coupling. The main idea of the DG formalism is to approximate the solution by piece-wise continuous polynomial functions, and to constrain weakly the compatibility between elements. The inter-elements weak enforcement of the continuity allows using discontinuous polynomial spaces of high degree and facilitates handling elements of different types and dynamic mesh modifications. Indeed, the possibility of using irregular and non conforming meshes in an algorithm makes it suitable for time dependent transient problems. It also allows having hanging nodes and different polynomial degrees at the interface, with a view to hp-adaptivity. These DG opportunities and their merits have been illustrated and discussed by Kaufmann et al. [32]. In addition, since the DG method allows discontinuities of the physical unknowns within the interior of the problem domain, it is a natural approach to capture the jumps across the material interface in coupled problems. Above all, DG methods are also characterized by their flexibility in terms of mesh design while keeping their high order accuracy [28] and their high scalability in parallel simulations [50, 8] while optimal convergence rates are still achieved. However, if not correctly formulated, discontinuous methods can exhibit instabilities, and the numerical results fail to approximate the exact solution. It is, therefore, important to develop a DG FEM which leads to reliable results for a wide

Title Suppressed Due to Excessive Length

3

variety of problems. By using an adequate inter element flux definition combined to stabilization techniques, the shortcomings of non-stabilized DG methods can be overcome [49, 52, 54]. Since the seminal work of Reed et al. [53], DG methods have been developed to solve hyperbolic, parabolic, and elliptic problems [13]. Most of DG methods for elliptic and parabolic problems rely on the Interior Penalty (IP) method. The main principle of IP, as introduced in [15, 69], is to constrain weakly the compatibility through the use of compatibility and/or stabilization terms, instead of building it into the finite element, which enables the use of discontinuous polynomial spaces of high degree. The interest in the symmetric interior penalty (SIPG) methods, in which the compatibility terms symmetrizes the formulation and which will be considered in this work, has been renewed by Wheeler [69] due to demands for optimality in the convergence rates with the mesh size hs , i.e. the rates of the convergence is k in the H1 -norm and k + 1 in the L2 -norm, where k is the polynomial approximation degree. However there exist different possible choices of traces and numerical fluxes as discussed by Arnold et al. [5], who have provided an analysis of a large class of discontinuous methods for second order elliptic problems with different numerical fluxes, and demonstrated that correctly formulated IP, NIPG (Non-Symmetric Interior Penalty), LDG (Local discontinuous Galerkin), and other DG methods are consistent and stables methods. In particular Arnold et al. [5] have proposed a framework for dealing with linear elliptic problems by means of DG methods and demonstrated that DG methods which are completely consistent and stable achieve optimal error estimates, and that the inconsistent DG methods like the pure penalty methods can still achieve optimal error estimates provided they are super-penalized. Besides, Georgoulis [19] has derived anisotropic hp-error bounds for linear second order elliptic diffusion convection reaction using Discontinuous Galerkin finite element methods (SIPG and NIPG), on shape-regular and anisotropic elements, and for isotropic and anisotropic polynomial degrees for the element bases. He has also observed optimal order of convergence in the L2 -norm for the SIPG formulation when a uniform mesh size refinement for different values of k is employed. Moreover, he has shown that the solution of the adjoint problem suffers from sub-optimal rates of convergence when a NIPG formulation is used. Yadav et al. [71] have extended the DG methods from a linear self-adjoint elliptic problem to a second order nonlinear elliptic problem. The non-linear system resulting from DG methods is then analyzed based on a fixed point argument. They have also shown that the error estimate in the L2 -norm for piece-wise polynomials of degree k ≥ 1 is k + 1. They have also provided numerical results to illustrate the theoretical results. Gudi et al. [23] have proposed an analysis for the most popular DG schemes such as SIPG, NIPG, and LDG methods for one dimension linear and non-linear elliptic problems, and the error estimate has been studied for each of these methods by reformulating the problems in a fixed point form. In addition, according to Gudi et al. [22], optimal errors in the H1 -norm and in L2 -norm are proved for SIPG for polynomial degrees larger or equal to 2, and a loss in the optimality in the L2 -norm is observed for NIPG and LDG. In that work a deterioration in the order of convergence in the mesh size hs is noted when linear polynomials are used. DG FEM have been widely developed to solve mechanical problems such as non-linear solid mechanics [49, 62, 61], (nearly) incompressible elasticity [25, 39], strain-gradient elasticity [16], strain gradient plasticity [43] strain gradient damage models [67], plate equations [24, 16, 6, 66], shell equations [48], as a non-exhaustive

4

Lina Homsi, Ludovic Noels

list. Truster et al. [63, 64] have derived a Variational Multiscale Discontinuous Galerkin (VMDG) method to account for geometric and material non-linearities in which computable expressions emerge during the course of the derivation for the stability tensor and numerical flux weighting tensors. Recently, DG has been used to solve coupled problems. For instance Wheeler and Sun [60] have proposed a primal DG method with interior penalty (IP) terms to solve coupled reactive transport in porous media. Liu et al. [40] have developed an incomplete IP method for convection-dominated thermo-poro-mechanics. Furthermore, Zheng et al. [73] have proposed a DG method to solve thermo-elastic coupling due to temperature and pressure dependent thermal contact resistance. In that work the DG method is used to simulate the temperature jump, and the mechanical sub-problem is solved by the DG finite element method with a penalty function. In [27], a coupled nonlinear Electro-Thermal DG method, has been derived by the authors in terms of energetically conjugated fields gradients and fluxes. This conjugated pair has been obtained by a particular choice of the test functions (δfT = δ( T1 ), δfV = δ( VT )) and of the trial functions (fT = T1 , fV = −V T ), where T is the temperature and V is the electrical potential [42, 72, 38], which has allowed developing a stable nonlinear DG formulation with optimal convergence rate. The main aim of this work is to derive a consistent and stable DG FEM for ETM coupling analyzes, which, to the authors knowledge, has not been introduced yet. To this end, the constitutive equations governing the ETM coupling are formulated in Section 2 as a function of the displacement u , the electric potential 1 u, −V V and the temperature T , under the form f (u T , T ). In Section 3, the DG method is first formulated in finite deformations and finite fields variations, resulting into a set of non-linear coupled equations, and then implemented within a three-dimensional finite element code. In particular, the parallel feature of the finite element method relies on the ghost elements method developed for DG formulation [8, 70]. Afterwards, the uniqueness and optimal numerical properties are derived for Electro-Thermo-Elasticity stated in a small deformation setting in Section 4. In particular, the convergence rates of the error in both the energy and L2 -norms are shown to be optimal with respect to the mesh size hs in terms of the polynomial degree approximation k (respectively in order k and k + 1). This section concludes with a numerical test supporting the developed theory, including the scalability property of the implementation. Finally the methodology is applied to study the Electro-Thermo-Mechanical behavior of SMPC in Section 5, in which a simple transversely isotropic hyperelastic formulation is used to model carbon fibers and an elasto-viscoplastic large deformation formulation is considered to model the SMP. In particular, the shape memory effect of SMPC unit cells electrically activated is studied in the large deformation regime.

2 Strong formulation of the Electro-Thermo-Mechanical problem

In this section an overview of the basic equations that govern Electro-ThermoMechanical coupled phenomena is presented. The body in its reference configuration Ω0 ∈ Rd , where d is the spatial dimension, whose Dirichlet boundary ∂D Ω0 and Neumann boundary ∂N Ω0 are the outer boundaries ∂Ω0 of the domain, is

Title Suppressed Due to Excessive Length

5 𝜕N Ω0h

Ω0

𝜕N Ωe0

𝜕D Ω

𝜕D Ω0 𝒙 = 𝝋(𝑿)

𝑿

𝒏 𝑬𝑍

𝑵

𝜕D Ωe0

Ω

𝜕I Ω0

𝑿−

𝜕D Ω0h

𝜕N Ω

𝜕N Ω0

Ωe+ 0

𝑵− 𝒙

𝑵+ Ωe− 0

s

𝑿+ 𝜕Ωe+ 0

𝜕Ωe− 0

𝑬𝑌 𝜕I Ω0h

𝑬𝑋

(a) Reference and current configurations

(b) Body discretization

Fig. 1 Definition of (a) the reference and current configurations of the body, (b) the body discretization.

subjected to a deformation mapping x = ϕ (X) defining the current configuration Ω ∈ Rd , see Fig. 1(a). The deformation gradient tensor is defined as F =

∂x F), = x ⊗ ∇0 with J = det(F ∂X

(1)

∂ is the gradient with respect to the reference its Jacobian, and where ∇0 = ∂X configuration. The displacement field reads u = x − X. Finally, the material properties may in general depend on the position X.

2.1 Balance equations The first balance equation is the equation of motion which is the balance of linear momentum in the absence of (inertial and external) body forces with respect to the reference frame ∇0 · P T = 0 ∀X ∈ Ω0 , (2) where P is the first Piola-Kirchhoff tensor, which is expressed in terms of the Cauchy stress as P = σ · F −T J . (3) The second balance equation is the electrical contribution which is the conservation of the electric current density flow. Defining the flow of electric current density with respect to the reference surface Je = je · F −T J ,

(4)

where je is the flow of electric current density with respect to the current surface, the second balance equation reads ∇0 · Je = 0 ∀X ∈ Ω0 .

(5)

The third balance equation is the conservation of the energy flux. The energy flux with respect to the reference surface reads Jy = Q + V Je = jy · F −T J = (q + V je ) · F −T J ,

(6)

6

Lina Homsi, Ludovic Noels

where Q and q are respectively the heat flux per unit surface in the reference and current configurations, V is the electrical potential, and where jy is the energy flux with respect to the current surface. Then the conservation of energy flux is stated as ∇0 · Jy = −ρ0 y˙ + F¯ ∀X ∈ Ω0 , (7) where y is the internal energy per unit mass and F¯ represents all the body energy sources per unit reference volume.

2.2 Constitutive models The first Piola-Kirchhoff stress tensor is evaluated under the generic form F , F˙ , V, T ; ξ (τ < t)) , P = P(F

(8)

where T is the temperature, and ξ (τ < t) is the set of internal variables evaluated at time τ lower than the current time t. In the applications, we will consider thermo-mechanical transverse isotropic and SMP thermo-mechanical constitutive behaviors, see A. In this work, we consider the thermo-electric coupling written in the current configuration under the matrix form      l αll −∇V je , (9) = j= −∇T jy V l + αTll k + αV l + α2 Tll where k denotes the symmetric tensor of thermal conductivity coefficients, l denotes the symmetric tensor of electric conductivity coefficients, and α is the Seebeck coefficient. The coefficients can be temperature and electric potential depen∂ dent. Moreover ∇ = ∂x is the gradient with respect to the current configuration. Using Eqs. (4) and (6), this last set of equations is rewritten      L L αL −∇0 V Je J= , (10) = L K + αV L + α2 TL L −∇0 T Jy V L + αTL where F , T, V ) = F −1 · k (T, V ) · F −T J and K (F F , T, V ) = F −1 · l (T, V ) · F −T J , L (F

(11)

are the coefficient tensors expressed in the reference configuration. However, in the relation (9), the vectors J and (∇0 V T ∇0 T T )T are not energetically conjugated, see the discussions in [42, 72, 38, 27]. Therefore we define the fields M = (fV fT )T with fV = − VT and fT = T1 , and the gradients of the fields vector in the reference frame are defined by ∇0M , a 2d × 1 vector in terms of    1   ∇0 fV ∇0 V − T I TV2 I ∇0 M = = . (12) ∇0 fT ∇0 T 0 − T12 I As a result the set of Eqs. (10) is rewritten as F , fV , fT ) ∇0M , J = Z 0 (F

(13)

Title Suppressed Due to Excessive Length

7

with the coefficient matrix expressed in the reference configuration  F , fV , fT ) = Z 0 (F

F , fT ) L 2 (F F , fT , fT ) L 1 (F F , fT , fT ) J y 1 (F F , fV , fT ) L 2 (F

 ,

(14)

where F , fT ) = L 1 (F F , fV , fT ) = J y 1 (F

1 f 1 F , fT , fT ) = − V L ; L 2 (F L + α 2 L ; and 2 fT fT fT 2 fV fV 1 2 1 K − 2α L + α L + 2 3 3 3L. fT fT fT fT

(15)

From Eq. (13), it can be seen that the symmetric coefficients matrix Z 0 is positive −1 1 definite since L 1 and J y 1 − L T 2 · L 1 · L 2 = f 2 K are positive definite for T > 0. T

2.3 Boundary conditions The body boundary ∂Ω0 , see Fig. 1(a), is divided into a Neumann part ∂N Ω0 on which the surface traction P ·N , the electric current density Je ·N , and the energy density Jy · N per unit reference surface, are respectively constrained to T¯ , J¯e , and J¯y , and into a Dirichlet part ∂N Ω0 , on which the displacement field u, the ¯ f¯V , and f¯T . In these fV -field, and the fT -field are respectively constrained to u, relations, N is the outward unit normal in the reference configuration. Note that to simplify the notations of the equations we have assumed that the Neumann and Dirichlet parts coincide for the three fields, but in all generalities the developed methodology remains applicable if they are different.

2.4 Strong form summary The set of governing equations (2, 5, and 7) is thus stated as finding u, M ∈  2 d + H (Ω0 ) × H2 (Ω0 ) × H2 (Ω0 ) such that ( ∇0 · P T = 0 J) = I i ∇T 0 (J

∀X ∈ Ω0 , ∀X ∈ Ω0 ,

(16)

where ∇0 is also used to represent a vector operator in the reference configuration, T I i = 0 − ρ0 y˙ + F¯ represents the contributions of the internal energy rate and +

the body energy sources, and where the notation H2 (Ω0 ) is used to consider only positive values in the Hilbert space H2 . This set of governing equations is completed by the Neumann boundary conditions  ¯ ∀X ∈ ∂N Ω0 ,  P · N = T !T (17) N 0  J = J¯ ∀X ∈ ∂N Ω0 , J N = 0 N

8

Lina Homsi, Ludovic Noels

where J¯ = (J¯e J¯y )T and J N is the flux per unit reference surface, and by the Dirichlet boundary conditions (

∀X ∈ ∂D Ω0 , ∀X ∈ ∂D Ω0 ,

¯ u=u ¯ M =M

(18)

¯ = (f¯V f¯T )T . where M Finally, the constitutive Eqs. (13, 8) read ( F , F˙ , M ; ξ (τ < t)) , P = P(F F , M ) ∇0M . J = Z 0 (F

(19)

3 Weak Discontinuous Galerkin formulation In this section we first derive the weak DG formulation of the problem stated under the strong form (16-19). Then we introduce a finite-element discretization of the test and trial functions, before briefly summarizing the resolution process.

3.1 Weak DG form derivation Let Ω0h be a shape regular family of discretization of Ω0 , such that Ω0h = ∪e Ω0e , see Fig. 1(b), with hs = maxΩ0e ∈Ω0h diam(Ω0e ) for Ω0e ∈ Ω0h with ∂Ω0e = ∂N Ω0e ∪ ∂D Ω0e ∪ ∂I Ω0e , and where ∂I Ω0h = ∪e ∂I Ω0e \ ∂Ω0h , is the intersecting boundary of the finite elements. Finally (∂DI Ω0 )s is a face either on ∂I Ω0h or on ∂D Ω0h , with P s s (∂DI Ω0 ) = ∂I Ω0h ∪ ∂D Ω0h . Since, at the interface between two elements, Fig. 1(b), each interior edge − (∂I Ω0 )s is shared by two elements − and + , where (∂I Ω0 )s ⊂ ∂I Ω0e and (∂I Ω0 )s ⊂   + ∂I Ω0e , we can define two useful operators, the jump operator J·K = ·+ − ·− that computes the discontinuity between the elements and the average operator h·i = − 1 + which is the mean between two element values. Those two operators 2 · +· can be extended on the Dirichlet boundary ∂D Ω0h as h·i = ·, J·K = (−·). The discontinuous Galerkin method results from the integration by parts on the elements of the governing equations multiplied by discontinuous test functions. Let us multiply the governing equations (16) by discontinuous test functions δu M, and integrate on Ω0h , yielding and δM 0=

XZ

h id ∀δu ∈ Πe H1 (Ω0e ) ;

(20)

Ω0e

e

XZ e

P (F F , M ) · ∇0 ) · δudΩ0 (P

Ω0e

MdΩ0 = IT i δM

XZ e



T F , M , ∇0M ) δM MdΩ0 ∇T 0 J (F

Ω0e

M ∈ Πe H1 (Ω0e ) × Πe H1 (Ω0e ) . ∀δM

(21)

Title Suppressed Due to Excessive Length

9

By performing an integration by parts on each element, and by using the divergence theorem, these equations become Z

δu · T¯ dS0 = ∂N Ω0h Z

Z F , M ) : ∇0 δudΩ0 + P (F Ω0h

∂I Ω0h

F , M )K · N − dS0 Jδu · P (F

Z

h id F , M) · N dS0 ∀δu ∈ H1 (Ω0e ) ;(22) δu · P (F ∂D Ω0h Z Z Z T MdΩ0 + MTJ¯dS0 = M)T J (F F , M , ∇0M )dΩ0 + − I i δM δM (∇0 δM Ω0h ∂N Ω0h Ω0h Z Z r z T T MN F , M , ∇0M ) dS0 − MN F , M , ∇0M )dS0 δM J (F δM J (F −

∂I Ω0h

∂D Ω0h

M ∈ Πe H1 (Ω0e ) × Πe H1 (Ω0e ) .(23) ∀δM  where the vector MN =

N- 0 0 N-

 M has been introduced for simplicity, where −

N − is defined as the reference outward unit normal of the minus element Ω0e , whereas N + is the reference outward unit normal of its neighboring element, with N − , and where we have used the Neumann boundary conditions (17). N + = −N However, when considering the jump operators, consistency is satisfied as long as the jump on the test function is accounted for, while the stress terms can be substituted by a numerical flux such as the average value. The set of Eqs. (22-23) thus becomes Z Z ¯ F , M ) : ∇0 δudΩ0 + δu · T dS0 = P (F ∂N Ω0h Ω0h Z id h P (F F , M )i · N − dS0 ∀δu ∈ Πe H1 (Ω0e ) ;(24) JδuK · hP ∂I Ω0h ∪∂D Ω0h Z Z Z M)T J (F F , M , ∇0M )dΩ0 + MTJ¯dS0 = MdΩ0 + (∇0 δM δM IT − i δM Ω0h ∂N Ω0h Ω0h Z r z T MN J(F F , M , ∇0M )i dS0 ∀δM M ∈ Πe H1 (Ω0e ) × Πe H1 (Ω0e ) ,(25) δM hJ ∂I Ω0h ∪∂D Ω0h

where we have considered the definitions of the operators on the Dirichlet boundary. In order to define the compatibility and stability terms, we define, on the one P F = I , M = M 0 ) and, on the other hand, hand, the four-order tensor H = ∂P F (F ∂F P F the second order tensor α th such that α th : H = − ∂P ∂T (F = I , M = M 0 ), with + M0 ∈ Πe H1 (Ω0e ) × Πe H1 (Ω0e ) the initial values of M. As a result both H and α th : H are constant during the simulation as it has been shown that this leads to accurate results in elasto-plasticity [50]. However, evolving terms could also be considered as discussed in [63, 64]. Moreover, we define the d × d × 2 matrix P ∂P M) = ∂M Y (M M . Assuming the stress tensor depends only on fT and not on fV , we 1 M)δM M = α th : H 2 δfT , with Y 0 = Y (M M0 ). have Y (M fT

10

Lina Homsi, Ludovic Noels

The compatibility of the displacement fields on the element interfaces and on the Dirichlet boundary is thus added to the weak form (24-25), leading to

Z

Z H : ∇0 δu) · N dS0 − ¯ · (H δu · T¯ dS0 − u ∂N Ω0h ∂D Ω0h Z  ¯ )M ¯ − Y (M ¯ 0 )M ¯ 0 · N dS0 = δu · Y (M ∂D Ω0h Z Z F , M ) : ∇0 δudΩ0 + P (F F , M )i · N − dS0 + JδuK · hP P (F Ω0h ∂I Ω0h ∪∂D Ω0h Z H : ∇0 δui · N − dS0 + JuK · hH ∂I Ω0h ∪∂D Ω0h Z h id Y (M M)M M − Y (M M0 )M M0 i · N - dS0 ∀δu ∈ Πe H1 (Ω0e ) (26) JδuK · hY ; ∂D Ω0h Z Z Z  T T¯ ¯N ¯ )∇0 δM M M M dS0 = − IT δM dΩ + δM J dS − M Z 0 (F , M 0 0 i Ω0h ∂N Ω0h ∂D Ω0h Z Z r z T MN J(F F , M , ∇0M )i dS0 M)T J (F F , M , ∇0M )dΩ0 + δM hJ (∇0 δM Ω0h ∂I Ω0h ∪∂D Ω0h Z Z z z r r T

T ¯ )∇0 δM F, M M dS0 Z0 (F F, M )∇0 δM Mi dS0 + Z 0 (F hZ MN + MN ∂I Ω0h

∂D Ω0h

M ∈ Πe H1 (Ω0e ) × Πe H1 (Ω0e )(27) ∀δM .

Some remarks arise from this formulation

– In the compatibility terms, the operator H has been considered as constant but not the operator Z 0 . As a result, the optimality of the convergence rate will only be demonstrated when linearizing the mechanical equations, and when stating the equations in a small deformation setting. This is due to the complexity of the mechanical behavior (19a) which is history-dependent, preventing to write a non-linear equation under a form similar to the considered electro-thermal coupling (19b). P (F F , M )i· – The last term of Eq. (26) appears from the fact that the term in JδuK·hP N − dS0 acting on the Dirichlet boundary condition prevents a symmetric operator to be obtained and as such would prevent the optimal convergence rate in the L2 -norm to be achieved, as it will be shown in Section 4.

What remains to be done is to stabilize the method through quadratic terms weighted by a stabilization parameter B and the mesh size hs , leading to state the

Title Suppressed Due to Excessive Length

11

 d + problem as finding u, M ∈ Πe H1 (Ω0e ) × Πe H1 (Ω0e ) × Πe H1 (Ω0e ) such that Z

δu · T¯ dS0 −

Z

H : ∇0 δu) · N dS0 + ¯ · (H u   Z HB ¯⊗N : : δu ⊗ N dS0 − u hs ∂D Ω0h Z Z  ¯ )M ¯ − Y (M ¯ 0 )M ¯ 0 · N dS0 = F , M ) : ∇0 δudΩ0 + δu · Y (M P (F ∂D Ω0h Ω0h Z P (F F , M )i · N − dS0 + JδuK · hP ∂I Ω0h ∪∂D Ω0h   Z HB JuK ⊗ N : : JδuK ⊗ N dS0 + hs ∂I Ω0h ∪∂D Ω0h Z H : ∇0 δui · N − dS0 + JuK · hH ∂I Ω0h ∪∂D Ω0h Z h id Y (M M)M M − Y (M M0 )M M0 i · N - dS0 ∀δu ∈ Πe H1 (Ω0e ) ; (28) JδuK · hY ∂D Ω0h Z Z Z  T T ¯N ¯ )∇0 δM MdΩ0 + MTJ¯dS0 − M dS0 + − I i δM δM M Z 0 (F , M Ω0h ∂N Ω0h ∂D Ω0h   Z Z B T ¯) M ¯ N dS0 = F ,M M)T J (F F , M , ∇0M )dΩ0 + MN Z 0 (F (∇0 δM δM hs Ω0h ∂D Ω0h  Z r z B T MN F, M ) JMN K dS0 + δM Z 0 (F hs ∂I Ω0h  Z r z B T ¯ ) JMN K dS0 + MN F, M δM Z 0 (F hs ∂ Ω Z D 0h z r T J(F F , M , ∇0M )i dS0 + MN hJ δM ∂I Ω0h ∪∂D Ω0h Z Z z z r r T

T ¯ )∇0 δM F, M M dS0 Z0 (F F, M )∇0 δM Mi dS0 + Z 0 (F hZ MN MN ∂N Ω0h

∂D Ω0h

∂D Ω0h

∂I Ω0h

M ∈ Πe H1 (Ω0e ) × Πe H1 (Ω0e ) ,(29) ∀δM  N 0 ¯ =N ¯ MM ¯ has been introduced. M 0 N The manifold of the trial functions G T = (uT M T ) is defined as

¯N = where the notation M

( X(+) s

=



)  d (+) G ∈ L2 (Ω0h ) × L2 (Ω0h ) × L2 (Ω0h ) such that . (+) G |Ω0e ∈ [Hs (Ω0e )]d × Hs (Ω0e ) × Hs (Ω0e ) ∀Ω0e ∈ Ω0h (+)

(30)

For the future use, we define X(+) as X2 and X+ the manifold such that fT > 0, while X is the manifold for which fT Q 0, with X+ ⊂ X. It should be noted that the trial functions in the previous equations of the weak formulation belong to  1 e d + H (Ω0 ) × H1 (Ω0e ) × H1 (Ω0e ); however for the numerical analysis, we will need  2 e d + to be in H (Ω0 ) × H2 (Ω0e ) × H2 (Ω0e ), as it will be discussed in Section 4.

12

Lina Homsi, Ludovic Noels

Therefore the weak form (28-29) is reformulated as finding G ∈ X+ such that ( F , M ; δu) = B(δu) A(F R G ∈ X, ∀δG (31) F , M ; δM M) = D(F F ; δM M) − Ω δM MTI i dΩ0 C(F 0h F , M ; δu) (B(δu)) is directly deduced from the right-hand (left-hand) where A(F F , M ; δM M) (D(F F ; δM M)) is directly deduced from the side of Eq. (28), and where C(F right-hand (left-hand) side of Eq. (29).

3.2 The finite element discretization of the coupled problem  T In the computational model we consider the approximation G T h = uh fVh fTh of the trial function defined in a finite dimensional space of real valued piece-wise polynomial functions. The following manifold is thus introduced ( )  d (+) G h ∈ L2 (Ω0h ) × L2 (Ω0h ) × L2 (Ω0h ) such that k(+) X = , (32) (+) G h |Ω0e ∈ [Pk (Ω0e )]d × Pk (Ω0e ) × Pk (Ω0e ) ∀Ω0e ∈ Ω0h +

where Pk (Ω0e ) is the space of polynomial functions of order up to k and Pk means that the polynomial approximation remains positive. The discretization of the system is carried out using the discontinuous Galerkin Finite element method by introducing the same shape functions for the trial funcM, which are thus interpolated tions u and M , and for the test functions δu and δM as a uh = Nua u a , M h = NM Ma , (33) ua , δuh = Nua δu

a M h = NM Ma , δM δM

(34)

ua and (δ)M Ma denote Mh where (δ)u nodal  the  values of respectively (δ)uh and (δ)M NfaV 0 a at node a, and where NM = is a matrix of the shape functions. 0 NfaT The finite element approximation of the weak form (31) is thus stated as finding + G h ∈ Xk such that ( F h , M h ; δuh ) = B(δuh ) A(F R Gh ∈ Xk . ∀δG (35) F h , M h ; δM Mh ) = D(F F h ; δM Mh ) − Ω δM MT C(F I ih dΩ0 h 0h 3.3 The system resolution in parallel The set of Eqs. (35) can be rewritten under the form       F aext G b = F aint G b + F aI G b ,

(36)

where G b is a (d + 2) × 1 vector gathering the unknowns at node b with G b = T T F h ; δM Mh ) at ub M b )T , F aext corresponds to the contributions of B(δuh ) and D(F (u F h , M h ; δuh ), C(F F h, node a, F aint corresponds to the volume contributions of A(F Mh ), and I ih , and where F aI corresponds to the interface contributions of M h ; δM

Title Suppressed Due to Excessive Length

13

F h , M h ; δuh ) and C(F F h , M h ; δM Mh ). The expressions of these forces are detailed A(F in B.1. The non-linear Eqs. (36) are linearized by means of an implicit formulation and solved using the Newton-Raphson scheme using as initial guess the last solution. To this end, the forces are written in a residual form. The predictor at iteration 0, reads G c = G c0 , and the residual at iteration i reads  i  i  i  i F aext G c − F aint G c − F aI G c = R a G c .

(37)

At iteration i, the first-order Taylor development yields the system to be solved, i.e.  a   i Fext Faint FaI ∂F ∂F ∂F Gb = −R Ra G c . (38) − − |G =G ci ∆G G b b b G G G ∂G ∂G ∂G The explicit expression of the tangent matrix of the coupled Electro-ThermoFa Fa Fa ∂F ∂F ∂F ab ext int I Mechanical system KG = ∂G Finally, the resolution Gb − ∂G Gb − ∂G Gb is given in B.2.  Gb = G b − G bi . of the system (38) yields the correction ∆G

Ω𝐼0h

Ω𝐼0h

Ωe0

Ωe0

g𝐼𝐼𝐼

g𝐼𝐼𝐼

g𝐼𝐼

Ω0

Ω0

Ω0

g𝐼

Ω0

g𝐼

Ω0

g𝐼𝐼𝐼 Ω0

𝐆𝒂

g𝐼𝐼 Ω0

Ω𝐼𝐼𝐼 0h

Ω𝐼𝐼 0h

(a)

Ω𝐼𝐼𝐼 0h

Ω𝐼𝐼 0h

(b)

i with Fig. 2 Parallel implementation using ghost elements with (a) the different partitions Ω0h j their ghost elements Ω0gj corresponding to the elements of the partition Ω0h having a common j interface, and (b) the exchange of the nodal field G b from the element in partition Ω0h to the gj i ghost element Ω0 in partition Ω0h .

The DG method has been implemented in Gmsh [20] in parallel using the ghost elements method suggested in [8, 70]. In this approach, the discretization Ω0h is i divided in partitions Ω0h using the METIS library [31], see Fig. 2(a). On top of its e i own bulk elements Ω0 , the partition Ω0h also owns the ghost elements Ω0gj which j correspond to the bulk elements of the partition Ω0h sharing a common partitions interface, see Fig. 2(a). The resolution steps are thus

14

Lina Homsi, Ludovic Noels

i – Step #1: Each partition Ω0h evaluates the bulk nodal forces F aint (36) and stiffFa ∂F e int ness matrix ∂G Gb (38) by performing a quadrature rule on its own elements Ω0 . There is no need to evaluate the bulk force and stiffness matrix contributions of the ghost elements Ω0gj . i – Step #2: Each partition Ω0h evaluates the interface nodal forces F aI (36) and Fa ∂F I stiffness matrix ∂G Gb (38) by performing a quadrature rule on an interface i shared either by two elements Ω0e belonging to this partition Ω0h or by an gj e element Ω0 belonging to this partition and a ghost element Ω0 . – Step #3: The system (38) is then solved using the MUMPS [3] library. – Step #4: To evaluate correctly the interface forces at partitions boundaries at Step #2 of the next iteration, the ghost elements Ω0gj need to be in the updated deformation state, which requires their nodal values G a to be communicated j from the original element Ω0e lying in partition Ω0h , see Fig. 2(b). This communication is achieved before each new iteration through the network via Message Passing Interface (MPI), whichs is the only communications required by the method.

4 Numerical Properties The demonstration of the numerical properties for Electro-Thermo-Mechanics coupling is derived following closely the approach developed by Gudi et al. [23] for nonlinear problems under the assumption d = 2, and under the assumptions of temperature independent thermo-mechanical material properties, (however J y1 , L 1 , L 2 remain temperature and electric potential dependent but C (the matrix form using Voigt notations of the material constant tensor H ), and α th are temperature and electric potential in-dependent), and in the absence of the heat source, such that the term F¯ in Eq. (16) is equal to zero. Finally, we also require a framework in small deformation and linear thermo-elasticity in order to demonstrate the stability and convergence rates. Indeed, as explained before the thermo-mechancal response in the non-linear range is history-dependent, preventing to write a nonlinear equation under a form similar to the onr of the electro-thermal coupling (19b).

4.1 Strong form in the small deformation setting Under the discussed assumptions, we can state the weak form (31) under a matrix form with the vector of the unknown fields G . In addition, we can  introduce  the C 0 0 coefficient matrix v of size (5d − 3) × (5d − 3) such that v =  0 l 1 l 2 . In 0 l 2 j y1 a small deformation setting, we have l 1 ' L 1 , l 2 ' L 2 , and j y1 ' J y1 , and thus Z ' Z 0 , the coefficient   matrix in the current configuration. Therefore v can also be C 0 written as v = . Moreover, assuming that the stress tensor does not depend 0Z

Title Suppressed Due to Excessive Length

15

 C α thc −C 0 0 2   fT  and on the electrical potential, we define the matrices o (fT ) =   0 0 0 00 0 o 0 = o (fT0 ) of size (5d − 3) × (d + 2). In these matrices, α thc is a vector of size T (3d − 3) × 1 with, for d = 3, αthc = αth αth αth 0 0 0 , and C αthc a vector of T C α thc ) = 3Kαth 3Kαth 3Kαth 0 0 0 size (3d − 3) × 1 and given for d = 3 by (C for isotropic materials, with K the bulk modulus and αth the thermal dilatation.  00 0 Finally we define h a matrix of size (d+2)×(d+2) with, for d = 3, h =  0 0 0 . 0 0 ρcv G, ∇G G) = v (G G)∇G G. We can now introduce the matrix w of size (5d−3)×1 as w (G G is the (5d − 3) × 1 vector of the gradient of the unknown fields, In this relation ∇G G = (∇)G G and is written for d = 3, using Voigt rules for the which is defined as ∇G mechanical contribution, as 

  ∂ εxx ∂x  εyy   0     εzz   0    ∂  2εxy    ∂y    2εxz   ∂   ∂z   2εyz     0  G) =  ∂fV  =  (∇G  ∂x   0  ∂fV    ∂y   0  ∂fV    ∂z   0  ∂fT    ∂x     ∂fT   0  ∂y   0 ∂fT 0 ∂z 

0 ∂ ∂y

0 0

0

∂ ∂z

∂ ∂x

0

∂ ∂x ∂ ∂ ∂z ∂y

0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0 0 0 0 0



        ux    uy      uz  . ∂     fV  ∂x  ∂  fT ∂y  ∂  ∂z ∂  0 ∂x  ∂   0 ∂y ∂ 0 ∂z

(39)

 Similarly, Σ T = σxx σyy σzz τxy τxz τyz jex jey jez jyx jyy jyz represents the energy conjugated stress. Using these definitions, the constitutive relations (19) are rewritten as G)∇G G + o (G G)G G − o 0G 0 , Σ = v (G

(40)

where G 0 is the vector of the initial values G T 0 = (ux0 uy0 uz0 fV0 fT0 ). The set of  d + governing Eqs. (16) is now stated as finding G ∈ H2 (Ω0 ) × H2 (Ω0 ) × H2 (Ω0 ) such that ˙ w(G G, ∇G G) + o (G G)G G − o 0G 0 ] = h G −∇T [w

∀x ∈ Ω .

(41)

This set of equations is completed with the Dirichlet BC (18) rewritten as ¯ G =G

∀ x ∈ ∂D Ω ,

(42)

and with the Neumann BC (17) rewritten as w + oG − o 0G 0 ) = w ¯ n¯T (w

∀ x ∈ ∂N Ω,

(43)

16

Lina Homsi, Ludovic Noels

where, for d = 3, 

nx  0   0   ny   nz   0 n¯ =   0   0   0   0   0 0

0 ny 0 nx 0 nz 0 0 0 0 0 0

0 0 nz 0 nx ny 0 0 0 0 0 0

0 0 0 0 0 0 nx ny nz 0 0 0

 0 0   0   0   0   0  , 0   0   0   nx   ny  nz

(44)

 ¯T = u ¯ T f¯V f¯T represents the unit outward normal in the current configuration, G  ¯ T = t¯T ¯jy ¯je gathers the constrained gathers the constrained fields, and where w fluxes. G)G G) consists of zero components and It can be noticed that the gradient of (oo(G Cα thc ) C C (C of the gradient of (− f 2 fT ). As the variation δ(− (C αf 2thc ) fT ) = (C αf 2thc ) δfT , the T T T G) can be rearranged in a new form o˜(G G) of size (d + 2) × (5d − 3), such matrix o (G that G)G G) = −˜ G)∇G G, ∇T (oo(G o (G

(45)

with, for d = 3, 

th 0 0 0 0 0 0 0 0 0 − 3Kα f2 T

0   G) =  0 o˜(G  0 0

00000000

0

0



th − 3Kα f2

0

0 0 0

th − 3Kα 2 fT 0 0

   .  

0 T

00000000 00000000 00000000

0 0 0

(46)

The operator ˜· can be seen as the transpose operator that accounts for the definition of the ∇ operator in the matrix form. Therefore Eq. (41) becomes ˙ w(G G, ∇G G)) + o˜(G G)∇G G = hG −∇T (w

∀x ∈ Ω .

(47)

4.2 Weak DG form in the small deformation setting Assuming a quasi-static process from now on, the associated DG form for the Electro-Thermo-Elasticity problem (47) is defined as finding G ∈ X+ such that G, δG G) = b(δG G), a(G

G ∈ X, ∀δG

(48)

Title Suppressed Due to Excessive Length

17

with Z

G)Tw(G G, ∇G G)dΩ + (∇δG

G, δG G) = a(G

Z

Ωh

G)∇G GdΩ + GTo˜(G δG

Ωh

Z

r

GnT δG

z

w(G G, ∇G G)i dS + hw Z r z r z

¯ )∇δG v(G G)∇δG Gi dS + G dS + GnT hv GnT v (G

∂I Ωh ∪∂D Ωh

Z ∂I Ωh

∂D Ω h

Z

r

GnT

 z  v (G G)B

∂I Ωh

Z ∂I Ωh ∪∂D Ωh

Z

Z

r

Gn K dS + JδG

hs ∂D Ω h E D T G)G G − o0G0 K dS + Gn Jo(G δG

GnT

z  v (G ¯ )B 

hs

Gn K dS − JδG

G)G G − o 0G 0 )dS , GTn¯T (oo(G δG

(49)

∂N Ω h

and Z G) = b(δG

¯ dS − G Tw δG

¯ nTv (G ¯ )∇δG G dS + G

∂D Ω h

∂N Ω h

Z

Z

¯

v (G )B ¯ GnT δG Gn hs ∂D Ωh

Z

¯ )G ¯ − o 0G 0 )dS , GnT (oo(G δG

dS +

(50)

∂D Ω h

where Gn = n¯−G is a (5d − 3) × 1 vector, which is defined by building n¯− using the unit outward normal of the “-” element, n− , following Eq. (44). Using the identity JabK = JaK hbi + hai JbK on ∂I Ωh , the relation (45), and performing an integration by parts lead to Z

Z T G o ˜ G G GT ∇T (oo(G G)G G − o0G0 )dΩ = δG (G )∇G dΩ = − δG Ωh Ωh XZ XZ G)T (oo(G G)G G − o 0G 0 )dΩ − GnT (oo(G G)G G − o 0G 0 )dS = (∇δG δG Ωe

e

Z

G)T (oo(G G)G G − o 0G 0 )dΩ − (∇δG Ωh

GnT (oo(G G)G G − o 0G 0 )dS δG

(51)

∂N Ωh

Z − Z

∂Ω e

e

Z

GnT (oo(G G)G G − o 0G 0 )dS + δG

∂D Ω h

∂I Ωh

D

Z

r

z

GnT hoo(G G)G G − o 0G 0 i dS+ δG

∂I Ωh

E

GnT Jo (G G)G G − o 0G 0 K dS. δG

Therefore, Eq. (48) can be rewritten as G, δG G) = b0 (δG G), a0 (G

G ∈ X, ∀δG

(52)

18

Lina Homsi, Ludovic Noels

with Z

0

G, δG G) = a (G

Z

T

G)T (oo(G G)G G − o 0G 0 )dΩ + (∇δG Ω Ωh Z h Z r z r z GnT hw w(G G, ∇G G)i dS + v(G G)∇δG Gi dS + δG GnT hv ∂I Ωh ∪∂D Ωh ∂I Ωh  Z Z r z

r z  v (G G)B ¯ )∇δG G dS + Gn K dS + JδG GnT v (G GnT hs ∂D Ω h ∂I Ωh Z r z  v (G ¯ )B  Gn K dS + JδG GnT hs ∂D Ω h Z r z GnT hoo(G G)G G − o 0G 0 i dS δG ∂I Ωh ∪∂D Ωh Z E D G)G G − o 0G 0 K dS , GnT Jo (G (53) − δG G) w (G G, ∇G G)dΩ + (∇δG

∂D Ω h

and G) = b0 (δG

Z

¯ dS − G Tw δG

Z ∂D Ωh

¯ nTv (G ¯ )∇δG G dS G

∂D Ωh

∂N Ωh

+

Z

GnT δG

¯ )B v (G ¯ n dS + G hs

Z

(54) ¯ )G ¯ − o 0G 0 )dS . GnT (oo(G δG

∂D Ωh

Henceforth, using Eq. (40), it is shown that Eq. (52), which is derived from Eq. (48), corresponds to the weak form (31). Unlike the usual case in DG, where the interface term involves o in the average operator h i, Eq. (49) shows that o is rather involved in the jump operator J K. This comes from the integration by parts in Eq. (51), in which o is G dependent. However, this allows the volume and consistency terms in Eq. (52) to be directly expressed in terms of the stress Σ = w + (ooG − o 0G 0 ), which is convenient when dealing with a non-linear formulation as in Eqs. (31).

4.3 Consistency  d To prove the consistency of the method, the exact solution G e ∈ H2 (Ω) × +

H2 (Ω) × H2 (Ω) of the problem stated by Eq. (47) is considered. This implies ¯ = −G wi = w , Jo (G Ge )G Ge − o 0G 0 K = 0 on ∂I Ωh , and JG e K = −G Ge , JG e K = 0, hw e e e e e ¯ ¯ ¯ G )G G − o 0G 0 K = −oo(G )G + o 0G 0 , hw wi = v (G )∇G G = v (G G )∇G G , and v (G G) = Jo (G ¯ ) = v(G v(G Ge ) on ∂D Ωh . Therefore, Eq. (48) becomes:  R R R ¯ )G ¯ − o 0G 0 dS + ¯ nTv (G ¯ )∇δG GnT o (G G Tw ¯ dS − ∂ Ω G G dS + ∂ Ω δG δG ∂N Ωh D h D h R R ¯ )B ¯ G T Ge Ge GnT v (G δG hs Gn dS = Ωh (∇δG ) w (G , ∇G )dΩ + ∂D Ω h r z R R GTo˜(G Ge )∇G Ge dΩ + ∂ Ω δG GnT hw w(G Ge , ∇G Ge )i dS − δG Ωh I h R R R T ¯ )G ¯ )∇δG Gne dS + GnTw (G Ge , ∇G Ge )dS − ∂ Ω Gne v (G GdS + ∂ Ω δG GnT hBs v (G δG ∂D Ωh D h D h R R GnT (oo(G Ge )G Ge − o 0G 0 ) dS + ∂ Ω δG GTn¯T (oo(G Ge )G Ge − o 0G 0 ) dS ∀δG G ∈ X. δG ∂D Ωh N h (55)

Title Suppressed Due to Excessive Length

19

Integrating the first term of the right hand side by parts leads to XZ XZ G)Tw (G Ge , ∇G Ge )dΩ = − GT ∇Tw (G Ge , ∇G Ge )dΩ+ (∇δG δG e

Ωe

e

XZ e

Ωe

GnTw (G Ge , ∇G Ge )dS, δG

∂Ω e

and Eq.(55) becomes  R R R ¯ nT v (G ¯ )∇δG ¯ )G ¯ − o 0G 0 )dS + G Tw ¯ dS − ∂ Ω G G dS + ∂ Ω δG GnT (oo(G δG ∂N Ω h D h D h   R R ¯) G ¯ n dS = − GT ∇Tw (G Ge , ∇G Ge )dΩ + GnT hBs v (G δG δG Ωh ∂D Ωh R R R T ¯ )∇δG GnTw (G Ge , ∇G Ge )dS + Ω δG GTo˜(G Ge )∇G Ge dΩ − ∂ Ω Gne v (G GdS + δG ∂N Ωh h D h R R ¯ )G Gne dS + ∂ Ω δG GnT (oo(G Ge )G Ge − o 0G 0 )dS + GnT hBs v (G δG ∂D Ωh D h R GTn¯T (oo(G Ge )G Ge − o 0G 0 )dS ∀δG G ∈ X. δG (56) ∂N Ω h The arbitrary nature of the test functions and the use of Eq. (45) lead to recover the set of conservation laws, Eqs. (41), and the boundary conditions, Eqs. (42-43).

4.4 Second order non-self-adjoint elliptic problem In this part, we will assume that ∂D Ωh = ∂Ωh . This assumption is not restrictive but simplifies the demonstrations. G), which is a symmetric and Starting from the definition of the matrix v (G positive definite matrix since its components C and Z are positive definite matrices, G) as λ(G G) let us define the minimum and maximum eigenvalues of the matrix v (G G); then for all ξ ∈ R(5d−3) and Λ(G 0 G)|ξ|2 ≤ ξiv ij (G G)ξj ≤ Λ(G G)|ξ|2 . 0 < λ(G

(57)

Also by assuming that k G kW1∞ ≤ α, then there is a positive constant Cα such that G) . 0 < Cα < λ(G Let us define    G ∈ (L2 (Ω0h ))(5d−3) |∇G Y = ∇G G

(58)

 |Ω e ∈

(H1 (Ω0e ))

(5d−3)

∀Ω0e ∈Ω0h

.

(59)

In the following analysis, we use the integral form of the Taylor’s expansions of G, ∇G G) = v (G G)∇G G, for (V V, ∇Q Q) ∈ X+ × Y in terms of (G G, ∇G G) ∈ X+ × Y : w (G V, ∇Q Q) − w (G G, ∇G G) = − wG (G G, ∇G G)(G G − V ) − w ∇G G)(∇G G − ∇Q Q)+ w (V G (G ¯ G − V , ∇G G − ∇Q Q) Rw (G

(60)

¯ G (G G, ∇G G)(G G − V) − w ¯ ∇G G)(∇G G − ∇Q Q) , =−w G (G G, ∇G G) = vG (G G)∇G G is the partial derivative of w with respect to G where wG (G G) the partial derivative of v with respect to G ), w ∇G G) = v (G G) is the (with vG (G G (G

20

Lina Homsi, Ludovic Noels

G expressed in the matrix form (and using partial derivative of w with respect to ∇G ¯ ¯G, w ¯ ∇G Voigt notations), and where w , G and Rw are the remainder terms, with in particular ¯ w (G G − V , ∇G G − ∇Q Q) =(G G − V )Tw ¯ GG (V V, ∇Q Q)(G G − V )+ R G − V )Tw ¯ G ∇G V)(∇G G − ∇Q Q) . 2(G G (V

(61)

¯ G, w ¯ ∇G ¯ GG , and w ¯ G∇G The other remainder terms w G, w G are given in C.1. Using the (5d−3) ¯G, w ¯ GG ∈ L ∞ (Ω × R(d+1) × R+ definition of w , if fT ≥ c2 > 0, then w ) 0 ×R + ∞ (d+1) ¯ ∇G ¯ G∇G and w (Ω × R × R0 ). G, w G ∈L G, ∇G G) = o˜(G G)∇G G a (d + 2) × 1 vector, and For future use, let us introduce d (G G) = o˜(G G) a (d + 2) × (5d − 3) matrix, dG (G G, ∇G G) = its partial derivatives d ∇G G (G G)∇G G a (d + 2) × (d + 2) matrix, dGG (G G, ∇G G) = o˜GG (G G)∇G G a (d + 2) × (d + o˜G (G G) = o˜G (G G) a (d + 2) × (5d − 3) × (d + 2) matrix. 2) × (d + 2) matrix, and d ∇G GG (G Similarly to Eqs. (60-61), one has V, ∇Q Q) − d (G G, ∇G G) = − dG (G G, ∇G G)(G G − V ) − d ∇G G)(∇G G − ∇Q Q)+ d (V G (G ¯ G − V , ∇G G − ∇Q Q) Rd (G

(62)

G, ∇G G)(G G − V ) − d¯∇G G)(∇G G − ∇Q Q) , = − d¯G (G G (G where ¯ d (G V, ∇Q Q)(G G − V )+ G − V , ∇G G − ∇Q Q) =(G G − V )Td¯GG (V R T¯ V)(∇G G − ∇Q Q) . G − V ) dG ∇G 2(G G (V

(63)

¯ ¯ G are given in C.1. Using the The other remainder terms d¯G , d¯∇G G , dGG , and dG ∇G 2 (5d−3) ¯ ¯ definition of d , if fT ≥ c > 0, then dG , dGG ∈ L ∞ (Ω × R(d+1) × R+ ) 0 ×R + ∞ (d+2) ¯ ¯ and d ∇G (Ω × R × R0 ). G , dG ∇G G, ∈ L G) = o (G G)G G and its first and Finally, we also define the (5d − 3) × 1 vector p (G G) of size (5d − 3) × (d + 2) and pGG (G G) of size (5d − 3) × second derivatives pG (G (d + 2) × (d + 2) respectively. The Taylor’s expansion of this last vector reads ¯ p (G V) − p (G G) = −p pG (G G)(G G − V) + R G − V ) = −¯ G)(G G − V) , p (V pG (G

(64)

where ¯ p (G G − V ) = (G G − V )Tp¯GG (V V)(G G − V) , R

(65)

with the other remainder terms p¯G and p¯GG given in C.1. Using the definition of p , if fT ≥ c2 > 0, then p¯G , p¯GG , ∈ L ∞ (Ω × R(d+1) × R+ 0 ). However, one has GTo T (G G)) oT (G G) ∂(G T T ∂o T G), which once computed explicitly as to pG = = G + o (G G G ∂G ∂G G), and Eq. (65) becomes derive Eq. (45) leads to pG = −oo(G ¯ p (G G − V ) = −(G G − V )¯ V)(G G − V) . R oG (V

(66)

2

Since for fT ≥ c > 0, w , o , and d are twice continuously differentiable function with all the derivatives through the second order locally bounded in a ball around G ∈ [R]d × R × R+ 0 as it will be shown in Section 4.5, we denote by Cy n Cy =max k w , d kW2 (Ω×R(d+1) ×R+ ×R(5d−3) ) , ∞

0

¯G, w ¯ GG , d¯G , d¯GG kL∞ (Ω×R(d+1) ×R+ ×R(5d−3) ) k w

(67)

0

o

¯ GG k ∞ ¯ ∇G ¯ G ∇G ¯G , d¯∇G k w G, w G, o G , d ∇GG L (Ω×R(d+1) ×R+ ) . 0

Title Suppressed Due to Excessive Length

21

 d + Let us define the solution G e ∈ H2 (Ω) × H2 (Ω) × H2 (Ω) of the strong e e e ¯ n on Gne = −G form stated by Eqs. (41-43). Since JGn K = 0 on ∂I Ω and JGn K = −G e ∂D Ω , and since Eq. (48) satisfies the consistency condition, we have Z Z e T e e G G G w G G GTo˜(G Ge )∇G Ge dΩ+ a(G , δG ) = (∇δG ) (G , ∇G )dΩ + δG Ωh Ωh Z Z r z T e e T Gn hw w(G G , ∇G G )i dS − Gn w (G Ge , ∇G Ge )dS− δG δG ∂I Ω h ∂D Ω h Z Z (68) T T v (G Ge )B Gn dS+ Ge )∇δG GdS + δG Gne v (G Gne hs ∂ Ω ∂D Ω h Z D h T e e Gn (oo(G G )G G − o 0G 0 )dS = b(δG G) ∀δG G ∈ X, δG ∂D Ωh

with Z

 ¯ nT v (G ¯ )∇δG G dS + G ∂D Ωh   Z B ¯ ¯ n dS. GnT v (G ) G δG hs ∂D Ω h

Z

¯ )G ¯ − o 0G 0 )dS+ GnT (oo(G δG

G) = − b(δG

∂D Ω h

(69)

G = δG Gh ∈ Xk in Eq. (68), subtracting the DG discretizaTherefore, using δG tion (48) E r Tsuccessively r Tfrom Eq.z (68), then adding and subtracting z D the two terms R R e T e T e B w G G Ge ) JδG Ghn K dS G − G G − G hw (G )∇δG i dS, and G G (G h ∇G n n hn hn hs w ∇G ∂ Ω ∂ Ω I

h

I

h

Ge )G Ge − o 0G 0 K = 0, JGne K = 0 on ∂I Ωh , and to this last relation, and using Jo (G e e e e ¯ ¯ )G ¯ + o 0G 0 on ∂D Ωh , one gets G )G G − o 0G 0 K = −oo(G Gn = −Gn , Jo (G JGn K = −G Z Gh )T (w w(G Ge , ∇G Ge ) − w (G Gh , ∇G Gh )) dΩ 0= (∇δG Ωh Z Ge )∇G Ge − o˜(G Gh )∇G Gh ) dΩ GT o (G + δG h (˜ Ωh Z z r w(G Ge , ∇G Ge ) − w (G Gh , ∇G Gh )i dS GT + δG hn hw ∂ Ω ∪∂ Ω Z I h D hr z T w∇G Ge )∇δG Gh i dS + Gne − G T G (G hn hw ∂I Ωh ∪∂D Ωh Z r T z T Ge ) − G T Gh ) hδG Ghn i dS − G e o T (G h o (G ∂ Ω ∪∂ Ω Z I hrD h z T w∇G Ge ) − w ∇G Gh )) ∇δG Gh i dS − Gne − G T G (G G (G hn h(w ∂I Ωh   Z r T z B w∇G Ge ) JδG Gne − GT Ghn K dS + G (G hn h s ∂I Ωh ∪∂D Ωh  Z r T z B e w G G Ghn K dS ∀δG Gh ∈ Xk . JδG − Gne − G T (w (G ) − w (G )) G G h ∇G ∇G hn h s ∂I Ωh (70)

Therefore, by applying the Taylor’s expansion (60-66) to the first, second, third, and fifth terms of Eq. (70), see C.2 for details, the latter expression is rewritten

22

Lina Homsi, Ludovic Noels +

as finding G h ∈ Xk such that: Ge ; G e − G h , δG Gh ) + B(G Ge ; G e − G h , δG Gh ) = N (G Ge , G h ; δG Gh ) A(G

Gh ∈ Xk . ∀δG (71)

ω ∈ X the In this last expression, we have defined for given ψ ∈ X+ , ω ∈ X and δω following forms: Z Z r z ψ ; ω , δω ω) = ω Tw ∇ψ ψ ω ωnT hw w∇ψ ψ ) ∇ω ω i dS+ A(ψ ∇δω (ψ )∇ω dΩ + δω ψ ψ (ψ Ωh ∂I Ωh ∪∂D Ωh Z r z w∇ψ ψ )∇δω ω i dS+ ω nT hw ψ (ψ ∂I Ωh ∪∂D Ωh  Z r z B ψ ) Jδω ωn K dS , w ∇ψ ω nT ψ (ψ hs ∂I Ωh ∪∂D Ωh (72) Z ψ ; ω , δω ω) = B(ψ

ω T (w wψ (ψ ψ , ∇ψ ψ )ω ω ) dΩ + ∇δω

ω Td ∇ψ ψ )∇ω ω dΩ δω ψ (ψ Ωh Z r z ωnT hw wψ (ψ ψ , ∇ψ ψ )ω ω i dS + ω Tdψ (ψ ψ , ∇ψ ψ )ω ω dΩ+ δω δω

Ωh

Z + Z

Z

∂I Ωh ∪∂D Ωh

r

Ωh

z

ψ ) hδω ω i dS , ω nTd T ψ (ψ ∇ψ

∂I Ωh ∪∂D Ωh

(73) ψ ; ., .) and the form which gather the terms such that for fixed ψ , the form A(ψ ψ ; ., .) are bi-linear. The non-linear terms have been gathered in N (G Ge , G h ; δG Gh ) B(ψ as Z ¯ w (G Gh )T (R Ge − G h , ∇G Ge − ∇G Gh ))dΩ Ge , G h ; δG Gh ) = (∇δG N (G Ωh Z r z

¯ w (G GT Ge − G h , ∇G Ge − ∇G Gh ) dS + δG R hn ∂ Ω ∪∂ Ω Z I hrD h z T w∇G Ge ) − w ∇G Gh )) ∇δG Gh i dS + Gne − G T G (G G (G hn h(w ∂I Ωh   Z r T z B (74) e e T w∇G G ) − w ∇G Gh )) JδG Ghn K dS (w + Gn − G hn G (G G (G h s ∂ Ω Z I h r z e T G − G h )To¯G Gh )(G Ge − G h ) hδG Ghn i dS + (G (G ∂I Ωh ∪∂D Ωh Z ¯ Ge − G h , ∇G GT Ge − ∇G Gh )dΩ + δG h Rd (G Ωh

= I1 + I2 + I3 + I4 + I5 + I6 . Compared with the fixed form from Gudi et al. [23] for non-linear elliptic probψ ; ., .) in lems, the formulations A and B are similar, except the last term of B(ψ ψ ) appears in the J·K operator instead of the h·i operator. Nevertheless, which d ∇ψ ψ (ψ this term becomes identical with the one in Gudi et al. [23] for fixed ψ . However the N is different in the fifth and sixth terms, i.e. I5 and I6 , so they will require

Title Suppressed Due to Excessive Length

23

a different treatment. Therefore in the following, we report the methodology developed by Gudi et al. [23] without demonstration, except when these two terms I5 and I6 explicitly appear, in which case the related demonstration is reported in Appendix. 4.5 Solution uniqueness Let us first define the mesh dependent norms, which will be considered in the following analysis, for G ∈ X |k G k|2∗ =

X

Gk2L2 (Ω e ) + k∇G

X e

e

|k G k|2 =

X

Gk2H1 (Ω e ) + kG

X e

e

2 h−1 s k JG n K kL2 (∂Ω e ) ,

(75)

2 h−1 s k JGn K kL2 (∂Ω e ) ,

(76)

and |k G k|21 =

X e

Gk2H1 (Ω e ) + kG

X

hs k G k2H1 (∂Ω e ) +

X e

e

2 h−1 s k JGn K kL2 (∂Ω e ) . (77) +

e

Let us first assume η = IhG − G ∈ X, with IhG ∈ Xk the interpolant of G e + in Xk . The last relation (71) thus becomes Ge ; IhG − G h , δG Gh ) + B(G Ge ; IhG − G h , δG Gh ) = A(G Ge ; η , δG Gh ) + B(G Ge ; η , δG Gh ) A(G Ge , G h ; δG Gh ) + N (G

Gh ∈ Xk . ∀δG (78)

Now in order to prove the existence of a solution G h of the problem stated by Eq. (70), which corresponds to the DG finite element discretization (48), we state the + + problem in the fixed point formulation and we define a map Sh : Xk → Xk as k+ k+ y) = Gy ∈ X , such that follows: for a given y ∈ X , find Sh (y Ge ; IhG − Gy , δG Gh ) + B(G Ge ; IhG − Gy , δG Gh ) = A(G Ge ; η , δG h ) + B(G Ge ; η , δG Gh ) A(G Ge , y ; δG Gh ) + N (G

Gh ∈ Xk . ∀δG (79)

The existence of a unique solution G h of the discrete problem (48) is equivalent to the existence of a fixed point of the map Sh , see [23]. For the subsequent analysis, we denote by Ck , a positive generic constant which is independent of the mesh size, but which does depend on the polynomial approximation degree k. We also use several Lemmata reported in C.3. Lemma 1 (Lower bound) For B larger than a constant, which depends on the polynomial approximation, there exist two constants C1k and C2k , such that Ge ; δG Gh , δG Gh ) + B(G Ge ; δG Gh , δG Gh ) ≥ C1k |k δG Gh k|2∗ −C2k k δG Gh k2L2 (Ω) A(G Gh ∈ Xk , ∀δG Ge ; δG Gh , δG Gh ) + B(G Ge ; δG Gh , δG Gh ) ≥ C1k |k δG Gh k|2 −C2k k δG Gh k2L2 (Ω) A(G Gh ∈ Xk . ∀δG

(80)

(81)

24

Lina Homsi, Ludovic Noels

Using the bounds (58) and (67), the Cauchy-Schwartz’ inequality, the trace inequality on the finite element space, see Lemma 7, the trace inequality, see Lemma 6, and the inverse inequality, see Lemma 8, the ξ-inequality –ξ > 0 : |ab| ≤ ξ 2 1 2 4 a + ξ b , as in Wheeler et al. [69] and Prudhomme et al. [52] analyzes, yields to prove this Lemma. The two constants C1k , C2k are independent of the mesh size, but require the stability parameter B to be larger than a constant, which depends on the polynomial degree, B > C k , to be positive . Lemma 2 (Upper bound) There exist C > 0 and Ck > 0 such that Ge ; m , δG G) + B(G Ge ; m , δG G) | ≤ C |k m k|1 |k δG G k|1 ∀m m , δG G ∈ X, | A(G

(82)

Ge ; m , δG Gh ) + B(G Ge ; m , δG Gh ) | ≤ C k |k m k|1 |k δG Gh k| ∀m m ∈ X , δG Gh ∈ Xk , | A(G (83) Ge ; m h , δG Gh ) + B(G Ge ; m h , δG Gh ) | ≤ C k |k m h k| |k δG Gh k| ∀m mh , δG Gh ∈ Xk . | A(G (84) Ge ; m, Applying the H¨ older’s inequality and the bound (67) on each term of A(G G)+B(G Ge ; m , δG G), then making use of the Cauchy-Schwartz’ inequality, lead to δG the relation (82). Finally, Eqs. (83) and (84) are easily obtained from the relation between the energy norms on the finite element space, Lemma 9. Using Lemma 1 and Lemma 2, the stability of the method is demonstrated through the following Lemmata. Lemma 3 (Auxiliary problem) We consider the following auxiliary problem, d+2 with φ ∈ L2 (Ω) : w∇G Ge )∇ψ ψ + wG (G Ge , ∇G Ge )ψ ψ ) + d ∇G Ge )∇ψ ψ + dG (G Ge , ∇G Ge )ψ ψ = φ on Ω, −∇T (w G (G G (G ψ = 0 on ∂Ω. (85) Assuming regular ellipticity of the operators and that wG and dG satisfy the weak  d minimum principle [21, Theorem 8.3], there is a unique solution ψ ∈ H2 (Ω) × H2 (Ω) × H2 (Ω) to the problem stated by Eq. (85) satisfying the elliptic property k ψ kH2 (Ωh ) ≤ C k φ kL2 (Ωh ) .

(86)

 d Moreover, for a given ϕ ∈ L2 (Ωh ) × L2 (Ωh ) × L2 (Ωh ) there exists a unique φ h ∈ Xk such that XZ Ge ; δG Gh , φ h ) + B(G Ge ; δG Gh , φ h ) = Gh dΩ ∀δG G h ∈ Xk , A(G ϕ T δG (87) e

Ωe

and there is a constant Ck such that : |k φ h k|≤ Ck k ϕ kL2 (Ωh ) .

(88)

Title Suppressed Due to Excessive Length

25

The proof of the first part is given in [21], by combining [21, Theorem 8.3] to [21, Lemma 9.17]. The second part was demonstrated by Gudi et al. [23]: The use of Lemma 1 Gh = φ h allows bounding |k φ h k| in terms of k ϕ kL2 (Ωh ) and and Eq. (87) with δG k φ h kL2 (Ωh ) ; The term k φ h kL2 (Ωh ) is then evaluated by using φ = φ h ∈ Xk in Eq. (85), multiplying the result by φ h and integrating it by parts on Ωh , resulting into k Ge ; ψ , φ h )+B(G Ge ; ψ , φ h ); Inserting the interpolant Ihψ in these last φ h k2L2 (Ωh ) = A(G 2 e G ; ψ −Ihψ , φ h )+B(G Ge ; ψ −Ihψ , φ h )+A(G Ge ; Ihψ , φ h )+ terms, i.e. k φ h kL2 (Ωh ) = A(G Ge ; Ihψ , φh ), bounding the last two terms by considering δG Gh = Ihψ in Eq. (87), B(G and bounding the first two terms by making successive use of Lemma 2 and of the energy bound of the interpolant error, Lemma 10, and using the regular ellipticity Eq. (86) lead to the bound k φ h kL2 (Ωh ) ≤ Ck k ϕ kL2 (Ωh ) , and thus to the proof of the solution uniqueness. Theorem 1 (Solution uniqueness to the problem (79)) The solution Gy to + the problem stated by Eq. (79) is unique for a given y ∈ Xk with Sh (yy) = Gy . +

In order to prove that the solution Gy is unique for a given y ∈ Xk , let us assume that there are two distinct solutions Gy 1 , Gy 2 , such that Ge ; IhG − Gy 1 , δG Gh ) + B(G Ge ; IhG − Gy 1 , δG Gh ) A(G Ge ; IhG − Gy 2 , δG Gh ) + B(G Ge ; IhG − Gy 2 , δG Gh ) ∀ δG Gh ∈ Xk . = A(G

(89)

For fixed G e , since A and B are bi-linear, this last relation thus becomes Ge ; Gy 1 − Gy 2 , δG Gh ) + B(G Ge ; Gy 1 − Gy 2 , δG Gh ) = 0 ∀ δG Gh ∈ Xk . A(G

(90)

Gh = Gy 1 −G Gy 2 ∈ Xk in Lemma 3 consists in stating that there Considering ϕ = δG is a unique Φ h ∈ Xk solution of the problem Eq. (87), with Ge ; Gy 1 − Gy 2 , Φ h ) + B(G Ge ; Gy 1 − Gy 2 , Φ h ) =k Gy 1 − Gy 2 k2L2 (Ωh ) , A(G

(91)

Gh = Φ h in Eq. (90) and with |k Φ h k|≤ Ck k Gy 1 − Gy 2 kL2 (Ωh ) . The choice δG y) = Gy is leads to k Gy 1 − Gy 2 kL2 (Ωh ) = 0, demonstrating that the solution Sh (y unique. + We will now show that Sh maps a ball Oσ (IhG ) ⊂ Xk into itself and is continuous in the ball. To this end, we define the ball Oσ with radius σ and centered at the interpolant IhG of G e as n o + Oσ (IhG ) = y ∈ Xk such that |k IhG − y k|1 ≤ σ , (92) 1 |k IhG − G e k|1 , 0 0 and detF F p(α) = 1 , F = F e(α) · F p(α) with detF

(123)

where F e(α) is the elastic distortion and F p(α) is the inelastic distortion with F p(α) (X, 0) = I . In these equations we have considered the possibility to account for several mechanisms α = 1, 2, 3. Moreover, the elastic decomposition of the deformation gradient can be written as F e(α) = R e(α) · U e(α) ,

(124)

with the elastic right and left Cauchy-Green strain tensors respectively equal to 2

C e(α) = U e(α) = F e(α)T · F e(α) and B e(α) = F e(α) · F e(α)T .

(125)

A.2.2 Elasto-visco-plasticity The material may be idealized to be isotropic. Accordingly, all constitutive functions are presumed to be isotropic in character. Let us assume that the free energy has the separable form ψR =

X

ψ (α) (ΦC e(α) , T ),

(126)

α

where Φce(α) represents a list of the principle invariants of C e(α) and T is the temperature. P The Cauchy stress is decomposed in terms of the mechanisms σ = (α) σ (α) with σ (α) =

1 e(α) e(α) e(α)T 1 F S (α) F T = F S F , J J

(127)

where S e(α) is the symmetric elastic second Piola-Kirchhoff stress S e(α) = 2

∂ψ (α) (Φc e(α) , T ) C e(α) ∂C

.

(128)

Moreover, the first Piola-Kirchhoff stress tensor can be computed following P (α) = J σ (α)F −T = F e(α) S e(α) F p(α)−T = F F p(α)−1 S e(α) F p(α)−T .

(129)

46

Lina Homsi, Ludovic Noels

The driving stress of the plastic flow is the symmetric Mandel stress, which is defined as M e(α) = J R e(α)Tσ (α)R e(α) = U e(α)S e(α)U e(α) = C e(α)S e(α) ,

(130)

where M e(α) is the elastic Mandel stress, R e(α) is the rotation matrix, and where it has been assumed that C e(α) and S e(α) permute. The corresponding equivalent shear stress is given by 1 e(α) M 0 |, τ¯(α) = √ |M 2 e(α)

where M 0

(131)

is the deviatoric part of the Mandel stress e(α)

M0

1 M e(α) , = M e(α) + pII , with p = − trM 3

(132)

q e(α) e(α) e(α) M 0 |= M 0 is the norm of the deviatoric part of the Mandel stress. and |M : M0 In order to account for the major strain-hardening and softening characteristics of polymeric materials observed during visco-plastic deformation, Srivastava et al. [59] have introduced macroscopic internal variables ξ (α) to represent important aspects of the microstructural resistance to plastic flow. The plastic flow follows p(α) F˙ = D p(α)F p(α) ,

(133)

where each F p(α) is to be regarded as an internal variable part of ξ (α) , and which is defined as a solution of the differential equation e(α)

M0 ), (134) 2¯ τα √ D p(α) | is the equivalent plastic shear where D p is the plastic stretching tensor, and ˙p(α) = 2|D strain rate. C e(α) , B p(α) , ξ (α) , T ) a list of constitutive variables, Therefore for given τ¯(α) and Λ (α) = (C the equivalent plastic shear strain rate ˙p(α) is obtained by solving a scalar strength relation such as Λ(α) , ˙p(α) ), τ¯(α) = Υ (α) (Λ (135) D p(α) = ˙p(α) (

Λ(α) , ˙p(α) ) is an isotropic function of its arguments. where the strength function Υ (α) (Λ

A.2.3 Partial differential governing equations In order to complete Eq. (7), y, the internal energy per unit mass, is defined as y = cv T , where the volumetric heat capacity per unit mass is a function of the glass transition temperature, and is defined as follows  c0 − c1 (T − Tg ) if T ≤ Tg (136) cv = c0 if T > Tg . Moreover, F¯ , the body source of heat, is expressed as F¯ = Qr +

X

τ¯(α) ˙p(α) + T

α

∂ 2 ψ e(α) e(α) : C˙ , C e(α) ∂T ∂C

(137)

where Qr is the scalar heat supply measured per unit reference volume and the last term of the right hand side is the thermo-elastic damping term which is neglected. Instead it is assumed that only a fraction v of the rate of the plastic dissipation contributes to the temperature change X F¯ = Qr + v τ¯(α) ˙p(α) , (138) α

where 0 ≤ v ≤ 1.

Title Suppressed Due to Excessive Length

47

The shear strain rate √ glass transition in amorphous polymers depends on the equivalent D 0 | to which the material is subjected, where D 0 = D − 31 trD D I denotes the total ˙ = 2|D deviatoric stretching tensor, with D=

1 ˙ −1 T (F F + F −TF˙ ) . 2

Eventually, the glass transition Tg is calculated from the following expression   Tr if ˙ ≤ r , ˙ Tg =  Tr + n log ( ) if ˙ > r , r

(139)

(140)

where Tr is the reference glass transition temperature at low strain rate, ˙ is the shear strain rate, and r is the reference strain rate.

A.2.4 The first micromechanism (α = 1) The first micromechanism (α = 1) represents an elastic resistance due to intermolecular energetic bond-stretching and a dissipation due to the thermally-activated plastic flow following chain segment rotation and relative slippage of the polymer chains between neighboring crosslinkage points. The following simple generalization of the classical strain energy function of infinitesimal isotropic elasticity is considered, and uses a logarithmic measure E e(1) =

1 ln C e(1) , 2

(141)

of the finite elastic strain [4]. The form of the elastic free energy is thus defined as e(1) 2

E0 ψ e(1) = G|E

| +

2   1  E e(1) − 3K trE E e(1) αth (T − T0 ) + f˜(T ) , K trE 2

(142)

where the deviatoric part of the logarithmic strain is denoted by E e0 , f˜(T ) is an entropy contribution to the free energy related to the temperature dependent specific heat of the material, and where the temperature dependent parameters G(T ), K(T ), αth (T ) are respectively the shear modulus, bulk modulus, and the coefficient of thermal expansion. The Mandel stress is thus obtained from C e(1) M e(1) = 2C

E e(1) , T) ∂ψ e(1) (E C ∂C

e(1)

=

E e(1) , T ) ∂ψ e(1) (E , E e(1) ∂E

(143)

if C e(1) and M e(1) permute. It should be noted that in this work E e(1) is computed by using a Taylor series approximation of Eq. (141), and not through the eigenvalue decomposition. e(1) e(1) e(1) E 0 |= E 0 : E 0 Substituting Eq. (142) in Eq. (143), as |E one can get directly M e(1) as e(1)

E0 M e(1) = 2GE

  E e(1) I − 3Kαth (T − T0 )II . + K trE

(144)

The thermal expansion is taken to have a bilinear temperature dependence, with the slope αth = αr above the glass transition temperature and the slope αth = αgl below it. Moreover, the evaluation equation for F p(1) follows Eqs. (133-134) with the thermallyactivated relation for the equivalent plastic strain rate following  if τ e(1) ≤ 0, 0 p(1) e(1) 0 (145) ˙ = Q(T ) τ V 1/m(1) 1  .0(1) exp (− ) exp (− )[sinh( )] if τ e(1) > 0, ξ KB T 2 KB T . (1)

where ˙p(1) is the plastic strain rate, the parameter 0 is a pre-exponential factor with units of 1/time, KB is Boltzmann’s constant, V 0 is an activation volume, m(1) is the sensitive

48

Lina Homsi, Ludovic Noels

parameter for the strain rate, Q(T ) is the temperature dependent activation energy, and τ e(1) denotes a net shear stress for the thermally activated flow τ e(1) = τ¯(1) − (Sa + Sb + αp p¯),

(146)

1 with αp > 0 a parameter introduced to account for the pressure sensitivity. The term exp (− ) ξ in Eq. (145) represents a concentration of flow defects, with  ξ=

ξgl ξgl + d(T − Tg )

if if

T ≤ Tg , T > Tg .

(147)

For the first micromechanism, besides the plastic strain gradient, the list ξ 1 = (ϕ, Sa , Sb ) of internal variables consists of three positive scalars, where the variable ϕ ≥ 0 and Sa ≥ 0 are introduced to model the yield peak which is observed in the stress-strain response of glassy polymers and Sb ≥ 0 is introduced to model the isotropic hardening at high strain. The evolution equations of S˙ a and ϕ˙ are governed by S˙ a = ha [b(ϕ∗ − ϕ) − Sa ] ˙p(1) ϕ˙ = g(ϕ∗ − ϕ)˙p(1) with

ϕ∗

with initial value Sa = Sa0 ,

with initial value ϕ = ϕ0 ,

(148) (149)

as

   p(1) ˙ T r   ) + hg ( )s  z (1 − Tg r ϕ∗ (˙p(1) , T ) =  ˙p(1) s   z hg ( ) r

if T ≤ Tg and ˙p(1) > 0 , (150) if T > Tg and ˙p(1) > 0 ,

which represents the temperature and strain rate dependency of ϕ, where z, r, hg , and s are taken to be constants. In particular we have introduced hg to get a small value for ϕ∗ when T > Tg , instead of 0 in order to improve the convergence of the numerical model. Then the evolution of Sb is governed by ¯ − 1)a , λ ¯= Sb = Sb0 + Hb (λ

p C /3, trC

(151)

¯ is an effective stretch which increases or decreases as the overall stretch increases or where λ decreases, and where Hb (T ) is a temperature dependent hardening parameter. The temperature evolutions of Hb , Q(T ), G(T ), and of the Poisson ratio ν(T ) follow a law ·(T ) =

1 1 1 (·gl + ·r ) − (·gl − ·r ) tanh( (T − Tg )) − L· (T − Tg ) , 2 2 ∆

(152)

where ·gl and ·r are the values in glassy and rubbery regions, and where L· represents the slope of the temperature variation of ·, and takes the value of L· = L·gl if T ≤ Tg and L· = L·r if T > Tg . The temperature dependence of the bulk modulus K(T ) is then obtained 2(1+ν(T )) by using the standard relation for isotropic materials K(T ) = G(T ) 3(1−2 ν(T )) .

A.2.5 The second micromechanism (α = 2) The second micromechanism (α = 2) represents the molecular chains between mechanical crosslinks. At temperatures below Tg the polymer exhibits a significant amount of mechanical crosslinking which disintegrates when the temperature is increased above Tg . Only deviatoric contributions are considered in the free energy function ¯ ψ (2) = ψ¯(2) (C

e(2)

, T) = −

¯ e(2) − 3 trC 1 (2) (2) ), µ Im ln(1 − (2) 2 Im

(153)

Title Suppressed Due to Excessive Length

49

¯ e(2) = F¯ e(2)TF ¯ e(2) = J− 23 C e(2) denotes the distrotional (or volume preserving) right where C (2) Cauchy strain tensor, the parameter Im is taken to be temperature constant, and where µ(2) is the rubbery shear modulus, which follows µ(2) = µg exp(−N (T − Tg )) ,

(154)

µ(2)

with µg the value of at the glass transition temperature, and N a parameter that represents the slope of temperature variation on a logarithmic scale. The corresponding Mandel stress is evaluated from Eq. (130) and Eq. (153) as M e(2) = C e(2) 2

¯ e(2) − 3 ∂ ψ¯(2) trC ¯ e(2) , )−1C = µ(2) (1 − 0 (2) e(2) C ∂C Im

(155)

e(2)

¯ e(2) − 1 trC ¯ e(2)I is the deviatoric part of C e(2) the right Cauchy Green tensor. =C 3 e(2) ¯ ¯ e(2) permute as well. Clearly, as C and C e(2) permute, M e(2) and C For the second mechanism, the equivalent plastic strain rate follows ¯ where C 0

(2)

˙p(2) = ˙0 (

1 τ¯(2) ) m(2) , S (2)

(156)

(2)

where ˙0 is a reference plastic shear strain rate, m(2) is the positive valued strain rate sensitivity parameter, and S (2) is a temperature dependent parameter which follows (152).

A.2.6 The third micromechanism (α = 3) The third micromechanism (α = 3) introduces the molecular chains between chemical crosslinks and represents the resistance due to changes in the free energy upon stretching of the molecular chains between the crosslinks. ¯ = F¯ TF¯ = J − 32 C , and is given by The free energy is a function of the deviatoric tensor C a deviatoric form ¯ −3 trC (3) ¯ ) = − 1 µ(3) Im ψ (3) = ψ¯(3) (C ln(1 − ), (157) (3) 2 Im (3)

where the material constants µ(3) > 0 and Im > 0 are assumed to be temperature-independent. The free energy (157) yields the corresponding second Piola stress S (3) as S (3) = 2

¯ ¯ −3 2 1 trC ∂ ψ¯(3) ∂C ¯ −1 ] . ¯ )C )−1 [II − (trC = J − 3 µ(3) (1 − (3) ¯ : ∂C C 3 ∂C Im

(158)

A.2.7 Finite increment form of the Shape Memory Polymer constitutive law The constitutive laws are formulated in a finite strain setting and solved following the predictorcorrector scheme during the time interval [tn ; tn+1 ], where we use the subscript n for the previous time tn and n + 1 for the current time tn+1 . The formulation can be summarized as follows: – Prediction step: The plastic deformation gradient is initialized to the value at the previous p(α) p(α) step F (pr) = F n , and the elastic deformation follows e(α)

p(α)−1

F n+1 = F n+1F n

.

(159)

– Correction step: In this step we solve the system of equations that has been presented for each mechanism. To extract the plastic increment using the evaluation equation of the plastic deformation gradient during the time step between the configurations n and n + 1, we consider p(α) p(α) Dp(α) )F Fn , F n+1 = exp(∆D (160) with ! e(α) M e(α) p(α) p(α) M D p(α) = (n+1 − n ) (α) = ∆p(α) . (161) ∆D 2¯ τ 2¯ τ (α) More details about the predictor-corrector algorithm and the stiffness computation can be found in [26].

50

Lina Homsi, Ludovic Noels

B The finite element formulation Using the interpolations (33-34), the gradients are readily obtained by: a Ma , ∇0 M h = ∇0 N M

∇0 uh = u a ⊗ ∇0 Nua , a

where ∇0 Nua

a ∇0 N M

∇0 Nua

(162)

a

M , Mh = u ⊗ δM (163) , ∇0 δM ∇0 δuh = δu   a ∇ N 0 0 fV a = and ∇0NM are the gradients of the shape functions at 0 ∇0 Nfa T

node a.

B.1 Nodal forces The expressions of the nodal forces (36) are obtained by substituting the interpolations (33-34) and (162-163) in the weak formulation (31). First the mechanical contribution reads XZ XZ ¯ dS0 − Fuaext = Nua T (¯ u ⊗ N : H ) · ∇0 Nua dS0 (∂N Ω0 )s

s

+

¯ ⊗N : u (∂D Ω0 )s

s

Z

(∂D Ω0 )s

s



XZ

HB hs



· N Nua dS0

 ¯ 0 )M ¯ 0 · N N a dS0 , ¯ )M ¯ − Y (M Y (M u



(164)

∂D Ω0h

Fuaint = Fua± I =

XZ

F h , M h ) · ∇0 Nua dΩ0 , and P (F

e Ω0 e a± Fu I1 + Fua± I2

(165)

+ Fua± I3 ,

(166)

with the three mechanical contributions to the interface forces related to the degrees of freedom of the nodes a± on each side of the interface elements reading1 XZ Fua± P (F F h , Mh )i · N − dS0 , (±Nua± ) hP (167) I1 = (∂I Ω0 )s

s

Fua± I2 Fua± I3

Z 1X = Juh K ⊗ N − : H ± · ∇0 Nua± dS0 , 2 s (∂I Ω0 )s   XZ HB · N − (±Nua± )dS0 . = (Juh K ⊗ N − ) : s h s (∂ Ω ) 0 I s

Secondly, the electro-thermal contributions read XZ XZ a a ¯ FM = NM J dS0 − ext (∂N Ω0 )s

s

+

a± FM = I

XZ

T

a ¯T ¯) B M ¯ N dS0 , NM NM Z 0 (F F h, M hs (∂D Ω0 )s XZ aT F h , M h , ∇M Mh )dΩ0 + ∇0 N M J (F

e a± FM I1

e Ω0

e

a± a± + FM + FM , I2

I3

(169)

a ¯ )M ¯ N dS0 F h, M ∇0 N M Z 0 (F

XZ s

a FM = int

s

(∂D Ω0 )s

(168)

e Ω0

(170) T

a NM I i dΩ0 , and

(171) (172)

1 The contributions on ∂ Ω D 0h can be directly deduced by removing the factor (1/2) accord¯) F h, M ingly to the definition of the average flux on the Dirichlet boundary and by using Z 0 (F a± Z F M F instead of 0 (F h , h ). However, there is one more additional term in u I1 in the Dirichlet P R Y (M M)M M − Y (M M0 )M M0 ) · N − dS0 . boundary, which is s (∂ Ω )s (−Nua ) (Y D

0

Title Suppressed Due to Excessive Length

51

with the three electric contributions to the interface forces1  T XZ a± a± ¯− NM J(F F h , M h , ∇M Mh )i dS0 , FM = (±N ) N hJ M I1

XZ

1 2

a± FM = I2

I3

 q y a±T ± F h , M h ) M hN dS0 , ∇0 N M Z 0 (F

(174)

  T  y B q a± ¯− NM F h , Mh ) (±N ) N Z 0 (F M hN dS0 . M hs (∂I Ω0 )s

(175)

(∂I Ω0

s

XZ

a± FM =

(173)

(∂I Ω0 )s

s

s

)s

B.2 Tangent operators ab = In order to derive the tangent matrix KG



Kuu KuM KMu KMM



Fa ∂F ext Gb ∂G

u ∆u M ∆M



Fa ∂F int Gb ∂G



 =−



Fa ∂F I , Gb ∂G

u, M ) Ru (u u, M ) RM (u

the system (38) is rewritten  .

(176)

The stiffness matrix has been decomposed into four sub-matrices as shown in Eq. (176) with respect to the discretization of the five independent field variables (d for displacement u , and 2 for M (for fV , and fT )), and can be obtained in a straighforward way from the internal forces, see details in [26].

C Derivation of the numerical properties C.1 Taylor’s remainders V − G ), ∇Q Qt = The remainder terms of Eqs. (60-61) are obtained by defining V t = G + t(V G + t(∇Q Q − ∇G G). They can thus be evaluated by ∇G 1

Z ¯ G (G G, ∇G G) = w

Vt , ∇Q Qt )dt, wG (V

1

Z ¯ ∇G G) = w G (G

0

Vt )dt , w ∇G G (V

(177)

0

and by 1

Z ¯ GG (V V, ∇Q Q) = w

wGG (V Vt , ∇Q Qt )dt , w ¯ G ∇G V) = (1 − t)w G (V

1

Z

0

wG ∇G Vt )dt , (1 − t)w G (V

(178)

0

G, ∇G G) = vGG (G G)∇G G, and wG ∇G G) = vG (G G) of w (G G, ∇G G), with the partial derivatives wGG (G G (G G) = 0. since w ∇G G∇G G (G The remainder terms of Eqs. (62-63) read ¯G (G G, ∇G G) = d

1

Z

Vt , ∇Q Qt )dt, dG (V

¯∇G G) = d G (G

1

Z

0

Vt )dt , d ∇G G (V

(179)

0

and ¯GG (V V, ∇Q Q) = d

Z

1

¯G∇G dGG (V Vt , ∇Q Qt )dt , d V) = (1 − t)d G (V

0

1

Z

dG∇G Vt )dt . (1 − t)d G (V

(180)

0

Finally, the remainder terms of Eqs. (64-65) read 1

Z ¯G (G G) = p 0

Vt )dt, pG (V

1

Z ¯GG (V V) = p 0

pGG (V Vt )dt . (1 − t)p

(181)

52

Lina Homsi, Ludovic Noels

C.2 Application of the Taylor’s expansion The first term of Eq. (70) is rewritten, using the Taylor’s expansion defined in Eq. (60), as Z

Gh )T (w w(G Ge , ∇G Ge ) − w (G Gh , ∇G Gh ))dΩ = (∇δG

Z

Ωh

Gh )T (w wG (G Ge , ∇G Ge )(G Ge − G h ))dΩ (∇δG Ωh

Z

Gh )T (w w∇G Ge )(∇G Ge − ∇G Gh ))dΩ − (∇δG G (G

+

Z

Ωh

¯ w (G Gh )T (R Ge − G h , ∇G Ge − ∇G Gh ))dΩ . (∇δG Ωh

(182) Similarly, the second term of Eq. (70) is rewritten, using the Taylor’s expansion defined in Eq. (62), as Z Ωh

GT Ge )∇G Ge − o˜(G Gh )∇G Gh ) dΩ = δG o (G h (˜

Z

GT Ge , ∇G Ge )(G Ge − G h ))dΩ + δG h dG (G

Z

Z = Ωh

Z − Ωh

Ωh

Ωh

GT d(G Ge , ∇G Ge ) − d (G Gh , ∇G Gh )) dΩ δG h (d GT Ge )(∇G Ge − ∇G Gh ))dΩ δG G (G h d ∇G

(183)

¯ Ge − G h , ∇G GT Ge − ∇G Gh )dΩ . δG h Rd (G

Likewise, the third term is rewritten, using the Taylor’s expansion defined in Eq. (60), as Z

r

z GT w(G Ge , ∇G Ge ) − w (G Gh , ∇G Gh )i dS = δG hn hw

r

z GT wG (G Ge , ∇G Ge )(G Ge − G h )i dS+ δG hn hw

r

z GT w∇G Ge )(∇G Ge − ∇G Gh )i dS− δG G (G hn hw

r

GT δG hn

∂I Ωh ∪∂D Ωh

Z ∂I Ωh ∪∂D Ωh

Z ∂I Ωh ∪∂D Ωh

Z ∂I Ωh ∪∂D Ωh

(184)

z

¯ w (G Ge − G h , ∇G Ge − ∇G Gh ) dS. R

G) = G To T (G G) and using The fifth term of Eq. (70) is developed by using the definition of p T (G the Taylor’s expansion defined in Eq. (64), leading to Z − ∂I Ωh ∪∂D Ωh

Z − ∂I Ωh ∪∂D Ωh

r T z T Ge ) − G T Gh ) hδG Ghn i dS = G e o T (G h o (G Z r z T T Ge − G T pG Ge ) hδG Ghn i dS + (G (G h )p

∂I Ωh ∪∂D Ωh

r z e ¯ T (G Ghn i dS. R p G − G h ) hδG (185)

T = −o oT (G G), using Eq. (66), this last term is also rewritten as However, as pG

Z − ∂I Ωh ∪∂D Ωh

Z ∂I Ωh ∪∂D Ωh

r T z T Ge ) − G T Gh ) hδG Ghn i dS = G e o T (G h o (G

r z T Ge − G T oT (G Ge ) hδG Ghn i dS (G h )o

Z

r

− ∂I Ωh ∪∂D Ωh

z T Ge − G h )To¯G Gh )(G Ge − G h ) hδG Ghn i dS. (G (G

(186)

Title Suppressed Due to Excessive Length

53

G0 )δG Gn = GnTo˜T (G G0 )δG G2 , and Eq. (186) Finally, using the definition of the ˜· operator G To T (G is rewritten as Z r T z T Ge ) − G T Gh ) hδG Ghn i dS = − G e o T (G h o (G ∂I Ωh ∪∂D Ωh

Z

r

∂I Ωh ∪∂D Ωh

z T Ge ) hδG Gh i dS Gne − G T o T (G (G hn )˜

Z

r

− ∂I Ωh ∪∂D Ωh

(187)

z T Gh )(G Ge − G h ) hδG Ghn i dS. Ge − G h )To¯G (G (G

C.3 General properties of the finite element method and Hilbert spaces The following properties will be used in the numerical properties derivation. Lemma5 (Interpolation inequality) For all G ∈ (Hs (Ω e ))n there exists a sequence G h ∈ n Pk (Ω e ) and a positive constant CkD depending on s and k but independent of G and hs , such that 1. for any 0 ≤ n ≤ s k G − G h kHn (Ω e ) ≤ CkD hµ−n k G kHs (Ω e ) , s 2. for any 0 ≤ n ≤ s − 1 +

2 r 2 µ−n−1+ r

k k G − G h kWn e ≤ CD hs r (Ω )

3. for any s > n +

(188)

(189)

k G kHs (Ω e ) , if d = 2,

1 2 1 µ−n− 2

k G − G h kHn (∂Ω e ) ≤ CkD hs

(190)

k G kHs (Ω e ) ,

where µ = min {s, k + 1}. The proof of the first and third properties can be found in [7], and the proof of the second property in the particular case of d = 2 can be found in [1, 2], see also the discussion by [23]. Remarks i) The approximation property in (2) is still valid for r = ∞, see [1]. G|Ω e ), which means ii) For G ∈ Xs , let us define the interpolant Ih G ∈ Xk by Ih G|Ω e = Gh (G IhG satisfies the interpolant inequality property provided in Lemma 5 on Ωh , see [28]. n Lemma 6 (Trace inequality) For all G ∈ Hs+1 (Ω e ) , there exists a positive constant CT , such that   1 s+1 k G krWs (∂Ω e ) ≤ CT k G krWs (Ω e ) + k G kr−1 k ∇ G k (191) 2 e s L (Ω ) , W2r−2 (Ω e ) r r hs where s = 0, 1 and r = 2, 4, or in other words   1 G kL2 (Ω e ) , k G k2L2 (Ω e ) + k G kL2 (Ω e ) k ∇G k G k2L2 (∂Ω e ) ≤ CT hs   1 4 G kL2 (Ω e ) . k G kL4 (∂Ω e ) ≤ CT k G k4L4 (Ω e ) + k G k3L6 (Ω e ) k ∇G hs

(192)

The first equation, s = 0 and r = 2, is proved in [52], and the second one, r = 4 and s = 0, is proved in [22]. 2 3K 20 fT

G0 )δG Gn = GnTo˜T (G G0 )δG G = − 3K In this case G To T (G f α n− δux − 20 T th x fT αth n− z δuz

fT

3K 20 fT

fT αth n− y δuy −

54

Lina Homsi, Ludovic Noels

Lemma 7 (Trace inequality on the finite element space) For all G h ∈ Pk (Ω e ) exists a constant CkK > 0 depending on k, such that −1 2

k ∇l Gh kL2 (∂Ω e ) ≤ CkK hs

Gh k2 2 hs k∇G L (∂Ω e ) G h ∈(Pk (Ω e ))n k∇G Gh k2 2 e L (Ω )

where CkK = sup

polynomial approximation only with hs =

k ∇l Gh kL2 (Ω e )

l = 0, 1,

n

there

(193)

is a constant which depends on the degree of the

|Ω e | , |∂Ω e |

see [25] for more details.

n and r ≥ 2, there exists CkI > 0, such Lemma 8 (Inverse inequality) For G h ∈ Pk (Ω e ) that d−d 2

k G h kLr (Ω e ) ≤ CkI hsr

d−1 − d−1 r 2

k G h kLr (∂Ω e ) ≤ CkI hs Gh kL2 (Ω e ) ≤ k ∇G

CkI h−1 s

(194)

k G h kL2 (Ω e ) , k G h kL2 (∂Ω e ) ,

k G h kL2 (Ω e ) .

(195) (196)

The proof of these properties can be found in [12, Theorem 3.2.6]. Note that Eqs. (194-195) involve the space dimension d. Lemma 9 (Relation between energy norms on the finite element space) From [69], for G h ∈ Xk , there exists a positive constant Ck , depending on k, such that |k G h k|1 ≤ Ck |k G h k| .

(197)

The demonstration directly follows by bounding the extra terms

P

e

hs k G k2H1 (∂Ω e ) of the

norm defined by Eq. (77), in comparison to the norm defined by Eq. (76), using successively the trace inequality, Eq. (192), and the inverse inequality, Eq. (196), for the first term, and the trace inequality on the finite element space, Eq. (193), for the second term. Lemma 10 (Energy bound of interpolant error) Let G e ∈ Xs , s ≥ 2, and let IhG ∈ Xk , be its interpolant. Therefore, there is a constant Ck > 0 independent of hs , such that |k G e − IhG k|1 ≤ Ck hsµ−1 k G e kHs (Ωh ) ,

(198)

with µ = min {s, k + 1}. The proof follows from Lemma 5, Eq. (188), and Eq. (190), applied on the mesh dependent norm (77).

Ge , y ; δG Gh ) C.4 The bound of the non-linear term N (G C.4.1 Intermediate bounds Ge , y ; δG Gh ), we have recourse to the following intermediate bounds, To bound the terms of N (G which are derived for the particular case of d = 2. Gh ∈ Xk , η = Ge − Ih G ∈ X and Lemma 11 (Intermediate bounds) Let ξ = Ih G − y, δG ζ = ξ + η . The terms in ξ can be bounded by recourse to the trace inequality, see Lemma 6, and to the inverse inequality, see Lemma 8, while bounding the terms in η makes use of the trace inequality, see Lemma 6 and of the interpolation inequality, see Lemma 5 for d = 2. Using the definition of the ball (92-93) thus leads to the different contributions, see [23, 26] for details, !1 2 X 2 k ζ kL2 (Ω e ) ≤ C k σ ≤ C k hsµ−1−ε k G e kHs (Ωh ) , (199) e

!1

4

X e



k4L4 (Ω e )

−1 2

≤ C k hs

3 −ε µ− 2

σ ≤ C k hs

k G e kHs (Ωh ) ,

(200)

Title Suppressed Due to Excessive Length

55

!1 2

X

k G e kHs (Ωh ) , ≤ C k σ ≤ C k hµ−1−ε s

k ∇ζζ k2L2 (Ω e )

(201)

e

!1 4

X

−3 4

k ζ k4L4 (∂Ω e )

≤ C k hs

µ− 7 −ε 4

σ ≤ C k hs

k G e kHs (Ωh ) ,

(202)

k G e kHs (Ωh ) ,

(203)

e

!1

4

X

k Jζζ K

k4L4 (∂Ω e )

k ∇ζζ

k4L4 (∂Ω e )

e

1

!1

4

X

µ− 3 −ε 4

≤ C k hs4 σ ≤ C k hs −3 4

≤ C k hs

µ− 7 −ε 4

σ ≤ C k hs

k G e kHs (Ωh )

(204)

e

with µ = min {s, k + 1}. Moreover, using the inverse inequality, see Lemma 8, one has  −1 k δG Gh kH1 (Ω e ) , Gh kW1 (Ω e ) ≤ CIk hs 2 k δG 4

−1 2

≤ CIk hs

| δG Gh |W1 (Ω e ) 4

(205)

Gh |H1 (Ω e ) . | δG

C.4.2 Bounds of the different contributions Ge , y ; δG Gh ) follows from the argumentation reported in [23] and is nominated The bound of N (G by the term with the largest bound. Ge , y ; δG Gh ), defined in Eq. (74), can be expanded using Eq. (61) as The first term of N (G Z XZ ¯ w (ζζ , ∇ζζ )dΩ = G)T Gh )T (ζζ Tw ¯ GG (y y, ∇y y)ζζ )dΩ (∇δG I1 = R (∇δG h Ωh

+2

e

XZ Ωe

e

Ωe

(206)

  T G)T ¯ G ∇G y)∇ζζ dΩ = I11 + 2I12 . (∇δG G (y h ζ w

The two term of the right hand side of Eq. (206) are bounded by using the generalized H¨ older’s inequality, the generalized Cauchy-Schwartz’ inequality, the definition of Cy in Eq. (67), and the bounds (199, 200, 201, and 205) as X Gh kL4 (Ω e ) | I11 | ≤ Cy k ζ kL4 (Ω e ) k ζ kL2 (Ω e ) k ∇δG e

!1

!1

4

X

≤ Cy



!1

2

X

k4L4 (Ω e )



4

X

k2L2 (Ω e )

e

e

k

(207)

Gh k4L4 (Ω e ) ∇δG

e

Gh |H1 (Ωh ) k G e kHs (Ωh ) , ≤ C k Cy hµ−2−ε σ | δG s | I12 | ≤ Cy

X

Gh kL4 (Ω e ) k ζ kL4 (Ω e ) k ∇ζζ kL2 (Ω e ) k ∇δG

e

!1

!1

4

≤ Cy

X



!1

2

X

k4L4 (Ω e )

e

k ∇ζζ

k2L2 (Ω e )

e

4

X

k

Gh k4L4 (Ω e ) ∇δG

(208)

e

Gh |H1 (Ωh ) k G e kHs (Ωh ) . ≤ C k Cy hµ−2−ε σ | δG s Combining the above results leads to Gh |H1 (Ωh ) k G e kHs (Ωh ) . | I1 |≤ C k Cy hµ−2−ε σ | δG s

(209)

Ge

Gh ), defined in Eq. (74), becomes by using Eq. (61), The second term of N (G G , y ; δG Z r zD E GT ¯ GG (y y, ∇y y)ζζ dS I2 = δG ζ Tw hn ∂I Ωh ∪∂D Ωh

Z +2 ∂I Ωh ∪∂D Ωh

r

GT δG hn

zD E ¯ G ∇G y)∇ζζ dS = I21 + 2I22 . ζ Tw G (y

(210)

56

Lina Homsi, Ludovic Noels

The two terms of the right hand side of Eq. (210) are bounded by using the generalized H¨ older’s inequality, the generalized Cauchy-Schwartz’ inequality, the definition of Cy in Eq. (67), and the bounds (202, 204), yielding !1 !1 2 2 r z 1 X X T 2 G h−1 k δG k | I21 | ≤ Cy hs2 k ζ k4L4 (∂Ω e ) s hn L2 (∂Ω e ) e

e

(211)

!1

2

k

e

≤ C Cy k G

kHs (Ωh ) hµ−2−ε σ s

X

h−1 s

e

!1

!1

4

1

X

| I22 | ≤ Cy hs2

Ghn K k JδG 4

X

k ζ k4L4 (∂Ω e )

e

X

k ∇ζζ k4L4 (∂Ω e )

e

k2L2 (∂Ω e )

,

r z 2 GT h−1 k δG s hn kL2 (∂Ω e )

!1 2

e

!1

2

X

≤ C k Cy k G e kHs (Ωh ) hµ−2−ε σ s

Ghn K k2L2 (∂Ω e ) h−1 k JδG s

e

. (212)

We now substitute Eqs. (211, 212) in Eq. (210) to obtain the final bound of the second term Ge , y ; δG Gh ) as of N (G !1 2 X 2 −1 G K k k JδG σ h | I2 |≤ C k Cy k G e kHs (Ωh ) hµ−2−ε . (213) hn s s L2 (∂Ω e ) e

Ge , y ; δG Gh ) as decomposed in Eq. (74), using Taylor’s Furthermore, for the third term of N (G expansion as in Eqs. (60-61), the generalized H¨ older’s inequality, the generalized CauchySchwartz’ inequality, the definition of Cy in Eq. (67), and the bounds (202, 203), leads to r z  X Z ¯ ∇G y)∇δG Gh dS | | I3 | ≤ | ζ nT ζ Tw GG (y ∂I Ω e

e

!1

!1

4

−1 2

X

≤ Cy hs

e

!1

4

X

k Jζζ K k4L4 (∂Ω e )

k ζ k4L4 (∂Ω e )

2

X

e

Gh k2L2 (∂Ω e ) hs k ∇δG

e

!1 2

≤ Cy C k k G e kHs (Ωh ) hµ−2−ε σ s

X

Gh |2H1 (∂Ω e ) hs | δG

.

e

(214) Ge

Gh ) defined in Eq. (74) is bounded using a Taylor’s Likewise, the fourth term of N (G G , y ; δG expansion as in Eqs. (60-61), the generalized H¨ older’s inequality, the generalized CauchySchwartz’ inequality, the definition of Cy in Eq. (67), and the bounds (202, 203) leading to  r z B X Z ¯ ∇G y) JδG Ghn K dS | | I4 | ≤ | ζ nT ζ Tw GG (y hs ∂I Ω e e !1 !1 !1 4 4 2 X X X −1 4 4 −1 2 2 ζ ζ G ≤ Cy hs k Jζ K kL4 (∂Ω e ) k kL4 (∂Ω e ) hs k JδG hn K kL2 (∂Ω e ) e

e

e

!1

2

≤ C k Cy k G e kHs (Ωh ) hµ−2−ε σ s

X e

Ghn K k2L2 (∂Ω e ) h−1 k JδG s

. (215)

Ge , y; δG Gh ) defined in Eq. (74) is derived using Eq. Then the bound of the fifth term of N (G (231) developed in C.4.3, following r z X Z Ge − G h )δG Ghn dS | Ge − G h )T I (G | I5 | ≤ 2Cy | (G e

X 1 | + Cy 8 e

∂I Ω e ∪∂D Ω e

Z ∂I Ω e

r z Ge − G h )T I JG Ge − G h K JδG Ghn K dS |=| I51 | + | I52 | , (G

(216)

Title Suppressed Due to Excessive Length

57

T . Using the generalized H¨ where I is a matrix of unit norm and has the same size of o¯G older’s inequality, the generalized Cauchy-Schwartz’ inequality, and the bounds (202, 203) one has r z X Z Ghn ) dS | | I51 | ≤ 2Cy | ζ T I (ζζ δG e

∂Ω e

−1 2

X

!1

!1 4

≤ 2Cy hs

e

!1

4

X

k Jζζ K k4L4 (∂Ω e )

2

X

k ζ k4L4 (∂Ω e )

e

Gh k2L2 (∂Ω e ) hs k δG

e

!1

2

X

≤ 2Cy C k k G e kHs (Ωh ) hµ−2−ε σ s

Gh k2L2 (∂Ω e ) hs k δG

,

e

(217) | I52

r z X Z 1 Ghn K dS | | ζ T I Jζζ K JδG | ≤ Cy e 8 ∂Ω e !1 !1 4 4 1 X X 1 4 4 2 ≤ Cy hs k Jζζ K kL4 (∂Ω e ) k Jζζ K kL4 (∂Ω e ) 8 e e ≤

1 σ Cy C k k G e kHs (Ωh ) hµ−ε s 8

!1 2

X

h−1 s

e

Gh K k JδG

k2L2 (∂Ω e )

!1

2

X e

Gh K k2L2 (∂Ω e ) k JδG h−1 s

. (218)

Combining Eqs. (217 and 218) leads to the final bound !1

2

k

| I5 | ≤ 2Cy C k G

e

X

kHs (Ωh ) hµ−2−ε σ s

hs k

Gh k2L2 (∂Ω e ) δG

e

(219)

!1

2

1 + Cy C k k G e kHs (Ωh ) hµ−ε σ s 8

X e

Gh K k2L2 (∂Ω e ) h−1 k JδG s

.

Ge , y ; δG Gh ) defined in Eq. (74), we rewrite it using Finally to bound the last term of N (G Eq. (63) as Z XZ ¯GG (y ¯ ζ , ∇ζζ )dΩ = GT GT ζ Td y, ∇y y)ζζ )dΩ δG δG I6 = h Rd (ζ h (ζ Ωh

+2

e

Ωe

e

XZ Ωe

(220)

¯G ∇G GT ζ Td y)∇ζζ )dΩ = I61 + 2I62 . δG G (y h (ζ

The two contributions are bounded using the generalized H¨ older’s inequality, the generalized Cauchy-Schwartz’ inequality, the definition of Cy in Eq. (67), the bounds (200, 201), and the inverse inequality of Lemma 8, following !1

!1

4

| I61 | ≤ Cy

X

!1

4

X

k ζ k4L4 (Ω e )

e

2

X

k ζ k4L4 (Ω e )

e

Gh k2L2 (Ω e ) k δG

(221)

e

Gh kL2 (Ωh ) k G e kHs (Ωh ) . ≤ C k Cy hµ−2−ε σ k δG s !1

!1

4

| I62 | ≤ Cy

X e



k4L4 (Ω e )

!1

2

X

k ∇ζζ

k2L2 (Ω e )

e

4

X

k

Gh k4L4 (Ω e ) δG

e

(222)

Gh kL2 (Ωh ) k G e kHs (Ωh ) . ≤ C k Cy hµ−2−ε σ k δG s Substituting Eqs. (221, 222) in Eq. (220) leads to Gh kL2 (Ωh ) k G e kHs (Ωh ) . | I6 | ≤ C k Cy hµ−2−ε σ k δG s

(223)

58

Lina Homsi, Ludovic Noels

Combining Eqs. (209, 213, 214, 215, 216, and 223), yields h Gh kH1 (Ωh ) + Ge , y ; δG Gh ) |≤ Ck Cy k G e kHs (Ωh ) hsµ−2−ε σ k δG | N (G !1  !1 2 2 X X 2 2 −1 . Ghn K k 2 Gh k 1 hs k δG + hs k JδG e e H (∂Ω )

e

(224)

L (∂Ω )

n

e

Finally, using the definition of the energy norm (77), this results yields the bound (95) of Ge , y ; δG Gh ). N (G

Ge , y ; δG Gh ) C.4.3 Declaration related to the fifth term of N (G Using the identities JabK = JaK hbi + hai JbK and hai hbi = habi − 41 JaK JbK on ∂I Ωh , the term q e y T (G G − G h )To¯G Gh )(G Ge − G h ) hδG Ghn i can be rewritten with an abuse of notations on the (G product operator as E D z T T Gh )δG Ghn JG Ge − G h K Ge − G h )To¯G Gh )(G Ge − G h ) hδG Ghn i = (G Ge − G h )To¯G (G (G (G z r z 1r e T T Gh ) hG Ge − G h i hδG Ghn i . Ge − G h )To¯G Gh ) JG Ge − G h K JδG Ghn K + (G G − G h )To¯G (G − (G (G 4 (225)

r

q e y T (G G − G h )To¯G Gh ) , where o¯G (G Gh ) corresponds Now, we need to solve explicitly the term (G Gh ) defined in Eq. (181) with pGG (V Vt ) = −o oG (V Vt ), yielding to −¯ pGG (G 1

Z Gh ) = o¯G (G

oG (V Vt )dt, (1 − t)o

(226)

0

Gh − G e ). As oG only involves terms in with V t = G e + t(G component. 1

Z

we compute α ¯ the nonzero

2 )dt. [feT + t(fT − feT )]3

(227)

1 2 dt = 2 . [feT + t(fT − feT )]3 fT feT

(228)

(1 − t)(

α ¯ = 3Kαth

2 , f3 T

0

For simplicity, let us define λ as Z

1

(1 − t)

λ= 0

q e y T (G G − G h )To¯G Gh ) , we need λ(feT − fT ) which reads It can be noticed that to evaluate (G λ(feT − fT ) =

1 2 fT feT

(feT − fT ) =

1 1 1 ( − e ), feT fT fT

(229)

and the jump of the last result is

Jλ(feT

− fT )K =

  

1 ( 1+ − f1e − fe T fT T 1 e 2 (fT − fT ) e fT fT

1 f− T

=

1 ) fe T − 1e2 fT fT

+

q y 1 =− e + fT − feT on ∂I Ωh fT fT f− q yT fT − feT on ∂D Ωh .

(230)

Title Suppressed Due to Excessive Length

59

Hence considering this equation in the matrix form, then substituting it in Eq. (225), and using the definition of Cy in Eq. (67), lead to z r XZ T Gh )(G Ge − G h ) hδG Ghn i dS | Ge − G h )To¯G (G | (G e

∂I Ω e ∪∂D Ω e

X Z ≤ Cy | ∂I Ω e

e

X 1 Cy | 8 e Z X + Cy |

Ge − Gh )T I JG Ge − Gh K δG Ghn dS | (G

Z

r

+

e

∂I

Ωe

r

∂I

Ωe

X Z + Cy | e

z Ge − G h )T I (G Ge − G h )δG Ghn dS | (G

r ∂D Ω e

z Ge − G h )T I JG Ge − G h K JδG Ghn K dS | (G

z Ge − G h )T I (G Ge − G h )δG Ghn dS | , (G

T. where I is a matrix of unit norm and has the same size of o¯G

(231)