# A Priori Model Reduction Method for Thermo-mechanical Simulations

David Ryckelynck

► To cite this version:

David Ryckelynck. A Priori Model Reduction Method for Thermo-mechanical Simulations. International Journal of Forming Processes, 2002, 5 (2-3-4), pp.477-491. 10.3166/ijfp.5.477-491. hal-03760221

HAL Id: hal-03760221

https://hal.science/hal-03760221v1

Submitted on 2 Mar 2023

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 CC BY 4.0 - Attribution - International License
A Priori Model Reduction Method
for Thermo-mechanical Simulations

David Ryckelynck

Laboratoire de mécanique des systèmes et des procédes

UMR CNRS - ENSAM - ESEM

Ecole Nationale Supérieure d'Arts et Métiers

151 boulevard de l'hôpital

75013 Paris

France

david.ryckelynck@paris.ensam.fr

ABSTRACT. We developed an algorithm to extract reduced-order model from Finite Element model, in the framework of thermo-mechanical simulations. This is an a priori approach in so far as we do not use the state forecasted by the Finite Element model. The basis functions of the reduced-order model are defined for large evolutions of the state variables of the studied system. An iterative non-incremental algorithm, a LATIN algorithm, constructs them. As the basis functions are improved during the iterative procedure, an approximate state evolution is also bettered. Thanks to the reduced-order model, an accurate overview of the variables evolutions is quickly obtained over the entire time interval. This way we hope to forecast process defaults as quickly as they are significant. An important advantage of this model reduction method is that we do not need a special reflection of what could be the main phenomena that happen during the case studied.

KEYWORDS: model reduction, Karhunen-Loève expansion, Krylov subspace, non-incremental approach, contact modelling.
# 1. Introduction

Thanks to the Finite Element Method, we are able to describe accurately the thermo-mechanical transformations that occur during processes like casting, forming and cutting. It is quite easy to construct a very complex finite element model. It is more difficult to forecast the state of the studied system during the process, by solving a time-dependant and non-linear thermo-mechanical problem. Generally a good knowledge of what are the main significant phenomena allows to simplify the finite element model. Obviously, the simpler the finite element model is, the faster we forecast the state of the system.

An interesting approach to simplify a finite element model is to build a reduced-order model, thanks to an algorithm. Finite element modelling uses generally a large number of variables to describe the state of the studied system. The aim of model reduction is to define a transformation matrix to get few basis functions, from the Finite Element shape functions, to describe the spatially distributed state. In the framework of non-linear time-dependant problems we can distinguish an a posteriori approach based on Karhunen-Loève expansion [SIR 91] [PAR 01] and an a priori approach based on the Krylov subspace [KNO 95] [WEI 95] [SUN 01]. The first kind of approach is based on the knowledge of the evolution of the state variables of the system, the second one is not based on this knowledge. In the Karhunen-Loève method, state snapshots, taken at various time instants, span a small subspace. The principal directions of this subspace are chosen as basis functions. Hence, the Karhunen-Loève expansion allows to extract the most significant part of the evolution of the state snapshots. The Krylov subspace allows to forecast the state evolution of the system. Using this subspace, only few unknowns are needed during the incremental computation of the state. At a given time step, the basis functions are built with the residual of the non-linear equations of the finite element problem and a matrix defined to correct the approximate state forecasted at the end of the time step. The basis functions can be used to simplify the computation of the state during several time increments [KNO 95]. In practice, the same Krylov subspace is used only for few time increments. This subspace had to be adapted if there is large evolution of the state variables after the time step when it was constructed. In this paper we propose an a priori reduction method because we want to use it as a computational strategy.

The purpose of our approach is to define residuals over a time interval involving a large state evolution. To do so, we need to improve the computation of an approximate state over the same time interval, thanks to an iterative algorithm. Obviously this approximate state is described thanks to the basis functions of the reduced-order model. A simple reduced-order model is supposed to be known and several correction stages are used to improve the basis functions thanks to the residuals and to improve also the approximate state over the entire time interval. So, this computational strategy is necessarily a non-incremental one. In order to improve the approximate state at each stage by a linear problem, we use the LATIN method [LAD 85] [LAD 96] to define our algorithm. Because the correction stage is a linear stage, we need basis functions for both displacements and stresses. The system Σ is split into sub-domains. Basis
functions are defined for each sub-domain $\Omega$ of $\Sigma$. Moreover, in the framework of the LATIN method, the basis functions for stresses satisfy some equilibrium conditions. Therefore we need in fact couples of basis functions to describe stresses inside each sub-domain $\Omega$ and normal stresses over the boundary of $\Omega$. Even if the approximate solution is far from the solution of the finite element problem, the residuals give efficient indications to construct the reduced-order model. Maybe it is due to equilibrium conditions satisfied when the corrections are computed. Several correction stages of the approximate state must be considered in order to obtain residuals that span a sufficient large subspace to define the reduced-order model. The subspace spanned by the residuals is called the LATIN subspace. During the correction stages the basis functions of the reduced-order model are extracted from the LATIN subspace thanks to the Karhunen-Loève expansion, by selecting the most significant patterns of the known approximate state. The subspace associated to these basis functions is called the LATIN-KL subspace. A 2D example illustrates the capability of the model reduction technique. It corresponds to a kind of casting problem. The part material is a kind of aluminium and the mould material is a kind of sand. We obtain a reduced-order model defined for large non-linear evolutions, during all the cooling of the part. Non-linear couplings are due to the thermal contraction of the part and thermal expansion of the mould, which modify the thermo-mechanical contact conditions between the part and the mould. An accurate approximate state evolution is obtained with few basis functions obtained during 100 correction stages.

## 2. The Finite Element model

In order to define the LATIN algorithm, we introduce here a convenient manner to describe the finite element model. The considered thermo-mechanical problem involves small strains, small displacements, elasticity, and contact producing thermo-mechanical coupling. Equations are given for each sub-domain $\Omega$. Boundary conditions are given for interfaces between a sub-domain $\Omega$ and a sub-domain $\Omega'$, and for interfaces between a sub-domain $\Omega$ and the boundary $\partial\Sigma$ of the system. A mesh is used to describe the displacement field $\underline{U}$ and the temperature field $T$ by introducing matrix of shape functions $\mathbf{N}_U$ and $\mathbf{N}_T$, and Finite Element degrees of freedom $q$ and $\theta$ respectively. The shape functions $\mathbf{N}_U$ and $\mathbf{N}_T$ define two subspaces $\mathcal{U}$ and $\bar{\mathcal{T}}$ respectively. We assume that these subspaces are large enough to represent accurately the evolution of the gradients of $\underline{U}$ and $T$. The time interval $[0, t_f]$ is divided into time steps $[t_j, t_{j+1}]$. The forward Euler scheme is used to solve the time differential equations.

### 2.1. Equations defined inside $\Omega \in \Sigma$

Inside each sub-domain $\Omega$ the following conditions must be satisfied. There is initial conditions for temperatures and stresses $\sigma$. Stresses, normal stresses $\underline{F}$ over $\partial\Omega$, temperatures, heat transfer $\underline{g}$, internal energy $e$ normal heat transfer $\phi$ over $\partial\Omega$
and heat dissipation $r$, must satisfy equilibrium conditions at any time instant $t_{j+1}$ (forward Euler scheme):

$$ \int_{\Omega} \mathrm{Tr} [\varepsilon (\underline{U}^*) \sigma] d\Omega = \int_{\partial\Omega} \underline{U}^* \cdot \underline{F} dS \quad \forall \underline{U}^* \in \mathcal{U} \quad [1] $$

$$ \begin{aligned} & \int_{\Omega} T^* \dot{e} d\Omega - \int_{\Omega} \underline{\mathrm{grad}}(T^*) \underline{g} d\Omega = \\ & - \int_{\partial\Omega} T^* \phi dS + \int_{\Omega} T^* r d\Omega \quad \forall T^* \in \mathcal{T} \end{aligned} \quad [2] $$

The thermo-mechanical constitutive relations concern the elastic strain rate $\dot{\varepsilon}^e$, the thermal expansion $\dot{\varepsilon}^{th}$, the heat dissipation, the internal energy and the heat transfer. The heat dissipation is a fraction of the irreversible mechanical rate of work.

$$ \varepsilon (\dot{\underline{U}}) = \dot{\varepsilon}^e + \dot{\varepsilon}^{th} \quad [3] $$

$$ \dot{\sigma} = \mathbf{K}(T) \dot{\varepsilon}^e \qquad \dot{\varepsilon}^{th} = \alpha(T) \dot{T} \quad [4] $$

$$ \dot{e} = \rho c(T) \dot{T} \qquad \underline{g} = -k(T) \underline{\mathrm{grad}}(T) \quad [5] $$

## 2.2. Boundary conditions over $\partial\Omega$

On each integration point on the boundary of $\Omega$ there is an interface between an other sub-domain $\Omega'$ or the considered point belongs to the boundary of $\Sigma$. Two kinds of interfaces between $\Omega$ and $\Omega'$ are considered : a perfect interface and a contact interface. In case of perfect interface the displacement and the temperature must be continuous and some equilibrium conditions must be satisfied. We assume that heat conduction decreases with the gap $j$ such that: $h(j) = h_o e^{-\gamma j}$.

# 3. The non-incremental approach to construct the reduced-order model

The purpose of the LATIN method is to build an algorithm to compute corrections of an approximate state thanks to stages, both global and linear, defined over the entire time interval. The residuals we want to define, must characterise what could be the mean value over $[0, t_f]$ of the corrections by using the Finite Element subspaces $\mathcal{U}$ and $\mathcal{T}$. But, to avoid expensive computations we do not want to compute these corrections by using the Finite Element subspaces. In order to define the residuals, we present the LATIN algorithm as if we used the Finite Element subspaces. Hence we modify the formulation of the linear correction stage to employ the reduced-order model for the computation of the corrections.
## 3.1. A LATIN algorithm

The Finite Element problem presented in section 2 is both global and non-linear. But, we chose a formulation such that the global equations are linear and the non-linear equations are local equations. The LATIN method is based on three principles : P1, P2 and P3 [LAD 85] [LAD 96]:

- P1 : in order to split the difficulties (global and non-linear equations), two groups of equations are created from the reference problem : a group of local equations and a group of linear equations;

- P2 : the algorithm is an iterative procedure that provides a solution for each group of equations at each iteration, these solutions are defined over the structure and over the entire time interval;

- P3 : an appropriate space-time representation has to be used to solve the global equations.

The reference problem is the Finite Element problem. The linear equations are: the initial conditions, the equilibrium conditions (1) (2) and the strain definition (3). Necessarily, the other equations belong to the group of the local equations.

The third principle of the LATIN method allows to exploit various kind of time an space representation of the variables ($U, \sigma, F, T, \phi$). In our approach we use basis functions associated to the reduced-order model. Different approaches were proposed to define a convenient representation of the corrections ($U - U_n, \sigma - \sigma_n, ...$), expressed linearly in terms of several space functions [LAD 96] [ART 92]. P. Ladevèze has proposed in [LAD 96] an approach very closed to the Karhunen-Loève expansion to construct these space functions at each correction stages. To employ the previously computed space functions during a given correction stage, P. Bussy [BUS 90] proposed to compute a part of the corrections while using the known space functions as basis functions. In practice, this last approach is efficient enough to provide the entire correction of the state variables [PEL 00]. The proposed approach can be understood as an extension of the approach proposed by P. Bussy, in order to construct a reduced-order model of a process, and not only to simplify the computation of the correction stage of a LATIN algorithm. The most simple approximate state ($U_o, \sigma_o, F_o, T_o, \phi_o$) is obtained from the initial state ($U_{ini}, \sigma_{ini}, F_{ini}, T_{ini}, \phi_{ini}$).

## 3.2. The Finite Element formulation of the correction stage

In accordance with the principle P2, the purpose of the correction stage is to construct a solution of the linear equations (1), (2) and (3). To do so, we have to replace the local equations by linear equations defined thanks to search directions. Let's consider that an approximate state ($U_n, \sigma_n, F_n, T_n, \phi_n$) satisfying (1), (2) and (3), is known. A fixed-point approach allows to replace the local equations defined over $\Omega$ by:
$$ \dot{\boldsymbol{\sigma}} = \mathbf{K}(T_n) \dot{\boldsymbol{\epsilon}}^e \qquad \dot{\boldsymbol{\epsilon}}^{th} = \boldsymbol{\alpha}(T_n) \dot{T}_n \qquad [6] $$

$$ r = \mathrm{Tr} \left[ \left( \boldsymbol{\varepsilon} \left( \dot{\underline{U}}_n \right) - \dot{\boldsymbol{\epsilon}}_n^e \right) \boldsymbol{\sigma}_n \right] \qquad [7] $$

$$ \dot{e} = \rho c(T_n) \dot{T} \qquad \underline{g} = -k(T_n) \underline{\mathrm{grad}}(T) \qquad [8] $$

On the boundary of each sub-domain, search directions have been defined in [CHA 96] for mechanical problems involving contact conditions. Similar search directions are used at the correction stage, for interfaces between $\Omega$ and $\Omega'$:

$$ \dot{\underline{F}} = \dot{\underline{F}}_{n+1/2} - \hat{k}^{me} \left( \dot{\underline{U}} - \dot{\underline{U}}_{n+1/2} \right) \qquad [9] $$

$$ \phi = \phi_{n+1/2} + \hat{h}^{th} (T - T_{n+1/2}) \qquad [10] $$

where $\underline{F}_{n+1/2}$, $\underline{U}_{n+1/2}$, $\underline{F}'_{n+1/2}$, $\underline{U}'_{n+1/2}$, $T_{n+1/2}$, $T'_{n+1/2}$, $\phi_{n+1/2}$ and $\phi'_{n+1/2}$ are defined by the contact conditions or the perfect interface conditions and these equations:

$$ \dot{\underline{F}}_{n+1/2} = \dot{\underline{F}}_n + \hat{k}^{me} \left( \dot{\underline{U}}_{n+1/2} - \dot{\underline{U}}_n \right) \\ \dot{\underline{F}}'_{n+1/2} = \dot{\underline{F}}'_n + \hat{k}^{me} \left( \dot{\underline{U}}'_{n+1/2} - \dot{\underline{U}}'_n \right) \qquad [11] $$

$$ \phi_{n+1/2} = \phi_n - \hat{h}^{th} (T_{n+1/2} - T_n) \\ \phi'_{n+1/2} = \phi'_n - \hat{h}^{th} (T'_{n+1/2} - T'_n) \qquad [12] $$

To find the variables $\underline{F}_{n+1/2}$, $\underline{U}_{n+1/2}$, $\underline{F}'_{n+1/2}$, $\underline{U}'_{n+1/2}$, $T_{n+1/2}$, $T'_{n+1/2}$, $\phi_{n+1/2}$ and $\phi'_{n+1/2}$, in case of contact condition, on an integration point, we assume that there is no contact, if it's not correct hence there is contact and we assume there is adherence. At last, if it's not correct there is sliding contact. In case of an integration point on $\partial\Sigma$, the variables $\underline{F}_{n+1/2}$, $\underline{U}_{n+1/2}$, $T_{n+1/2}$ and $\phi_{n+1/2}$ must satisfy the boundary conditions. In practice, the Young modulus $E$ and a characteristic length $d$ of $\Omega$ are used to define the parameters of the algorithm such that:

$$ \hat{k}^{me} = \frac{E(T_n)}{d} \qquad \hat{h}^{th} = \frac{k(T_n)}{d} $$
### 3.3. The LATIN subspace

Basis functions have to be defined for $\underline{U}$, $\underline{\sigma}$, $\underline{E}$, $T$ and $\phi$. Let's note respectively $\underline{W}_U$, $\underline{W}_\sigma$, $\underline{W}_F$, $W_T$ and $W_\phi$ the residuals defining these basis functions. The three residuals $\underline{W}_U$, $\underline{W}_\sigma$ and $\underline{W}_F$ must characterise respectively $\underline{U} - \underline{U}_n$, $\underline{\sigma} - \underline{\sigma}_n$ and $\underline{E} - \underline{E}_n$, over $[0, t_f]$. To avoid several resolutions of Finite Element problems during a correction stage, we introduce two given scalar time functions $\lambda^{me}$ and $\lambda^{th}$ such that:

$$ 
\begin{aligned}
\underline{U} &= \underline{U}_n + \lambda^{me}(t) \underline{W}_U \\
\underline{\sigma} &= \underline{\sigma}_n + \lambda^{me}(t) \underline{W}_\sigma \\
\underline{F} &= \underline{F}_n + \lambda^{me}(t) \underline{W}_F \\
T &= T_n + \lambda^{th}(t) W_T \\
\phi &= \phi_n + \lambda^{th}(t) W_\phi
\end{aligned}
$$

The couple $(\underline{W}_\sigma, \underline{W}_F)$ must satisfy the equilibrium conditions (1). Because the residuals are space functions, only a weak form of the equations of the correction stage can be satisfied. The main equations of this weak form are:

$$ \int_{0}^{t_f} \dot{\lambda}^{me} \dot{\underline{\sigma}} d\tau = \int_{0}^{t_f} \dot{\lambda}^{me} \underline{K}(T_n) (\dot{\underline{\varepsilon}} - \underline{\alpha}(T_n) \dot{\underline{T}}_n) d\tau \quad [13] $$

$$ \int_{0}^{t_f} \dot{\lambda}^{me} \underline{\dot{F}} d\tau = \int_{0}^{t_f} \dot{\lambda}^{me} \underline{\dot{F}}_{n+1/2} - \hat{k}^{me} (\underline{\dot{U}} - \underline{\dot{U}}_{n+1/2}) d\tau \quad [14] $$

Hence $\underline{W}_U$ is computed to satisfy a linear system:

$$ 
\begin{aligned}
\underline{W}_U &= \underline{N}_U \underline{q}_W \\
\underline{K}_W^{me} \underline{q}_W &= \underline{F}_W^{me}(t_f)
\end{aligned}
$$

The residuals $\underline{W}_\sigma$ and $\underline{W}_F$ are deduced from $\underline{W}_U$ by using the equations (13) and (14). The sign of $\dot{\lambda}^{me}$ is chosen such that $\|\underline{F}_W^{me}(t)\|$ is an increasing time function, like it has been proposed in [PEL 00].

The same approach is used to compute $W_T$ and $W_\phi$:

$$ 
\begin{aligned}
W_T &= \underline{N}_T \underline{\theta}_W \\
\underline{K}_W^{th} \underline{\theta}_W &= \underline{F}_W^{th}(t_f) \\
W_\phi &= \int_{0}^{t_f} \dot{\lambda}^{th} \underline{N}_T^T (\phi_{n+1/2} - \phi_n + \hat{h}^{th} (T_n - T_{n+1/2})) d\tau
\end{aligned}
$$
The residuals can be used as basis functions only if they have significant values. If it is the case they are normalised.

### 3.4. The LATIN-KL subspace

As the number of correction stages increases, the size of the LATIN subspace increases also. More over the fields $W_U$, $W_σ$, $W_F$, $W_T$ and $W_φ$ do not have any physical meaning. Since we know an approximate state evolution ($\tilde{U}_n$, $\tilde{\sigma}_n$, $\tilde{F}_n$, $\tilde{T}_n$, $\tilde{\phi}_n$) we can use the Karhunen-Loève expansion [SIR 91] to define a subspace with a better physical sense. Therefore we propose to use the Karhunen-Loève method to extract a subspace from the LATIN subspace. This subspace is the LATIN-KL subspace. For each variable $\tilde{x}$ represented with the basis functions of the LATIN subspace, we use the coefficients of the linear combination $\tilde{x}(t)$ of basis functions placed in a matrix $X$ and a particular state evolution $\tilde{x}_o$:

$$ \tilde{x}(t) = X \tilde{x}(t) + \tilde{x}_o $$

An approximate value of each variable is known at the end of each time step. Hence we know $\tilde{x}$ over the entire time interval. To construct the LATIN-KL subspace, we apply the Karhunen-Loève expansion to the different variables $\tilde{x}(t)$. The values of $\tilde{x}$ on each time step define a smaller subspace than the LATIN subspace. The size of this subspace can't be higher than the number of time steps. For example, if there is only one time step the size of the LATIN-KL subspace is one. The Karhunen-Loève expansion is used to choose the principal directions $\tilde{v}_i$ of the subspace built with the values $\tilde{x}(t)$. These principal directions are defined by the following eigensystem:

$$ M \cdot \tilde{v}_i = \mu_i \tilde{v}_i \quad \text{with} \quad M = \int_0^{t_f} \dot{\tilde{x}}(t) \dot{\tilde{x}}^T(t) dt \quad [15] $$

Let's note $\mu_{max}$ the highest eigenvalue of $M$. Only a part of the principal vectors is saved in a matrix $V$. These vectors are corresponding to the eigenvalues greater than $10^{-6} \mu_{max}$. Hence we obtain a matrix $X_v$ containing the basis functions of the LATIN-KL subspace such that:

$$ X_v = X V \quad [16] $$

In practice, only the basis functions of the LATIN-KL subspace are saved. For each sub-domain $\Omega$, the different matrix of basis functions are: $A_U$ for the displacements, $A_σ$ for the stresses, $A_F$ for the normal stresses over $\partial\Omega$, $A_T$ for the temperatures and $A_φ$ for the normal heat transfer over $\partial\Omega$. Therefore the known approximate state is represented with few degrees of freedom $a_{U n}$, $a_{\sigma n}$, $a_{T n}$ and $a_{\phi n}$, to define the linear combinations:
$$ 
\begin{aligned}
\underline{U}_n &= \mathbf{A}_U \underline{a}_U n(t) + \underline{U}_o \quad \sigma_n = \mathbf{A}_\sigma \underline{a}_\sigma n(t) + \sigma_o \quad \underline{F}_n = \mathbf{A}_F \underline{a}_\sigma n(t) + \underline{F}_o \\
T_n &= \mathbf{A}_T \underline{a}_T n(t) + T_o \quad \phi_n = \mathbf{A}_\phi \underline{a}_\phi n(t) + \phi_o
\end{aligned}
$$

The same degrees of freedom $\underline{a}_\sigma n(t)$ are used for the stresses and the normal stresses on the boundary in order to obtain a couple $(\sigma_n, \underline{F}_n)$ that satisfy the equilibrium conditions (1) at any time instant. So, for stresses, the matrix $\mathbf{V}$ in equation (16) is obtained from $\underline{a}_\sigma n$. But this matrix modifies both $\mathbf{A}_\sigma$ and $\mathbf{A}_F$. After each correction stage, $\underline{W}_U$, $\mathbf{W}_\sigma$, $\underline{W}_F$, $W_T$ and $W_\phi$ are added to the convenient matrix $\mathbf{A}_x$, if we keep a set of independent space functions. The Karhunen-Loève simplification is done only once after three correction stages, before adding the last computed vectors $W_x$. When the basis functions have been modified, we need to adapt their linear combination. To do so, we use the Moore-Penrose inverse of $\mathbf{V}$:

$$ \tilde{\underline{x}}_v(t) = (\mathbf{V}^T \mathbf{V})^{-1} \mathbf{V}^T \tilde{\underline{x}}(t) \quad \forall t \in [0, t_f] $$

### 3.5. Improvement of the approximate state over $[0, t_f]$

To find a correction of the approximate state with the basis functions of the reduced-order model we use a classical weak formulation of the equations of the correction stage. The main variables of the correction stages are the displacements and the temperatures. From the corrections of these variables we deduce the corrections of the stresses and the heat transfer. During the correction stage, the correction $(\delta\underline{U}, \delta\sigma, \delta\underline{F}, \delta T, \delta\phi)$ belong to the known LATIN-KL subspace. So the few variables of this global stage are $\delta\underline{a}_U$, $\delta\underline{a}_\sigma$, $\delta\underline{a}_T$, and $\delta\underline{a}_\phi$ such that:

$$ 
\begin{aligned}
\delta\underline{U} &= \mathbf{A}_U \delta\underline{a}_U(t) \quad \delta\sigma = \mathbf{A}_\sigma \delta\underline{a}_\sigma(t) \quad \delta\underline{F} = \mathbf{A}_F \delta\underline{a}_\sigma(t) \\
\delta T &= \mathbf{A}_T \delta\underline{a}_T(t) \quad \delta\phi = \mathbf{A}_\phi \delta\underline{a}_\phi(t)
\end{aligned}
$$

A weak form of the equations of the correction stage provides the four following incremental problems:

$$ 
\begin{cases}
\delta \underline{a}_U(0) = 0 \\
\mathbf{G}_U^{me}(t_{j+1}) \delta \underline{a}_U(t_{j+1}) = \underline{N}_U^{me}(t_{j+1}) \quad \forall t_{j+1}
\end{cases}
$$
$$ 
\begin{cases}
\delta \underline{a}_T(0) = 0 \\
\mathbf{G}_T^{th}(t_{j+1}) \delta \underline{a}_T(t_{j+1}) = \underline{N}_T^{th}(t_{j+1}) \quad \forall t_{j+1}
\end{cases}
$$
$$ 
\begin{cases}
\delta \underline{a}_\sigma(0) = 0 \\
\mathbf{G}_\sigma^{me}(t_{j+1}) \delta \underline{a}_\sigma(t_{j+1}) = \underline{N}_\sigma^{me}(t_{j+1}) \quad \forall t_{j+1}
\end{cases}
$$
$$ \left\{ \begin{array}{l} \delta \underline{a}_{\phi}(0) = 0 \\ \mathbf{G}_{\phi}^{th}(t_{j+1}) \delta \underline{a}_{\phi}(t_{j+1}) = \underline{N}_{\phi}^{th}(t_{j+1}) \quad \forall t_{j+1} \end{array} \right. $$

### 3.6. Algorithm of the linear correction stage

The local computation of $(\underline{U}_{n+1/2}, \underline{F}_{n+1/2}, T_{n+1/2}, \phi_{n+1/2})$ is both necessary for the computation of the residuals and for the computation of the correction of the approximate state. More over, this local computation and the global computations are incremental ones. Hence we use the following algorithm:

- For each sub-domain $\Omega$: we do a copy of the boundary variables $(\underline{U}'_n, \underline{F}'_n, T'_n, \phi'_n)$ of the neighbour sub-domain $\Omega'$, over the entire time interval $[0, t_f]$

- For each sub-domain $\Omega$

- Initial conditions

- For each time step $[t_j, t_{j+1}]$

For each integration point inside $\Omega$: local contributions to

$$ \begin{array}{l} \mathbf{G}_{\mathrm{U}}^{\mathrm{me}}(t_{j+1}), \underline{N}_{\mathrm{U}}^{\mathrm{me}}(t_{j+1}), \mathbf{G}_{\sigma}^{\mathrm{me}}(t_{j+1}), \underline{N}_{\sigma}^{\mathrm{me}}(t_{j+1}), \mathbf{G}_{\mathrm{T}}^{\mathrm{th}}(t_{j+1}), \\ \underline{N}_{\mathrm{T}}^{\mathrm{th}}(t_{j+1}), \underline{\sigma}, \mathbf{K}_{\mathrm{W}}^{\mathrm{me}}, \mathbf{K}_{\mathrm{W}}^{\mathrm{th}} \end{array} $$

For each integration point on $\partial\Omega$: local contributions to

$$ \begin{array}{l} \mathbf{G}_{\mathrm{U}}^{\mathrm{me}}(t_{j+1}), \underline{N}_{\mathrm{U}}^{\mathrm{me}}(t_{j+1}), \mathbf{G}_{\sigma}^{\mathrm{me}}(t_{j+1}), \underline{N}_{\sigma}^{\mathrm{me}}(t_{j+1}), \mathbf{G}_{\mathrm{T}}^{\mathrm{th}}(t_{j+1}), \\ \underline{N}_{\mathrm{T}}^{\mathrm{th}}(t_{j+1}), \mathbf{G}_{\phi}^{\mathrm{th}}(t_{j+1}), \underline{N}_{\phi}^{\mathrm{th}}(t_{j+1}), \underline{F}, \mathbf{K}_{\mathrm{W}}^{\mathrm{me}}, \mathbf{K}_{\mathrm{W}}^{\mathrm{th}} \end{array} $$

Computation of: $\delta \underline{a}_{\mathrm{U}}(t_{j+1}), \delta \underline{a}_{\mathrm{T}}(t_{j+1})$

For each integration point inside $\Omega$: local contributions to

$$ \mathbf{G}_{\sigma}^{\mathrm{me}}(t_{j+1}), \underline{N}_{\sigma}^{\mathrm{me}}(t_{j+1}) $$

For each integration point on $\partial\Omega$: local contributions to

$$ \mathbf{G}_{\phi}^{\mathrm{th}}(t_{j+1}), \underline{N}_{\phi}^{\mathrm{th}}(t_{j+1}) $$

Computation of: $\delta \underline{a}_{\sigma}(t_{j+1}), \delta \underline{a}_{\phi}(t_{j+1})$

Choice of the sign of $\lambda^{\mathrm{me}}$ and $\lambda^{\mathrm{th}}$ for the contributions to

$$ \underline{F}_{\mathrm{W}}^{\mathrm{me}}(t_{j+1}), \underline{F}_{\mathrm{W}}^{\mathrm{th}}(t_{j+1}), W_{\phi} $$

- Computation of new space functions $\underline{W}_{\mathrm{U}}$, $\underline{W}_{\sigma}$, $\underline{W}_{\mathrm{F}}$ and $W_T$

- For each sub-domain $\Omega$: update the approximate state $\underline{a}_{x n+1} = \underline{a}_{x n} + \delta \underline{a}_{x}$ and eventually update the LATIN-KL subspace

There is no time synchronisation for the computation of the corrections on each sub-domain. This algorithm could be easily used on a parallel computer.

## 4. 2D examples

### 4.1. A first casting problem

Let's consider the Finite Element model on the Figure 1. 40 time steps are used to describe the evolution of the state variables. The two sub-domains $\Omega_1$ and $\Omega_2$ are
elastic bodies, respectively in aluminium and sand. There is no initial gap between the two sub-domains. This is a plan strain problem.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Meshes on each sub-domain</figcaption></figure>

The way $\Omega_1$ is cooled down depends on the contact conditions between the two sub-domains. More over these contact conditions are defined by the thermal expansion of each sub-domain. The thermo-mechanical coupling is strong and similar to the one encountered in casting problems.

## 4.2. Few stages to obtain the main part of the state evolution

To forecast the defaults in a process with a simulation, it is interesting to obtain a quick overview of the evolution of the state variables. After 20 correction stages a good overview of the displacements and temperature fields is obtained (Figure 2).

Even accurate local results are quickly available, with few corrections stages. On Figure 3 there is the example of the temperature evolution (Figure 1) on the point P obtained after different correction stages.

## 4.3. Discussion about the reduced-order model and the LATIN-KL subspace

The LATIN subspace is very closed to the subspace proposed by P. Bussy to represent the corrections of the state variables. The Karhunen-Loève expansion allows a very nice reduction of the number of basis functions of the reduced-order model. On Figure 4 we can observe, for the displacement defined over $\Omega_2$, the difference between the size of the LATIN-KL subspace and the size of the LATIN subspace, during the correction stages.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** An interesting overview of the displacement and temperature field after 20 correction stages</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Different local evolutions of the temperature at the point P</figcaption></figure>

On Figure 5 there is the first basis functions of $A_T$ and the corresponding values of the degree of freedom $a_{T,n}$ after 100 correction stages (n = 100). These basis functions allow to represent 80% of the temperature field over the entire time interval.

In Table 1, we summarised the number of unknowns of the reduced-order model obtained after 100 correction stages, and the number of unknowns of the Finite Element model.

A more complex geometrie have been studied for the same kind of problem. The system is split into four sub-domains (Figure 6). The numbers of unknowns of the reduced-order model for each sub-domain are given in Table 2.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Sizes of the LATIN subspace and LATIN-KL subspace during the correction stages</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Basis functions for the temperature fields and corresponding degrees of freedom (n=100)</figcaption></figure>

Despite a representation of the displacement and the stresses is needed for our approach, only few unknowns are necessary to obtain an accurate approximation of the state evolution.

**5. Conclusion**

Thanks to the proposed reduction method we can quickly obtain an approximate evolution of the state of the studied system. Therefore, if there is any important default on the process, we can see it very quickly without any specific reflection to construct a simple model to describe what is happening in the case studied. The main phenomena are described automatically by the reduction method. So we suggest that very fine Finite Element model should be realised to describe what could happen during the transformation. The reduction method will automatically extract the basis functions to describe the evolution of the state with few unknowns. To perform accurate simulations, we just have to continue correction stages as long as it is necessary. This ap-
<table><caption>**Table 1.** Number of unknowns for each variables over each sub-domain, for the reduced-order model (n=100) / for the Finite Element model on Figure 1</caption><thead><tr><th></th><th>U</th><th>σ</th><th>T</th><th>φ</th></tr></thead><tbody><tr><th>Ω<sub>1</sub></th><td>7 / 636</td><td>8 / 1647</td><td>10 / 318</td><td>6 / 47</td></tr><tr><th>Ω<sub>2</sub></th><td>12 / 926</td><td>9 / 2532</td><td>11 / 463</td><td>6 / 47</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** A more complex problem : Meshes and sub-domains</figcaption></figure>

proach should be interesting to study processes involving multiphysic coupling which are difficult to describe with a simple model.

## 6. References

[ART 92] ARTZ M., COGNARD J.-Y., LADEVÈZE P., “An efficient computational method for complex loading histories”, *Comp. Plas. Fund. and Appl.*, Editions Owen, Onate and Hinton, Barcelona, 1992, p. 225-236.

[BUS 90] BUSSY P., ROUGÉE P., VAUCHEZ P., “The Large Time INcrement method for numerical simulation of metal forming processes”, *Proc. NUMETA*, Editions Elsevier, 1990, p. 102-109.

[CHA 96] CHAMPANEY L., COGNARD J.-Y., DUREISSEIX D., LADEVÈZE P., “Numerical experimentations of parallel strategies in structural non-linear analysis”, *Calc. parallèles.*, vol. 8, num. 2, 1996, p. 245-249.

[KNO 95] KNOLL D.-A., MCHUGH P.-R., “Newton-Krylov methods applied to a system of convection-diffusion-reaction equations”, *Comput. phys. commun.*, vol. 88, 1995, p. 141-160.

[LAD 85] LADEVÈZE P., “Sur une famille d’algorithmes en mécanique des structures”, *CRAS, Paris, 300, série II*, vol. 2, 1985, p. 41-44.

[LAD 96] LADEVÈZE P., *Mécanique non linéaire des Structures*, Etudes en mécanique des matériaux et des structures, Editions Hermès, Paris, 1996.
<table><caption>**Table 2.** Number of unknowns for each variables over each sub-domain, for the reduced-order model (n=100) / for the Finite Element model on Figure 6</caption><thead><tr><th></th><th>U</th><th>σ</th><th>T</th><th>φ</th></tr></thead><tbody><tr><td>Ω<sub>1</sub></td><td>15 / 720</td><td>30 / 2468</td><td>20 / 360</td><td>19 / 94</td></tr><tr><td>Ω<sub>2</sub></td><td>26 / 1422</td><td>32 / 5200</td><td>17 / 711</td><td>22 / 242</td></tr><tr><td>Ω<sub>3</sub></td><td>16 / 1162</td><td>32 / 4144</td><td>23 / 583</td><td>11 / 154</td></tr><tr><td>Ω<sub>4</sub></td><td>19 / 884</td><td>30 / 3084</td><td>26 / 442</td><td>15 / 134</td></tr></tbody></table>

[OCA 89] O'CALLAHAN J., AVITABILE P., RIEMER R., "System Equivalent Reduction Expansion Process", *7th International Modal Analysis Conference*, Las Vegas, 1989, p. 29-37.

[PAR 01] PARK H.-M., JUNG W.-S., "The Karhunen-Loève Galerkin method for the inverse natural convection problems", *Int. j. heat mass transfer*, vol. 44, 2001, p. 155-167.

[PEL 00] PELLE J-P., RYCKELYNCK D., "An Efficient Adaptive Strategy to Master the Global Quality of Viscoplastic Analysis", *Computers & Structures*, Editions Elsevier Science, vol. 78, num. 1-3, 2000, p. 169-184.

[RYC 00] RYCKELYNCK D., "A LATIN Approach to Solve Thermomechanical Problems", *Techniques and Developments*, Editions B.-H.-V. Topping, Civil-Comp Press, 2000, p. 237-248.

[SIR 91] SIROVICH L., "Empirical Eigen functions and low dimensional systems", *New Perspectives in Turbuence*, vol. 5, num. 139, 1991.

[SUN 01] SUNDAR S., BHAGAVAN B.-K., PRASAD S., "Newton-Preconditioned Krylov Subspace Solvers for System of Nonlinear Equations: A Numerical Experiment", *Applied Mathematics Letters*, vol. 14, 2001, p. 195-200.

[WEI 95] WEISS R., "A theoretical overview of Krylov subspace methods", *Appl. numer. math.*, vol. 19, 1995, p. 207-233.