Finite element convergence analysis for the thermoviscoelastic Joule heating problem

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.


Introduction
Consider the following system of coupled equations: with initial conditions over the convex polygonal or polyhedral domain Ω ⊂ R d with d ≤ 3. Together with appropriate boundary conditions, to be specified later, these equations describe the evolution of the temperature θ , electric potential φ and deformation u of a conducting body. Here A, B are constant tensors of order 4, describing the viscosity and elasticity of the body, and M is a constant matrix describing the thermal expansion of the body. The vector f consists of external forces and σ (θ) denotes the electrical conductivity, which here depends on the temperature. In addition, we have used the notation for the linearized strain tensor and : for the Frobenius inner product. The coupling of electricity and temperature through (1.1)-(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 (1.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]23,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 σ ≥ σ min > 0. See also [30] for the degenerate case where σ = 0 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,24] 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 O(h + k) instead of O(h 2 + k).
The main contribution of this article is therefore an error analysis for a fully discrete discretization applied to the problem (1.1)-(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 semiimplicit 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 φ andu in terms of the error in θ , in the spirit of [11,22]. The latter error then fulfills an equation similar to (1.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 k ≤ h d/r by performing the analysis in two steps, where the first considers only the discretization in time, cf. [22]. Finally, in order to produce theu 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, for t ∈ [0, T ] and x ∈ ∂Ω. 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 Sect. 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 Sect. 3, where we also perform the error analysis. In Sect. 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 Sect. 3.2. These results are confirmed by the numerical experiments presented in Sect. 4, and conclusions and future work is summarized in Sect. 5.

Weak formulation and discretization
In order to present a weak formulation of the problem, we introduce the spaces as well as the space of symmetric matrices, The idea here is that θ and φ − φ b belong to V , u ∈ V and ε(u) ∈ Q. On Q, we have the inner product which gives rise to the norm · Q . To simplify some notation, we use the inner product on V instead of the usual one. The norm · V induced by this inner product is equivalent to · H 1 (Ω) d 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 · B , which arises from the elasticity operator through as well as the norm · A+kB defined analogously for a small positive constant k. Under Assumption 3.1 in the next section, both of these norms are equivalent to the V -norm. In the following, we will omit the specification of Ω and simply write L 2 or L 2 . Additionally, the L 2 -and L 2 -norms will both simply be denoted by · and the corresponding inner products by (·, ·), where no confusion can arise. By multiplying the Eqs. (1.1), (1.2) with the test function χ ∈ V , Eq. (1.3) with χ ∈ V and then using Green's formula we get for all χ ∈ V and χ ∈ V , respectively. In (2.3), we have made use of the identity (ε(u), ∇v) = (ε(u), ε(v)) as well as the similar identities (Aε(u), ∇v) = (Aε(u), ε(v)) and (Bε(u), ∇v) = (Bε(u), ε(v)). The latter two hold because we assume A and B 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 t ∈ (0, T ] for a given T . We now discretize the time interval [0, T ] using a constant temporal step size k, which results in the grid t n = nk with n = 1, 2, . . . , N and N k = T . We will abbreviate function evaluations at these times by sub-scripts, so that θ n = θ(t n ), φ n = φ(t n ), u n = u(t n ) and f n = f (t n ).
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 ≈ φ n and U n ≈ u n . Additionally, we denote by D t the first-order backward difference quotient, i.e.
With this notation given, we now consider the following semi-implicit temporal discretization of Eqs. (1.1)-(1.3), where D 2 t = D t D t , and its corresponding weak form, for n = 1, . . . , N and for all χ ∈ S h and χ ∈ S h , respectively. The initial conditions are the same as in the continuous case: Θ 0 = θ 0 , U 0 = u 0 and D t U 0 = v 0 . (We use a fictitious point U −1 to define D t U 0 .) Note that this discretization results in a decoupling of the equations; we solve first for Θ n using (2.4) then use this to find Φ n from (2.5) and U n from (2.6). 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 S h ⊂ V and S h ⊂ V . These consist of continuous, piecewise linear functions with zero trace on ∂Ω, defined on a quasi-uniform mesh with mesh-width h. Then the fully discrete problem we are interested in is given by for n = 1, . . . , N and for all χ ∈ S h and χ ∈ S h , respectively. Here, the approxima- defined on all of Ω.) As initial conditions, we take U 0 h = 0, D t U 0 h = 0 and Θ 0 h = I h θ 0 , 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.1), see [11,22].

Error analysis
Our main goal is to estimate the errors Θ n h − θ n , Φ n h − φ n and U n h − u n . In order to do this, we will generalize the analysis of [22] (cf. also [11]) for the case with no deformation. This consists of first showing that the time-discrete approximations are O(k)-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 O(h 2 )-close. The main problem here is the nonlinear term σ (θ)|∇φ| 2 , which is handled in a two-step fashion: first using Ch and then using this to estimate ∇(Φ n h − Φ n ) in a stronger norm. In our case, the temperature Eq. (1.1) contains the extra term M : ε(u), so our idea is to also bound the error inu by the error in the temperature. Then we show that the approximations U n 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, C denotes a generic constant independent of k, h and n but possibly depending on T , that may differ from line to line.

Assumption 3.1
The viscosity and elasticity tensors A = (a i jkl ) and B = (b i jkl ) are symmetric, and both yield Lipschitz continuous and strongly coercive bilinear forms. That is, and there are positive constants

Assumption 3.2
The electrical conductivity σ belongs to C 1 (R) and there are positive constants σ min , σ max and σ max such that for all θ ≥ 0 we have By [20], these assumptions guarantee the existence of a weak solution to the problem, i.e functions (θ, φ, u) satisfying (2.1)-(2.3) with the time derivatives interpreted in a weak sense. Thus for example θ ∈ L 2 (0, T ; V ) andθ ∈ L 2 (0, T ; V ) . 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 H 2 = H 2 (Ω) d : The assumptions on θ and φ are essentially the same as in the non-deformable situation given in [22], while the assumptions on u and f are new. We note that for the nondeformable case, the existence of solutions with similar regularity properties was shown in [11] when d ≤ 2, 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:

The constant C is independent of k, h and n, but may depend on the final time T = N k and the problem data.
To abbreviate expressions like the above in the following, we introduce e n θ = Θ n − θ n , e n φ = Φ n − φ n and e n u = U n − u n as well as

The time-discrete case
We start by considering the semi-discrete case, and first provide a bound for D t e n u in terms of e n θ .
due to the regularity assumptions on u. We note that for any sequence {g n } we have where · B is the norm induced by the inner product Bε(·), ε(·) . Thus by choosing χ = D t e n u and using the Cauchy-Schwarz inequality as well as Young's inequality, Proof To begin with, we see that the error e n φ satisfies Multiplying this equation by e n φ and integrating directly yields by the regularity assumptions. This inequality for e n φ corresponds to Lemma 3.1 for e n u . Further, we see that the error e n θ satisfies where R n θ = σ (θ n−1 ) − σ (θ n ) |∇φ n−1 | 2 + σ (θ n ) ∇φ n−1 + ∇φ n · ∇φ n−1 − ∇φ n + M : ε(u n −u n−1 ) + M : ε(u n−1 − D t u n−1 ).
is bounded by R n θ ≤ Ck, again by the regularity assumptions. After multiplying by e n θ and integrating, we therefore get The last term of this expression can be shown to be bounded by C( e n θ 2 + e φ 2 H 1 ), see [22, p. 627], and for the second term we observe that for a generic u ∈ V , As a completely analogous calculation holds also for (∇u) T and M is symmetric, we thus have Canceling the last term, summing up and using Eq. From these preliminary bounds, we may deduce the desired regularity of Θ n and Φ n and then test (3.2) with −Δe n θ to acquire e n θ 2 For details, we refer to [22,Theorem 3.1]. Let us instead investigate the remaining questions of the regularity of U n and the pointwise bound for D t e n u in the V -norm. By the defining equation, we have that ∇ · Aε(D t e n u ) + Bε(e n u ) = D 2 t e n u + ∇ · MΘ n + D t 2 u(t n ) −ü(t n ) +∇ · Aε(D t u(t n ) −u(t n )) , (3.5) where the right-hand side is in L 2 since D 2 t e n u ≤ k −1 ( D t e n u + D t e n−1 u ) ≤ C. Let us denote it by g n . Then we can rewrite the previous equation as ∇ · Aε(D t e n u ) + kBε(D t e n u ) = g n + ∇ · Bε(e n−1 u ) . An application of Grönwall's lemma thus shows that D t e n u H 2 ≤ C, which also implies that e n u , U n and D t U n are all in H 2 . We may now multiply (3.5) by ∇ · (A + kB)ε(D t e n u ) and integrate to get But the first term in the right-hand side is bounded by Ck 2 and in the second term we may again use that D t e j u 2 ≤ C, which implies the stated regularity for U n .

The fully discrete case
We now turn to the fully discretized case and first prove an analogue to Lemma 3.1. with the constant C independent of k, h and n. Remark 3.1 In the case of a first-order equation, one would typically first add and subtract the Ritz projection of e n u 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 (Aε(·), ε(·)) 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 Aand B-terms into account simultaneously, i.e. it is a projection of C 1 (0, T ; V )-functions rather than of elements in V . In the present situation, we need of course to consider a discretized version, but it nevertheless simplifies matters.
Proof Subtracting (2.9) from (2.12), we see that for all χ ∈ S h . Now let e n u,h = η n + ρ n , where η n = U n h − W n ∈ S h and ρ n = W n − U n , with the discrete Ritz-Volterra projection W n of U n satisfying W 0 = U 0 = 0 and for all χ ∈ S h . We note that Eq. (3.6) may also be stated as Aε D t ρ n + Bε(ρ n ), ε(χ ) = 0, and that since D t U 0 = 0, also D t W 0 = 0. Additionally, we need the Ritz projection R h given by the viscosity term. For a generic u ∈ V , this is defined by for all χ ∈ S h , and we have the inequality We start by estimating the V -norms of D t ρ n and ρ n . To this end, we observe that for a generic u, we have We therefore take ϕ ∈ C ∞ 0 (Ω) d with ϕ = 1 and let Ψ ∈ V be the solution to Then where the last term is bounded by Moreover, since D t W n ∈ S h , the first term is bounded by By expressing ρ n in terms of D t ρ j and noting that ρ 0 = 0, we thus have and under the step size restriction Ck < 1 we can eliminate the last term of the sum and apply Grönwall's lemma. This shows that By using the regularity shown in Theorem 3.2 and then summing over n, we see that Using these bounds we may now estimate ρ also in the L 2 -norm, by instead letting Ψ ∈ V be the solution to Then as before, For R 4 , we note that Ψ H 2 ≤ C ϕ ≤ C, so that by using integration by parts and observing that both ρ n and Ψ are zero on ∂Ω we get, Hence similarly to the calculation for the V -norm, Grönwall's lemma implies that To bound η n , we also need a bound on the second derivative of ρ n . For this, we apply where R φ contains terms related to the potential φ. Choosing χ = e n h , we know from [22] that We additionally know that e 0 Hence ifCh 5/2 ≤ 1 we have that e n h ≤ h 1/2 . Thus by induction e n h ≤ h 1/2 holds for all n such that 0 ≤ n ≤ N . But then also the other calculations just performed are valid for 1 ≤ n ≤ N , so in fact e n h ≤ h 3/2 . This preliminary bound may be used as in [22, p. 631] to show e n φ,h ≤ Ch and to improve the bound of the quadratic potential term to Hence, This proves e n θ,h ≤ Ch 2 , and from [22] we find e n φ,h + h e n φ,h H 1 ≤ Ch 2 . Applying Lemma 3.2 gives D t e n u,h ≤ Ch 2 . Finally, by inverse inequalities we find also that e n θ,h H 1 + D t e n u,h V ≤ Ch.
Proof (of Theorem 3.1) This follows directly from Theorem 3.2 and Theorem 3.3 upon observing that, e.g.
where the last term is bounded in the proper way due to the regularity assumptions on the solution to the continuous system.

Numerical experiments
We have implemented both the method based on (2.10)-(2.12) and the corresponding fully implicit method based on implicit Euler, using FEniCS (see e.g. [4,26]). These implementations were then used to verify our theoretical results by applying them to the following test examples.

Problem 1
First consider the two-dimensional problem with Ω = (0, 1) 2 , M = I , f = [0, 0] T and the viscosity and elasticity tensors given in Voigt notation by We take the electrical conductivity to be given by which has a rather steep slope close to θ = 2. The initial conditions are given by θ 0 (x, y) = 0 and u 0 (x, y) = v 0 (x, y) = [0, 0] T . These functions also define the Dirichlet boundary conditions for θ and u, while for φ they are given by φ b (x, y) = 5(1 − x). We discretize Ω by first subdividing it into squares and then dividing each square into four triangles. With N x squares in each dimension, each triangle has diameter h = 1/N x and the full grid has 4N 2 x triangles. We take N x ∈ {4, 8, 16, 32, 64}. Since the error should be O(h 2 + k), we choose the number of time steps to be N t = N 2 x /2. With the final time T = 1, this gives k = 2h 2 . We emphasize here that the time steps could be taken much larger than this, but illustrating the error is then less straightforward. Finally, because the exact solution of the problem is not available we cannot compute the exact errors. Instead, we compare the different approximations to a reference approximation computed by the implicit Euler scheme with N x = 128 and N t = 8192. Figure 1 shows the errors for the different discretizations on a logarithmic scale, for both the semi-implicit method (left) and the method based on implicit Euler (right). These clearly exhibit the expected error behaviour predicted by Theorem 3.3, except for the first points where the grid is very coarse. We also note that the errors are very similar in size, which means that the semi-implicit method is much more efficient. A peculiar effect in this case is that the semi-implicit errors in θ and φ are actually less than the implicit Euler errors, though this does not hold for the error in u.

Problem 2
In the second experiment, we investigated the influence of the viscosity on the errors.
To this end, we employ the same data as presented in Sect. 4.1 except for the viscosity operator which we set to A = γ ⎡ ⎣ 1 1 0 1 1 0 0 0 1 ⎤ ⎦ (in Voigt notation). In this case, we used N x ∈ {4, 8, 16, 32} with N t = N 2 x /4 and took N x = 64, N t = 1024 for the reference approximation. We only used the semi-implicit scheme here. The first observation is that varying γ has essentially no effect on the errors in θ and φ. This is to be expected, as the influence of u on θ is not so large. We therefore omit the plots of these errors, and instead present the error in u for different values of γ in Fig. 2.
We observe that the error clearly increases as γ is decreased, which is to be expected. Indeed, an inspection of the convergence proof indicates that the L 2 -error should be inversely proportional to the coercivity constant of A, and thus also of γ . This is, however, in the worst case. In the current situation, Fig. 2 indicates that even γ = 0 would be perfectly feasible, though smaller step sizes might be necessary to enter the asymptotic regime.

Problem 3
For our last numerical experiment, we consider a 3D problem arising from an engineering application, inspired by [16,17]. We let Ω be as in Fig. 3, which also shows a typical spatial tetrahedral discretization. This represents a micro-electro-mechanical system (MEMS) used for precise positioning on small scales. When an electric current is passed through the device from the upper-left connector to the lower-left connector, it heats up. This causes a deformation, which due to the asymmetrical design of the component makes the tip move downwards.
We employ homogeneous Neumann boundary conditions everywhere except for at the left-most edge of the two connectors. These correspond to the component being insulated and stress-free. On the left-most edge we choose the Dirichlet boundary conditions corresponding to the component being clamped and having a potential difference applied between the two connectors. The equations, including physical constants, are Here, ρ denotes the density, c the specific heat capacity, K = kI the thermal conductivity matrix, M = mI the thermal expansion matrix and σ the electrical conductivity. Additionally, θ indicates the deviation from the ambient temperature Θ 0 = 293.15 K. We choose the elasticity and viscosity operators to be given on Lamé parameter form: Aε(u) = 2η 1 ε(u) + η 2 tr ε(u)I and Bε(u) = 2με(u) + λ tr ε(u)I, are given in terms of Poisson's ratio ν and Young's modulus E, and η 1 , η 2 are corresponding viscosity parameters. Here, tr denotes the trace of a matrix; tr τ = τ 11 + τ 22 .
The parameter values we have used, similar to the material properties of silicon, are listed in Table 1. In addition to this, we take f = [0, 0, 0] T and choose the electrical conductivity as σ (θ) = 38 · 10 6 27 3000 + 550 We solve the problem until the time T = 0.1 using the semi-implicit method for different spatial and temporal discretizations. The maximum sizes h of the tetrahedrons that were used and the corresponding number of vertices are listed in Table 2. The time steps were again taken proportional to h 2 but modified slightly to yield an integer number of steps. Since the temporal grids thus generated are not refinements of each other, we measured the error as the sum of the errors at only the points t j = j · 10 −2 for j = 1, . . . , 10. These errors are listed in Table 2, and also plotted in Fig. 4. While we cannot apply Theorem 3.3 directly, due to the mixed boundary conditions and the non-convexity of the domain, we observe that we still acquire almost O(h 2 + k) convergence. The curves wiggle because k = Ch 2 is only approximately satisfied, and the different magnitudes of the errors reflect the relative sizes of the solution components. The larger error in θ for the coarsest mesh indicates that it violates either the k < k 0 or h < h 0 mesh size limitations.   The approximation to the solution of the problem defined in Sect. 4.3 at t = T and with the finest spatial and temporal discretization. In the right-most plot, the grid has been deformed according to the computed displacement and then super-imposed over the original mesh to illustrate the deformation. We note that the grid is never deformed in the actual computations (this figure is in color in the electronic version of the article) Finally, Fig. 5 shows the approximations Θ N h , Φ N h and U N h at T , viewed from the side. At this point in time the solutions have just reached their steady state, and we see that the body deforms in the expected fashion.

Conclusions and outlook
We have presented a fully discrete numerical method for the fully coupled thermoviscoelastic thermistor problem (1.1)-(1.3) and proved optimal convergence orders in both space and time. These theoretical results are validated by experimental results.
We reiterate that mixed boundary conditions and re-entrant corners might lead to order reductions. In that case an adaptive mesh refinement strategy may be used, which requires a good a posteriori error estimate. It is possible that the ideas in [3] regarding this can be extended to the present, deformable case.
As illustrated by Sect. 4.3, a typical thermistor is not convex, so a further item that could be improved in the analysis is therefore the shape of the computational domain itself. In this direction we note that the stationary version of the non-deformable problem has been studied in [17,19] for very general domains. It is our ambition to extend these ideas to the time-dependent deformable case in the future.
Finally, a similar analysis would apply also for higher-order methods both in time and space. See e.g. [24] for a Crank-Nicolson-approach to the non-deformable Joule heating problem. However, such an analysis would require extra regularity assumptions that are unfeasible in real-world engineering applications.