Robust simulation of buckled structures using reduced order modeling

Lightweight metallic structures are a mainstay in aerospace engineering. For these structures, stability, rather than strength, is often the critical limit state in design. For example, buckling of panels and stiffeners may occur during emergency high-g maneuvers, while in supersonic and hypersonic aircraft, it may be induced by thermal stresses. The longstanding solution to such challenges was to increase the sizing of the structural members, which is counter to the ever present need to minimize weight for reasons of efficiency and performance. In this work we present some recent results in the area of reduced order modeling of post- buckled thin beams. A thorough parametric study of the response of a beam to changing harmonic loading parameters, which is useful in exposing complex phenomena and exercising numerical models, is presented. Two error metrics that use but require no time stepping of a (computationally expensive) truth model are also introduced. The error metrics are applied to several interesting forcing parameter cases identified from the parametric study and are shown to yield useful information about the quality of a candidate reduced order model. Parametric studies, especially when considering forcing and structural geometry parameters, coupled environments, and uncertainties would be computationally intractable with finite element models. The goal is to make rapid simulation of complex nonlinear dynamic behavior possible for distributed systems via fast and accurate reduced order models. This ability is crucial in allowing designers to rigorously probe the robustness of their designs to account for variations in loading, structural imperfections, and other uncertainties.


Introduction
Aircraft structural components operating in complex loading environments can experience various limit state conditions. Examples of these limit state conditions include, thermal buckling due to spatially varying temperature fields, and snap-through dynamic instabilities caused by intense aeroinduced loading. Furthermore, excessive panel deformations could negatively influence the aerodynamic performance. While it is certainly possible to simulate the structural dynamic response under these extreme conditions using finite element models (FEMs), the large computational cost can be prohibitive, especially during the first stages of preliminary design of the aircraft structural components.
The computational challenges discussed in the previous paragraph have led to the development of reduced order models (ROMs) for the computation of the nonlinear structural dynamic response. These ROMs have extended modal analysis methods of linear dynamical systems to problems where the severity of the loading environment leads to important nonlinear effects. A recent review of these methods given by Mignolet et al. [1] discusses different approaches that have been developed and validated. These approaches lead to ROMs that are parametric and differ in the type of basis functions used to represent the structural response. Validation of these ROMs has been traditionally performed by comparing the ROM results to a "truth" model such as the response obtained from the FEM or data from experiments. In this case, the error introduced by the model reduction can be easily determined as the difference between the "truth" model and the ROM. Based on this error, the analyst can decide whether to enrich the ROM basis or not. In reality, however, loading conditions along an aircraft trajectory do not remain fixed, and comparing ROM results to the FEM is computationally expensive and defeats the purpose of the ROM. This highlights the importance in the development of error metrics that can provide information on the robustness and accuracy of the ROM solution. This in turn will lead to more robust ROMs and is an important step towards the development of adaptive nonlinear ROMs. In the context of a different field of study than the one from the present paper, Drohmann and Carlberg [2] introduced a technique for the modeling of errors introduced by ROMs of stationary problems. They use supervised machine learning methods to construct a mapping from ROM error indicators, such as the residual of the ROM solution, to a statistic distribution of the error. The approach is applied to the solution of a stationary heat transfer problem.
The present work is a preliminary effort in the development of error measures that can be used to shed light on the confidence of the ROM solution. A post-buckled beam (see schematic in Fig. 1) is used as an analogue of the geometrically curved or thermally buckled thin aircraft panels.
The outline of this paper is as follows. In Section 2, background on the ROM formulation is presented. Section 3 presents a description of the error metrics and Section 4 presents results of the application of the error metrics to the ROM solution of the nonlinear dynamic response of a thin beam. Section 5 summarizes the work presented in this paper.

ROM Formulation
The structural ROMs considered in this study are based on a representation of the nonlinear geometric response in terms of a set of basis vectors, where u(t) represents the vector of displacements of the finite element degrees of freedom, ψ (n) are specified, constant basis, and η n (t) are the time dependent generalized coordinates.
where M FE is the mass matrix, D FE is the damping matrix, F int is the restoring force vector, which includes linear and nonlinear components, and F ext is the vector with the external load at the nodes. As discussed in [1], the normal modes of vibration of the structure provide an excellent basis. The ROM equations are shown in Eq. (3) where ij M denotes the elements of the mass matrix,   K are the linear, quadratic, and cubic stiffness coefficients and F i are the modal forces. A detailed derivation of Eq. (2) is given by Kim et al. [3].
There are two important aspects in the construction of a ROM. First, the selection of the basis vectors, which following the discussions of Kim et al. [3], Hollkamp and Gordon [4], Muravyov and Rizzi [5] and Mignolet et al. [1], the core of the basis is formed by the linear modes of the structure, hereafter referred to as out-of-plane modes, which would appear significantly in the linear response. The other key aspect is the computation of the ROM parameters. The terms i F can conveniently be obtained from a finite element computation in which all degrees of freedom are blocked and in which the external loading is applied. Then, i F are obtained as the projections on the basis vector ) (m i  of the reaction forces followed by a change of sign. Finally, the modal mass and damping matrices are obtained as in the linear case: The identification of the stiffness coefficients ( ijlp K is also performed using the FEM. In this work, the identification procedure of Hollkamp and Gordon [4] is used, where the coefficients are estimated in a least squares approach using nonlinear static finite element solutions.

Parameter space study
A study of the forcing parameter and initial condition (IC) space was conducted for several ROMs of the post-buckled beam in order to identify regions where the proposed error metric could be evaluated. The forcing parameters of interest were the frequency and amplitude of a spatially uniform harmonic force. Wiebe and Spottswood conducted a similar parameter space evaluation [6] using the single degree-of-freedom Duffing equation because it allowed for high-resolution analysis of the parameter and IC space while still providing a qualitative comparison with a clamped, post-buckled beam experiment. High-dimension models, of even simple structures like the post-buckled clamped beam, make the parameter space study intractable due to the simulation expense. In the present study, nonlinear ROMs like those described previously, provided a mechanism to reduce the higherdimensional models into simpler ones thus allowing for such a detailed study. Figure denotes the parameter space response mapping for a 3-and 7-mode ROM. Even the 3-mode ROM of the postbuckled beam displays a sufficiently complex response for the selected parameters. Figure 3 shows the single-well and cross-well time histories for the highlighted points in the maroon and blue regions respectively in Figure . Figure was generated for a 160 x 160 grid of amplitude and frequency pairs. Each amplitude/frequency pair was simulated 100 times to gauge the sensitivity to initial conditions. Each grid point used different position and velocity "snap-shots" which were selected at random (uniformly in time) from prior-generated simulations of a chaotic response of the ROM. Specifically, a database of 1000 position and velocity "snap-shots" of the beam for each ROM was generated from simulation and then position and velocity pairs were randomly selected as initial conditions. The alternative approach of selecting spatially uncorrelated values as initial conditions for the beam was ruled out as this may have yielded very high internal forces and a poorly-conditioned tangent stiffness matrix. The 160 x 160 grid provided for a broad range of both loading amplitude and relevant frequencies as defined by the aforementioned experimental studies [6]. The response types are organized in the same manner as those used to characterize and explore the Duffing response in [6]: where the digits, in order, denote the existence (1) or nonexistence (0) of single-well non-chaotic, chaos, snap-through/cross well P n (n > 1), and cross-well P 1 response. Chaos was estimated in a similar fashion [6] by calculating the sign of the largest Lyapunov exponent using a renormalization procedure.   Sensitivity to initial conditions, a classic characteristic of nonlinear systems, is also present in this system in two different contexts. First, as shown in Fig. 4 (denoted by the white dot in the green region of Fig. 2), the system presents chaos, i.e., sensitivity to initial conditions for a time history. Second, the parameter space plots in Fig. 2 shows that beyond sensitivity to forcing parameters (the boundaries between response types appear fractal), many of the initial conditions yield coexisting responses. The former case makes it difficult to quantitatively determine the accuracy of a ROM. This is because even the 'truth' model will present quantitative differences in any two time series for changes in time step, initial conditions, or any other algorithmic change no matter how small. This requires further study and will not be discussed further. The issue of co-existing responses is particularly concerning when using ROMs, for example, the 3-mode ROM completely misses a coexisting snap-through response in the region highlighted by the white oval in Fig. 2 (a). It is hoped that in the future, the proposed error metrics may help uncover such potentially critical errors. In order to allow for a clear presentation of the method, however, the rest of the paper will remain in less troubled waters, i.e. the focus will be on the periodic responses in Fig. 3.

State-space angle and residual
The parameter space study presented in the previous subsection offers a qualitative approach to determining the robustness of a ROM to a wide array of loading parameters. The two approaches presented here add a quantitative component. Figure 4 shows a schematic of the method. The time stepping (using Newmark's method in this case, but any method could be used) is performed on the low-order ROM (hereafter called M 1 ) with N modes. The real state space is higher than 2-dimensional and cannot be plotted, however, the concept is equally as valid in any dimension. At each time step, the state space velocity vector of M 1 is given by where the term  is equal to the acceleration of the generalized coordinates of M 1 .
A higher-dimensional model, the 'truth model', is then used to approximate (in many cases an exact transformation is possible) the state of M 1 . This step may be done in a variety of ways, depending on what is used for the truth model. If a FEM is used, it is a trivial step to determine the physical/nodal coordinates of the FEM using the ROM coordinates and shape functions. Depending on the shape functions of the FEM, the transformation from M 1 to the FEM may be exact or approximate where ( ̅ , ̅ ) is the state of the lower-order ROM represented in the state space of the truth model. Because they are both ROMs, M 1 and M 2 may both be represented by Eq. (4). Overbars will be used to indicate that Eq. (3) or equivalently Eq. (4) is being applied to the truth model, and no overbars will be used when it is the lower-order ROM. Thus, the equivalent state in Eq. (5) can then be used to determine the equivalent state space acceleration given by and the overall equivalent state space velocity vector given by The new state-space velocity vector z  represents the direction that the M 2 system will take if the displacement and velocity of the M 1 system are imposed. Therefore, the angle between the M 2 statespace velocity (the truth model) and z  is a measure of how close the solution of the lower-order ROM (M 1 ) is to the response predicted by the higher-order ROM (M 2 ). At a given time instance, the solution of the ROM equations is sought to minimize the force residual: Computing the force residual using the higher-fidelity ROM and the state of the lower-fidelity ROM could also serve as an indication of the accuracy of the lower-fidelity ROM solution. Therefore, the other error metric that will be used in this work is the norm of the force residual computed with the parameters of the higher-fidelity ROM and the state of the lower-fidelity one.

Application of Error Metrics
The error metric techniques described above were assessed on an isotropic flat beam subjected to a uniform temperature and a load proportional to its first mode shape. The material properties of the beam are given in Table 1. The beam was modeled using an in-house finite element code. The geometry was discretized using 40 beam elements and the two ends of the beam were clamped. The temperature of the beam was set equal to 11.6°C above the reference temperature; this is equal to approximately three times the buckling temperature of the beam.

Results
The beam response resulting from the parameters selected in the study from Section 3.1 will be used in this section. The 7-mode ROM will be considered the "truth" model and the error metrics discussed in the previous section will be used to assess the error in the response obtained from a 3mode and a 5-mode ROM. The response predicted by all ROMs to the two sets of parameters selected before are qualitatively the same. The periodic single-well response corresponds to a forcing frequency of 70Hz and amplitude of 2.1g's while the periodic cross well behavior was obtained with a forcing frequency of 138Hz and forcing amplitude equal to 18g's (see Fig. 3 for time series). Shown in Fig. 6 are the results for the single well response at the center of the beam. Figure 6 shows projected state-space plots (for the displacement and velocity at the (a) ¼ point and (b) midspan) of one period of the stationary response of the beam. Note that the initial rise and vibration amplitudes at the quarter point are smaller, however, in both it is easy to see a very large discrepancy between the 3-mode ROM and the 7-mode ROM, while the difference between the 5-mode ROM and 7-mode ROM is smaller. This observation can be made as well by looking at the norm of the force residual (Fig. 7) which shows a large reduction between the 3-mode and the 5-mode ROMs. The second error metric, that is the angle between the state space velocity vectors of the 'truth' model and the approximate models (θ 3 for the 3-mode ROM and θ 5 for the 5-mode ROM) has also been shown at a few selected points on the orbit. For the points shown, the 5-mode model shows a smaller angle and hence a reduced error in the direction of the state space velocity vector. Note that when components of the acceleration vector reverse sign, it is common to see, for a very brief period of time, a very large angle error. That is because the angle will show a ~180 degree error in any component for which the acceleration of the approximate model does not reverse sign at the exact same time as the truth model. Thus some filtering is required to eliminate very brief spikes in the angle error. The points shown in Fig. 6 are near the region of maximum acceleration in order to ensure this does not occur. A similar series of plots is shown for the cross-well behavior Fig. 8 and Fig. 9. Similar to the single-well response, the 5-mode ROM matches the response predicted by the 7-mode ROM better than the 3-mode model, however, the improvement is less drastic than it is for the single-well response. It is interesting, however, that the angle error at regions of high acceleration appears to be similar or slightly better for the 3-mode model. The magnitudes of the error, however, are 1-2 orders of magnitude lower than they are for the single-well response, suggesting both models are doing significantly better than they did for the single-well response. This would indicate that the cross-well response engages fewer modes, and therefore has a lower spatial complexity than the single-well response. For the case of the cross-well response it is known that the first asymmetric mode is necessary to capture correctly the snap-through events. For illustration purposes, Figs. 10 and 11 show the values of the error metrics for a 2-mode ROM built using the first two symmetric modes (i.e., no asymmetric mode is present). The discrepancy with the 7-mode ROM is quite large and the difference in the error metrics between these results and the ones shown previously for the 3-mode and 5-mode ROM is more significant.

Summary
The work presented in this paper is a preliminary investigation on the development of error metrics capable of providing an indication on the level of the approximation error in a ROM solution. This error is the result of representing the displacement field of the response in terms of a truncated basis. The first error metric explored was the angle between the state-space velocity of a low-fidelity ROM and an equivalent state-space velocity obtained by imposing the state of the low-fidelity ROM on the equations of motion of a high-fidelity ROM. The second metric was the norm of the force residual obtained by imposing the state of the low-fidelity ROM on the equations of motion of the high-fidelity ROM. The first method yields information about the error in the magnitude in the direction of the state space velocity vector of the ROM, and the second about the error in the magnitude. When combined these two error metrics indicate the effect of the 'artificial stiffening' of the ROM, something which is inevitable when truncating the DOFs of a structure. Although the error metrics rely of a 'truth model' (in this paper a higher-order ROM), they have the advantage that no time stepping of the higher-cost model is required.
The results of a parameter and initial condition space study provided a qualitative measure of the accuracy of a 3-mode ROM of a clamped-clamped buckled beam. The study classified the response of the beam as single-well non-chaotic, chaos, snap-through/cross well P n (n > 1) and cross-well P 1 . Two sets of parameters were selected for which the same type of response was obtained with a 3-mode and a 7-mode ROM, one corresponding to a periodic single-well response and the other set to a periodic cross-well response. Even though qualitatively the types of the response predicted by the two ROMs were the same, the norm of the force residual metric indicated a larger error in the prediction of the 3mode ROM than for a 5-mode ROM. The state space velocity angle error metric also showed that the 5-mode ROM reduced the error for the single-well response, but did very little, or even slightly increased the error for the cross-well response. It is expected that response somewhere in between, linger near unstable equilibria (i.e. orbits that almost enough, or barely have enough energy to snapthrough) will present the most spatially-complex responses that may be useful in exercising the methods. Unfortunately, these types of responses frequently exhibit chaos, for which quantitative comparisons are difficult and often meaningless. More research is required to investigate this case.
It should be noted, that while the ROMs in this study gave good results with relatively few modes, more complex structures such as plates and shells may require far more modes making the benefit of a ROM much larger. It is hoped that with further study proposed error metrics will allow for similar studies to be done using the ROM alone without reliance on a purpose-defeating computationallyexpensive FEM running simultaneously.