Computer simulations are routinely performed in today’s technological world for modeling and response prediction of almost any physical phenomenon. The ever-increasing demand for realistic simulations leads to a higher level of detail during the modeling phase, which in turn increases the complexity of the models and results in a bigger problem size. Typically, such physical processes are mathematically modeled using partial differential equations (PDEs), which are discretized (e.g. using Finite Elements, Finite differences, Finite volumes methods etc.) to obtain problems with a finite (but usually large) number of unknowns. Despite the tremendous increase in computational power over the past decades, however, the time required to solve high-dimensional discretized models remains a bottleneck towards efficient and optimal design of structures. Model order reduction (MOR) aims to reduce the computational efforts in solving such large problems.
The classical approach to model reduction involves a linear projection of the full system onto a set of basis vectors. This linear projection is characterized by a matrix whose columns span a suitable low-dimensional subspace. Various techniques have been applied to high-dimenisional systems to obain such a reduction basis, including the Proper Orthogonal Decomposition (POD)POD1 ; POD2 ; POD3
(also knowns as Singular Value Decomposition (SVD), Karhunen-Loeve Decomposition), Linear Normal Modes (LNM) and Krylov subspace projectionkrylov . Once a suitable basis is chosen, the reduced-order model (ROM) is then obtained using Galerkin projection. Similar linear projection techniques have been devised for component-mode synthesis (CMS), such as the Craig-Bampton method Craig Bampton . An implicit assumption to all linear projection techniques is that the full system dynamics evolves in a lower-dimensional linear invariant subspace of the phase space of the system. While such linear subspaces do exist for linear systems (linear modal subspaces), they are generally non-existent in nonlinear systems. This results in a priori unknown and potentially large errors for linear projection-based reduction methods, necessitating the verification of the accuract of the reduction procedure on a case-by-case basis.
More recent trends in model reduction account for this issue by constructing the reduced solution over nonlinear manifolds qm ; qm2 ; NLDR . The seminal idea of Shaw and Pierre shawpierre93 is to construct assumed nonlinear invariant surfaces (nonlinear normal modes) that act as continuations of linear normal mode families under the addition of nonlinear terms near an equilibrium (see kerschen2014
for a recent review). More heuristic reduction procedures in structural dynamics include the static condensation approaches, where the fast or stiff variables in the system are intuitively identified, and statically enslaved to the slow or flexible ones (cf.mignolet for a review). Guyan reduction Guyan is a classic example of this reduction philosophy at the linear level.
Most classic reduction techniques, coming from a intuive and heuristic standpoint, do not provide an a priori estimate of their accuracy or even validity for a given system. To this effect, Haller and Ponsioenhaller16sf proposed requirements for mathematically justifiable and robust model reduction in a non-linear mechanical system. These requirements ensure not only that the lower dimensional attracting structure (manifold) in the phase space is robust under perturbation, but also that the full system trajectories are attracted to it and synchronize with the reduced model trajectories at rates that are faster than typical rates within the manifold.
Recent advances in nonlinear dynamics enable such an exact model-reduction using the slow-fast decomposition (SFD) haller16sf and spectral submanifold (SSM) based reduction haller16 . SFD is a general procedure to identify if a mechanical system exhibits a global partitioning of degrees of freedom into slow (flexible) and fast (stiff) components such that the fast variables can be enslaved to the slow ones. This results in a global ROM containing only the slow degrees of freedom of the full system. SSMs, on the other hand, are the smoothest nonlinear extensions of the linear modal subspaces near an asymptotically stable equilibrium. Neither SFD nor SSM-based reduction has been applied to problems with high numbers of degrees of freedom. A beam model is often the first step in showing the potential of a new technique for reduction of high dimensional systems (cf. ilbeigi16 ; craig89 ; nnmbeam ). To this effect, we combine here, for the first time, the application of these techniques on a finite-element discretized nonlinear von Kármán beam model.
We first show that the beam model satisfies the requirements of SFD. The corresponding ROM is subsequently obtained on a slow manifold defined over the transverse degrees of freedom of the beam. This SFD-reduced ROM posseses no clear spectral gaps but a spectral quotient analysis nonetheless enables a further reduction to an SSM using the formulas given by Szalai et. al Szalai16 . Importantly, our spectral quotient analysis returns the number of modes relevant for reduction, instead of postulating or deducing this number from numerical experimentation. In the end, our two-step reduction results in an exact ROM with a drastic reduction in the number of degrees of freedom of the system.
In the next section, we start by reviewing the main steps involved in the derivation of the governing PDEs for the von-Karman beam and the non-dimensionalization we performed upon them. The finite-element disretized equations obtained from the resulting nondimensionalized PDEs are then presented. This system of equations is then first reduced using the SFD in Section 3. The SSM-based reduction applied to the SFD-reduced system in Section 4. The conclusions, along with scope for further work, are presented in Section 5.
|Length of beam (m)|
|Height of beam (m)|
|Width of beam (m)|
|Area of cross section (m) =|
|Young’s Modulus (Pa)|
|Viscous damping rate of material (Pa s)|
We briefly summarize the main steps leading to the derivation of the governing partial differential equations (PDEs) for the von Kármán beam (see, e.g., Reddy reddy10 for a detailed derivation). We consider a straight 2D beam aligned initially with the axis, as shown in Figure 1. The motion of the beam takes place in the plane.
Assuming the Euler-Bernoulli hypothesis for the kinematics of bending (i.e., that the lines initially straight and perpendicular to the beam axis remain so after deformation), we obtain the displacement field:
where denote the displacements of a material point with coordinates , and are the displacements of any material point lying on the reference line given by the axis. The von Kármán strain approximation of the Green-Lagrange strain for moderate rotations is given by
Using the virtual work principle, we can formulate the equations of motion in terms of the primary unknowns as
where and are external body forces per unit length of the beam in the and directions, respectively; is the area of the cross-section; and is the component of the Cauchy stress. We choose the Kelvin-Voigt model for viscoelasticity as constitutive law to relate the stress to the von Kármán strain (2) and to the corresponding strain rate as
Here denotes Young’s modulus and is the rate of viscous damping for the material. Further assumptions include a uniform rectangular cross section with the reference axis passing through the centroid of the cross-section, and a non-dimensionalization of variables as ( have been introduced as scaling factors for transverse and axial loading respectively). These lead to the following dimensionless PDEs governing the beam behavior:
Here denotes , is the length to thickness ratio and is a dimensionless constant resulting from non-dimensionalization. We have chosen the displacements to be scaled with respect to the thickness of the beam, and not the length. Moreover, the applied loading has been nondimensionalized and scaled with respect to a load which leads to transverse displacements in the order of the thickness of the beam. This is because the von Kármán kinematic approximations are justified for displacements and forces in this range. Furthermore, the time has been nondimensionalized with respect to a time period that is representatitve of natural frequency of the oscillation of the beam, since the structural response is generally studied at such time scales.
Upon finite-element discretization of the non-dimensional system (3) with cubic shape functions for and linear shape functions for (see, e.g., Crisfield crisfield96 ), we obtain the finite-dimensional discretized version of (3) as
where are the finite dimensional (discretized) counterparts of the unknowns respectively ( being the number of unknowns dependent on the finite-element discretization), and and are the corresponding mass and stiffness matrices. Here, (a bilinear function of the form ), (a cubic function of the form ) and (a quadratic function of the form ) correspond to the nonlinear elastic forces in the beam; (a quadratic function of the form ), (a linear function of the form ) and (a linear function of the form ) correspond to the nonlinear contributions resulting from the viscoelastic material damping.
3 Slow-Fast Decomposition
3.1 Verification of assumptions for the application of SFD
In the finite-element discretized system (4), the variables (representing the axial displacement components) are stiffer and hence faster than the variables (representing the transverse displacement components). Such a global difference of speeds indicates the possible existence of a lower-dimensional slow manifold, as described in mathematical terms by the geometric singular perturbation theory of Fenichel fenichel79 . If such a slow manifold is robust and attracts nearby solutions, then the dynamics on this manifold provides and exact reduced-order model with which all nearby solutions synchronize exponentially fast.
For general finite-dimensional mechanical systems characterized by such a dicotomy of time scales, Haller and Ponsioen haller16sf deduced conditions under which positions and velocities in the fast degrees of freedom (,) can be expressed as a graph over their slow counterparts (,), resulting in a globally exact model reduction. If these conditions for a Slow-Fast Decomposition (SFD) are satisfied, then all trajectories of the full system (close enough to the slow manifold in the phase space) synchronize with the reduced model trajectories at rates faster than those within the slow manifold. To check these conditions, we take , the thickness to length ratio of the beam, as the required non-dimensional small parameter. The system (4) can then be reformulated in terms of as
The mass-normalized forcing terms are defined in terms of the new variable as
The main conditions of the SFD procedure, as derived by Haller and Ponsioen haller16sf , are the following.
Nonsingular extension to : and should possess smooth (in fact ) extension to their respective limits, which is the case here by (6).
Existence of a critical manifold: The algebraic equation should be solvable for on an open, bounded domain. Indeed, for any open and bounded, the set defined by
satisfies for all , forming a critical manifold in the language of singular perturbation theory. Clearly, the critical manifold is independent of in our current setting. If, however, we assumed that (i.e., the beam is excited in the axial direction by forces that are an order of magnitude larger than those in the transverse direction), then we would obtain a time-dependent critical manifold.
For clarity, we denote any expression evaluated over the critical manifold by
Asymptotic stability of the critical manifold: With the matrices
the equilibrium solution of the unforced, constant-coefficient linear system
should be asymptotically stable for all fixed parameter values . This is again satisfied in our setting since are positive definite matrices and .
3.2 Global reduced-order model from SFD
As shown by Haller and Ponsioen haller16sf , assumptions (A1)-(A3) guarantee that the critical manifold perturbs into a nearby attracting slow manifold for small enough. On this slow manifold, the discretized beam system (4) admits an exact reduced order model given by
constitute the higher-order terms in the equations describing the slow manifold as
In the context of the beam example considered here, the exact reduced-order model in (9) can be explicitly written out in the following form:
Note that the reduced order model (11) is conservative (contains only inertial and elastic force terms) at leading order, whereas the full system (4) is dissipative due to viscoelastic damping. The terms in the ROM account for these damping contributions and hence cannot be ignored. Apart from the damping contributions, there also exists a conservative correction at the level which includes the static response of the variables to the corresponding loading (cf. the expression for in (12)).
3.3 Specific results from SFD
We now consider a specific beam with geometric and material parameters as follows: length m; thickness to length ratio in the range from ; Young’s modulus GPa; density Kg/m; material viscous damping rate Pa s. We use a spatially uniform load on the beam in the axial and in the transverse direction, given by , , . Here is the constant used to nondimensionalize time and is the loading frequency (chosen to be the first natural frequency of the beam in this case, cf. Figure 3). Using these parameters, we obtain .
The graphs presented above show the response of a single degree of freedom in time. In order to check the performance of SFD globally, we use the following measure,
where is the full vector of generalized displacements at time , obtained from the full nonlinear solution; is the solution based on the reduced model; and is the set of time instants at which the error is recorded. The error recorded for different cases is shown in Table 2. We observe that for of order , the ROM at the leading order (which is conservative) provides a close approximation for the full system. The terms, however, become important when is of order , because they account for damping which becomes more significant in this-range. We further calculated the terms in the ROM to compare changes in accuracy (see A for general calculation of the terms). We observe that the terms only provide a marginal improvement in terms of accuracy in our case. A further increase in values to the order of renders the reduced model highly inaccurate, thus higher-order terms in are necessary to achieve the desired accuracy. Note, however, the values around are rather unphysical; the parameter , denoting the thickness-to-length ratio of a beam, is not expected attain such high values.
|SFD -||SFD - up to||SFD - up to|
4 SSM reduction
Though the SFD reduction is robust and globally valid in the phase space, the reduced system on the slow manifold still contains degrees of freedom, which is generally a large number. Specifically, due to the choice of the shape functions for discretization, the reduced model (11) obtained from SFD still contains two-thirds of the total number of unknowns of the full system (4
). The following eigenvalue analysis of the linearized reduced system, however, reveals the existence of a further separation in time scales within the slow manifold near the origin, showing potential for a further reduction via spectral submanifolds (cf. Haller and Ponsioenhaller16 )
4.1 Separation in damping: spectral quotient analysis
The eigenvalue problem for the linearization of the reduced system (11) can be formulated as
where and are the
eigenvalue and eigenvector for any. The negative real part of the eigenvalue
represents the exponential rate of decay of trajectories towards the equilibrium position along the corresponding two-dimensional eigenspace of the linearized system. The smoothest local extension of such an invariant subspace is known as the spectral submanifold (SSM), a notion introduced by Haller and Ponsioenhaller16 . Tangent to a slow subspace (i.e., the subspace spanned by eigenvectors corresponding to eigenvalues with the lowest-magnitude real parts), such an SSM offers an opportunity for a further, drastic reduction from the SFD-based slow manifold to a two-dimensional invariant manifold.
Estimating the real parts of the eigenvalues arising from (14) is sensitive to the choice of the eigensolver algorithm, especially for large values. An accurate first-order approximation to the eigenvalues, however, can be obtained using the undamped eigenvalues from their defining equation
Based on this formula, the modal frequencies, real parts of the eigenvalues, and the ratios of the subsequent real parts are shown as a function of the mode number in Figure 6.
Note that the first eigenspace is about times slower in terms of its exponential decay rate than the second slowest eigenspace, which renders the first mode an optimal choice for reduction using its corresponding SSM. We will, therefore, use this mode to perform a single-mode SSM reduction for the beam model, as discussed below.
4.2 SSM-based reduction
In order to perform a further reduction of the reduced system (11) obtained from SFD, we now compute the slowest single-mode SSM, whose reduced dynamics is given by a two-dimensional ODE as a final reduced-order model. Szalai et al. Szalai16 have obtained general formulae for such one-mode SSMs and their associated backbone curves for autonoumous dynamical systems (no external forcing in our setting). The reduced-order model obtained in this fashion is given in terms of an internal parametrization of the SSM, rather than a projection of the SSM on its underlying modal subspace. This construct allows the SSM to be recovered more globally, even if it develops a fold over the underlying modal subspace.
For this general formulation, we require the autonomous counterpart of system (11) to be in the phase space form
where is a class function representing the nonlinear terms in the ROM (11). By inspection of (11), we find that is a strictly cubic polynomial in and hence belongs to the class of analytic functions. The system (16) can be diagonalized as
where are the eigenvalues of the matrix ordered such that . The matrix
of the linear transformationleading to (17), contains the corresponding eigenvectors of .
The notion of an SSM was introduced by Haller and Ponsioen haller16 , more formally defined as follows:
A spectral submanifold (SSM), , corresponding to a spectral subspace of is an invariant manifold of the dynamical system (17) such that
is tangent to at the origin and has the same dimension as ;
is strictly smoother than any other invariant manifold satisfying (i).
When is a function, then under the assumptions
the relative spectral quotient satisfies ,
there are no resonances up to order between and the rest of the spectrum of ,
the existence and uniqueness of a class SSM is guaranteed by the main theorem of Haller and Ponsioen haller16 , deduced from the more general results of Cabré et al. cabre2003 . The theorem states that the SSM can be viewed as an embedding of an open set into the phase space of system (17), described by a map . Furthermore, there exists a quadratic polynomial mapping , such that the reduced dynamics on the SSM can be written as
where are the coordinates for the mode (for which the SSM construction is performed).
In our beam setting, we obtain from the eigenvalue approximation in (15). This shows that the relative spectral quotient depends on the chosen discretization and monotonically increases with the number of discretization variables. This can be expected physically from a proportionally damped mechanical structure in which the damping for the modes is expected to increase monotonically with the corresponding oscillation frequency. Since here, we obtain that (B1) is satisfied. Furthermore, we assume generic parameter values under which the non-resonance requirement (B2) is satisfied by the beam system. The theorem of Haller and Ponsioen haller16 then applies to the system (17), and the formulae obtained by Szalai et al. Szalai16 yield the coefficients of and the quadratic and cubic Taylor coefficients of . Specifically, let be expanded in indicial notation as
and let be of the polynomial form
where denotes the order terms of .
The first-order terms in this expansion are given by
is an all-zero matrix except for two non-zero entries given by. The component of the quadratic terms can be written as
is an all-zero 3-tensor sincehas no quadratic components. Finally, using the formulae derived by Szalai et al. Szalai16 , we can write the cubic terms as
where is a sparse 4-tensor with nonzero entries given as
as obtained from the general formulae of Szalai et al. Szalai16 , reproduced in the B. Note that the Einstein summation convention has not been followed in the above expressions. We finally express the reduced dynamics in the polar coordinates using a transformation such that , given by
We now consider the SFD-reduced beam system (11) with the physically relevant parameter value , and further reduce it to its slowest, two-dimensional SSM with . Let the modal coordinates be partitioned such that represents the displacement of master modes, and represents the rest of the slow modes obtained previously from the SFD. For the two-dimensional SSM over longer time-scales, free oscillations of the beam are expected to decay towards the origin. To illustrate this over longer time scales, we consider three sets of intial conditions.
Initial condition on the SSM: Starting with an initial condition on the two-dimensional SSM, we observe in Figure (a)a that the full solution indeed stays on the SSM, thereby showing that the SSM is indeed invariant. Since the computed is only a third-order appoximation of the SSM, an initial condition far enough from the the fixed point may not stay over this approximate surface (cf. Figure (b)b). Indeed, a higher-order approximation to would be required to verify the invariance of numerically for larger initial conditions.
Figure 7: The SSM plotted in the modal coordinates . Full solution trajectory is shown for an initial condition taken on the SSM. Modal amplitude of the mode () is plotted against the master modal variables () with . The full solution trajectory (a) stays approximately on the computed SSM, (b) does not stay exactly on the computed SSM for an initialization far enough from the fixed point.
Initial condition off the SSM but still on the slow manifold: We illustrate that the SSM is attracting nearby trajectories within the SFD-based slow manifold by launching a trajectory away from the SSM but still inside the slow manifold . As shown in Figure 8, a trajectory of the full solution which starts at an arbitrarily chosen point in the phase space of system (16), close enough to the computed SSM, quickly converges towards the SSM and synchronizes with the flow on the SSM. Figure 9 shows that this rate of attraction towards the SSM is indeed faster than the typical decay rate within the SSM.
Figure 8: The SSM plotted in the modal coordinates . Full solution trajectory is shown for an initial condition taken off the SSM but still inside the slow manifold. Modal amplitude of the mode () is plotted against the master modal unknows (). The full solution trajectory (a) quickly decays onto the SSM, and (b) synchronises with the dynamics on the SSM.
Initial condition off the slow manifold: Although the SSM is computed for the SFD-reduced system (11), the reduced dynamics on the SSM also captures the asymptotic behavior of the full system by construction. Indeed, as seen in Figure 10, a trajectory initialized off the slow manifold in the phase space of system (4) (thus off the SSM as well) is also attracted towards the SSM, and synchronizes with the corresponding flow on the SSM. This synchronization can be seen in more detail in Figure 11, where the time histories of slow and fast variables are compared between the slow and full solution. As expected, the response of the fast degree of freedom, in particular, shows that the full solution performs rapid oscilations around, and stabilizes on, the reduced solution obtained from SFD. Finally, the full solution approaches the SSM-reduced solution.
Figure 11: The comparison of full and different reduced solutions for (a) slow and (b) fast degrees of freedom. Note the two timescales in the zoom-in for the dynamics of the fast DOF. The full solution quickly decays to the slow manifold (with SFD-reduced solutions) after which it decays to the SSM.
Next, we study the rate of decay of the enslaved variables towards the computed SSM. From Figure 12, we see that the initial decay rate of a trajectory towards the SSM is approximately which is actually dominated by the decay rate suggested by the second slowest eigenvalue, i.e., . The slow dynamics within the SSM is expected to occur at the rate suggested by the first eigenvalue, i.e., . The final decay-rate for the full solution is, however, a decay rate of which is between the slowest and second slowest rates, still being an order of magnitude faster than the dynamics within the SSM, as expected from the underlying theory.
We have demonstrated two exact model reduction techniques, the slow-fast decomposition (SFD) and spectral-submanifold (SSM)-based reduction, on the finite-element model of a von Kármán beam. The SFD enabled us to express the fast axial variables of the beam as a graph over the slow transverse variables, yielding a reduced-order model only in terms of the slow variables. This SFD-reduced model possesses a gap in its spectral quotient disctribution, indicating an opportunity for further reduction using the SSM. Subsequently, we carried out a single-mode SSM-based reduction on the SFD-reduced model. This two-stage, exact reduction procedure resulted in a drastic reduction of the model dimension, and ensured that the full system trajectories synchronized with the reduced model trajectories at rates much faster than typical decay rates within the manifold.
The application of the above mentioned exact reduction techniques to our beam model confirms their potential for truly high-dimensional systems. However, significant work remains to be done to enable application of these techniques to simulations of realistic structures, especially on the computational implementation side. In particular, these techniques, in their present form, require the nonlinear coefficients of the system to be known apriori. For finite-element-based applications, often these nonlinear coefficients are embedded in the software at the element level and their full counterparts are never calculated during a simulation (for reasons of computational efficiency). This is certainly a computational challenge for application of SFD/SSM-based reduction techniques to real-world applications and shall be the focus of our future endeavours.
The current work features a single-mode SSM. This was justified in our beam example, because decay rates in the first mode were exceptionally slower than the rest. For more general damped-mechanical systems, however, we expect higher-dimensional SSMs (i.e., smoothest nonlinear continuations of high-dimensional linear modal subspaces) to be more relevant for reduction. Furthermore, the polynomial expansion of the SSM mapping contains terms only up to the third order in the current work. As also discussed in Section 4.3, higher-order terms would play an important role in capturing trajectories initialized farther from the fixed point. A MATLAB computational toolbox to compute multi-mode SSMs with terms of up to arbitrary order is currently under development (cf. Ponsioen and Haller Ponsioen17 ).
The current work covers the autonomous (non-forced) beam model. Existence and uniqueness results for SSMs are also available for non-autonomous, damped mechanical systems as discussed by Haller and Ponsioen haller16 . For (quasiperiodically) forced systems, the SSM appears as a quasiperiodically deforming invariant surface in the phase space, which is again an ideal tool for model reduction. The computation of such non-autonomous SSMs will find useful application in non-autonomous mechanical systems, and is subject to ongoing work (cf. Breunung and Haller Breunung17 ).
Appendix A terms in the slow manifold
As shown by Haller and Ponsioen haller16sf , the enslavement of the fast variables to the slow ones along the slow manifold is given by the functions
where are as shown in (12) and is given by
In order to obtain the terms in the reduced model, we have calculated the general expressions for as