FEM analysis for a thermoviscoelastic system

Finite element convergence analysis for the thermoviscoelastic Joule heating problem

Axel Målqvist Mathematical Sciences
Chalmers University of Technology and the University of Gothenburg
, SE-412 96 Göteborg
Sweden.
and  Tony Stillfjord
Februart 20, 2017
Abstract.

We consider a system of equations that model the temperature, electric potential and deformation of a thermoviscoelastic body. A typical application is a thermistor; an electrical component that can be used e.g. as a surge protector, temperature sensor or for very precise positioning. We introduce a full discretization based on standard finite elements in space and a semi-implicit Euler-type method in time. For this method we prove optimal convergence orders, i.e. second-order in space and first-order in time. The theoretical results are verified by several numerical experiments in two and three dimensions.

Key words and phrases:
Partial differential equations, thermoviscoelastic, Joule heating, thermistor, convergence analysis, finite elements
2010 Mathematics Subject Classification:
65M12, 65M60, 74D05, 74H15
Both authors were supported by the Swedish Research Council under grant 2015-04964

A. Målqvist]axel@chalmers.se

T. Stillfjord]tony.stillfjord@gu.se

1. Introduction

Consider the following system of coupled equations:

 ˙θ =Δθ+σ(θ)|∇ϕ|2−M:ϵ(˙u), (1) 0 =∇⋅(σ(θ)∇ϕ), (2) ¨u =∇⋅(Aϵ(˙u)+Bϵ(u)−Mθ)+f, (3)

with initial conditions

 θ(0,x)=θ0(x),u(0,x)=u0(x)and˙u(0,x)=v0(x),

over the convex polygonal or polyhedral domain with . Together with appropriate boundary conditions, to be specified later, these equations describe the evolution of the temperature , electric potential and deformation of a conducting body. Here , are constant tensors of order 4, describing the viscosity and elasticity of the body, and is a constant matrix describing the thermal expansion of the body. The vector consists of external forces and denotes the electrical conductivity, which here depends on the temperature. In addition, we have used the notation

 ϵ(u)=12(∇u+(∇u)T)

for the linearized strain tensor and for the Frobenius inner product.

The coupling of electricity and temperature through (1)–(2) is commonly known as Joule heating and is typically used to model thermistors, see e.g. [5, 9]. These are electrical components used for example as surge protectors or temperature sensors. The inclusion of thermoviscoelastic effects through (3) allows us to also model their use as actuators on the micro-scale, cf. [16].

We note that the Joule heating problem, both stationary and time-dependent, has been considered extensively in different contexts. For discussions on existence and uniqueness, see e.g. [2, 5, 6, 8, 9, 17, 18, 19, 24, 31] and the references therein. For the fully coupled, deformable problem the literature is less extensive. We refer mainly to [20] for the non-degenerate case that we consider here, with . See also [30] for the degenerate case where is allowed; this requires a more generalized solution concept.

However, to our knowledge there exists no numerical analysis for methods applied to the fully coupled case. Many authors have analyzed methods for similar problems. For example, [12] considers the quasi-static version where the -term is ignored, [1], [11] and [22] considers the non-deformable case, [13, 14] treat the purely thermoviscoelastic case (no ) with nonlinear constituent law, etc. Additionally, in the deformable case a common theme seems to be suboptimal convergence orders, i.e. errors of the form instead of .

The main contribution of this article is therefore an error analysis for a fully discrete discretization applied to the problem (1)–(3), which shows optimal convergence orders in both time and space. For the spatial discretization we consider standard finite elements, and for the temporal discretization a semi-implicit Euler-type method. Our approach also allows us to analyze e.g. the implicit Euler method, but the semi-implicit method benefits from a greatly decreased computational cost while the errors are comparable.

The central idea of our proof is to bound the errors in and in terms of the error in , in the spirit of [11] and [23]. The latter error then fulfills an equation similar to (1), to which we may apply a Grönwall inequality after properly handling the quadratic potential term. We note that we avoid any time step restrictions of the form by performing the analysis in two steps, where the first considers only the discretization in time, cf. [23]. Finally, in order to produce the error bound, we extend the concept of Ritz-Volterra projections for damped wave equations (see [25]) to the discrete and vector-valued viscoelasticity case.

For simplicity, we consider Dirichlet boundary conditions,

 θ(t,x)=0,ϕ(t,x)=ϕb(t,x)andu(t,x)=0

for and . This is a simplified case of the ideal situation with an arbitrary polygon and mixed boundary conditions, corresponding to where the body is clamped and insulated. As is well known (see e.g. [15]) the solutions to such a problem would typically suffer from a lack of regularity in the vicinity of re-entrant corners and boundary condition transitions, which leads to suboptimal convergence orders for finite-element based numerical methods. We therefore restrict ourselves to the simplified model, and will indicate possible generalizations by our numerical experiments.

A brief outline of the article is as follows. In Section 2 we write the problem on weak form and discretize it in both time and space. The assumptions on the data and solutions to the continuous problem are given in Section 3, where we also perform the error analysis. In Subsection 3.1, the time-discrete system is shown to be first-order convergent, and then the full discretization is shown to be second-order convergent to the time-discrete system in Subsection 3.2. These results are confirmed by the numerical experiments presented in Section 4, and conclusions and future work is summarized in Section 5.

2. Weak formulation and discretization

In order to present a weak formulation of the problem, we introduce the spaces

 V:=H10(Ω)⊂L2(Ω),and V:=H10(Ω)d⊂L2(Ω)d=:L2(Ω),

as well as the space of symmetric matrices,

 Q={ξ=(ξij)di,j=1⊂L2(Ω)d×d;ξji=ξij,1≤i,j≤d}.

The idea here is that and belong to , and . On , we have the inner product

 (ξ,ζ)Q:=∫Ωξ(x):ζ(x)dx=d∑i,j=1(ξij,ζij)L2(Ω).

which gives rise to the norm . To simplify some notation, we use the inner product

 (u,v)V=(ϵ(u),ϵ(v))Q

on instead of the usual one. The norm induced by this inner product is equivalent to by Korn’s inequality, see e.g. [10, Chapter III, Theorems 3.1, 3.3] and [27]. We will on several occasions make use also of the norm , which arises from the elasticity operator through

 ∥u∥2B=(Bϵ(u),ϵ(u))Q,

as well as the norm defined analogously for a small positive constant . Under Assumption 3.1 in the next section, both of these norms are equivalent to the -norm. In the following, we will omit the specification of and simply write or . Additionally, the - and -norms will both simply be denoted by and the corresponding inner products by , where no confusion can arise.

By multiplying the equations (1), (2) with the test function , Equation (3) with and then using Green’s formula we get

 (˙θ,χ)+(∇θ,∇χ) =(σ(θ)|∇ϕ|2,χ)−(M:ϵ(˙u),χ), (4) (σ(θ)∇ϕ,∇χ) =0, (5) (¨u,χ)+(Aϵ(˙u)+Bϵ(u),ϵ(χ))Q =(Mθ,ϵ(χ))Q+(f,χ), (6)

for all and , respectively. In (6), we have made use of the identity as well as the similar identities and . The latter two hold because we assume and to be symmetric; see Assumption 3.1 in the next section. Note also that we have omitted the time parameter here and in the original equation; both are supposed to hold for all times for a given .

We now discretize the time interval using a constant temporal step size , which results in the grid with and . We will abbreviate function evaluations at these times by sub-scripts, so that

 θn=θ(tn),ϕn=ϕ(tn),un=u(tn)andfn=f(tn).

The approximations of these solution values should belong to the same spaces as in the continuous case, and we will denote them by capital letters and superscripts:

 Θn≈θn,Φn≈ϕnandUn≈un.

Additionally, we denote by the first-order backward difference quotient, i.e.

 DtΘn=Θn−Θn−1k.

With this notation given, we now consider the following semi-implicit temporal discretization of Equations (1)–(3),

 DtΘn =ΔΘn+σ(Θn−1)|∇Φn−1|2−M:ϵ(DtUn−1), (7) 0 =∇⋅(σ(Θn)∇Φn), (8) D2tUn =∇⋅(Aϵ(DtUn)+Bϵ(Un)−MΘn)+fn, (9)

where , and its corresponding weak form,

 (DtΘn,χ)+(∇Θn,∇χ)=(σ(Θn−1)|∇Φn−1|2,χ)−(M:ϵ(DtUn−1),χ), (10) (σ(Θn)∇Φn,∇χ)=0, (11) (12)

for and for all and , respectively. The initial conditions are the same as in the continuous case: , and . (We use a fictitious point to define .) Note that this discretization results in a decoupling of the equations; we solve first for using (7) then use this to find from (8) and from (9). This implies a significant decrease in computational effort compared to the fully coupled case arising from e.g. the implicit Euler discretization.

For the spatial discretization, we introduce the finite element spaces and . These consist of continuous, piecewise linear functions with zero trace on , defined on a quasi-uniform mesh with mesh-width . Then the fully discrete problem we are interested in is given by

 (DtΘnh,χ)+(∇Θnh,∇χ)=(σ(Θn−1h)|∇Φn−1h|2,χ)−(M:ϵ(DtUn−1h),χ), (13) (σ(Θnh)∇Φnh,∇χ)=0, (14) (D2tUnh,χ)+(Aϵ(DtUnh)+Bϵ(Unh),ϵ(χ))Q=(MΘnh,ϵ(χ))Q+(fn,χ), (15)

for and for all and , respectively. Here, the approximations satisfy , and . (We assume that is defined on all of .) As initial conditions, we take , and , the Lagrangian interpolant of the exact initial condition.

Remark 2.1.

We assume the domain to be a convex polygon or polyhedron in order that the standard interpolation and regularity estimates for linear elliptic problems are satisfied, see [7, Section 3.2]. Similarly, the quasi-uniformity of the mesh guarantees that the standard inverse inequalities are satisfied. These are needed to handle the nonlinear potential term in (1), see [11, 23].

3. Error analysis

Our main goal is to estimate the errors , and . In order to do this, we will generalize the analysis of [23] (cf. also [11]) for the case with no deformation. This consists of first showing that the time-discrete approximations are -close to the solutions of the continuous system, and also proving that these approximations exhibit a certain regularity. The key part here is to express the error in the potential in terms of the error in the temperature, and then only working with the temperature equation. With the given regularity, the time-discrete and fully discrete approximations can then be compared and shown to be -close. The main problem here is the nonlinear term , which is handled in a two-step fashion: first using that to show that in fact and then using this to estimate in a stronger norm.

In our case, the temperature equation (1) contains the extra term , so our idea is to also bound the error in by the error in the temperature. Then we show that the approximations possess certain regularity, which may be used to also express the fully discrete deformation errors in terms of the fully discrete temperature errors. The key part in the latter step is to utilize the concept of Ritz-Volterra projections [25], which we here generalize to the vector-valued viscoelasticity case, as well as to discrete time.

Before we perform this extended analysis, we state the general assumptions on the given data. In these, as well as throughout the rest of the paper, denotes a generic constant independent of , and but possibly depending on , that may differ from line to line.

Assumption 3.1.

The viscosity and elasticity tensors and are symmetric, and both yield Lipschitz continuous and strongly coercive bilinear forms. That is,

 aijkl=ajikl=aklij,bijkl=bjikl=bklij,

and there are positive constants such that for all we have

 max((Aϵ(u),ϵ(v))Q,(Bϵ(u),ϵ(v))Q) ≤C1∥u∥V∥v∥Vand min((Aϵ(u),ϵ(u))Q,(Bϵ(u),ϵ(u))Q) ≥C2∥u∥2V.
Assumption 3.2.

The electrical conductivity belongs to and there are positive constants , and such that for all we have

 0<σmin≤σ(θ)≤σmaxand|σ′(θ)|≤σ′max.
Assumption 3.3.

The function , and is regular enough that

 ∥ϕb∥L∞(0,T;W2,12/5)+∥˙ϕb∥L2(0,T;H1)+∥∇ϕb∥L∞(0,T;L∞)≤C.

By [20], these assumptions guarantee the existence of a weak solution to the problem, i.e functions satisfying (4)–(6) with the time derivatives interpreted in a weak sense. Thus for example and . For optimal convergence orders more regularity is required, and explicit conditions on the data that guarantees such regularity is currently unknown. We therefore also make the following regularity assumption, where :

Assumption 3.4.

There exist solutions to (4)–(6) over the time interval which are regular enough that

 ∥θ∥L∞(0,T;H2)+∥˙θ∥L∞(0,T;L2)+∥˙θ∥L2(0,T;H2)+∥¨θ∥L1(0,T;L2) ≤C, ∥ϕ∥L∞(0,T;W2,12/5)+∥˙ϕ∥L2(0,T;H1)+∥ϕ∥L∞(0,T;W1,∞) ≤C, ∥˙u∥L∞(0,T;H2)+∥¨u∥L∞(0,T;H2)+∥u(3)∥L1(0,T;L2) ≤C

The assumptions on and are essentially the same as in the non-deformable situation given in [23], while the assumptions on and are new. We note that for the non-deformable case, the existence of solutions with similar regularity properties was shown in [11] when , with weak requirements on the initial values. In the general elliptic/parabolic case, the absence of reentrant corners in the convex domain makes such regularity plausible, see e.g. [15, Chapters 3,4] and [28, Chapter 19]. In the displacement equation the viscosity term acts as damping, and we expect regular solutions to be present also there, see e.g. [21]. We are not aware of any regularity results for the fully coupled system, but we note that our numerical experiments with smooth data suggest that Assumption 3.4 is satisfied in practice.

The following main theorem will be proved in the next two subsections:

Theorem 3.1.

Let Assumptions 3.1-3.4 be satisfied and let and be solutions to the equations (4)–(6) and (13)–(15), respectively. Then there are positive constants and such that if and we have for that

 ∥Θnh−θn∥+∥Φnh−ϕn∥+∥DtUnh−˙un∥≤C(h2+k),

and

 ∥Θnh−θn∥H1+∥Φnh−ϕn∥H1+∥DtUnh−˙un∥V≤C(h+k).

The constant is independent of , and , but may depend on the final time and the problem data.

To abbreviate expressions like the above in the following, we introduce

 enθ=Θn−θn,enϕ=Φn−ϕnandenu=Un−un

as well as

 enθ,h=Θnh−Θn,enϕ,h=Φnh−Φnandenu,h=Unh−Un.

3.1. The time-discrete case

We start by considering the semi-discrete case, and first provide a bound for in terms of .

Lemma 3.1.

Let Assumptions 3.1-3.4 be satisfied and let and be solutions to the equations (4)–(6) and (10)–(12), respectively. Then we have

 ∥Dtenu∥2+∥enu∥2V+kn∑j=1∥Dteju∥2V≤Ck2+Ckn∑j=1∥ejθ∥2,

for , with the constant independent of and .

Proof.

By equations (6) and (12), we see that the error satisfies

 (D2tenu,χ)+(Aϵ(Dtenu)+Bϵ(enu),ϵ(χ)) =(Menθ,ϵ(χ))+(¨u(tn)−D2tu(tn),χ) +(Aϵ(˙u(tn)−Dtu(tn)),ϵ(χ)) ≤C∥enθ∥∥χ∥V+Ck∥χ∥+Ck∥χ∥V

due to the regularity assumptions on . We note that for any sequence we have

 2(Dt2gn,Dtgn)≥Dt∥Dtgn∥2and2(Bϵ(gn),ϵ(Dtgn))≥Dt∥gn∥2B,

where is the norm induced by the inner product . Thus by choosing and using the Cauchy–Schwarz inequality as well as Young’s inequality, , we get

 Dt∥Dtenu∥2+2C2∥Dtenu∥V+Dt∥enu∥2B≤Ck2+C∥enθ∥2+C2∥Dtenu∥2V.

Canceling the final term, summing over and modifying the constants then yields

 ∥Dtenu∥2+kn∑j=1∥Dteju∥V+∥enu∥2B≤Ck2+Ckn∑j=1∥ejθ∥2,

and the Lemma follows from the equivalence between the - and -norms.

Theorem 3.2.

Let Assumptions 3.1-3.4 be satisfied and let and be solutions to the equations (1)–(3) and (7)–(9), respectively. Then there is a positive constant such that if then

 ∥enθ∥2H1+∥enϕ∥2H1+∥Dtenu∥2V≤Ck2,

for , with the constant independent of and . In addition, the approximations have the following regularity:

 ∥Θn∥2H2+∥DtΘn∥2+kn∑j=1∥DtΘj∥2H2 ≤C, ∥Φn∥W2,12/5+∥Φn∥W1,∞ ≤C, ∥DtUn∥2H2+∥D2tUn∥2V+kn∑j=1∥D2tUj∥2H2 ≤C.
Proof.

To begin with, we see that the error satisfies

 −∇⋅(σ(Θn)∇enϕ))=∇⋅((σ(Θn)−σ(θn))∇ϕn).

Multiplying this equation by and integrating directly yields

 ∥∇enϕ∥2≤C∥∇ϕn∥L∞∥enθ∥∥∇enϕ∥,

so that

 ∥∇enϕ∥≤C∥enθ∥ (16)

by the regularity assumptions. This inequality for corresponds to Lemma 3.1 for . Further, we see that the error satisfies

 Dtenθ−Δenθ =(σ(Θn−1)−σ(θn−1))|∇ϕn−1|2+σ(Θn−1)(∇Φn−1+∇ϕn−1)⋅∇en−1ϕ (17) −M:ϵ(Dten−1u)+Rnθ,

where

 Rnθ =(σ(θn−1)−σ(θn))|∇ϕn−1|2+σ(θn)(∇ϕn−1+∇ϕn)⋅(∇ϕn−1−∇ϕn) +M:ϵ(˙un−˙un−1)+M:ϵ(˙un−1−Dtun−1).

is bounded by , again by the regularity assumptions. After multiplying by and integrating, we therefore get

 Dt∥enθ∥2+2∥∇enθ∥2 ≤C∥en−1θ∥∥enθ∥∥∇ϕn−1∥L∞+(M:ϵ(Dten−1u),enθ)+Ck∥enθ∥ (18) +(σ(Θn−1)(∇Φn−1+∇ϕn−1)enθ,∇en−1ϕ).

The last term of this expression can be shown to be bounded by , see [23, p.627], and for the second term we observe that for a generic ,

 (M:(∇u),χ)L2=(∇u,Mχ)Q=−(u,∇⋅(Mχ))L2=−(u,M∇χ)L2.

As a completely analogous calculation holds also for and is symmetric, we thus have

 (M:ϵ(u),χ)=−(u,M∇χ)≤C∥u∥∥∇χ∥. (19)

This implies that (18) reduces to

 Dt∥enθ∥2+2∥∇enθ∥2≤C(k2+∥en−1θ∥2+∥enθ∥2+∥en−1ϕ∥2H1+∥Dten−1u∥2)+∥∇enθ∥2.

Canceling the last term, summing up and using Equation (16) and Lemma 3.1 thus yields

 ∥enθ∥2+kn∑j=1∥∇ejθ∥2≤Ck2+Ckn∑j=1∥ejθ∥2.

Under the step size restriction , we can eliminate the last term of the sum. An application of Grönwall’s lemma then shows that the left-hand side is bounded by . Using Equation (16) and Lemma 3.1 again, we see that in fact

 ∥enθ∥2+kn∑j=1∥∇ejθ∥2+∥∇enϕ∥2+∥Dtenu∥2+∥enu∥2V+kn∑j=1∥Dteju∥2V≤Ck2

From these preliminary bounds, we may deduce the desired regularity of and and then test (17) with to acquire

 ∥enθ∥2H1+kn∑j=1∥Δejθ∥2≤Ck2.

For details, we refer to [23, Theorem 3.1]. Let us instead investigate the remaining questions of the regularity of and the pointwise bound for in the -norm. By the defining equation, we have that

 ∇⋅(Aϵ(Dtenu)+Bϵ(enu)) =D2tenu+∇⋅(MΘn)+Dt2u(tn)−¨u(tn) (20) +∇⋅(Aϵ(Dtu(tn)−˙u(tn))),

where the right-hand side is in since . Let us denote it by . Then we can rewrite the previous equation as

 ∇⋅(Aϵ(Dtenu)+kBϵ(Dtenu))=gn+∇⋅(Bϵ(en−1u)).

Now since both and induce bounded and coercive inner products on , we see that

 ∥Dtenu∥2H2 ≤C∥∇⋅(Aϵ(Dtenu)+kBϵ(Dtenu))∥2 ≤C∥gn∥2+C∥en−1u∥2H2

But since , we can estimate the second term by Cauchy–Schwarz as

 ∥en−1u∥2H2≤kn−1∑j=1∥Dteju∥2H2.

An application of Grönwall’s lemma thus shows that

 ∥Dtenu∥H2≤C,

which also implies that , and are all in . We may now multiply (20) by and integrate to get

 (Dtϵ(Dtenu),(A+kB)ϵ(Dtenu))+∥∇⋅((A+kB)ϵ(Dtenu))∥2≤C∥enθ∥2H1+C∥en−1θ∥2H2,

where we have used the Cauchy-Schwarz and Young inequalities and canceled a term . The first term on the left-hand side can be estimated from below by , so summing up and using the equivalence of the - and -norms, we get

 ∥Dtenu∥2V+kn∑j=1∥Dteju∥2H2≤Ckn−1∑j=1∥ejθ∥2H1+Ckn−1∑j=1∥Dteju∥2H2.

But the first term in the right-hand side is bounded by and in the second term we may again use that . Defining

 wn=∥Dtenu∥2V+kn∑j=1∥Dteju∥2H2,

we thus have

 wn≤Ck2+Ckn−1∑j=1wj,

and an application of Grönwall’s lemma shows that . This yields the final desired error bound, and additionally shows that , which implies the stated regularity for .

3.2. The fully discrete case

We now turn to the fully discretized case and first prove an analogue to Lemma 3.1.

Lemma 3.2.

Let Assumptions 3.1-3.4 be satisfied and and be solutions to equations (10)–(12) and (13)–(15), respectively. Then there is a positive constant such that if we have for that

 ∥enu,h∥2+∥Dtenu,h∥2 ≤Ch4+Ckn∑j=1∥ejθ,h∥2and ∥enu,h∥2V+kn∑j=1∥Dteju,h∥2V ≤Ch2+Ckn∑j=1∥ejθ,h∥2,

with the constant independent of , and .

Remark 3.1.

In the case of a first-order equation, one would typically first add and subtract the Ritz projection of in order to work only in the finite element space. This approach is viable also in the second-order case, if one defines the Ritz projection using the inner product. We refer to [29] for the scalar-valued case. However, we choose to instead work with a Ritz-Volterra projection, see [25] for the scalar-valued case. Such a projection takes both the - and -terms into account simultaneously, i.e. it is a projection of -functions rather than of elements in . In the present situation, we need of course to consider a discretized version, but it nevertheless simplifies matters.

Proof.

Subtracting (12) from (15), we see that

 (D2tenu,h,χ)+(Aϵ(Dtenu,h)+Bϵ(enu,h),ϵ(χ))=(Menθ,h,ϵ(χ))

for all . Now let , where

 ηn=Unh−Wn∈Shand