# Variational h-adaption for coupled thermomechanical problems

Rohit Pethe, Thomas Heuzé, Laurent Stainier

► To cite this version:

Rohit Pethe, Thomas Heuzé, Laurent Stainier. Variational h-adaption for coupled thermomechanical problems. Engineering Computations, 2019, 10.1108/EC-05-2019-0243. hal-02379535

HAL Id: hal-02379535

https://hal.science/hal-02379535

Submitted on 25 Nov 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Variational h-adaption for coupled thermo-mechanical problems

Rohit Pethe¹, Thomas Heuzé¹, and Laurent Stainier¹

¹Research Institute in Civil and Mechanical Engineering (GeM, UMR 6183 CNRS/ECN/UN), École Centrale de Nantes, 1 rue de la Noë, F-44321 Nantes, France,

## Abstract

In this paper, a variational mesh h-adaption approach is presented for strongly coupled thermo-mechanical problems. The mesh is adapted by local subdivision controlled by an energy criterion. Thermal and thermo-mechanical problems are of interest here. In particular, steady and transient purely thermal problems, transient strongly coupled thermo-elasticity and thermo-plasticity problems are investigated. Different test cases are performed to test the robustness of the algorithm for the problems listed above. Since the algorithm is based on a set of tolerance parameters, parametric analyses and a study of their respective influence on the mesh adaption is carried out. This detailed analysis is performed on uni-dimensional problems, and a final example is provided in two dimensions.

**Keywords** mesh adaption, variational approach, thermo-mechanics, nonlinear coupled problems, multi-physics

## 1 Introduction

In a number of transient problems (purely thermal, mechanical or thermo-mechanical), zones of high gradients of fields of interest evolve with time and loading. It is therefore interesting to use a dynamical mesh adaption algorithm to capture the solution in zones of high gradients in order to maintain the required precision. Many methods of mesh adaption have been proposed in the literature based on error-estimation. In these methods, the strategy is to adapt the mesh to minimize an error estimate, typically an upper bound, among all meshes of fixed size; or by recursive application of local refinement steps (Verfürth 1996) (Ainsworth and Oden 2000). But these methods have certain limitations. Rigorous estimates can be derived for linear constitutive models (for example elasticity), but it becomes more complex when non-linear constitutive models are used. Moreover, admissible fields need to be reconstructed (Ladeveze, Pelle, and Rougeot 1991) (Zienkiewicz and Zhu 1987). In addition, standard error bounds require a certain regularity of the solution for their validity (Ciarlet 1988). Therefore, it can be difficult and costly to use this approach for complex problems involving non-linear constitutive models and/or large deformation. In addition, methods based on global remeshing of the domain of interest require to transfer internal variables between meshes, which can lead to artificial diffusion of the latter unless specific methods are used (Barthold, Schmidt, and Stein 1998) (Brancherie, Villon, and Ibrahimbegovic 2008).

Variational formulations allow us to express finite element problems as problems of minimization (or maximization) of an energy-like potential. This holds true for non-linear problems as well (Dacorogna 1989). In some cases, the energy functional is evident, whereas, in some cases it
needs careful formulation. For instance, in inelastic problems and dynamical problems, minimum principles can be obtained by careful time discretization (Radovitzky and Ortiz 1999)(Ortiz and Stainier 1999)(Yang, Stainier, and Ortiz 2006)(Ortiz and Repetto 1999). In these cases, the energy-like functional is incremental and incorporates the free energy, inertia and kinetics of material.

An alternative approach of mesh adaption for purely mechanical problems was recently proposed (Mosler and Ortiz 2007)(Mosler and Ortiz 2009), based on the variational approach of (Ortiz and Stainier 1999). This technique uses an error indicator rather than an error estimator. In a variational approach, an energy-like potential is to be minimized (or maximized), and the gain in this scalar value associated to a given mesh adaption indicates the level of approximation, following the minimum (or maximum) criterion. No error estimates are then used at any stage of the algorithm. It allows mesh adaption in presence of large deformations and non-linear constitutive behavior. In addition, it was shown in (Mosler and Ortiz 2007) that variational *h*-adaption could be combined with variational *r*-adaption, at least for hyperelastic behavior. Indeed, *r*-adaption would involve remapping in the presence of internal variables, and was not considered by these authors for dissipative behaviors. Recently, *hp*-adaptive energy minimisation has also been treated in (Houston and Wihler 2016) for linear problems, the *h*- part of the adaptive process being close to former ones (Mosler and Ortiz 2007)(Mosler and Ortiz 2009). In these, the authors addressed isothermal, steady state mechanical problems. The extension of the algorithm to thermo-mechanical problems raises some additional difficulties, which are addressed in the present work.

In this work, an *h*-adaption algorithm for problems in multi-physics is developed. The variational energy-like potential value is used to construct an error indicator and the variational principle itself drives mesh refinement and coarsening. Similar to (Mosler and Ortiz 2007), the algorithm presented is also based on a variational approach (Ortiz and Stainier 1999) but with an extended functional that admits a saddle point (Yang, Stainier, and Ortiz 2006; Stainier 2013). This algorithm now accounts for heat conduction, transient thermal and thermo-mechanical coupling effects. For problems in multiphysics, the different physics have different temporal and spatial scales. Different meshes are used for each physics in order to account for its own spatial scale, these multiple meshes are adapted sequentially, the solution scheme relying on staggered algorithms (Farhat, Park, and Y. 1991)(Armero and Simo 1992). This permits to accurately capture the different spatial scales associated with each physics while maintaining the cost effectiveness of the approach. For constitutive behaviors involving internal variables, these are stored at integration points of each mesh, and no transfer (remapping) is necessary. Transfers between meshes only involve interpolation of external fields (displacement and temperature in this case). Though this paper focuses on 1-D problems in order to better illustrate and analyze the method, this algorithm can easily be extended to 2-D and 3-D problems using different subdivision schemes such as longest edge propagation path (LEPP) bisection algorithm of Rivara (Rivara 1991)(Rivara and Levin 1992)(Bänsch 1991)(Rivara and Inostroza 1997)(Rivara 1997). In order to demonstrate this, first results from our ongoing work in 2-D is shown on a transient thermal test case. The edge bisection technique is used for refinement and coarsening the mesh. It enables to keep the same mesh topology so that a simple interpolation can be used between data located at integration points of the two meshes, avoiding a costly and diffusive projection of fields.

The structure of the article is as follows. In the second section, the variational formulation which is the base for our mesh adaption strategies is first presented. Then the mesh adaption algorithm is introduced. For the sake of clarity, the algorithm is first explained for steady state (thermal) problems followed by transient (thermal) problems and finally for strongly coupled (thermo-mechanical) problems. In the third section, studies on different test cases are presented. For each test case, cost and parametric analyses of the variational mesh adaption algorithm are carried out and comparison is made with respect to a uniform refinement mesh technique. At the end of this section, first results of a 2-D test case are presented to emphasize the extendability of the presented approach to 2D (and 3D) problems.
## 2 Formulations and methods

In this section, we recall the outlines of the variational formulation of coupled thermo-mechanical problems initially proposed in (Yang, Stainier, and Ortiz 2006) (see also (Stainier 2013)). The presentation is specialized to the 1D case for simplicity, but the framework is actually very general, and is not limited to this case.

### 2.1 Steady state thermal problem

Assuming homogeneous properties, the local uni-dimensional thermal equilibrium reads

$$k \frac{d^2 T}{dx^2} + r = 0 \quad \forall x \in ]0, L[, \qquad (1)$$

where k is the thermal conductivity, r is some external heat source and T is the temperature field. In addition, let us consider the Dirichlet boundary conditions:

$$\begin{aligned} T(x = 0) &= T_0 \\ T(x = L) &= T_1 \end{aligned} \qquad (2)$$

This can be reformulated as a variational problem. The solution of the steady thermal problem is then found by minimizing the following convex potential:

$$\phi(T(x)) = \frac{1}{2} \int_{0}^{L} k \left(\frac{dT}{dx}\right)^2 dx - \int_{0}^{L} rT dx \qquad (3)$$

Indeed, the first variation of equation (3) yields

$$\langle D\Phi, \delta T \rangle = \int_{0}^{L} k \frac{dT}{dx} \frac{d(\delta T)}{dx} dx - \int_{0}^{L} r \delta T dx = 0 \quad \forall \delta T : \delta T(0) = \delta T(L) = 0, \qquad (4)$$

which is the weak form of equation (1). Since Φ is convex, any approximated, discretized temperature field T_h will lead to the following inequality:

$$\Phi(T_h) \geq \Phi(T_{\text{analytical}}) \qquad (5)$$

The variational mesh adaption algorithm presented here exploits directly this property plus the fact that, as a scalar energy, Φ is the sum of all the elementary contributions Φ(T_h) = Σ_e Φ_e(T_h^e) in a finite element discretisation framework. Indeed, the latter property allows us to define local patches of elements, refine them to decrease their values of Φ_e, then add them to those of other patches. One local patch of elements may be isolated of the rest of the mesh by fixing the known temperature field on its boundary during the refinement procedure. Assuming a 1D thermal problem and patch made of one linear element, the mesh refinement procedure consists first in fixing the temperature at its end nodes, and in adding an additional node at the middle of this element. The procedure then amounts to minimize the value of energy-like potential Φ_e over the patch. If the local improvement in the potential is considered as significant, the node is added to the global mesh. Otherwise, the mesh (in the local patch) is kept unchanged. This procedure is shown in figure 1. After looping over all the elements, a new global mesh is defined. However, since the boundary temperature of all the patches has been fixed for the refinement procedure based on the previous known temperature field, a global thermal problem needs to be solved on the refined mesh. This defines an iterative refinement procedure until the required precision is obtained.
<figure><img src="image_2.png" /><figcaption>Figure 1: Mesh adaption algorithm on 1D problem.</figcaption></figure>

Thus, a given local patch may be refined or derefined, depending on the energy difference associated with mesh subdivision, the energy being itself related to the local regularity of the solution. Two tolerance parameters $Tol_r$ and $Tol_d$ associated with the refinement and derefinement bounds respectively are then introduced in the variational mesh adaption algorithm. A third tolerance parameter $Tol_0$ enables to stop the overall iterative procedure of adapted meshes. In one iteration of the procedure, two successive loops on patches are carried out in the sequence for the refinement and derefinement, followed by a new solution on the new mesh.

1. We begin with an arbitrary coarse mesh and solve our problem on that mesh. Get energy-like potential $\phi_{G_1}$.

2. **while** $\frac{\phi_{G_2} - \phi_{G_1}}{\phi_{G_1}} \le Tol_0$ **do**

3. $\phi_{G_1} = \phi_{G_2}$

4. Division of our full geometry $\Omega$ into different patches $\Omega_i$.

5. **for** $\Omega_i = \text{First Patch to} \ \Omega_i = \text{Last Patch do}$

6. Refine the current patch locally and solve a local problem on this small patch with the temperature field on the boundary of the patch imposed (given by the complete solution we calculated in earlier iteration).

7. Calculate the local energy-like potential $\phi_{L_2}$.

8. **if** $\frac{\phi_{L_2} - \phi_{L_1}}{\phi_{L_1}} < Tol_r$ **then**

9. Add patch to list of unrefined patches.

10. **else**

11. Refine the patch in global mesh.

12. Add the patch to the set of refined patches.

13. **end if**

14. **end for**

15. **for** first unrefined patch **to** last unrefined patch **do**

16. Locally derefine the mesh on the patch or locally merge the patch with adjacent unrefined patch.

17. Calculate the values of fields on deleted nodes by interpolation.
18: Calculate local values of energy-like potentials with interpolated values $\phi_{L_2}$ and calculated values $\phi_{L_1}$.

19: if $\frac{\phi_{L_2} - \phi_{L_1}}{\phi_{L_1}} < Told$ then

20: Derefine the patch in the global mesh (or merge two unrefined patches in global mesh).

21: end if

22: end for

23: Solve the problem again on the new mesh thus obtained and calculate $\phi_{G_2}$.

24: end while

## 2.2 Transient thermal problem

Transient thermal problems are governed by the heat equation, written in 1D as

$$c \frac{\partial T}{\partial t} + \frac{\partial}{\partial x} \left( -k \frac{\partial T}{\partial x} \right) = r \quad \forall x \in [0, L] \quad \forall t \in [0, t_{end}], \quad (6)$$

plus a set of initial and boundary conditions. Here, c is the heat capacity per unit volume, k is the thermal conductivity, r is internal heat source density, and T is the temperature field. Assuming that k is homogeneous in $\Omega$, we get

$$c \frac{\partial T}{\partial t} - k \frac{\partial^2 T}{\partial x^2} = r. \quad (7)$$

The above equation can be written in a discretized time setting (backward-Euler) as follows

$$c \frac{\Delta T}{\Delta t} - k \frac{\partial^2 T_{n+1}}{\partial x^2} = r_{n+1}. \quad (8)$$

Let's now define $T_{ref}$ as the reference temperature, and $\theta$ as the variation around it so that $T = T_{ref} + \theta$. Putting this in above equation, one gets

$$\frac{c}{T_{ref}} \frac{\Delta \theta}{\Delta t} - \frac{k}{T_{ref}} \frac{\partial^2 \theta}{\partial x^2} = \frac{r}{T_{ref}}. \quad (9)$$

The Helmholtz free energy can be defined as

$$W(\theta) = -\frac{c}{T_{ref}} \frac{\theta^2}{2}, \quad (10)$$

so that the entropy reads

$$\eta(\theta) = -\frac{dW}{d\theta} = c \frac{\theta}{T_{ref}}. \quad (11)$$

Introducing equations (10) and (11) in equation (9), one gets

$$\frac{\Delta \eta}{\Delta t} - \frac{k}{T_{ref}} \frac{\partial^2 \theta}{\partial x^2} = \frac{r}{T_{ref}}. \quad (12)$$

Let's now define the following incremental energy-like potential:

$$\Phi(\theta_{n+1}) = \int_{0}^{L} (W(\theta_{n+1}) - W(\theta_{n}) + \eta_{n}(\theta_{n+1} - \theta_{n}))dx \\ - \Delta t \int_{0}^{L} \frac{1}{2} \frac{k}{T_{ref}} \left( \frac{\partial \theta_{n+1}}{\partial x} \right)^2 dx + \Delta t \int_{0}^{L} \frac{r}{T_{ref}} (\theta_{n+1} - \theta_{n})dx \quad (13)$$
It is important to note that the first variation of equation (13) gives the weak form of equation (12). The incremental variational principle reads:

$$ \theta_{n+1} = \arg \sup_{\theta} \Phi(\theta_{n+1}) \qquad (14) $$

Although presented here in a heuristic fashion, it can easily be verified that this variational principle is a specific case of the more general variational formulation derived in (Yang, Stainier, and Ortiz 2006).

In the case of transient problems, the mesh is adapted at the first time step as explained in section 2.1. From the second time step onward, the mesh adaption procedure is started with the final adapted mesh obtained at the previous time step as an initial mesh. This avoids a complete re-meshing of the domain at each time step. Two techniques may also be used in order to maintain a low cost of the algorithm without compromising the precision of the solution obtained.

The first one is to avoid to adapt the mesh at every time step. If domains of high gradients of fields of interest do not change much between two time steps, the mesh does not need to be adapted. After the refinement loop over the patches of the procedure, if no major changes have been found in the mesh, adaption at that time step can be bypassed.

The second technique pertains to the values of tolerance parameters in our algorithm. The current solution uses the final adapted mesh from the previous time step as initial mesh. Therefore, well-adapted meshes are required at the early time steps. Accordingly, small values of all of our parameters ($Tol_r, Tol_d, Tol_0$ and $Tol_u$) are hence used at the first time steps. Then, these values can increase with time steps till a given time, after which the values can be set constant.

## 2.3 Thermo-elasticity problem

### 2.3.1 Time discretization

Following (Armero and Simo 1992), the thermoelasticity coupled equations read:

$$ \rho \frac{\partial^2 u}{\partial t^2} = \frac{\partial}{\partial x} (E(\varepsilon - \alpha\theta)) + \rho b \qquad (15) $$

$$ \tilde{c} \frac{\partial \theta}{\partial t} = \frac{\partial}{\partial x} \left( \tilde{k} \frac{\partial \theta}{\partial x} \right) - \alpha E \frac{\partial \varepsilon}{\partial t} + \frac{r}{T_{ref}} \qquad (16) $$

where $\rho$ is the mass density, $u$ denotes the displacement field, $E$ is the Young's modulus, $\varepsilon$ is the linearized strain tensor, $\alpha$ is the coefficient of thermal expansion, $\theta$ is the temperature variation such that $T = T_{ref} + \theta$. External loads $\rho b$ and $r$ represent the body force vector and the external heat source density respectively. One also sets $\tilde{c} = \frac{c}{T_{ref}}$ and $\tilde{k} = \frac{k}{T_{ref}}$, where $c$ is the heat capacity and $k$ is the thermal conductivity. An incremental energy-like potential, convex with respect to the displacement $u$ and concave with respect to the temperature $\theta$, defined at time step $n+1$ is by:

$$ \Phi(u_{n+1}, \theta_{n+1}) = \int_0^L \left\{ \frac{\rho}{2\beta} \left( \frac{\Delta u}{\Delta t} \right)^2 - \frac{\rho}{\beta} \left( v_n + \left( \frac{1}{2} - \beta \right) \Delta t a_n \right) \frac{\Delta u}{\Delta t} \right. \\ + W(\epsilon_{n+1}, \theta_{n+1}) - W(\epsilon_n, \theta_n) + \eta(\theta_{n+1} - \theta_n) - \Delta t \psi(\theta_{n+1}) \\ \left. + \Delta t \frac{r(t_{n+1})}{T_{ref}} (\theta_{n+1} - \theta_n) - b(t_{n+1})(u_{n+1} - u_n) \right\} dx \\ - \left. (\bar{f}(t_{n+1})(u_{n+1} - u_n)) \right|_{x \in \Gamma_t} - \Delta t \left. \left( \bar{q}(t_{n+1}) \frac{\theta_{n+1} - \theta_n}{T_{ref}} \right) \right|_{x \in \Gamma_q} \qquad (17) $$
as proposed in (Yang, Stainier, and Ortiz 2006) and (Radovitzky and Ortiz 1999). The incremental problem then amounts to solving the saddle-point problem

$$u_{n+1}, \theta_{n+1} = \arg\inf_{u} \sup_{\theta} \Phi(u, \theta). \qquad (18)$$

Here, $\beta$ and $\gamma$ are the Newmark time integration parameters, $u$, $v$, $a$ are the displacement, velocity and acceleration fields respectively, $\psi(\nabla\theta)$ is the heat conduction (Biot) potential, $\eta$ the entropy, and $W(\varepsilon, \theta)$ the Helmholtz free energy. The boundary on which external forces $\bar{f}$ are applied is $\Gamma_t$, whereas, $\Gamma_q$ is the boundary on which the heat flux $\bar{q}$ is imposed. The acceleration $a$ and velocity $v$ fields can be obtained from the displacement field $u$ using Newmark's formulas:

$$\begin{aligned} a_{n+1} &= \frac{u_{n+1} - u_n}{\beta \Delta t^2} - \frac{v_n}{\beta \Delta t} - \frac{\frac{1}{2} - \beta}{\beta} a_n \\ v_{n+1} &= v_n + (1 - \gamma) \Delta t a_n + \gamma \Delta t a_{n+1} \end{aligned} \qquad (19)$$

The Newmark's parameters are set at $\beta = \frac{1}{4}$ and $\gamma = \frac{1}{2}$ so that scheme can be unconditionally stable with respect to the time step size $\Delta t$. The Helmholtz free energy can be given as follows:

$$W(\varepsilon, \theta) = \frac{1}{2} E \varepsilon^2 - \alpha E \varepsilon \theta - \frac{1}{2} \frac{c}{T_{ref}} \theta^2 \qquad (20)$$

so that, the entropy $\eta$ reads

$$\eta = -\frac{\partial W}{\partial \theta} = \frac{c\theta}{T_{ref}} + \varepsilon E \alpha. \qquad (21)$$

In addition, the heat conduction potential can be given as

$$\psi(\nabla\theta) = \frac{1}{2} \frac{k}{T_{ref}} \left(\frac{\partial\theta}{\partial x}\right)^2. \qquad (22)$$

### 2.3.2 Algorithm

An adiabatic staggered algorithm is used to solve the problem. This partition of the thermo-mechanical operator is known to preserve numerical stability (Armero and Simo 1992). The mechanical part is solved by

$$u_{n+1} = \arg\inf_{u} \sup_{\theta_{ad}} \Phi_{ad}(u, \theta_{ad}), \qquad (23)$$

where $\Phi_{ad}$ is built from the potential $\Phi$, defined by equation (17), by removing the conduction ($\psi(\nabla\theta)$) and the prescribed heat flux ($\int_{\Gamma_q}$) terms. The stationary condition of $u$ leads to the following system of equations

$$\left[ \begin{matrix} \frac{1}{\beta \Delta t^2} [M] + [E] & [B] \\ [B]^T & -[\tilde{C}] \end{matrix} \right] \left\{ \begin{matrix} \{u_{n+1}\} \\ \{\theta_{ad}\} \end{matrix} \right\} = \left\{ \begin{matrix} \frac{[M]}{\beta \Delta t^2} \{u_n\} + \frac{[M]}{\beta \Delta t} \{\{v_n\} + (\frac{1}{2} - \beta)\Delta t \{a_n\}\} + \{b\} + \{\bar{f}\} \\ -[\tilde{C}] \{\theta_n\} + [B]^T \{u_n\} \end{matrix} \right\} \qquad (24)$$

Here, $[M]$, $[E]$, $[B]$, $[\tilde{C}]$ and $[\tilde{k}]$ are the classical mass, elastic stiffness, coupling, normalized capacity and normalized conduction matrices respectively. Similarly, the thermal part is solved by

$$\theta_{n+1} = \arg\sup_{\theta} \Phi(u_{n+1}, \theta). \qquad (25)$$
Therefore, one gets

$$[-[\tilde{C}] - \Delta t[\tilde{k}]] \{\theta_{n+1}\} = \left\{-[\tilde{C}]\{\theta_n\} + [B]^T\{u_n\} - \frac{\{r\}}{T_{ref}} + \frac{\Delta t}{T_{ref}}\{\bar{q}\} - [B]^T\{u_{n+1}\}\right\} \quad (26)$$

This staggered algorithm consists first in solving the mechanical problem by equation (24) assuming adiabatic conditions. Then, a thermal step at constant geometry is performed solving equation (26). It has been shown (Armero and Simo 1992) that this staggered algorithm is unconditionally stable, provided a fixed time step (Adam 2003). This separation into two steps also allows us to use different meshes for the mechanical and the thermal fields.

Therefore, in one time step, first, the mechanical mesh is adapted according to the adaption procedure explained in section 2.2 to get the mesh that best describes the displacement solution field using a minimum number of elements. During this, thermal fields are interpolated onto the mechanical mesh. Second, the thermal mesh is also adapted to get a mesh that best describes temperature field using a minimum number of elements. During this, the mechanical fields are interpolated on the thermal mesh. The adapted mechanical and thermal meshes at a given time step serve as initial meshes for the mechanical and thermal steps for the next time step.

For unidimensional problems, doing interpolation or projection of fields associated with one mesh onto another mesh won't make big difference in terms of cost and precision. However, doing interpolation rather than projection becomes important for 2D problems (and even more for 3D ones) since the gain of computational cost is significant, while it preserves a good accuracy.

## 3 Mathematical analysis

The variational functional $\Phi$ in equation (17) is equivalent to an $H_1$ norm. This is shown below for transient purely thermal problem for simplicity. Equation (13) can first be rewritten as follows

$$\Phi(\theta_{n+1}) = \int_{0}^{L} \left\{ -\frac{c}{T_{ref}} \frac{\theta_{n+1}^2}{2} + \frac{c\theta_{n+1}\theta_n}{T_{ref}} - \frac{c\theta_n^2}{2T_{ref}} - \Delta t \frac{k}{2T_{ref}} \nabla\theta \cdot \nabla\theta \right\} dx. \quad (27)$$

Defining the following dimensionless field and variable $\bar{\theta} = \frac{\theta}{T_{ref}}$, $\bar{x} = \frac{x\Delta t k}{c}$, and the functional $J = -\frac{\Phi(\bar{\theta})}{cT_{ref}}$, one gets

$$J(\bar{\theta}) = -\frac{\Phi(\bar{\theta})}{cT_{ref}} = \int_{0}^{\bar{L}} \left\{ \frac{1}{2}(\bar{\theta} - \bar{\theta}_n)^2 + \frac{1}{2} \nabla \bar{\theta} \cdot \nabla \bar{\theta} \right\} d\bar{x}. \quad (28)$$

The functional $J(\bar{\theta})$ can be put on the form

$$J(\bar{\theta}) = \frac{1}{2}a(\bar{\theta}, \bar{\theta}) - l(\bar{\theta}) + c_1, \quad (29)$$

where

$$a(\bar{\theta}, \bar{\theta}) = \|\bar{\theta}\|_1^2 = \int_{\bar{\Omega}} (\bar{\theta}^2 + \|\nabla \bar{\theta}\|^2) d\Omega \quad (30)$$

$$l(\bar{\theta}) = \int_{\bar{\Omega}} \bar{\theta} \bar{\theta}_n d\Omega \quad (31)$$

$$c_1 = \int_{\bar{\Omega}} \frac{1}{2} \bar{\theta}_n^2 d\Omega \quad (32)$$
denote the $\mathcal{H}_1$ (squared) norm, a linear form and a constant depending on the solution at the previous time step, respectively. The associated weak form is obtained by minimizing the functional $J(\bar{\theta})$:

$$a(\bar{\theta}, \delta\bar{\theta}) - l(\delta\bar{\theta}) = 0 \quad \forall \delta\bar{\theta} \in \mathcal{V}, \qquad (33)$$

where $\mathcal{V}$ is the space of continuous real valued functions in the domain $\Omega$, satisfying homogeneous Dirichlet boundary conditions. Let $\bar{\theta}_h$ be the finite element solution and $\bar{\theta}$ the exact one. The error in the potential $J(\bar{\theta})$ due to the finite element discretization reads

$$J(\bar{\theta}) - J(\bar{\theta}_h) = \frac{1}{2} (a(\bar{\theta}, \bar{\theta}) - a(\bar{\theta}_h, \bar{\theta}_h)) + l(\bar{\theta}_h) - l(\bar{\theta}). \qquad (34)$$

Provided the $\mathcal{H}_1$ (squared) norm of the error reads

$$\frac{1}{2} \|\bar{\theta} - \bar{\theta}_h\|_1^2 = \frac{1}{2} a(\bar{\theta}, \bar{\theta}) + \frac{1}{2} a(\bar{\theta}_h, \bar{\theta}_h) - a(\bar{\theta}, \bar{\theta}_h), \qquad (35)$$

the combination of equations (35), (34) and (33) yields:

$$\frac{1}{2} \|\bar{\theta} - \bar{\theta}_h\|_1 = J(\bar{\theta}_h) - J(\bar{\theta}). \qquad (36)$$

The difference in the variational potential is directly related to the $\mathcal{H}_1$ (squared) norm of the interpolation error. This result is also found in (Houston and Wihler 2016), although shown differently.

# 4 Results and discussion

## 4.1 Steady state

### 4.1.1 Analytical solution

Consider the following particular boundary conditions and heat source density:

$$\begin{aligned} r &= x^m \\ T_0 &= 0 \\ T_1 &= 0 \end{aligned} \qquad (37)$$

where $m$ is a constant. Therefore, the analytical solution of the problem for a bar of length $L$ reads

$$T = \frac{L^{m+1}x}{(m+1)(m+2)} - \frac{x^{m+2}}{(m+1)(m+2)}. \qquad (38)$$

The energy potential is given by:

$$\Phi(T) = \frac{k}{2} L^{2m+3} \left[ \frac{1}{(m+1)^2(2m+3)} - \frac{1}{(m+1)^2(m+2)^2} \right] - L^{2m+3} \left[ \frac{1}{(m+1)(m+2)^2} - \frac{1}{(m+1)(m+2)(2m+3)} \right]. \qquad (39)$$
### 4.1.2 Numerical solution

The algorithm is started with a coarse initial mesh of two elements. Figure 2 shows the analytical solution as well as the variationally adapted temperature solutions on each iterated mesh, computed with P1-finite elements. These solutions are obtained with $m = 51$, $L = 10$ m, the thermal conductivity and the cross-section area are set to unity. After few iterations and with less than 50 elements, the sharp temperature gradient generated by the large value of $m$ is well captured. Refined elements have been introduced close to that sharp gradient, and few elements are sufficient to represent the remaining part of the solution.

<figure><img src="image_3.png" /><figcaption>Figure 2: Analytical solution and variationally adapted temperature fields.</figcaption></figure>

### 4.1.3 Cost Analysis

In order to assess the usefulness of this algorithm, the error between the computed solution and the analytical one can be plotted as a function of the number of nodes of the mesh. Three cases are considered. In the first case, the plot is made for uniform refinement of the mesh. This can be used as a reference. In the second case, the error at each refinement iteration in the variational adaptive mesh algorithm is plotted with respect to the number of nodes of the mesh. However, since the mesh adaption is done in several iterations, a consistent comparison between a uniform refinement and the variational one should account for path of refinement followed during mesh adaption. One way to accomplish this is to account for a cumulated number of nodes associated with all the calculations performed during the mesh adaption process. Therefore, the error at each refinement iteration is also plotted with respect to the cumulative number of nodes. In the third case, a comparison is performed with the Superconvergent Patch Recovery method (Zienkiewicz and Zhu 1992a; Zienkiewicz and Zhu 1992b) (also denoted ZZ2). The above method provides the global error estimator

$$ \|e_q\|^2 = \sum_{i=1}^{N} \|e_q\|_i^2, \qquad (40) $$

computed from elementary contributions defined in 1D as

$$ \|e_q\|_i^2 = \int_{\Omega_i} \frac{e_q^2}{k} dx, \qquad (41) $$
provided

$$e_q = q^* - q_h, \qquad (42)$$

where $q_h = -k \frac{dT_h}{dx}$ is the finite element reconstruction of the heat flux vector, which is here constant elementwise, and $q^*$ is the nodally reconstructed $ZZ2$ solution from patches (see (Zienkiewicz and Zhu 1992a)). Elements are split into two if the relative value of the local estimator $\|e_q\|_i^2/\|q_h\|_i^2$ is lower than some tolerance, set here identical to that of the variational mesh adaption algorithm.

Figure 3 shows the adapted temperature solutions on each iterated mesh adapted with the $ZZ2$ error estimator as well as the analytical solution. A close path of mesh adaption is observed between the variational (figure 2) and the $ZZ2$ (figure 3) ones, though the final meshes differ.

<figure><img src="image_4.png" /><figcaption>Figure 3: Analytical solution and adapted temperature fields with the $ZZ2$ error estimator.</figcaption></figure>

Figure 4a shows the $L_2$ norm of the error on the temperature field, and figure 4b shows some energy norm of the error. More precisely, the potential energy $\Phi$ (3) is computed with the difference between the approximated and the analytical temperature fields. From figures 4a and 4b, several points can be emphasized. First, the curves of adaptive meshing algorithms are below that linked to uniform mesh refinement, which is quite expectable. Second, the curves of adaptive meshing techniques plotted as a function of the cumulated number of nodes cross from above the uniform mesh one, showing as expected that there is a number of nodes beyond which adaptive remeshing techniques are more performant and more cost effective than a uniform mesh technique. At last, the variational mesh adaption technique appears slightly more performant than the $ZZ2$ one on the range of error computed, both for cumulated and non-cumulated nodes, and both for $L_2$ error and the energy norm.
<figure><img src="image_5.png" /><figcaption>Figure 4: Error norms of the energy-like potential as a function of the number of nodes.</figcaption></figure>

**Remark 1.** In the variational mesh adaption approach, the cost of the refinement step of the mesh essentially consists of the computation of local solutions on refined patches. These problems are very small and cheap ones, because they consist in evaluating the solution at a single node. For thermal analysis, it yields one linear equation per patch to be solved since the temperature is fixed on the patch boundary. The computation of the error indicator for the whole domain thus have a complexity of $O(N)$, where N is the number of nodes. Once the mesh has been adapted, a new solution is computed on the updated mesh, whose complexity is $O(N^3)$ if a direct solver is used. Hence, the cost of the refinement step remains far smaller than the computation of the solution.

**Remark 2.** The nodally reconstructed heat flux vector $q^*$ in the ZZ2 approach requires the solution of a small system of linear equations on each patch, whose dimension equals that of the polynomial basis used for the reconstruction. Though the dimension of the basis is usually small, the solution of a linear system on each patch yields a higher complexity of the refinement step than that achieved by the variational mesh adaption approach.

### 4.1.4 Parametric Analysis

The algorithm explained in section 2 exploits the additive property of the energy-like potential $\Phi$ by summing its local values over all elements. Improving the local value of $\Phi$ on a patch allows to reduce its global value on the mesh, and hence leads to a reduction of error. Three tolerance parameters have been introduced in the algorithm which influence its performance. However, the parameters $Tol_0, Tol_r$ and $Tol_d$ are not independent. In order to study the sole effect of each parameter, others are set to a constant value. The results of this parametric study are shown below.

**Effect of $Tol_0$:** Figures 5 and 6 show the influence of the tolerance parameter $Tol_0$, while fixing $Tol_r$ and $Tol_d$ to a particular value.
<figure><img src="image_6.png" /><figcaption>**Figure 5:** *Effect of the parameter Tol₀, Tol_r and Tol_d being fixed to 10⁻⁴ (figure shows L₂ error vs. the number of nodes).*</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6:** *Effect of the parameter Tol₀, Tol_r and Tol_d being fixed to 0.5 (figure shows L₂ error vs. the number of nodes).*</figcaption></figure>

The parameter *Tol*₀ allows to decide when to stop the algorithm. It doesn't have any effect on the path followed. Therefore, one can observe that, as it decreases, the number of iterations increases. Therefore, the parameter *Tol*₀ should be selected such that the algorithm stops when a solution of a required precision (with respect to the current *Tol*<sub>r</sub>) has been obtained. For example as shown in figure 6, since the *Tol*<sub>r</sub> and *Tol*<sub>d</sub> parameters are set at 0.5, a value of *Tol*₀ ranging between 10³ and 10² is enough.

**Effect of Tol<sub>r</sub>:** Figures 7 and 8 show the influence of the tolerance parameter *Tol*<sub>r</sub>, while fixing *Tol*<sub>0</sub> and *Tol*<sub>d</sub> to a particular value.

<figure><img src="image_8.png" /><figcaption>**Figure 7:** *Effect of the parameter Tol<sub>r</sub>, Tol<sub>0</sub> is set at 10<sup>-2</sup> and Tol<sub>d</sub> at 10<sup>-4</sup> (figure shows L<sub>2</sub> error vs. the number of nodes).*</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8:** *Effect of the parameter Tol<sub>r</sub>, Tol<sub>0</sub> and Tol<sub>d</sub> being fixed to 10<sup>-2</sup> (figure shows L<sub>2</sub> error vs. the number of nodes).*</figcaption></figure>

The parameter *Tol*<sub>r</sub> drives the precision of the converged solution. As shown in figure 7, when the *Tol*<sub>r</sub> ranges between 5 × 10<sup>-3</sup> and 10<sup>-2</sup> the error of the converged solution is of the order of 10<sup>-3</sup>. In the graphs, one can observe that the algorithm carries out a few more iterations after convergence. This occurs because the value of *Tol*<sub>0</sub> has been set to a constant value in order to study the sole effect of *Tol*<sub>r</sub>. Whereas in normal circumstances, the value of *Tol*<sub>0</sub> is changed according to that of *Tol*<sub>r</sub>, so that the algorithm stops immediately after reaching convergence.
**Effect of $Tol_d$:** The results are shown in figures 9 and 10. Convergence and stability of the algorithm depend on the parameter $Tol_d$. The latter should be less than or equal to $Tol_r$, otherwise the algorithm will keep on refining and derefining the same patch entering in an unending loop. For example in figure 10, $Tol_r$ is fixed to $10^2$. All the curves that correspond to the values of $Tol_d$ less than or equal to $10^2$ converge to the solution, whereas all the other curves diverge. This effect can also be observed in one of the curves in figure 8.

### 4.1.5 Improved Algorithm

In problems involving sharp gradients of the main field, many iterations of this iterative adaption process may be performed before convergence occurs, particularly if the initial mesh is coarse. Hence, it could be interesting to accelerate the refinement procedure by dividing an element in more than two elements. An application of this idea is shown in figures 11 and 12. The refinement procedure is so that:

1. **if** $\frac{\phi_{L_2} - \phi_{L_1}}{\phi_{L_1}} > Tol_u$ **then**

2. Subdivide 1 element in 4 elements, where $Tol_u > Tol_r$.

3. **else if** $Tol_u > \frac{\phi_{L_2} - \phi_{L_1}}{\phi_{L_1}} > Tol_r$ **then**

4. Subdivide 1 element in 2 elements.

5. **else if** $\frac{\phi_{L_2} - \phi_{L_1}}{\phi_{L_1}} < Tol_d$ **then**

6. Consider derefinement.

7. **end if**

<figure><img src="image_10.png" /><figcaption>**Figure 9:** Effect of the parameter $Tol_d$, $Tol_0$ is set at $10^{-2}$ and $Tol_r$ at $10^{-3}$ (figure shows $L_2$ error vs. the number of nodes).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10:** Effect of the parameter $Tol_d$, $Tol_0$ is set at $10^{-2}$ and $Tol_r$ at $10^{-2}$ (figure shows $L_2$ error vs. the number of nodes).</figcaption></figure>

**Parametric Analysis:** Figures 13 and 14 show the results of a parametric analysis carried out for the parameter $Tol_u$ while keeping $Tol_r$ constant at $10^4$. All algorithms give equivalent results after convergence and at the beginning. However, there is a big difference in the path followed to reach the converged state.
<figure><img src="image_12.png" /><figcaption>Figure 11: Comparison of the L₂ error of the temperature field between the original and improved algorithms.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Comparison of the energy error between the original and improved algorithms.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Parametric analysis of the parameter TOLu, represented in terms of the L₂ error of the temperature.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14: Parametric analysis of the parameter TOLu, represented in terms of the energy error.</figcaption></figure>

## 4.2 Thermo-elasticity

Consider a bar with homogeneous Dirichlet thermal and mechanical boundary conditions at its two ends:

$$ T(0, t) = T(L, t) = 0; \quad u(0, t) = u(L, t) = 0 \quad \forall t \qquad (43) $$

along with a sinusoidal initial velocity:

$$ u(x, 0) = 0; \quad v(x, 0) = \sin\left(\frac{\pi x}{L}\right); \quad T(x, 0) = 0 \quad \forall x \in [0, L] \qquad (44) $$

This test case has been introduced in (Armero and Simo 1992). With these conditions, the bar is expected to vibrate, though damped through thermal dissipation.

### 4.2.1 Numerical solution fields

An adiabatic staggered scheme is used for the solution, as well as the algorithm of mesh adaption explained in section 2.3.2. The time step is set at 1 second in this test case. The problem is solved
on a very fine mesh (4097 nodes) to obtain a reference solution, which has also been compared with the results obtained in (Armero and Simo 1992) to ensure correctness.

<figure><img src="image_16.png" /><figcaption>(a) Displacement field at time = 1 second.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Temperature field at time = 1 second.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) Displacement field at time = 50 seconds.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(d) Temperature field at time = 50 seconds.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(e) Displacement field at time = 301 seconds.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(f) Temperature field at time = 301 seconds.</figcaption></figure>

<caption>Figure 15: Displacement and temperature fields at different times.</caption>

Figures 15a, 15b, 15c, 15d, 15e and 15f show the displacement and temperature fields at times 1, 50 and 301 seconds respectively. In some parts of the bar, the algorithm has introduced more nodes even where the solution field does not vary much. Two items may explain this behavior. First, our criterion for the mesh refinement is not directly related to the smoothness of the solution profile, but to the value of the energy-like potential. Recall that its value also allows to account for the variation of the solution field with respect to time, hence in a sense it generalizes the quasi-static
approach of the ZZ2 error estimator. Second as explained in section 2.2, the mesh is not adapted at each time step in order to achieve a better cost effectiveness. Therefore, more nodes are required to capture the solution at different time steps for which the same mesh is used. This maintains the accuracy of the solution field and also the cost-effectiveness of the algorithm. However, it is evident from these figures that a very good solution field is captured at all the time steps.

### 4.2.2 Cost analysis

A cost analysis is performed on both thermal and mechanical meshes. The $L_2$ error of the displacement field and of the temperature field are computed on the mechanical and the thermal meshes respectively. The results are shown in figures 16a, 16b, 16c, 16d, 16e, and 16f. It is evident that the introduced mesh adaption algorithm is almost always more cost-effective with respect to a simple uniform mesh for both thermal and mechanical meshes.
<figure><img src="image_22.png" /><figcaption>(a) Cost analysis of the mechanical mesh at time=1 second.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) Cost analysis of the thermal mesh at time=1 second.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(c) Cost analysis of the mechanical mesh at time=50 seconds.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(d) Cost analysis of the thermal mesh at time=50 seconds.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(e) Cost analysis of the mechanical mesh at (f) Cost analysis of the thermal mesh at time=301 seconds.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 16:** Cost analysis of the mechanical and thermal meshes different times</figcaption></figure>

Recall that obtaining cost-effectiveness and a good accuracy of the solution at the first time step of the calculation is crucial, because the following adapted meshes directly depend on the previous ones. For example, at time steps 50 and 301, no mesh adaption is carried out because the meshes used at the previous time steps are good enough to represent the solution.

In this test case, the solution fields do not vary sharply with respect to time and space. Therefore, the algorithm drives the mesh adaption at quite few time steps, from time 1 to 301 seconds. The mechanical mesh is adapted at only 4 time steps, whereas the thermal mesh is adapted only at 3
time steps.

## 4.3 Thermo-elasto-plasticity

### 4.3.1 Numerical solution fields

<figure><img src="image_28.png" /><figcaption>Figure 17: Geometry of the test case (Heuzé et al. 2014).</figcaption></figure>

Heuzé et al. (Heuzé et al. 2014) have extended the well-known viscometer test case to thermo-elastic-plastic solid behaviors in small and large strains. In this test case, the sole mechanical part acts on the thermal part, so that the mechanical problem is solved independently of the thermal problem. The mechanical problem is first solved, followed by the thermal one taking into account the effect of the mechanical solution. The geometry of the problem is shown in figure 17. The gap between the two cylinders is discretized by a radial 1D mesh. Zero displacement is prescribed on the inner cylinder while a driven rotation is prescribed on the outer cylinder. Temperature of external and internal cylinders are fixed to zero. Therefore, the boundary conditions of the problem can be stated as:

$$ 
\begin{align} 
u_{\theta}(r = a) &= 0 \nonumber \\ 
u_{\theta}(r = b) &= u_{\theta}(b) \nonumber \\ 
T(r = a) &= 0 \nonumber \\ 
T(r = b) &= 0 
\end{align} 
\tag{45} $$

where a and b denote the inner and outer radii respectively, and $u_{\theta}$ the curved arc length swept since the finite strain framework is assumed. However the following differences arise between the test case given in (ibid.) and the present one. A hyperelastic-plastic constitutive law is considered here, whereas a hypo-elastic-plastic constitutive law was used in (ibid.). The analytical solution developed in (ibid.) relies on certain assumptions, that is: dilatation effects are neglected, thermal and mechanical parameters are fixed independently of the temperature and additional terms linked to the objective derivative are neglected. The solution developed in small strains is extended to the
large strains in a straightforward manner, but its validity remains bounded from above when the rotations and hence the objective derivative become important. We solve the problem on a very fine mesh (5000 elements) with the numerical data of (Heuzé et al. 2014) and use that solution as a reference one.

<figure><img src="image_29.png" /><figcaption>**Figure 18:** Displacement solution (Heuzé et al. 2014).</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 19:** Equivalent plastic strain reference solution (Heuzé et al. 2014).</figcaption></figure>

The reference solutions in displacement and equivalent plastic strains are plotted in figures 18 and 19 respectively. According to (ibid.), the thermal solution is valid once the viscometer is completely elastic-plastic. As seen in figure 19, at rotation of $\theta = 3°$ of the outer cylinder, the viscometer is completely elastic-plastic. Therefore, a coupled mechanical problem is solved starting at a rotation of outer cylinder of 3 degrees, provided the initial temperature being given by the analytical solution at that rotation, and a rotation evolution of the outer cylinder prescribed so that the plastic crown radius varies exponentially in time (eq.(35) of (ibid.)), consistently with eq.(24) of (ibid.). The reference solution in temperature is shown in figure 20. As seen from figure 18, the displacement field does not vary much but the thermal field presents the interest of a strong temperature gradient close to the inner cylinder. Therefore, it will not be very interesting to use adaptive meshing technique on the displacement mesh. Therefore, we solve our problem by adapting only the thermal mesh and keeping the mechanical mesh constant.

<figure><img src="image_31.png" /><figcaption>**Figure 20:** Temperature reference solution.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Figure 21:** Superposed plot of the analytical (Heuzé et al. 2014) and numerical equivalent plastic strain.</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>(a) Rotation of $\theta = 4^\circ$ of the outer cylinder.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(b) Rotation of $\theta = 5^\circ$ of the outer cylinder.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(c) Rotation of $\theta = 6^\circ$ of the outer cylinder.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(d) Rotation of $\theta = 9^\circ$ of the outer cylinder.</figcaption></figure>

<caption>**Figure 22:** Solution on the adapted mesh for different rotations $\theta$ of the outer cylinder.</caption>

Figures 22a, 22b, 22c and 22d show solutions on adapted mesh at rotations $\theta = 4^\circ, 5^\circ, 6^\circ, 9^\circ$ of the outer cylinder respectively. Figure 21 shows analytical and numerical plastic strain distributions. One can observe that the numerical solution is very close to the analytical solution. The small differences between the numerical solution and the analytical one can be attributed to the different formulations of the mechanical constitutive models in large strains adopted in these two solutions. However, it is harmless for the mesh adaption purpose we are interested in here.

### 4.3.2 Analysis

As seen in figure 20, the domains of interest of the solutions field (domains with high gradients of temperature) do not evolve much in time. Therefore, our algorithm adapts the mesh only at the first time step, and then decides to use the same mesh for the following time steps.
<figure><img src="image_37.png" /><figcaption>Figure 23: $L_2$ error analysis at different rotations $\theta$ of the outer cylinder.</figcaption></figure>

Figures 23a, 23b, 23c, and 23d show the $L_2$ error of the adapted meshes. One can observe that mesh adaption has taken place only at a rotation of 4 degrees, all other time steps use the same mesh. The adaptive meshing still appears more economical than a uniform mesh.

## 4.4 Bidimensional transient thermal test case

In order to demonstrate the extendibility of this algorithm to 2D problems, the results of one test case obtained from our ongoing work in 2D is here presented. Let's consider a rectangular geometry whose boundary temperature is prescribed to zero. An external heat source is introduced in the domain, which follows a circular path in time centered within the rectangle. The heating area at one instant is generated by an arc length of 1 degree and a length of 1m in the radial direction. Figures 24a, 24b, 24c, 24d, 24e and 24f show the solution fields on the adapted meshes at different time steps. A strong mesh adaption is performed in this test case because the location of strong temperature gradients moves with the prescribed heat source. The mesh coarsening upstream from the heat source appears as efficient as the mesh refinement where the heat source is located. It is therefore evident that this adaption strategy works well and is more cost effective than using a simple uniform mesh.
<figure><img src="image_38.png" /><figcaption>**Figure 24:** *Temperature field at different time steps.*</figcaption></figure>

# 5 Conclusion

In this paper, a strategy for mesh adaption based on a variational approach for multiphysics problems has been proposed, in particular attention has been paid to thermo-mechanical problems. The variational approach uses an error indicator to adapt the mesh, based on the optimality property
of the solution. The geometry is divided into patches and according to the level of improvement of
the local value of an energy-like potential, refinement or de-refinement of the patch is performed.

This strategy was first tested on simple steady state and transient thermal problems, for which
a complete parametric analysis was performed. The effect of each parameter was studied and
the strategies of selection of these parameters were discussed. Then the algorithm was applied
to a strongly coupled problem of thermo-elasticity, using an adiabatic staggered algorithm and
two different meshes for the thermal and the mechanical solution fields. Finally, the strategy was
successfully tested on a more complicated thermo-elasto-plasticity benchmark test case. In this
weakly coupled problem, the mechanical solution field was calculated on a fixed mesh, whereas the
thermal mesh was adapted.

In all these test cases, it has been demonstrated that the developed strategy is reliable, eco-
nomical and more effective than using a simple uniform mesh. The first perspective of this method
in current progress is to extend this strategy to 2-D and 3-D problems, using different subdivision
schemes, for example, the mesh subdivision scheme by (Rivara and Inostroza 1997) and (Rivara
1997).

References

Adam, L. (2003). "Modélisation du comportement thermo-élasto-viscoplastique des métaux soumis à grandes déformations. Application au formage superplastique". PhD thesis. Université de Liège.

Ainsworth, M. and JT. Oden (2000). "A Posterior Error Estimation in Finite Element Analysis." In: *Wiley: New York*.

Armero, F. and J.C. Simo (1992). "A new unconditionally stable fractional step method for non-linear coupled thermomechanical problems". In: *International Journal for Numerical Methods in Engineering* 35, pp. 737-766.

Bänsch, E. (1991). "Local mesh refinement in 2 and 3 dimensions". In: *Impact of Computing in Science and Engineering* 3, pp. 181-191.

Barthold, F.-J., M. Schmidt, and E. Stein (1998). "Error indicators and mesh refinements for finite-element-computations of elastoplastic deformations". In: *Computational Mechanics* 22.3, pp. 225-238.

Brancherie, D., P. Villon, and A. Ibrahimbegovic (2008). "On a consistent field transfer in non linear inelastic analysis and ultimate load computation". In: *Computational Mechanics* 42.2, pp. 213-226.

Ciarlet, P. (1988). "Mathematical Elasticity." In: *Three-dimensional Elasticity, North-Holland* 1, pp. 1482-1496.

Dacorogna, B. (1989). "Direct Methods in the Calculus of Variations." In: *Springer: New York*.

Farhat, C., K.C. Park, and Dubois-Pelerin, Y. (1991). "An unconditionally stable staggered algorithm for transient finite element analysis of coupled thermoelastic problems". In: *Computer Method in Applied Mechanics and Engineering* 85.3, pp. 349-365.

Heuzé, T. et al. (2014). "Benchmark tests based on the couette viscometer - II: Thermo-elasto-plastic solid behaviour in small and large strains". In: *Computers and mathematics with applications* 67, pp. 1482-1496.

Houston, P. and T. Wihler (2016). "Adaptive energy minimisation for hp-finite element methods". In: *Computers and Mathematics with Applications* 71, pp. 977-990.

Ladeveze, P., J-P. Pelle, and Ph. Rougeot (1991). "Error estimation and mesh optimization for classical finite elements". In: *Engineering Computations* 8, pp. 69-80.

Mosler, J. and M. Ortiz (2007). "Variational h-adaptation in finite deformation elasticity and plasticity". In: *International Journal of Numerical Methods in Engineering* 72, pp. 505-523.
Mosler, J. and M. Ortiz (2009). "An error-estimate-free and remapping-free variational mesh refinement and coarsening method for dissipative solids at finite strains". In: *International Journal of Numerical Methods in Engineering* 77, pp. 437-450.

Ortiz, M. and EA. Repetto (1999). "Nonconvex energy minimisation and dislocation in ductile single crystals". In: *Journal of the Mechanics and Physics of Solids* 47, pp. 397-462.

Ortiz, M. and L. Stainier (1999). "The variational formulation of viscoplastic constitutive updates". In: *Computational Methods for Applied Mechanical Engineering* 171, pp. 419-444.

Radovitzky, R. and M. Ortiz (1999). "Error estimation and adaptive meshing in strongly nonlinear dynamic problems". In: *Computer Methods in Applied Mechanics and Engineering* 172, pp. 419-444.

Rivara, M. (1991). "Local modification of meshes for adaptive and/or multigrid finite-element methods". In: *Journal of Computational and Applied Mathematics* 36, pp. 79-89.

— (1997). "New longest edge algorithm for the refinement and/or improvement of unstructured triangulations". In: *International Journal for Numerical Methods in Engineering* 40, pp. 3313-3324.

Rivara, M. and P. Inostroza (1997). "Using longest-side bisection techniques for the automatic refinement of delaunay triangulations". In: *International Journal for Numerical Methods in Engineering* 40, pp. 581-597.

Rivara, M. and Ch. Levin (1992). "refinement algorithm suitable for adaptive and multi-grid techniques". In: *Communications in Applied Numerical Methods* 8, pp. 281-290.

Stainier, I. (2013). "A Variational Approach to Modeling Coupled Thermo-Mechanical Nonlinear Dissipative Behaviors". In: ed. by StÃlphane P.A. Bordas. Vol. 46. Advances in Applied Mechanics. Elsevier, pp. 69-126.

Verfürth, R. (1996). "A Review of A Priori Error Estimation and Adaptive Mesh-refinement Techniques". In: *Wiley/Teubner: New York/Stuttgart*.

Yang, Q., L. Stainier, and M. Ortiz (2006). "A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids". In: *Journal of the Mechanics and Physics of Solids* 54, pp. 401-424.

Zienkiewicz, O.C. and J.Z. Zhu (1987). "A simple error estimator and adaptive procedure for practical engineering analysis". In: *International Journal of Numerical Methods in Engineering* 24, pp. 337-357.

Zienkiewicz, Olgierd Cecil and Jian Zhong Zhu (1992a). "The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique". In: *International Journal for Numerical Methods in Engineering* 33.7, pp. 1331-1364.

— (1992b). "The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity". In: *International Journal for Numerical Methods in Engineering* 33.7, pp. 1365-1382.