We present an efficient method to significantly reduce the offline cost associated with the construction of training sets for hyper-reduction of geometrically nonlinear, finite element (FE)-discretized structural dynamics problems. The reduced-order model is obtained by projecting the governing equation onto a basis formed by vibration modes (VMs) and corresponding modal derivatives (MDs), thus avoiding cumbersome manual selection of high-frequency modes to represent nonlinear coupling effects. Cost-effective hyper-reduction is then achieved by lifting inexpensive linear modal transient analysis to a quadratic manifold (QM), constructed with dominant modes and related MDs. The training forces are then computed from the thus-obtained representative displacement sets. In this manner, the need of full simulations required by traditional, proper orthogonal decomposition (POD)-based projection and training is completely avoided. In addition to significantly reducing the offline cost, this technique selects a smaller hyper-reduced mesh as compared to POD-based training and therefore leads to larger online speedups, as well. The proposed method constitutes a solid alternative to direct methods for the construction of the reduced-order model, which suffer from either high intrusiveness into the FE code or expensive offline nonlinear evaluations for the determination of the nonlinear coefficients.

## Introduction

The use of finite element (FE)-based structural dynamic models to simulate structural response has become a standard practice in many fields of application, such as mechanical engineering, aerospace engineering, and bio-mechanics. These models often tend to possess large number of degrees-of-freedom (DOFs) in order to account for extremely detailed geometric features and material distribution. However, routine simulations to explore different load scenarios, geometric layouts and material choice for such large models carry prohibitive computational costs. In this context, reduced-order models (ROMs) offer a reprieve from such extreme computational costs and allow effective design and optimization activities.

In a broad sense, ROMs are low-dimensional counterparts of the original model, often referred to as high-fidelity model (HFM). Classically, this reduction is achieved through a linear projection of the full system of equations onto a reduction basis, which has been precomputed *offline*. This basis spans a low-dimensional subspace suitable for capturing the HFM solution. As a result, the ROM is characterized by only a few DOFs. The construction cost of the reduced operators in the projected equations, however, scales with the size of the HFM and not with that of the ROM. This constitutes a serious bottleneck for achieving significant speed-ups, as the assembly of the nonlinear reduced operators dominates the cost of the time integration. Hence, such ROMs are effective only if these reduced operators can be precomputed offline, which is not the case for nonlinear systems in general.

More often than not, nonlinear modeling is essential for design and analysis of realistic structures, even in the preliminary stages. Thin-walled structural components, for instance, are typically employed in the industry when high stiffness-to-weight and strength-to-weight ratios must be achieved. Among other nonlinear effects, the geometric nonlinearities are particularly important in modeling their behavior, including peculiar effects like bending-stretching coupling; buckling; snap-through; mode jumping, etc., due to finite rotations [1]. This has given rise to a pressing need for fast and effective reduction of large nonlinear dynamical systems. To this effect, various techniques have been developed over the recent years, which have made the *online* computation of the reduced nonlinear operators in ROMs much cheaper [2,3–6]. These are commonly referred to as *hyper-reduction* techniques.

Generally, hyper-reduction techniques alleviate the computational costs of the nonlinear terms by optimally selecting a small set of nodes (or elements) in the mesh over which the nonlinearity is evaluated. The nonlinearity is then cheaply interpolated over the rest of the mesh. This selection process is performed with the help of *training* vectors, which are usually obtained from the solution of the HFM. Such *full*-solution vectors are often also used for the construction of the reduction basis in projection-based ROMs, for example, using the proper orthogonal decomposition (POD) [7–9]. The use of these full-solution *snapshots* for training and reduction purposes is a computationally expensive affair, which can be unaffordable, especially at the preliminary design stage of structures.

Many techniques enable the construction of a ROM or a reduction basis for nonlinear problems without the need of a full solution (cf., Refs. [10–12]). Furthermore, for a certain class of problems—characterized by mild geometric nonlinearities—the use of vibration modes (VMs) and modal derivatives (MDs) has been shown to be effective for the construction of a modal-based reduction basis [13–15]. Indeed, the robustness and effectiveness of full-simulation vectors in general nonlinear reduction scenarios (e.g., using POD) cannot be discounted. Nonetheless, such modal-based reduction techniques, albeit being suboptimal, find use when the available time and computational resources do not allow for creating a database of full simulation runs. By avoiding full solutions, these techniques go a long way in reducing offline costs incurred during construction of an effective reduction basis.

The generation of training vectors for hyper-reduction of nonlinear terms in the ROM, however, still leads to tremendous offline costs in the form of full-simulations. To this effect we propose a modal-based training-set generation technique, which completely avoids the HFM simulations, thereby, making the reduction procedure truly *simulation-free.* As an alternative to hyper-reduction of nonlinearity, techniques which precompute the reduced nonlinear operators in the form of polynomial tensors up to cubic degree are also common in literature [16–20]. These techniques, however, are usually based upon an ad-hoc selection of VMs and are limited in their scope (see Sec. 7 for a discussion).

The modal derivatives have been classically used to form a reduction basis along with VMs in order to systematically capture the geometrically nonlinear behavior. Recently, these were also used for reduction via a quadratic manifold (QM) [13], where a linear subspace, formed by a truncated set of VMs, captures the linearized dynamics near the equilibrium and the corresponding MDs provide the necessary nonlinear (quadratic) extension to this subspace. In this work, we use this notion of a QM to generate meaningful training vectors from a linear modal superposition of the underlying linearized system. This unified approach builds upon the linear modal signature, which is rather cheaply available and essential for the analysis of any structural system. We have tested this approach on a simple, illustrative example as well as on a realistic structural model. To demonstrate effectiveness of the approach, we have compared the offline costs involved in reduction and hyper-reduction of these examples with that of other established techniques. This is done by taking into account all the effort needed to construct such a simulation-free hyper-reduced ROM.

This paper is organized as follows: in Sec. 2, we start by reviewing the concepts of projection-based model reduction, which results in the reduction of the dimensionality of the HFM. The concept of hyper-reduction is reviewed in Sec. 3, where we propose the use of a stability-preserving, FE-based hyper-reduction technique, known as energy conserving sampling and weighing (ECSW) [2]. The main idea of modal-based training-set generation using the QM is presented in Sec. 5. The results for the numerical examples, along with comparisons with other techniques, are presented in Sec. 6. Section 7 discusses the state-of-the-art alternatives in the spirit of simulation-free reduction. Finally, the conclusions are given in Sec. 8.

## Dimensionality Reduction

**u**of the structure, and $g(t)\u2208\mathbb{R}n$ is the time dependent external load vector. In our case, the nonlinear term $f(u)$ models the effect of geometric nonlinearities, arising from large deflections and rotations. Though the techniques presented in this work are general, von Kármán kinematics has been used to model geometric nonlinearities for shell-based structures. As discussed in Sec. 1, the system Eq. (1) is referred to as the HFM. The response of the HFM can be extremely time consuming to compute if the dimension

*n*of the system is large. The classical notion of model reduction aims to reduce this dimensionality by introducing a linear mapping on to a suitable low-dimensional subspace $V$ as

**V**is known as the reduction basis since its columns form a basis for $V$. The ROM is then obtained using Galerkin projection as

*m*. This is, however, not the case for the computation of the nonlinear term $f\u0303(q(t))$. For FE-based applications, this evaluation is usually carried out online in the following manner:

where $fe(ue)\u2208\mathbb{R}Ne$ is the contribution of the element *e* toward the vector $fnl(u)$ (*N _{e}* being the number of DOFs for the element

*e*), $Ve$ is the restriction of

**V**to the rows indexed by the DOFs corresponding to

*e*, and

*n*is the total number of elements in the structure. Since the reduced nonlinear term $f\u0303(q)$ is evaluated in the space of full variables, the computational cost associated with its evaluation does not scale with

_{e}*m*alone. Indeed, Eq. (3) shows that this cost scales linearly with the number of elements in the structure and can, hence, be high for large systems. Thus, despite the reduction in dimensionality achieved in Eq. (2), the evaluation of the reduced nonlinear term $f\u0303(q)$ emerges as a new bottleneck for the fast prediction of system response using the ROM. Hyper-reduction techniques help mitigate these high computational costs by approximation of the reduced nonlinear term in a computationally affordable manner.

## Hyper-Reduction

Different hyper-reduction techniques exist in literature, all aiming at reducing the evaluation of subspace-projected nonlinear terms. The empirical interpolation method [6] approximates a nonaffine function using a set of basis functions over a continuous, bounded domain, and hence finds applications for reduced basis representations of parameterized partial differential equations. The discrete empirical interpolation method (DEIM) [3] achieves this on a spatially discretized domain. For FE applications, the unassembled DEIM or UDEIM [4] was shown to be more efficient, where the unassembled element-level nonlinear force vector is used in the DEIM algorithm to return a small set of elements (instead of nodes) over which the nonlinearity is evaluated. For conservative/Hamiltonian systems, however, the DEIM (or UDEIM) lead to a loss of numerical stability during time integration. Though a symmetrized version of DEIM has been proposed recently [5] to avoid this issue, its potential for FE-based applications remains unexplored. To this end, the ECSW hyper-reduction method [2], which directly approximates the reduced nonlinear term $f\u0303(q)$ without the loss of numerical stability [21], is remarkable.

*E*of the structure ($|E|\u226ane$) to cheaply approximate $f\u0303(q)$ as (cf., Eq. (3))

where $\xi e\u2208\mathbb{R}+$ is a positive weight attached to each element $e\u2208E$, which is empirically chosen to ensure a good approximation of the summation in Eq. (3). Clearly, if $|E|\u226ane$, then the evaluation of the approximation in Eq. (4) would come at a fraction of the original computational cost associated with Eq. (3). In doing so, ECSW approximates the virtual work done by the internal forces on the set of vectors in the basis **V**. As a consequence, the ECSW preserves the structure and the stability properties of the underlying full model (cf., Ref. [21]).

*n*training vectors in the set with $u(i)$ representing the

_{t}*i*th vector, then corresponding reduced unknowns $q(i)$ can be easily calculated using least squares as

**G**as follows:

*E*used in Eq. (4) as

An optimally sparse solution to $(P1)$ is NP-hard to obtain. However, a greedy-approach-based algorithm [22], which finds a suboptimal solution, has been found to deliver an effective reduced mesh *E* [2]. For the sake of completeness, we have included its pseudo-code in Algorithm 1, where $\zeta E$ and $GE$ denote, respectively, the restriction of $\zeta \u2208\mathbb{R}ne$ and column-wise restriction of $G$ to the elements in the *active* subset *E*. The set *Z* is the disjoint *inactive* subset, which contains the zero entry-indices of $\xi $ and $\zeta $.

Input:$G,b,\tau $ |

Output:$\xi \u2208\mathbb{R}ne$ sparse, $E\u2282{1,\u2026,ne}$ |

1: $E\u2190\u2205,Z\u2190{1,\u2026,ne},\xi \u21900\u2208\mathbb{R}ne$ |

2: while$\Vert G\xi \u2212b\Vert 2>\tau \Vert b\Vert 2$do |

3: $\mu \u2190GT(b\u2212G\xi )$ |

4: $[|\nu |,e]=max(\mu )$ ▹ $max$: returns maximum value in a vector followed by its location (index) |

5: $E\u2190E\u222a{e},Z\u2190Z\{e}$ |

6: while true do |

7: $\zeta E\u2190GE\u2020b$ ▹$\u2020$ represents pseudo-inverse |

8: $\zeta Z\u21900$ |

9: if$\zeta E>0$then |

10: $\xi \u2190\zeta $ |

11: break |

12: $\eta =mink\u2208E\xi k/(\xi k\u2212\zeta k)$ |

13: $\xi \u2190\xi +\eta (\zeta \u2212\xi )$ |

14: $Z\u2190{i|\xi i=0}$ |

15: $E\u2190{1,\u2026,ne}\Z$ |

Input:$G,b,\tau $ |

Output:$\xi \u2208\mathbb{R}ne$ sparse, $E\u2282{1,\u2026,ne}$ |

1: $E\u2190\u2205,Z\u2190{1,\u2026,ne},\xi \u21900\u2208\mathbb{R}ne$ |

2: while$\Vert G\xi \u2212b\Vert 2>\tau \Vert b\Vert 2$do |

3: $\mu \u2190GT(b\u2212G\xi )$ |

4: $[|\nu |,e]=max(\mu )$ ▹ $max$: returns maximum value in a vector followed by its location (index) |

5: $E\u2190E\u222a{e},Z\u2190Z\{e}$ |

6: while true do |

7: $\zeta E\u2190GE\u2020b$ ▹$\u2020$ represents pseudo-inverse |

8: $\zeta Z\u21900$ |

9: if$\zeta E>0$then |

10: $\xi \u2190\zeta $ |

11: break |

12: $\eta =mink\u2208E\xi k/(\xi k\u2212\zeta k)$ |

13: $\xi \u2190\xi +\eta (\zeta \u2212\xi )$ |

14: $Z\u2190{i|\xi i=0}$ |

15: $E\u2190{1,\u2026,ne}\Z$ |

## Classical (Hyper) Reduction

*best*spans the vectors in this ensemble can be obtained by the solution to the following minimization problem:

**V**are the left singular vectors of

**U**, obtained by the singular value decomposition (SVD) of

**U**as

^{2}

**U**is factorized into unitary matrices $L=[l1,l2,\u2026,ln]\u2208\mathbb{R}n\xd7n$ (containing the left singular vectors) and $R\u2208\mathbb{R}nt\xd7nt$ (containing the right singular vectors); and the diagonal (rectangular) matrix $S\u2208\mathbb{R}n\xd7nt$ (containing corresponding singular values on the diagonal). These singular values represent the relative importance of corresponding vectors of

**L**in forming the basis

**V**. If the singular values (and the corresponding singular vectors) are arranged in a descending order $S11\u2265S22\u2265\cdots \u22650$, it can be shown that

Thus, the left singular vectors $li$ corresponding to the highest singular values are the most relevant for constructing a reduction basis using POD. Since an SVD can be performed over any solution snapshots of any general nonlinear problem, POD is seen as a versatile method. It should be noted, however, that such a reduction basis is optimal only for capturing a solution which is characteristic of the applied loading, and the procedure needs to be repeated (or a database of full solutions need to be created in the first place) to take other types of loading into account.

Once a reduction basis is obtained, the same ensemble of full solution vectors can be used as training vectors for hyper-reduction (e.g., using ECSW) to accelerate the computation of nonlinearity during the online stage. The steps involved in a POD-based ECSW, using full-solution snapshots for training, have been outlined in Algorithm 2.

Input:$u(i),i\u2208{1,\u2026,nt}$ snapshots from full solution |

Output:$\xi \u2208\mathbb{R}ne$ sparse, $E\u2282{1,\u2026,ne}$ for ECSW-POD |

1: $U\u2190[u(1),\u2026,u(nt)]$ |

2: $[V]\u2190POD(U,M)$ |

3: $[G,b]\u2190CONSTRUCTGb(U,V)$^{a} |

4: $[\xi ,E]\u2190SNNLS(G,b)$^{b} |

Input:$u(i),i\u2208{1,\u2026,nt}$ snapshots from full solution |

Output:$\xi \u2208\mathbb{R}ne$ sparse, $E\u2282{1,\u2026,ne}$ for ECSW-POD |

1: $U\u2190[u(1),\u2026,u(nt)]$ |

2: $[V]\u2190POD(U,M)$ |

3: $[G,b]\u2190CONSTRUCTGb(U,V)$^{a} |

4: $[\xi ,E]\u2190SNNLS(G,b)$^{b} |

CONSTRUCTGb constructs matrix **G** and vector **b** using the training sets in *U* and reduction basis *V*, according to Eq. (5).

SNNLS performs sparse non-negative least squares solution according to Algorithm 1.

## Simulation-Free Reduction

In this section, we propose a systematic procedure for obtaining a reduction basis, as well as the necessary training vectors for hyper-reduction, while completely avoiding full simulation runs.

### Modal-Based Reduction Basis.

**V**as

**V**is obtained after an orthogonalization of the matrix $\Psi $ in Eq. (8), $orth$ represents any routine to orthogonalize any general matrix such that the result is a orthonormal matrix with a full column rank (e.g., the Gram–Schmidt orthogonalization). The matrix $\Psi $ contains the relevant VMs and MDs, $\varphi 1,\u2026,\varphi M$ represent a truncated set of VMs obtained from the solution of the undamped eigenvalue problem

^{3}

*η*of the VM $\varphi j$, after

_{j}**K**is replaced by the tangent stiffness matrix $(\u2202f(u)/\u2202u)$. The singular system Eq. (11) requires an additional normalization condition to solve for $\theta ij$. We use the mass normalization of the VMs, which results in the following normalization constraint:

where $\theta ijstatic$ were first introduced in Ref. [15]), and termed as static MDs (SMDs) in Ref. [13].

*M*, VMs are included. This could potentially increase the reduction basis significantly as

*M*increases. Techniques to effectively select MDs (cf., Refs. [13] and [26]) relevant for a given load scenario are useful in such situations. The SMDs have been used in this work for construction of the QM. We have used the following simple criterion (called maximum modal interaction in Ref. [13]) to assign weight

*W*to an SMD $\theta ijstatic$:

_{ij}where $\eta i(t)$ and $\eta j(t)$ are the modal amplitudes for VMs $\varphi i$ and $\varphi j$, respectively, obtained during linear modal superposition for the applied loading. Along with the *M* VMs, a total of *M* most significant SMDs have been used to construct the reduction basis. Thus, after orthogonalization Eq. (9), an effective reduction basis with size $m\u22642M$ can be obtained with minimal offline costs using VMs and MDs.

### Quadratic Manifold-Based Training Set Generation.

Starting with the idea of modal superposition, we propose a QM-based approach [13] to obtain relevant training vectors without the need of these full solution vectors. The linear modal superposition using a few significant VMs is a well-established technique to obtain the reduced solution of a linear dynamical system. However, when the nonlinearities become significant, the modal basis can be augmented with MDs to effectively capture the response.

where $M\u2282{1,\u2026,n}$ is the set of indices corresponding to the most significant modes in modal superposition 15.

**z**using a quadratic mapping as

^{4}An arbitrary vector of modal coordinates is denoted by $z\u2208\mathbb{R}M$. As discussed in Ref. [13], the tangent space of the QM at equilibrium ($u=0$) is the modal subspace represented by $\Phi $, which is corrected using the MD information to account for nonlinear behavior upon departure from the equilibrium. Furthermore, the MDs capture the inherent bending/torsion-stretching coupling arising from geometric nonlinearities. These components arise naturally as second components of the QM. This provides a straight-forward method to generate physically meaningful training vectors using the snapshots of the modal amplitudes $\eta (t)$ in Eq. (17), as

where $ti\u2208T$ (with $T$ being the simulation time-span) are the time instants at which the modal snapshots are captured; *n _{t}* is the number of such snapshots chosen for training. Thus, the idea involves the lifting of the linear modal snapshot to the QM to account for the nonlinearity, as illustrated in Fig. 1.

This simple criterion implicitly assumes that the linear behavior is correct up to a first-order and the essential nonlinear bending/torsion-stretching coupling effects (captured using the QM) are of second order. This does not necessarily imply that the nonlinear forces are small. In fact, the range of deflections we are interested in makes the linear and nonlinear forces comparable. This is particularly in accordance with the von Kármán kinematics. By construction, however, this approach is not restricted to von Kármán nonlinearities. The main advantage of this method lies in the fact that the linear modal solution is available for any given system at practically no cost, resulting in physically relevant training vectors without the need of expensive high-fidelity simulation. The training vectors thus obtained could, in principle, be used for training any hyper-reduction technique. However, for our application and testing purposes, we use the ECSW-based hyper-reduction technique. We have shown the essential steps for a simulation-free hyper-reduction using the ECSW in Algorithm 3.

Input: HFM model as described by Eq. (1), tolerance for ECSW-based training τ |

Output:$\xi \u2208\mathbb{R}ne$ sparse, $E\u2282{1,\u2026,ne}$ for ECSW-Modal |

1: $[\Phi ,\eta (t)]\u2190ModalSuperposition(model,g(t),T)$^{a} |

2: $M\u2190size(\Phi ,2)$^{b} |

Calculation of (S)MDs to construct$\Omega $in Eq. (17) |

3: $\Omega \u2190$ all-empty 3-tensor$\u2208\mathbb{R}n\xd7M\xd7M$, $\Psi \u2190\Phi $ |

4: for$j\u21901\u2009to\u2009m$do |

5: Compute $(\u22022f(u)/\u2202u\u2202u)|u=0\varphi j$ |

6: for$i\u21901\u2009to\u2009M$do |

7: $\omega \u2190\theta ij$ or $\theta ijstatic$^{c} |

8: $\Omega (:,i,j)\u2190\omega $, $\Psi \u2190[\Psi ,\omega ]$ ▹ matlab notation |

Perform selection of MDs in$\Psi $, if needed |

9: if MD selection = true then |

10: $\Psi \u2190MDSelection(\Psi ,\eta (t))$^{d} |

11: $[V]\u2190orth(\Psi )$^{e} |

Generate training vectors |

12: Build function $QM$: $\Gamma (z)\u2190QM(z)$^{f} |

13: for$i\u2208{1,2,\u2026,nt}$do |

14: $u(i)=QM(\eta (ti))$^{g} |

15: $U\u2190[u(1),\u2026,u(nt)]$ |

Perform training |

16: $[G,b]\u2190CONSTRUCTGb(U,V)$^{h} |

17: $[\xi ,E]\u2190SNNLS(G,b)$^{i} |

Input: HFM model as described by Eq. (1), tolerance for ECSW-based training τ |

Output:$\xi \u2208\mathbb{R}ne$ sparse, $E\u2282{1,\u2026,ne}$ for ECSW-Modal |

1: $[\Phi ,\eta (t)]\u2190ModalSuperposition(model,g(t),T)$^{a} |

2: $M\u2190size(\Phi ,2)$^{b} |

Calculation of (S)MDs to construct$\Omega $in Eq. (17) |

3: $\Omega \u2190$ all-empty 3-tensor$\u2208\mathbb{R}n\xd7M\xd7M$, $\Psi \u2190\Phi $ |

4: for$j\u21901\u2009to\u2009m$do |

5: Compute $(\u22022f(u)/\u2202u\u2202u)|u=0\varphi j$ |

6: for$i\u21901\u2009to\u2009M$do |

7: $\omega \u2190\theta ij$ or $\theta ijstatic$^{c} |

8: $\Omega (:,i,j)\u2190\omega $, $\Psi \u2190[\Psi ,\omega ]$ ▹ matlab notation |

Perform selection of MDs in$\Psi $, if needed |

9: if MD selection = true then |

10: $\Psi \u2190MDSelection(\Psi ,\eta (t))$^{d} |

11: $[V]\u2190orth(\Psi )$^{e} |

Generate training vectors |

12: Build function $QM$: $\Gamma (z)\u2190QM(z)$^{f} |

13: for$i\u2208{1,2,\u2026,nt}$do |

14: $u(i)=QM(\eta (ti))$^{g} |

15: $U\u2190[u(1),\u2026,u(nt)]$ |

Perform training |

16: $[G,b]\u2190CONSTRUCTGb(U,V)$^{h} |

17: $[\xi ,E]\u2190SNNLS(G,b)$^{i} |

Modal superposition performs modal superposition on the model using Eq. (15) and return the truncated set of $M=|M|$ significant VMs in $\Phi $ and their corresponding amplitudes $\eta (t)$ according to Eq. (16).

size$(\Phi ,2)$ returns the number of columns in the matrix $\Phi $ (same as matlab$size$ function).

Use MD selection techniques (cf., Eq. (13), Refs. [13] and [26]) to retain few relevant MDs in $\Psi $.

orth computes a full column rank orthonormal matrix **V** from $\Psi $ shown in (same as matlab$orth$ function).

Using VMs and MDs, build function $QM$ with argument **z**, which returns $\Gamma $ according to Eq. (17).

QM produces training vectors from modal snapshots $\eta (ti)$ using Eq. (17).

CONSTRUCTGb constructs matrix **G** and vector **b** using the training sets in **U** and reduction basis **V**, according to Eq. (5).

SNNLS performs sparse non-negative least squares solution according to Algorithm 1.

## Numerical Results

### Setup.

where $u(t)\u2208\mathbb{R}n$ is the vector of generalized displacements at the time *t*, obtained from the HFM solution, $u\u0303(t)\u2208\mathbb{R}n$ is the solution based on the (hyper) reduced model, and *S* is the set of time instants at which the error is recorded. The mass matrix **M** provides a relevant normalization for the generalized displacements, which could be a combination of physical displacements and rotations, as is the case in the shell models shown here.

*T*

_{full}and $Tsim\u22c6$ represent the computation time taken during the time integration of $full$ and (hyper)reduced solution respectively. The superscript $\u22c6$ denotes the reduction technique being used. The speed-up defined in this manner takes only the online costs into account. For a fairer comparison, we take offline costs into account by defining an $effective$ speed-up as

where $Toff\u22c6$ represents the computational time spent $offline$ for setting up a (hyper-)reduced model; $con,coff\u2208[0,1]$ represent the relative weights given to the online and offline costs, respectively, such that $con+coff=1$. The HFM simulation is assumed to carry zero offline costs and, thus, $Seff=S=1$ for a full HFM solution. A high value of *S*_{eff} would make a given (hyper-) reduction method favorable, with *S*_{eff} < 1 implying that the corresponding reduction technique is effectively more expensive than a full solution run.

One method of choosing *c*_{on} and *c*_{off} may be based on the number of simulations *n*_{sim} performed using a ROM as $con=(nsim/nsim+1),\u2009coff=(1/nsim+1)$. In this work, we have given an equal weightage to the online and offline costs, i.e., $con=coff=0.5$_{,} which leads to a rather conservative estimate of the effective speed-up.

During hyper-reduction using ECSW, the convergence tolerance *τ* for finding a sparse solution to the NNLS problem Eq. (6), as required in Algorithm 1, is chosen to be 0.01. This is within the practically recommended range as reported in Ref. [2]. The results from the following reduction techniques are compared in both examples considered in this work (model-I and model-II):

*Simulation-free reduction* (no hyper-reduction): These techniques involve the simulation of the Galerkin ROM Eq. (2), where the reduction basis **V** is obtained without the use of full simulation snapshots. We consider the following two subtypes:

**L1**: Here, a few significant VMs (say*M*in number) of the structure are selected based on linear modal superposition for the given load. These VMs, along with the corresponding SMDs (which would be $(M(M+1)/2)$ in number), are used to construct**V**. Thus, the size of the basis would be $m=(M(M+3)/2)$.**L2**: A similar basis as that in L1 is used here, except for the fact that*M*of most significant SMDs are selected according to the maximum modal interaction criterion [13], originally proposed in Ref. [26] (cf., Eq. (13)). Thus, the basis contains*M*VMs as well as SMDs, making the basis size $m=2M$.*Classical reduction*(no hyper-reduction): Here, the nonlinear system using is reduced using a POD basis.**POD-1**: The reduction basis of the same size as in L1 is constructed through the SVD of the full-solution-snapshots-matrix, as explained in Sec. 4. Thus, $m=(M(M+3)/2)$ left singular vectors with the highest singular values are included in**V**.**POD-2**: Same as POD-1, except the basis size is matched with that in L2. Thus, $m=2M$ most significant singular vectors are used.*Classical hyper-reduction*.**ECSW-POD**: The ECSW method is used with the same basis in POD-1, to hyper-reduce the nonlinearity in the ROM. The HFM simulation snapshots are used in the offline stage for training to obtain the reduced mesh.*Simulation-free hyper-reduction*: The proposed QM-based training set generation is used to reduce offline costs during hyper-reduction. Thus, a truly simulation-free reduction is demonstrated with the following approaches:**Modal-ECSW-L1**: Here, the previously described ROM in L1 is hyper-reduced using ECSW with the help of QM-based training vectors as described in Algorithm 3.**Modal-ECSW-L2**: Same as Modal-ECSW-L1, except the ROM from L2 is hyper-reduced, instead of that from L1.

### Shallow-Arch Structure.

where **l** is a constant load vector corresponding to a uniform pressure distribution of 1 N/mm^{2}. Here, *p*(*t*)—termed as the dynamic load factor—determines the time-dependency of the external load. The results are shown for a quasi-periodic choice for *p*(*t*) as given by Eq. (19), where *ω* is a typical loading frequency chosen to be the second eigenfrequency of the linearized system. The amplitude of loading is kept large enough to trigger nonlinear behavior, as shown in Fig. 3.

A linear modal analysis shows that the first and the fifth modes of the linearized model are sufficient for capturing the linear dynamics using modal superposition. We consider these modes along with the corresponding SMDs for reduction in L1.^{5} Thus, number of modes *M* = 2 and the reduction basis size $m=(M(M+3)/2)=5$ for the technique L1. The L2 approach uses *M* = 2 most significant SMDs in the reduction basis along with the 2 VMs, introducing $m=2M=4$ reduced unknowns. POD-1 and POD-2 techniques involve the use of full-solution vectors for generating reduction bases of size same as L1 and L2, respectively. As described in Sec. 6.1, the ROM POD-1 is hyper-reduced using the classical ECSW-POD method, while the ROMs L1 and L2 are hyper-reduced using the simulation-free hyper-reduction techniques Modal-ECSW-L1 and Modal-ECSW-L2, respectively. The results for hyper-reduction, compared with their nonhyper-reduced counterparts, are shown in Figs. 4 and 5. The relative error estimates, and (effective) speed-ups for different reduction techniques are reported in Table 1.

Type | Reduction technique | Basis size (m) | # Elements | GRE (%)_{M} | $S\u22c6$ | $Seff\u22c6$ |
---|---|---|---|---|---|---|

Classical reduction | POD-1 | 5 | 400 | 2.86 | 2.85 | 0.74 |

Classical hyper-reduction | ECSW-POD | 5 | 19 | 3.04 | 50.4 | 0.94 |

Simulation-free reduction | L1 | 5 | 400 | 3.36 | 2.88 | 2.85 |

Simulation-free hyper-reduction | Modal-ECSW-L1 | 5 | 11 | 4.95 | 79.9 | 15.9 |

Type | Reduction technique | Basis size (m) | # Elements | GRE (%)_{M} | $S\u22c6$ | $Seff\u22c6$ |
---|---|---|---|---|---|---|

Classical reduction | POD-1 | 5 | 400 | 2.86 | 2.85 | 0.74 |

Classical hyper-reduction | ECSW-POD | 5 | 19 | 3.04 | 50.4 | 0.94 |

Simulation-free reduction | L1 | 5 | 400 | 3.36 | 2.88 | 2.85 |

Simulation-free hyper-reduction | Modal-ECSW-L1 | 5 | 11 | 4.95 | 79.9 | 15.9 |

We report the comparison of accuracy and computational speed gains for different methods in Table 1. The computations were performed for POD-1, POD-2, L1, L2, Modal-ECSW-L1, and Modal-ECSW-L2. Out of these, we report here only the details for POD-1, L1 and Modal-ECSW-L1. For this small problem, the left out cases give practically the same as the ones reported in Table 1, which suffice in demonstrating the main point. The other cases are relevant for the wing example discussed in Sec. 6.3.

Here, the POD-1 and POD-2 technique show an effective speed-up $Seff\u22c6<1$, as expected from the very definition of $Seff\u22c6$ (cf., Eq. (18)). Furthermore, the hyper-reduction of the POD-based ROMs (ECSW-POD) increases the value the effective speed-up, but these values are still poor. Indeed, this is due to the use of full solution vectors, whether for construction of the reduction basis **V**, or for generation of training sets for hyper-reduction. The simulation-free reduction techniques (L1 and L2), on the other hand, show $Seff\u22c6>1$. The simulation-free hyper-reduction techniques (Modal-ECSW-L1 and Modal-ECSW-L2) significantly increase $Seff\u22c6$, making them an ideal choice for effective reduction with minimal offline costs. It is interesting to note that the POD-2 approach leads to a loss of accuracy in the reduced solution, whereas the simulation-free L2 approach with the same size of the reduction basis performs much better in terms of the accuracy (cf., Figs. 4 and 5).

### Wing Structure.

where *H*(*t*) is the Heaviside function and *ω* chosen as the average of the first- and second-natural frequency of vibration. The load amplitude is chosen so that the linear and nonlinear internal forces have magnitudes of similar order (see Fig. 7, cf., Ref. [13] for linear and nonlinear response of this model). This is in agreement with the range of applicability of von-Kármán-kinematics-based shell elements, used here.

The solution to the linearized system can be accurately reproduced using modal superposition with the first five VMs of the structure. These VMs (*M* = 5), along with the corresponding SMDs are used to construct reduction basis to perform simulation-free reduction in L1 (*m* = 20) and L2 (*m* = 10), as discussed in Sec. 6.1. The ROMs L1 and L2 are further equipped with the QM-based training-set generation technique to perform simulation-free hyper-reduction using ECSW. The results for these techniques are compared with the classical (hyper-)reduction approaches and reported in Table 2.

Type | Reduction technique | Basis size (m) | # Elements | GRE (%)_{M} | $S\u22c6$ | $Seff\u22c6$ |
---|---|---|---|---|---|---|

Classical reduction | POD-1/POD-2 | 20/10 | 49,968 | 0.31/0.98 | 2.77/2.51 | 0.73/0.71 |

Classical hyper-reduction | ECSW-POD | 20 | 278 | 0.38 | 394.6 | 0.97 |

Simulation-free reduction | L1/L2 | 20/10 | 49,968 | 1.65/1.84 | 2.72/2.74 | 2.69/2.72 |

Simulation-free hyper-reduction | Modal-ECSW-L1 | 20 | 156 | 1.40 | 638.7 | 38.48 |

Modal-ECSW-L2 | 10 | 70 | 1.61 | 1355 | 44.38 |

Type | Reduction technique | Basis size (m) | # Elements | GRE (%)_{M} | $S\u22c6$ | $Seff\u22c6$ |
---|---|---|---|---|---|---|

Classical reduction | POD-1/POD-2 | 20/10 | 49,968 | 0.31/0.98 | 2.77/2.51 | 0.73/0.71 |

Classical hyper-reduction | ECSW-POD | 20 | 278 | 0.38 | 394.6 | 0.97 |

Simulation-free reduction | L1/L2 | 20/10 | 49,968 | 1.65/1.84 | 2.72/2.74 | 2.69/2.72 |

Simulation-free hyper-reduction | Modal-ECSW-L1 | 20 | 156 | 1.40 | 638.7 | 38.48 |

Modal-ECSW-L2 | 10 | 70 | 1.61 | 1355 | 44.38 |

The following observations can be made from the results in Table 2:

The general need for hyper-reduction is apparent from the results for classical reduction as well as simulation free reduction approaches. Even for a model with truly high number of DOFs (as Model-II here), the online speed-up ($S\u22c6$) shows a value between 2 and 3 for these techniques. Indeed, the evaluation and projection of nonlinearity is a major bottleneck for saving computational time.

For the conservative choice of equal weights assigned to online and offline costs, the Classical reduction methods (POD-1/2) show an effective speed-up $Seff\u22c6<1$, even when equipped with hyper-reduction (ECSW-POD). This makes such techniques heavily dependent on an expensive database of full solution runs to obtain any reasonable effective speed-up for a range of load cases.

The simulation-free reduction techniques, on the other hand, lead to effective speed-ups $Seff\u22c6>1$. Indeed, the reduction basis constructed using modes and SMDs is much cheaper to obtain than a POD basis, which requires the HFM solution vectors.

The simulation-free hyper-reduction using the QM-based training-set generation, when applied to the ROMs L1 and L2, results in effective speed-ups which are orders of magnitude higher than any other methods. Specifically, the online speed-up for Modal-ECSW-L2 is approximately four times that of ECSW-POD. This is due to the fact that the sNNLS routine selected a quarter the number of elements for Modal-ECSW-L2 in comparison with ECSW-POD (online speed-up $S\u22c6$ for ECSW is inversely proportional to the number of elements selected in sNNLS routine). Similar arguments can be made for Modal-ECSW-L1, indicating that the QM-based training set generation results in a smaller (and, hence, effective) set of elements for ECSW.

## Discussion on Alternative Simulation-Free Reduced-Order Models Constructions

Model reduction based on vibration modes for geometrically nonlinear systems has been a topic of extensive research. As an alternative to simulation-free hyper-reduction, several methods approach the problem by assuming that the reduced nonlinear force is an at most cubic function of the reduced variables, and determining the nonlinear ROM coefficients offline ([16–20], cf., Ref. [11] for a review). The evaluation of these reduced nonlinear coefficients can be done

- (1)
*In a nonintrusive fashion*: using a so-called stiffness evaluation procedure, whereby the structure is prescribed with numerous static displacements which are carefully selected, and a large set of nonlinear algebraic equations need to be solved to be able to determine these coefficients, - (2)
*Intrusively*: for low order polynomial nonlinearities, a direct method involves the evaluation these coefficients at the element level, exactly and intrusively (see, e.g., Ref. [19]).

Model reduction based on nonlinear coefficient evaluation and natural modes has some drawbacks. It often involves the selection of a few VMs in the reduction basis to capture geometrically nonlinear effects in a rather ad hoc manner. Typically, even bending dominated response of geometrically nonlinear structures requires membrane modes to be included in the reduction basis. To the best of the authors' knowledge, there are no methods which directly give a set of vibration modes relevant for the reduction of geometrically nonlinear problems. Rather, as done for instance in Refs. [17–19] and [27], membrane-dominated modes need to be manually selected from the high-frequency spectrum. One must therefore compute many modes before finding the required additional modes. More importantly, it is often very difficult to judge which modes are required to capture the nonlinear coupling, especially for structures with complex geometries for which the membrane/bending distinction is not well defined. Implicit condensation and expansion (cf., e.g., Ref. [28]), which involves the identification of a set of membrane modes using a series of static solutions, does away with this issue to a large extent. However, selecting relevant static fields remains a challenging and cumbersome task to ensure success with this approach. To this end, the use of MDs, along with a few low-frequency modes easily identifiable from linear analysis, is a much straightforward way of capturing the geometrically nonlinear coupling (cf., Refs. [13], [15], [29], and [30]).

The indirect computation of the nonlinear coefficients is prone to several numerical challenges. Consider, for instance, a nonlinear ROM constructed using 30 modes with up to cubic nonlinearities. As many as 34,280 nonlinear static load cases are needed to completely specify the nonlinear ROM. As the underlying HFM is assumed to be high-dimensional, the computational burden associated with these nonlinear solutions/evaluations^{6} can be overwhelming [31]. Although a significant reduction in the number of coefficients that need to be identified was achieved in Ref. [27], by taking into account the tangent matrix of the nonlinear identification problems, one still needs load cases coming from well-selected static displacement fields to identify the nonlinear coefficients. In addition, appropriate modal amplitudes that trigger the nonlinearity to the right amount need to be prescribed, making the entire procedure somewhat difficult.

The direct computation of tensor coefficients in an intrusive manner does away with much of these numerical challenges, since the tensors can be calculated exactly in a systematic manner (see, e.g., [19] and [30]). However, it is strictly applicable for polynomial nonlinearities in the full system. In particular, direct tensor computation cannot be applied to other types of elements where analytical polynomial expressions for internal force are not available, e.g., corotational finite elements.

Furthermore, nonlinear coefficients evaluation restricts the nonlinearities in reduced equations to low-degree polynomials, usually at most cubic. One reason for this is the computation of higher-order tensors becoming intractable very quickly with the size of the reduction basis, both in terms of memory requirements as well as speed [30]. This restriction is acceptable when the full model has only up to cubic nonlinearities, e.g., in the case of von Kármán nonlinearities. Thus, for other cases with higher degrees of nonlinearity: while it is convenient to look only at lower order terms, useful physical effects might be ignored. Hyper-reduction methods address all the above issues to a large extent.

In a broad sense, Hyper-reduction is a semi-intrusive^{7} method, which facilitates very fast approximation of arbitrary nonlinearities. In ECSW-based hyper-reduction, this is done by evaluating the nonlinear terms over a reduced element set and assigning appropriate weights to the elements in such a set. Training vectors are required to obtain this reduced set of elements and corresponding weights in a systematic manner. Hence, hyper-reduction does not face the same issues encountered in nonlinear coefficient evaluation. However, the training vectors, which usually come from full simulations, carry a huge and undesirable offline costs. We have attempted to alleviate this bottleneck in the context of geometrical nonlinearities.

Though the technique proposed here was demonstrated for von-Kármán kinematics-based shell structures, the method is fairly general to include arbitrary geometrical nonlinearities and it will be a subject of future efforts. Moreover, the use of the quadratic manifold has been shown to be effective in capturing the geometrically nonlinear response of solid finite element-based structures as well [14].

In the context of reducing offline costs for hyper-reduction, some very recent efforts have been made in Ref. [32], where a sequence of Krylov force vectors are used to somehow approximate the subspace spanned by the inertial forces associated with the full system. The basis vectors in this Krylov force subspace are then amplified with random Gaussian amplitudes. Finally, these force realizations are applied to the system in order to solve for nonlinear static displacements, which are then used as training vectors for hyper-reduction. The following preliminary observations can be made about the approach in Ref. [32] in comparison with the QM-based lifting approach proposed here:

- (1)
Since a set of random force realizations are used for training, this method is shown to be capable of handling large displacements, unlike the method proposed here.

- (2)
It is not known (a priori or a posteriori) how many moments (or basis vectors) in the Krylov force subspace would be required to suitably approximate the inertial forces in the system. The same holds true for the standard deviation associated to the Gaussian distribution required for generating physically relevant random amplitudes and the number of random realizations needed to ensure good training. Though the approach has a wide scope, it seems that the user would need to adjust these parameters on a case-by-case basis.

- (3)
Since the solution of nonlinear static displacements for each random realization can be a fairly expensive affair, one might incur fairly high offline costs in the process. Indeed, even after use of parallelization in the offline stage, the hyper-reduction based on the approach in Ref. [32] delivered an effective speed-up factor of only about 2.15 on a realistic problem with size comparable to the wing example presented in this work. Note that such an effective speed-up is even lower than that obtained using simulation-free reduction methods without hyper-reduction (cf., Table 2).

- (4)
Finally, by taking into account randomly generated amplitudes, one might over-train for hyper-reduction, resulting in a relatively large set of selected elements. This leads to limited online performance as can be deduced from the results in Ref. [32], where an online speed-up factor of less than 5 is obtained for a realistic example with size comparable to model-II (cf., Table 2 for the corresponding online speedups).

## Conclusions

We have introduced a new method to generate training sets for hyper-reduction of geometrically nonlinear structural dynamics problems, without the need of full solution snapshots thereby reducing offline costs significantly. This method essentially involves the lifting of snapshots obtained from the (inexpensive) linear modal superposition solution, on to a quadratic manifold, which is tangent to the corresponding linear modal subspaces and captures the second-order nonlinear effects. As discussed in Sec. 5.2, this modal-based technique results in physically meaningful training vectors, which capture the essential bending-stretching coupling, characteristic of geometrically nonlinear thin-walled structures. The advantage of this method lies in the fact that the linear modal solution is available for any given system at practically no costs.

As remarked earlier, the QM-based training sets still rely on the underlying linear footprint of the model and correct for the relevant nonlinear effects up to the second-order. It is well known that von Kármán kinematics is valid up to deflections that produce linear and nonlinear internal forces of comparable magnitude (cf., Fig. 7). In the experience of the authors, the proposed method has been found to perform well for such range of deflections. However, for higher displacements, the full solution may not evolve over a quadratic manifold, as also noted in Ref. [14]. The quadratic manifold is not expected to return physically relevant training vectors for hyper-reduction in such cases. Especially for resonant loading of lightly damped systems, the modal amplitudes are expected to grow significantly for a linearized system (until steady-state is reach over long times). In such situations, the system response may be significantly more nonlinear than predicted by the quadratic manifold lifting.

The computational speed-ups calculated in the numerical results conservatively assume that the online and offline stages of the simulation carry equal costs. However, it should also be noted that, in general, the full-solution vectors used for training and reduction purposes can be expected to be optimal only for the load case from which these full-solution vectors are initially obtained. Thus, practically, a database of full-solution runs for different load cases would be needed before any meaningful benefits of model-reduction could be observed. We contend the use of simulation-free reduction techniques to ease preliminary geometrically nonlinear analysis of structures, where such an expensive database of full solutions is unavailable or unaffordable. Thus, we think that the naively defined effective speed-ups used here are still indicative.

The QM-based training set generation is a physically and intuitively inspired approach. The range of applicability of this technique and the associated error-bounds can yet not be a priori established. An abstract and analytical approach would be needed to come up with such a priori estimates, which forms a part of our future endeavors.

## Acknowledgment

The authors are thankful to anonymous reviewers, whose feedback was very valuable for this work.

## Funding Data

Air Force Office of Scientific Research, Air Force Material Command, USAF (FA9550-16-1-0096).

## Nomenclature

- DEIM =
discrete empirical interpolation method

- DOF =
degrees-of-freedom

- ECSW =
energy conserving sampling and weighing

- EIM =
empirical interpolation method

- FE =
finite element

- GRE =
global relative error

- HFM =
high-fidelity model

- MD =
modal derivative

- NNLS =
non-negative least-squares

- ODE =
ordinary differential equation

- POD =
proper orthogonal decomposition

- QM =
quadratic manifold

- ROM =
reduced-order model

- SMD =
static modal derivative

- SVD =
singular value decomposition

- UDEIM =
unassembled discrete empirical interpolation method

- VM =
vibration mode

Without loss of generality, we assume a zero mean of the training samples in **U**.

Here we neglect the damping contribution in eigenvalue problem to avoid complex eigenvalues and vectors. Note that for damped linear systems with low damping or $modal/Rayleigh$ damping as explained in Ref. [25], the eigenvectors for an undamped system are a good approximation for the damped counterpart and still form a good basis for linear modal superposition. Such a damping is very popular in structural dynamics and the theory is illustrated in this context. Under large damping, the real and imaginary parts of the most significant complex eigenvectors would be added to the matrix $\Psi $.

The mapping Eq. (17) can be written using the Einstein summation convention in the indicial notation as: $\Gamma I=\Phi Iizj+12\Omega Iijzizj\u2003I\u2208{1,\u2026,n},\u2003i,j\u2208{1,\u2026,M}.$

Note that due to the spatially uniform nature of the loading, the second, third and fourth VMs are not excited (since they is anti-symmetric bending and twisting modes), even though the applied loading is at the second natural frequency of the model. Thus, the loading frequency is not a resonant frequency for the ROM.

Depending on whether a force based or displacement based coefficient evaluation is used, cf., Ref. [11].

We refer to hyper-reduction as semi-intrusive because although the element level contributions of nonlinear force residual and Jacobian are required for hyper-reduction, the actual expressions for these element-level nonlinearities are not needed. In other words, the element is treated as a black box.