Model Order Reduction for Temperature-Dependent Nonlinear Mechanical Systems: A Multiple Scales Approach

TL;DR

This paper introduces a multiple scales approach for model order reduction in temperature-dependent nonlinear mechanical systems, effectively capturing slow thermal dynamics and improving accuracy over traditional methods.

Contribution

It develops an adaptive reduction basis that varies with temperature, providing a systematic way to reduce complex thermo-mechanical models.

Findings

Achieves better accuracy than standard Galerkin projection.

Reduces the number of unknowns in the models.

Demonstrates effectiveness on linear and nonlinear beam examples.

Abstract

The thermal dynamics in thermo-mechanical systems exhibits a much slower time scale compared to the structural dynamics. In this work, we use the method of multiple scales to reduce the thermo-mechanical structural models with a slowly-varying temperature distribution in a systematic manner. In the process, we construct a reduction basis that adapts according to the instantaneous temperature distribution of the structure, facilitating an efficient reduction in the number of unknown. As a proof of concept, we demonstrate the method on a range of linear and nonlinear beam examples and obtain a consistently better accuracy and reduction in the number of unknowns than standard the Galerkin projection using a constant basis.

Tables

Table 1. Table 1 : Values for geometrical parameters and physical properties used in the beam models

Parameters (symbol)	Value [unit]
Length of beam ($L$)	0.1 [m]
Thickness of beam ($h$)	1 [mm]
Width of beam ($b$)	10 [mm]
Curvature parameter for the curved-
beam ($w$; cf. Figure 4a; $w=0$ for straight beam)	5 [mm]
Young’s Modulus ($E$)	70 [GPa]
Material damping modulus ($\kappa$)	0.1 (GPa s)
Density ($\rho$)	2700 [kg/m3]
Coefficient of linear expansion (${\alpha }_T$)	23.1$\times {10}^{-6}$ [K-1]

Table 2. Table 2 : The uniform-in-time error estimate ( 51 ) for the different reduction approaches considered in the curved beam with geometric nonlinearities (cf. Figure  12 ).

Reduction method	Basis size ($m$)	Error
Modal	60	9.47%
Modal-POD	20	136.75%
Modal MMS $𝒪\left(1\right)$	20	3.69%
Modal MMS $𝒪\left(ϵ\right)$	20	3.19%

Equations

$$M\ddot{x}+C\dot{x}+K(T)x+=f(t)+g(T),$$

$$M\ddot{x}+C\dot{x}+K(T(t))x=f(t)+g(T(t)),$$

$$u=V(p(t))q(t)$$

$$\mathbf{M}\ddot{\mathbf{u}}+\mathbf{C}\dot{\mathbf{u}}+\mathbf{f}\left(\mathbf{u},\mathbf{T}\right)$$

$$\mathbf{u}(0)$$

$$\dot{\mathbf{u}}(0)$$

$${\mathbf{M}}_T{\mathbf{T}}^{\prime }+{\mathbf{f}}_T\left(\mathbf{T}\right)=\mathbf{h}(\tau ),$$

$$\dot{\mathbf{T}}=ϵ{\mathbf{T}}^{\prime }.$$

$$\mathbf{u}=\mathbf{u}(t,\tau ).$$

$$\mathbf{u}(t,\tau )={\mathbf{u}}_0(t,\tau )+ϵ{\mathbf{u}}_1(t,\tau )+{ϵ}^2{\mathbf{u}}_2(t,\tau )+\dots ,$$

$$\dot{\mathbf{u}}=$$

$$\ddot{\mathbf{u}}=$$

$$\mathbf{M}\left(\frac{{\partial }^2\mathbf{u}}{\partial t^2}+2ϵ\frac{{\partial }^2\mathbf{u}}{\partial t\partial \tau }+{ϵ}^2\frac{{\partial }^2\mathbf{u}}{\partial {\tau }^2}\right)+\mathbf{C}\left(\frac{\partial \mathbf{u}}{\partial t}+ϵ\frac{\partial \mathbf{u}}{\partial \tau }\right)+\mathbf{f}\left(\mathbf{u}(t,\tau ),\mathbf{T}(\tau )\right)=\mathbf{p}\left(t;ϵ\right),$$

$$\text{O}(1):\mathbf{M}\frac{{\partial }^2{\mathbf{u}}_0}{\partial t^2}+\mathbf{C}\frac{\partial {\mathbf{u}}_0}{\partial t}+\mathbf{f}\left({\mathbf{u}}_0,\mathbf{T}(\tau )\right)=\mathbf{p}\left(t,0\right).$$

$$\text{O}(ϵ):\mathbf{M}\frac{{\partial }^2{\mathbf{u}}_1}{\partial t^2}+\mathbf{C}\frac{\partial {\mathbf{u}}_1}{\partial t}+\left[\frac{\partial \mathbf{f}}{\partial \mathbf{u}}\left({\mathbf{u}}_0,\mathbf{T}(\tau )\right)\right]{\mathbf{u}}_1$$

$${\mathbf{u}}_1(0)=\frac{\partial {\mathbf{u}}_1(0)}{\partial t}=0.$$

$$\mathbf{M}\ddot{\mathbf{u}}+\mathbf{C}\dot{\mathbf{u}}+\mathbf{K}({\mathbf{T}}_0)\mathbf{u}+\mathbf{b}({\mathbf{T}}_0)=\mathbf{g}\left(t\right),$$

$${\mathbf{u}}_{eq}$$

$$\mathbf{M}\ddot{\mathbf{y}}+\mathbf{C}\dot{\mathbf{y}}+\mathbf{K}({\mathbf{T}}_0)\mathbf{y}=\mathbf{g}\left(t\right),$$

$$\left(\mathbf{K}({\mathbf{T}}_0)-{\omega }_i^2\mathbf{M}\right){\mathbf{\varphi }}_i=0,i=1,\dots ,n$$

$$\mathbf{y}(t)=\sum _{i=1}^n{\mathbf{\varphi }}_i{q}_i(t)\approx \sum _{i\in M}{\mathbf{\varphi }}_i{q}_i(t),$$

$${\mathbf{\Phi }}^{\top }\mathbf{M\Phi }\ddot{\mathbf{q}}+{\mathbf{\Phi }}^{\top }\mathbf{C\Phi }\dot{\mathbf{q}}+{\mathbf{\Phi }}^{\top }\mathbf{K\Phi q}={\mathbf{\Phi }}^{\top }\mathbf{g}\left(t\right),$$

$$⟹\mathbf{u}(t)\approx {\mathbf{u}}_{eq}+\sum _{i\in M}{\mathbf{\varphi }}_i{q}_i(t)={\mathbf{u}}_{eq}+\mathbf{\Phi }\mathbf{q}(t).$$

$$\mathbf{M}\ddot{\mathbf{u}}+\mathbf{C}\dot{\mathbf{u}}+\mathbf{K}\left(\mathbf{T}(t)\right)\mathbf{u}+\mathbf{b}(\mathbf{T}(t))=\mathbf{g}\left(t\right),$$

$$\mathbf{M}\ddot{\mathbf{u}}+\mathbf{C}\dot{\mathbf{u}}+\mathbf{K}\left(\mathbf{T}(t)\right)\mathbf{u}=\tilde{\mathbf{g}}\left(t\right),$$

$$\mathbf{u}(t)\approx \mathbf{Vq}(t).$$

$$\mathbf{M}\frac{{\partial }^2{\mathbf{u}}_0}{\partial t^2}+\mathbf{C}\frac{\partial {\mathbf{u}}_0}{\partial t}+\mathbf{K}\left(\mathbf{T}(\tau )\right)\mathbf{u}+\mathbf{b}(\mathbf{T}(\tau ))=\mathbf{p}\left(t,0\right).$$

$${\mathbf{u}}_{eq}(\tau )=-\mathbf{K}{\left(\mathbf{T}(\tau )\right)}^{-1}\mathbf{b}(\mathbf{T}(\tau )).$$

$$\left[\mathbf{K}(\mathbf{T}(\tau ))-{\omega }_i^2(\tau )\mathbf{M}\right]{\mathbf{\varphi }}_i(\tau )=0$$

$${\mathbf{u}}_0(t,\tau )\approx {\mathbf{u}}_{eq}(\tau )+\mathbf{\Phi }(\tau ){\mathbf{q}}_0(t).$$

$$\mathbf{\Phi }(\tau {)}^{\top }\left[\mathbf{M}\mathbf{\Phi }(\tau )\frac{{\partial }^2{\mathbf{q}}_0}{\partial t^2}+\mathbf{C}\mathbf{\Phi }(\tau )\frac{\partial {\mathbf{q}}_0}{\partial t}+\mathbf{K}\left(\mathbf{T}(\tau )\right)\left[{\mathbf{u}}_{eq}(\tau )+\mathbf{\Phi }(\tau ){\mathbf{q}}_0(t)\right]+\mathbf{b}(\mathbf{T}(\tau ))\right]$$

$$\mathcal{O}(1):\underset{{\mathbf{M}}_{\mathbf{\Phi }}(\tau )}{\underbrace{\mathbf{\Phi }(\tau {)}^{\top }\mathbf{M}\mathbf{\Phi }(\tau )}}\frac{{\partial }^2{\mathbf{q}}_0}{\partial t^2}+\underset{{\mathbf{C}}_{\mathbf{\Phi }}(\tau )}{\underbrace{\mathbf{\Phi }(\tau {)}^{\top }\mathbf{C}\mathbf{\Phi }(\tau )}}\frac{\partial {\mathbf{q}}_0}{\partial t}+\underset{{\mathbf{K}}_{\mathbf{\Phi }}(\tau )}{\underbrace{\mathbf{\Phi }(\tau {)}^{\top }\mathbf{K}\left(\mathbf{T}(\tau )\right)\mathbf{\Phi }(\tau )}}{\mathbf{q}}_0(t)=\mathbf{\Phi }(\tau {)}^{\top }\mathbf{p}\left(t,0\right)$$

Full text

Institute for Mechanical Systems, ETH Zürich

Leonhardstrasse 21, 8092 Zürich, Switzerland

1 Introduction

Model reduction plays an important role in identifying low dimensional features in high-dimensional dynamical systems, whose numerical simulation would otherwise be either extremely computationally-intensive or entirely infeasible. In such cases, model reduction leads to a small number of variables that describe the system dynamics evolving over low-dimensional subspaces (or manifolds), making computations easier. Structural dynamics applications such as aerospace structures, civil structures, micro-electro-mechanical systems (MEMS) etc., often involve such high-dimensional dynamical systems, which are obtained after spatial discretization of the governing equations of motion. Such structures also exhibit changes in mechanical behavior under the influence of external environmental effects such as thermal, acoustic etc., which manifest as coupling terms in the equations of motion. Model reduction in such coupled domains is especially challenging. We propose a multiple-scales-based technique to reduce structural dynamics equations varying geometrical nonlinear structural dynamics. In this work, we focus on structural systems featuring geometric nonlinearities and thermal stresses induced by temperature fields varying slowly in time and space, and we propose a multiple-scales-based technique to reduce the temperature-dependent equations of motion.

For modeling mechanical structures under the influence of temperature changes, a one-way coupling is often used. This means that the temperature changes affect the structural dynamics but not vice-versa. This is a reasonable assumption because roughly speaking, the thermal dynamics evolves over a much slower time-scale in comparison with the structural dynamics and though the structure can respond to temperature-changes, the same is not true for thermal response to structural motion. More specifically, this slow thermal time scale is reflected in the large time constant of temperature evolution law in comparison with the time-period of oscillation of the structure. Due to such a slow nature of temperature evolution, structural dynamics equations are sometimes even modeled with a constant temperature over a given structural time span. When the temperature does not change during a structural dynamic simulation, reduced-order models (ROM) can be constructed using many standard techniques in literature, e.g., using vibration modes (VMs) [1], modal derivatives [2], proper orthogonal decomposition (POD) [3, 4], implicit condensation and expansion (ICE) [5], among others, cf. ref. [6] for a review.

Nonlinear ROMs are often local in nature, i.e, these are constructed to approximate dynamical behavior locally in the vicinity of an equilibrium, as opposed to globally in the phase space. This is, especially, the case when ROMs are constructed from vibration modes and modal derivatives (cf. refs. [2, 7]), which are computed around an equilibrium and span the relevant subspaces. In the context of temperature dependence, consider a simple linear mechanical system such as

[TABLE]

where $M,C,K$222we introduce these equations/symbols informally for the purpose of the discussion in this introductory section. A more precise treatment follows from Section 2 onwards. are mass, damping and stiffness matrices; $t$ is the time; and $T$ represents temperature variable(s). Usually, model reduction is performed by assuming the temperature $T$ to be a static parameter and a reduction basis (comprising, e.g., a few VMs) is computed for a given value of this static temperature parameter. Interestingly, in specific cases, a basis comprising of VMs of the unheated structure, i.e., cold modes, in combination with some carefully selected dual modes is found to be suitable for capturing the behavior of heated structures [8, 9]. It was further reported in ref. [9] that the use of hot modes in the basis leads to an improvement in results. In general, however, a practical strategy would be to construct a database of reduced-order bases for a set of static temperature fields and perform interpolation among these precomputed bases [10, 11, 12] to come up with ROMs for other values of temperature fields.

Now, if the temperature field is indeed static for a given structural dynamics problem, these approaches may be expected to perform well. However, when the temperature varies dynamically during a structural simulation (e.g., due to change in thermal boundary conditions or heat source), i.e.,

[TABLE]

the aforementioned techniques are not directly applicable for reducing structural dynamics equations due to two main reasons:

Non-existence of an equilibrium point for local reduction: As mentioned above, most reduction techniques usually construct relevant subspaces which are attached to an equilibrium point. However, in the case of dynamic temperature dependence, there is no such (unique) equilibrium. Even in the absence of mechanical loading, the internal force generated due to the dynamically varying temperature leads to a continuously “changing” equilibrium. Then, any local reduction technique would not be applicable since there is no unique equilibrium around which the dynamics can be approximated. 2. 2.

No invariant subspaces for reduction: Each temperature configuration corresponds to a new set of VMs which form relevant invariant subspaces for reduction, if temperature is kept static. However, due to a continuous change in temperature, we cannot find a fixed invariant subspace, even for the linear system (2).

Despite these issues, model reduction for structural equations has been successfully performed with dynamic temperature-dependence in some cases. One simple approach is to obtain a reference temperature configuration by averaging temperature variation over time for each point in space. This results in a static temperature field, around which a reduction basis can be calculated to reduce the problem where the temperature is changing dynamically. While this approach has been demonstrated to perform well in specific situations [13], it is easy to envision scenarios where this simple technique can lead to incorrect predictions. Indeed, consider a beam under the influence of a spatial temperature pulse which is oscillating along the length of the beam (modeling thermal effects of an oscillating shock), as examined in ref. [14]. In this case, the above-mentioned averaging approach, would lead to a spatially uniform temperature distribution, around which the reduction basis would be constructed. Clearly, such an approximation would be inconsistent with the physics of the problem.

Furthermore, the above-mentioned technique of complimenting a basis of cold modes with dual modes would not be theoretically justifiable for reducing (2) due to the two issues above. However, these techniques have been found to work well in specific cases [15, 14] where the structural dynamics is reduced under the influence of time-varying temperature distribution.

In the context of dynamic parameter dependence, another technique [16] involves the construction of a smooth parameter-dependent basis to reduce linear problems with application to moving external loads or boundary conditions, encountered during gear meshing. In this method, a time-dependent mapping from a low-dimensional space is introduced into the full system of unknowns such that

[TABLE]

where $V$ is a basis matrix, parameterized with respect to the parameters $p$, which are changing in time; and $q(t)$ are the reduced set of unknowns. In this setting, a set of bases which are constructed for static parameter values are interpolated to obtain a family of bases $V(p)$, which smoothly depends on the parameters $p$. The time derivative of eq. (3) leads to convective terms in the corresponding ROM, whose evaluation can pose a numerical challenge. It is important to note that unlike the modal subspaces which are invariant for static parameters, eq. (3) describes an invariance relationship which does not hold in general, even for linear systems such as (2).

Indeed, consider a simple two-degree-of-freedom linear mechanical oscillator, as shown in Figure 1. The spring and damping constant in this example are assumed to be temperature-dependent such that the eigen-frequencies of the system remain unchanged and well-separated for any variation in temperature. The eigenvectors, however, change upon varying temperature. While keeping the temperature constant, we force the system near its first natural frequency. This prompts us to reduce the system using the first eigenmode, resulting in a single-degree-of-freedom ROM at any fixed temperature. The same is depicted in Figure 2 for a few of values of temperature.

Since the first mode calculated at different temperatures is found to capture the dynamics of the system in Figure 1, the following questions arise in the case of dynamic temperature variation:

• •

Is it reasonable that a single-mode basis which adapts instantaneously to the dynamic temperature change such as in eq. (3) would be effective for model reduction?

• •

And if yes, when can we expect such a reduction approach to produce good results?

The answer to these question lies in an interesting observation, depicted in Figure 3. We see that an instantaneously-adapting basis with a single mode is effective in capturing system dynamics under a slowly changing temperature. However, the same ROM is not effective when temperature varies rapidly. This shows that a slow variation of parameters is necessary (but not sufficient) to guarantee effective model reduction using interpolation of bases obtained from static parameter values in the context of dynamic parameter dependence.

To this end, we apply the method of multiple scales (MMS) (see, e.g., ref. [17]) to the temperature-dependent structural dynamics equations, where the temperature is taken to be evolving at a different, slower time scale than that of the structural motion. This results in a set of partial differential equations (PDEs) in the slow thermal time and the fast structural time. Due to this treatment, the slowly varying temperature unknowns can be treated as static parameters in a mathematically justifiable manner. Interestingly, this procedure automatically avoids any convective terms arising from the time derivative of the mappings such as (3), which are usually cumbersome to evaluate [16].

Using the MMS, we obtain an asymptotic set of PDEs in the structural time scale, where the temperature can be seen as a slowly varying parameter. This allows for the calculation of a parameter-dependent equilibrium around which the aforementioned local reduction techniques can be applied using the corresponding subspaces in an adaptive manner. Thus, the use of MMS resolves both the issues mentioned above in this slow-fast, thermo-mechanical setting.

After applying the MMS to temperature-dependent structural equations in Section 2, we first show the relevance of MMS on a simple single-degree-of-freedom thermo-mechanical oscillator. Next, we propose the use of a basis which instantaneously adapts to the slowly varying temperature in Section 3. In doing, we preserve minimal number of unknowns in the basis due to changing temperature. In Section 3, we first consider linear systems, reduced using temperature-dependent VMs, followed by geometrically nonlinear systems, reduced using VMs and modal derivatives in an adaptive basis. As a proof-of-concept, we test these propositions in Section 4 on straight and curved beam examples in the linear as well as the geometrically nonlinear setting. Finally, some conclusion with scope for further work are laid down in Section 6.

2 Multiple time scale dynamics

After spatial discretization using finite elements (FE) of partial differential equations (PDEs) governing temperature-dependent momentum balance, we obtain a set of ordinary differential equations (ODEs). Along with the initial conditions for generalized displacements and velocities, these ODEs govern the temperature-dependent response of the underlying structure in the form of an initial value problem as:

[TABLE]

where $\mathbf{u}(t)\in\mathbb{R}^{n}$ is the vector of generalized displacements; $\mathbf{M}\in\mathbb{R}^{n\times n}$ is the mass matrix; $\mathbf{C}\in\mathbb{R}^{n\times n}$ is the damping matrix; $\mathbf{T}\in\mathbb{R}^{n_{T}}$ represents the spatially discretized and possibly time-dependent temperature field prescribed on the mechanical structure, $\mathbf{f}:\mathbb{R}^{n}\times\mathbb{R}^{n_{T}}\mapsto\mathbb{R}^{n}$ gives nonlinear elastic internal force as a function of the displacement $\mathbf{u}$ and temperature $\mathbf{T}$ of the structure; and $\mathbf{g}(t)\in\mathbb{R}^{n}$ is the time-dependent external load vector, $\dot{(\bullet)}=\frac{d(\bullet)}{dt}$ denotes the derivative with respect to time $t$. The system (4) is sometimes also referred to as the high fidelity model (HFM) in the literature.

As mentioned in the Introduction, the temperature field over structures is expected to evolve much slower as compared to the structural time scale $t$. The evolution of the temperature field over the structure is usually modeled using the heat equation, which upon spatial discretization (using, e.g., finite elements) gives

[TABLE]

where $\tau$ is the thermal time-scale; $(\bullet)^{\prime}$ denotes the derivative with respect to $\tau$; $\mathbf{M}_{T}\in\mathbb{R}^{n_{T}\times n_{T}}$ is the thermal mass matrix; $\mathbf{f}_{T}:\mathbb{R}^{n_{T}}\mapsto\mathbb{R}^{n_{T}}$ gives the heat flux through the structure due to physical effects such as conduction, convection, radiation, or a combination there-of; $\mathbf{h}(\tau)$ represents the externally-applied heat source to the structure.

The slowness of the thermal dynamics with respect to the structural one can be mathematically modeled with the relation $\tau=\epsilon t$, where $0<\epsilon\ll 1$ is a small physical parameter. In practice, $\epsilon$ can be quantified as, e.g., the ratio of the time period of the fundamental mode of oscillation of the structural system (4) after linearization and the thermal time constant obtained from the linearized spectral analysis of (5). Now, if the temperature field is a function of the slow time scale $\tau$, i.e., $\mathbf{T}=\mathbf{T}(\tau)$, then we get

[TABLE]

Equation (6) suggests that if $\mathbf{T}$ naturally evolves over the timescale $\tau$, i.e., $\mathbf{T}$ varies at $\mathcal{O}(1)$ speed with respect to $\tau$, then it varies at a much smaller $\mathcal{O}(\epsilon)$ speed over the structural time scale $t$. Physically, the solution to the HFM (4) can be assumed to have a fast component due structural dynamics, as well as a slow component under the influence of slowly changing temperature. In other words, the solution evolves over two timescales as

[TABLE]

The method of multiple scales (cf., e.g., ref. [17]) then allows us to express the solution in terms of an expansion in $\epsilon$ as

[TABLE]

where $\mathbf{u}_{0},\mathbf{u}_{1},\mathbf{u}_{2}\dots$ are the solution components at different orders in $\epsilon$, which can be recursively obtained in the following manner. Using chain rule for the time derivative of eq. (8), we obtain

[TABLE]

Upon substituting (8) and (9) into the governing equations (4), we obtain a PDE in $t$ and $\tau$ given by

[TABLE]

where we have introduced the dependence of the external mechanical force $\mathbf{p}$ on $\epsilon$. We note that there is no reason for the external forcing to be dependent on the physical parameter $\epsilon$. In this perturbative setting, however, we treat $\epsilon$ as scaling constant that could be used to distinguish between components of a given forcing that have different magnitudes. Thus, $\mathbf{p}$ can be seen as an appropriately-scaled version of $\mathbf{g}$ in (4), e.g., if $\mathbf{g}$ is $\mathcal{O}(\epsilon)$ small in comparison with other terms of the governing equations (4), then we can simply write $\mathbf{g}(t)=\mathbf{p}(t,\epsilon)=\epsilon\mathbf{l}(t),$ where $\mathbf{l}(t)$ is $\mathcal{O}(1)$. According to the MMS, eq. (10) must be satisfied at all orders in $\epsilon$. Using the expansion (8) and after Taylor expansion of the nonlinear internal force $\mathbf{f}$ around $\epsilon=0$, we obtain the leading order terms in the PDE 10 as:

[TABLE]

The solution to eq. (11) gives us the leading-order component $\mathbf{u}_{0}$ in eq. (8). Equation (11) can be integrated in time $t$ using usual time integration schemes adopted for eq. (4), e.g., Newmark’s scheme, with the same initial conditions, i.e., $\mathbf{u}_{0}(0)=\mathbf{x}_{0},\frac{\partial\mathbf{u}_{0}(0)}{\partial t}=\mathbf{v}_{0}$. Collecting the $\mathcal{O}(\epsilon)$ terms in eq. (10), we obtain

[TABLE]

Upon solving for $\mathbf{u}_{0}$ from eq.(11), it is easy to see that eq. (12) represents a linear system of equations in $\mathbf{u}_{1}$, the solution to which gives us an $\mathcal{O}(\epsilon)$ correction to $\mathbf{u}_{0}$ as given in eq. (8). Physically, (12) represents the system linearized around the leading-order solution $\mathbf{u}_{0}$ with a forcing term on the right-hand-side (dependent on $\mathbf{u}_{0}$). Since the initial conditions of the HFM (4)) are independent of $\epsilon$, they are included in the leading-order problem (11), shown above. Consequently, the initial conditions for solving eq. (12) should be identically zero, i.e.,

[TABLE]

Thus, using the method of multiple scales, we obtain a series of problems which–as shall be shown later–can be used to systematically reduce the HFM at different orders in $\epsilon$. Equations (11) and (12) give only the equations for $\mathbf{u}_{0}$ and $\mathbf{u}_{1}$. Higher-order terms in eq. (8) can be found by following a similar recursive procedure by collecting higher-order terms in eq. (10). Furthermore, as is the case with the $\mathcal{O}(\epsilon)$ problem (12), it easy to see that all higher problems would be linear in their corresponding unknowns.

3 Adaptive basis for model reduction

3.1 Linear systems: reduction using vibration modes

For the sake of consistency, we first briefly review the concept of modal super-position which is the classic choice for modal reduction of linear mechanical system.

3.1.1 Static temperature dependence

It is well known that modal subspaces (i.e., subspaces spanned by VMs) are invariant and useful in reduction of linear systems. This is also referred to as modal superposition in literature (cf., e.g., Géradin & Rixen [1]). Consider the case of linear structural dynamics independent of temperature or when the temperature configuration remains fixed during a dynamic simulation, i.e., $\mathbf{f}(\mathbf{u},\mathbf{T})=\mathbf{K}(\mathbf{T}_{0})\mathbf{u}+\mathbf{b}(\mathbf{T_{0}})$ in eq. (4) such that

[TABLE]

where $\mathbf{T}_{0}$ is a fixed temperature configuration. In the absence of dynamic excitation $\mathbf{g}(t)$, the above system has a unique equilibrium given as

[TABLE]

Upon shifting the origin to $\mathbf{u}_{eq}$, we get

[TABLE]

where $\mathbf{y}=\mathbf{u}-\mathbf{u}_{eq}$. The VMs of system (14) form a basis for the generalized displacements $\mathbf{y}$ around the equilibrium $\mathbf{u}_{eq}$. For proportionally damped structures, the undamped VMs $\boldsymbol{\phi}_{i}$ that are obtained from the solution of the eigenvalue problem,

[TABLE]

lead to decoupling of the system (14) and are used in a modal truncation approximation as

[TABLE]

where $\mathcal{M}\subset\{1,2,\dots,n\}$ is a set containing the indices of the VMs relevant for modal truncation (cf. Géradin & Rixen  [1]), such that $m:=|\mathcal{M}|\ll n$.

The corresponding ROM is obtained via the classic Galerkin projection as

[TABLE]

where $\boldsymbol{\Phi}\in\mathbb{R}^{n\times m}$ is a matrix containing the VMs indexed by the set $\mathcal{M}$. The solution to the original system (13) can then be recovered as

[TABLE]

3.1.2 Presence of dynamic temperature-dependence

If we allow temperature dependence in an arbitrary fashion, such that

[TABLE]

we see that even in the absence of external forcing $\mathbf{g}(t)$, the system (19) has no equilibrium point due to the presence of dynamic thermal excitation $\mathbf{b}(\mathbf{T}(t))$. Thus, due to the absence of equilibrium and dynamic temperature variation, one cannot use specific VMs to reduce the system, even in a linear system such as (19), as we discussed in the Introduction.

This problem is commonly tackled in the literature (cf. refs. [8, 9, 13, 15, 14]) by simply treating the temperature-dependent term $\mathbf{b}(\mathbf{T}(t))$ as yet another time-dependent forcing term similar to $\mathbf{g}(t)$. Specifically, we obtain the system

[TABLE]

where $\tilde{\mathbf{g}}\left(t\right):={\mathbf{g}}\left(t\right)-\mathbf{b}(\mathbf{T}(t))$. Note that the linear system (20) now exhibits a unique trivial equilibrium point, i.e., $\mathbf{u}_{eq}=\mathbf{0}$ in the absence of forcing $\tilde{\mathbf{g}}\left(t\right)$. This manipulation leads to the routinely performed Galerkin projection about the origin via the mapping

[TABLE]

The reduction basis $\mathbf{V}$ is usually constructed in the literature by an ad hoc selection of modes (hot/cold/dual modes) that are expected to account for the mechanical forcing as well as the dynamic temperature dependence, as discussed in the Introduction.

However, if it is a priori known that the thermal dynamics is slow, as is usually the case for thermo-mechanical systems, we can apply the method of multiple scales, as described in Section 2, to obtain the leading order problems as follows (cf. eq. (11)):

[TABLE]

In this PDE in time $t$, the slow time $\tau$ can be seen as parameter, dependent upon which equilibria can be calculated for the unforced systems as

[TABLE]

Furthermore, the solution to the eigenvalue problem

[TABLE]

results in the temperature-dependent VMs $\boldsymbol{\phi}_{i}(\tau)$ around this manifold of equilibria. The solution to the leading order problem (22) is then approximated using a parameterized basis $\mathbf{\boldsymbol{\Phi}}(\tau)$ containing a truncated set of VMs as

[TABLE]

A reduced-order model can be obtained by projecting the leading order system (22) on to the slowly varying basis as

[TABLE]

and using eq. (23) we obtain the following reduced order model at the leading order

[TABLE]

Furthermore, the $\mathcal{O}(\epsilon)$ correction according to eq. (11) is given by

[TABLE]

Using the same modal basis to reduce this linear problem and observing that origin is the unique fixed point for this system, we get

[TABLE]

which upon Galerkin projection yields the reduced-order model at $\mathcal{O}(\epsilon)$

[TABLE]

Thus, in the linear case, the reduced operators are the same for reduced problems at the leading order, as well as at $\mathcal{\mathcal{O}}(\epsilon)$. In fact, it is easy to see that these linear operators remain the same for higher order corrections, as well.

We would like to emphasize that although this formal procedure can be carried out even if the underlying thermal dynamics is not slow, i.e., cases where $\epsilon$ is not small enough, it would not be mathematically justifiable and can very easily lead to spurious results, as shown in the example in Section 1.

3.2 Geometrically nonlinear systems

For reduction of linear systems, a basis comprising of a truncated set of VMs can be easily identified by the examination of the external forcing and the linear spectrum. However, such a basis is usually not suitable for reduction of geometrically nonlinear systems because of membrane effects, which become significant in the nonlinear regime due to the bending-stretching-torsion coupling. In such cases, modal derivatives have emerged as a simple and effective model-reduction tool (cf. refs. [2, 18, 7, 19, 20]). We discuss the use of modal derivatives in the context of dynamic temperature dependence.

The temperature-dependent equilibrium $\mathbf{u}_{eq}\left(\tau\right)$ is given by the solution to static problem

[TABLE]

We linearize the nonlinear internal force $\mathbf{f}$ around the configuration $\mathbf{u}_{eq}\left(\tau\right)$ to obtain the tangent stiffness matrix $\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\left(\mathbf{u}_{eq}(\tau),\mathbf{T}(\tau)\right)$ evaluated at the temperature-dependent equilibrium $\mathbf{u}_{eq}\left(\tau\right)$. Since $\tau$ is taken as a slowly varying parameter, we can compute the temperature-dependent VMs of the structure similar to the linear case as

[TABLE]

The static modal derivatives (cf. ref. [7]), which are relevant for capturing the membrane effects in thin-walled structures are given as

[TABLE]

where $\boldsymbol{\theta}_{ij}$ is the static modal derivative and physically represents the change in mode $\boldsymbol{\phi}_{i}$ as the structure is perturbed in the direction of mode $\boldsymbol{\phi}_{j}$. Note that here the modal derivatives are calculated around a non-trivial, temperature-dependent equilibrium. This is in contrast with the case of no temperature dependence usually treated in literature, where the modal derivatives are simply constructed around the origin.

The reduction basis is formed by combining a truncated set of modes with the corresponding modal derivatives. As observed in ref. [7], a set of $k$ VMs results in $k(k+1)/2$ modal derivatives. These modal derivatives may be selected to further reduce the basis size (cf. refs. [22, 28]).

[TABLE]

where $\mathbf{V}(\tau)\in\mathbb{R}^{n\times m}$ is the reduction basis. Analogous to the linear case (cf. (25)), we express the leading-order reduced solution as

[TABLE]

Upon substituting (35) into the leading order MMS expansion (11) and performing Galerkin projection, we get

[TABLE]

The $\mathcal{O}(\epsilon)$ correction $\mathbf{u}_{1}$ given in (12) can be reduced using the same adaptive basis as

[TABLE]

And again, by means of Galerkin projection, we obtain a reduced version of (12) which is linear in the reduced unknowns $\mathbf{q}_{1}$ as

[TABLE]

In the special case of a small $\mathcal{O}(\epsilon)$ mechanical forcing, e.g. $\mathbf{p}(t,\epsilon)=\mathbf{\epsilon}\mathbf{l}(t)$ we have

[TABLE]

and the leading-order ROM is given by

[TABLE]

which is the system response to thermal loading only. This would be a good approximation of the system response for small enough $\epsilon$ as long $\mathbf{u}_{eq}$ is a stable equilibrium. Indeed, the mechanical system would be performing small-amplitude oscillations around the applied thermal loading. These can be captured by the $\mathcal{O}(\epsilon)$ ROM as

[TABLE]

3.3 Basis interpolation and mode veering

Though the chosen subspace spanned by basis vectors at each temperature configuration would be uniquely defined, the same is not true for the basis representation of such a subspace. Each subspace would be defined by an equivalence class of orthogonal bases. There is a possibility of mode veering and rotation between bases calculated at two different temperature configurations if the representative basis for a given temperature configuration is chosen arbitrarily. This would pose a serious issue in adaptive model reduction, as proposed in this work. Indeed, the modal amplitudes in ROM equations (36)),(38) must change continuously in time. This is generally not true if the corresponding basis vectors are changing discontinuously, as would be expected when mode veering and/or rotation occurs between bases. A method to avoid this issue was proposed in ref. [10] by identifying a congruence transformation which transforms a basis such that its vectors are consistent with those in a selected reference basis. We reproduce the procedure from ref. [10] in Algorithm 1.

4 Numerical examples

We demonstrate the proposed method on a set of beam examples, with displacements in 2-dimensional $(x,z)$ plane. The beam is modeled using finite elements with a linear elastic material (see Table 1 for geometrical parameters and physical properties). Furthermore, we consider geometric nonlinearities in the model via the von-Kármán strain approximation. The thermal effects are modelled by assuming that the temperature distribution is uniform throughout the thickness of the beam. Across all examples, the beam is clamped at both ends and discretized using 60 elements.

A pulse-shaped temperature field is applied to the beam, as shown in Figure 4. It is modeled as using the following equation:

[TABLE]

where $x_{0}=x_{c}-p/2$, $x_{c}$ is the pulse-center location, $p$ is the width of the pulse, $T_{c}$ is the pulse height, $H$ is the Heaviside step function. The center of the pulse $x_{c}$ moves slowly across the beam according to the relation

[TABLE]

where $x_{0}$ is the initial location of the center of the pulse and $A$ is the amplitude of temperature-pulse-oscillation along the length of the beam.

We consider two geometries of the doubly-clamped beam model: one which is initially straight and the other with a small initial curvature. The FE model of the beam contains a temperature-dependent linear stiffness matrix. As a result, the VMs are also temperature-dependent. The variation in the linear VMs and natural frequencies of the beam for both models as the center of the temperature profile ($x_{c}$) changes location are depicted in Figure LABEL:fig:modes.

Along with the imposed thermal dynamics on the beam, we also excite the structure mechanically with a load oscillating at a much faster rate than that of the thermal distribution, the details which are described in the following examples. Across all examples, we perform numerical time integration starting with the static equilibrium at the initial thermal configuration.

Model reduction strategy: During model reduction using the proposed MMS-based technique, a basis comprising of temperature-dependent VMs (and modal derivatives in case of nonlinear problems) is interpolated as the center $x_{c}$ slowly changes according to eq. (43). In order to perform this interpolation, we use a database of $n_{d}=19$ bases, each computed around the equilibrium when the thermal pulse center is at location

[TABLE]

We ensure that all the bases in this database are consistent using Algorithm 1. As the temperature profile moves along the length of the beam, we obtain a basis relevant for the instantaneous temperature configuration by performing linear interpolation between the bases in the database. In general, such an interpolation is not guaranteed to preserve orthogonality of the bases and can even lead to ill-conditioning. This can be avoided by interpolating the bases over a Stiefel manifold, as proposed in ref. [11]. However, in the authors’ experience, this approach is not suitable for performing interpolation online, as it can be computationally intensive and a linear interpolation is found to work just as well in structural dynamics applications.

Comparison with a stacking-based approach: A simple and robust approach for reduction would constitute stacking the modes relevant to different temperature configurations in a single basis and performing orthogonalization, e.g., using the Gram-Schmidt method or singular value decomposition (SVD) (cf. Golub & van Loan [27]) . Specifically, let $\mathbf{V}^{(j)}\in\mathbb{R}^{n\times m}$ be an orthonormal basis for the $j^{\mathrm{th}}$ temperature configuration in the database of $n_{d}$ configurations. We stack these bases in a matrix

[TABLE]

to perform SVD as

[TABLE]

where $\mathbf{L}\in\mathbb{R}^{n\times n}$ is an orthonormal matrix containing the left singular vectors, $\mathbf{\Sigma}\in\mathbb{R}^{n\times n_{d}m}$ is a diagonal matrix containing the corresponding singular values, and $\mathbf{R}\in\mathbb{R}^{n_{d}m\times n_{d}m}$ is an orthonormal matrix containing the right singular vectors. The columns of $\mathbf{L}$ corresponding to the non-zero singular values in $\mathbf{\Sigma}$ form an orthonormal basis spanning the columns of $\boldsymbol{\mathcal{V}}$. Following this procedure, the matlab command $\texttt{orth}(\boldsymbol{\mathcal{V}})$ economically computes and returns the required orthonormal basis that we have used in this work.333Arguably, such a basis is conservative and might contain spurious directions due to numerical tolerances. A workaround to avoid such spurious modes would be to consider the modes associated to the non-zero singular values greater than a chosen numerical threshold, which is considerably higher than the machine zero value.

Remark.

The usual approach of plotting the singular value decay to cut-off the mode selection after appearance of a “knee” in such a plot is particularly justifiable in the context of POD, where a dominant subspace or principal components are identified from simulation snapshots or data. Such simulation snapshots carry trajectory information in contrast to the reduced bases that already represent dominant direction at different parameter configurations in the matrix $\boldsymbol{\mathcal{V}}$. In other words, one needs a basis that spans all the columns of $\boldsymbol{\mathcal{V}}$ as opposed to obtaining ”dominant directions” that approximately span the data contained in $\boldsymbol{\mathcal{V}}$.

The orthogonalized basis is expected to capture behavior across a range of temperature configurations in a more robust manner in comparison with other approaches mentioned in the Introduction. However, the size of such a basis can very quickly become unmanageable as more bases at different temperature configurations are added to the database, as shown by Figure 6. Furthermore, note that our simple examples feature a single parameter for basis description, i.e., the location $x_{c}$ of the center of a temperature pulse. This simple stacking-based approach may lead to a further increase in basis size as the parameter-dimensionality increases. To avoid such high-dimensionality for a reduced-basis, we adopt the following methods to limit the size of the reduction basis used for comparison purposes, when necessary:

Modal-POD: Here, we use $m$ modes associated to the highest singular values obtained from the SVD (46) of the matrix $\boldsymbol{\mathcal{V}}$ given in eq. (45). This results in the best constant basis of size $m$ that can be obtained from the temperature configurations in the database. Note that the size $m$ is chosen to be the same as that of the reduction basis used in the MMS-based reduction approach to facilitate a fair comparison in terms of the number of unknowns. 2. 2.

Modal: One might wish to avoid the Modal-POD type reduction in view of the above-mentioned remark. In that case, one could select a subset of temperature configurations at random from the $n_{d}$ configurations available in the database. Thus, lesser data is used to construct the matrix $\boldsymbol{\mathcal{V}}$ and a smaller basis size is obtained after orthogonalization.

4.1 Straight beam

First we consider a linear model of a straight beam, whereby the nonlinear terms arising from the von-Kármán assumption are neglected. Note, however, that the model still contains a temperature-dependent linear stiffness matrix along with a thermal load vector. In addition, the beam is mechanically excited as follows:

[TABLE]

where $\omega_{s}$ is the loading frequency, chosen to be the average of the first and second natural frequency (calculated for temperature configuration $x_{c}=L/2$), and $\mathbf{l}$ is load amplitude vector, representing a spatially uniform load in the transverse direction with a density of $10^{4}$N/m. The temperature profile traverses the beam span according to the relation (43) with $x_{0}=0.5L,A=0.3L$.

Using the stacking approach described above that involves stacking all the modes in the database together, we obtain a basis of size (rank) 93, which is quite large for full system with 177 DOFs. It turns out that this simple problem can be reduced using a constant basis of five modes according to the Modal-POD described above. We compare the solution using the MMS-based reduction approach proposed in this work using first five VMs which are adapted to the instantaneous temperature configuration. We test our approach for $\epsilon=0.01,0.001$.

The results show that the former approach using a constant basis is just as effective in capturing transverse displacements of the beam as the multiple-scales approach using an adaptive basis (cf. Figure 7a). Note, however, that the axial displacements are totally missed by this approach, as shown in Figure 7b. This is not unexpected since the first five VMs contain only low-frequency bending fields. Such a basis would thus be ineffective in capturing the membrane effects. On the other hand, the MMS-based approach systematically takes care of this issue while still using a similar basis which contains bending modes only. This is achieved by making use of the scale separation when axial components are systematically captured using a temperature-dependent equilibrium according to eq. (25). For $\epsilon=0.01$, while the MMS-based approach qualitatively captures the axial response, it is not as accurate as in capturing the transverse displacements (cf. Figure 7). However, the full system response, including the axial displacements is expected to be captured accurately for sufficiently small $\epsilon$. Indeed, this is the case for a smaller value of $\epsilon=0.001$, as shown in Figure 8.

Instead of using the proposed MMS-based reduction, the axial contributions can also be captured by appending a set of carefully selected axial modes in the original constant basis (Modal-POD). Arguably, this would be a better way to reduce the problem because even though we increase the basis size, it completely avoids the basis interpolation required in the MMS-based approach. While selecting such axial modes is relatively easy in this illustrative straight beam example, this could already be challenging if we add curvature to the beam, thereby coupling the bending and membrane DOFs. The MMS-based approach is still as effective and efficient in such situations, as we show in the next example.

4.2 Curved beam

Next, we consider a beam with the same properties as the previous example but curvature added (cf. Table 1). We consider the following mechanical excitation applied to the beam

[TABLE]

where $\omega_{c}$ is a typical loading frequency, once again taken as the average of the first and second natural frequency of the system (calculated for the temperature configuration $x_{c}=L/2$); $\mathbf{l}_{0}$ is the spatial load vector at the leading order whose shape represents a uniform pressure at the magnitude of $10^{3}$ N/m (cf. Figure 9a for its shape). In addition to this uniform load, we have a random $\mathcal{O}(\epsilon)$ perturbation that is modeled as follows. The perturbation shape $\tilde{\mathbf{l}}_{1}$ (see Figure 9b) is constructed by a linear combination of the first five VMs of the beam (computed for $x_{c}=L/2$) such that the coefficients for this combination are pseudo-random values from the standard uniform distribution in the interval $(0,1)$. The time variation of this perturbation is given by the scalar function $a(t)$ which is modeled by choosing pseudo-random values at each time instant from the uniform distribution in the interval (0,1) and thereafter filtering out the frequency contents higher than the third natural frequency of the beam. The resulting amplitude has time history as shown in  Figure 9c.

Finally, we assume that the random perturbation $\tilde{\mathbf{l}}_{1}$ models noise and has small amplitude in comparison with the leading-order forcing vector $\mathbf{l}_{0}$. As discussed in Section 2, we scale this given component of the forcing as $\tilde{\mathbf{l}}_{1}=\epsilon\mathbf{l}_{1}$, where $\mathbf{l}_{1}:=\tilde{\mathbf{l}}_{1}/\epsilon$. Our choice of loading amplitudes ensures that $\mathbf{l}_{0}$ and $\mathbf{l}_{1}$ are of similar magnitudes after this scaling, and in particular $\|\mathbf{l}_{1}\|_{2}=\|\mathbf{l}_{0}\|_{2}$. We obtain

[TABLE]

4.2.1 Linear model

We first consider the linear example of the curved beam with the mechanical forcing (49) along with a time-varying temperature distribution according to relation (43), as before. Since the spectral content of the applied mechanical forcing is contained within the first five modes of the systems, these modes constitute an ideal choice for reduction. Thus, similar to the previous straight-beam example, we use the MMS-based reduction with a 5-mode basis that adapts instantaneously to the temperature configuration of the beam; see Figure 10 for results when $\epsilon=10^{-3}$. The presence of the $\mathcal{O}(\epsilon)$ (pseudo-random) perturbation to the leading-order forcing, however, leads to more interesting results in this case. Note that the leading order ROM is uninfluenced by the $\mathcal{O}(\epsilon)$ component of the forcing and misses the corresponding fluctuating components. Nonetheless, these components are successfully captured by the $\mathcal{O}(\epsilon)$-ROM, as shown in Figure 10.

On the other hand, the stacking-based approach described earlier results in a basis of size 95, which is again not feasible for reducing a 177 DOF system. To alleviate this:

Modal ($m=15$): we apply the stacking approach to a subset of the sampled temperature configurations. The first five modes are available around the temperature configurations $x_{c}^{(j)}$ described in eq. (44). We select 3 of these temperature configurations at random ($j=2,12,18$ in this instance) and stack the corresponding modes in a matrix and perform orthogonalization to obtain a reduction basis of size 15. 2. 2.

Modal-POD ($m=5$): we obtain the $m=5$ modes associated to the highest singular values of the matrix $\boldsymbol{\mathcal{V}}$ which has $n_{d}=19$ temperature configurations (cf. eqs. (44),(45)).

A conventional reduction is then performed via Galerkin projection using these bases. The corresponding reduced solutions show a significantly worse accuracy (see Figure 10) when compared to the MMS-based reduction at leading order, as well as at $\mathcal{O}(\epsilon)$.

Upon increasing $\epsilon$, which quantifies the speed of thermal loading, we observe that the multiple scales assumption starts to deviate from the true behavior of the system, as expected. The same is demonstrated in the introduction with the simple examples in Figure 3. In particular, for $\epsilon=0.01$, we see in Figure 11 that the MMS-based reduction approach is unable to capture the Full system response, even on this linear example. The same example also demonstrates the failure of the conventional Galerkin projection-based techniques to reduce the Full system using a constant basis.

4.2.2 Nonlinear model

Finally, we consider the above-described curved beam model, but this time with the presence of geometric nonlinearities. As discussed in Section 3.2, it is well-known that a basis consisting solely of a few linear VMs is insufficient to capture the nonlinear response of a system, and the enrichment of the reduction basis with modal derivatives is an effective way to capture geometrically nonlinear response. Accordingly, we augment the database of 5-mode bases used in the linear case with the corresponding static modal derivatives. Since $k=5$ VMs result in $k(k+1)/2=15$ modal derivatives, we obtain a basis of size $m=20$ at each temperature configuration in the database.

To compare the results with another reduction approach, we again use the stacking-based approach in the following two ways:

Modal ($m=60$): Similar to the linear case, we again use a subset of randomly selected temperature configuration ($x_{c}^{(j)},j=4,7,13$) from the database. In contrast to the linear case, here we obtain a reduction basis of size 60 (since each of the three temperature configurations correspond to a basis of size 20). 2. 2.

Modal-POD ($m=20$) approach, where we use the $m=20$ modes associated to the highest singular values of the matrix $\boldsymbol{\mathcal{V}}$ which contains the first five VMs and the corresponding 15 modal derivatives at $n_{d}=19$ temperature configurations (cf. eq. (44),(45)).

Figures 12(a) and (b) show the $z$ and $x$-direction displacements at the quarter span of the beam for $\epsilon=0.001$. To provide a more global picture, we compare the error in the vector of generalized displacements relative to the full nonlinear solution $\mathbf{u}(t)$ at each time instant (cf. Figure 12c). This error is defined as

[TABLE]

where $\mathbf{u}_{red}(t)$ denotes the reduced solution obtain from the different approaches. It is interesting that despite have a basis size of 60, a constant basis reduction provides a worse accuracy than the MMS-based reduction with a basis of size 20 that adapts to the instantaneous temperature configuration.

The improvement in accuracy between the $\mathcal{O}(1)$ and the $\mathcal{O}(\epsilon)$ reduced solutions from the MMS is not very apparent from Figure 12c. For this reason, we use the following global error estimate which is uniform over time interval $[0,I]$, where $I$ is the total time over which the simulation is performed

[TABLE]

where the approximation is perform using a one-point numerical integration over the uniformly sampled (over $[0,I]$) time instants in the set $\mathcal{I}$. According to Table 2, this estimate of error shows a marginal improvement in accuracy using the $\mathcal{O}(\epsilon)$ solution over the $\mathcal{O}(1)$ solution, as expected.

5 Discussion

The numerical examples above serve as proof of concept for effective reduction using the MMS-based approach, proposed in this work. Although these examples are considerably simpler than the targeted applications of model reduction techniques, we do not envision issues with the scalability of this systematic approach. This is because the adaptive interpolation ensures that the number of reduced variables remains uniform as the temperature is dynamically changing during a simulation. Nonetheless, we discuss some general sources of computational bottleneck and possible ways to address them in the following.

Arbitrarily varying temperature fields.

For simplicity, we have only considered externally-imposed temperature variation on the beam in this work. However, this variation in structural systems is usually obtained by solving the heat equation (5), which is given in its linearized form by the system

[TABLE]

where the heat source $\mathbf{h}$ is externally applied, and the temperature $\mathbf{T}$ is indeed obtained as an output. For the numerical examples treated in this work, a single parameter ($x_{c}$) was sufficient to model the variation in temperature distribution over the beam. Thus, the database of reduction bases (used in adaptive basis selection) was obtained for a range of single-parameter values (44). General temperatures fields, however, require $n_{T}$ (number of DOFs in eq. (5)) parameters for description, which could be potentially large and would require a huge database of reduction bases to be computed. Thus, the memory requirements for the adaptive basis interpolation could be a potential bottleneck for adaptive basis selection in case of multi-dimensional parameters. One simple approach to reduce the number of temperature field parameters is as follows.

The temperature distribution of the structure can be spectrally decomposed into a basis of thermal modes obtained from the eigenvalue analysis of the linear operators in $\eqref{eq:lin_heat}$ as

[TABLE]

where $\mathbf{\boldsymbol{\psi}}_{j}$ is the $j^{\mathrm{th}}$ eigenvector related to the thermal problem (5) and $\lambda_{j}$ is the corresponding eigenvalue. The temperature configuration can then be suitably approximated using a few thermal modes as

[TABLE]

where $\boldsymbol{\Psi}\in\mathbb{R}^{n_{T}\times m_{T}}$ is a matrix containing the thermal modes of the structure and the slowly varying amplitudes $\mathbf{q}_{T}(\tau)\in\mathbb{R}^{m_{T}}$ of the thermal modes determine the temperature distribution across the structure at a given time. It is then easy to see that the reduction basis for structural dynamics can be parameterized by the thermal mode amplitudes $\mathbf{q}_{T}$, thereby reducing the number of required parameters to $m_{T}\ll n_{T}$.

Computational time and hyper-reduction.

In the case of linear systems such as eq. (13) with time-invariant coefficient matrices, the reduced operators can be precomputed at once when a constant reduction basis is used. However, with time-dependence in the coefficient matrices such as in eq. (19)), the reduced operators need to be updated online at each time step during time integration. This operation poses computational bottlenecks in reduction of such systems, and a reduction in the number of unknowns does not lead to a tantamount reduction in computational time. Note that this issue affects reduction techniques using a constant basis, as well as the MMS-based reduction technique proposed here. A possible solution to this issue is to precompute the reduced operators $\mathbf{M}_{\boldsymbol{\Phi}}$, $\mathbf{C}_{\boldsymbol{\Phi}}$ and $\mathbf{K}_{\boldsymbol{\Phi}}$ in eq. (26) along with the reduction bases at each temperature configuration in the database, during the offline stage. Thereafter, this reduced information can be directly interpolated online (cf. ref. [10]) during time integration. This completely avoids the interpolation of reduction bases and projection at each time step at the cost of precomputation of reduced information.

For nonlinear reduced-order models, the evaluation of reduced nonlinear operators by projection onto the reduction basis can also be a computational bottleneck [21] during online time integration. Indeed, due to these extra operations, the computation of a reduced solution becomes possibly more expensive than that for a full system solution [22]. Hyper-reduction techniques provide reprieve in such cases by a fast approximation of the reduced nonlinear operators using training and selection procedures. In the finite-element context, the energy conserving sampling and weighing method (ECSW) [23] is found to be effective. Based on offline training, this method provides a single reduced finite-element mesh for a given set of reduction bases and parameter values, as also discussed in ref. [24]. Furthermore, the hyper-reduction of geometrically nonlinear structural systems using ECSW with VMs and modal derivatives basis has already been affirmed in previous works [25, 26]. These features make it possible to train a ROM across a range of temperature configurations and obtain a reduced mesh so as to accelerate computation for adaptive basis reduction using the MMS-based approach proposed in this work.

6 Conclusion

In this work, we focused on model reduction for temperature-dependent structural dynamics equations with time varying temperature fields. We discussed that a slow variation in the structural temperature is not only physically relevant but also essential to justify model reduction using an adaptive reduction basis (cf. Section 1). Subsequently, we proposed the systematic use of the method of multiple scales to exploit this slow temperature dependence in the structural dynamics equations. We consistently reduced the equations of motion using a temperature-dependent basis that slowly adapts to the instantaneous temperature configuration of the structure (cf. Section 3).

We treated numerical examples in the linear case using an adaptive basis containing vibration modes and in the geometrically nonlinear setting, the basis was enriched using the corresponding modal derivatives. In the process, we also concluded (cf. Section 4) that the use constant reduction basis obtained by systematically combining reduction bases obtained from different temperature configurations is not a scalable strategy, as this results in a potentially large number of reduced variables. On the other hand, the MMS-based reduction approach turned out to be very efficient in terms of number of variables and showed consistently better accuracy.

While we addressed general computational bottlenecks in the reduction procedure and discussed possible strategies to mitigate them (cf. Section 5 ), the main focus of this work was efficient dimensionality reduction of the governing equations. Implementation of these strategies in line with the discussion in Section 5 is essential to obtain the necessary and sufficient computational speed. This forms part of our future efforts. In combination with these measures, we expect the MMS-based adaptive reduction would enable us to reduce realistic structures under the influence of thermal environments in a robust and efficient manner.

Acknowledgments.

The authors would like to thank the anonymous reviewers of this work for their very instrumental feedback. The support of the Air Force Office of Scientific Research, Air Force Material Command, USAF under Award No.FA9550-16-1-0096 is acknowledged.