# An enhanced parametric nonlinear reduced order model for imperfect
structures using Neumann expansion

Jacopo Marconi  Paolo Tiso<EMAIL_ADDRESS>Davide E. Quadrelli  Francesco
Braghin  Department of Mechanical Engineering, Politecnico di Milano, Via La
Masa 1, 20156 Milan, Italy  Institute for Mechanical Systems, ETH Zürich,
Leonhardstrasse 21, 8092 Zürich, Switzerland

###### Abstract

We present an enhanced version of the parametric nonlinear reduced order model
for shape imperfections in structural dynamics we studied in a previous work
[1]. The model is computed intrusively and with no training using information
about the nominal geometry of the structure and some user-defined displacement
fields representing shape defects, i.e. small deviations from the nominal
geometry parametrized by their respective amplitudes. The linear superposition
of these artificial displacements describe the defected geometry and can be
embedded in the strain formulation in such a way that, in the end, nonlinear
internal elastic forces can be expressed as a polynomial function of both
these defect fields and the actual displacement field. This way, a tensorial
representation of the internal forces can be obtained and, owning the
reduction in size of the model given by a Galerkin projection, high simulation
speed-ups can be achieved. We show that by adopting a rigorous deformation
framework we are able to achieve better accuracy as compared to the previous
work. In particular, exploiting Neumann expansion in the definition of the
Green-Lagrange strain tensor, we show that our previous model is a lower order
approximation with respect to the one we present now. Two numerical examples
of a clamped beam and a MEMS gyroscope finally demonstrate the benefits of the
method in terms of speed and increased accuracy.

###### keywords:

Nonlinear Modeling , Reduced Order Models , Parametric , Geometric
Nonlinearities , Defects

††journal: arXiv

##  1 Introduction

The Finite Element (FE) method has long been a fundamental analysis and design
tool in many areas of science and engineering. In structural mechanics it is
almost mandatory to use FE models to investigate the behavior of complex
systems, which often have many geometric details that would be difficult to
handle with alternative approaches, such as lumped parameter or analytical
models [2]. However, large FE simulations would often require considerable
computational resources and time, so in some cases designers may prefer to
perform real experiments rather than numerical ones. On the one hand, this
need for fast and affordable FE simulations has given rise to numerical
techniques to improve computational efficiency: domain decomposition and
substructuring [3, 4] and FE Tearing and Interconnecting (FETI, [5]) are just
a few examples. On the other hand, model order reduction methods have emerged,
consisting in the construction of a Reduced Order Model (ROM), whose number of
degrees of freedom (dofs) is much smaller than that of the reference Full
Order Model (FOM). The use of linear ROMs also in industrial contexts is
nowadays well established as the theory underlying them. Guyan reduction [6]
and modal analysis [7] are two well-known examples in mechanical statics and
dynamics, respectively, where FOM’s static deformations and Vibration Modes
(VMs, also known as eigenmodes or natural modes of the linear system) are used
to construct a Reduced Basis (RB) that projects the governing equations onto a
lower dimension subspace. Linear ROMs were also successfully coupled with
substructuring techniques in the Craig-Bampton and Rubin methods [8, 9], which
are available in many commercial software.

For nonlinear FE studies, where the demand for reduction is dire, many
solutions have been proposed over the last decades, but none of them seems to
have prevailed over the others, as each of them offers certain advantages,
requires certain costs and/or targets specific problems. Overall, however, the
literature is mature enough to provide the analyst with many different options
in several practical applications, ranging from bolted joints [10], gears
[11], contacts [12, 13], friction [14] and viscoplasticity [15] to flexible
multi-body dynamics with geometric nonlinearities [16] and substructuring
[17].

Nonlinear ROMs can be classified according to (i) whether they are RB-
projection based or not, (ii) whether they are data- or model-driven and (iii)
their (non-)intrusiveness. In the following we consider mostly projection
approaches, as the one adopted in this work; alternatively, one could resort
to different strategies, such as normal form theory or Spectral Submanifolds.
The most recent contributions in this sense include [18] and [19, 20]. In
(ii), for data-driven ROMs we usually refer to ROMs constructed using previous
FOM simulation data (or experimental data, [21]), as opposed to model-driven
methods that rely on some intrinsic properties of the model itself for ROM
construction, such as modal approaches [22, 23, 24, 25]. As for intrusiveness,
we usually denote a ROM as non-intrusive [26] if it can be used with routines
and solvers of commercial FE software and, conversely, as intrusive a method
requiring dedicated routines. Specifically, intrusive methods require access
and manipulation to element-level quantities, as for instance nonlinear
generalized forces and jacobians. Other distinctions can be made in terms of
the types of nonlinearities that a given model can handle and the way
nonlinear functions are evaluated [27]. All these differences ultimately
affect the two phases that all ROMs have in common: the offline phase, in
which the ROM is constructed, and the online phase, in which the simulation
responses are retrieved. As the main goal of ROMs is to reduce computational
effort and time, a key aspect to keep in mind when choosing a method is the
overhead cost to pay in the offline phase; in the case of data-driven methods,
this cost can be as high as the cost associated to the solution of the FOM
[12]. Generally speaking then, data-driven methods (usually based on Proper
Orthogonal Decomposition, or POD, strategies [28]) are used in scenarios where
the high cost associated to the data generation can be amortized: typically,
this is the case of multi-query analysis. In this sense, although not as
versatile and generally applicable as data-driven POD-based approaches, model-
driven strategies in structural dynamics are desirable, for no FOM simulation
is required a priori. Rayleigh-Ritz procedures [29], dual modes [26] and Modal
Derivatives (MDs) [30, 31, 32] are some popular examples.

One way to mitigate the offline overhead costs of all the aforementioned
methods, but especially the data-driven ones, is to resort to (nonlinear)
parametric ROMs, (NL-)pROMs. Also in this context, the literature on linear
systems is quite well developed and consolidated. An extensive survey and
comparison of these methods can be found in [33, 34]. The reduction of
nonlinear parametric Partial Differential Equations (PDEs) is instead still an
active research topic, which has attracted increasing interest in various
disciplines over the years. Interestingly, the vast majority of nonlinear
parametric model order reduction methods is data-driven, POD-based. Some
recent examples include non-intrusive interpolation methods for evaluating
nonlinear functions with hypersurfaces [35, 36] and use of Gaussian Processes
and machine learning for error evaluation and refinement of the pROM [37] or
interpolation on the Grassman manifold via tangent spaces [38]. Alternatively,
many of these methods approximate the nonlinear function using hyper-reduction
methods as the Discrete Empirical Interpolation Method (DEIM) [39, 40] to
speed up the evaluation, and in this sense online basis selection and adaptive
algorithms were studied [41, 42]. However, as mentioned above, POD (and DEIM)
needs a number of FOM simulations to construct the ROM. For this reason, [43]
implemented a Multi-Fidelity strategy in which the parametric dependence was
reconstructed using a large number of low-fidelity models and a minimal number
of high-fidelity evaluations. Other approaches exploit machine learning to
construct an input-output relationship, with convolutional neural networks
[44] and autoencoders [45], which require the training of a network, again,
using preexisting data. Note that most of the above methods lead to pROMs that
are only evaluated in the online phase, i.e. no simulation is actually
performed111By simulation we refer to the solution of a set of equations
describing a system in any kind of analysis setting (e.g. in time or frequency
domain)., but the solutions at the known parameter locations are
“interpolated” to obtain the result.

Although model-driven NL-pROMs seem to be less popular, they offer the
undeniable advantage of being simulation-free, thus considerably cutting down
the offline costs. Interesting recent examples are loosely based on the
extension of methods for linear systems, such as the Non-Linear Moment
Matching (NLMM) scheme [46, 47, 48]. In Ref. [49], a non-parametric ROM is
constructed with NLMM and DEIM for each parameter instance sampled from the
parameter space. These models are then “adjusted” onto a common subspace where
they are interpolated to produce the pROM. This strategy, however, requires
the solution of a set of nonlinear algebraic equations on the FOM at different
time instances, for different signal generators, and at each point on the
parameter grid. For large systems, the computational effort could still be
significant, although lower than that of POD methods.

In this paper we propose a NL-pROM for geometric nonlinearities and
parametrized shape defects to study the behavior of imperfect structures. This
is motivated by the fact that, as it is observed in many engineering
applications, even small imperfections can significantly change
characteristics and performances of a system, as for instance in the case of
MEMS devices [50, 51] and mistuning of gas turbine blades [14]. Other ROMs
have already been developed in this sense [52, 53], but limited to localized
defects. Regarding geometric nonlinearities, we recall that in the case of
continuum finite elements with linear elastic constitutive law and Total
Lagrangian formulation, as in our study, the nonlinear elastic forces are a
polynomials which can be represented using tensors, so that
qualitatively222Due to memory limitations, third and fourth order stiffness
tensors cannot be computed for the FOM, but they can be constructed in reduced
form directly operating at element level [27]. the FOM governing equations
write333$\otimes$ denotes the outer product, : and  $\boldsymbol{\vdots}$ the
double and triple contraction operations.

$\textbf{M}\ddot{\textbf{u}}^{F}+\textbf{C}_{d}\dot{\textbf{u}}^{F}+\tensor[_{2}]{\textbf{K}}{}^{F}\textbf{u}^{F}+\tensor[_{3}]{\textbf{K}}{}^{F}:(\textbf{u}^{F}\otimes\textbf{u}^{F})+\tensor[_{4}]{\textbf{K}}{}^{F}\scalebox{0.75}{
$\boldsymbol{\vdots}$
}(\textbf{u}^{F}\otimes\textbf{u}^{F}\otimes\textbf{u}^{F})=\textbf{f}_{ext}(t)$
(1)

where M, $\textbf{C}_{d}\in\mathbb{R}^{n\times n}$ are the mass and damping
matrices,
$\textbf{u}^{F},\,\dot{\textbf{u}}^{F},\,\ddot{\textbf{u}}^{F}\in\mathbb{R}^{n}$
the displacement, velocity and acceleration vectors, and
$\textbf{f}_{ext}(t)\in\mathbb{R}^{n}$ an external forcing, being $n$ the FOM
number of dofs. $\tensor[_{2}]{\textbf{K}}{}^{F}\in\mathbb{R}^{n\times n}$,
$\tensor[_{3}]{\textbf{K}}{}^{F}\in\mathbb{R}^{n\times n\times n}$ and
$\tensor[_{4}]{\textbf{K}}{}^{F}\in\mathbb{R}^{n\times n\times n\times n}$ are
the stiffness tensors for the linear, quadratic and cubic elastic internal
forces.

Conceptually, the method retrace the one we presented in [1], but it is based
on a different deformation scheme (of which our earlier work resulted to be a
sub-case). An overview of the individual steps of the method is shown in Fig.
1. The user defines as input data the nominal structure (in terms of geometry,
material properties and FE mesh) and a number $m$ of displacement fields
representing the shape defects, which are intended as small deviations from
the nominal geometry (Fig. 1a). These can be known analytically, from
experimental measurements or previous simulations, and finally they can be
discretized with displacement field vectors $\textbf{U}_{i}$ and collected in
a matrix $\textbf{U}=[\textbf{U}_{1},...,\textbf{U}_{m}]$. Each defect can
then be leveraged in amplitude by the parameter vector
$\boldsymbol{\xi}=[\xi_{1},...,\xi_{m}]^{T}$ (Fig. 1b) so that the final
defected geometry represented by the model is given by the global defect
displacement field $\textbf{u}_{d}=\textbf{U}\boldsymbol{\xi}$, i.e. a linear
superposition of the selected defects (Fig. 1c). With this information about
the nominal structure and shape defects, we assemble the RB using a modal
approach with VMs, MDs and Defect Sensitivities (DSs). We then construct the
reduced stiffness tensors, once and for all, projecting the element-level
tensors with the selected RB (Fig. 1d). In this way, linear, quadratic and
cubic elastic forces can be evaluated directly with respect to the reduced
coordinates _and_ shape defect magnitudes without switching between the full
and reduced order space when evaluating the nonlinear function. Our strategy
can then be classified as model-driven (simulation-free). Finally, in the
online phase, the simulation is performed with the reduced governing equations
(Fig. 1e). Notice that the model is used to run a simulation, not to evaluate
a solution as in interpolation-like techniques: as such, different forcing
terms and also different analysis types (e.g. transient, frequency response)
are possible.

Figure 1: Overview of the proposed method, schematically illustrated for a
pinned beam.

All of this is possible thanks to the modified definition of the Green-
Lagrange strain tensor we use. Specifically, our strain tensor embeds two
subsequent transformations: (i) the one from nominal to defected geometry
(which, at the end, will be parametrized), and (ii) the one from the defected
configuration to the deformed/final one. The deformation produced by the
latter is the one we measure, so no strain/stresses are introduced by the
presence of the defect in (i); however, the deformation of (ii) will depend on
(i). The formulation we obtain however contains rational terms which cannot be
used for a tensorial representation (which can describe polynomials only).
Given the assumed small entity of the shape defects, we advocate the use of a
Neumann expansion to approximate the Green-Lagrange tensor, obtaining again a
polynomial form. Applying standard FE procedures, we finally get to the
expression of the reduced elastic internal forces, which will parametrically
depend on the defect amplitudes $\boldsymbol{\xi}$. In this framework, we show
that the model in [1] (whose deformation formulation was based on [54])
corresponds to a lower order Neumann expansion with integrals evaluated on the
nominal volume, and that the higher order approximation we propose here leads
to better accuracy and to a larger applicability range.

The work is organized as follows: the modified strain formulation is given in
Section 2 and approximated using Neumann expansion in Section 3. In Section 4
the FE discretization is developed and then used in Section 5 to construct the
reduced order stiffness tensors. The choice and computation of the RB is
described in Section 6. Finally, numerical studies in Sections 7 and 8
demonstrate the effectiveness of the proposed approach on a 2D FE clamped beam
and on a MEMS gyroscope and computational times are discussed.

##  2 Strain formulation: a two-steps deformation approach

Figure 2: Scheme for the considered deformation setting. A nominal structure,
of coordinates $\textbf{x}_{0}$, undergoes a deformation described by the
transformation map $\Phi_{1}$. The structure is now in the deformed
configuration (coordinates $\textbf{x}_{d}$). A second transformation
$\Phi_{2}$ and the displacement u describe the deformation from the defected
configuration to the final one.

Let us consider the scheme depicted in Fig. 2. A nominal body of coordinates
$\textbf{x}_{0}=\\{x_{0},\,y_{0},\,z_{0}\\}$ undergoes a first deformation
described by the map $\boldsymbol{\Phi}_{1}$, which brings the body in a new
configuration with coordinates
$\textbf{x}_{d}=\\{x_{d},\,y_{d},\,z_{d}\\}=\boldsymbol{\Phi}_{1}(\textbf{x}_{0})$.
The displacement corresponding to this operation is
$\textbf{u}_{d}=\\{u_{d},\,v_{d},\,w_{d}\\}=\textbf{x}_{d}-\textbf{x}_{0}$. We
will refer to this second configuration as the defected configuration. As it
will be detailed later, in our method $\textbf{u}_{d}$ will be a user-defined
displacement field representing a shape defect which, superimposed to the
nominal geometry, defines the configuration with respect to which we measure
deformation. Let us now consider a second deformation, described by the map
$\boldsymbol{\Phi}_{2}$, from the defected configuration to the final one,
with coordinates
$\textbf{x}(t)=\\{x(t),\,y(t),\,z(t)\\}=\boldsymbol{\Phi}_{2}(\textbf{x}_{d},t)$.
We will refer to the latter as to the deformed or final configuration, whose
displacement is given by
$\textbf{u}=\\{u,\,v,\,w\\}=\textbf{x}-\textbf{x}_{d}$.

Considering the infinitesimal line segment $\text{d}\textbf{x}_{0}$ in the
nominal geometry, we can define the line segments d$\textbf{x}_{d}$ and dx in
the defected and deformed configurations as

$\displaystyle\text{d}\textbf{x}_{d}=\textbf{F}_{1}\text{d}\textbf{x}_{0},$
(2a)
$\displaystyle\text{d}\textbf{x}=\textbf{F}_{2}\text{d}\textbf{x}_{d}=\textbf{F}_{2}\textbf{F}_{1}\text{d}\textbf{x}_{0},$
(2b)

where the deformation gradients $\textbf{F}_{1}$ and $\textbf{F}_{2}$ are
given by

$\displaystyle\textbf{F}_{1}=\nabla_{0}\textbf{x}_{d}=\frac{\partial\textbf{x}_{d}}{\partial\textbf{x}_{0}}=\textbf{I}+\nabla_{0}\textbf{u}_{d}=\textbf{I}+\frac{\partial\textbf{u}_{d}}{\partial\textbf{x}_{0}}=\textbf{I}+\textbf{D}_{d},$
(3a)
$\displaystyle\textbf{F}_{2}=\nabla_{d}\textbf{x}=\frac{\partial\textbf{x}}{\partial\textbf{x}_{d}}=\textbf{I}+\nabla_{d}\textbf{u}=\textbf{I}+\frac{\partial\textbf{u}}{\partial\textbf{x}_{d}}=\textbf{I}+\textbf{D}_{2}.$
(3b)

and where $\textbf{D}_{d}$ and $\textbf{D}_{2}$ are the displacement
derivative matrices of the first and second transformations, respectively.
Using the chain rule, we can also define

$\textbf{D}=\frac{\partial\textbf{u}}{\partial\textbf{x}_{0}}=\frac{\partial\textbf{u}}{\partial\textbf{x}_{d}}\frac{\partial\textbf{x}_{d}}{\partial\textbf{x}_{0}}=\textbf{D}_{2}\textbf{F}_{1},$
(4)

so that $\textbf{D}_{2}=\textbf{D}\textbf{F}_{1}^{-1}$ can be referred to the
nominal coordinates.

Using Eqs. (2)–(4), the stretch between deformed and defected configurations
writes

$\begin{split}\mathcal{S}=\text{d}\textbf{x}^{T}\text{d}\textbf{x}-\text{d}\textbf{x}_{d}^{T}\text{d}\textbf{x}_{d}&=\text{d}\textbf{x}_{0}^{T}\textbf{F}_{1}^{T}(\textbf{F}_{2}^{T}\textbf{F}_{2}-\textbf{I})\textbf{F}_{1}\text{d}\textbf{x}_{0}\\\
&=\text{d}\textbf{x}_{0}^{T}(\textbf{D}+\textbf{D}^{T}+\textbf{D}^{T}\textbf{D}+\textbf{D}_{d}^{T}\textbf{D}+\textbf{D}^{T}\textbf{D}_{d})\text{d}\textbf{x}_{0}.\end{split}$
(5)

Measuring the deformation with respect to the defected configuration, the
second order Green-Lagrange strain tensor $\textbf{E}_{2}$ is defined as

$\mathcal{S}=2\text{d}\textbf{x}^{T}_{d}\textbf{E}_{2}\text{d}\textbf{x}_{d}=2\text{d}\textbf{x}^{T}_{0}\textbf{F}_{1}^{T}\textbf{E}_{2}\textbf{F}_{1}\text{d}\textbf{x}_{0},$
(6)

which, rearranged, leads to

$\textbf{E}_{2}=\frac{1}{2}\textbf{F}_{1}^{-T}(\textbf{D}+\textbf{D}^{T}+\textbf{D}^{T}\textbf{D}+\textbf{D}_{d}^{T}\textbf{D}+\textbf{D}^{T}\textbf{D}_{d})\textbf{F}_{1}^{-1}.$
(7)

Looking at Eqs. (5) and (7), it can be easily verified that $\textbf{E}_{2}$
correctly satisfies the minimum requirements for a strain measure to vanish
under a rigid body translation ($\textbf{F}_{2}=\textbf{I}$) and/or rotation
($\textbf{F}_{2}^{T}\textbf{F}_{2}=\textbf{R}^{T}\textbf{R}=\textbf{I}$, with
R orthonormal rotation matrix), for any $\textbf{F}_{1}$. Eq. (7) is indeed an
exact expression for the strains from defected to final configuration. Notice
however that in this form all the quantities are computed with respect to the
nominal coordinates $\textbf{x}_{0}$.

##  3 Strain approximations

The introduced strain measure, being referred to the nominal geometry only,
paves the way for the pre-computation of the stiffness tensors, as it will be
shown in the following sections. However, as mentioned in the introduction, a
tensorial formulation can be applied only when the internal forces display a
polynomial dependence on the displacements, which in the present case include
both $\textbf{u}_{d}$ and u. The inverse of the deformation gradient
$\textbf{F}_{1}$ in Eq. (7) entails a rational dependence on $\textbf{u}_{d}$,
and therefore needs some attention. Let us consider the following known
result:

#### Neumann expansion

If P is a square matrix and the Neumann series
$\sum_{n=0}^{+\infty}\textbf{P}^{n}$ is convergent, we have that

$(\textbf{I}-\textbf{P})^{-1}=\sum_{n=0}^{+\infty}\textbf{P}^{n}$ (8)

A spectral norm444The spectral norm of a matrix $\textbf{A}_{1}$ is defined as
the square root of the largest singular value of
$\textbf{A}_{1}^{*}\textbf{A}_{1}$, being $\textbf{A}_{1}^{*}$ the conjugate
transpose of $\textbf{A}_{1}$, that is:
$\left\lVert\textbf{A}_{1}\right\rVert_{2}=\sqrt{\lambda_{max}(\textbf{A}_{1}^{*}\textbf{A}_{1})}$
$\varepsilon=\left\lVert\textbf{P}\right\rVert_{2}<1$ is a sufficient
condition for the convergence of the Neumann series. Moreover, it can be shown
[55] that truncating the sum to order $N$ the norm error is bounded as

$\delta_{N}=\left\lVert\left(\textbf{I}-\textbf{P}\right)^{-1}-\sum_{n=0}^{N}\textbf{P}^{n}\right\rVert_{2}\leq\dfrac{\varepsilon^{N+1}}{1-\varepsilon}=\delta_{lim}.$
(9)

Letting $\textbf{P}=-\textbf{D}_{d}$, we can expand $\textbf{F}_{1}^{-1}$
using the Neumann series as

$\textbf{F}_{1}^{-1}=(\textbf{I}+\textbf{D}_{d})^{-1}\approx\sum_{n=0}^{N}(-\textbf{D}_{d})^{n}.$
(10)

Under the assumption of small defects (i.e.
$\left\lVert\textbf{D}_{d}\right\rVert_{2}\ll 1$) the series is guaranteed to
converge. Moreover, we can truncate the sum in Eq. (10) to $N=1$, obtaining:

$\textbf{E}_{2,N}=\frac{1}{2}(\textbf{I}-\textbf{D}_{d})^{T}(\textbf{D}+\textbf{D}^{T}+\textbf{D}^{T}\textbf{D}+\textbf{D}_{d}^{T}\textbf{D}+\textbf{D}^{T}\textbf{D}_{d})(\textbf{I}-\textbf{D}_{d}),$
(11)

which, solving the product, can be rewritten as:

$\textbf{E}_{2,N}=\frac{1}{2}\Big{[}\left(\textbf{D}+\textbf{D}^{T}+\textbf{D}^{T}\textbf{D}+\cancel{\textbf{D}_{d}^{T}\textbf{D}}+\cancel{\textbf{D}^{T}\textbf{D}_{d}}\right)+\\\
-\left(\cancel{\textbf{D}_{d}^{T}\textbf{D}}+\textbf{D}_{d}^{T}\textbf{D}^{T}+\textbf{D}_{d}^{T}\textbf{D}^{T}\textbf{D}+\textbf{D}_{d}^{T}\textbf{D}_{d}^{T}\textbf{D}+\cancel{\textbf{D}_{d}^{T}\textbf{D}^{T}\textbf{D}_{d}}\right)+\\\
-\left(\textbf{D}\textbf{D}_{d}+\cancel{\textbf{D}^{T}\textbf{D}_{d}}+\textbf{D}^{T}\textbf{D}\textbf{D}_{d}+\cancel{\textbf{D}_{d}^{T}\textbf{D}\textbf{D}_{d}}+\textbf{D}^{T}\textbf{D}_{d}\textbf{D}_{d}\right)+\\\
+\left(\cancel{\textbf{D}_{d}^{T}\textbf{D}\textbf{D}_{d}}+\cancel{\textbf{D}_{d}^{T}\textbf{D}^{T}\textbf{D}_{d}}+\textbf{D}_{d}^{T}\textbf{D}^{T}\textbf{D}\textbf{D}_{d}+\textbf{D}_{d}^{T}\textbf{D}_{d}^{T}\textbf{D}\textbf{D}_{d}+\textbf{D}_{d}^{T}\textbf{D}^{T}\textbf{D}_{d}\textbf{D}_{d}\right)\Big{]}.$
(12)

Finally, neglecting the terms $\mathcal{O}(\textbf{D}_{d}^{2})$, i.e. assuming
that the first transformation $\Phi_{1}$ is linear, Eq. (12) reduces:

$\textbf{E}_{2,N1}=\frac{1}{2}\left(\textbf{D}+\textbf{D}^{T}+\textbf{D}^{T}\textbf{D}-\textbf{D}_{d}^{T}\textbf{D}^{T}-\textbf{D}\textbf{D}_{d}-\textbf{D}_{d}^{T}\textbf{D}^{T}\textbf{D}-\textbf{D}^{T}\textbf{D}\textbf{D}_{d}\right)$
(13)

The modified Green-Lagrange strain tensor $\textbf{E}_{2,N1}$ is a polynomial
function of the derivatives of the displacement fields u and $\textbf{u}_{d}$,
and can be thus used to compute a ROM using tensors.

Remark 1 (on Budiansky approximation). The strain formulation in [54], used by
Budiansky to study buckling in presence of defects, was obtained by
subtracting the strain that a defect would produce on the nominal structure
from the strain of the deformed structure measured with respect to the nominal
configuration. It can be shown that truncating the Neummann series to the
zero-th order (i.e. setting $N=0$, so that $\textbf{F}_{1}^{-1}=\textbf{I}$)
and using Eq. (7) and (10), the strain writes:

$\textbf{E}_{2,N0}=\frac{1}{2}\left(\textbf{D}+\textbf{D}^{T}+\textbf{D}^{T}\textbf{D}+\textbf{D}_{d}^{T}\textbf{D}+\textbf{D}^{T}\textbf{D}_{d}\right),$
(14)

which is the same strain tensor we adopted in [1] following Budiansky’s
approximation.

##  4 Finite Element formulation

In this section we derive the elastic internal forces (at element-level) for
the FE discretization of the full order model based on the strain as defined
in Eq. (13). We remark that this full model represents just an approximation
of the reference full order model FOM-d (where the defect is embedded directly
in the mesh by shifting the position of the nodes). Although not offering any
direct advantage over FOM-d, this full model will allow us to compute the
parametric ROM, as it will be explained in Section 5.

First, it is convenient to switch to Voigt notation. Let
$\boldsymbol{\theta}=\\{u_{,x}\,u_{,y}\,u_{,z}\,v_{,x}\,v_{,y}\,v_{,z}\,w_{,x}\,w_{,y}\,w_{,z}\\}^{T}$
be the vectorized form of D and
$\boldsymbol{\theta}_{d}=\\{u_{d,x}\,u_{d,y}\,u_{d,z}\,v_{d,x}\,v_{d,y}\,v_{d,z}\,w_{d,x}\,w_{d,y}\,w_{d,z}\\}^{T}$
be the vectorized form of $\textbf{D}_{d}$555For the derivatives of the
displacement components, we use the notation $u_{,x}=\frac{\partial
u}{\partial x_{0}}$ and $u_{d,x}=\frac{\partial u_{d}}{\partial x_{0}}$
(similar definition for $v$, $w$ and the other spatial coordinates).. Eq. (13)
rewrites:

$\textbf{E}_{1,N1}=\left(\textbf{H}+\frac{1}{2}\textbf{A}_{1}(\boldsymbol{\theta})+\textbf{A}_{2}(\boldsymbol{\theta}_{d})+\textbf{A}_{3}(\boldsymbol{\theta}_{d})\textbf{A}_{1}(\boldsymbol{\theta})\right)\boldsymbol{\theta}$
(15)

where

$\textbf{H}=\small\begin{bmatrix}1&0&0&0&0&0&0&0&0\\\ 0&0&0&0&1&0&0&0&0\\\
0&0&0&0&0&0&0&0&1\\\ 0&1&0&1&0&0&0&0&0\\\ 0&0&1&0&0&0&1&0&0\\\
0&0&0&0&0&1&0&1&0\\\ \end{bmatrix},\normalsize\hskip
5.69046pt\textbf{A}_{1}(\boldsymbol{\theta})=\small\begin{bmatrix}u_{,x}&0&0&v_{,x}&0&0&w_{,x}&0&0\\\
0&u_{,y}&0&0&v_{,y}&0&0&w_{,y}&0\\\ 0&0&u_{,z}&0&0&v_{,z}&0&0&w_{,z}\\\
u_{,y}&u_{,x}&0&v_{,y}&v_{,x}&0&w_{,y}&w_{,x}&0\\\
u_{,z}&0&u_{,x}&v_{,z}&0&v_{,x}&w_{,z}&0&w_{,x}\\\
0&u_{,z}&u_{,y}&0&v_{,z}&v_{,y}&0&w_{,z}&w_{,y}\\\ \end{bmatrix},\\\
\normalsize\hskip 5.69046pt\text{s.t.: }\hskip
5.69046pt\textbf{D}+\textbf{D}^{T}+\textbf{D}^{T}\textbf{D}\longleftrightarrow(2\textbf{H}+\textbf{A}_{1}(\boldsymbol{\theta}))\boldsymbol{\theta},$
(16)
$\textbf{A}_{2}(\boldsymbol{\theta}_{d})=(-1)\small\begin{bmatrix}u_{d,x}&v_{d,x}&w_{d,x}&0&0&0&0&0&0\\\
0&0&0&u_{d,y}&v_{d,y}&w_{d,y}&0&0&0\\\ 0&0&0&0&0&0&u_{d,z}&v_{d,z}&w_{d,z}\\\
u_{d,y}&v_{d,y}&w_{d,y}&u_{d,x}&v_{d,x}&w_{d,x}&0&0&0\\\
u_{d,z}&v_{d,z}&w_{d,z}&0&0&0&u_{d,x}&v_{d,x}&w_{d,x}\\\
0&0&0&u_{d,z}&v_{d,z}&w_{d,z}&u_{d,y}&v_{d,y}&w_{d,y}\\\ \end{bmatrix},\\\
\normalsize\hskip 5.69046pt\text{s.t.: }\hskip
5.69046pt-\textbf{D}^{T}_{d}\textbf{D}^{T}-\textbf{D}\textbf{D}_{d}\longleftrightarrow
2\textbf{A}_{2}(\boldsymbol{\theta}_{d})\boldsymbol{\theta},$ (17)
$\textbf{A}_{3}(\boldsymbol{\theta}_{d})=\left(-\frac{1}{2}\right)\small\begin{bmatrix}2u_{d,x}&0&0&v_{d,x}&w_{d,x}&0\\\
0&2v_{d,y}&0&u_{d,y}&0&w_{d,y}\\\ 0&0&2w_{d,z}&0&u_{d,z}&v_{d,z}\\\
2u_{d,y}&2v_{d,x}&0&u_{d,x}+v_{d,y}&w_{d,y}&w_{d,x}\\\
2u_{d,z}&0&2w_{d,x}&v_{d,z}&u_{d,x}+w_{d,z}&v_{d,x}\\\
0&2v_{d,z}&2w_{d,y}&u_{d,z}&u_{d,y}&v_{d,y}+w_{d,z}\end{bmatrix},\normalsize\\\
\hskip 5.69046pt\text{s.t.: }\hskip
5.69046pt-\textbf{D}_{d}^{T}\textbf{D}^{T}\textbf{D}-\textbf{D}^{T}\textbf{D}\textbf{D}_{d}\longleftrightarrow
2\textbf{A}_{3}(\boldsymbol{\theta}_{d})\textbf{A}_{1}(\boldsymbol{\theta})\boldsymbol{\theta}.$
(18)

Let us now define $\textbf{u}^{e}$ and $\textbf{u}_{d}^{e}$ as the nodal
displacement vectors of a (continuum) finite element. Calling G the shape
function derivatives matrix, such that
$\boldsymbol{\theta}=\textbf{G}\textbf{u}^{e}$ and
$\boldsymbol{\theta}_{d}=\textbf{G}\textbf{u}_{d}^{e}$, and exploiting the
property by which
$\textbf{A}_{1}(\boldsymbol{\theta})\delta\boldsymbol{\theta}=\textbf{A}_{1}(\delta\boldsymbol{\theta})\boldsymbol{\theta}$,
the virtual variation of the strain in Eq. (15) writes

$\delta\textbf{E}_{1,N1}=\left(\textbf{H}+\textbf{A}_{1}+\textbf{A}_{2}+2\textbf{A}_{3}\textbf{A}_{1}\right)\textbf{G}\delta\textbf{u}^{e}=\textbf{B}_{N1}\delta\textbf{u}^{e},$
(19)

where B is the strain-displacements matrix and where we dropped the explicit
dependencies on $\boldsymbol{\theta}_{d}$ and $\boldsymbol{\theta}$ to ease
the notation. The virtual work of internal forces is given by

$W_{int}=\int_{V_{d}}\delta\textbf{E}_{1,N1}^{T}\textbf{S}_{1}\text{
d}V_{d}=(\delta\textbf{u}^{e})^{T}\int_{V_{d}}\textbf{B}_{N1}^{T}\textbf{S}_{1}\text{
d}V_{d},$ (20)

where $\textbf{S}_{1}=\textbf{C}\textbf{E}_{1}$ is the Piola-Kirchhoff stress
in Voigt notation, being C the linear elastic constitutive matrix, and where
$V_{d}$ is the volume of the defected configuration. The expression for the
internal forces $\textbf{f}_{int}$ follows from the virtual work:

$\textbf{f}_{int}=\int_{V_{d}}\textbf{B}_{N1}^{T}\textbf{C}\textbf{E}_{1,N1}\text{
d}V_{d}.$ (21)

Finally, the tangent stiffness matrix can be computed as usual taking the
virtual variation of the internal forces (see B). Equations (15) and (21) can
be used to perform tests and/or simulations of the full model and to compare
the results to the corresponding FOM-d in order to assess the quality of the
approximation before the reduction of the model. In the next section, the
DpROM derived from this formulation is presented.

##  5 Tensorial representation of internal forces

###  5.1 Element–level tensors

Under the hypothesis that for small defects $V_{d}\approx V_{o}$, Eq. (21) in
full can be written as

$\textbf{f}_{int}=\int_{V_{o}}\textbf{G}^{T}\left(\textbf{H}+\textbf{A}_{1}+\textbf{A}_{2}+2\textbf{A}_{3}\textbf{A}_{1}\right)^{T}\textbf{C}\left(\textbf{H}+\frac{1}{2}\textbf{A}_{1}+\textbf{A}_{2}+\textbf{A}_{3}\textbf{A}_{1}\right)\textbf{G}\textbf{u}^{e}\,\text{d}{V_{o}}.$
(22)

To “extract” the displacement vectors u and $\textbf{u}_{d}$ from matrices
$\textbf{A}_{1}$, $\textbf{A}_{2}$ and $\textbf{A}_{3}$, we can write:

$\displaystyle\textbf{A}_{1}$
$\displaystyle={\textbf{L}_{1}}\cdot\boldsymbol{\theta}={\textbf{L}_{1}}\cdot(\textbf{G}\textbf{u}^{e}),$
(23a) $\displaystyle\textbf{A}_{2}$
$\displaystyle=\textbf{L}_{2}\cdot\boldsymbol{\theta}_{d}=\textbf{L}_{2}\cdot(\textbf{G}\textbf{u}_{d}^{e}),$
(23b) $\displaystyle\textbf{A}_{3}\textbf{A}_{1}$
$\displaystyle=(\textbf{L}_{3}\cdot\boldsymbol{\theta}_{d})\cdot\boldsymbol{\theta}=(\textbf{L}_{3}\cdot(\textbf{G}\textbf{u}_{d}^{e}))\cdot(\textbf{G}\textbf{u}^{e}),$
(23c)

where ${\textbf{L}_{1}},\textbf{L}_{2}\in\mathbb{R}^{6\times 9\times 9}$ and
$\textbf{L}_{3}\in\mathbb{R}^{6\times 9\times 9\times 9}$ are constant sparse
matrices (see A).

Dropping for convenience the integral operation over the volume $V_{o}$
(implicitly assumed), we can separate the contributions in Eq. (22) as

$\displaystyle\textbf{f}_{1}$
$\displaystyle=\textbf{G}^{T}\left(\textbf{H}^{T}\textbf{C}\textbf{H}+\textbf{H}^{T}\textbf{C}\textbf{A}_{2}+\textbf{A}_{2}^{T}\textbf{C}\textbf{H}+\textbf{A}_{2}^{T}\textbf{C}\textbf{A}_{2}\right)\textbf{G}\textbf{u}^{e},$
(24a)
$\displaystyle\begin{split}\textbf{f}_{2}&=\textbf{G}^{T}\left(\frac{1}{2}\textbf{H}^{T}\textbf{C}\textbf{A}_{1}+\textbf{A}_{1}^{T}\textbf{C}\textbf{H}+\frac{1}{2}\textbf{A}_{2}^{T}\textbf{C}\textbf{A}_{1}+\textbf{A}_{1}^{T}\textbf{C}\textbf{A}_{2}+\right.\\\
&\hskip
38.69576pt2\textbf{A}_{1}^{T}\textbf{A}_{3}^{T}\textbf{C}\textbf{H}+\textbf{H}^{T}\textbf{C}\textbf{A}_{3}\textbf{A}_{1}+2\textbf{A}_{1}^{T}\textbf{A}_{3}^{T}\textbf{C}\textbf{A}_{2}+\textbf{A}_{2}^{T}\textbf{C}\textbf{A}_{3}\textbf{A}_{1}\Big{)}\textbf{G}\textbf{u}^{e},\end{split}$
(24b) $\displaystyle\textbf{f}_{3}$
$\displaystyle=\textbf{G}^{T}\left(\frac{1}{2}\textbf{A}_{1}^{T}\textbf{C}\textbf{A}_{1}+\textbf{A}_{1}^{T}\textbf{C}\textbf{A}_{3}\textbf{A}_{1}+\textbf{A}_{1}^{T}\textbf{A}_{3}^{T}\textbf{C}\textbf{A}_{1}+2\textbf{A}_{1}^{T}\textbf{A}_{3}^{T}\textbf{C}\textbf{A}_{3}\textbf{A}_{1}\right)\textbf{G}\textbf{u}^{e},$
(24c)

where $\textbf{f}_{1}$, $\textbf{f}_{2}$ and $\textbf{f}_{3}$ are the linear,
quadratic and cubic terms in the displacement u, respectively. These can be
recasted in tensorial form as

$\displaystyle\textbf{f}_{1}$
$\displaystyle=\left[\tensor[_{2n}]{\textbf{K}}{}+\tensor[_{3d}]{\textbf{K}}{}\cdot\textbf{u}_{d}^{e}+\tensor[_{4dd}]{\textbf{K}}{}:(\textbf{u}_{d}^{e}\otimes\textbf{u}_{d}^{e})\right]\cdot\textbf{u}=\tensor[_{2}]{\textbf{K}}{}(\textbf{u}_{d}^{e})\cdot\textbf{u}^{e},$
(25a) $\displaystyle\textbf{f}_{2}$
$\displaystyle=\left[\tensor[_{3n}]{\textbf{K}}{}+\tensor[_{4d}]{\textbf{K}}{}\cdot\textbf{u}_{d}^{e}+\tensor[_{5dd}]{\textbf{K}}{}:(\textbf{u}_{d}^{e}\otimes\textbf{u}_{d}^{e})\right]:(\textbf{u}\otimes\textbf{u})=\tensor[_{3}]{\textbf{K}}{}(\textbf{u}_{d}^{e}):(\textbf{u}^{e}\otimes\textbf{u}^{e}),$
(25b) $\displaystyle\textbf{f}_{3}$
$\displaystyle=\left[\tensor[_{4n}]{\textbf{K}}{}+\tensor[_{5d}]{\textbf{K}}{}\cdot\textbf{u}_{d}^{e}+\tensor[_{6dd}]{\textbf{K}}{}:(\textbf{u}_{d}^{e}\otimes\textbf{u}_{d}^{e})\right]\scalebox{0.75}{
$\boldsymbol{\vdots}$
}(\textbf{u}\otimes\textbf{u}\otimes\textbf{u})=\tensor[_{4}]{\textbf{K}}{}(\textbf{u}_{d}^{e})\scalebox{0.75}{
$\boldsymbol{\vdots}$
}(\textbf{u}^{e}\otimes\textbf{u}^{e}\otimes\textbf{u}^{e}).$ (25c)

###  5.2 Reduced tensors and internal forces

We now derive the reduced internal forces and tensors via Galerkin projection.
Let us assume the following reduction for the displacement vectors
$\textbf{u}^{e}$ and $\textbf{u}_{d}^{e}$:

$\displaystyle\textbf{u}^{e}$
$\displaystyle\approx\textbf{V}^{e}\boldsymbol{\eta},\text{ with
}\textbf{V}^{e}\in\mathbb{R}^{n_{e}\times
m},\,\boldsymbol{\eta}\in\mathbb{R}^{m},$ (26a)
$\displaystyle\textbf{u}_{d}^{e}$
$\displaystyle=\textbf{U}^{e}\boldsymbol{\xi},\text{ with
}\textbf{U}^{e}\in\mathbb{R}^{n_{e}\times
m_{d}},\,\boldsymbol{\xi}\in\mathbb{R}^{m_{d}}$ (26b)

being $\textbf{V}^{e}$ and $\textbf{U}^{e}$ the partitions of the RBs
($\textbf{V}\in\mathbb{R}^{n\times m}$ and $\textbf{U}\in\mathbb{R}^{n\times
m_{d}}$, with $n$ number of dofs of the full order model) relative to the
element, $\boldsymbol{\eta}$ and $\boldsymbol{\xi}$ the reduced coordinates.
$m\ll n$ is thus the number of vectors included in the RB V while $m_{d}$
represents the number of the assumed shape defects, collected column-wise in
U. Introducing $\Gamma=\textbf{G}\textbf{V}^{e}$ and
$\Upsilon=\textbf{G}\textbf{U}^{e}$, we can directly define the reduced order
tensors using Einstein’s notation as

$\displaystyle\tensor[_{2n}]{Q}{{}_{IJ}}$
$\displaystyle=\Gamma_{iI}H_{ji}C_{jk}H_{kl}\Gamma_{kJ},$ (27a)
$\displaystyle\tensor[_{3d}]{Q}{{}_{IJ\underline{K}}}$
$\displaystyle=\Gamma_{iI}\left(H_{ji}C_{jk}{L_{2}}_{kla}\Upsilon_{aK}+{L_{2}}_{jia}\Upsilon_{aK}C_{jk}H_{kl}\right)\Gamma_{lJ},$
(27b) $\displaystyle\tensor[_{4dd}]{Q}{{}_{IJ\underline{KL}}}$
$\displaystyle=\Gamma_{iI}{L_{2}}_{jia}\Upsilon_{aK}C_{jk}{L_{2}}_{klb}\Upsilon_{bL}\Gamma_{lJ},$
(27c) $\displaystyle\tensor[_{3n}]{Q}{{}_{IJK}}$
$\displaystyle=\Gamma_{iI}\left(\frac{1}{2}H_{ji}C_{jk}{L_{1}}_{kla}\Gamma_{aK}+{L_{1}}_{jia}\Gamma_{aK}C_{jk}H_{kl}\right)\Gamma_{lJ},$
(27d)
$\displaystyle\begin{split}\tensor[_{4d}]{Q}{{}_{IJK\underline{L}}}&=\Gamma_{iI}\left(\frac{1}{2}{L_{2}}_{jia}\Upsilon_{aL}C_{jk}{L_{1}}_{klb}\Gamma_{bK}+{L_{1}}_{jia}\Gamma_{aK}C_{jk}{L_{2}}_{klb}\Upsilon_{bL}\right.+\\\
&+2{L_{3}}_{jiab}\Upsilon_{bL}\Gamma_{aK}C_{jk}H_{kl}+H_{ji}C_{jk}{L_{3}}_{klab}\Upsilon_{bL}\Gamma_{aK}\Big{)}\Gamma_{lJ},\end{split}$
(27e) $\displaystyle\tensor[_{5dd}]{Q}{{}_{IJK\underline{LM}}}$
$\displaystyle=\Gamma_{iI}\left(2{L_{3}}_{jiab}\Upsilon_{bL}\Gamma_{aK}C_{jk}{L_{2}}_{klc}\Upsilon_{cM}+{L_{2}}_{jia}\Upsilon_{aL}C_{jk}{L_{3}}_{klbc}\Upsilon_{cM}\Gamma_{bK}\right)\Gamma_{lJ},$
(27f) $\displaystyle\tensor[_{4n}]{Q}{{}_{IJKL}}$
$\displaystyle=\frac{1}{2}\,\Gamma_{iI}{L_{1}}_{jia}\Gamma_{aK}C_{jk}{L_{1}}_{klb}\Gamma_{bL}\Gamma_{lJ},$
(27g) $\displaystyle\tensor[_{5d}]{Q}{{}_{IJKL\underline{M}}}$
$\displaystyle=\Gamma_{iI}\left({L_{1}}_{jia}\Gamma_{aK}C_{jk}{L_{3}}_{klbc}\Upsilon_{cM}\Gamma_{bL}+{L_{3}}_{jiab}\Upsilon_{bM}\Gamma_{aK}C_{jk}{L_{1}}_{klc}\Gamma_{lL}\right)\Gamma_{lJ},$
(27h) $\displaystyle\tensor[_{6dd}]{Q}{{}_{IJKL\underline{MN}}}$
$\displaystyle=2\,\Gamma_{iI}{L_{3}}_{jiab}\Upsilon_{bM}\Gamma_{aL}C_{jk}{L_{3}}_{klcd}\Upsilon_{dN}\Gamma_{cK}\Gamma_{lJ}.$
(27i)

where, for convenience, tensor dimensions of size $m$ are denoted by capital
letter subscripts, dimensions of size $m_{d}$ by underlined capital letter
ones. So, for example, $\tensor[_{3d}]{\textbf{Q}}{}\in\mathbb{R}^{m\times
m\times m_{d}}$. The global reduced tensors of the full structure can then be
computed directly summing up the element contributions, a procedure which is
highly parallelizable. Reduced (global) internal forces can therefore be
defined as
$\textbf{f}_{int,r}=\textbf{f}_{1r}+\textbf{f}_{2r}+\textbf{f}_{3r}$, where

$\displaystyle\textbf{f}_{1r}$
$\displaystyle=\left[\tensor[_{2n}]{\textbf{Q}}{}+\tensor[_{3d}]{\textbf{Q}}{}\cdot\boldsymbol{\xi}+\tensor[_{4dd}]{\textbf{Q}}{}:(\boldsymbol{\xi}\otimes\boldsymbol{\xi})\right]\cdot\boldsymbol{\eta}=\tensor[_{2}]{\textbf{Q}}{}(\boldsymbol{\xi})\cdot\boldsymbol{\eta},$
(28a) $\displaystyle\textbf{f}_{2r}$
$\displaystyle=\left[\tensor[_{3n}]{\textbf{Q}}{}+\tensor[_{4d}]{\textbf{Q}}{}\cdot\boldsymbol{\xi}+\tensor[_{5dd}]{\textbf{Q}}{}:(\boldsymbol{\xi}\otimes\boldsymbol{\xi})\right]:(\boldsymbol{\eta}\otimes\boldsymbol{\eta})=\tensor[_{3}]{\textbf{Q}}{}(\boldsymbol{\xi}):(\boldsymbol{\eta}\otimes\boldsymbol{\eta}),$
(28b) $\displaystyle\textbf{f}_{3r}$
$\displaystyle=\left[\tensor[_{4n}]{\textbf{Q}}{}+\tensor[_{5d}]{\textbf{Q}}{}\cdot\boldsymbol{\xi}+\tensor[_{6dd}]{\textbf{Q}}{}:(\boldsymbol{\xi}\otimes\boldsymbol{\xi})\right]\scalebox{0.75}{
$\boldsymbol{\vdots}$
}(\boldsymbol{\eta}\otimes\boldsymbol{\eta}\otimes\boldsymbol{\eta})=\tensor[_{4}]{\textbf{Q}}{}(\boldsymbol{\xi})\scalebox{0.75}{
$\boldsymbol{\vdots}$
}(\boldsymbol{\eta}\otimes\boldsymbol{\eta}\otimes\boldsymbol{\eta}),$ (28c)

and the reduced tangent stiffness matrix $\tensor[_{t}]{\textbf{Q}}{}$ can be
written simply as

$\tensor[_{t}]{Q}{{}_{IJ}}=\tensor[_{2}]{Q}{{}_{IJ}}+(\tensor[_{3}]{Q}{{}_{IJj}}+\tensor[_{3}]{Q}{{}_{IjJ}})\eta_{j}+(\tensor[_{4}]{Q}{{}_{IJij}}+\tensor[_{4}]{Q}{{}_{IiJj}}+\tensor[_{4}]{Q}{{}_{IijJ}})\eta_{i}\eta_{j}.$
(29)

Finally, the equations of motion for the reduced system write:

$\textbf{M}_{r}\ddot{\boldsymbol{\eta}}(t)+\textbf{C}_{r}\dot{\boldsymbol{\eta}}(t)+\textbf{f}_{int,r}(\boldsymbol{\eta}(t),\boldsymbol{\xi})=\textbf{f}_{ext,r}(t)$
(30)

where $\textbf{M}_{r}=\textbf{V}^{T}\textbf{M}\textbf{V}$ and
$\textbf{C}_{r}=\textbf{V}^{T}\textbf{C}_{d}\textbf{V}$ are the reduced mass
and damping matrices,
$\textbf{f}_{ext,r}(t)=\textbf{V}^{T}\textbf{f}_{ext}(t)$ the reduced external
forces acting on the system. Notice that, in accordance to the hypothesis made
for the internal forces, also these matrices must be integrated over the
nominal volume $V_{o}$.

Remark 2 (on tensor computation). Equations (27) give directly the stiffness
tensors in reduced form, and this is in general highly desirable as their
integration over the element volume takes multiple evaluations (e.g. through
Gauss quadrature). Since the computational complexity highly depends on the
number of dofs of the tensor, it is preferable to integrate directly the
reduced ones as long as the number of reduced coordinates $m$ is lower than
the number of element’s dofs $n_{e}$ (e.g. $n_{e}=60$ for a serendipity
hexahedron with quadratic shape functions). In case $m>n_{e}$, it is
computationally more efficient to compute the element tensors first (using Eqs
(27) and replacing both $\Gamma$ and $\Upsilon$ with G) for Gauss integration,
then project the element tensors using V and U accordingly. A similar
reasoning can be done for $m_{d}$, but under the very likely hypothesis that
$m_{d}\ll n_{e}$ it results almost always convenient to adopt the reduced
form.

###  5.3 Volume integration

The tensors in Eqs. (28) can be computed once and for all integrating on the
nominal volume $V_{o}$, under the said hypothesis that the defect does not
change the defected volume $V_{d}$ significantly. When this hypothesis cannot
be made, one can adopt the following approximation. Let $\textbf{Q}_{e}$ be
the generic expression of a tensor to be integrated over the volume $V_{d}$
(element level). We can compute the final reduced tensor Q as:

$\textbf{Q}=\sum_{e=1}^{N_{el}}\int_{V_{d}}\textbf{Q}_{e}\text{
d}V_{d}=\sum_{e=1}^{N_{el}}\int_{V_{o}}\textbf{Q}_{e}\text{det}(\textbf{F}_{1})\text{
d}V_{o}.$ (31)

where $N_{el}$ is the total number of elements. The determinant of
$\textbf{F}_{1}$ can now be approximated retaining only first order terms. To
the purpose of illustration, let us consider the following 2D example where
the global defect is given by the linear superposition of two shape defects,
that is:

$\textbf{u}_{d}(\textbf{x}_{0},\boldsymbol{\xi})=\begin{Bmatrix}u_{d}\\\
v_{d}\end{Bmatrix}=\begin{bmatrix}f_{u}^{(1)}(\textbf{x}_{0})&f_{u}^{(2)}(\textbf{x}_{0})\\\
f_{v}^{(1)}(\textbf{x}_{0})&f_{v}^{(2)}(\textbf{x}_{0})\end{bmatrix}\begin{Bmatrix}\xi_{1}\\\
\xi_{2}\end{Bmatrix}=\left[\textbf{f}^{\,(1)},\,\textbf{f}^{\,(2)}\right]\boldsymbol{\xi}.$
(32)

where we denote with $\textbf{f}^{\,(i)}=[f_{u}^{(i)},\,f_{v}^{(i)}]^{T}$ the
vector of the functions describing the i-th shape-defect for the
x-displacement $u_{d}$ and the y-displacement $v_{d}$, respectively. We can
approximate the determinant of $\textbf{F}_{1}$ as

$\displaystyle\text{det}(\textbf{F}_{1})=$
$\displaystyle\,1+u_{d,x}+v_{d,y}+u_{d,x}v_{d,y}-u_{d,y}v_{d,x}$
$\displaystyle=$
$\displaystyle\,1+\xi_{1}\left(f_{u,x}^{(1)}+f_{v,y}^{(1)}\right)+\xi_{2}\left(f_{u,x}^{(2)}+f_{v,y}^{(2)}\right)+\xi_{1}^{2}\left(f_{u,x}^{(1)}f_{v,y}^{(1)}-f_{u,y}^{(1)}f_{v,x}^{(1)}\right)+$
$\displaystyle+\xi_{1}\xi_{2}\left(f_{u,x}^{(1)}f_{v,y}^{(2)}+f_{u,x}^{(2)}f_{v,y}^{(1)}-f_{u,y}^{(1)}f_{v,x}^{(2)}-f_{u,y}^{(2)}f_{v,x}^{(1)}\right)+\xi_{2}^{2}\left(f_{u,x}^{(2)}f_{v,y}^{(2)}-f_{u,y}^{(2)}f_{v,x}^{(2)}\right)$
$\displaystyle\approx$
$\displaystyle\,1+\xi_{1}\left(f_{u,x}^{(1)}+f_{v,y}^{(1)}\right)+\xi_{2}\left(f_{u,x}^{(2)}+f_{v,y}^{(2)}\right)$

where we neglected higher order terms. Generalizing this result for $m_{d}$
defects, we can write

$\text{det}(\textbf{F}_{1})\approx 1+\sum_{i=1}^{m_{d}}\xi_{i}\left(\text{div
}\textbf{f}^{\,(i)}\right)$ (33)

so that Eq. (31) can be approximated as:

$\textbf{Q}\approx\sum_{e=1}^{N_{el}}\int_{V_{o}}\textbf{Q}_{e}\text{
d}V_{0}+\sum_{i=1}^{m_{d}}\left(\xi_{i}\sum_{e=1}^{N_{el}}\int_{V_{o}}\textbf{Q}_{e}\left(\text{div
}\textbf{f}^{\,(i)}\right)\text{
d}V_{o}\right)=\textbf{Q}^{\prime}+\sum_{i=1}^{m_{d}}\xi_{i}\textbf{Q}^{\prime\prime}_{i}$
(34)

where $\textbf{Q}^{\prime}$ is the tensor evaluated on the nominal volume and
$\textbf{Q}^{\prime\prime}_{i}$ is the contribution of the i-th defect, which
can be computed again once and for all offline, referring to the nominal
volume. The additional computational burden to compute
$\textbf{Q}^{\prime\prime}_{i}$ grows less than linearly with the number of
defects, since in a quadrature integration scheme we can use the
$\textbf{Q}_{e}$ evaluated at integration points (using Eqs. (27)) both for
$\textbf{Q}^{\prime}_{i}$ and for $\textbf{Q}^{\prime\prime}_{i}$. The
additional computations therefore involve only scalar by tensor
multiplications and tensor sums, so that most of the added computational time
is merely due to memory access management. Notice that one could also compute
all the additional tensors needed to describe det($\textbf{F}_{1}$) with no
approximation (even tough this is in most cases unnecessary, for h.o.t. do not
improve accuracy significantly). However, the first order approximation we
presented has the advantage to introduce only one new term for every
additional defect. Finally, if adopting this correction one should use the
mass matrix computed over the defected volume $V_{d}$ in accordance with the
present formulation. This way a new mass matrix $\textbf{M}_{d}$ must be
computed for each new parameter realization, however, being the matrix
constant during the analysis, this additional cost is negligible.

Remark 3 (on computational efficiency). The corrective terms
$\textbf{Q}^{\prime\prime}_{i}$ in Eq. (34) are null for an isochoric
transformation between nominal and defected domain (det($\textbf{F}_{1}$)=1).
In practice, one can set up a procedure to avoid the computation of these
terms to speed up the construction of the reduced tensors.

##  6 Reduction Basis

To construct the system described so far it is necessary to select the bases V
and U. The latter is simply a collection of user-defined displacement vectors,
each representing one specific defect, so that the final defected shape is
given by a linear superposition (see Eq. (26b)). The (properly said) RB is V,
whose choice may not be trivial, as it must correctly represent the system
response over a range of parameters without the possibility to be changed
(since a change would require to recompute the stiffness tensors). As
previously done in [1], our choice is to use a modal-based approach including
VMs, MDs and Defect-Sensitivities (DSs) [56] in the RB, as this solution
offers a way to construct a basis in a direct way, that is without convoluted
basis selection strategies, the need of computing all (or an excessively high
number of) eigenvectors or the need for previously computed simulations. We
remark however that, in principle, one could use also other RBs, as long as
they are valid over the parameter space.

Let us consider the following eigenvalue problem

$\left(\textbf{K}_{t}-\omega_{i}^{2}\textbf{M}\right)\boldsymbol{\Phi}_{i}=\textbf{0}$
(35)

where $\textbf{K}_{t}=\textbf{K}_{t}(\textbf{u},\textbf{u}_{d})$ is the
tangent stiffness matrix, M the mass matrix, $\omega_{i}$ the i-th
eigenfrequency and $\boldsymbol{\Phi}_{i}$ the corresponding eigenvector.
Static Modal Derivatives $\boldsymbol{\theta}_{ij}$ (MDs) are computed
neglecting the mass term, by taking the derivative of Eq. (35) with respect to
$\eta_{j}$ and evaluating the resulting expression at equilibrium (i.e.
$\eta_{j}=0$) and for $\boldsymbol{\xi}=\textbf{0}$:

$\boldsymbol{\theta}_{ij}=\frac{\partial\boldsymbol{\Phi}_{i}}{\partial\eta_{j}}\bigg{|}_{0}=-\textbf{K}_{0}^{-1}\frac{\partial\textbf{K}_{t}(\boldsymbol{\Phi}_{j}\eta_{j},\textbf{0})}{\partial\eta_{j}}\bigg{|}_{0}\boldsymbol{\Phi}_{i}.$
(36)

Defect Sensitivities (DSs) $\boldsymbol{\Xi}_{i,j}$ can be obtained following
a similar procedure, differentiating with respect to $\xi_{j}$:

$\boldsymbol{\Xi}_{i,j}=\frac{\partial\boldsymbol{\Phi}_{i}}{\partial\xi_{j}}\bigg{|}_{0}=-\textbf{K}_{0}^{-1}\frac{\partial\textbf{K}_{t}(\textbf{0},\textbf{U}_{j}\xi_{j})}{\partial\xi_{j}}\bigg{|}_{0}\boldsymbol{\Phi}_{i}.$
(37)

Expressions for the tangent stiffness derivatives are given in B.

Table 1: Acronyms for the different models considered in the numerical
studies.

Model | Description
---|---
FOM-d | Full Order Model with defect included by shifting the mesh nodes from the nominal configuration
| (no approximation). It is the reference model.
ROM-d | Reduced Order Model computed from FOM-d. Its reduction basis comprises VMs and MDs.
DpROM-N0 | Defect parametric Reduced Order Model, based on the 0-th order Neumann expansion (see Eq. (14)).
| Its reduction basis comprises VMs, MDs and DSs. Tensors are up to the 4-th
order (see [1]).
DpROM-N1 | Defect parametric Reduced Order Model, based on the 1-st order Neumann expansion (see Eq. (13)).
| Its reduction basis comprises VMs, MDs and DSs. Tensors are up to the 6-th
order (see Eqs. (27)).
DpROM-N1t | Truncated version of DpROM-N1, where terms of order $\mathcal{O}(\textbf{D}_{d}\textbf{D}^{2})$ in Eq. (13) have been neglected.
| As a consequence, $\tensor[_{5d}]{\textbf{Q}}{}$,
$\tensor[_{5dd}]{\textbf{Q}}{}$, $\tensor[_{6dd}]{\textbf{Q}}{}$ and the last
two terms in Eq. (27e) are null.
-v (suffix) | as the corresponding DpROM, but with the volume integration correction described in Section 5.3.

Remark 4 (on higher derivatives). Given the higher accuracy of the model,
larger defect magnitudes can be considered as compared to [1]. To fully
exploit the increased applicability range, a richer RB might be necessary,
reason why we here introduce second Defect-Sensitivities (DS2s) and MDs-
Sensitivities (MDSs). Let us take the derivative of Eq. (36) with respect to
the k-th defect amplitude $\xi_{k}$. We define the MDS
$\boldsymbol{\theta}_{ij,k}$ as:

$\boldsymbol{\theta}_{ij,k}=\frac{\partial\boldsymbol{\theta}_{ij}}{\partial\xi_{k}}\bigg{|}_{0}=-\textbf{K}_{0}^{-1}\left(\frac{\partial\textbf{K}_{t}}{\partial\xi_{k}}\bigg{|}_{0}\boldsymbol{\theta}_{ij}+\frac{\partial^{2}\textbf{K}_{t}}{\partial\eta_{j}\partial\xi_{k}}\bigg{|}_{0}\boldsymbol{\Phi}_{i}+\frac{\partial\textbf{K}_{t}}{\partial\eta_{j}}\bigg{|}_{0}\boldsymbol{\Xi}_{i,k}\right).$
(38)

Notice that $\boldsymbol{\theta}_{ij,k}\neq\boldsymbol{\theta}_{ji,k}$. In the
same manner, the second Defect Sensitivities (DS2s) with respect to $\xi_{k}$
write:

$\boldsymbol{\Xi}_{i,jk}=\frac{\partial\boldsymbol{\Xi}_{i,j}}{\partial\xi_{k}}\bigg{|}_{0}=-\textbf{K}_{0}^{-1}\left(\frac{\partial\textbf{K}_{t}}{\partial\xi_{k}}\bigg{|}_{0}\boldsymbol{\Xi}_{i,j}+\frac{\partial^{2}\textbf{K}_{t}}{\partial\xi_{j}\partial\xi_{k}}\bigg{|}_{0}\boldsymbol{\Phi}_{i}+\frac{\partial\textbf{K}_{t}}{\partial\xi_{j}}\bigg{|}_{0}\boldsymbol{\Xi}_{i,k}\right).$
(39)

It is evident that the blind inclusion of DS2s and/or MDSs in the RB would add
an unacceptable number of unknowns, especially when considering MDSs.
Depending on the type of the analysis (linear/nonlinear), on the kind of the
defect (i.e. affecting the linear or the nonlinear dynamics) and on the entity
of the defect itself (large/small), one can decide whether to include some
vectors or not. Pre-selection strategies to reduce the basis size, as the one
presented in [57] and [32], are beyond the scopes of this work and are not
treated hereafter.

##  7 Numerical tests – I

Figure 3: (a) Model-I: Nominal mesh and defected mesh with $\xi=1$ (and a x5
scale factor). (b) Frequency Responses and backbone curves for different
defect amplitudes $\xi$ using the Harmonic Balance method (7 harmonics) for:
ROM-d ( ), DpROM-N1 ( ), DpROM-N0 ( , only backbones). Backbones have been
computed also for the FOM-d ( ) using the shooting method for validation. The
vertical displacement of the mid-span of the beam is shown (first harmonic
coefficient of the Fourier series, normalized over the beam thickness). For
each plot, the detuning parameter $\sigma$ is referred to the corresponding
FOM-d first eigenfrequency. The bottom-right figure collects the backbone
curves for comparison.

We consider now a FE model of an aluminum beam, of length $l_{x}=2\,m$,
thickness $t_{y}=50\,mm$ and width $w_{z}=0.2\,m$, clamped at both ends. We
use a 2D plain strain model, with a mesh of 80 quadrilateral elements with
quadratic shape functions (630 dofs). A Rayleigh damping matrix
$\textbf{C}_{d}=\alpha\textbf{M}_{d}+\beta\textbf{K}_{0}$ is introduced by
imposing a quality factor $Q_{1}=Q_{2}=100$ on the first and second modes of
the linear system ($\alpha=3.1,\,\beta=6.3\cdot 10^{-6}$). A nodal load $F$ is
applied to the mid-span of the beam (with $F=1$ kN and $F=4$ kN for the forced
responses). A shape defect defined as the vertical translation of the nodes
$v_{d}(x,\xi)=\xi t_{y}\sin(\pi/l_{x}x)$ is imposed, deforming the nominal
geometry of the straight beam into a shallow arch. Notice that this kind of
defect represents an isochoric transformation (see Section 5.3).

Table 1 reports the acronyms used for the models of this and the next
numerical studies. For $\xi=\\{0,\,0.25,\,0.5,\,0.75,\,1\\}$, backbones and
Frequency Responses (FR) are computed for ROM-d, DpROM-N1 and DpROM-N0,
constructed using 5 VMs, 15 MDs, and 5 DSs (only for DpROMs). The Harmonic
Balance (HB) method was used (with 7 harmonics) using the NLvib Matlab tool
[58] (slightly modified to adapt the direct use of tensors) and our in-house
Matlab FE code. To validate the results of the ROMs, the Shooting Method is
used to find the backbones of the corresponding FOM-d. Results are shown in
Fig. 3. Computations were carried out in Matlab 2020a on a local machine
equipped with an Intel(R) Xeon(R) Silver 4214 CPU @2.20 GHz and 256 GB RAM
@2666 MHz. Tensors were built in a Julia sub-routine, called by the main
Matlab code, which uses the TensorOperations package [59] for the tensor
contraction. At present, the tensor construction is implemented serially,
therefore leaving space for possible future speed-ups exploiting parallel
computing, as remarked earlier.

Table 2: Average computational times. The data of the two DpROMs, being very similar, are clustered togeter. For the FR, the times refer to the higher forcing ($F=4$ kN). ROM-d, DpROM-N1/0 and FOM-d have 20, 25 (due to the 5 additional DSs) and 630 dofs, respectively. Notice that, being the size of the FOM-d very small, no significant conclusions in terms of speedups can be drawn from this data (refer to the next section for a detailed discussion). Model–I |  | ROM-d | DpROM-N1/0 | FOM-d
---|---|---|---|---
Harmonic Balance (HB) | Frequency Response | 649 s | 673 s | –
| Backbone | 237 s | 273 s | –
Shooting | Backbone | 31 min | 35 min | 83 h 18 min
ROM construction |  | 0.97 s | 6.5 s / 2 s | –

As it can be observed, the shift from hardening to softening behavior is well
captured by all the models, with a minor loss of accuracy of the DpROMs as
$\xi$ increases. In particular, DpROM-N0 shows a significant frequency offset
of the first linear eigenfrequency which remains constant throughout the
backbone curve (the same happens for the FRs, but we omitted to plot them for
the sake of figures clarity). The main goal of the present test was to assess
the accuracy of the method verifying the results against the FOM and over a
range of frequencies. However, computational times are collected in Table 2
for completeness. Run times for the shooting method with the (Dp)ROMs are
included for comparison. These figures, however, must be taken just as an
indication, first because of the difference between FOM and ROMs in terms of
convergence during continuation (ROMs are less likely to incur into numerical
artifacts) and, secondly, because speed and convergence of this kind of
analysis is highly sensitive to several parameters and finding the best
combination by trial-and-error usually leads to sub-optimal performances. Last
but not least, the size of the FOM in this case is too low to really
appreciate the savings in terms of ROM construction.

##  8 Numerical tests – II

###  8.1 MEMS gyroscope

The last example we present is a prototype MEMS mono-axial gyroscope, shown in
Fig. 4(a). The device consists in a mass suspended by four S-shaped springs,
connected to the ground on the bottom of the anchors. It is a monolithic
piece, produced via Deep Reactive Ion Etching (DRIE), a process which removes
material from a plane silicon wafer to obtain the final geometry. The etching
procedure is the main cause of production shape defects of MEMS devices, as it
will be detailed later. In operative conditions, the mass is kept in motion by
comb finger electrodes at the natural frequency of the drive mode (i.e. a mode
featuring motion mainly in the x-direction), so that in presence of an
external angular rate $\Omega$ (along the y-axis) a vertical displacement
$w_{sense}$ arises due to Coriolis effect along the z-axis (sense). The latter
is then converted into an electrical signal through the parallel plate
electrode placed on the ground below the mass, providing the measure for the
angular rate. In general, a defect or a combination of them may create a
coupling between the x- and z-axes so that the drive motion generates an
additional out-of-plane displacement which superimposes to the Coriolis
displacement to be measured. This is usually referred to as quadrature error
since, being proportional to the drive displacement, it is in phase quadrature
with the Coriolis signal, proportional to the drive velocity. Though it is
possible to tell apart the two contributions, this is highly undesirable as it
requires dedicated, over-sized electronics to accommodate the larger
displacements. Ultimately, this results in higher power consumption.

(a)

(b)

Figure 4: (a) Model–II, a MEMS gyroscope. Drive and sense direction are
indicated by arrows. (b) Meshed model, with 14,920 quadratic hexahedra, for a
total of 87,767 nodes and 261,495 dofs. Approximate dimension are given.

###  8.2 FE model, defects and simulation details

The FE model is shown in Fig. 4(b) and describes in detail the geometry and
mesh of the device, counting 14,920 quadratic hexahedral elements for a total
of 261,495 dofs. For the present study we selected two typical defects
occurring in the production of MEMS devices, namely the wall angle (shown in
Fig. 5(a)) and a restriction of the cross section of the beams (Fig. 5(b)).
The first is generated by the fact that the plasma beam of the DRIE process
might be not perfectly perpendicular to the working plane, while the second
one typically comes from an overexposure to the chemical attacks (over-
etching). In the spirit of our method, we can describe the global defects as
the superposition of these two displacement fields (see Eq. (26b)), letting
$\textbf{U}=[\textbf{U}_{1},\,\textbf{U}_{2}]$ with the associated amplitude
parameter vector $\boldsymbol{\xi}=[\xi_{1},\,\xi_{2}]^{T}$. The wall angle
shape defect
$\textbf{U}_{1}=[\textbf{u}_{d1},\,\textbf{v}_{d1},\,\textbf{w}_{d1}]^{T}$ is
given by

$u_{d1}(\xi_{1},z)=\xi_{1}\tan(\alpha_{y})z,$ (40)

and $v_{d1}=w_{d1}=0$. The tapering of the beams
$\textbf{U}_{2}=[\textbf{u}_{d2},\,\textbf{v}_{d2},\,\textbf{w}_{d2}]^{T}$ is
defined as

$v_{d2}(\xi_{2},x,y)=\xi_{2}\sin\left(\frac{\pi}{L_{b}}(x-x_{off})\right)(y-y_{mid})\frac{2}{W_{b}},$
(41)

and $u_{d2}=w_{d2}=0$, where $L_{b}$ and $W_{b}$ are the length and the width
of the beam, $x_{off}$ an offset depending on the location of each beam and
$y_{mid}$ is the y-coordinate corresponding to the middle line of each beam.
To ease the interpretation of the amplitude parameters, in the following
$\xi_{1}$ is reported in degrees to represent the physical wall angle coming
from the product $\xi_{1}\tan(\alpha_{y})$ in Eq. (40), while $\xi_{2}$ is
reported as a percentage of the beam thickness.

(a)

(b)

Figure 5: (a) First defect shape, $\textbf{U}_{1}$: constant wall angle
$\alpha_{y}$ (only one beam is shown). The colormap indicates x-displacement.
(b) Second defect shape, $\textbf{U}_{2}$: tapering of the suspension beams,
3D and top views (only one beam is shown). The colormap indicates
y-displacement (absolute value).

We compute the FR of the MEMS gyroscope using the NLvib Matlab tool and our
in-house Matlab FE code. We used a reduction basis with 3 VMs, the
corresponding 6 MDs and 3 DSs per defect (only for the DpROMs), and we used
the HB method with $H=5$ harmonics (with $N_{s}=3H+1$ time samples per period,
the minimum number of samples by which no sampling error is introduced in the
harmonics up the the H-th order when considering polynomial nonlinearities up
to the third order [60], as in our case). Given the size of the model, we take
as reference the results of ROM-d, as it would be prohibitively time and
memory consuming to compute the frequency response for FOM-d. Apart from the
practical issues, we justify this choice considering on the one side the good
results obtained for lower dimensional models (as the one presented in the
previous section), and on the other side considering that, ultimately, our
DpROMs will be at best as good as ROM-d, which is not parametric and not
approximated in its formulation.

The frequency response was obtained forcing the system in the center of the
suspended mass with a nodal load directed along the x-direction, with
amplitude $F=0.4\,\mu$N, and using a Rayleigh damping matrix with $\alpha=105$
and $\beta=0$. Figure 6 reports the FRs around the first eigenfrequency of the
system for the x-displacement $u$ (drive direction) and the z-displacement $w$
(sense direction) for all the combinations of
$\xi_{1}=\\{0^{\circ},\,0.5^{\circ},\,1^{\circ}\\}$ and
$\xi_{2}=\\{0\%,\,0.5\%,\,1\%,\,1.5\%,\,2\%\\}$.

###  8.3 Tested models

For the present study, we considered also a truncated version of DpROM-N1
(named DpROM-N1t) where we considered negligible in Eq. (11) the strains of
order $\mathcal{O}(\textbf{D}_{d}\textbf{D}^{2})$. This further assumption
allows us to drop all the terms multiplying $\textbf{A}_{3}$ or, equivalently,
$\textbf{L}_{3}$, so that $\tensor[_{5d}]{\textbf{Q}}{}$,
$\tensor[_{5dd}]{\textbf{Q}}{}$, $\tensor[_{6dd}]{\textbf{Q}}{}$ and the last
two terms in Eq. (27e) are null. Although introducing a new approximation,
DpROM-N1t is significantly cheaper to construct and, as it will be shown, does
not introduce any appreciable accuracy loss in our studies. For each of the
presented pROMs then, we test the same models with the volume-integration
correction described in Section 5.3. We will address to these models appending
the suffix “-v” to the name of the model itself (e.g. DpROM-N1t-v). To recap,
the results for a total of 6 models are reported in the following:
DpROM-N0(-v), DpROM-N1t(-v) and DpROM-N1(-v). Again, see Table 1 for a quick
reference.

(a) Drive response.

(b) Sense response (out of plane).

Figure 6: Frequency Responses and backbone curves for different defect
amplitudes $\xi_{1}$ (shown in degrees) and $\xi_{2}$ (reported as percentage
of the beam width) using the Harmonic Balance method (5 harmonics) for: ROM-d
( ), DpROM-N0 ( ), DpROM-N1t ( ), DpROM-N1 ( ), DpROM-N0-v ( ), DpROM-N1t-v (
) and DpROM-N1-v ( ). The displacements $u$ and $w$ of the center of the mass
is shown (first harmonic coefficient of the Fourier series). For each plot,
the percent detuning parameter $\sigma$ is referred to the corresponding FOM-d
first eigenfrequency.

###  8.4 Results

Considering first the effect of the wall angle defect only, it is apparent how
DpROM-N0 performances quickly degrade as soon as the parameter $\xi_{1}$ is
increased. This can be seen both in the error on the linear eigenfrequency and
especially in the overestimated $w$-response, approximately one order of
magnitude higher than the reference. This may be due to the fact that the
S-shaped beams are specifically designed to minimize the cross-coupling
between the drive (x-) and sense (z-) axes created by the wall angle, so that
the $w$-response is so small (2 orders of magnitude lower than the
$u$-response) that it cannot be accurately captured by DpROM-N0. The same
observations can be made for DpROM-N0-v, as the wall angle defect by itself
represent an isochoric transformation. The responses of all the other tested
DpROMs show instead a perfect match with the reference when $\xi_{2}=0\%$.

If the tapering defect only is considered (i.e. with $\xi_{1}=0^{\circ}$), we
observe that DpROM-N0/N1t/N1 have similar responses, with an error on the
eigenfrequency that translates the whole response by an approximately constant
$\Delta f$. The models with the volume-integration correction were then
tested. If on the one hand DpROM-N0-v still presents relevant errors, on the
other hand DpROM-N1t-v and DpROM-N1v show very accurate results in the full
range of the tested $\xi_{2}$. Such an improvement was expected, as the volume
changed by this defect affects the suspension beams dimensions, to which the
eigenfrequencies of the system are very sensitive.

For the remaining cases, the trends observed for the parameters $\xi_{1}$ and
$\xi_{2}$ individually mix. Notice that looking at some results (e.g.
$u$-response for $\xi_{1}=1^{\circ}$, $\xi_{2}=0.5\%$), it may seem that
DpROM-N0 gives better results than DpROM-N1. This is however just a
coincidence, as for DpROM-N0 the first defect shifts the first eigenfrequency
to lower frequencies while the second defect to higher frequencies, so that
the two errors in this case cancel out. Indeed, when the volume correction is
used in DpROM-N0-v, only the first effect is observed, and the frequencies are
shifted to the left.

Figure 7: Transient response of the center of the suspended mass of the
gyroscope for: ROM-d ( ), DpROM-N0-v ( ), DpROM-N1t-v ( ) and DpROM-N1-v ( ).
The forcing is harmonic with period T0. Case with $\xi_{1}=1^{\circ}$,
$\xi_{2}=2\%$.

In Fig. 7 we also show the transient response of the forced node for ROM-d and
DpROMs-v (case with $\xi_{1}=1^{\circ}$, $\xi_{2}=0.5\%$). Each model is
forced at its own first resonance frequency $f_{0}$ (as it is usually the case
for MEMS gyroscopes) with a harmonic forcing, taking 100 samples per period
and for a time span equal to 10 times T${}_{0}=1/f_{0}$, with $F=50\mu$N. The
integration was carried out in Matlab with our in-house code, using a Newmark
integration scheme. Looking at the responses along the three axes, we observe
that the three DpROMs yield correct results but for DpROM-N0-v along the sense
z-direction ($w$ component). Also, considering the z-response, we can see that
DpROM-N1 is slightly better DpROM-N1t, fact that was not very visible in the
FRs.

###  8.5 Computational times

Table 3 reports the average time for the FR analyses and for the construction
of the different models. To compare in terms of time ROM-d and the DpROMs, it
is convenient to consider the variable costs (Tvar), i.e. the ones that have
to be sustained for each new parameter realization, and the overhead costs
(Toh), i.e. the ones sustained once and for all independently from the number
of realizations. In the case of ROM-d we have that
$\text{T}_{var}=\text{T}_{constr}+\text{T}_{sim}$, being $\text{T}_{constr}$
the time to construct the model (i.e. RB and tensors computation) and
$\text{T}_{sim}$ the time for one simulation, while T${}_{oh}=0$. For ROM-d
indeed, there are no common overhead costs, but a new model must be
constructed for each new realization of the parameters. In the case of DpROMs
instead, we have that $\text{T}_{var}^{p}=\text{T}_{sim}^{p}$ and
$\text{T}_{oh}^{p}=\text{T}_{constr}^{p}$ (we use the superscript “p” to
distinguish the parametric models from ROM-d). For the parametric models we
have in fact to pay upfront the cost of model construction, which is generally
more expensive than the one for ROM-d, but thereafter only the simulation cost
must be sustained for each new case. The first trivial conclusion is then that
there exist a number $\bar{N}$ of parameter realization above which DpROMs
become convenient, that is:

$N\geq\bar{N}=\bigg{\lceil}\frac{T_{oh}^{p}}{T_{var}-T_{var}^{p}}\bigg{\rceil}.$
(42)

For $\bar{N}$ to be positive and finite, it follows that

$T_{var}>T_{var}^{p}\,\,\,\longleftrightarrow\,\,\,T_{constr}+T_{sim}>T_{sim}^{p}.$
(43)

From Eq. (43) it can be seen how the convenience of the parametric model over
the non-parametric one depends on the relative weight between the simulation
and construction times of the latter and the simulation time of the former, as
it can be observed in Table 3 looking at the different speedups666Speedups are
computed considering the variable costs only, with respect to ROM-d. for the
FR and transient analyses.

That said, it is clearly difficult to draw general and definitive conclusions
on the benefits of the two solutions, ROM-d and DpROMs, time-wise. In the
experience of the authors, transient analysis offer the best gains, as
simulation speed is very high, grows almost linearly with the simulation time
span, and is less sensitive to the number of dofs than other kind of analysis,
as the ones requiring continuation methods. When continuation is required, one
could potentially find greater benefits in using a model with a low number of
dofs, so that ROM-d could actually become the best choice. We remark however
that for ROM-d we have to take into account also the construction cost as a
variable cost, and that for large FE models the sole computation of structural
eigenmodes can already take several minutes, making this cost very high.

Table 3: Average computational times for the FR with the HB method and for the construction of each ROM (comprising the time to compute VMs, MDs, DSs and the tensors). ROM-d counts 9 dofs, while DpROMs 15 dofs. Apart from the construction costs, items for the DpROMs are clustered and averaged. Notice that the construction time for ROM-d is sustained for each new parameter realization and contributes to the variable costs. DpROM names are reported by their respective suffix. Model–II | ROM-d | N0 | N1t | N1 | N0-v | N1t-v | N1-v
---|---|---|---|---|---|---|---
ROM construction | 209 s | 335 s | 333 s | 816 s | 353 s | 357 s | 1,063 s
HB FR | 21.7 s | 49.2 s
Transient analysis | 0.43 s | 0.51 s
Overhead cost | – | 335 s | 333 s | 816 s | 353 s | 357 s | 1,063 s
Variable cost (FR/transient) | 230.7 s / 209.4 s | 49.2 s / 0.51 s
Speedup (FR/transient) | – / – |  $4.7\times$ / $410.6\times$

##  9 Conclusions

We presented a ROM for geometric nonlinearities that can parametrically
describe a shape imperfection with respect to the nominal (blueprint) design,
named for brevity DpROM. The imperfection is given by the superposition of
user-defined defect shapes, whose amplitudes are parameters of the model and
can be changed without reconstructing the model itself. This result has been
made possible thanks to a polynomial representation of the internal forces
resulting from a two-step deformation process (which brings the nominal
geometry into the defected one and then into the deformed one) and from the
approximation of the strains obtained by a Neumann expansion. The latter
allowed to eliminate rational expressions under the hypothesis of small
defects, so that the elastic internal forces are written as simple polynomials
both with respect to the displacement field representing the defect and with
respect to the actual displacement field. Using a Galerkin projection and a
modal-based approach for selecting the RB, the reduced internal forces have
been recast in tensorial form, where the linear, quadratic and cubic stiffness
tensors are found to be functions of a parameter vector collecting the
amplitudes of the defects imposed on the structure. Within this framework we
tested different versions of the DpROM for different degrees of approximation.
In particular, we have shown that the model we had previously developed using
Budiansky’s approach corresponds to the 0th-order expansion of our model,
without volume integration correction (i.e. DpROM-N0). Finally, in the
numerical studies we showed that the higher order approximation DpROM-N1
effectively leads to more accurate results and that for volume-changing
defects a large improvement can be achieved by approximating the tensor
integral over the real volume of the defective geometry (DpROM-N1-v). The
truncated version DpROM-N1t(-v) was also presented, which has almost the same
accuracy as its complete counterpart, but without the need to construct
tensors with dimensionality higher than four. The computational costs were
then critically discussed, taking into account different types of analysis. In
particular, we showed that in transient studies we can usually expect very
high speedups from the parametric models. In the case of FR analysis, which we
used to assess the quality of the solutions over a range of frequencies as an
alternative to multiple time analyses, the gains will be more contained. In
this context, to reduce the dofs of both the parametric and non-parametric
ROMs and make FR analysis faster and thus closer to transient analysis in
terms of time and speedups, we think that an a priori selection of the RB
vectors and hyperreduction strategies would actually be very beneficial, and
they can actually constitute the spur for future investigation.

## Acknowledgements

The authors want to express their gratitude to the AMS division at
STMicroelectronics in Cornaredo (MI, Italy), for their support in this work.

##  Appendix A $\textbf{L}_{1}$, $\textbf{L}_{2}$ and $\textbf{L}_{3}$
matrices

We report in Tables 4 and 5 the expressions for the matrices $\textbf{L}_{1}$,
$\textbf{L}_{2}$ and $\textbf{L}_{3}$ defined in Eqs. 23.

Table 4: Elements $L^{(1)}_{ijk}$, $L^{(2)}_{ijk}$ and $L^{(3)}_{ijkl}$ of the
sparse $3\times 4\times 4$ matrices $\textbf{L}_{1}$, $\textbf{L}_{2}$ and of
the sparse $3\times 4\times 4\times 4$ matrix $\textbf{L}_{3}$, respectively,
in the 2D case (name subscripts are moved to superscripts to avoid confusion
with the indexes).

$L^{(1)}_{111}=1$, |  $L^{(1)}_{321}=1$, |  $L^{(1)}_{312}=1$, |  $L^{(1)}_{222}=1$, |  $L^{(1)}_{123}=1$, |  $L^{(1)}_{343}=1$, |  $L^{(1)}_{324}=1$, |  $L^{(1)}_{244}=1$.
---|---|---|---|---|---|---|---
$L^{(2)}_{111}=1$, |  $L^{(2)}_{331}=1$, |  $L^{(2)}_{312}=1$, |  $L^{(2)}_{232}=1$, |  $L^{(2)}_{123}=1$, |  $L^{(2)}_{343}=1$, |  $L^{(2)}_{324}=1$, |  $L^{(2)}_{244}=1$.
$L^{(3)}_{1111}=1$, |  $L^{(3)}_{3211}=\frac{1}{2}$, |  $L^{(3)}_{3121}=\frac{1}{2}$, |  $L^{(3)}_{1331}=1$, |  $L^{(3)}_{3431}=\frac{1}{2}$, |  $L^{(3)}_{3341}=\frac{1}{2}$, |  $L^{(3)}_{3112}=1$, |  $L^{(3)}_{2212}=\frac{1}{2}$,
$L^{(3)}_{2122}=\frac{1}{2}$, |  $L^{(3)}_{3332}=1$, |  $L^{(3)}_{2432}=\frac{1}{2}$, |  $L^{(3)}_{2342}=\frac{1}{2}$, |  $L^{(3)}_{1213}=\frac{1}{2}$, |  $L^{(3)}_{1123}=\frac{1}{2}$, |  $L^{(3)}_{3223}=1$, |  $L^{(3)}_{1433}=\frac{1}{2}$,
$L^{(3)}_{1343}=\frac{1}{2}$, |  $L^{(3)}_{3443}=1$, |  $L^{(3)}_{3214}=\frac{1}{2}$, |  $L^{(3)}_{3124}=\frac{1}{2}$, |  $L^{(3)}_{2224}=1$, |  $L^{(3)}_{3434}=\frac{1}{2}$, |  $L^{(3)}_{3344}=\frac{1}{2}$, |  $L^{(3)}_{2444}=1$.

Table 5: Elements $L^{(1)}_{ijk}$, $L^{(2)}_{ijk}$ and $L^{(3)}_{ijkl}$ of the
sparse $6\times 9\times 9$ matrices $\textbf{L}_{1}$, $\textbf{L}_{2}$ and of
the sparse $6\times 9\times 9\times 9$ matrix $\textbf{L}_{3}$, respectively,
in the 3D case (name subscripts are moved to superscripts to avoid confusion
with the indexes).

$L^{(1)}_{111}=1$, |  $L^{(1)}_{421}=1$, |  $L^{(1)}_{531}=1$, |  $L^{(1)}_{412}=1$, |  $L^{(1)}_{222}=1$, |  $L^{(1)}_{632}=1$, |  $L^{(1)}_{513}=1$, |  $L^{(1)}_{623}=1$, |  $L^{(1)}_{333}=1$,
---|---|---|---|---|---|---|---|---
$L^{(1)}_{144}=1$, |  $L^{(1)}_{454}=1$, |  $L^{(1)}_{564}=1$, |  $L^{(1)}_{445}=1$, |  $L^{(1)}_{255}=1$, |  $L^{(1)}_{665}=1$, |  $L^{(1)}_{546}=1$, |  $L^{(1)}_{656}=1$, |  $L^{(1)}_{366}=1$,
$L^{(1)}_{177}=1$, |  $L^{(1)}_{487}=1$, |  $L^{(1)}_{597}=1$, |  $L^{(1)}_{478}=1$, |  $L^{(1)}_{288}=1$, |  $L^{(1)}_{698}=1$, |  $L^{(1)}_{579}=1$, |  $L^{(1)}_{689}=1$, |  $L^{(1)}_{399}=1$.
$L^{(2)}_{111}=1$, |  $L^{(2)}_{441}=1$, |  $L^{(2)}_{571}=1$, |  $L^{(2)}_{412}=1$, |  $L^{(2)}_{242}=1$, |  $L^{(2)}_{672}=1$, |  $L^{(2)}_{513}=1$, |  $L^{(2)}_{643}=1$, |  $L^{(2)}_{373}=1$,
$L^{(2)}_{124}=1$, |  $L^{(2)}_{454}=1$, |  $L^{(2)}_{584}=1$, |  $L^{(2)}_{425}=1$, |  $L^{(2)}_{255}=1$, |  $L^{(2)}_{685}=1$, |  $L^{(2)}_{526}=1$, |  $L^{(2)}_{656}=1$, |  $L^{(2)}_{386}=1$,
$L^{(2)}_{137}=1$, |  $L^{(2)}_{467}=1$, |  $L^{(2)}_{597}=1$, |  $L^{(2)}_{438}=1$, |  $L^{(2)}_{268}=1$, |  $L^{(2)}_{698}=1$, |  $L^{(2)}_{539}=1$, |  $L^{(2)}_{669}=1$, |  $L^{(2)}_{399}=1$.
$L^{(3)}_{1111}=1$, |  $L^{(3)}_{4211}=\frac{1}{2}$, |  $L^{(3)}_{5311}=\frac{1}{2}$, |  $L^{(3)}_{4121}=\frac{1}{2}$, |  $L^{(3)}_{5131}=\frac{1}{2}$, |  $L^{(3)}_{1441}=1$, |  $L^{(3)}_{4541}=\frac{1}{2}$, |  $L^{(3)}_{5641}=\frac{1}{2}$, |  $L^{(3)}_{4451}=\frac{1}{2}$,
$L^{(3)}_{5461}=\frac{1}{2}$, |  $L^{(3)}_{1771}=1$, |  $L^{(3)}_{4871}=\frac{1}{2}$, |  $L^{(3)}_{5971}=\frac{1}{2}$, |  $L^{(3)}_{4781}=\frac{1}{2}$, |  $L^{(3)}_{5791}=\frac{1}{2}$, |  $L^{(3)}_{4112}=1$, |  $L^{(3)}_{2212}=\frac{1}{2}$, |  $L^{(3)}_{6312}=\frac{1}{2}$,
$L^{(3)}_{2122}=\frac{1}{2}$, |  $L^{(3)}_{6132}=\frac{1}{2}$, |  $L^{(3)}_{4442}=1$, |  $L^{(3)}_{2542}=\frac{1}{2}$, |  $L^{(3)}_{6642}=\frac{1}{2}$, |  $L^{(3)}_{2452}=\frac{1}{2}$, |  $L^{(3)}_{6462}=\frac{1}{2}$, |  $L^{(3)}_{4772}=1$, |  $L^{(3)}_{2872}=\frac{1}{2}$,
$L^{(3)}_{6972}=\frac{1}{2}$, |  $L^{(3)}_{2782}=\frac{1}{2}$, |  $L^{(3)}_{6792}=\frac{1}{2}$, |  $L^{(3)}_{5113}=1$, |  $L^{(3)}_{6213}=\frac{1}{2}$, |  $L^{(3)}_{3313}=\frac{1}{2}$, |  $L^{(3)}_{6123}=\frac{1}{2}$, |  $L^{(3)}_{3133}=\frac{1}{2}$, |  $L^{(3)}_{5443}=1$,
$L^{(3)}_{6543}=\frac{1}{2}$, |  $L^{(3)}_{3643}=\frac{1}{2}$, |  $L^{(3)}_{6453}=\frac{1}{2}$, |  $L^{(3)}_{3463}=\frac{1}{2}$, |  $L^{(3)}_{5773}=1$, |  $L^{(3)}_{6873}=\frac{1}{2}$, |  $L^{(3)}_{3973}=\frac{1}{2}$, |  $L^{(3)}_{6783}=\frac{1}{2}$, |  $L^{(3)}_{3793}=\frac{1}{2}$,
$L^{(3)}_{1214}=\frac{1}{2}$, |  $L^{(3)}_{1124}=\frac{1}{2}$, |  $L^{(3)}_{4224}=1$, |  $L^{(3)}_{5324}=\frac{1}{2}$, |  $L^{(3)}_{5234}=\frac{1}{2}$, |  $L^{(3)}_{1544}=\frac{1}{2}$, |  $L^{(3)}_{1454}=\frac{1}{2}$, |  $L^{(3)}_{4554}=1$, |  $L^{(3)}_{5654}=\frac{1}{2}$,
$L^{(3)}_{5564}=\frac{1}{2}$, |  $L^{(3)}_{1874}=\frac{1}{2}$, |  $L^{(3)}_{1784}=\frac{1}{2}$, |  $L^{(3)}_{4884}=1$, |  $L^{(3)}_{5984}=\frac{1}{2}$, |  $L^{(3)}_{5894}=\frac{1}{2}$, |  $L^{(3)}_{4215}=\frac{1}{2}$, |  $L^{(3)}_{4125}=\frac{1}{2}$, |  $L^{(3)}_{2225}=1$,
$L^{(3)}_{6325}=\frac{1}{2}$, |  $L^{(3)}_{6235}=\frac{1}{2}$, |  $L^{(3)}_{4545}=\frac{1}{2}$, |  $L^{(3)}_{4455}=\frac{1}{2}$, |  $L^{(3)}_{2555}=1$, |  $L^{(3)}_{6655}=\frac{1}{2}$, |  $L^{(3)}_{6565}=\frac{1}{2}$, |  $L^{(3)}_{4875}=\frac{1}{2}$, |  $L^{(3)}_{4785}=\frac{1}{2}$,
$L^{(3)}_{2885}=1$, |  $L^{(3)}_{6985}=\frac{1}{2}$, |  $L^{(3)}_{6895}=\frac{1}{2}$, |  $L^{(3)}_{5216}=\frac{1}{2}$, |  $L^{(3)}_{5126}=\frac{1}{2}$, |  $L^{(3)}_{6226}=1$, |  $L^{(3)}_{3326}=\frac{1}{2}$, |  $L^{(3)}_{3236}=\frac{1}{2}$, |  $L^{(3)}_{5546}=\frac{1}{2}$,
$L^{(3)}_{5456}=\frac{1}{2}$, |  $L^{(3)}_{6556}=1$, |  $L^{(3)}_{3656}=\frac{1}{2}$, |  $L^{(3)}_{3566}=\frac{1}{2}$, |  $L^{(3)}_{5876}=\frac{1}{2}$, |  $L^{(3)}_{5786}=\frac{1}{2}$, |  $L^{(3)}_{6886}=1$, |  $L^{(3)}_{3986}=\frac{1}{2}$, |  $L^{(3)}_{3896}=\frac{1}{2}$,
$L^{(3)}_{1317}=\frac{1}{2}$, |  $L^{(3)}_{4327}=\frac{1}{2}$, |  $L^{(3)}_{1137}=\frac{1}{2}$, |  $L^{(3)}_{4237}=\frac{1}{2}$, |  $L^{(3)}_{5337}=1$, |  $L^{(3)}_{1647}=\frac{1}{2}$, |  $L^{(3)}_{4657}=\frac{1}{2}$, |  $L^{(3)}_{1467}=\frac{1}{2}$, |  $L^{(3)}_{4567}=\frac{1}{2}$,
$L^{(3)}_{5667}=1$, |  $L^{(3)}_{1977}=\frac{1}{2}$, |  $L^{(3)}_{4987}=\frac{1}{2}$, |  $L^{(3)}_{1797}=\frac{1}{2}$, |  $L^{(3)}_{4897}=\frac{1}{2}$, |  $L^{(3)}_{5997}=1$, |  $L^{(3)}_{4318}=\frac{1}{2}$, |  $L^{(3)}_{2328}=\frac{1}{2}$, |  $L^{(3)}_{4138}=\frac{1}{2}$,
$L^{(3)}_{2238}=\frac{1}{2}$, |  $L^{(3)}_{6338}=1$, |  $L^{(3)}_{4648}=\frac{1}{2}$, |  $L^{(3)}_{2658}=\frac{1}{2}$, |  $L^{(3)}_{4468}=\frac{1}{2}$, |  $L^{(3)}_{2568}=\frac{1}{2}$, |  $L^{(3)}_{6668}=1$, |  $L^{(3)}_{4978}=\frac{1}{2}$, |  $L^{(3)}_{2988}=\frac{1}{2}$,
$L^{(3)}_{4798}=\frac{1}{2}$, |  $L^{(3)}_{2898}=\frac{1}{2}$, |  $L^{(3)}_{6998}=1$, |  $L^{(3)}_{5319}=\frac{1}{2}$, |  $L^{(3)}_{6329}=\frac{1}{2}$, |  $L^{(3)}_{5139}=\frac{1}{2}$, |  $L^{(3)}_{6239}=\frac{1}{2}$, |  $L^{(3)}_{3339}=1$, |  $L^{(3)}_{5649}=\frac{1}{2}$,
$L^{(3)}_{6659}=\frac{1}{2}$, |  $L^{(3)}_{5469}=\frac{1}{2}$, |  $L^{(3)}_{6569}=\frac{1}{2}$, |  $L^{(3)}_{3669}=1$, |  $L^{(3)}_{5979}=\frac{1}{2}$, |  $L^{(3)}_{6989}=\frac{1}{2}$, |  $L^{(3)}_{5799}=\frac{1}{2}$, |  $L^{(3)}_{6899}=\frac{1}{2}$, |  $L^{(3)}_{3999}=1$.

##  Appendix B Tangent stiffness matrix derivatives

The virtual variation wrt u of the internal elastic forces as defined in Eq.
(22) writes

$\begin{split}\delta\textbf{f}_{int}&=\int_{V_{o}}\left[\textbf{G}^{T}\left(\textbf{H}+\textbf{A}_{1}+\textbf{A}_{2}+2\textbf{A}_{3}\textbf{A}_{1}\right)^{T}\textbf{C}\left(\textbf{H}+\frac{1}{2}\textbf{A}_{1}+\textbf{A}_{2}+\textbf{A}_{3}\textbf{A}_{1}\right)\textbf{G}\delta\textbf{u}+\right.\\\
&+\textbf{G}^{T}\left(\textbf{H}+\textbf{A}_{1}+\textbf{A}_{2}+2\textbf{A}_{3}\textbf{A}_{1}\right)^{T}\textbf{C}\left(\frac{1}{2}\delta\textbf{A}_{1}+\textbf{A}_{3}\delta\textbf{A}_{1}\right)\textbf{G}\textbf{u}+\\\
&+\left.\textbf{G}^{T}\left(\delta\textbf{A}_{1}+2\textbf{A}_{3}\delta\textbf{A}_{1}\right)^{T}\textbf{C}\left(\textbf{H}+\frac{1}{2}\textbf{A}_{1}+\textbf{A}_{2}+\textbf{A}_{3}\textbf{A}_{1}\right)\textbf{G}\textbf{u}\right]\text{
d}V_{o}.\end{split}$ (44)

Recalling that
$\textbf{A}_{1}\delta\boldsymbol{\theta}=\delta\textbf{A}_{1}\boldsymbol{\theta}$,
we can write

$\begin{split}\delta\textbf{f}_{int}&=\int_{V_{o}}\left[\textbf{G}^{T}\left(\textbf{H}+\textbf{A}_{1}+\textbf{A}_{2}+2\textbf{A}_{3}\textbf{A}_{1}\right)^{T}\textbf{C}\left(\textbf{H}+\textbf{A}_{1}+\textbf{A}_{2}+2\textbf{A}_{3}\textbf{A}_{1}\right)\textbf{G}\delta\textbf{u}+\right.\\\
&+\textbf{G}^{T}\delta\textbf{A}_{1}^{T}\underbrace{\left(\textbf{I}+2\textbf{A}_{3}^{T}\right)\textbf{C}\left(\textbf{H}+\frac{1}{2}\textbf{A}_{1}+\textbf{A}_{2}+\textbf{A}_{3}\textbf{A}_{1}\right)\textbf{G}\textbf{u}}_{\textbf{N}(\textbf{u},\textbf{u}_{d})}\Big{]}\text{
d}V_{o},\end{split}$ (45)

where the second term on the right-hand side can be rewritten to put in
evidence the displacement virtual variation $\delta\textbf{u}$ as

$\delta f_{I}^{\prime\prime}=\int_{V_{o}}G_{iI}L_{1jik}G_{kl}\delta
u_{l}N_{j}\text{ d}V_{o},$ (46)

where Einstein notation was used for convenience. The tangent stiffness matrix
therefore writes:

$\textbf{K}_{t}=\textbf{K}^{\prime}+\textbf{K}^{\prime\prime}$ (47)

where

$\textbf{K}^{\prime}=\int_{V_{o}}\textbf{G}^{T}\left(\textbf{H}+\textbf{A}_{1}+\textbf{A}_{2}+2\textbf{A}_{3}\textbf{A}_{1}\right)^{T}\textbf{C}\left(\textbf{H}+\textbf{A}_{1}+\textbf{A}_{2}+2\textbf{A}_{3}\textbf{A}_{1}\right)\textbf{G}\text{
d}V_{o},$ (48)
$K_{IJ}^{\prime\prime}=\int_{V_{o}}G_{iI}L_{1jik}G_{kJ}N_{j}\text{ d}V_{o}.$
(49)

Substituting $\textbf{u}=\boldsymbol{\Phi}_{i}\eta_{i}$ and
$\textbf{u}_{d}=\textbf{U}_{j}\xi_{j}$ in Eq. (47), taking the derivative wrt
either $\eta_{i}$ and/or $\xi_{j}$ and evaluating the resulting expressions at
equilibrium and with zero defect amplitudes, as required by Eqs. (36)–(39), we
can write the derivatives of $\textbf{K}_{t}$ as:

$\displaystyle\frac{\partial\textbf{K}_{t}}{\partial\eta_{j}}\bigg{|}_{0}$
$\displaystyle=\int_{V_{o}}\left[\textbf{G}^{T}\left(\textbf{H}^{T}\textbf{C}{\textbf{A}_{1}}_{j}+{\textbf{A}_{1}}_{j}^{T}\textbf{C}\textbf{H}\right)\textbf{G}+\textbf{G}\cdot_{11}\left[\left({\textbf{L}_{1}}\cdot\textbf{G}\right)\cdot_{11}\left(\textbf{C}\textbf{H}\textbf{G}\boldsymbol{\Phi}_{i}\right)\right]\right]\text{
d}V_{o},$ (50a)
$\displaystyle\frac{\partial\textbf{K}_{t}}{\partial\xi_{j}}\bigg{|}_{0}$
$\displaystyle=\int_{V_{o}}\textbf{G}^{T}\left(\textbf{H}^{T}\textbf{C}{\textbf{A}_{2}}_{j}+{\textbf{A}_{2}}_{j}^{T}\textbf{C}\textbf{H}\right)\textbf{G}\text{
d}V_{o},$ (50b)
$\displaystyle\begin{split}\frac{\partial^{2}\textbf{K}_{t}}{\partial\eta_{j}\partial\xi_{k}}\bigg{|}_{0}&=\int_{V_{o}}\left[\textbf{G}^{T}\left({\textbf{A}_{2}}_{k}^{T}\textbf{C}{\textbf{A}_{1}}_{j}+{\textbf{A}_{1}}_{j}^{T}\textbf{C}{\textbf{A}_{2}}_{k}+2{\textbf{A}_{1}}_{j}^{T}{\textbf{A}_{3}}_{k}^{T}\textbf{C}\textbf{H}+2\textbf{H}^{T}\textbf{C}{\textbf{A}_{3}}_{k}{\textbf{A}_{1}}_{j}\right)\textbf{G}+\right.\\\
&+\textbf{G}\cdot_{11}\left[\left({\textbf{L}_{1}}\cdot\textbf{G}\right)\cdot_{11}\left(\left(\textbf{C}\textbf{A}_{2}+2\textbf{A}_{3}^{T}\textbf{C}\textbf{H}\right)\textbf{G}\boldsymbol{\Phi}_{i}\right)\right]\Big{]}\text{
d}V_{o},\end{split}$ (50c)
$\displaystyle\frac{\partial^{2}\textbf{K}_{t}}{\partial\xi_{j}\partial\xi_{k}}\bigg{|}_{0}$
$\displaystyle=\int_{V_{o}}\textbf{G}^{T}\left({\textbf{A}_{2}}_{j}^{T}\textbf{C}{\textbf{A}_{2}}_{k}+{\textbf{A}_{2}}_{k}^{T}\textbf{C}{\textbf{A}_{2}}_{j}\right)\textbf{G}\text{
d}V_{o},$ (50d)

where, recalling that $\eta_{i}$ and $\xi_{j}$ are scalars, we used
$\textbf{A}_{1}(\textbf{G}\boldsymbol{\Phi}_{i}\eta_{i})={\textbf{A}_{1}}_{i}\eta_{i}$
and
$\textbf{A}_{2}(\textbf{G}\textbf{U}_{j}\xi_{j})={\textbf{A}_{2}}_{j}\xi_{j}$
(same for $\textbf{A}_{3}$) to avoid a cumbersome notation, and where
$\cdot_{ij}$ denotes the contraction of the i-th dimension of the first term
with the j-th dimension of the second term.

## References

  * Marconi et al. [2020] J. Marconi, P. Tiso, F. Braghin,  A nonlinear reduced order model with parametrized shape defects,  Computer Methods in Applied Mechanics and Engineering 360 (2020) 112785\.
  * Belytschko et al. [2013] T. Belytschko, W. K. Liu, B. Moran, K. Elkhodary, Nonlinear finite elements for continua and structures, John wiley & sons, 2013.
  * Toselli and Widlund [2005] A. Toselli, O. B. Widlund, Domain Decomposition Methods — Algorithms and Theory, volume 34 of Springer Series in Computational Mathematics, Springer Berlin Heidelberg, Berlin, Heidelberg, 2005. doi:10.1007/b137868.
  * Klerk et al. [2008] D. D. Klerk, D. J. Rixen, S. N. Voormeeren,  General Framework for Dynamic Substructuring: History, Review and Classification of Techniques,  AIAA Journal 46 (2008) 1169–1181.
  * Farhat and Roux [1991] C. Farhat, F.-X. Roux,  A method of finite element tearing and interconnecting and its parallel solution algorithm,  International Journal for Numerical Methods in Engineering 32 (1991) 1205–1227.
  * Guyan [1965] R. J. Guyan,  Reduction of stiffness and mass matrices,  AIAA Journal 3 (1965) 380–380.
  * He and Fu [2001] J. He, Z.-F. Fu, Modal Analysis, Elsevier, 2001. doi:10.1016/B978-0-7506-5079-3.X5000-1.
  * Craig and Bampton [1968] R. Craig, M. Bampton,  Coupling of Substructures for Dynamic Analyses 6 (1968) 1313–1319.
  * Rubin [1975] S. Rubin,  Improved Component-Mode Representation for Structural Dynamic Analysis,  AIAA Journal 13 (1975) 995–1006.
  * Pichler et al. [2017] F. Pichler, W. Witteveen, P. Fischer,  Reduced-Order Modeling of Preloaded Bolted Structures in Multibody Systems by the Use of Trial Vector Derivatives,  Journal of Computational and Nonlinear Dynamics 12 (2017) 051032.
  * Blockmans et al. [2015] B. Blockmans, T. Tamarozzi, F. Naets, W. Desmet,  A nonlinear parametric model reduction method for efficient gear contact simulations,  International Journal for Numerical Methods in Engineering 102 (2015) 1162–1191.
  * Balajewicz et al. [2015] M. Balajewicz, D. Amsallem, C. Farhat,  Projection-based model reduction for contact problems,  International Journal for Numerical Methods in Engineering 106 (2015) 644–663.
  * Géradin and Rixen [2016] M. Géradin, D. J. Rixen,  A ‘nodeless’ dual superelement formulation for structural and multibody dynamics application to reduction of contact problems,  International Journal for Numerical Methods in Engineering 106 (2016) 773–798.
  * Mehrdad Pourkiaee and Zucca [2019] S. Mehrdad Pourkiaee, S. Zucca,  A Reduced Order Model for Nonlinear Dynamics of Mistuned Bladed Disks with Shroud Friction Contacts,  Journal of Engineering for Gas Turbines and Power 141 (2019) 1–13.
  * Ghavamian et al. [2017] F. Ghavamian, P. Tiso, A. Simone,  POD–DEIM model order reduction for strain softening viscoplasticity,  Computer Methods in Applied Mechanics and Engineering 317 (2017) 458–479.
  * Wu et al. [2019] L. Wu, P. Tiso, K. Tatsis, E. Chatzi, F. van Keulen,  A modal derivatives enhanced Rubin substructuring method for geometrically nonlinear multibody systems,  Multibody System Dynamics 45 (2019) 57–85.
  * Wu et al. [2018] L. Wu, P. Tiso, F. van Keulen,  Interface Reduction with Multilevel Craig–Bampton Substructuring for Component Mode Synthesis,  AIAA Journal (2018) 1–15.
  * Vizzaccaro et al. [2020] A. Vizzaccaro, L. Salles, C. Touzé,  Comparison of nonlinear mappings for reduced-order modelling of vibrating structures: normal form theory and quadratic manifold method with modal derivatives,  Nonlinear Dynamics (2020).
  * Jain et al. [2018] S. Jain, P. Tiso, G. Haller,  Exact nonlinear model reduction for a von Kármán beam: Slow-fast decomposition and spectral submanifolds,  Journal of Sound and Vibration 423 (2018) 195–211.
  * Ponsioen et al. [2020] S. Ponsioen, S. Jain, G. Haller,  Model reduction to spectral submanifolds and forced-response calculation in high-dimensional mechanical systems,  Journal of Sound and Vibration 488 (2020) 115640.
  * Perez et al. [2017] R. Perez, G. Bartram, T. Beberniss, R. Wiebe, S. M. Spottswood,  Calibration of aero-structural reduced order models using full-field experimental measurements,  Mechanical Systems and Signal Processing 86 (2017) 49–65.
  * Touzé et al. [2014] C. Touzé, M. Vidrascu, D. Chapelle,  Direct finite element computation of non-linear modal coupling coefficients for reduced-order shell models,  Computational Mechanics 54 (2014) 567–580.
  * Hollkamp and Gordon [2008] J. J. Hollkamp, R. W. Gordon,  Reduced-order models for nonlinear response prediction: Implicit condensation and expansion,  Journal of Sound and Vibration 318 (2008) 1139–1153.
  * Kuether et al. [2015] R. J. Kuether, B. J. Deaner, J. J. Hollkamp, M. S. Allen,  Evaluation of geometrically nonlinear reduced-order models with nonlinear normal modes,  AIAA Journal 52 (2015) 3273–3285.
  * Amabili [2013] M. Amabili,  Reduced-order models for nonlinear vibrations, based on natural modes: the case of the circular cylindrical shell,  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 371 (2013) 20120474.
  * Mignolet et al. [2013] M. P. Mignolet, A. Przekop, S. A. Rizzi, S. M. Spottswood,  A review of indirect/non-intrusive reduced order modeling of nonlinear geometric structures,  Journal of Sound and Vibration 332 (2013) 2437–2460.
  * Jain [2015] S. Jain, Model Order Reduction for Non-linear Structural Dynamics, 2015. URL: http://resolver.tudelft.nl/uuid:cb1d7058-2cfa-439a-bb2f-22a6b0e5bb2a.
  * Lu et al. [2019] K. Lu, Y. Jin, Y. Chen, Y. Yang, L. Hou, Z. Zhang, Z. Li, C. Fu,  Review for order reduction based on proper orthogonal decomposition and outlooks of applications in mechanical systems,  Mechanical Systems and Signal Processing 123 (2019) 264–297.
  * Noor and Peterst [1980] A. K. Noor, J. M. Peterst,  Reduced Basis Technique for Nonlinear Analysis of Structures,  AIAA Journal 18 (1980) 455–462.
  * Idelsohn and Cardona [1985] S. R. Idelsohn, A. Cardona,  A reduction method for nonlinear structural dynamic analysis,  Computer Methods in Applied Mechanics and Engineering 49 (1985) 253–279.
  * Sombroek et al. [2018] C. S. Sombroek, P. Tiso, L. Renson, G. Kerschen,  Numerical computation of nonlinear normal modes in a modal derivative subspace,  Computers and Structures 195 (2018) 34–46.
  * Jain et al. [2017] S. Jain, P. Tiso, J. B. Rutzmoser, D. J. Rixen,  A quadratic manifold for model order reduction of nonlinear structural dynamics,  Computers and Structures 188 (2017) 80–94.
  * Benner et al. [2015] P. Benner, S. Gugercin, K. Willcox,  A Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems,  SIAM Review 57 (2015) 483–531.
  * Baur Peter Benner Bernard Haasdonk Christian Himpe Immanuel Maier Mario Ohlberger and Dynamik Komplexer [2017] U. Baur Peter Benner Bernard Haasdonk Christian Himpe Immanuel Maier Mario Ohlberger, F. Dynamik Komplexer, Comparison of methods for parametric model order reduction of instationary problems, 2017. URL: http://www.mpi-magdeburg.mpg.de/preprints/.
  * Xiao et al. [2015] D. Xiao, F. Fang, A. G. Buchan, C. C. Pain, I. M. Navon, A. Muggeridge,  Non-intrusive reduced order modelling of the Navier-Stokes equations,  Computer Methods in Applied Mechanics and Engineering 293 (2015) 522–541.
  * Xiao et al. [2017] D. Xiao, F. Fang, C. C. Pain, I. M. Navon,  A parameterized non-intrusive reduced order model and error analysis for general time-dependent nonlinear partial differential equations and its applications,  Computer Methods in Applied Mechanics and Engineering 317 (2017) 868–889.
  * Xiao [2019] D. Xiao,  Error estimation of the parametric non-intrusive reduced order model using machine learning,  Computer Methods in Applied Mechanics and Engineering 355 (2019) 513–534.
  * Zimmermann [2019] R. Zimmermann,  Manifold interpolation and model reduction (2019) 1–36.
  * Barrault et al. [2004] M. Barrault, Y. Maday, N. C. Nguyen, A. T. Patera,  An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations,  Comptes Rendus Mathematique 339 (2004) 667–672.
  * Chaturantabut and Sorensen [2010] S. Chaturantabut, D. C. Sorensen,  Nonlinear Model Reduction via Discrete Empirical Interpolation,  SIAM Journal on Scientific Computing 32 (2010) 2737–2764.
  * Phalippou et al. [2020] P. Phalippou, S. Bouabdallah, P. Breitkopf, P. Villon, M. Zarroug,  ‘On-the-fly’ snapshots selection for Proper Orthogonal Decomposition with application to nonlinear dynamics,  Computer Methods in Applied Mechanics and Engineering 367 (2020) 113120\.
  * Cho et al. [2020] H. Cho, S. J. Shin, H. Kim, M. Cho,  Enhanced model-order reduction approach via online adaptation for parametrized nonlinear structural problems,  Computational Mechanics 65 (2020) 331–353.
  * Kast et al. [2020] M. Kast, M. Guo, J. S. Hesthaven,  A non-intrusive multifidelity method for the reduced order modeling of nonlinear problems,  Computer Methods in Applied Mechanics and Engineering 364 (2020) 112947\.
  * Hesthaven and Ubbiali [2018] J. S. Hesthaven, S. Ubbiali,  Non-intrusive reduced order modeling of nonlinear problems using neural networks,  Journal of Computational Physics 363 (2018) 55–78.
  * Maulik et al. [2020] R. Maulik, B. Lusch, P. Balaprakash,  Reduced-order modeling of advection-dominated systems with recurrent neural networks and convolutional autoencoders,  arXiv (2020).
  * Astolfi [2008] A. Astolfi,  Model reduction by moment matching for nonlinear systems,  in: 2008 47th IEEE Conference on Decision and Control, volume 2015-Febru, IEEE, 2008, pp. 4873–4878. URL: http://ieeexplore.ieee.org/document/7039956/http://ieeexplore.ieee.org/document/4738791/. doi:10.1109/CDC.2008.4738791.
  * Astolfi [2010] A. Astolfi,  Model Reduction by Moment Matching for Linear and Nonlinear Systems,  IEEE Transactions on Automatic Control 55 (2010) 2321–2336.
  * Ionescu and Astolfi [2016] T. C. Ionescu, A. Astolfi,  Nonlinear moment matching-based model order reduction,  IEEE Transactions on Automatic Control 61 (2016) 2837–2847.
  * Rafiq and Bazaz [2020] D. Rafiq, M. A. Bazaz,  A framework for parametric reduction in large-scale nonlinear dynamical systems,  Nonlinear Dynamics 102 (2020) 1897–1908.
  * Acar and Shkel [2008] C. Acar, A. Shkel, MEMS Vibratory Gyroscopes, Springer, 2008.
  * Izadi et al. [2018] M. Izadi, F. Braghin, D. Giannini, D. Milani, F. Resta, M. F. Brunetto, L. G. Falorni, G. Gattere, L. Guerinoni, C. Valzasina,  A comprehensive model of beams’ anisoelasticity in MEMS gyroscopes, with focus on the effect of axial non-vertical etching,  5th IEEE International Symposium on Inertial Sensors and Systems, INERTIAL 2018 - Proceedings (2018) 1–4.
  * Wang et al. [2018a] X. Q. Wang, G. P. Phlipot, R. A. Perez, M. P. Mignolet,  Locally enhanced reduced order modeling for the nonlinear geometric response of structures with defects,  International Journal of Non-Linear Mechanics 101 (2018a) 1–7.
  * Wang et al. [2018b] X. Q. Wang, P. J. O’Hara, M. P. Mignolet, J. J. Hollkamp,  Reduced Order Modeling with Local Enrichment for the Nonlinear Geometric Response of a Cracked Panel,  AIAA Journal 57 (2018b) 421–436.
  * Budiansky [1967] B. Budiansky,  Dynamic Buckling of Elastic Structures: Criteria and Estimates,  in: Proceedings of an International Conference Held at Northwestern University, Evanston, Illinois, Pergamon Press Ltd, 1967. URL: http://linkinghub.elsevier.com/retrieve/pii/B9781483198217500107. doi:10.1016/B978-1-4831-9821-7.50010-7.
  * Wang et al. [2013] X. Wang, S. Cen, C. Li,  Generalized neumann expansion and its application in stochastic finite element methods,  Mathematical Problems in Engineering 2013 (2013).
  * Hay et al. [2010] A. Hay, J. Borggaard, I. Akhtar, D. Pelletier,  Reduced-order models for parameter dependent geometries based on shape sensitivity analysis,  Journal of Computational Physics 229 (2010) 1327–1352.
  * Tiso [2011] P. Tiso,  Optimal second order reduction basis selection for nonlinear transient analysis,  in: Proceedings of the 29th IMAC A Conference on Structural Dynamics 2011, volume 3, 2011, pp. 27–39. URL: http://link.springer.com/10.1007/978-1-4419-9299-4http://link.springer.com/10.1007/978-1-4419-9299-4{_}3. doi:10.1007/978-1-4419-9299-4_3.
  * Krack and Gross [2019] M. Krack, J. Gross, Harmonic Balance for Nonlinear Vibration Problems, 2019\. doi:10.1007/978-3-030-14023-6.
  * Jutho et al. [2019] Jutho, getzdan, S. Lyon, M. Protter, M. P. S, Leo, J. Garrison, F. Otto, E. Saba, D. Iouchtchenko, A. Privett, A. Morley, Jutho/tensoroperations.jl: v1.1.0, 2019. URL: https://doi.org/10.5281/zenodo.3245497. doi:10.5281/zenodo.3245497.
  * Woiwode et al. [2020] L. Woiwode, N. N. Balaji, J. Kappauf, F. Tubita, L. Guillot, C. Vergez, B. Cochelin, A. Grolet, M. Krack,  Comparison of two algorithms for Harmonic Balance and path continuation,  Mechanical Systems and Signal Processing 136 (2020) 106503.