Finite element convergence analysis for the thermoviscoelastic Joule heating problem

BIT Numer Math (2017) 57:787–810 DOI 10.1007/s10543-017-0653-1 Finite element convergence analysis for the thermoviscoelastic Joule heating problem Axel Målqvist1 · Tony Stillfjord1 Received: 2 November 2016 / Accepted: 2 March 2017 / Published online: 15 March 2017 © The Author(s) 2017. This article is published with open access at Springerlink.com 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 firstorder in time. The theoretical results are verified by several numerical experiments in two and three dimensions. Keywords Partial differential equations · Thermoviscoelastic · Joule heating · Thermistor · Convergence analysis · Finite elements Mathematics Subject Classification 65M12 · 65M60 · 74D05 · 74H15 Communicated by Rolf Stenberg. Both authors were supported by the Swedish Research Council under Grant 2015-04964. B Tony Stillfjord Axel Målqvist 1 Mathematical Sciences, Chalmers University of Technology and the University of Gothenburg, 412 96 Göteborg, Sweden 123 788 A. Målqvist, T. Stillfjord 1 Introduction Consider the following system of coupled equations: θ̇ = Δθ + σ (θ )|∇φ|2 − M : ε(u̇),   0 = ∇ · σ (θ )∇φ ,   ü = ∇ · Aε(u̇) + Bε(u) − Mθ + f, (1.1) (1.2) (1.3) with initial conditions θ (0, x) = θ0 (x), u(0, x) = u 0 (x) and u̇(0, x) = v0 (x), over the convex polygonal or polyhedral domain Ω ⊂ Rd 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 ε(u) =  1 ∇u + (∇u)T 2 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–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 123 Finite element convergence analysis for the… 789 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 φ and u̇ 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 the u̇ 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) and u(t, x) = 0 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. 2 Weak formulation and discretization In order to present a weak formulation of the problem, we introduce the spaces V := H01 (Ω) ⊂ L 2 (Ω), and V := H01 (Ω)d ⊂ L 2 (Ω)d =: L2 (Ω), as well as the space of symmetric matrices,   Q = ξ = (ξi j )i,d j=1 ⊂ L 2 (Ω)d×d ; ξ ji = ξi j , 1 ≤ i, j ≤ d . 123 790 A. Målqvist, T. Stillfjord The idea here is that θ and φ − φb belong to V , u ∈ V and ε(u) ∈ Q. On Q, we have the inner product  (ξ, ζ ) Q := Ω ξ(x) : ζ (x) dx = d    ξi j , ζi j L 2 (Ω) . i, j=1 which gives rise to the norm · Q . To simplify s (...truncated)