# Efficient weakly coupled projection basis for the reduction of thermo-mechanical models

Ph. Nachtergaele, Daniel Rixen, A.M. Steenhoek

► To cite this version:

Ph. Nachtergaele, Daniel Rixen, A.M. Steenhoek. Efficient weakly coupled projection basis for the reduction of thermo-mechanical models. Journal of Computational and Applied Mathematics, 2010, 234 (7), pp.2272-2278. 10.1016/j.cam.2009.08.116. hal-01517953

HAL Id: hal-01517953

https://hal.science/hal-01517953

Submitted on 4 May 2017

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
# Efficient weakly coupled projection basis for the reduction of thermo-mechanical models

Ph. Nachtergaele<sup>a</sup>, D.J. Rixen<sup>b</sup>, A.M. Steenhoek<sup>b,a</sup>

Open Engineering, Rue des Chasseurs Ardennais, 8, B-4031 Liège, Belgium

<sup>b</sup> Delft University of Technology, Mekelweg 2, 2628CD Delft, The Netherlands

In the work presented in this contribution we have investigated the reduction of dynamical thermo-mechanical systems and we propose a new reduction method that can be seen as an extension of the common Craig-Bampton method Craig and Bampton (1968) [2] where multi-physics is now implicitly included in the projection basis. The efficiency of this new approach has been evaluated in multiple configuration of a representative thermo-mechanical model.

Keywords:

Reduction

Craig-Bampton

Thermo-mechanics

Coupling

## 1. Introduction

Reducing complex models of multi-physical models is an important step in the design and optimization of high-tech systems. For instance, one is interested in representing the dynamical behaviour of components with low order models in order to integrate these models in a functional analysis of the whole system or to apply optimization strategies. The principle of model reduction has been established long time ago [1,2] and substructuring is widely used in structural mechanics. However, novel applications require the reduction methods to be generalized and enhanced in order to efficiently address multidisciplinary problems. In this paper we will discuss an efficient way to build a reduction basis for thermo-mechanical systems. First we shortly summarize the concepts of model reduction and component mode synthesis. Then we outline the thermo-mechanical equations and discuss a special enrichment of the internal modes in order to efficiently represent the dynamic response of a structure to a thermal load. Finally we investigate the efficiency of the proposed basis on simple model of heat actuator.

## 2. Concept of model reduction

Let us consider a finite element model characterized by the second order differential matrix equation in which $q$ represents the vector of degrees of freedom (DOF), $K$, $D$ and $M$ respectively the stiffness, damping and inertia matrices, and $f$ the vector of forces

$$M\ddot{q} + D\dot{q} + Kq = f. \qquad (1)$$

The driving idea of reduction methods is to make an approximation on the model DOF and express them in terms of a linear combination of representative modes. It is clear that in order to obtain a reduced model, the number of modes must be much smaller than the number of DOF in the model. The previous system is then rewritten in terms of "generalized" DOF that
physically correspond to the amplitude of the considered modes. Mathematically this consists in building an approximation of the DOF vector according to a "trial" subspace defined by the modes, assumed to be linearly independent.

$$q = R\hat{q}. \qquad (2)$$

In this approximation, R is a matrix which columns contain the vectors spanning the approximation space. Replacing this expression in the full set of equations leads to the following equation

$$MR\ddot{q} + DR\dot{q} + KR\hat{q} = f + r. \qquad (3)$$

The residual load appears because, due to the approximation on the DOFs, the equation cannot be satisfied exactly in general (unless the approximation space represented by R spans the exact solution, which is never the case in practice). It can be noted that once an approximate solution is computed (as described below) the equation above can always be used to evaluate the associated error in the exact problem. Clearly, there are more equations to satisfy than degrees of freedom in the reduced space. Obviously the idea is to find the set of reduced DOF that minimizes the residue. In the Petrov-Galerkin method one wants to ensure that the residue is null when projected in a given subspace of the problem. Hence in order to get an approximation of the problem, the equation is also projected onto a "test" subspace

$$\tilde{R}^T \{ MR\ddot{q} + DR\dot{q} + KR\hat{q} \} = \tilde{R}^T \{ f + r \} \qquad (4)$$

and the equation governing the behaviour of the reduced model is obtained by requesting the projection of the residue to be zero

$$\hat{M}\ddot{q} + \hat{D}\dot{q} + \hat{K}\hat{q} = \hat{f} \qquad (5)$$

$$\hat{M} = \tilde{R}^T MR \qquad \hat{D} = \tilde{R}^T DR \qquad \hat{K} = \tilde{R}^T KR \qquad \hat{f} = \tilde{R}^T f. \qquad (6)$$

The choice of representative trial and test bases is the key of the method and is directly related to the quality of the approximation. Indeed, if the trial base defines a subspace that contains the exact solution of the problem, then the residual term is null and the reduced model is exact. On the other hand one needs to use a test bases such that the error remaining after solving the reduced problem (the error being in the space orthogonal to $\tilde{R}$) does not significantly alter the solution. Often (such as in symmetric problems) the test and trial bases are taken as identical. This case corresponds to a Galerkin technique (equivalent to a Rayleigh-Ritz method for dynamic problems) and can be interpreted as a virtual work principle (see [3]).

## 2.1. Static and dynamic condensation

The well-known Guyan-Irons static condensation algorithm [1], can be seen as a projection method using static deformation shapes as trial and test bases. Nevertheless, in the case of dynamic analysis, the dynamic behaviour of the internal domain of the model, that means the condensed DOF, is not taken into account. The improvement proposed in [2] consists in adding to the trial and test bases a given number of modes directly related to the dynamic behaviour of the internal part of the model. Craig proposed to consider the internal eigen modes $\phi$ of the model, computed by constraining the retained DOF to zero. These modes indeed represent the dynamic neglected in the quasi-static modes. In a general way, the modes added to the trial and test bases, respectively $\Phi$ and $\tilde{\Phi}$, can be different. Nevertheless, when using internal eigen modes, the common practice consists in considering the same modes in the two bases. In the Craig-Bampton method the following trial and test bases are considered

$$q = \begin{pmatrix} q_r \\ q_c \end{pmatrix} = R \begin{pmatrix} q_r \\ \gamma \end{pmatrix} = R\hat{q} \qquad (7)$$

$$R = \begin{bmatrix} I & 0 \\ -K_{cc}^{-1} K_{cr} & \Phi \end{bmatrix} \qquad \tilde{R}^T = \begin{bmatrix} I & -K_{rc} K_{cc}^{-1} \\ 0 & \tilde{\Phi}^T \end{bmatrix} \qquad (8)$$

where we denote with a subscript c the DOF that are condensed (i.e. the ones internal to a substructure) and by r the remaining ones (the DOF on an interface).

## 3. Reduction of thermo-mechanical models

In the case of thermo-mechanical problems, the DOF representing the model contains structural and thermal DOF, respectively noted $\ddot{u}$ and $\dot{u}$, and reordering the equations according to the physical nature of the DOF leads to the classical formulation of thermo-mechanical problems

$$\begin{bmatrix} M_{uu} & 0 \\ 0 & 0 \end{bmatrix} (\ddot{u} \quad \ddot{\theta}) + \begin{bmatrix} D_{uu} & 0 \\ D_{\theta u} & D_{\theta \theta} \end{bmatrix} (\dot{u} \quad \dot{\theta}) + \begin{bmatrix} K_{uu} & K_{u\theta} \\ 0 & K_{\theta\theta} \end{bmatrix} (u \quad \theta) = \begin{pmatrix} f_u \\ f_\theta \end{pmatrix}. \qquad (9)$$
In the previous expression, $M_{uu}$, $D_{uu}$ and $K_{uu}$ correspond to the structural inertia, damping and stiffness matrices, that are symmetric. The $D_{\theta\theta}$ and $K_{\theta\theta}$ matrices are symmetric, and, in the context of second order equations in time, can be seen as damping and stiffness terms. Physically however they describe respectively the heat capacitance and conduction effects. The two remaining matrices $D_{\theta u}$ and $K_{\theta u}$ characterize the thermo-mechanical coupling. Due to these terms, the complete system is unsymmetric. Partitioning again the DOF vector in terms of the retained part $q_r$ and condensed part $q_c$, each of them containing both structural and thermal DOF, the equations of the complete system can thus be reordered and rewritten

$$ \begin{bmatrix} M_{rr}^{uu} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & M_{cc}^{uu} & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} \begin{pmatrix} \ddot{q}_r^u \\ \ddot{q}_r^\theta \\ \ddot{q}_c^u \\ \ddot{q}_c^\theta \end{pmatrix} + \begin{bmatrix} D_{rr}^{uu} & 0 & D_{rc}^{uu} & 0 \\ D_{rr}^{\theta u} & D_{rr}^{\theta \theta} & D_{rc}^{\theta u} & D_{rc}^{\theta \theta} \\ D_{cr}^{uu} & 0 & D_{cc}^{uu} & 0 \\ D_{cr}^{\theta u} & D_{cr}^{\theta \theta} & D_{cc}^{\theta u} & D_{cc}^{\theta \theta} \end{bmatrix} \begin{pmatrix} \ddot{q}_r^u \\ \ddot{q}_r^\theta \\ \dot{q}_c^u \\ \dot{q}_c^\theta \end{pmatrix} + \begin{bmatrix} K_{rr}^{uu} & K_{rr}^{u\theta} & K_{rc}^{uu} & K_{rc}^{u\theta} \\ 0 & K_{rr}^{\theta u} & 0 & K_{rc}^{\theta \theta} \\ K_{cr}^{uu} & K_{cr}^{u\theta} & K_{cc}^{uu} & K_{cc}^{u\theta} \\ 0 & K_{cr}^{\theta \theta} & 0 & K_{cc}^{\theta \theta} \end{bmatrix} \begin{pmatrix} q_r^u \\ q_r^\theta \\ q_c^u \\ q_c^\theta \end{pmatrix} = \begin{pmatrix} f_r^u \\ f_r^\theta \\ f_c^u \\ f_c^\theta \end{pmatrix}. \quad (10) $$

### 3.1. Static projection modes

We want to stress that, unlike for purely structural models, the static modes in the trial and test matrices are not identical when handling thermo-mechanical models. Indeed, in (8) the block $-K_{cc}^{-1}K_{cr}$ is not simply the transpose of $K_{cr}K_{cc}^{-1}$ since $K_{cc}$ is not symmetric and since $K_{rc}$ is not the transpose of $K_{cr}$. Also a particular attention must be paid to the structure of the matrices in order to implement an efficient reduction algorithm. Indeed, it is preferable to factorize both $K_{cc}^{uu}$ and $K_{cc}^{\theta\theta}$ rather than factorize the entire $K_{cc}$ matrix. It must also be noted that despite the symmetry of the mass matrix in the complete model, the reduced mass matrix of a thermo-mechanical system is non-symmetric. This comes from the fact that the test and trial spaces are not identical.

### 3.2. Internal projection modes

The choice for the internal modes could be based on the state-space form of the internal problem, similarly to what was proposed for highly damped structures [4]. However, our main concern was to deal with systems composed of a structural and a thermal part that have well separate time constants, namely the period of the lowest pure structural frequencies is much smaller than the time constant of the thermal diffusion in the system. In this case, a quasi-static behaviour of the structure can be assumed as a good approximation. Therefore, we conjectured that the dynamic coupling between thermal and structural problem can be neglected in the process of internal mode computation. This last one is thus performed by solving two decoupled eigen problems. The structural symmetric eigen problem characterizes the dynamic response of the system to structural time dependent loads.

$$ K_{cc}^{uu} \Phi_{u}^{u} = M_{cc}^{uu} \Phi_{u}^{u} \Omega_{u}. \quad (11) $$

The “thermal-driven” internal modes are computed in two steps. First, the decoupled thermal eigen problem, also symmetric, that is representative of the thermal response to thermal loading is solved

$$ K_{cc}^{\theta\theta} \Phi_{\theta}^{\theta} = D_{cc}^{\theta\theta} \Phi_{\theta}^{\theta} \Omega_{\theta}. \quad (12) $$

By doing this, the coupling between the two physics is not considered, leading to poor quality approximations. Therefore, in accordance with the assumption of quasi-static behaviour of the structure with respect of the thermal excitations, a correction is proposed. This new approach, which is innovative regarding the current practice, consists in adding a coupling term $\Phi_{\theta}^{u}$ to the thermal internal modes that corresponds to the static structural response to the thermal eigen modes.

$$ \Phi_{\theta}^{u} = -(K_{cc}^{uu})^{-1} K_{cc}^{u\theta} \Phi_{\theta}^{\theta}. \quad (13) $$

According to the quasi-static approximation, the trial and test basis are defined by the following projection operators in which the last column of the matrices will be called the “thermal-driven modes”

$$ R = \begin{bmatrix} I & 0 & 0 & 0 \\ 0 & I & 0 & 0 \\ -[K_{cc}^{uu} & K_{cc}^{u\theta} & K_{cr}^{uu} & K_{cr}^{u\theta} \\ 0 & K_{cc}^{\theta\theta} & K_{cr}^{\theta\theta} & K_{cr}^{\theta\theta} \end{bmatrix}^{-1} \begin{bmatrix} K_{cc}^{uu} & K_{cr}^{u\theta} \\ 0 & K_{cr}^{\theta\theta} \end{bmatrix} \begin{pmatrix} \Phi_{u}^{u} \\ 0 \end{pmatrix} \begin{pmatrix} \Phi_{\theta}^{u} \\ \Phi_{\theta}^{\theta} \end{pmatrix} \quad (14) $$

$$ \tilde{R}^T = \begin{bmatrix} I & 0 & 0 & 0 \\ 0 & I & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} - \begin{bmatrix} K_{rc}^{uu} & K_{rc}^{u\theta} & K_{cc}^{uu} & K_{cc}^{u\theta} \\ 0 & K_{rc}^{\theta\theta} & K_{cc}^{\theta\theta} & K_{cc}^{\theta\theta} \end{bmatrix} \begin{bmatrix} K_{cc}^{uu} & K_{cc}^{u\theta} \\ 0 & K_{cc}^{\theta\theta} \end{bmatrix}^{-1} \quad (15) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic view of the thermal actuator model.</figcaption></figure>

<table><caption>**Table 1**

Mechanical and thermal time constants of the model.</caption><thead><tr><td>Model length</td><td>Heat capacitance</td><td>Mechanical time constant</td><td>Thermal time constant</td><td>Ratio</td></tr></thead><tbody><tr><td>200 m</td><td>1.64 × 10<sup>6</sup></td><td>1.22929 × 10<sup>1</sup></td><td>5.5022 × 10<sup>8</sup></td><td>10<sup>7</sup></td></tr><tr><td>200 m</td><td>1.64 × 10<sup>3</sup></td><td>1.22929 × 10<sup>1</sup></td><td>5.5022 × 10<sup>5</sup></td><td>10<sup>4</sup></td></tr><tr><td>200 µm</td><td>1.64 × 10<sup>6</sup></td><td>1.22929 × 10<sup>-5</sup></td><td>5.5022 × 10<sup>-4</sup></td><td>10<sup>1</sup></td></tr><tr><td>200 µm</td><td>1.64 × 10<sup>3</sup></td><td>1.22929 × 10<sup>-5</sup></td><td>5.5022 × 10<sup>-7</sup></td><td>10<sup>-2</sup></td></tr></tbody></table>

# 4. Applications

The efficiency of this approach has been evaluated using a simplified model of heat actuator, illustrated on Fig. 1. The model is made of thermo-mechanical volume elements and a symmetry condition has been applied on the lower face so that the motion of the actuator is in the plane of the structure. The two anchors are clamped and an iso-temperature constraint is applied on each of them.

The whole structure is thermally loaded using volume heat source to simulate a Joule heating effect. To simulate the physics, a different value of the heat generation rate is considered in the two arms of the structure. In terms of model reduction, two mechanical DOF and two thermal DOF are retained. The mechanical retained DOF correspond to the in-plane displacements of the actuation point located at the corner of the structure tip, in the symmetry plane. The retained thermal DOF, named T1 and T2, are the anchor's temperatures. Note that the thermal DOF retained for the reduction will no longer be unknowns when the reduced model is used for the analysis since the anchor temperatures are then imposed.

## 4.1. Projection modes

The projection modes corresponding to our reduction approach are presented on Fig. 2. The first two modes correspond to the retained mechanical DOF modes ("static modes"). As these modes represent a static mechanical behaviour, they do not involve any thermal component.

The internal thermal-driven projection modes correspond to the eigen modes of the thermal condensed DOF to which the corresponding quasi-static mechanical response is added. The reduced model involves ten internal thermal-driven modes even if only the four first modes are illustrated.

## 4.2. Transient response

In order to evaluate the efficiency of the proposed reduction basis, transient analyses has been performed for different configuration of the model. A geometrical scaling of the model and two values of the heat capacitance of the material have been considered. Table 1 summarizes the mechanical and thermal time constants of the model in the different configurations of analyses.

For each studied configuration, the lateral displacement of the actuation point is compared for the initial 3D model, the reduced model generated using thermal-driven internal modes and a reduced model generated using pure thermal internal modes (no structural modes are used here). The simulation of the models in meter demonstrates that, even if the ratio of time constants is huge, the thermal-driven approach provides significantly more accurate results only when a low value of heat capacitance is considered, as illustrated on Fig. 3.

The simulation of the models in micron for real value of heat capacitance, presented on Fig. 4, shows that the two reduction methods lead to accurate results, even if the ratio of time constants is small. However, in the same geometrical configuration, if the heat capacitance is significantly reduced the thermal-driven approach provides again more accurate results, as illustrated on Fig. 5.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Projection modes (temperature on deformed shape).</figcaption></figure>

Firstly, these results demonstrate that, even if the ratio between mechanical and thermal time constants correspond to a relevant property in terms of thermoelastic damping [5], this value is not a significant indicator in terms of model reduction strategy. Secondly, for usual values of heat capacitance, the two reduction methods provide similar results very close to the reference values.

However, for low values of heat capacitance, the thermal-driven approach reveal significantly more accurate than the reduction strategy based on uncoupled pure thermal internal modes. Therefore, taking into account the low additional computation cost required by the thermal-driven approach, we suggest using thermal-driven modes when building reduced models of thermo-mechanical systems.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Transient analysis of 200 m model.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Transient analysis of 200 µm model with real heat capacitance.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Transient analysis of 200 µm model with low heat capacitance.</figcaption></figure>
## 5. Conclusion

In this paper we have proposed an extension of the classical Craig-Bampton reduction method to thermo-mechanical systems. This new algorithm introduces an enrichment of the reduction basis allowing a more accurate approximation of the thermo-mechanical response. The idea consists in adding to the pure thermal modes their corresponding quasi-static structural response. It has been shown using a simple model of heat actuator that the quasi-static structural enrichment of the thermal modes can significantly improves the fidelity of the reduced model dynamic response. However, it also appears that the efficiency of this improvement depends on the system configuration that is not only characterized by the ratio between structural and thermal time constants. Further investigations are required on this subject to understand in which configuration the thermal-driven mode reduction method is particularly well suited.

## References

[1] R.J. Guyan, Reduction of stiffness and mass matrices, AIAA Journal 3 (2) (1965) 380.

[2] R.R. Craig, M.C.C. Bampton, Coupling of substructures for dynamic analysis, AIAA Journal 6 (7) (1968) 1313-1319.

[3] M. Géradin, D. Rixen, Mechanical Vibrations: Theory and Application to Structural Dynamics, Wiley & Sons, Chichester, 1997.

[4] R.R. Craig, Y.T. Chung, Generalized substructure coupling procedure for damping systems, AIAA Journal 20 (3) (1982) 442-444.

[5] R. Lifshitz, M.L. Roukes, Thermoelastic damping in micro- and nano-mechanical systems, Physical Review B 61 (8) (2000) 5600-5609.