# Non-incremental strategies for simulating thermomechanical models with uncertainty

Francisco Chinesta, Etienne Prulière, Amine Ammar, Elías Cueto

► To cite this version:

Francisco Chinesta, Etienne Prulière, Amine Ammar, Elías Cueto. Non-incremental strategies for simulating thermomechanical models with uncertainty. International Journal of Material Forming, 2009, 2 (S1), pp.563-566. 10.1007/s12289-009-0422-z. http-02424650

HAL Id: hal-02424650

https://hal.science/hal-02424650

Submitted on 27 Dec 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
# NON-INCREMENTAL STRATEGIES FOR SIMULATING THERMOMECHANICAL MODELS WITH UNCERTAINTY

F. Chinesta¹, E. Prulière¹*, A. Ammar², E. Cueto³

¹ Ecole Centrale de Nantes, EADS Corporate International Chair, Nantes, France

² Université Joseph Fourier, Laboratoire de Rhéologie, UMR CNRS, Grenoble, France

³ I3A, Universidad de Zaragoza, Spain.

**ABSTRACT:** Description of complex materials involves numerous computational challenges. Thus, the analysis of transient models needs for intensive computation. In some of our former works we proposed a technique based on the separated representation of the unknown field able to circumvent the curse of dimensionality in the treatment of highly multidimensional models. In this work, we are addressing the challenge related to the transient behaviour. For this purpose we propose a separated representation of transient models leading to a non-incremental strategy, allowing to impressive CPU time savings. The use of separated representations allows computing the solution of transient parametric models (the parameters are assumed as additional coordinates). The parameters include thermal coefficients but also thermal sources and/or initial conditions. The resulting curse of dimensionality is circumvented by using the separated representation that implies a complexity that scales linearly with the dimension of the space.

**KEYWORDS:** Reduced representation, Parametric model, Multi-dimensional models, Separated representation.

## 1 INTRODUCTION

The complexity of coupled thermomechanical models defined in domains involving very fine space and time meshes induce computational issues. On the other hand, some times the parametric nature of such models (because of different uncertainties) induces additional difficulties. The use of explicit incremental methods is restricted to very short time steps because of the stability requirements. On the other hand, implicit incremental methods need the solution of a linear system at each time step in the general non-linear case. In any case, the treatment of parametric models requires the solutions of many direct models. To circumvent these difficulties one could use reduced approximation bases within the incremental framework. Another possibility lies in the use of an alternative non-incremental technique. This work explores the use of a particular non-incremental technique based on the separated representation of the unknown field.

## 2 SOLVING PROBLEMS WITH UNCERTAINTY

Let's consider the general heat equation:

$$ \frac{\partial u(t, \mathbf{x})}{\partial t} - \frac{\partial}{\partial \mathbf{x}} \left( k \frac{\partial}{\partial \mathbf{x}} u(t, \mathbf{x}) \right) = f(t, \mathbf{x}) \quad (1) $$

with the following initial condition:

$$ u(t=0, x) = u_0(x) \quad (2) $$

where $u$ is the temperature field and $k$ the thermal conductivity. In general the thermal conductivity depends on the temperature, but in what follows we are assuming that the conductivity is constant but unknown.

### 2.1 ON A SEPARATED REPRESENTATION

To solve multidimensional problems with a low calculation cost, an efficient method has been recently developed. This technique was successfully applied for treating some multi-dimensional models encountered in the kinetic theory description of complex fluids [1] [2] [3] and for treating high resolution thermal homogenization [5]. For the sake of simplicity the method is described here in the 1D transient case defined by:

$$ \frac{\partial u(t, x)}{\partial t} - k \frac{\partial^2 u(t, x)}{\partial x^2} = 0 \quad (3) $$

For that purpose let us consider the weak formulation of the heat equation:

$$ \int_{\Omega} u^* \left( \frac{\partial u}{\partial t} - k \frac{\partial^2 u}{\partial x^2} \right) = 0 \quad (4) $$

<footnote>* Corresponding author: 1 rue de la Noë, BP 92101, 44321 Nantes - France, (33) 2 40 37 16 00, Etienne.Pruliere@ec-nantes.fr</footnote>
$$ \text{where } \Omega = ]0, t_{\max}] \times ]x_{\min}, x_{\max}] = [\Omega_t \times \Omega_x $$

The main idea is the assumption that the solution can be written in the following separated form:

$$ u(t,x) \approx \sum_{i=1}^{N} T_i(t) X_i(x) + u_0 \quad (5) $$

The introduction of $u_0$ in the previous expression enables to enforce an homogeneous initial condition, i.e.:

$$ T_i(t=0) = 0 \quad \forall i \le N \quad (6) $$

Building-up such a solution needs an iterative procedure. At the iteration $n$ we assume that the function $T_i$ and $X_i$ are known $\forall i \le n$. Then we are looking for:

$$ u \approx \sum_{i=1}^{n} T_i(t) X_i(x) + R(t)S(x) + u_0 \quad (7) $$

where $R$ and $S$ have to be determined. We assume the weighting field:

$$ u^* = R^*(t)S(x) + R(t)S^*(x) \quad (8) $$

Thus, the resulting weak formulation writes:

$$ \int_{\Omega} (R^*S + RS^*) \left[ \sum_{i=1}^{n} \left( \frac{dT_i}{dt} X_i - kT_i \frac{d^2 X_i}{dx^2} \right) + \frac{dR}{dt} S - kR \frac{d^2 S}{dx^2} - k \frac{d^2 u_0}{dx^2} \right] d\Omega = 0 \quad (9) $$

that defines a non-linear problem which needs an iterative process to be solved. An alternate directions strategy has given excellent results in our former studies [1],[2],[6]. This method is performed in two steps:

**Step 1:** $R$ being known, we are looking for $S$.

**Step 2:** $S$ being known, we are looking for $R$.

Starting with an arbitrary tentative function $R$, we perform step 1 and then step 2, and again both steps until reaching convergence. Only step 1 will be described now, step 2 being very similar. Thus, $R$ is assumed known and then $R^* = 0$. It results:

$$ \sum_{i=1}^{n} \left( \int_{\Omega_i} R \frac{dT_i}{dt} \int_{\Omega_x} S^* X_i - \int_{\Omega_i} kRT_i \int_{\Omega_x} S^* \frac{d^2 X_i}{dx^2} \right) + \quad (10) $$

$$ \int_{\Omega_i} R \frac{dR}{dt} \int_{\Omega_x} S^* S - \int_{\Omega_i} kR^2 \int_{\Omega_x} S^* \frac{d^2 S}{dx^2} - \int_{\Omega_i} kR \int_{\Omega_x} S^* \frac{d^2 u_0}{dx^2} = 0 $$

The integrals on $\Omega_t$ can be computed numerically. We denote by $\alpha_i$, $\beta_i$ and $\gamma$ these integrals, leading to:

$$ \alpha_0 \int_{\Omega_x} S^* S - \beta_0 \int_{\Omega_x} S^* \frac{d^2 S}{dx^2} = \gamma \int_{\Omega_x} S^* \frac{d^2 u_0}{dx^2} - \sum_{i=1}^{n} \left( \alpha_i \int_{\Omega_x} S^* X_i - \beta_i \int_{\Omega_x} S^* \frac{d^2 X_i}{dx^2} \right) \quad (11) $$

This is a 1D steady-state weak formulation (3D in the general case) which can be solved for instance with the

finite elements method, but if we extract the strong form, then it could be solved using any discretization technique.

## 2.2 ADDRESSING PARAMETRIC MODELS

### 2.2.1 Introducing new coordinates

In this section, we consider a model in which there is an uncertainty on one or many parameters. For illustration purposes, we want to solve the heat equation where the thermal conductivity $k$ is a badly known, but constant, parameter. A way to overcome this difficulty is to take $k$ as a new coordinate belonging to the uncertainty interval. The transient solution for a particular value of the conductivity can then be computed by restricting the general solution to each particular value of this extra-coordinate. Obviously, the price to be paid is the increase of the model dimensionality. However, this is not a serious issue when one proceeds within the separated representation framework just described. The solution is now assumed as:

$$ u(t, x, k) \approx \sum_{i=1}^{N} T_i(t) X_i(x) K_i(k) + u_0 \quad (12) $$

The iterative method described in the previous section is then used to compute the solution.

### 2.2.2 Results

We consider the heat problem in the domain $\Omega_x = ]-1, 1[$ with, for the sake of simplicity, homogeneous essential boundary conditions. The initial condition is: $u_0(x) = 1 - x^2$ and we want to compute the temperature until $t_{max}=0.5s$. and for $1 \le k \le 10$.

The temperature field versus the thermal conductivity is depicted in Figure 1 at time $t=0.1s$. and in Figure 2 at time $t=0.5s$. A reduced number of sums allow computing the solution with a reasonable accuracy. Obviously, the accuracy can be improved by using more sums in the separated approximation. In the limit case, the computed solution approaches the one computed by using the standard finite element method.

<figure><img src="image_2.png" /><figcaption>Figure 1: Temperature versus thermal conductivity and position at time $t=0.1s$.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Temperature versus thermal conductivity and position at time t=0.5s.</figcaption></figure>

# 3 TOWARDS FURTHER REDUCTION

## 3.1 REDUCED SEPARATED REPRESENTATION

Though the separated representation allows considerable time savings, one could envisage additional CPU time savings. For instance the heat equation defined in large 2D or 3D domains needs the solution of large linear systems at step 1 of the iterative strategy presented above. We analyze the solution of the space step:

$$S(x) = \sum_{i=1}^{n_x} B_i(x) s_i \quad (13)$$

where $B_i$ are the basis functions whose number is $n_x$ and $s_i$ are the weights associated with these basis functions. With this approximation, the virtual field yields:

$$S^*(x) = \sum_{i=1}^{n_x} B_i(x) s_i^* \quad (14)$$

Introducing these notations in Eq. (11) it results:

$$\begin{aligned} & \alpha_0 \sum_{i,j=1}^{n_x} s_i^* \left( \int_{\Omega_x} B_i B_j \right) s_j - \beta_0 \sum_{i,j=1}^{n_x} s_i^* \left( \int_{\Omega_x} B_i \frac{d^2 B_j}{dx^2} \right) s_j = \\ & = \gamma \sum_{i=1}^{n_x} s_i^* \left( \int_{\Omega_x} B_i \frac{d^2 u_0}{dx^2} \right) - \\ & \quad - \sum_{i=1}^{n_x} s_i^* \left[ \sum_{j=1}^{n_x} \left( \alpha_j \int_{\Omega_x} B_i X_j - \beta_j \int_{\Omega_x} B_i \frac{d^2 X_j}{dx^2} \right) \right] \end{aligned} \quad (15)$$

that leads to the linear system:

$$Ms = V \quad (16)$$

where $M$ is a $n_x \times n_x$ matrix. If $n_x$ is smaller than the number of nodes distributed in the space, then significant CPU time savings can be accomplished.

**Remark:** Such an approximation may also be performed on the time space for the function $R(t)$.

## 3.2 BUILDING-UP THE APPROXIMATION

The main remaining question is how to build-up an optimal approximation basis. A solution is to build the basis functions by invoking the proper orthogonal decomposition as was performed in [4]. The main issue is the necessity to perform some previous simulations. We have opted for an “a priori” strategy. The idea is to check (at each iteration of the computation of $S$ (respect. $R$)) if the basis needs to be enriched, and if that is the case, to propose an optimal enrichment.

### 3.2.1 Checking the accuracy of the basis

We assume that the basis functions are known for $i \le n_x$. The computation of $S$ requires the resolution of Eq. (15). The weights $s_i$ in that basis are assumed known. To check the accuracy of the basis, we compute the residual function that can be easily derived from Eq. (15):

$$\begin{aligned} \mathrm{Res} = & \alpha_0 \sum_{i=1}^{n_x} B_i s_i - \beta_0 \sum_{i=1}^{n_x} \frac{d^2 B_i}{dx^2} s_i - \\ & \gamma \frac{d^2 u_0}{dx^2} + \sum_{j=1}^{n_x} \left( \alpha_j X_j - \beta_j \frac{d^2 X_j}{dx^2} \right) \end{aligned} \quad (17)$$

If the residual norm is greater than a given small enough value:

$$\|\mathrm{Res}\| = \sqrt{\int_{\Omega_x} |\mathrm{Res}|^2} > \varepsilon \|S\| \quad (18)$$

the basis must be enriched. On the contrary, if the norm of the residual is small enough, we continue to the next iteration keeping the same reduced basis.

### 3.2.2 Enrichment of the basis

The main aim is to define the enriched function $B'$ and the corresponding weight $s'$ in order to cancel the residual. From Eq. (17) we can write:

$$\begin{aligned} & \alpha_0 \left( \sum_{i=1}^{n_x} B_i s_i + B' s' \right) - \beta_0 \left( \sum_{i=1}^{n_x} \frac{d^2 B_i}{dx^2} s_i + \frac{d^2 B'}{dx^2} s' \right) \\ & - \gamma \frac{d^2 u_0}{dx^2} + \sum_{j=1}^{n_x} \left( \alpha_j X_j - \beta_j \frac{d^2 X_j}{dx^2} \right) = 0 \end{aligned} \quad (19)$$

that taking into account Eq. (17) leads to:

$$\alpha_0 B' s' - \beta_0 \frac{d^2 B'}{dx^2} s' = - \mathrm{Res} \quad (20)$$

which is an ordinary differential equation, that could be solved under the normality constraint $\|B'\| = 1$.

Obviously, the solution of that equation implies the solution of a linear system whose size corresponds to the number of nodes used in the discretization, however, this solution is only performed when the reduced basis needs for an enrichment, and not at each iteration as in the standard scheme. As soon as the new function $B'$ is
computed, it is added to the approximation basis, $s'$ is added to the weights list and $n_x$ is increased of one unit. In the worst case, in which an enrichment must be performed at each iteration, the CPU time would be the same than the standard strategy. However, the analyzed cases prove that the CPU time can be reduced in one order of magnitude.

### 3.3 EXAMPLES

Let's consider the same example as in section 2.2.2 but with $t_{max}=1$ and $k$ varying linearly in time such as $k=2t$. Both, the time and the physical spaces are approximated employing 1000 nodes. The basis enrichment criterion in Eq. (18) is set to $\epsilon = 10^{-3}$. The computed solution is depicted in Figure 3.

<figure><img src="image_4.png" /><figcaption>**Figure 3: Temperature versus time and position computed using a reduced separated representation**</figcaption></figure>

A reference solution $u_{ref}$ was computed using an implicit finite differences scheme. The finite difference scheme needed 120 seconds to achieve the results, whereas the standard separated representation required 13 seconds. The reduced separated representation only needed 7 seconds. The error (in the usual Lebesgue norm) is depicted Figure 4. In general the error is lower than 1% except in the boundaries neighbourhood where the solution being close to 0, the error in percent is artificially higher.

<figure><img src="image_5.png" /><figcaption>**Figure 4: Error of the reduced separated representation method versus time and position**</figcaption></figure>

In this case the CPU time saving is not significant because the space domain is 1D, but the higher is the dimension of the space, the more significant will be the expected CPU time savings when one proceeds within the reduced separated representation framework.

## 4 CONCLUSIONS

The separated representation has been successfully applied as a general solver for parametric linear and non-linear transient models. A further reduction was accomplished by using reduced bases for accounting each model coordinate. This technique opens new possibilities to perform efficient simulations of complex thermomechanical models.

## REFERENCES

[1] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings. A New Family of Solvers for Some Classes of Multidimensional Partial Differential Equations Encountered in Kinetic Theory Modeling of Complex Fluids. J. Non-Newtonian Fluid Mech. 139:153-176, 2006.

[2] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings. A New Family of Solvers for Some Classes of Multidimensional Partial Differential Equations Encountered in Kinetic Theory Modeling of Complex Fluids. Part II: Transient Simulation Using Space-Time Separated Representation. J. Non-Newtonian Fluid Mech. 144:98-121, 2007.

[3] B. Mokdad, E. Prulière, A. Ammar, F. Chinesta. On the Simulation of Kinetic Theory Models of Complex Fluids Using the Fokker-Planck Approach. Applied Rheology.

[4] D. Ryckelynck, F. Chinesta, E. Cueto, A. Ammar. On the A Priori Model Reduction: Overview and Recent Developments. Archives of Computational Methods in Engineering. 13(1):91-128, 2006.

[5] F. Chinesta, A. Ammar, F. Lemarchand, P. Beauchene, F. Boust, Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization, Comput. Methods Appl. Mech. Engrg. 197:400-413, 2008.

[6] E. Prulière, A. Ammar, N. El Kissi, F. Chinesta. Recirculating flows involving short fiber suspensions: numerical difficulties and efficient micro-macro solvers. Archives of Computational Methods in Engineering, to be published.