## Abstract

The development of reduced-order models remains an active research area, despite advances in computational resources. The present work develops a novel order-reduction approach that is designed to incorporate isolated regions that contain, for example, nonlinearitites or accumulating damage. The approach is designed to use global modes of the overall system response, which are then naturally coupled to the response in the isolated region of interest. Two examples are provided to demonstrate both the accuracy and the computational efficiency of the proposed approach. The performance of this approach is compared to the exact response corresponding to a finite element simulation for the chosen problems. In addition, the accuracy and computational efficiency are shown relative to a standard Galerkin reduction based on the linear normal modes. It is found that the proposed reduction offer computational efficiency comparable to a Galerkin reduction, but more accurately represents the response of the system when both are compared to the finite element simulation.

## 1 Introduction

Despite recent advances in computing power, large-scale finite element analysis (FEA) can still present computational difficulties simply arising from, for example, the overall size of the problems considered or large collections of analyses required for uncertainty quantification. Therefore, reduction techniques originating when computational resources were a fraction of current capabilities continue to retain their importance and usefulness. Component mode synthesis (CMS), based on the concept of substructuring, is one of the most widely used reduction methods in structural dynamics. Originally based on the the work of Hurty [1] and extended by Craig and Bampton [2], CMS decomposes the response of an individual component of a larger structure into component and constraint modes [3,4]. The dominant modes of the component can be identified and retained, thus significantly reducing the overall order of the system, while the constraint modes can be used to couple adjacent components together. Such an approach is particularly attractive when the response of the overall structure is localized to only a few of the components of the system, but can be used to develop reduced-order models of the response of complex structural systems.

However, there are numerous examples where the response of the overall system can be characterized as global in nature, so that the response of the system is more closely represented by the structural modes of the overall system rather than localized within only a few components. Jointed structures provide a canonical example of such systems [5–10], where the response of the overall system with joints is closely related to the modes of the overall structure in the absence of the interfaces, known as the monolithic system. In such a case, the structural modes of the monolithic system are significantly different from the vibration modes of the individual components. As identified by Segalman [11–13], in such cases, it is desirable to describe the response and ultimately develop reduced-order models in terms of the structural modes of the overall system, rather than a collection of component-level modes. However, in such a system, the dynamics introduced by the presence of the interfaces, which are internal to the structure, nonetheless significantly influence the response of the overall system [6].

More generally, the response of system with isolated features, including nonlinearitites, could in principle be represented by CMS, so that the isolated effects are introduced through coupling between the constraint modes, while the exterior linear regions were decomposed into their individual component modes. However, the same issue remains—the response of the overall system may be better represented by the modes of the overall system rather than those of the individual (linear) components. In the present work, we develop a novel order-reduction method for systems with isolated nonlinearitites. In contrast to CMS, which utilizes the modes of the individual components to achieve the reduction in the overall model order, the present method employs the modes of the global system, obtained from the corresponding monolithic system in the absence of the nonlinearitites.

## 2 Global Model Reduction

*reference system*. A mixed formulation for the equations of motion, similar to that presented in Quinn [14], is developed here but extended to a form appropriate for the equations of motion arising from a finite element discretization of the equations of motion for a general structural system.

### 2.1 Mixed Formulation.

*N*degree-of-freedom mechanical system given in Eq. (1) can be expanded in the form

The closure of this system requires the deviatoric response $(z_\alpha ,z_\beta )$ in the boundary regions be determined outside of this system of equations. In addition, Eq. (10) cannot be used to determine the response within the interior region defined by the boundary and isolated regions $(C\alpha ,Cc,C\beta )$. However, if $(z_\alpha ,z_\beta )$ can be determined then the resulting response in the exterior regions is exact. In Ref. [14], the deviatoric force $Q_$ was specified as a constitutive equation, based on a reduced-order model for an interface as developed by Ref. [15].

*unreduced isolated system*.

In summary, the nonlinear equations of motion for the reference system given in Eq. (1) are written in terms of $x_\u2208RN$, with *N* = *N*_{1} + *N*_{r} + *N*_{2}. These have been decomposed into a set of linear equations in terms of the mixed variable $w_\u2208RN$ given in Eq. (10) together with a system of nonlinear equations given in Eq. (13), written in terms of $z_\u2208RNr$. At first glance, this appears to be a step backwards, as the original *N* degree-of-freedom system has been replaced by a system with *N* + *N*_{r} degrees-of-freedom. However, the linear equations for $w_$ may be able to be reduced using existing model reduction techniques, so that if *N* ≫ *N*_{r}, a significant computational reduction may be achieved and, as illustrated below, a significant increase in accuracy as compared to traditional Galerkin reduction in terms of the eigenmodes of the linear system.

### 2.2 Modal Analysis.

^{2}. However, such a eigenanalysis can be performed in the linear system with $N_\u22610_$. More generally, various order reduction techniques could be applied to this linear system, including projection based methods [16]. Moreover, assuming proportional damping, the undamped modes $\varphi _$ defined as

*q*

_{k}(

*t*) the

*k*

^{th}modal coordinate, defined as

*δ*

_{mn}is the Kronecker delta function. Other expansions, including those based on Proper Orthogonal Decompositions [17–19], could be employed as well. In this, $\varphi _kTf_(t)$ describes the modal forcing while $\varphi _kTQ_$ describes the effect of the isolated nonlinearitites, both acting on the

*k*

^{th}mode. However, as described above, the deviatoric forcing $Q_$ must be determined outside of Eq. (10). In particular, the deviatoric displacement $z_$ is driven by the mixed response of the interior region as defined by

*M*mode truncation then the deviatoric response in the interior region as defined by $z_$ can be determined from Eq. (13), which can then be used to determine the reduced-order modal response of the system as given in Eq. (15). With the modal reduction, the resulting system has

*M*+

*N*

_{α}+

*N*

_{c}+

*N*

_{β}degrees-of-freedom. While

*M*≤

*N*, typically the retained modes offer a significant reduction in the overall degree-of-freedom for the system, so that

*M*≪

*N*

_{1}+

*N*

_{2}. This formulation nonetheless couples together all of the retained modes through the reconstruction of the mixed response in the interior region from the modal expansion, as determined from Eq. (18). The resulting response of the system can then be reconstructed from the modal coordinates $q_(t)$ and the interior coordinates $z_(t)$ as

## 3 Example: Finite Element Analysis Rod Elements

*m*and ℓ, respectively, while the axial rigidity is

*EA*, then for the the uniform rod with

*N*elements the corresponding mass and stiffness matrices are

*N*. The spatial position of the rod, denoted

*s*lies in the interval 0 ≤

*s*≤ ℓ, and if the interior region is located in the interval

*s*∈ (

*s*

_{1},

*s*

_{2}), then the boundary nodes are assumed to correspond to $\alpha =\u230as1N\u2113\u230b$ and $\beta =\u2308s2N\u2113\u2309$, so that

*i*and

*j*in the isolated region are assumed to depend on the relative displacement between node pairs, defined as

*v*

_{i,j}(

*t*) ≡

*x*

_{i}(

*t*) −

*x*

_{j}(

*t*), so that individual forces are assumed to take the form

*h*[

*v*(

*t*)] is a hysteretic operator acting on the relative displacement

*v*(

*t*). These additional terms used in the interior region are representative of nonlinearities that are often seen in structural systems, including linear detuning and cubic stiffness and damping terms, as well as a hysteretic force commonly arising from an interface [20]. These terms were chosen to highlight the ability of the proposed approach to incorporate such common terms. Therefore, the nonlinear force acting between nodes

*i*and

*j*is defined as

*h*[

*v*(

*t*)] is a force component arising from hysteretic damping and is determined by a differential constitutive equation. In the case of a Bouc model for the hysteretic damping, the equations of motion are augmented by constitutive relations of the form

*N*degrees-of-freedom, the system is simulated with the nominal parameters

*ρ*, ℓ,

*EA*) do not restrict the system dynamics. In the internal region the parameters associated with the nonlinearities are scaled by the number of elements

*N*, so that the equivalent stiffness and damping across the internal interval

*s*∈ (

*s*

_{1},

*s*

_{2}) = (0.25, 0.35) remains (approximately) constant as the number of elements varies.

*N*= 80, and for this discretization, the number of elements in each region is

For this comparison, the initial conditions of the system are chosen so that the rod is statically displaced by a unit force applied at *s* = 0.30 and then released from rest. Note that this initial state excites every mode in the structure, and the response of the system at *s* = 0.50 is shown in Fig. 2. In Fig. 2(a), the response of the system is shown for the original equations of motion, while in Fig. 2(b), the response is shown at the same location for the modally reduced isolated model, retaining *M* = 10 structural modes. Note that this system has nine degrees-of-freedom associated with the interior region (c.f., Eq. (12)), so that the total degree-of-freedom of the reduced model is 19. Finally, for comparison, the response for a 19 degrees-of-freedom Galerkin reduction is shown in Fig. 2(c).

In Fig. 3(a) the error as defined by Eq. (28) is shown for both the localized reduced-order model as well as the corresponding Galerkin reduction. Although not shown, when the system is simulated with the unreduced isolated formulation given in Eqs. (10) and (13), the response is identical to the simulation with the original equations up to numerical accuracy of the integrator used. The nonlinearities are localized in the interior region *s* ∈ (*s*_{1}, *s*_{2}) = (0.25, 0.35). As a result, accurate resolution of the dynamics in this region is critical to capture the effect of the nonlinearities, including the energy loss due to dissipation. In Fig. 3(b), the error in the response of the interior region is shown as a function of time for the modally reduced isolated model as well as for the Galerkin reduction, both retaining 19 total degrees-of-freedom. As illustrated, the isolated model is more accurate in the interior region as compared to the Galerkin reduction.

*M*, the number of structural modes retained, varies, while in Figs. 4(b), the computational efficiency is shown, defined as the ratio of the computational time of the reduced system relative to that of the original equations of motion, that is

*T*is measured using the matlab’s built-in toc function. As observed in Fig. 4(a), there is a minimum degree-of-freedom required for the modally reduced isolated model due to the fixed number of degrees-of-freedom associated with the interior region, which with

*N*= 80 is 9. The error in the reduced isolated system quickly decreases as the number of retained spatial modes increases, while the error in the Galerkin reduced model decreases more slowly. However, because of the additional degrees-of-freedom associated with the interior region in the reduced isolated system, it requires a minimum number of spatial modes to match the accuracy of the Galerkin approximation, in this case with

*M*= 6. Above this spatial resolution, however, the reduced isolated model with retaining

*M*> 6 spatial modes is more accurate than the corresponding Galerkin approximation with

*M*+ 9 spatial modes. Moreover, as seen in Fig. 4(b), the computational effort associated with the reduced isolated model is comparable to that of the Galerkin reduction of the same size. Note that the simulation uses the matlab adaptive integrator ode23t, so that the computational time depends on a number of factors in additional to model size, including internal integrator tolerances and estimated numerical error.

## 4 Example: Finite Element Analysis Beam Elements

*E*

*I*and

*ρA*are considered to be unity and the system is modeled using Euler-Bernoulli beam elements. Each node contains both rotational and transverse degrees-of-freedom, so an

*P*node discretization has

*N*= 2

*P*degrees-of-freedom. In addition, proportional damping is used with

*ζ*= 0.02. The system is assumed to have contain cubic stiffness and damping nonlinearities over the interval

*s*∈ (0.40, 0.50). External transverse loading is applied to the beam locations at

*s*=

*s*

_{1}= 0.33 and

*s*=

*s*

_{2}= 0.67, of the form

*a*)

*b*)

*ω*= 3.516 of the beam. Finally, the beam is initially at rest in its undeformed position.

With *P* = 35 nodes for the system, the resulting system has *N* = 70 degrees-of-freedom. As illustrated in Fig. 5(a), the average error as defined in Eq. (28) depends on the number of modes retained in the expansion of the structural system, and as *M* increases the error decreases, while as shown in Fig. 5(b), the computational time ratio is reduced for *M* < *N*, and approaches unity as *M* → *N*.

## 5 Conclusions

A novel order-reduction approach has been developed specifically for systems that contain isolated features, including nonlinearitites. The approach assumes linear exterior regions connected by localized regions that could, for example, include nonlinearities such as those arising from joints and interfaces. In contrast to a familiar component mode synthesis, the overall system is reduced using the global modes of the system, which then couple to the deviatoric effects introduced by the isolated region. The result is a reduced-order modeling approach that accurately represents the response in the isolated region and therefore can accurately describe the nonlinearities, yet is nonetheless written in terms of the global modes of the system. Two examples are considered, the transient response of a longitudinal rod and the forced response of an Euler-Bernoulli beam. The former is compared against both a standard FEA analysis as well as a Galerkin expansion using the linear mode shapes; the latter is compared against FEA analysis. Note however, that the response within the isolated region as described by $z_$ has not been reduced. Future efforts can seek to develop appropriate reduction methods focused on this interior region, that can then be combined with the framework presented here.

## Footnote

Note that $N_$ need not necessarily be nonlinear, but instead could, for example, include some linear loss of stiffness due to increasing damage. While in this instance a modal decomposition could in principle be undertaken in the damaged state, it is not utilized here. In the present work any additional terms in $N_$ will be referred to as “nonlinear” regardless of the specific form of the terms.

## Acknowledgment

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE–NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper. Data are provided by a third party listed in Acknowledgment.

## References

*Dynamics and Bifurcations of Non-Smooth Mechanical Systems*