# Variational h-adaption for problems in multiphysics

Rohit Pethe, Laurent Stainier, Thomas Heuzé

► To cite this version:

Rohit Pethe, Laurent Stainier, Thomas Heuzé. Variational h-adaption for problems in multiphysics. 13e colloque national en calcul des structures, Université Paris-Saclay, May 2017, Giens, Var, France. hal-01899354

HAL Id: hal-01899354

https://hal.science/hal-01899354

Submitted on 19 Oct 2018

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 problems in multiphysics.

R. Pethe<fnref n="1" />, L. Stainier<fnref n="1" />, T. Heuzé<fnref n="1" />

**Résumé** — A 2D variational mesh adaption approach is presented for problems in multiphysics. The mesh is adapted by local subdivision controlled by an energy criterion that serves as an error indicator. Thermal and thermo-mechanical problems are studied here through different test cases to test the algorithm. Then, parametric analyses are carried out, that is the effect of different parameters of our algorithm is studied for mesh adaption.

**Mots clés** — mesh adaption, variational approach, thermomechanics, error indicator, nonlinear coupled problems, multiphysics.

## 1 Introduction

In number of transient problems (purely mechanical and thermal or thermo-mechanical), zones of high gradients of fields of interest evolve with time and loading. Therefore, a dynamic mesh adaption algorithm able to capture the solution in the zones of high gradients is necessary to maintain the required precision. Many methods of mesh adaption are proposed in the literature based on error-estimation. In these methods, the strategy is to adapt the mesh to minimize an error bound among all meshes of fixed size; or by recursive application of local refinement steps [1]. But, certain limitations are experienced by these methods. They work well with linear constitutive models (for example elasticity), but become highly complex when non-linear constitutive models are used. Also, reconstruction of admissible fields is necessary. In addition, for their validity, standard error bounds require the solution to be regular. Therefore, it is difficult and costly to use this approach for complex problems.

However, using variational formulations finite element problems can be expressed as problems of minimization (or maximization) of an energy like potential. This holds true for non-linear structure of the problem as well. However, in some cases, a careful formulation of energy like potential can be needed. For example, in inelastic problems and dynamic problems minimum principle can be obtained by careful time discretization [2, 3]. 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 [4], based on the variational approach of [2]. This technique uses an error indicator rather than an error estimator. It allows mesh adaption in presence of large deformations and non-linear constitutive behavior. In variational approach, an energy like potential is used which indicates error. No error estimates are used at any stage of the algorithm.

In this work, we develop h-adaption algorithm for problems in multi-physics. The variational energy-like potential value is used as an error indicator and the variational principle itself drives mesh refinement and coarsening. This approach can be combined with different geometric refinement procedures. First, the edge bisection technique is used for refinement and coarsening of triangular finite elements. This technique allows to use the same integration points for new refined or coarsened elements and avoid the problem of loss of accuracy by projection of fields on different meshes. Second, the extendibility to general geometric adaptions like RIVARA's LEPP bisection algorithm [5, 6] is demonstrated. Moreover, for problems in multiphysics, each physic has its own temporal and spatial scale. Therefore, different meshes are used to approximate fields associated to different physics and are adapted in the same timestep. This work is the 2-D extension of the technique developed earlier in 1-D [7, 8].

<footnote n="1">¹ GeM, Ecole Centrale Nantes, [pethe.rohit,laurent.stainier,thomas.heuze]@ec-nantes.fr</footnote>
# 2 Algorithm

## 2.1 Simple uniphysics

Assuming homogeneous properties, the local unidimensional thermal equilibrium equation reads :

$$ \mathrm{div}(k \mathbf{grad}(T)) + r = 0 \quad \forall x \in \Omega \quad (1) $$

where k is the thermal conductivity, r is some external heat source and T is the temperature field. Equation (1) should be complemented with appropriate Dirichlet and Neumann boundary conditions. This can be reformulated as a variational problem. The solution for steady thermal problem is then found by minimizing the following convex potential :

$$ \Phi(T(x)) = \int_{\Omega} \frac{1}{2} k \nabla T \cdot \nabla T d\Omega - \int_{\Omega} r T d\Omega - \int_{\Gamma_q} \bar{q} T d\Gamma \quad (2) $$

where Ω is the domain and Γq is the boundary with applied heat flux q̄. Since Φ is convex, any approximated, discretized temperature field Tq will lead to following inequality :

$$ \Phi(T_h) \geq \Phi(T_{\text{analytical}}) \quad (3) $$

The variational mesh adaption algorithm presented here exploits directly this property plus the fact that Φ is the sum of all the elementary contributions Φ(Th) = Σe Φ(The) in a finite element discretization framework. Indeed, the latter property allows us to work on local patches of elements, independently of the remaining of the mesh. 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 thermal problem and a 2D patch made of two triangular linear elements as shown in figure 1, the mesh refinement procedure consists in adding a middle node to the common edge of the two triangles, hence defining four triangles; the potential Φ can be recomputed on that patch provided fixed boundary temperature. This procedure amounts to minimize the value of energy like potential Φ over the patch. If the local improvement in the potential is significant, the node is added to the global mesh. However, since the boundary temperature of all patches has been fixed for the refinement procedure based on the previous known temperature field, a global finite element 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 – Isolation and testing of patches.</figcaption></figure>

One local patch may be refined or derefined, depending on the local regularity of the solution. Thus, two tolerance parameters Tolr and Told associated to the refinement and derefinement bounds respectively are introduced in the variational mesh adaption algorithm. A third tolerance parameter Tol0 enables to stop the overall iterative procedure of adapted meshes. The whole procedure is summarized below : The algorithm starts with an arbitrary coarse mesh and the problem is solved on that mesh. Second, the whole geometry is divided in different patches Ωi such that ∪iΩi = Ω. Third, each patch is refined locally and a local problem is solved. Fourth, comparing the values of energy like potential for the refined patch and unrefined one (from complete solution calculated earlier), refinement is effectively performed (or not performed). A global refined mesh is obtained by doing this procedure on all the patches. Then again, global solution is calculated. This procedure is repeated until convergence.

In case of transient problems, the initial mesh is a general coarse mesh and we adapt the mesh at first time step as explained above. From second time step onwards, 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. Transient thermal algorithms also involve a mesh coarsening step which is in the same spirit as that of refinement step. The expression of incremental energy like potential for transient thermal problems is given by :

$$ \Phi(\theta_{n+1}) = \int_x (W(\theta_{n+1}) + \eta_n \theta_{n+1}) dx - \Delta t \int_x \frac{1}{2} \frac{k}{T_{ref}} \left( \frac{\partial \theta_{n+1}}{\partial x} \right)^2 dx + \Delta t \int_x \frac{r}{T_{ref}} \theta_{n+1} dx \quad (4) $$

where, $\theta = T - T_{ref}$ is the temperature variation with respect to the reference temperature $T_{ref}$, $W(\theta)$ is the Helmholtz free energy, and $\eta$ is the entropy.

## 2.2 Multiphysics

The algorithm for problems in multiphysics uses staggered approach. This approach allows us to use different meshes for different physics. At each resolution step, the corresponding mesh is adapted to represent the field of interest precisely. In the following, this approach is explained for thermo-elasticity problems. Following [9], the thermoelasticity coupled equations read :

$$ \rho \ddot{u} = \nabla (\mathbf{E} : (\boldsymbol{\epsilon} - \alpha \boldsymbol{\theta})) + \rho \mathbf{b} + \mathbf{f} \quad (5) $$

$$ \tilde{c} \dot{\theta} = \operatorname{div}(\tilde{\mathbf{k}} \cdot \nabla \theta) - \alpha : \mathbf{E} : \dot{\boldsymbol{\epsilon}} + \frac{r}{T_{ref}} \quad (6) $$

where $\rho$ is the mass density, $\mathbf{u}$ is displacement field, $\mathbf{E}$ is the elasticity tensor, $\boldsymbol{\epsilon}$ is the strain tensor, $\alpha$ is the coefficient of thermal expansion tensor. External loads $\mathbf{b}$ and $r$ denote body forces and external heat source density respectively. The heat capacity is $c$ and $\tilde{c} = \frac{c}{T_{ref}}$. Similarly, $k$ is the thermal conductivity and $\tilde{k} = \frac{k}{T_{ref}}$. An incremental energy like potential, convex with respect to the displacement $\mathbf{u}$ and concave with respect to the temperature variation $\theta$ defined at time step $n+1$ is given as follows :

$$ \begin{aligned} I(\mathbf{u}_{n+1}, \theta_{n+1}) = & \int_{\Omega} \left\{ \frac{\rho}{2\beta} \left(\frac{\Delta \mathbf{u}}{\Delta t}\right) \cdot \left(\frac{\Delta \mathbf{u}}{\Delta t}\right) - \frac{\rho}{\beta} (\mathbf{v}_n + (\frac{1}{2} - \beta)\Delta t \mathbf{a}_n) \cdot \frac{\Delta \mathbf{u}}{\Delta t} + \right. \\ & \left. W(\epsilon_{n+1}, \theta_{n+1}) - W(\epsilon_n, \theta_n) + \eta(\theta_{n+1} - \theta_n) - \Delta t \psi(\theta_{n+1}) + \right. \\ & \left. \Delta t \frac{r}{T_{ref}} (\theta_{n+1} - \theta_n) - \mathbf{b} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) \right\} d\Omega \\ & - \int_{\Gamma_t} \mathbf{t}_{n+1} \cdot (\mathbf{u}_{n+1} - \mathbf{u}_n) dS - \Delta t \int_{\Gamma_q} q_{n+1} \frac{\theta_{n+1} - \theta_n}{T_{ref}} dS \end{aligned} \quad (7) $$

Here, $\beta$ and $\gamma$ are the Newmark time integration parameters, $\mathbf{u}, \mathbf{v}, \mathbf{a}$ are the displacement, velocity and acceleration fields respectively, $\psi(\theta)$ is the heat conduction (Biot) potential, $\eta$ the entropy, and $W(\boldsymbol{\epsilon}, \boldsymbol{\theta})$ the Helmholtz free energy. An adiabatic staggered approach is used to solve the above problem. The solution of the mechanical part is associated to a saddle point :

$$ u_{n+1} = \inf_{u} \sup_{\theta_{ad}} I_{ad} \quad (8) $$

where $I_{ad}$ is $I$ in equation (7) without conduction ($\psi(\theta)$) and heat flux boundary ($\int_{\Gamma_q}$) terms. Therefore, one gets the following discrete 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\}\} + \{t_{n+1}\} \\ -[\tilde{C}] \{\theta_n\} + [B]^T \{u_n\} \end{matrix} \right\} \quad (9) $$

Similarly, thermal part is solved by

$$ \theta_{n+1} = \sup_{\theta} I \quad (10) $$
that reads through the following discrete system :

$$[-[\tilde{C}] - \Delta t[\tilde{k}]] \{\theta_{n+1}\} = \{-[\tilde{C}]\{\theta_{ad}\} - \frac{\{r\}}{T_{ref}} + \frac{\Delta t}{T_{ref}}\{q_{n+1}\}\} \quad (11)$$

This staggered algorithm consists first in solving the mechanical problem by equation (9) assuming an adiabatic thermal state. Then, a thermal step at constant geometry is performed solving equation (11). It has been shown [9] that this staggered algorithm is unconditionally stable provided a fixed time step. This separation into two steps also allows us to use different meshes for the mechanical and the thermal fields.

Therefore, in one time step, we first solve the mechanical step by interpolating thermal fields onto the mechanical mesh. Second, we adapt the mechanical mesh to get the mesh that best represents the displacement solution field using minimum number of elements. Third, we solve the thermal problem on the thermal mesh by interpolating mechanical fields on the thermal mesh. Fourth, we adapt the thermal mesh to get a mesh that best represents temperature field using minimum number of elements. Then, we start the next time step with the adapted mechanical and thermal meshes at previous time step as initial meshes of mechanical and thermal part respectively.

# 3 Numerical test cases

## 3.1 Steady state thermal problem

Consider a rectangular plate of width W and height H subjected to following boundary conditions :

$$ 
\begin{aligned}
T &= T_1 & \text{at } x = 0, x = W, y = 0 \\
T &= T_2 & \text{at } y = H
\end{aligned}
\qquad (12) $$

Here, the boundary conditions impose a jump in temperature from $T_1$ to $T_2$ at the two corners of the plate. Because of this jump, the energy like potential (2) should indicate drastic improvement in the solution with each refinement on the patches near the two corners associated with temperature jump. Therefore, expected adapted mesh by using the mesh adaption algorithm explained earlier should be very fine near these corners due to this singularity. This simulation of mesh adaption was carried out using a starting mesh of four 6-Node triangular elements as shown in figure 2. The final adapted mesh and the solution field is shown in figure 3. The adapted mesh obtained fits well what was expected. This simulation was carried out using edge bisection technique.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Solution field on initial mesh along with the boundary conditions.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Representation of numerical solution on adapted mesh.</figcaption></figure>

Note that here the obtained mesh is symmetric and well structured. This is because of the structure of this problem and the usage of quadratic elements. But, in general this is not true. One often obtains anisotropic adapted meshes when a simple edge bisection technique is used. The algorithm can be combined with other patching strategies to produce meshes with desired geometric characteristics.

The analytical solution for this test case can be found in [10, Chapter 3]. It is given as follows :

$$ \frac{T - T_1}{T_2 - T_1} = \frac{2}{\pi} \sum_{n=1}^{\infty} \sin\left(\frac{n\pi x}{W}\right) \frac{\sinh(n\pi y/W)}{\sinh(n\pi H/W)} \quad (13) $$
The efficiency of the algorithm is analysed by comparing the error in the numerical solution with respect to the analytical solution (13), as a function of the number of nodes of the mesh.

Three cases are considered. In the first case, plot is made for the uniform meshes. This will be used as a reference. In the second case, the error at each refinement iteration in the 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 to all the calculations performed during the mesh adaption process. Therefore, in the third case, we plot the error at each refinement iteration in the variational adaptive mesh algorithm with respect to the cumulative number of nodes.

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Analysis of the algorithm. Number of nodes on X axis and L₂ error in temperature on Y axis.</figcaption></figure>

Figure 4 shows this plot for L₂ norm of error in temperature field with respect to the analytical solution (13). One can observe that even the curve using cumulated number of nodes crosses the line for uniform meshes. In other words using variational mesh adaption, more precise solution is obtained with less computational cost.

## 3.2 Transient thermal problem

Consider a rectangular plate with homogeneous dirichlet boundary conditions imposed on its boundary $Γ$. Also consider an external heat source that moves along a circular path with time and heats the plate. The center of the circular path is $(x_c, y_c)$, whereas inner and outer radii are $a$ and $b$ respectively. The heat source density is given as :

$$ \begin{aligned} \bar{q} &= Q & \text{if } a \le r \le b \text{ and } \theta < \text{time} < \theta + 1 \\ \bar{q} &= 0 & \text{otherwise.} \end{aligned} \qquad (14) $$

here, $r = \sqrt{(x-x_c)^2 + (y-y_c)^2}$ and $\theta = \tan^{-1}(\frac{y-y_c}{x-x_c})$ is expressed in degrees. This is an example in which the important domains change drastically with time because of movement of external heat source. Therefore, while simulating with adaptive meshing technique, one expects active mesh refinement and coarsening at each time step.

This test case was simulated by using the mesh adaption algorithm explained earlier. At each time step, coarsening precedes the process of refining. The patches are constructed in such a way that we construct meshes with RIVARA's LEPP algorithm [5, 6]. These patches will be called as RIVARA patches in this article.

The results shown in figures 5, 6, 7, 8, 9, 10 demonstrate the efficiency of the process of mesh adaption. The efficiency with respect to a simple uniform mesh is evident in this test case because the mesh in overall domain is uniform except for the important domains of high gradients and external heat source locations.
<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Adapted mesh and solution field at time 24 seconds.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6 – Adapted mesh and solution field at time 72 seconds.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7 – Adapted mesh and solution field at time 120 seconds.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 8 – Adapted mesh and solution field at time 234 seconds.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9 – Adapted mesh and solution field at time 264 seconds.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 10 – Adapted mesh and solution field at time 345 seconds.</figcaption></figure>

### 3.3 Thermo-elasticity problem

Consider a rectangular plate with a circular hole. The plate is submitted to compression through imposed displacement in its Y direction. Making use of the symmetries, only one quarter of the plate is modelled. Zero heat flux boundary condition is prescribed for the thermal part whereas for the mechanical part, negative displacement is imposed on top with symmetry conditions on planes x = 0 and y = 0. Figure 11 shows the computational domain and the boundary conditions.

<figure><img src="image_12.png" /><figcaption>FIGURE 11 – Geometry and boundary conditions for thermo-elasticity problem.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12 – Mechanical adapted mesh obtained using simple edge bisection technique.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 13 – Mechanical adapted mesh obtained using RIVARA’s technique.</figcaption></figure>

Since stress concentration is expected near hole, accordingly the mesh adaption should produce a finer mesh around the hole at the first time step. Unlike the transient thermal test case explained in earlier section, severe mesh coarsening and refinement should not be observed at each timestep because important domains in the field of interest do not change with respect to time.

Figures 12 and 13 show adapted meshes for mechanical part obtained by using simple patch of two traingles and by using RIVARA patch respectively. Observe that in both cases, we have fine mesh near the hole and coarser mesh far from the hole. Therefore, both these meshes represent a good solution field. However, one can observe that in figure 12, mesh has elongated elements (anisotropic mesh). This difference is even more evident in the temperature field shown in figures 14 and 15. Observe the four elongated elements in figure 14. This anisotropy is avoided in RIVARA’s algorithm by finding the LEPP
<figure><img src="image_15.png" /><figcaption>FIGURE 14 – Thermal adapted mesh obtained using simple edge bisection technique.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 15 – Thermal adapted mesh obtained using RIVARA’s technique.</figcaption></figure>

and refining backwards as shown in figure 15. It is important to note that the mesh adaption is performed only at the first time step and no major changes in the mesh are observed after the first time step. The high computational cost associated with using a uniform fine mesh is avoided by using adaptive algorithm irrespective of the nature of patch used.

# 4 Conclusion

A strategy for mesh adaption using a variational approach was extended in this work for 2D multi-physics problems. The variational approach uses an error indicator to adapt the mesh and not an error estimator. The geometry is divided into patches and according to the level of improvement in the local value of energy like potential, decision is made whether to refine or de-refine the patch.

This strategy was first tested on a simple steady state thermal problem using a simple edge bisection technique. Quadratic elements were used in this test case and cost effectiveness with respect to a simple uniform mesh was demonstrated. In the second transient thermal test case, different type of patches are used manifesting the extendibility and compatibility of the algorithm to use different types of geometric refinement techniques. Ability to effectively refine and coarsen the mesh in order to capture precise solution fields maintaining the low computational cost was established. The third test case examined the approach of using different meshes in case of multiphysics. It also compared meshes obtained by using RIVARA’s algorithm and simple patch of two triangles.

# Références

[1] R. Verfürth. A review of a priori error estimation and adaptive mesh-refinement techniques. *Wiley/Teubner : New York/Stuttgart*, 1996.

[2] M. Ortiz and L. Stainier. The variational formulation of viscoplastic constitutive updates. *Computational Methods for Applied Mechanical Engineering*, 171:419–444, 1999.

[3] Q. Yang, L. Stainier, and M. Ortiz. A variational formulation of the coupled thermo-mechanical boundary-value problem for general dissipative solids. *Journal of the Mechanics and Physics of Solids*, 54:401–424, 2006.

[4] J. Mosler and M. Ortiz. Variational h-adaptation in finite deformation elasticity and plasticity. *International Journal of Numerical Methods in Engineering*, 72:505–523, 2007.

[5] M. Rivara and P. Inostroza. Using longest-side bisection techniques for the automatic refinement of delaunay triangulations. *International Journal for Numerical Methods in Engineering*, 40:581–597, 1997.

[6] M. Rivara. New longest edge algorithm for the refinement and/or improvement of unstructured triangulations. *International Journal for Numerical Methods in Engineering*, 40:3313–3324, 1997.

[7] R. Pethe, T. Heuzé, and L. Stainier. Variational approach for thermal mesh adaptation. In *Proceedings of 12ème colloque national en calcul des structures*, CSMA 2015, Giens (Var), France, May 2015.

[8] R. Pethe, T. Heuzé, and L. Stainier. Variational approach for thermo-mechanical mesh adaptation. In *Proceedings of Congrès Français de Méchanique*, CFM 2015, Lyon, France, July 2015.

[9] F. Armero and J.C. Simo. A new unconditionally stable fractional step method for non-linear coupled thermomechanical problems. *International Journal for Numerical Methods in Engineering*, 35:737–766, 1992.

[10] J.P. Holman. *Heat Transfer*. McGraw-Hill Companies, 1996.