A coupled electro-thermal Discontinuous Galerkin

Journal of Computational Physics 348 (2017) 231–258

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A coupled electro-thermal Discontinuous Galerkin method L. Homsi a , C. Geuzaine b , L. Noels a,∗ a

Computational & Multiscale Mechanics of Materials (CM3), Department of Aerospace and Mechanical Engineering, University of Liège, Quartier Polytech 1, Allée de la Découverte 9, B-4000 Liège, Belgium b Department of Electrical Engineering and Computer Science, University of Liège, Quartier Polytech 1, Allée de la Découverte 10, B-4000 Liège, Belgium

a r t i c l e

i n f o

Article history: Received 27 August 2016 Received in revised form 10 April 2017 Accepted 12 July 2017 Available online 18 July 2017 Keywords: Discontinuous Galerkin method Electro-thermal coupling Nonlinear elliptic problem Error estimates

a b s t r a c t This paper presents a Discontinuous Galerkin scheme in order to solve the nonlinear elliptic partial differential equations of coupled electro-thermal problems. In this paper we discuss the fundamental equations for the transport of electricity and heat, in terms of macroscopic variables such as temperature and electric potential. A fully coupled nonlinear weak formulation for electro-thermal problems is developed based on continuum mechanics equations expressed in terms of energetically conjugated pair of fluxes and fields gradients. The weak form can thus be formulated as a Discontinuous Galerkin method. The existence and uniqueness of the weak form solution are proved. The numerical properties of the nonlinear elliptic problems i.e., consistency and stability, are demonstrated under specific conditions, i.e. use of high enough stabilization parameter and at least quadratic polynomial approximations. Moreover the prior error estimates in the H1 -norm and in the L2 -norm are shown to be optimal in the mesh size with the polynomial approximation degree. © 2017 Elsevier Inc. All rights reserved.

1. Introduction Electro-thermal materials received a significant interest in recent years due to their capability to convert electricity directly into heat and vice versa, which promises a wide range of applications in energy and environment fields. The main interest of this work is to derive a consistent and stable Discontinuous Galerkin (DG) method for two-way electro-thermal coupling analyzes considering electro-thermal effects such as Seebeck and Peltier effects, but also Joule heating. These effects describe the direct conversion of the difference in electric potential into a temperature difference within the system (Peltier effect), and vice versa (Seebeck effect). This is typical of thermo-electric cells which could work in two ways: electric generations [1,2] and heat pumps which operate in cool or heat modes [3]. Electro-thermal continuum has extensively been developed in the literature [3–6]. For example, as a non-exhaustive list, Ebling et al. [1] have implemented thermo-electric elements into the finite element method and have validated it by analytical and experimental results for the figure of merit values. Liu [6] has developed a continuum theory of thermo-electric bodies. He has applied it to predict the effective properties of thermo-electric composites. However he has considered that the temperature and voltage variations are small, which leads to a linear system of partial differential equations. Pérez-

*

Corresponding author. E-mail addresses: [email protected] (L. Homsi), [email protected] (C. Geuzaine), [email protected] (L. Noels).

http://dx.doi.org/10.1016/j.jcp.2017.07.028 0021-9991/© 2017 Elsevier Inc. All rights reserved.

232

L. Homsi et al. / Journal of Computational Physics 348 (2017) 231–258

Aparicio et al. [5] have proposed an electro-thermal formulation for simple configurations and have provided a comparison between analytical and numerical results. The constitutive equations that govern electro-thermal coupling can be formulated in term of ( −TV , 1T ) instead of (V, T), where T is the temperature and V is the electric potential. Such a formulation has been considered in the literature, e.g. by Mahan [4], Yang et al. [7], Liu [6], in order to obtain a conjugated pair of fluxes and fields gradients. Mahan [4] has provided a comparison between the different energy fluxes that have been developed and used by different researchers and concluded that all these different treatments result in the same equation. In this paper a discontinuous Galerkin weak formulation in terms of energetically conjugated fields gradients and fluxes is proposed. The main strengths of the DG method are the ease in handling elements of different types, the high order accuracy reached for higher polynomial orders and the high scalability in parallel simulations. Indeed the possibility of using irregular and non-conforming meshes in algorithm makes it suitable for time dependent transient problems as it allows for easy mesh modification dynamically with time. Above all, 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 electro-thermal coupled problems. However, if not correctly formulated, discontinuous methods can exhibit instabilities, and the numerical results fail to approximate the exact solution. For practitioners, it is important to have methods available which yield reliable results for a wide 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, e.g. [8–10]. Since the seminal work of Reed et al. [11], DG methods have been developed to solve hyperbolic, parabolic, and elliptic problems. The state of the art of DG methods and their developments can be found in [12]. Most of DG methods for elliptic and parabolic problems rely on the Interior Penalty (IP) method. The main principle of IP, as introduced by [13,14], is to constrain weakly the compatibility instead of building it into the finite element which makes it easier to use discontinuous polynomial spaces of high degree. The interest in the symmetric interior penalty (SIPG) methods, which will be considered in this paper, has been renewed by Wheeler [14] due to demands for optimality of 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 other possible choices of traces and numerical fluxes as discussed by Arnold et al. [15], who have provided an analysis of a large class of discontinuous methods for second order elliptic problems with different numerical fluxes, and declared that IP, NIPG (Non-Symmetric Interior Penalty), LDG (Local discontinuous Galerkin) and other DG methods are consistent and stable methods. In particular Arnold et al. [15] have proposed a framework for dealing with linear elliptic problems by means of DG methods and demonstrated that DG methods which are consistent, adjoint 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 [16] has derived anisotropic hp version 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. [17] have extended the DG methods from a linear self-adjoint elliptic problem to a second order nonlinear elliptic problem. The nonlinear system resulting from DG methods is then analyzed based on a fixed point argument. They have 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. [18] have proposed an analysis for the most popular DG schemes such as SIPG, NIPG and LDG methods for nonlinear elliptic problems, and the error estimates have been studied for each of these methods by reformulating the problems in a fixed point form. In addition, according to Gudi et al. [18], 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 L2 -norm is observed for NIPG and LDG. In that work a deterioration in the order of convergence in hs is noted when linear polynomials are used. Recently, DG has been used to solve coupled problems. For instance Wheeler and Sun [19] have proposed a primal DG method with interior penalty (IP) terms to solve coupled reactive transport in porous media. In that work, a cut off operator is used in the DG scheme to treat the coupling and achieve convergence. They have declared that optimal convergence rates for both flow and transport terms can be achieved if the same polynomial degree of approximation is used. However if they are different, the behavior for the coupled system is controlled by the part with the lowest degree of approximation, and the error estimate in L2 (H1 )-norm is nearly optimal in k with a loss of 12 when polynomials with different degrees are used. Furthermore, Zheng et al. [20] have proposed a DG method to solve the thermo-elastic coupling problems 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. To the authors knowledge there is no development related to DG methods for electro-thermal coupling, which is the aim of this paper. The main advantage of this work is the aptitude to deal with complex geometry and the capability of the formulation to capture the electro-thermal behavior for composite materials with high contrast: one phase has a high electric conductivity (e.g., carbon fiber) and other is a resistive material (e.g., polymers). The key point in being able to develop a stable DG method for electro-thermal coupling is to formulate the equations in terms of energetically conjugated


233

pairs of fluxes and fields gradients. Indeed, the use of energetically consistent pairs allows writing the strong form in a matrix form suitable to the derivation of a SIPG weak form as it will be demonstrated in this paper. This paper is organized as follows. Section 2 describes the governing equations of electro-thermal materials. In order to develop the DG formulation, the weak form is formulated in terms of a conjugated pair of fluxes and fields gradients, resulting in a particular choice of the test functions (δ( 1T ), δ( −TV )) and of the trial functions ( 1T , −TV ). A complete nonlinear coupled finite element algorithm for electro-thermal materials is then developed in Section 3 using the DG method to derive the weak form. This results into a set of non-linear equations which is implemented within a three-dimensional finite element code. Section 4 focuses on the demonstration of the numerical properties of the DG method under the assumption of a bi-dimensional problem, based on rewriting the nonlinear formulation in a fixed point form following closely the approach developed in [21,18]. The numerical properties of the nonlinear elliptic problem, i.e. consistency and the uniqueness of the solution, can then be demonstrated, and the prior error estimate is shown to be optimal in the mesh size for polynomial approximation degrees k > 1. In Section 5, several examples of applications in one, two, and three dimensions are provided for single and composite materials, in order to verify the accuracy and effectiveness of the electro-thermal DG formulation and to illustrate the algorithmic properties. We end by some conclusions, remarks, and perspectives in Section 6. 2. Governing equations In this section an overview of the basic equations that govern electro-thermal phenomena is presented for a structure characterized by a domain ⊂ Rd , with d = 2, or 3 the space dimension, whose external boundary is ∂. In particular we discuss the choice of the conjugated pair of fluxes and fields gradients that will be used to formulate the strong form in a matrix form. 2.1. Strong form The first balance equation is the electrical charge conservation equation. When assuming a steady state, the solution of the electrical problem consists in solving the following Poisson type equation for the electrical potential

∇ · je = 0

∀x ∈ ,

(1)

2

where j e [A/m ] denotes the flow of electrical current density vector, which is defined as the rate of charge carriers per unit area or the current per unit area. At zero temperature gradient, the current density j e is described by Ohm’s law which is the relationship between the electric potential V [V] gradient and the electric current flux per unit area through the electric conductivity l [S/m], with

j e = l · (−∇ V).

(2)

However when T [K] varies inside the body, an electromotive force (∇ V)s per unit length appears, and reads

(∇ V)s = −α ∇ T,

(3)

where α [V/K] is the Seebeck coefficient which is in general temperature dependent and defined as the derivative of the electric potential with respect to the temperature. By taking in consideration the Seebeck effect, Eq. (3), and adding it to Ohm’s Law, Eq. (2), for systems in which the particle density is homogeneous [4], the current density is rewritten as

j e = l · (−∇ V) + αl · (−∇ T).

(4)

The second balance equation is the conservation of the energy flux, which is a combination of the inter exchanges between the thermal and electric energies:

∇ · j y = −∂t y

∀ x ∈ .

(5)

The right hand side of this equilibrium equation is the time derivative of the internal energy density y [J/m3 ]

y = y0 + cv T,

(6)

which consists of the constant y0 independent of the temperature and of the electric potential, and of the volumetric heat capacity cv [J/(K · m3 )] multiplied by the absolute temperature T. Moreover the energy flux j y is defined as

jy = q + V je,

(7)

where q [W/m2 ] is the heat flux. On the one hand, at zero electric current density, j e = 0 (open circuit), the heat flux is given by the Fourier’s Law

q = k · (−∇ T),

(8)

234


where k [W/(K · m)] denotes the symmetric matrix of thermal conductivity coefficients, which may depend on the temperature. On the other hand, at zero temperature gradient, the heat flux is given by

q = βα j e = α T j e ,

(9)

where the coupling between the heat flux q and the electric current density j e is governed by the Peltier coefficient βα = α T. By superimposing the previous terms to the Fourier’s Law, Eq. (8), the thermal flux can be rewritten as:

q = k · (−∇ T) + α T j e = (k + α 2 Tll) · (−∇ T) + α Tll · (−∇ V).

(10)

The first term is due to the conduction and the second term corresponds to the Joule heating effect. Let us now define the Sobolev space Wsr (), with s a non-negative integer and r ∈ [1, ∞[, the subspace of all functions from the norm Lr () whose generalized derivatives up to order s exist and belong to Lr (), which is defined as

Wsr () = f ∈ Lr (), ∂ α f ∈ Lr () ∀ | α |≤ s, s ≥ 1 .

(11)

When r = 2, the spaces are Hilbert spaces equipped with the scalar product: Ws2 () = Hs (). For s = 0, the norm is the L2 norm. + Therefore the conservation laws are written as finding V, T ∈ H2 () × H2 () such that

∇ · je = 0

∀ x ∈ ,

(12)

∇ · j y = ∇ · q + j e · ∇ V = −∂t y

∀ x ∈ ,

(13)

2+

where T belongs to the manifold H , such that T is always strictly positive. These equations are completed by suitable boundary conditions, where the boundary ∂ is decomposed into a Dirichlet boundary ∂D and a Neumann boundary ∂N (i.e., ∂D ∪ ∂N = ∂, and ∂D ∩ ∂N = 0). On the Dirichlet BC, one has

T = T¯ > 0 ,

¯ V=V

∀ x ∈ ∂D ,

(14)

¯ are the prescribed temperature and electric potential respectively. The natural Neumann boundary conditions where T¯ and V are constraints on the secondary variables: the electric current for the electric charge equation and the energy flux for the energy equation, i.e. j e · n = ¯je ,

j y · n = ¯jy ∀ x ∈ ∂N ,

(15)

where n is the outward unit normal to the boundary ∂. For simplicity we consider the same boundary division into Neumann and Dirichlet parts for the both fields T and V. However in the general case this could be different. The set of Eqs. (12), (13) can be rewritten under a matrix form. First we rewrite Eqs. (4), (7), (10) under the form

j=

je jy

=

αl k + α Vll + α 2 Tll

l Vll + α Tll

−∇ V −∇ T

(16)

. +

The set of governing Eqs. (12), (13) thus becomes finding V, T ∈ H2 () × H2 () such that

div j =

0

−∂t y

= i,

(17)

where we have introduced i =

0

−∂t y

for a future use.

2.2. The conjugated driving forces First the weak form of the conservation of electric charge carriers, Eq. (1), is obtained by taking the inner product of this equation with a suitable scalar test function δ fV ∈ H1 ( ) over a sub-domain ⊂ , yielding

∇ · j e δ fV d = 0 ∀δ fV ∈ H1 ( ).

(18)

After a simple formal integration by parts and using the divergence theorem, we obtain

−

j e · ∇δ fV d +

j e · n δ fV dS = 0 ∀δ fV ∈ H1 ( ).

(19)

∂

Secondly, taking the inner product of the second balance equation, Eq. (13), with the test function δ fT ∈ H1 ( ), over the sub-domain ⊂ leads to


∇ · q δ fT d +

j e · ∇ Vδ fT d = −

∂t yδ fT d ∀δ fT ∈ H1 ( ).

q · ∇δ fT d =

(20)

Moreover by applying the divergence theorem, one obtains

235

∇ V · j e δ fT d +

q · n δ fT dS +

∂

∂t yδ fT d ∀δ fT ∈ H1 ( ).

(21)

By substituting the internal energy, Eq. (6), and the thermal flux, Eq. (10), this last equation reads

k · (−∇ T) + α Tjje · ∇δ fT d =

cV ∂t Tδ fT d +

∇ V · j e δ fT d +

k · (−∇ T) + α Tjje · n δ fT dS.

(22)

∂

In order to define the conjugated forces, let us substitute δ fV by − VT in Eq. (19). This results into

je · n − ∂

Substituting δ fT by

V

dS =

T

je · −

∇V

T

+

∇ T d . 2

V T

(23)

+

∈ H1 ( ) in Eq. (22) leads to:

1 T

(−∇ T)

(−∇ T) · k ·

T2

−α

je T

−∇ T cV je · ∇ T d = ( ∂t T)d + ∇ V · d + k·( ) + αj e · n dS. T

T

∂

T

(24)

By subtracting Eq. (23) from Eq. (24), one gets

=

cV T

∂t Td +

k·(

−∇ T

∂

T

) + αj e + j e ( ) · n dS V T

∇V ∇V (−∇ T) V je −j e · (−∇ T) · k · + je · − j e 2 · ∇ T d + − α · ∇ T d , 2 T

or

1 T

T

(cV ∂t T)d +

∂

T

1 q + j e V · n dS = T

T

T

(25)

−∇ T · j e V − k · ∇ T + αj e T d . 2 T

(26)

Henceforth, as j y = q + j e V, this last result is rewritten as

∂t yδ fT d +

j y · n δ fT dS =

∂

j y · ∇δ fT d .

(27)

By this way we recover the conservation equation of the energy flux, Eq. (5), which shows that j e , j y and ∇(− VT ), ∇( 1T ) is a conjugated pair of fluxes and fields gradients as derived in [6]. 2.3. Strong form in terms of the conjugated pair of fluxes and fields gradients

fV , with fV = − VT and fT = 1T , then the gradients of the fields fT vector ∇M , a 2d × 1 vector in terms of (∇ fV , ∇ fT ), is defined by Let us define a 2 × 1 vector of the unknown fields M =

∇M =

∇ fV ∇ fT

=

∇(− VT ) ∇( 1T )

=

− 1T I 0

V I T2 − 12 I T

∇V ∇T

(28)

.

Hence, the fluxes defined by Eq. (16) can be expressed in terms of fV , fT , yielding

j=

je jy

=

lT VTll + α T2l

VTll + α T2 l T2k + 2α T2 Vll + α 2 T3l + TV2l

∇ fV ∇ fT

.

(29)

The 2d × 1 fluxes vector j is the product of the fields gradients vector ∇ M , which derived from the state variables (fV , fT ), by a coefficients matrix Z (V, T) of size 2d × 2d, which is temperature and electric potential dependent. This formulation of the conjugated forces leads to a symmetric coefficients matrix Z (V, T) such that

236


Fig. 1. Interface between two elements (e+ ) and (e− ).

j = Z ∇ M.

(30)

−1 From Eq. (29), the symmetric coefficients matrix Z (V, T) is positive definite if Z 00 and Z 11 − Z T10Z 00 Z 01 are positive definite. −1 Z01 = kT2 is also positive definite, then Z(V, T) is a positive definite Since Z00 = lT is positive definite, and Z11 − ZT10 Z00 matrix. The coefficient matrix Z (V, T) in Eq. (29) could also be rewritten in terms of (fV , fT ) = (− VT , 1T ), since T = f1 , V = − ffV , as T

⎛ Z(fV , fT ) = ⎝

1 l fT

−

− fV2 l + α 12 l fT

k f2T

fT

fV l+ f2T

α f12 l T

− 2α fV3 l + α 2 13 l + fT

fT

T

⎞ f2V f3T

l

⎠.

(31)

Since the coefficients matrix is positive definite, the energy can be defined by

⎛

∇M Tj = ∇M TZ (fV , fT )∇M = ∇ fV ∇ fT ⎝

1 l fT

−

− fV2 l + α 12 l fT

fT

k f2T

fV l+ f2T

T

− 2α fV3 l + α 2 13 l + fT

⎞

α f12 l fT

f2V f3T

l

⎠

∇ fV ∇ fT

≥ 0.

(32)

Finally, the strong form (16), (17) can be expressed as

⎧ ⎪ ⎨div(j ) = i ¯ M =M ⎪ ⎩n¯ T j ¯ j = where n¯ =

n 0 0 n

∀ x ∈ , ∀ x ∈ ∂D , ∀ x ∈ ∂N ,

¯ ∈ L2 (∂D ) × L2+ (∂D ), and j¯ = ,M

(33)

¯je ¯jy .

3. Thermo-electrical analysis with the Discontinuous Galerkin (DG) finite element method Let the domain ⊂ Rd be approximated by a discretized domain h ⊂ Rd such that ≈ h = ∪e e , where a finite element in h is denoted by e . The boundary ∂h is decomposed into a region of Dirichlet boundary ∂D h , and a region e of Neumann boundary ∂N h . The intersecting boundary of the finite elements is denoted by ∂I h = ∪e ∂ \ ∂h as shown in the Fig. 1, with ∂N h = ∪e ∂N e , ∂D h = ∪e ∂D e , ∂h ∪ ∂I h = ∪e ∂e , and ∂I e = ∂e ∂I h . Within this finite element discretization, an interior face (∂I )s = ∂e+ ∩ ∂e− is shared by elements e+ and e− , and n − is the unit normal vector pointing from element e− toward element e+ , see Fig. 1. Similarly, an exterior Dirichlet edge (∂D )s = ∂e ∩ ∂D h is the intersection between the boundary of the element e and the Dirichlet boundary, and s − n = n is s used to represent the outward unit normal vector. Finally (∂DI ) is a face either on ∂I h or on ∂D h , with ) = ∂ ∪ ∂ . (∂ I h D h s DI 3.1. Weak discontinuous form The weak form of the matrix form Eq. (33) will be formulated by considering a two-field coupled problem defined in terms of the energetically conjugated pair of and fields gradients as defined in Section 2.3. fluxes fV , we can define the associated norm of the standard broken Sobolev space To account for the discontinuity in M = fT s Wr (h ) of order s and exponent r with 1 ≤ r < ∞. Starting from the Sobolev space norm and semi-norm


237

⎧ 1 ⎪ α f r dx + α f r dx r , ⎨ M Ws (e ) = ∂ ∂ T V |α |≤s e |α |≤s e r 1 ⎪ r ⎩| M | s e s r s r = ∂ f dx + ∂ f dx , T V Wr ( ) e e

(34)

the norm and semi-norm of the broken Sobolev space read

⎧ ⎪ ⎨ M Ws ( r

⎪ ⎩| M |

h)

Wsr (h )

1 r M rWs (e ) , r 1 r r = . e | M |Wr (e ) =

e

(35)

s

For the case r = ∞, the norm is defined as

M Ws∞ (h ) = maxe M Ws∞ (e ) , with M Ws∞ (e ) = max|α |≤s { ∂ α fV L∞ (e ) , ∂ α fT L∞ (e ) } . We can define the broken Sobolev spaces for the case r = 2 as (+)

Xs

(+) = M ∈ L2 (h ) × L2 (h ) |M

s(+)

M |e ∈Hs (e )×H

and

(36)

(37)

,

(e ) ∀e ∈h

Y = ∇ M ∈ (L2 (h ))3 × (L2 (h ))3 |

3 3 ∇M |e ∈ H1 (e ) × H1 (e ) ∀e ∈h

(38)

,

+ where we define X+ s as the manifold such that fT > 0 while Xs is the manifold for which fT 0, with Xs ⊂ Xs . For concise-

ness, we have used an abuse of notations in Eq. (37), in which “•(+) ” holds either for “•” or for “•+ ”. Still for conciseness, (+) we define X(+) = X2 . Let us derive the weak form of the governing equation (33) for electro-thermal coupling by multiplying it by a suitable test function δM ∈ X1 , performing a volume integral, and using the divergence theorem on each element e . This leads to state the problem as finding M ∈ X+ 1 such that

−

e

∇δM Tj (M , ∇M ) d +

e

e

δM Tn¯ Tj (M , ∇M ) dS =

∂e

δM Ti d ∀δM ∈ X1 .

(39)

h

The surface integral of this last equation is rewritten as

e

δM Tn¯ Tj (M , ∇M )dS =

e

∂e

δM Tn¯ Tj (M , ∇M )dS +

e

∂N e

δM Tn¯ T j (M , ∇M )dS.

(40)

∂I e ∪∂D e

The second term of the right hand side of Eq. (40) can be rewritten using

e

δM Tn¯ Tj (M , ∇M )dS =

∂I e

e

T T T T δM − n¯ − j − (M , ∇M ) + δM + n¯ + j + (M , ∇M ) dS, n¯ + = −n¯ − ,

∂I h

δMT n¯ T j(M, ∇ M)dS = −

∂D e

(41)

−δMT n¯ − j(M, ∇ M) dS and n¯ − = n¯ . T

∂D h −

+

In these equations we use superscript “− (+)” to refer to the value of element e (e ). Moreover we define δMn = the n− 0 n− 0 − δM and Mn = − M for future use.

0 n 0 n We can introduce trace operators to manipulate the numerical flux and obtain the primal formulation. On ∂I h , the average • and the jump J•K operators are defined as • = 12 (•+ + •− ), J K = (•+ − •− ). The definition of these two trace operators can be extended on the Dirichlet boundary ∂D h as • = •, J•K = (−•). Therefore using Eqs. (41), Eq. (39) becomes

e

e

∇δM Tj (M , ∇M )d + h

δM Ti d = ∂N h

δM Tn¯ Tj (M , ∇M )dS −

r

z δMnT j (M , ∇M ) dS.

(42)

∂I h ∪∂D h

Applying the Neumann boundary conditions specified in Eq. (33) allows this last result to be rewritten as finding M ∈ X+ 1 such that

238


δM Tj¯ dS =

∇δM Tj (M , ∇M )d +

h

∂N h

r

δM Ti d +

z δMnT j (M , ∇M ) dS ∀δM ∈ X1 .

(43)

∂I h ∪∂D h

h

Applying the mathematical identity JabK = JaK b + JbK a on ∂I h , and by neglecting the second term because only consistency of the test functions needs to be enforced, then the consistent flux related to the last term of Eq. (43) reads r z

δMnT

j (M , ∇M ) .

Moreover, on the one hand, due to the discontinuous nature of the trial functions in the DG weak form, the inter-element discontinuity is allowed, so the continuity of unknown variables is enforced weakly by using symmetrization and stabilization terms at the interior elements boundary interface ∂I h . On the other hand, the Dirichlet boundary condition (33) is also enforced in a weak sense by considering the same symmetrization and stabilization terms at the Dirichlet elements boundary interface ∂D h . By using the definition of the conjugated force, Eq. (29), the virtual energy flux δj (M ) reads

δj (M ) = Z (M )∇δM .

(44)

This last result allows formulating the symmetrization and quadratic stabilization terms so the weak form Eq. (43) becomes finding M ∈ X+ 1 such that:

T¯

δM j dS −

¯ )∇δM dS + ¯ T Z (M M n

∂N h

∂D h

=

∇δM Tj (M , ∇M )d +

h

r

MnT

+ ∂I h ∪∂D h

z

hs

r

δM Ti d + ∂I h ∪∂D h

h

B

δMnT

∂D h

r

Z (M )∇δM dS +

δMnT

¯ ¯ n dS Z (M ) M

δMnT

z

j (M , ∇M ) dS

z B hs

∂I h ∪∂D h

(45)

Z (M ) JMn K dS ∀δM ∈ X1 ,

n 0 ¯ . In this equation B is the stability parameter which has to be sufficiently high to guarantee stability, M 0 n as it will be shown in Section 4, and hs is the characteristic length of the mesh, which will also be defined in Section 4. The last two terms of the left hand side of Eq. (45) make sure that the Dirichlet boundary condition (33) is weakly enforced, as it will be shown in Section 4. The last three terms of the right hand side of Eq. (45) are the interface terms, which correspond to a SIPG method:

¯n = where M

1. The first term ensures consistency despite the discontinuity of the test function δM between two elements, and involves the consistent numerical flux which is here the traditional average flux. 2. The second term achieves symmetry of the weak form and thereby also of the stiffness matrix after FE discretization. It also ensures (weakly) the continuity of solution across element boundaries and the optimal convergence rate in the L2 -norm. 3. The last term ensures stability, as it is well known that the discontinuous formulation of elliptic problems requires quadratic terms. The stabilization term depends on a stability parameter, which is independent of mesh size and material properties, as it will be shown in Section 4. Eventually the weak form (45) is thus summarized as finding M ∈ X+ 1 such that:

¯ , δM ) − a(M , δM ) = b(M

δM Ti d ∀δM ∈ X1 ,

(46)

h

with

r

∇δM Tj (M , ∇M )d +

a(M , δM ) = h

r

∂I h ∪∂D h

z

δMnT

z

r

MnT Z (M )∇δM dS +

+ ∂I h ∪∂D h

and

¯ , δM ) = b(M

∂N h

j (M , ∇M ) dS

δMnT

z B hs

∂I h ∪∂D h

δM Tj¯dS −

∂D h

¯ )∇δM dS + ¯ Z (M M n T

∂D h

δMnT

B hs

(47)

Z (M ) JMn K dS,

¯) M ¯ n dS. Z (M

(48)


239

It should be noted that the test functions in the previous equations of the weak formulation belong to X 1 , however for the numerical analysis, we will need to be in X2 . 3.2. Finite element discretization In the finite element method, the trial function M is approximated by M h , which is defined over a finite element e using the interpolation concepts in terms of the standard shape function Na ∈ R at node a, see [22], yielding

M h = N a M a , ∇M h = ∇N a M a ,

(49)

where M a denotes the nodal value of M h at node a, N a =

a

N 0 0 Na

is a matrix of the shape functions and ∇N a =

∇ Na 0 is a matrix of the shape function gradients. In order to obtain a Galerkin formulation, the test functions are 0 ∇ Na approximated using the same interpolation, i.e.

δM h = N a δM a ,

∇δM h = ∇N a δM a .

(50)

In this work, we assume a constant mesh size on the elements, but the theory can be generalized by considering bounded |e | element sizes such as in [18]. We assume the discretization is shaped with a regular mesh of size hs defined as |∂e | . We also assume shape regularity of h so that there exist constants c1 , c2 , c3 , and c4 , independent of hs , such that

c1 diam ((∂DI )s ) ≤ hs ≤ c2 diam ((∂DI )s ), and c3 diam (e ) ≤ hs ≤ c4 diam (e ). + fVh ∈ Xk of the solution is thus defined in the space The finite discontinuous polynomial approximation M h =

(51)

fTh

(+)

Xk

(+) = M h ∈ L2 (h ) × L2 (h ) |M

(+)

k e k h |e ∈P ( )×P

(e ) ∀e ∈h

(52)

, +

where Pk (e ) is the space of polynomial functions of order up to k and Pk means that the polynomial approximation remains positive. In the numerical framework, the positive nature of the polynomial approximation is not directly enforced, but results naturally from the resolution of the well-posed finite element problem. + Using these definitions, the problem becomes finding M h ∈ Xk such that

¯ , δM h ) − a(M h , δM h ) = b(M

δM Thi d ∀δM h ∈ Xk .

(53)

h

The set of Eqs. (53) can be rewritten under the form:

F aext M b = F aint M b + F aI M b ,

(54)

where M b is the vector of the unknown fields at node b. The nonlinear Eqs. (54) are solved using the Newton–Raphson scheme. To this end, the forces are written in a residual form. The predictor, iteration 0, reads M b = M b0 , the residual at iteration i reads

F aext M bi − F aint M bi − F aI M bi = R ,

(55)

and at iteration i, the first order Taylor development yields the system to be solved, i.e.

∂ F aI b ∂ F aext ∂ F aint 0 =R + − − M − M bi . ∂Mb ∂Mb ∂Mb

a

In this last equation

F aext =

e

F aint =

N aj¯dS −

s

∂N e

T

(∂D )s

∇N a j (M h , ∇M h )d +

e

e a±

T ¯ )M ¯ n dS + ∇N a Z (M

F aI ± = F I1 + F aI2± + F aI3± ,

(56)

s

e

N ai d,

(∂D )s

N an¯ T

¯) Z (M

B hs

¯ n dS, M

(57)

(58)

e

(59)

240


with the three contributions to the interface forces on ∂I h 1

F aI1± =

s

F aI2± = F aI3± =

T ±N a± n¯ − j (M h , ∇M h ) dS,

(∂I )s

1 2

s

T ∇N a± Z ± (M h ) JM hn K dS,

(61)

T B ±N a± n¯ − Z (M h ) JM hn K dS.

(62)

(∂I )s

s

(60)

hs

(∂I )s

In these equations the symbol “±” refers to the node a± (with “+” for node a+ and “−” for node a− ). This system is solved by means of a Newton–Raphson method with the stiffness matrix computed in Appendix A, where the iterations continue until the convergence to a specified tolerance is achieved. 4. Numerical properties for DG method In this section, the numerical properties of the weak formulation stated by Eq. (46) are studied in steady state conditions (ii = 0), and under the assumption that d = 2. It is demonstrated that the framework satisfies two fundamental properties of a numerical method: consistency and stability. Moreover we show that the method possesses the optimal convergence rate with respect to the mesh size. 4.1. Discontinuous space and finite element properties In this part, we will assume that ∂D h = ∂h . This assumption is not restrictive but simplifies the demonstrations. Let us also define the norms, which appear in the analysis of the interior penalty method, for M ∈ X1

| M |2∗ =

e

| M |2 =

e

∇M 2L2 (e ) + M 2H1 (e ) +

1 2 h− s JMn K 2

L (∂e )

e

1 2 h− s JMn K 2

L (∂e )

e

,

(63)

,

(64)

and

| M |21 =

e

M 2H1 (e ) +

e

hs M 2 1 H

(∂e )

+

1 2 h− s JMn K 2

L (∂e )

e

,

(65)

where ∂e = ∂I e ∪ ∂D e . Eqs. (63)–(65) define norms as |M |∗ = 0 only when M = cst on h and is equal to 0 on ∂D h . Lemma 4.1 (Relation between energy norms on the finite element space). From [14], for M h ∈ Xk , there exists a positive constant Ck , depending on k, such that

| M h |1 ≤ Ck | M h |,

(66)

where Ck is a positive constant, independent of the mesh size. The demonstration directly follows by bounding the extra terms

e hs

M 22

L (∂e )

and

e hs

∇M 22

L (∂e )

of the norm

defined by Eq. (65), in comparison to the norm defined by Eq. (64), using successively the trace inequality, Eq. (B.5), and the inverse inequality, Eq. (B.9), for the first term, and the trace inequality on the finite element space, Eq. (B.6), for the second term. 4.2. Consistency +

2 To prove the consistency of the method, the exact solution M e ∈ H2 () × H problem stated by Eqs. (33) is () of the e e e ¯ ¯ ) = Z(Me ) M j j M M M j considered. This implies J K = 0, = on ∂I h , and J K = − = − , = j = Z(Me )∇ Me , and Z(M) = Z(M on ∂D h . Therefore, Eq. (46) becomes:

1 The contributions on ∂D h can be directly deduced by removing the factor (1/2) accordingly to the definition of the average flux on the Dirichlet boundary.


δM Tj¯dS −

∂N h

¯ )∇δM dS + ¯ Z (M M n T

∂D h

∇δM Tj (M e , ∇M e )d +

= h

r

δMnT

hs

∂D h

z δMnT j (M e , ∇M e )dS δMnT j (M e , ∇M e )dS −

∂D h

(67)

∂D h

T

241

B ¯ ¯ n dS Z (M ) M

∂I h

Mne Z (M e )∇δM dS +

−

δMnT

B hs

Z (M e )Mne dS ∀δM ∈ X .

∂D h

Integrating the first term of the right hand side by parts leads to

e

T

e

e

∇δM j (M , ∇M )d = −

e

e

T

e

e

δM ∇j (M , ∇M )d +

e

e

δMnT j (M e , ∇M e )dS,

(68)

∂e

and Eq. (67) becomes

δMT¯jdS −

∂N h

¯ )∇δM dS + ¯ Z(M M n T

∂D h

δMnT

B hs

¯) M ¯ n dS = − Z(M

∂D h

δMnT j (M e , ∇M e )dS −

+

∂N h

T

Mne Z (M e )∇δM dS +

∂D h

e

δMnT

δMT ∇ j(Me , ∇ Me )d

e

B hs

(69)

Z (M e )Mne dS ∀δM ∈ X.

∂D h

The arbitrary nature of the test functions leads to recover the set of conservation laws and the boundary conditions stated by Eqs. (33). 4.3. Stability of the DG formulation The demonstration of the stability follows closely the approach developed by [8,14,17,18] for linear and nonlinear elliptic problems. Since the problem is herein coupled, and since the elliptic operator is particularized, we report the modified main steps of the methodology, including the Lemmata, Theorems, and demonstrations, that was initially developed in [17,18] for d = 2. The main idea to prove the solution uniqueness and to establish the prior error estimates is to linearize the problem by reformulating the nonlinear problem in a fixed point form which is the solution of the linearized problem as proposed in [17,18,23]. Starting from the definition of the matrix Z (M ), Eq. (31), which is a symmetric and positive definite matrix, as we have proved in Section 2.2, let us define the minimum and maximum eigenvalues of the matrix Z (M ) as λ(M ) and (M ), then for all ξ ∈ R2d 0

0 < λ(M )|ξ |2 ≤ ξiZ ij (M )ξj ≤ (M )|ξ |2 .

(70)

Also by assuming that M W1 ≤ α , then there is a positive constant Cα such that ∞

0 < Cα < λ(M ).

(71)

In the subsequent analysis, we use the following integral form of the Taylor’s expansions of j , defined in Eq. (30), for

(V , ∇P ) ∈ X+ × Y in terms of (M , ∇M ) ∈ X+ × Y

j (V , ∇P ) − j (M , ∇M ) = −jM (M , ∇M )(M − V ) − j ∇M (M )(∇M − ∇P ) + R¯ j (M − V , ∇M − ∇P )

= −j¯M (M , ∇M )(M − V ) − j¯∇M (M )(∇M − ∇P ),

(72)

where jM is the partial derivative of j with respect to M , j ∇M is the partial derivative of j with respect to ∇M expressed in the matrix form, and where j¯M , j¯∇M , and R¯ j are the remainder terms. With V t = M + t(V − M ), ∇P t = ∇M + t(∇P − ∇M ), we have

j¯M (M , ∇M ) =

1 0

jM (V t , ∇P t )dt, j¯∇M (M , ∇M ) =

1 j ∇M (V t , ∇P t )dt, 0

R¯ j (M − V, ∇ M − ∇ P) = (M − V)T¯jMM (V, ∇ P)(M − V) + 2(M − V)T¯jM ∇M (V, ∇ P)(∇ M − ∇ P), and

(73) (74)

242


j¯MM (V , ∇P ) =

1

(1 − t)jMM (V t , ∇P t )dt, j¯M ∇M (V , ∇P ) =

0

1 (1 − t)jM ∇M (V t , ∇P t )dt.

(75)

0

2 Z Z Using the definition Eq. (30) of j , we have jM = ∂∂M ∇M , j ∇M = Z , jMM = ∂∂MZ2 ∇M , jM ∇M = j ∇MM = ∂∂M , j ∇M ∇M = 0. If fT ≥ + + ∞ d d ∞ ¯ ¯ ¯ ¯ ¯ fT0 > 0, then jM , jMM ∈ L ( × R × R0 × R × R ) and j ∇M , jM ∇M , j ∇MM ∈ L ( × R × R0 ). Since j is a twice continuously differential function with all the derivatives through the second order, which can be shown to be locally bounded in a ball around M ∈ R × R+ 0 following the argumentation of [17,18] for d = 2, we denote by Cy

Cy = max j W2 (×R×R+ ×Rd ×Rd ) , j¯M , j¯MM L∞ (×R×R+ ×Rd ×Rd ) , j¯∇M , j¯M ∇M , j¯∇MM L∞ (×R×R+ ) . ∞ 0 0 0

(76)

We can now study the weak form defined by Eq. (46) under the assumptions i = 0 and j¯ independent of M . The problem thus reads as finding M ∈ X+ such that

¯ , δM ) ∀δM ∈ X, a(M , δM ) = b(M

(77)

¯ , δM ) by Eq. (48). with a(M , δM ) defined by Eq. (47) and b(M 4.3.1. Derivation of the non-self-adjoint linear elliptic problem + Let us define M e ∈ H2 () × H2 () the solution of the strong form stated by Eq. (33). Thus as JM e K = 0 on ∂I e and as e e e ¯ JM K = −M = −M on ∂D , we have

a(M e , δM e ) =

T

∇δM e j (M e , ∇M e )d +

h

−

r

T

δMne

z

j (M e , ∇M e ) dS −

∂I h

¯ T Z (M e )∇δM e dS + M n

∂D h

δMne Tj (M e , ∇M e )dS

∂D h

δMne T

B

¯ n dS = b(M ¯ , δM e ) ∀δM e ∈ X, Z (M )M

(78)

e

hs

∂D h

as the weak form stated by Eq. (46) is consistent, see Section 4.2. Using the weak formulation (77), we state the Discontin+ uous Galerkin finite element method for the problem as finding M h ∈ Xk , such that

¯ , δM h ) ∀δ M h ∈ Xk ⊂ X. a(M h , δM h ) = b(M

(79)

Therefore, using δrM e = δM h in it from the DG " adding and subtracting succesz Eq. (78), subtracting r discretization z ! (79), then Bj T T eT e eT e M M j M M M M M sively ∂ ∪∂ )∇δ h dS and ∂ ∪∂ ( ) JδMhn K dS, yield ∇M ( n − n − hn hn h ∇M I

h

D

h

n

I

0 = a(M e , δM h ) − a(M h , δM h ) = h

D

h

n

h

s

∇δM Th j (M e , ∇M e ) − j (M h , ∇M h ) d

+ ∂I h ∪∂D h

r

δM Thn

z

r

T

z

r

T

z

r

T

z B

Mne − M Thn

+ ∂I h ∪∂D h

Mne − M Thn

− ∂I h ∪∂D h

Mne − MThn

− ∂I h ∪∂D h

+ ∂I h ∪∂D h

j (M e , ∇M e ) − j (M h , ∇M h ) dS

r

T

Mne − M Thn

j ∇M (M e )∇δM h dS

(80)

j ∇M (M e ) − j ∇M (M h ) ∇δM h dS

hs

z B hs

j∇M (Me ) − j∇M (Mh )

JδMhn K dS

j ∇ M (M e ) JδM hn K dS ∀δM h ∈ Xk . +

Using the Taylor series defined in Eq. (72) to rewrite the first two terms, the set of Eqs. (80) becomes finding M h ∈ Xk such that:

A(M e ; M e − M h , δM h ) + B(M e ; M e − M h , δM h ) = N (M e , M h ; δM h )

∀δM h ∈ Xk .

(81)


243

In this last equation, for given ψ ∈ X+ , ω ∈ X and δω ∈ X, we have defined the following forms:

r

ψ )∇ω d + ∇δω Tj ∇ψ ψ (ψ

ψ ; ω , δω ) = A(ψ h

r

+

T n

ω

∂I h ∪∂D h

z

z

r

ψ )∇δω dS + j ∇ψ ψ (ψ

ψ , ∇ψ ψ )ω d + ∇δω T jψ (ψ

ψ ) ∇ω dS j ∇ψ ψ (ψ

∂I h ∪∂D h

ψ ; ω , δω ) = B(ψ

δωnT

T n

ω

∂I h ∪∂D h

r

δωnT

z B hs

ψ ) Jδωn K dS, j ∇ψ ψ (ψ

z

(82)

ψ , ∇ψ ψ )ω dS. jψ (ψ

(83)

∂I h ∪∂D h

h

ψ ; ., .) and the form B (ψ ψ ; ., .) are bi-linear. The nonlinearity of Eq. (81) is thus gathered in the For fixed ψ , the form A(ψ term N (M e , M h ; δM h ) which reads

N (M e , M h ; δM h ) = h

∇δM Th (R¯ j (M e − M h , ∇M e − ∇M h ))d

r

+ ∂I h ∪∂D h

δM Thn

z

R¯ j (M e − M h , ∇M e − ∇M h ) dS

r

T

z

r

T

z B

Mne − M Thn

+ ∂I h ∪∂D h

Mne − M Thn

+

(84)

j ∇M (M e ) − j ∇M (M h ) ∇δM h dS

∂I h ∪∂D h

hs

j ∇M (M e ) − j ∇M (M h )

JδM hn K dS.

4.3.2. Solution uniqueness + + Let us first define η = IhM − M e ∈ X, with IhM ∈ Xk the interpolant of M e in Xk . The last relation (81) thus becomes

A(M e ; IhM − M h , δM h ) + B(M e ; IhM − M h , δM h ) = A(M e ; η , δM h ) + B(M e ; η , δM h ) + N (M e , M h ; δM h )

∀δM h ∈ Xk .

(85)

In order to prove the existence of a solution M h of the problem stated by Eq. (80), which corresponds to the DG finite + + element discretization (79), the problem is stated in the fixed point formulation and a map S h : Xk → Xk is defined as + + follows [18]: for a given y ∈ Xk , find Sh (y ) = My ∈ Xk , such that

A(M e ; IhM − My , δM h ) + B(M e ; IhM − My , δM h ) = A(M e ; η ,δδM h ) + B(M e ; η , δM h ) + N (M e , y ; δM h )

∀δM h ∈ Xk .

(86)

The existence of a fixed point of the map Sh is equivalent to the existence of a solution M h of the discrete problem (79), see [18]. For the following analysis, we denote by Ck , a positive generic constant which is independent of the mesh size, but may depend on CT , CkD , CkI , CkK , and on k which are defined in Appendix B, and can take different values at different places. To demonstrate the uniqueness, we have recourse to the following Lemmata, following [18]. Lemma 4.2 (Lower bound). For B larger than a constant, which depends on the polynomial approximation only, there exist two constants Ck1 and Ck2 , such that

A(M e ; δM h , δM h ) + B(M e ; δM h , δM h ) ≥ Ck1 | δM h |2∗ −Ck2 δM h 2L2 () ∀δM h ∈ Xk , e

e

A(M ; δM h , δM h ) + B(M ; δM h , δM h ) ≥

Ck1

| δM h |

2

−Ck2

δM h 2L2 ()

k

∀δM h ∈ X .

(87) (88)

Proceeding by using the bounds (71) and (76), the Cauchy–Schwartz’ inequality, the trace inequality on the finite element ξ space (B.6), the trace inequality, Eq. (B.4), and inverse inequality, Eq. (B.9), the ξ -inequality –ξ > 0 : |ab| ≤ 4 a2 + ξ1 b2 , as in Wheeler et al. [14] and Prudhomme et al. [8] analyzes, yields to prove this Lemma 4.2. The two positive constants Ck1 , Ck2 are independent of the mesh size, but do depend on k and B . The proof is detailled in [24]. In particular, for Ck1 to be positive, the following constraint on the stabilization parameter should be satisfied: B >

C2y

C2α

Ck . Therefore for the method to be stable,

the stabilization parameter should be large enough depending on the polynomial approximation.

244


Lemma 4.3 (Upper bound). There exist C > 0 and Ck > 0 such that

| A(M e ; u , δM ) + B(M e ; u , δM ) | ≤ C | u |1 | δM |1 ∀u , δM ∈ X, e

e

(89)

k

k

(90)

| A(M ; u h , δM h ) + B(M ; u h , δM h ) | ≤ C | u h | | δM h | ∀u h , δM h ∈ X .

(91)

| A(M ; u , δM h ) + B(M ; u , δM h ) | ≤ C | u |1 | δM h | ∀u ∈ X , δM h ∈ X , e

e

k

k

Applying the Hölder’s inequality, and the bound (76) on each term of A(M e ; u , δM ) + B (M e ; u , δM ) and then applying the Cauchy–Schwartz’ inequality, lead to relation (89). Therefore relations (90) and (91) are easily deduced from the relation between energy norms on the finite element space, Eq. (66). The proof is detailled in [24]. Lemma 4.4 (Energy bound). Let M e ∈ Xs , s ≥ 2, and let IhM ∈ Xk , s ≥ 2, be its interpolant. Therefore, there is a positive constant Ck > 0 independent of hs , such that

| M e − IhM |1 ≤ Ck hμ−1 M e Hs (h ) ,

(92)

with μ = min {s, k + 1}. The proof follows from applying the interpolation inequalities, Eq. (B.1) and Eq. (B.3), on the mesh dependent norm (65). Lemma 4.5 (Auxiliary problem). We consider the following auxiliary problem, with φ ∈ L2 () × L 2 ():

ψ + jM (M e , ∇M e )ψ ψ = φ on , −∇ T j ∇M (M e )∇ψ

(93)

ψ = 0 on ∂.

Assuming regular ellipticity of the operator, there is a unique solution ψ ∈ H2 () × H2 () to the problem stated by Eq. (93), which satisfies the elliptic property

ψ H2 (h ) ≤ C φ L2 (h ) .

(94)

Moreover, for a given ξ ∈ L2 (h ) × L2 (h ) there exists a unique φ h ∈ Xk such that

φ h ) + B(M e ; δM h ,φ φ h) = A(M e ; δM h ,φ

e

ξ T δM h d ∀δM h ∈ Xk ,

(95)

e

and there is a constant Ck such that:

| φ h |≤ Ck ξ L2 (h ) .

(96)

The proof of the first statement is given in [25], by combining [25, Theorem 8.3] to [25, Lemma 9.17]. The proof of the second statement follows the methodology described by Gudi et al. [18]: The use of Lemma 4.2 and Eq. (95) with δM h = φ h allows bounding | φ h | in terms of ξ L2 (h ) and φ h L2 (h ) ; φ h L2 (h ) is then estimated by considering φ = φ h ∈ Xk in Eq. (93), multiplying the result by φ h and integrating it by parts on h yielding φ h 22 = A(M e ; ψ , φ h ) + L (h )

B(M e ; ψ , φ h ). Inserting the interpolant Ihψ in these last terms, i.e. φ h 2L2 ( ) = A(M e ; ψ − Ih ψ , φ h ) + B(M e ; ψ − Ih ψ , φ h ) + h A(M e ; Ih ψ , φ h ) + B(G e ; Ihψ , φ h ), bounding the last two terms by considering δM h = Ihψ in Eq. (95), bounding the first two terms by making successive use of Lemmata 4.3 and 4.4, and using the regular ellipticity, Eq. (94), allow deriving the bound φ h L2 (h ) ≤ Ck ξ L2 (h ) , which results into the proof of (96). The existence of the solution of the discrete problem is demonstrated by proving, using these Lemmata [18], that the map Sh has a fixed point. Theorem 4.6 (Solution uniqueness to the problem stated by Eq. (86)). The solution My to the problem stated by Eq. (86) is unique for +

a given y ∈ Xk with Sh (y ) = My . Let us assume that there are two distinct solutions My 1 , My 2 , which result into

A(M e ; IhM − My 1 , δM h ) + B(M e ; IhM − My 1 , δM h ) = A(M e ; IhM − My 2 , δM h ) + B(M e ; IhM − My 2 , δM h )

∀ δM h ∈ Xk .

(97)

For fixed M e , A and B are bi-linear, therefore this last relation becomes

A(M e ; My 1 − My 2 , δM h ) + B(M e ; My 1 − My 2 , δM h ) = 0 ∀ δM h ∈ Xk .

(98)


245

Using Lemma 4.5, with ξ = δM h = My 1 − My 2 ∈ Xk results in stating that there is a unique h ∈ Xk solution of the problem Eq. (95), with

h ) + B(M e ; My 1 − My 2 , h ) = My 1 − My 2 2L2 ( ) , A(M e ; My 1 − My 2 ,

(99)

h

and that | h |≤ Ck My 1 − My 2 L2 ( ) . Choosing δM h as h in Eq. (98), we have My 1 − My 2 L2 ( ) = 0. Therefore, the h h solution Sh (y) = My is unique. +

It is now demonstrated that Sh maps a ball Oσ (IhM ) ⊂ Xk into itself and is continuous in the ball. The ball Oσ is defined with a radius σ and is centered at the interpolant IhM of M e as [18]

+

Oσ (IhM ) = y ∈ Xk such that | IhM − y |1 ≤ σ , with

σ=

| IhM − M e |1 ε

hs

1

, 0 n +

255

1 2

μ−n− 12

M − M h Hn (∂e ) ≤ CkD hs

(B.3)

M Hs (e ) ,

where μ = min {s, k + 1}. The proof of the first and third properties can be found in [26], and the proof of the second property in the particular case of d = 2 can be found in [27,28], see also the discussion by [18]. Remarks i) The approximation property in (2) is still valid for r = ∞, see [29]. ii) For M ∈ Xs , let us define the interpolant IhM ∈ Xk by IhM |e = M h (M |e ), which means IhM satisfies the interpolant inequality property provided in Lemma B.1, see [23]. Lemma B.2 (Trace inequality). For all M ∈ Hs+1 (e ) × Hs+1 (e ), there exists a positive constant CT , such that

M rWs (∂e ) ≤ CT r

1

M rWs (e ) + M rW−s1 (e ) ∇ s+1M L2 (e ) , r 2r−2 hs

(B.4)

where s = 0, 1 and r = 2, 4, or in other words

M 22

L (∂e )

≤ CT

M 4L4 (∂e ) ≤ CT

1 hs 1 hs

M 22

L (e )

+ M L2 (e ) ∇ M L2 (e ) ,

M 4L4 (e ) + M 3L6 (e ) ∇M L2 (e ) .

(B.5)

The first equation, s = 0 and r = 2, is proved in [8], and the second one, r = 4 and s = 0, is proved in [30]. Lemma B.3 (Trace inequality on the finite element space). For all M h ∈ Pk (e ) × Pk (e ) there exists a constant CkK > 0 depending on k, such that − 12

∇ lM h L2 (∂e ) ≤ CkK hs where CkK = supυ ∈PK (e )

∇ lM h L2 (e )

hs ∇M h 22

L (∂e )

∇M h 22

L (e )

| | |∂e | . e

l = 0, 1 ,

(B.6)

is a constant which depends on the degree of the polynomial approximation only, with hs =

We refer to [31] for more details. Lemma B.4 (Inverse inequality). For M h ∈ Pk (e ) × Pk (e ) and r ≥ 2, there exists CkI > 0, such that d

M h Lr (e ) ≤ CkI hsr

− d2 d −1

M h Lr (∂e ) ≤ CkI hs r

M h L2 (e ) , 1 − d− 2

M h L2 (∂e ) ,

1 ∇M h L2 (e ) ≤ CkI h− s M h L2 (e ) .

(B.7) (B.8) (B.9)

The proof of these properties can be found in [32, Theorem 3.2.6]. Note that Eqs. (B.7)–(B.8) involve the space dimension d. Appendix C. Sh maps the ball Oσ (Ih M) into itself +

Let y ∈ Oσ (IhM ) ⊂ Xk and Sh (y ) = My be a solution of the problem given by Eq. (86). Then using Lemma 4.2, Eq. (88), Lemma 4.3, Eq. (90), Lemma 4.7, Eq. (104), and the definition of the ball (100), we successively find that

256


Ck1 | IhM − My |2 −Ck2 IhM − My 22

L (h )

≤ A(M e ; IhM − My , IhM − My ) + B(M e ; IhM − My , IhM − My ) ≤ A(M e ; IhM − M e , IhM − My ) + B(M e ; IhM − M e , IhM − My ) + N (M e , y , IhM − My ) k

e

k

μ−2 −ε

e

≤ C | IhM − M |1 | IhM − My | +C Cy M Hs (h ) hs k

ε

k

μ−2 −ε

e

≤ (C hs + C Cy M Hs (h ) hs

(C.1)

σ | IhM − My |

)σ | IhM − My | .

Let us define Ck (Ck , Cy , CM ) a constant, that can depend on Ck , Cy and CM , then, as 0 < ε < rewritten for k ≥ 2:

1 , 4

the last expression can be

Ck1 | IhM − My |2 −Ck2 IhM − My 22

L (h )

≤ Ck σ hεs | IhM − My | .

(C.2)

Then, in order to estimate IhM − My L2 ( ) , we consider the auxiliary problem defined in Lemma 4.5. Choosing ξ = h δM h = IhM − My , there exists φ h such that, | φ h |≤ Ck IhM − My L2 () with

φ h ) + B(Me ; Ih M − My ,φ φ h) Ih M − My 2L2 ( ) = A(Me ; Ih M − My ,φ h

φ h ) + B(M e ; IhM − M e ,φ φ h ) + N (M e , y ;φ φ h) ≤ A(M e ; IhM − M e ,φ μ−2 −ε

≤ Ck | IhM − M e |1 | φ h | +Ck Cy M e Hs (h ) hs k

ε

k

μ−2 −ε

e

≤ (C σ hs + C Cy M Hs (h ) σ hs

σ | φ h |

(C.3)

) IhM − My L2 (h )

≤ Ck σ hεs IhM − My L2 (h ) if k ≥ 2. Substituting Eq. (C.3) in Eq. (C.2) gives

Ck1 | IhM − My |2 ≤ Ck σ hεs | IhM − My | +Ck2 IhM − My 22

L (h )

k

ε

≤ C σ hs | IhM − My |

+Ck2 Ck

(C.4)

ε

σ hs | IhM − My | if k ≥ 2.

Hence, we get

| IhM − My |≤ Ck σ hεs if k ≥ 2, and for a mesh size hs small enough and a given ball size

(C.5)

σ , IhM − My −→ 0, hence Sh maps Oσ (IhM ) to itself.

Appendix D. The continuity of the map Sh in the ball Oσ (Ih M) +

For y 1 , y 2 ∈ Oσ (IhM ) ⊂ Xk , let My 1 = Sh (y 1 ), My 2 = Sh (y 2 ) be the solutions of the linearized problem stated by Eq. (86), leading to

A(M e ; IhM − My 1 , δM h ) + B(M e ; IhM − My 1 , δM h ) = A(M e ; η ,δδM h ) + B(M e ; η , δM h ) + N (M e , y 1 ; δM h )

∀δM h ∈ Xk ,

(D.1)

and

A(M e ; IhM − My 2 , δM h ) + B(M e ; IhM − My 2 , δM h ) = A(M e ; η , δM h ) + B(M e ; η , δM h ) + N (M e , y 2 ; δM h )

∀δM h ∈ Xk , where

(D.2)

η = IhM − M e . By subtracting Eq. (D.1) from Eq. (D.2), we have A(M e ; My 2 − My 1 , δM h ) + B(M e ; My 2 − My 1 , δM h ) = N (M e , y 2 ; δM h ) − N (M e , y 1 ; δM h ).

Choosing ζ 1 = M e − y 1 ∈ X and ζ 2 = M e − y 2 ∈ X, the right hand side of Eq. (D.3) can be rewritten as follows:

(D.3)


N (M e , y 2 ; δM h ) − N (M e , y 1 ; δM h ) =

257

∇δM Th R¯ j (ζζ 2 , ∇ζζ 2 ) − R¯ j (ζζ 1 , ∇ζζ 1 ) d

h

r

+ ∂I h ∪∂D h

δM Thn

z

r

T

z

r

T

z

r

T

z B

M e − y T2n

+ ∂I h ∪∂D h

M e − y T1n

− ∂I h ∪∂D h

M e − y T2n

+ ∂I h ∪∂D h

r

T

M e − y T1n

−

R¯ j (ζζ 2 , ∇ζζ 2 ) − R¯ j (ζζ 1 , ∇ζζ 1 ) dS

j ∇M (M e ) − j ∇M (y 2 ) ∇δM h dS (D.4)

j ∇M (M e ) − j ∇M (y 1 ) ∇δM h dS

hs

z B hs

∂I h ∪∂D h

(j ∇M (M e ) − j ∇M (y 2 )) JδM hn K dS (j ∇M (M e ) − j ∇M (y 1 )) JδM hn K dS.

By applying Taylor series, Eqs. (72)–(75), to rewrite the right hand side, every term will be either in y 1 − y 2 or in ∇(y 1 − y 2 ). So proceeding as to establish Lemma 4.7 and using Lemma 4.1, Eq. (66), lead to [18] μ−2 −ε

| N (M e , y 2 ; δM h ) − N (M e , y 1 ; δM h ) | ≤ Ck Cy M e Hs (h ) hs

| y 1 − y 2 | | δM h | .

(D.5)

Choosing δM h = My 2 − My 1 , and using Eq. (88), Eq. (D.3) becomes:

Ck1 | My 2 − My 1 |2 −Ck2 My 2 − My 1 22

≤ A(M e ; My 2 − My 1 , My 2 − My 1 )

L (h )

+ B(M e ; My 2 − My 1 , My 2 − My 1 ) e

(D.6) e

≤ N (M , y 2 ; My 2 − My 1 ) − N (M , y 1 ; My 2 − My 1 ). Similarly, setting δM h = My 2 − My 1 in Eq. (D.5), Eq (D.6) becomes: μ−2 −ε

| My 2 − My 1 |2 ≤ Ck1 Cy M e Hs (h ) hs Since My 2 − My 1 22

L (h )

| y 2 − y 1 | | My 2 − My 1 | +Ck2 My 2 − My 1 2L2 ( ) . h

(D.7)

≤| My 2 − My 1 | My 2 − My 1 L2 (h ) , this last relation becomes μ−2 −ε

| My2 − My1 | ≤ Ck1 Cy M e Hs (h ) hs In order to estimate My 2 − My 1 22

L (h )

| y 2 − y 1 | +Ck2 My2 − My1 L2 (h ) .

(D.8)

, we consider ξ = My 2 − My 1 in Lemma 4.5. Therefore, there exists a unique φ h

satisfying Eq. (95) ∀δM h ∈ X . In particular for δM h = My 2 − My 1 , this implies k

φ h ) + B(M e ; My 2 − My 1 ,φ φ h) My 2 − My 1 2L2 ( ) = A(M e ; My 2 − My 1 ,φ h

φ h ) − N (M e , y 1 ;φ φ h) = N (M e , y 2 ;φ μ−2 −ε

| y 2 − y 1 || φ h |

μ−2 −ε

| y 2 − y 1 | My 2 − My 1 L2 (h ) ,

≤ Ck Cy M e Hs (h ) hs ≤ Ck Cy M e Hs (h ) hs

(D.9)

where we have used Eq (D.3), Eq. (D.5), and Eq. (96). Therefore, substituting Eq. (D.9) in Eq. (D.8) yields μ−2 −ε

| My 1 − My 2 | ≤ Ck Cy M e Hs (h ) hs

| y 1 − y 2 | .

(D.10)

References [1] D. Ebling, M. Jaegle, M. Bartel, A. Jacquot, H. Böttner, Multiphysics simulation of thermoelectric systems for comparison with experimental device performance, J. Electron. Mater. 38 (7) (2009) 1456–1461. [2] Y.Y. Hsiao, W.C. Chang, S.L. Chen, A mathematic model of thermoelectric module with applications on waste heat recovery from automobile engine, Energy 35 (3) (2010) 1447–1454. [3] J.L. Pérez-Aparicio, R. Palma, R.L. Taylor, Finite element analysis and material sensitivity of Peltier thermoelectric cells coolers, Int. J. Heat Mass Transf. 55 (4) (2012) 1363–1374. [4] G.D. Mahan, Density variations in thermoelectrics, J. Appl. Phys. 87 (10) (2000) 7326–7332.

258


[5] J.L. Pérez-Aparicio, R.L. Taylor, D. Gavela, Finite element analysis of nonlinear fully coupled thermoelectric materials, Comput. Mech. 40 (1) (2007) 35–45. [6] L. Liu, A continuum theory of thermoelectric bodies and effective properties of thermoelectric composites, Int. J. Eng. Sci. 55 (2012) 35–53. [7] Y. Yang, S. Xie, F. Ma, J. Li, On the effective thermoelectric properties of layered heterogeneous medium, J. Appl. Phys. 111 (1) (2012) 013510. [8] S. Prudhomme, F. Pascal, J. Oden, A. Romkes, Review of a Priori Error Estimation for Discontinuous Galerkin Methods, Tech. Rep., TICAM, UTexas, 2000. [9] L. Noels, R. Radovitzky, A general discontinuous Galerkin method for finite hyperelasticity. Formulation and numerical applications, Int. J. Numer. Methods Eng. 68 (1) (2006) 64–97. [10] A. Romkes, S. Prudhomme, J. Oden, A priori error analyses of a stabilized discontinuous Galerkin method, Comput. Math. Appl. 46 (8) (2003) 1289–1311. [11] W.H. Reed, T. Hill, Triangular mesh methods for the neutron transport equation, Los Alamos Report LA-UR-73-479 http://www.osti.gov/scitech/ servlets/purl/4491151, 1973. [12] B. Cockburn, G.E. Karniadakis, C.W. Shu, The Development of Discontinuous Galerkin Methods, Springer, 2000. [13] J. Douglas, T. Dupont, Interior Penalty Procedures for Elliptic and Parabolic Galerkin Methods, Springer Berlin Heidelberg, Berlin, Heidelberg, ISBN 978-3-540-37550-0, 1976, pp. 207–216. [14] M.F. Wheeler, An elliptic collocation-finite element method with interior penalties, SIAM J. Numer. Anal. 15 (1) (1978) 152–161. [15] D.N. Arnold, F. Brezzi, B. Cockburn, L.D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (5) (2002) 1749–1779. [16] E.H. Georgoulis, Discontinuous Galerkin Methods on Shape-Regular and Anisotropic Meshes, D. Phil. Thesis, University of Oxford, 2003. [17] S. Yadav, A. Pani, E.J. Park, Superconvergent discontinuous Galerkin methods for nonlinear elliptic equations, Math. Comput. 82 (283) (2013) 1297–1335. [18] T. Gudi, N. Nataraj, A.K. Pani, hp-Discontinuous Galerkin methods for strongly nonlinear elliptic boundary value problems, Numer. Math. 109 (2) (2008) 233–268. [19] S. Sun, M.F. Wheeler, Discontinuous Galerkin methods for coupled flow and reactive transport problems, Appl. Numer. Math. 52 (2) (2005) 273–298. [20] X.P. Zheng, D.H. Liu, Y. Liu, Thermoelastic coupling problems caused by thermal contact resistance: a discontinuous Galerkin finite element approach, Sci. China, Phys. Mech. Astron. 54 (4) (2011) 666–674. [21] S. Kesavan, Topics in Functional Analysis and Applications, John Wiley & Sons, 1989. [22] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method: Solid Mechanics, vol. 2, Butterworth-Heinemann, 2000. [23] P. Houston, J. Robson, E. Süli, Discontinuous Galerkin finite element approximation of quasilinear elliptic boundary value problems I: the scalar case, IMA J. Numer. Anal. 25 (4) (2005) 726–749. [24] L. Homsi, Development of Non-Linear Electro-Thermo-Mechanical Discontinuous Galerkin Formulations, Ph.D. thesis, University of Liege, Belgium, 2017. [25] D. Gilbarg, N.S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer, 1983. [26] I. Babuška, M. Suri, The hp version of the finite element method with quasiuniform meshes, RAIRO-Modél. Math. Anal. Numér. 21 (2) (1987) 199–238. [27] M. Ainsworth, D. Kay, The approximation theory for the p-version finite element method and application to non-linear elliptic PDEs, Numer. Math. 82 (3) (1999) 351–388. [28] M. Ainsworth, D. Kay, Approximation theory for the hp-version finite element method and application to the non-linear Laplacian, Appl. Numer. Math. 34 (4) (2000) 329–344. [29] D.K.M. Ainsworth, The approximation theory for the p-version finite element method and application to the nonlinear elliptic PDEs, Numer. Math. 82 (1999) 351–388. [30] T. Gudi, Discontinuous Galerkin Methods for Nonlinear Elliptic Problems, Ph.D. thesis, Indian Institute of Technology, Bombay, 2006. [31] P. Hansbo, M.G. Larson, Discontinuous Galerkin methods for incompressible and nearly incompressible elasticity by Nitsche’s method, Comput. Methods Appl. Mech. Eng. 191 (17) (2002) 1895–1908. [32] P. Ciarlet, Conforming Finite Element Methods for Second-Order Problems, Chap. 3, SIAM, 2002, pp. 110–173.