Variational h-adaption for coupled thermomechanical problems


Purpose
The purpose of this paper is to present a variational mesh h-adaption approach for strongly coupled thermomechanical problems.


Design/methodology/approach
The mesh is adapted by local subdivision controlled by an energy criterion. Thermal and thermomechanical problems are of interest here. In particular, steady and transient purely thermal problems, transient strongly coupled thermoelasticity and thermoplasticity problems are investigated.


Findings
Different test cases are performed to test the robustness of the algorithm for the problems listed above. It is found that a better cost-effectiveness can be obtained with that approach compared to a uniform refining procedure. Because the algorithm is based on a set of tolerance parameters, parametric analyses and a study of their respective influence on the mesh adaption are carried out. This detailed analysis is performed on unidimensional problems, and a final example is provided in two dimensions.


Originality/value
This work presents an original approach for independent h-adaption of a mechanical and a thermal mesh in strongly coupled problems, based on an incremental variational formulation. The approach does not rely on (or attempt to provide) error estimation in the classical sense. It could merely be considered to provide an error indicator. Instead, it provides a practical methodology to adapt the mesh on the basis of the variational structure of the underlying mathematical problem.



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 energylike 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.

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.

Steady state thermal problem
Assuming homogeneous properties, the local uni-dimensional thermal equilibrium reads 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: This can be reformulated as a variational problem. The solution of the steady thermal problem is then found by minimizing the following convex potential: Indeed, the first variation of equation (3) yields which is the weak form of equation (1). Since Φ is convex, any approximated, discretized temperature field T h will lead to the following inequality: 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 e h ) 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. 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 T ol r and T ol d associated with the refinement and derefinement bounds respectively are then introduced in the variational mesh adaption algorithm. A third tolerance parameter T ol 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 energylike potential φ G1 .

4:
Division of our full geometry Ω into different patches Ω i .

5:
for Ω i = F irst P atch to Ω i = Last P atch 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).
Add patch to list of unrefined patches. Refine the patch in global mesh.

12:
Add the patch to the set of refined patches. 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.
Derefine the patch in the global mesh (or merge two unrefined patches in global mesh).
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 Ω, we get The above equation can be written in a discretized time setting (backward-Euler) as follows Let's now define T ref as the reference temperature, and θ as the variation around it so that T = T ref + θ. Putting this in above equation, one gets The Helmholtz free energy can be defined as so that the entropy reads Introducing equations (10) and (11) in equation (9), one gets Let's now define the following incremental energy-like potential: It is important to note that the first variation of equation (13) gives the weak form of equation (12). The incremental variational principle reads: 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 (T ol r ,T ol d ,T ol 0 and T ol 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.

Time discretization
Following (Armero and Simo 1992), the thermoelasticity coupled equations read: where ρ is the mass density, u denotes the displacement field, E is the Young's modulus, ε is the linearized strain tensor, α is the coefficient of thermal expansion, θ is the temperature variation such that T = T ref + θ. External loads ρb and r represent the body force vector and the external heat source density respectively. One also setsc = c 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 θ, defined at time step n + 1 is by: 6 as proposed in (Yang, Stainier, and Ortiz 2006) and (Radovitzky and Ortiz 1999). The incremental problem then amounts to solving the saddle-point problem Here, β and γ are the Newmark time integration parameters, u, v, a are the displacement, velocity and acceleration fields respectively, ψ(∇θ) is the heat conduction (Biot) potential, η the entropy, and W (ε, θ) the Helmholtz free energy. The boundary on which external forcesf are applied is Γ t , whereas, Γ q is the boundary on which the heat fluxq is imposed. The acceleration a and velocity v fields can be obtained from the displacement field u using Newmark's formulas: The Newmark's parameters are set at β = 1 4 and γ = 1 2 so that scheme can be unconditionally stable with respect to the time step size ∆t. The Helmholtz free energy can be given as follows: so that, the entropy η reads In addition, the heat conduction potential can be given as (22)

Algorithm
An adiabatic staggered algorithm is used to solve the problem. This partition of the thermomechanical operator is known to preserve numerical stability (Armero and Simo 1992).The mechanical part is solved by where Φ ad is built from the potential Φ, defined by equation (17), by removing the conduction (ψ(∇θ)) and the prescribed heat flux ( Γq ) terms. The stationary condition of u leads to the following system of equations Therefore, one gets 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 unidimentional 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.

Mathematical analysis
The variational functional Φ 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 Defining the following dimensionless field and variableθ = θ T ref ,x = x∆tk c , and the functional The functional J(θ) can be put on the form where denote the 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 where V is the space of continuous real valued functions in the domain Ω, satisfying homogeneous Dirichlet boundary conditions. Letθ h be the finite element solution andθ the exact one. The error in the potential J(θ) due to the finite element discretization reads Provided the H 1 (squared) norm of the error reads the combination of equations (35), (34) and (33) yields: The difference in the variational potential is directly related to the H 1 (squared) norm of the interpolation error. This result is also found in (Houston and Wihler 2016), although shown differently.

Results and discussion
4.1 Steady state

Analytical solution
Consider the following particular boundary conditions and heat source density: where m is a constant. Therefore, the analytical solution of the problem for a bar of length L reads .
The energy potential is given by: 9

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.

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 computed from elementary contributions defined in 1D as where q h = −k 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 2 i / q h 2 i 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 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 Φ (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.  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.

Parametric Analysis
The algorithm explained in section 2 exploits the additive property of the energy-like potential Φ by summing its local values over all elements. Improving the local value of Φ 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 T ol 0 ,T ol r and T ol 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 T ol 0 : Figures 5 and 6 show the influence of the tolerance parameter T ol 0 , while fixing T ol r and T ol d to a particular value.  The parameter T ol 0 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 T ol 0 should be selected such that the algorithm stops when a solution of a required precision (with respect to the current T ol r ) has been obtained. For example as shown in figure 6, since the T ol r and T ol d parameters are set at 0.5, a value of T ol 0 ranging between 10 3 and 10 2 is enough.
Effect of T ol r : Figures 7 and 8 show the influence of the tolerance parameter T ol r , while fixing T ol 0 and T ol d to a particular value.  The parameter T ol r drives the precision of the converged solution. As shown in figure 7, when the T ol r ranges between 5 × 10 −3 and 10 −2 the error of the converged solution is of the order of 10 −3 . In the graphs, one can observe that the algorithm carries out a few more iterations after convergence. This occurs because the value of T ol 0 has been set to a constant value in order to study the sole effect of T ol r . Whereas in normal circumstances, the value of T ol 0 is changed according to that of T ol r , so that the algorithm stops immediately after reaching convergence.
Effect of T ol d : The results are shown in figures 9 and 10. Convergence and stability of the algorithm depend on the parameter T ol d . The latter should be less than or equal to T ol r , otherwise the algorithm will keep on refining and derefining the same patch entering in an unending loop. For example in figure 10, T ol r is fixed to 10 2 . All the curves that correspond to the values of T ol 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.

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: Subdivide 1 element in 2 elements. 5: else if Consider derefinement. 7: end if  Parametric Analysis: Figures 13 and 14 show the results of a parametric analysis carried out for the parameter T ol u while keeping T ol 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.

Thermo-elasticity
Consider a bar with homogeneous Dirichlet thermal and mechanical boundary conditions at its two ends: along with a sinusoidal initial velocity: 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.

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.   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 instroduced 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.

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.   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. Heuzé et al. (Heuzé et al. 2014) have extended the well-known viscometer test case to thermoelastic-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:

Thermo-elasto-plasticity
where a and b denote the inner and outer radii respectively, and u θ 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 19: Equivalent plastic strain reference solution (Heuzé et al. 2014).
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 θ = 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 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.

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.     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.

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.

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, economical 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).