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,36]. 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) [79]. 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. [1012]). 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 [1315]. 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 [1620]. 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

The partial differential equations for momentum balance in a structural continuum are first FE discretized along the spatial dimensions to obtain a system of second-order ordinary differential equations (ODEs). Along with the initial conditions for generalized displacements and velocities, these ODEs govern the response of the underlying structure. More specifically, this response can be described by the solution to an initial value problem of the following form:
$Mu¨(t)+Cu˙(t)+f(u(t))=g(t)u(t0)=u0, u˙(t0)=v0$
(1)
where the solution $u(t)∈ℝn$ is a high-dimensional generalized displacement vector, $M∈ℝn×n$ is the mass matrix, $C∈ℝn×n$ is the damping matrix, $f:ℝn↦ℝn$ gives the nonlinear elastic internal force as a function of the displacement u of the structure, and $g(t)∈ℝ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
$u(t)≈Vq(t) , V∈ℝn×m$
where $q(t)∈ℝm$ ($m≪n$) is the low-dimensional vector of reduced variables, and V is known as the reduction basis since its columns form a basis for $V$. The ROM is then obtained using Galerkin projection as
$VTMV︸M̃q¨(t)+VTCV︸C̃q˙(t)+VTf(Vq(t))=VTg(t)$
where $M̃, C̃∈ℝm×m$ are the reduced mass and damping matrices, respectively. Often, the internal force can be split into its linear and nonlinear contributions as $f(u)=Ku+fnl(u)$ to obtain a reduced stiffness matrix as well
$M̃q¨(t)+C̃q˙(t)+VTKV︸K̃q(t)+VTfnl(Vq(t))︸f̃(q(t))=VTg(t)$
(2)
It is easy to see that each of the reduced matrices $M̃, C̃$, and $K̃$ can be precomputed in the offline stage prior to time integration. Hence, during time integration (online phase), the computational cost associated with the evaluation of the linear terms in Eq. (2) scales only with the number of reduced variables m. This is, however, not the case for the computation of the nonlinear term $f̃(q(t))$. For FE-based applications, this evaluation is usually carried out online in the following manner:
$f̃(q)=VTfnl(Vq)=∑e=1neVeTfe(Veq)$
(3)

where $fe(ue)∈ℝNe$ is the contribution of the element e toward the vector $fnl(u)$ (Ne 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 ne is the total number of elements in the structure. Since the reduced nonlinear term $f̃(q)$ is evaluated in the space of full variables, the computational cost associated with its evaluation does not scale with 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̃(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̃(q)$ without the loss of numerical stability [21], is remarkable.

Essentially, ECSW aims to identify a small set of elements E of the structure ($|E|≪ne$) to cheaply approximate $f̃(q)$ as (cf., Eq. (3))
$f̃(q)=∑e=1neVeTfe(Veq)≈∑e∈EξeVeTfe(Veq)$
(4)

where $ξe∈ℝ+$ is a positive weight attached to each element $e∈E$, which is empirically chosen to ensure a good approximation of the summation in Eq. (3). Clearly, if $|E|≪ne$, 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]).

The elements and weights are determined to approximate virtual work over the chosen training sets, which generally come from full solution run(s). If there are nt training vectors in the set with $u(i)$ representing the ith vector, then corresponding reduced unknowns $q(i)$ can be easily calculated using least squares as
$q(i)=(VTV)−1VTu(i)$
and element level contribution of projected internal force for each of the training vectors can be assembled in a matrix G as follows:
$G=[g11…g1ne⋮⋱⋮gnt1…gntne]∈ℝmnt×ne, b=[b1⋮bnt]∈ℝmntgie=VeTfe(Veq(i)), bi=f̃(q(i))=∑e=1negie , ∀i∈{1,…,nt}, e∈{1,…,ne}$
(5)
The set of elements and weights is then obtained by a sparse solution to the following non-negative least-squares (NNLS) problem:
$(P1):ξ=argminξ̃∈ℝne,ξ̃≥0‖Gξ̃−b‖2$
(6)
A sparse solution to $(P1)$ returns a sparse vector $ξ$, the nonzero entries of which form the reduced mesh E used in Eq. (4) as
$E={e:ξe>0}$

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 $ζE$ and $GE$ denote, respectively, the restriction of $ζ∈ℝ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 $ξ$ and $ζ$.

Algorithm 1

Sparse NNLS for ECSW

 Input:$G,b,τ$ Output:$ξ∈ℝne$ sparse, $E⊂{1,…,ne}$ 1: $E←∅,Z←{1,…,ne},ξ←0∈ℝne$ 2: while$‖Gξ−b‖2>τ‖b‖2$do 3:  $μ←GT(b−Gξ)$ 4:  $[|ν|,e]=max(μ)$ ▹ $max$⁠: returns maximum value in a vector followed by its location (index) 5:  $E←E∪{e},Z←Z\{e}$ 6:  while true do 7:   $ζE←GE†b$ ▹$†$ represents pseudo-inverse 8:   $ζZ←0$ 9:   if$ζE>0$then 10:    $ξ←ζ$ 11:    break 12:   $η=mink∈Eξk/(ξk−ζk)$ 13:    $ξ←ξ+η(ζ−ξ)$ 14:    $Z←{i|ξi=0}$ 15:    $E←{1,…,ne}\Z$
 Input:$G,b,τ$ Output:$ξ∈ℝne$ sparse, $E⊂{1,…,ne}$ 1: $E←∅,Z←{1,…,ne},ξ←0∈ℝne$ 2: while$‖Gξ−b‖2>τ‖b‖2$do 3:  $μ←GT(b−Gξ)$ 4:  $[|ν|,e]=max(μ)$ ▹ $max$⁠: returns maximum value in a vector followed by its location (index) 5:  $E←E∪{e},Z←Z\{e}$ 6:  while true do 7:   $ζE←GE†b$ ▹$†$ represents pseudo-inverse 8:   $ζZ←0$ 9:   if$ζE>0$then 10:    $ξ←ζ$ 11:    break 12:   $η=mink∈Eξk/(ξk−ζk)$ 13:    $ξ←ξ+η(ζ−ξ)$ 14:    $Z←{i|ξi=0}$ 15:    $E←{1,…,ne}\Z$

## Classical (Hyper) Reduction

We briefly review the classical notion of reduction and hyper-reduction with the use of HFM solution vectors. Let $ui∈ℝn$, $i∈{1,…,nt}$ be training vectors obtained from the full solution of system Eq. (1). Let $U:=[u1,u2,…,unt]∈ℝn×nt$ be the ensemble of snapshots obtained from the full solution. A lower dimensional POD basis $V=[v1v2…vm]∈ℝn×m$ containing $m≪nt$ orthogonal vectors which best spans the vectors in this ensemble can be obtained by the solution to the following minimization problem:
$minvi∈ℝn∑j=1ns‖uj−∑i=1m(ujTvi)vi‖22$
This is a least squares problem and the vectors in V are the left singular vectors of U, obtained by the singular value decomposition (SVD) of U as2
$U=LSRT$
(7)
where is U is factorized into unitary matrices $L=[l1,l2,…,ln]∈ℝn×n$ (containing the left singular vectors) and $R∈ℝnt×nt$ (containing the right singular vectors); and the diagonal (rectangular) matrix $S∈ℝn×nt$ (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≥S22≥⋯≥0$, it can be shown that
$∑j=1ns‖uj−∑i=1m(ujTli)li‖22=∑i=m+1rSii2$

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.

Algorithm 2

ECSW-POD

 Input:$u(i),i∈{1,…,nt}$ snapshots from full solution Output:$ξ∈ℝne$ sparse, $E⊂{1,…,ne}$ for ECSW-POD 1: $U←[u(1),…,u(nt)]$ 2: $[V]←POD(U,M)$ 3: $[G,b]←CONSTRUCTGb(U,V)$a 4: $[ξ,E]←SNNLS(G,b)$b
 Input:$u(i),i∈{1,…,nt}$ snapshots from full solution Output:$ξ∈ℝne$ sparse, $E⊂{1,…,ne}$ for ECSW-POD 1: $U←[u(1),…,u(nt)]$ 2: $[V]←POD(U,M)$ 3: $[G,b]←CONSTRUCTGb(U,V)$a 4: $[ξ,E]←SNNLS(G,b)$b
a

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

b

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.

For geometrically nonlinear problems with separation in time scales between in-plane and out-of-plane behavior (beams, shells, thin plates, and their assemblage), it is well known that the VMs and MDs form a meaningful basis to represent the nonlinear displacement field [13,15,23,24]. We propose their use in construction of the required reduction basis V as
$Ψ=[ϕ1,…,ϕM,…,θij,…], i,j∈{1,…,M}$
(8)
$V=orth(Ψ)$
(9)
where V is obtained after an orthogonalization of the matrix $Ψ$ 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 $Ψ$ contains the relevant VMs and MDs, $ϕ1,…,ϕM$ represent a truncated set of VMs obtained from the solution of the undamped eigenvalue problem3
$(K−ωi2M)ϕi=0 ∀i∈{1,…,n}$
(10)
and $θij$ are the MDs, obtained from the solution to the following problem:
$(K−ωi2M)θij+(∂2f(u)∂u∂u|u=0ϕj−∂ωi2∂ηjM)ϕi=0, i,j∈{1,.…,M}$
(11)
This represents a derivative of the eigenvalue problem Eq. (10) with respect to the modal amplitude ηj of the VM $ϕj$, after K is replaced by the tangent stiffness matrix $(∂f(u)/∂u)$. The singular system Eq. (11) requires an additional normalization condition to solve for $θij$. We use the mass normalization of the VMs, which results in the following normalization constraint:
$ϕiTMθij=0 ∀i,j∈{1,.…,M}$
Various techniques exist to obtain MDs (cf., Ref. [13] for a brief review). These include a somehow effective static approximation of the MDs, obtained by neglecting the mass terms in Eq. (11) as
$θijstatic=−K−1[∂2f(u)∂u∂u|u=0ϕj]ϕi$
(12)

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

It is also easy to see that the size of the matrix $Ψ$ is $O(M2)$, if all MDs corresponding to 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 Wij to an SMD $θijstatic$:
$Wij=∫0T|ηi(t)ηj(t)|dt, i,j∈M$
(13)

where $ηi(t)$ and $ηj(t)$ are the modal amplitudes for VMs $ϕi$ and $ϕ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≤2M$ 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.

For proportionally damped systems, the linearized solution $ulin(t)$ of the HFM in response to the applied external loading $g(t)$ can be spectrally decomposed using linear modal superposition as (cf., Ref. [25])
$η¨i+βiη˙i(t)+ωi2η=γi(t)$
(14)
$ulin(t)=∑i=1nϕiηi(t)$
(15)
where $ηi(t)$ are the time varying modal amplitudes corresponding to the $M$-normalized VM $ϕi$; $βi:=ϕiTCϕi$ is the damping coefficient resulting from proportional damping assumption; and $γi(t):=ϕiTg(t)$ denotes the modal forcing. Note that Eq. (14) is a system of decoupled linear ODEs which can be very easily integrated. Depending on the spatial and temporal properties of the external load $g(t)$, the summation in Eq. (15) can be truncated to obtain
$ulin(t)≈∑i∈Mϕiηi(t)=Φη(t)$
(16)

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

Clearly, the snapshots of the linearized solution $ulin(t)$ would not be effective in training for nonlinear hyper-reduction. Indeed, consider, for example, a flat thin walled shell structure with some externally applied out-of-plane loading. The modal truncation would feature bending dominated modes which are out of plane for a linearized solution. The nonlinear response, however, is expected to feature in-plane effects due to bending/torsion-stretching coupling. These essential features would be missing in snapshots of $ulin(t)$, making them poor for training nonlinear hyper-reduction. The modal amplitudes, however, still contain important bending information of the structure and can be useful in generating effective training vectors by using the notion of a QM, which was introduced for nonlinear model reduction in Ref. [13]. There, the full unknowns are mapped to a lower dimensional set of variables z using a quadratic mapping as
$u(t)≈Γ(z(t)):=Φz(t)+12Ω:(z(t)⊗z(t))$
(17)
where $Φ∈ℝn×M$ containing a truncated set of $M(=|M|)$ significant VMs, resulting from the linear modal superposition; and $Ω∈ℝn×M×M$ is a third-order tensor containing the corresponding MDs.4 An arbitrary vector of modal coordinates is denoted by $z∈ℝM$. As discussed in Ref. [13], the tangent space of the QM at equilibrium ($u=0$) is the modal subspace represented by $Φ$, 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 $η(t)$ in Eq. (17), as
$u(i)=Γ(η(ti)), i∈{1,…,nt}$

where $ti∈T$ (with $T$ being the simulation time-span) are the time instants at which the modal snapshots are captured; nt 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.

Fig. 1
Fig. 1
Close modal

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.

Algorithm 3

Simulation-free hyper-reduction using ECSW

 Input: HFM model as described by Eq. (1), tolerance for ECSW-based training τ Output:$ξ∈ℝne$ sparse, $E⊂{1,…,ne}$ for ECSW-Modal 1: $[Φ,η(t)]←ModalSuperposition(model,g(t),T)$a 2: $M←size(Φ,2)$b Calculation of (S)MDs to construct$Ω$in Eq. (17) 3: $Ω←$ all-empty 3-tensor$∈ℝn×M×M$⁠, $Ψ←Φ$ 4: for$j←1 to m$do 5:  Compute $(∂2f(u)/∂u∂u)|u=0ϕj$ 6:  for$i←1 to M$do 7:   $ω←θij$ or $θijstatic$c 8:   $Ω(:,i,j)←ω$⁠, $Ψ←[Ψ,ω]$ ▹ matlab notation Perform selection of MDs in$Ψ$, if needed 9: if MD selection = true then 10:  $Ψ←MDSelection(Ψ,η(t))$d 11: $[V]←orth(Ψ)$e Generate training vectors 12: Build function $QM$⁠: $Γ(z)←QM(z)$f 13: for$i∈{1,2,…,nt}$do 14:  $u(i)=QM(η(ti))$g 15: $U←[u(1),…,u(nt)]$ Perform training 16: $[G,b]←CONSTRUCTGb(U,V)$h 17: $[ξ,E]←SNNLS(G,b)$i
 Input: HFM model as described by Eq. (1), tolerance for ECSW-based training τ Output:$ξ∈ℝne$ sparse, $E⊂{1,…,ne}$ for ECSW-Modal 1: $[Φ,η(t)]←ModalSuperposition(model,g(t),T)$a 2: $M←size(Φ,2)$b Calculation of (S)MDs to construct$Ω$in Eq. (17) 3: $Ω←$ all-empty 3-tensor$∈ℝn×M×M$⁠, $Ψ←Φ$ 4: for$j←1 to m$do 5:  Compute $(∂2f(u)/∂u∂u)|u=0ϕj$ 6:  for$i←1 to M$do 7:   $ω←θij$ or $θijstatic$c 8:   $Ω(:,i,j)←ω$⁠, $Ψ←[Ψ,ω]$ ▹ matlab notation Perform selection of MDs in$Ψ$, if needed 9: if MD selection = true then 10:  $Ψ←MDSelection(Ψ,η(t))$d 11: $[V]←orth(Ψ)$e Generate training vectors 12: Build function $QM$⁠: $Γ(z)←QM(z)$f 13: for$i∈{1,2,…,nt}$do 14:  $u(i)=QM(η(ti))$g 15: $U←[u(1),…,u(nt)]$ Perform training 16: $[G,b]←CONSTRUCTGb(U,V)$h 17: $[ξ,E]←SNNLS(G,b)$i
a

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

b

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

c

Use MDs from Eq. (11) or SMDs from Eq. (12) in the quadratic part.

d

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

e

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

f

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

g

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

h

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

i

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

## Numerical Results

### Setup.

We illustrate the performance of the proposed hyper-reduction techniques on two examples of thin-walled structures. One of them is an academically simple model of a shallow-arched rectangular plate, named model-I. The other example is the model of a NACA-airfoil-based wing model, with realistically high number of DOFs, introduced in Ref. [13], and referred to as model-II. The details of both models are given in Secs. 6.2 and 6.3, respectively. Both structures are modeled using flat, triangular-shell elements with six DOFs per node (i.e., 18 DOFs per element). For both the models, the accuracy of the results was compared to the corresponding full nonlinear solutions using a mass-normalized global relative error (GRE) measure, defined as
$GREM=∑t∈S(u(t)−ũ(t))TM(u(t)−ũ(t))∑t∈Su(t)TMu(t)×100$

where $u(t)∈ℝn$ is the vector of generalized displacements at the time t, obtained from the HFM solution, $ũ(t)∈ℝ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.

Fig. 2
Fig. 2
Close modal
Fig. 3
Fig. 3
Close modal
All computations were performed in matlab on the Euler cluster of ETH Zürich. Since the success of reduction techniques is often reported in terms of savings in simulation time, we define an online speedup $S⋆$ computed according to the following simple formula:
$S⋆=TfullTsim⋆$
where Tfull and $Tsim⋆$ represent the computation time taken during the time integration of $full$ and (hyper)reduced solution respectively. The superscript $⋆$ 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
$Seff⋆=conTfullcoffToff⋆+conTsim⋆$
(18)

where $Toff⋆$ represents the computational time spent $offline$ for setting up a (hyper-)reduced model; $con,coff∈[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 Seff would make a given (hyper-) reduction method favorable, with Seff < 1 implying that the corresponding reduction technique is effectively more expensive than a full solution run.

One method of choosing con and coff may be based on the number of simulations nsim performed using a ROM as $con=(nsim/nsim+1), coff=(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:

1. 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)$.

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$.

3. Classical reduction (no hyper-reduction): Here, the nonlinear system using is reduced using a POD basis.

4. 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.

5. 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.

6. 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.

7. 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:

8. 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.

9. Modal-ECSW-L2: Same as Modal-ECSW-L1, except the ROM from L2 is hyper-reduced, instead of that from L1.

### Shallow-Arch Structure.

We consider a simple arch structure, as shown in Fig. 2. A uniform pressure distribution is chosen to act normal to the plate surface as the external load. A time varying amplitude (load function) is used given by
$g(t)=p(t)lp(t)=p0[sin(ωt)+sin(πωt)]$
(19)

where l is a constant load vector corresponding to a uniform pressure distribution of 1 N/mm2. 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.

Fig. 4
Fig. 4
Close modal

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.

Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal
Table 1

Starting with a linear modal superposition with M = 2 modes (first and the fifth), reduction techniques as described in Sec. 6.1 are formulated. Global relative error, speed-ups and effective speed-ups (for $con=coff=0.5$) in model-I for these reduction techniques are tabulated. A total of nt = 200 training vectors, chosen uniformly from the solution time-span, are used to setup hyper-reduction in Modal-ECSW-L1/L2 as well as ECSW-POD. Computational time for the full solution run for model-I is given by Tfull = 66.34 s.

TypeReduction techniqueBasis size (m)# ElementsGREM (%)$S⋆$$Seff⋆$
Classical reductionPOD-154002.862.850.74
Classical hyper-reductionECSW-POD5193.0450.40.94
Simulation-free reductionL154003.362.882.85
Simulation-free hyper-reductionModal-ECSW-L15114.9579.915.9
TypeReduction techniqueBasis size (m)# ElementsGREM (%)$S⋆$$Seff⋆$
Classical reductionPOD-154002.862.850.74
Classical hyper-reductionECSW-POD5193.0450.40.94
Simulation-free reductionL154003.362.882.85
Simulation-free hyper-reductionModal-ECSW-L15114.9579.915.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⋆<1$, as expected from the very definition of $Seff⋆$ (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⋆>1$. The simulation-free hyper-reduction techniques (Modal-ECSW-L1 and Modal-ECSW-L2) significantly increase $Seff⋆$, 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.

For the application of the proposed simulation-free hyper-reduction methods to more realistic models, we consider the model of an NACA-airfoil-based wing structure, introduced in Ref. [13]. This model (referred to as Model-II hereafter) contains truly high number of DOFs, thereby allowing for the appreciation of obtained accuracy and computational speed-ups. We simulate the response of the structure to a low-frequency pulse load, applied as a spatially uniform pressure load on the highlighted area on the structure skin (cf., Fig. 6). This pressure load takes the shape of a pulse in time as shown in Fig. 7. The dynamic load function is given as
$p(t)=A sin2(ωt)[H(t)−H(πω−t)]$
(20)

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.

Fig. 7
Fig. 7
Close modal

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.

Table 2

Starting with a linear modal superposition with M = 5 modes (first five VMs), reduction techniques as described in Sec. 6.1 are formulated. Global relative error, speed-ups and effective speed-ups (for $con=coff=0.5$) in model-II for these reduction techniques are tabulated. A total of nt = 200 training vectors, chosen uniformly from the solution time-span, are used to setup hyper-reduction in Modal-ECSW-L1/L2 as well as ECSW-POD. Time for a full solution run $Tfull=3.808×104$ s. We see from the boldfaced values that the simulation-free hyper-reduction methods lead to selection of less number of elements in comparison with classical hyper-reduction. This results in better online performance as indicated by the respective S values. Furthermore, the effective speed-up Seff is much better for simulation-free hyper-reduction methods in comparison with the classical hyper-reduction, as expected.

TypeReduction techniqueBasis size (m)# ElementsGREM (%)$S⋆$$Seff⋆$
Classical reductionPOD-1/POD-220/1049,9680.31/0.982.77/2.510.73/0.71
Classical hyper-reductionECSW-POD202780.38394.60.97
Simulation-free reductionL1/L220/1049,9681.65/1.842.72/2.742.69/2.72
Simulation-free hyper-reductionModal-ECSW-L1201561.40638.738.48
Modal-ECSW-L210701.61135544.38
TypeReduction techniqueBasis size (m)# ElementsGREM (%)$S⋆$$Seff⋆$
Classical reductionPOD-1/POD-220/1049,9680.31/0.982.77/2.510.73/0.71
Classical hyper-reductionECSW-POD202780.38394.60.97
Simulation-free reductionL1/L220/1049,9681.65/1.842.72/2.742.69/2.72
Simulation-free hyper-reductionModal-ECSW-L1201561.40638.738.48
Modal-ECSW-L210701.61135544.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⋆$) 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⋆<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⋆>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⋆$ 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 ([1620], cf., Ref. [11] for a review). The evaluation of these reduced nonlinear coefficients can be done

1. (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. (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. [1719] 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/evaluations6 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-intrusive7 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. (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. (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. (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. (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 =

•
• ROM =

reduced-order model

•
• SMD =

static modal derivative

•
• SVD =

singular value decomposition

•
• UDEIM =

unassembled discrete empirical interpolation method

•
• VM =

vibration mode

2

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

3

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 $Ψ$.

4

The mapping Eq. (17) can be written using the Einstein summation convention in the indicial notation as: $ΓI=ΦIizj+12ΩIijzizj I∈{1,…,n}, i,j∈{1,…,M}.$

5

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.

6

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

7

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.

## References

1.
Riks
,
E.
,
1997
, “
Buckling Analysis of Elastic Structures: A Computational Approach
,”
,
34
, pp.
1
76
.
2.
Farhat
,
C.
,
Avery
,
P.
,
Chapman
,
T.
, and
Cortial
,
J.
,
2014
, “
Dimensional Reduction of Nonlinear Finite Element Dynamic Models With Finite Rotations and Energy-Based Mesh Sampling and Weighting for Computational Efficiency
,”
Int. J. Numer. Methods Eng.
,
98
(
9
), pp.
625
662
.
3.
Chaturantabut
,
S.
, and
Sorensen
,
D. C.
,
2010
, “
Nonlinear Model Reduction Via Discrete Empirical Interpolation
,”
SIAM J. Scientic Comput.
,
32
(
5
), pp.
2737
2764
.
4.
Tiso
,
P.
, and
Rixen
,
D. J.
,
2013
, “
Discrete Empirical Interpolation Method for Finite Element Structural Dynamics
,”
Topics in Nonlinear Dynamics, Volume 1
(Conference Society for Experimental Mechanics Series, Vol.
35
), G. Kerschen D. Adams, and A. Carrella, eds.,
Springer
,
New York
.
5.
Chaturantabut
,
S.
,
Beattie
,
C.
, and
Gugercin
,
S.
,
2016
, “
Structure-Preserving Model Reduction for Nonlinear Port-Hamiltonian Systems
,”
SIAM J. Sci. Comput.
,
38
(5), pp. B837–B865.
6.
Barrault
,
M.
,
,
Y.
,
Nguyen
,
N. C.
, and
Patera
,
A. T. A.
,
2004
, “
Empirical Interpolation” Method: Application to Efficient Reduced-Basis Discretization of Partial Differential Equations
,”
C. R. Math. Acad. Sci. Paris
,
339
(
9
), pp.
667
672
.
7.
Kosambi
,
D.
,
1943
, “
Statistics in Function Space
,”
J. Indian Math. Soc.
,
7
, pp.
76
78
.http://repository.ias.ac.in/99240/
8.
Amabili
,
M.
,
Sarkar
,
A.
, and
Paidoussis
,
M. P.
,
2003
, “
Reduced-Order Models for Nonlinear Vibrations of Cylindrical Shells Via the Proper Orthogonal Decomposition Method
,”
J. Fluids Struct.
,
18
(
2
), pp.
227
250
.
9.
Kerschen
,
G.
,
Golinval
,
J. C.
,
Vakakis
,
A. F.
, and
Bergman
,
L. A.
,
2005
, “
The Method of Proper Orthogonal Decomposition for Dynamical Characterization and Order Reduction of Mechanical Systems: An Overview
,”
Nonlinear Dyn.
,
41
(
1–3
), pp.
147
169
.
10.
Bai, Z.,
2002
, “
Krylov Subspace Techniques for Reduced-Order Modeling of Large-Scale Dynamical Systems
,”
Appl. Numer. Math.
,
43
(
1–2
), pp.
9
44
.
11.
Mignolet
,
M. P.
,
Przekop
,
A.
,
Rizzi
,
S. A.
, and
Spottswood
,
S. M.
,
2013
, “
A Review of Indirect/Non-Intrusive Reduced Order Modeling of Nonlinear Geometric Structures
,”
J. Sound Vib.
,
332
(
10
), pp.
2437
2460
.
12.
Jain
,
S.
,
Tiso
,
P.
, and
Haller
,
G.
,
2018
, “
Exact Nonlinear Model Reduction for a Von Karman Beam: Slow-Fast Decomposition and Spectral Submanifolds
,”
J. Sound Vib.
,
423
, pp.
195
211
.
13.
Jain
,
S.
,
Tiso
,
P.
,
Rixen
,
D. J.
, and
Rutzmoser
,
J. B.
,
2017
, “
A Quadratic Manifold for Model Order Reduction of Nonlinear Structural Dynamics
,”
Comput. Struct.
,
188
, pp.
80
94
.
14.
Rutzmoser
,
J. B.
,
Rixen
,
D. J.
,
Tiso
,
P.
, and
Jain
,
S.
,
2017
, “
Generalization of Quadratic Manifolds for Reduced Order Modeling of Nonlinear Structural Dynamics
,”
Comput. Struct.
,
192
, pp.
196
209
.
15.
Idelsohn
,
S. R.
, and
Cardona
,
A.
,
1985
, “
A Reduction Method for Nonlinear Structural Dynamic Analysis
,”
Comput. Methods Appl. Mech. Eng.
,
49
(
3
), pp.
253
279
.
16.
Muravyov
,
A. A.
, and
Rizzi
,
S. A.
,
2003
, “
Determination of Nonlinear Stiffness With Application to Random Vibration of Geometrically Nonlinear Structures
,”
Comput. Struct.
,
81
(
15
), pp.
1513
1523
.
17.
Amabili
,
M.
,
2013
, “
Reduced-Order Models for Nonlinear Vibrations, Based on Natural Modes: The Case of the Circular Cylindrical Shell
,”
Philos. Trans. R. Soc. A
,
371
(
1993
), p.
20120474
.
18.
Lazarus
,
A.
,
Thomas
,
O.
, and
De
,
J. F.
,
2012
, “
Finite Element Reduced Order Models for Nonlinear Vibrations of Piezoelectric Layered Beams With Applications to NEMS
,”
Finite Elem. Anal. Des.
,
49
(
1
), pp.
35
51
.
19.
Touz
,
C.
,
Vidrascu
,
M.
, and
Chapelle
,
D.
,
2014
, “
Direct Finite Element Computation of Non-Linear Modal Coupling Coefficients for Reduced-Order Shell Models
,”
Comput. Mech.
,
54
(
2
), pp.
567
580
.
20.
Kuether
,
R. J.
,
Deaner
,
B. J.
,
Hollkamp
,
J. J.
, and
Allen
,
M. S.
,
2015
, “
Evaluation of Geometrically Nonlinear Reduced-Order Models With Nonlinear Normal Modes
,”
AIAA J.
,
53
(
11
), pp.
3273
3285
.
21.
Farhat
,
C.
,
Chapman
,
T.
, and
Avery
,
P.
,
2015
, “
Structure-Preserving, Stability, and Accuracy Properties of the Energy-Conserving Sampling and Weighting Method for the Hyper Reduction of Nonlinear Finite Element Dynamic Models
,”
Int. J. Numer. Methods Eng.
,
102
(
5
), pp.
1077
1110
.
22.
Peharz
,
R.
, and
Pernkopf
,
F.
,
2012
, “
Sparse Nonnegative Matrix Factorization With l0-Constraints
,”
Neurocomputing
,
80
, pp.
38
46
.
23.
Slaats
,
P. M. A.
,
de Jongh
,
J.
, and
Sauren
,
A. A. H. J.
,
1995
, “
Model Reduction Tools for Nonlinear Structural Dynamics
,”
Comput. Struct.
,
54
(
6
), pp.
1155
1171
.
24.
Barbič
,
J.
, and
James
,
D. L.
,
2005
, “
Real-Time Subspace Integration for St. Venant-Kirchhoff Deformable Models
,”
ACM Trans. Graph.
,
24
(
3
), pp.
982
990
.
25.
,
M.
, and
Rixen
,
D.
,
1997
,
Mechanical Vibrations: Theory and Application to Structural Dynamics
, 2nd ed.,
Wiley
, Chichester, UK.
26.
Tiso
,
P.
,
2011
, “
Optimal Second Order Reduction Basis Selection for Nonlinear Transient Analysis
,”
Modal Analysis Topics, Volume 3
(Conference Society for Experimental Mechanics Series),
T.
Proulx
, ed.,
Springer
,
New York
, pp.
27
39.
27.
Perez
,
R.
,
Wang
,
X. Q.
, and
Mignolet
,
M. P.
,
2014
, “
Nonintrusive Structural Dynamic Reduced Order Modeling for Large Deformations: Enhancements for Complex Structures
,”
ASME J. Comput. Nonlinear Dyn.
,
9
(
3
), p.
031008
.
28.
Hollkamp
,
J. J.
, and
Gordon
,
R. W.
,
2008
, “
Reduced-Order Models for Nonlinear Response Prediction: Implicit Condensation and Expansion
,”
J. Sound Vib.
,
318
(
4–5
), pp.
1139
1153
.
29.
Sombroek
,
C. S. M.
,
Tiso
,
P.
,
Renson
,
L.
, and
Kerschen
,
G.
,
2018
, “
Numerical Computation of Nonlinear Normal Modes in a Modal Derivatives Subspace
,”
Comput. Struct.
,
195
, pp.
34
46
.
30.
Jain
,
S.
,
2015
, “
Model Order Reduction for Nonlinear Structural Dynamics
,”
Master's thesis
, Delft University of Technology, Delft, The Netherlands.https://repository.tudelft.nl/islandora/object/uuid:cb1d7058-2cfa-439a-bb2f-22a6b0e5bb2a
31.
Schoneman
,
J. D.
,
Allen
,
M. S.
, and
Kuether
,
R. J.
,
2017
, “
Nonlinear Modal Substructuring of Panel and Stiffener Assemblies Via Characteristic Constraint Modes
,”
Dynamics of Coupled Structures, Volume 4
(Conference Proceedings of the Society for Experimental Mechanics Series), Springer, Cham, Switzerland, pp.
307
326
.
32.
Rutzmoser
,
J. B.
, and
Rixen
,
D. J.
,
2017
, “
A Lean and Efficient Snapshot Generation Technique for the Hyper-Reduction of Nonlinear Structural Dynamics
,”
Comput. Methods Appl. Mech. Eng.
,
325
, pp.
330
349
.