Semi-implicit Integration and Data-Driven Model Order Reduction in Structural Dynamics with Hysteresis

A key idea in numerical integration of ordinary differential equations (ODEs) is outlined here for a general reader.

We consider ODE systems written in first order form, $\dot{\operatorname{\boldsymbol{y}}}=\boldsymbol{f}(\operatorname{\boldsymbol{y}},t)$, where $\operatorname{\boldsymbol{y}}$ is a state vector and $t$ is time. In single-step methods, which we consider in this paper, we march forward in time using some algorithm equivalent to

The specific form of $\boldsymbol{H}$ above is derived from the form of $\boldsymbol{f}(\operatorname{\boldsymbol{y}},t)$ and depends on the integration algorithm. The actual evaluation of $\boldsymbol{H}$ may involve one or multiple stages, but that is irrelevant: the method is single-step in time. For smooth systems, $\boldsymbol{H}$ is guided by the first few terms in Taylor expansions of various quantities about points of interest in the current time step. In such cases, as $h\rightarrow 0$, the error in the computed solution goes to zero like $h^{m}$ for some $m>0$. If $m>1$, the convergence is called superlinear. If $m=2$, the convergence is called quadratic. Values of $m>2$ are easily possible for smooth systems with moderate numbers of variables: see, e.g., the well known Runge-Kutta methods [3]. Unfortunately, for large structural systems, for the asymptotic $h^{m}$ rate to hold, $h$ may need to be impractically small. For example, the highest natural frequency of a refined finite element (FE) model for a structure may be, say, $10^{6}$ Hz. This structure may be forced at, say, 10 Hz. A numerical integration algorithm that requires time steps that resolve the highest frequency in the structure, i.e., time steps much smaller than $10^{-6}$ seconds in this example, is impractical. A high order of convergence that requires time steps smaller than $10^{-6}$ is of little use. We seek accurate results with time steps much smaller than the forcing period, but much larger than the time period of the highest natural frequency of the structure, e.g., $10^{-2}$ or $10^{-3}$ seconds. To develop such practically useful algorithms, we must consider the stability of the numerical solution for larger values of $h$, i.e., $10^{-6}<h<10^{-2}$, say.

Now consider the nature of $\boldsymbol{H}$. If $\boldsymbol{H}$ does not depend explicitly on $\operatorname{\boldsymbol{y}}(t+h)$, the algorithm is explicit. Otherwise it is implicit. For nonlinear systems, $\boldsymbol{H}$ is a nonlinear function of its $\operatorname{\boldsymbol{y}}$-arguments. Then each implicit integration step requires iterative solution for $\operatorname{\boldsymbol{y}}(t+h)$. For linear dynamics, with $\boldsymbol{H}$ linear in its $\operatorname{\boldsymbol{y}}$-arguments, the $\operatorname{\boldsymbol{y}}(t+h)$ can be moved over to the left hand side and integration proceeds without iteration, although usually with matrix inversion. The algorithm is still called implicit in such cases: implicitness and iteration are not the same. An algorithm can be neither fully implicit nor fully explicit, and we will present one such algorithm in detail.

We present a semi-implicit approach that can be used for high dimensional finite element models with distributed hysteresis. For simplicity, we adopt the popular Bouc-Wen model [4, 5, 6] as our damping mechanism. After presenting and validating our numerical integration method, we present a way to obtain useful lower order models for the structure, starting from a refined FE model. A key issue is that the refined or full FE model computes hysteretic variables at a large number of Gauss points in the structure, and a smaller subset needs to be used for the model order reduction to be practical.

Thus, our contribution is twofold. First, we present a simple semi-implicit algorithm for a structural FE model with distributed hysteresis and demonstrate its convergence and utility. Second, we use this algorithm to compute some responses of the structure and use those responses to construct accurate lower order models with reduced numbers of both vibration modes and hysteretic variables. These lower order models can be used later for quick simulations under similar initial conditions or loading.

We begin with a review of some numerical integration methods that are available in popular software or research papers.

Structural systems with Bouc-Wen hysteresis continue to be studied in research papers using algorithms that are not as efficient as the one we will present below. For example, as recently as 2019, Vaiana et al. [14] considered a lumped-parameter model and used an explicit time integration method from Chang’s family [15]. That method requires tiny time steps: the authors used one hundredth of the time period of the highest structural mode. Such small time steps are impractical for refined FE models. Our algorithm allows much larger time steps. Thus, in the area of FE models with distributed hysteresis, our paper makes a useful contribution.

Next, we acknowledge that the commercial finite element software Abaqus [7] allows users to specify complex material responses and also to choose between explicit and implicit numerical integration options. For many nonlinear and nonsmooth dynamic problems, explicit integration needs to be used in Abaqus. As outlined above, implicit or semi-implicit algorithms can be useful for somewhat simpler material modeling and in-house FE codes.

Considering general purpose software for numerical work, many analysts dealing with hysteresis may begin with Matlab [10]. Matlab’s built in function ode15s is designed for stiff systems but not specifically for systems with hysteresis. We have found that ode15s can handle ODEs from FE models with a modest number of elements and with hysteresis, but it struggles with higher numbers of elements because its adaptive time steps become too small. In this paper we will use ode15s to obtain numerically accurate solutions for systems of moderate size. Results obtained from our own algorithm will then be validated against ode15s results. Subsequently, for larger systems, our algorithm will continue to work although ode15s freezes and/or crashes.

For those programming their own integration routines in structural dynamics, the well known Newmark method [11] from 1959 remains popular (see, e.g., [9, 8]), although it cannot guarantee stability and may require tiny time steps as noted in, e.g., [12]. In that paper [12] of 1977, Hilber, Hughes and Taylor modified the Newmark method and showed unconditional stability for linear structural problems. However, their method (known as HHT-$\alpha$) can show spurious oscillations of higher modes even without hysteretic damping (see, e.g., Fig. 7 of the paper by Piché [31], which we will take up below).

An appreciation of the issues faced for full three-dimensional (3D) simulation with hysteretic damping can be gained, e.g., from the work of Triantafyllou and Koumousis [13]. Their formulation is actually based on plasticity in 3D; they compare their solutions with those from Abaqus; and they include dynamics through the Newmark algorithm. Their algorithm is rather advanced for a typical analyst to implement quickly. Note that our present application is easier than theirs because we have only one-dimensional Bouc-Wen hysteresis. Additionally, our primary application here is in model order reduction. And so we develop for our current use a numerical integration approach that is simpler than that in [13].

Finally, readers interested in hybrid simulations (with an experimental structure in the loop) may see, e.g., Mosqueda and Ahmadizadeh [16, 17] who used a modified version of the HHT-$\alpha$ for calculating the seismic response of a multi-degree-of-freedom structure. Purely simulation-based studies such as the present one can hopefully provide theoretical inputs into planning for such hybrid simulations in future work.

We close this section with a brief discussion of model oder reduction. While there are very many papers on the topic, we make special note of the condensation approaches of Guyan [18] and Irons [19], wherein finding normal modes of the full system was avoided. With growth in computational power, system eigenvectors become more easily available, and direct modal projection has become common. For a recent discussion of this and related methods, see [20]. We will use straightforward modal projection for our displacement degrees of freedom. However, we will also have a large number of internal hysteretic state variables. We will present a sequential approach to selecting a subset of these hysteretic variables, using an algorithm which is related closely to a method called “QR with column pivoting” in [21]. Minor differences from [21] here are that we work with rows and not columns, and our data matrix is not rank deficient: we separately impose a termination criterion.

Note that modal reduction, without the added complication of hysteretic damping, is well known. For example, Stringer et al.[22, 23] successfully applied modal reduction to rotor systems, including ones with gyroscopic effects. Other recent examples in dynamics and vibrations can be seen in [24, 25]. Proper Orthogonal Decomposition (POD) [26] based reduced order models are often used in fluid mechanics: a representative sample may be found in [27, 28, 29, 30].

We will begin by adopting an existing implicit scheme for linear structural dynamics without hysteresis that is both simple and significantly superior to both the Newmark and HHT-$\alpha$ methods. Our adopted scheme, due to Piché [31], is an L-stable method from the Rosenbrock family (see [32] and references therein) which is stable, implicit without iteration for linear structures, and second order accurate. For smooth nonlinear systems of the form

Piché suggests a one-time linearization of $\boldsymbol{f}$ at the start of each time step. We will not use that linearization step because hysteresis is nonsmooth.

Here we extend Piché’s formulation to include hysteretic damping in the otherwise-linear structural dynamics. In our method the stiffness and inertia terms are integrated implicitly (without iteration, being linear) and the nonsmooth hysteresis is monitored at Gauss points [33] and treated with explicit time-integration. Thus, overall, our method is semi-implicit. Our proposed approach, when applied to a refined FE model of an Euler-Bernoulli beam with distributed Bouc-Wen hysteretic damping, is easy to implement and can be dramatically faster than Matlab’s ode15s. In fact, it continues to work well beyond refinement levels where ode15s stops working.

We emphasize that ode15s is not really designed for such systems. We use it here because it has adaptive step sizing and error estimation: when it does work, it can be highly accurate. For this reason, to examine the performance of our algorithm for modest numbers of elements, we will use results obtained from ode15s with a tight error tolerance. With higher number of elements, ode15s fails to run because the time steps required become too small.

Having shown the utility of our proposed semi-implicit integration algorithm, we will turn to our second contribution of this paper, namely model order reduction. Such Reduced Order Models (ROMs) can save much computational effort. In the physical system, if only a few lower modes are excited, the displacement vector can be approximated as a time-varying linear combination of those modes only. A remaining issue here is that the number of Gauss points used for the hysteresis can be reduced too, and that is where we offer an additional contribution.

In recent work related to ours, for an Euler-Bernoulli beam with hysteretic dissipation, Maity et al. [33] used the first few undamped analytically obtained modes to approximate the solution and performed the virtual work integration using a few Gauss points chosen over the full domain. However, they used a different hysteresis model motivated by distributed microcracks. Furthermore, in our finite element model, virtual work integrations are performed over individual elements and not the whole domain. The number of Gauss points in the finite element model increases with mesh refinement. When we project the response onto a few modes, we can explicitly reduce the number of Gauss points retained as well, and a practical method of doing so is one of the contributions of this paper.

In what follows, we present the finite element model in section 2, outline the numerical integration algorithm in section 3, develop the approach for model order reduction in section 4, and present our results in section 5.

The elemental stiffness and inertia matrices are routine and given in many textbooks: each beam element has two nodes, and each node has one translation and one rotation. The hysteresis variable $z$ is a scalar quantity distributed in space. However, it enters the discretized equations through integrals, and so $z$ values only need to be known at Gauss points within each element. The evolution of $z$ at the Gauss points is computed using Eq. (3), and is thus driven by displacement variables. Some further details are given in appendix A.

After assembly we obtain a system of equations of the form

where $\mathbf{B}\in\mathbb{R}^{n_{\textrm{e}}n_{\textrm{g}}\operatorname{\times}N}$ and the different symbols mean the following:

1. 1.

   $\boldsymbol{q}$ is a column vector of nodal displacements and rotations for $n_{\textrm{e}}$ beam elements, with $N=2n_{\textrm{e}}$ for a cantilever beam,

2. 2.

   $\mathbf{M}$ and $\mathbf{K}$ are mass and stiffness matrices of size $2\,n_{\textrm{e}}\times 2\,n_{\textrm{e}}$,

3. 3.

   $\boldsymbol{z}$ is a column vector of length $n_{\textrm{g}}n_{\textrm{e}}$ from $n_{\textrm{g}}$ Gauss points per element,

4. 4.

   $(\cdot\circ\cdot)$ and $(\cdot)^{\circ\,(\cdot)}$ denote elementwise multiplication and exponentiation,

5. 5.

   $\mathbf{A}$ is a matrix of weights used to convert $\boldsymbol{z}$ values into virtual work integrals,

6. 6.

   $\boldsymbol{\chi}$ is a column vector of curvatures at the Gauss points,

7. 7.

   $\mathbf{B}$ maps nodal displacements and rotations $\boldsymbol{q}$ to curvatures $\boldsymbol{\chi}$ at the Gauss points, and

8. 8.

   $\boldsymbol{f}_{0}(t)$ incorporates applied forces.

In Eq. (7a) above we can include viscous damping by adding $\mathbf{C}\dot{\boldsymbol{q}}$, for some suitable $\mathbf{C}$, on the left hand side. As this paper focuses on hysteretic damping, we have taken ${\bf C}={\bf 0}.$

Piché’s algorithm uses a numerical parameter $1-\frac{1}{\sqrt{2}}$ which is written as $\gamma$ for compact presentation. The symbol should not be confused for Euler’s constant. We now proceed as follows. This subsection presents implicit integration for the structural part alone: the hysteretic variable vector $\boldsymbol{z}$ is not integrated in an implicit step: this compromise simplifies the algorithm greatly.

1. 1.

   Define

   $$\boldsymbol{z}(t_{0})=\boldsymbol{z}_{0},\quad\boldsymbol{F}_{0}=-\mathbf{A}\boldsymbol{z}_{0}+\boldsymbol{f}_{0}(t_{0}),$$

   $$\operatorname{\boldsymbol{y}}_{0}=\boldsymbol{q}(t_{0}),\quad\operatorname{\boldsymbol{v}}_{0}=\dot{\boldsymbol{q}}(t_{0}),\quad\dot{\boldsymbol{\chi}}_{0}=\mathbf{B}\operatorname{\boldsymbol{v}}_{0},$$

   $$\dot{\boldsymbol{z}}_{0}=\left(\bar{A}-\alpha\,\textrm{sign}(\dot{\boldsymbol{\chi}}_{0}\circ\boldsymbol{z}_{0})\circ\left\lvert{\boldsymbol{z}_{0}}\right\rvert^{\circ n_{\textrm{h}}}-\beta\left\lvert{\boldsymbol{z}_{0}}\right\rvert^{\circ n_{\textrm{h}}}\right)\circ\dot{\boldsymbol{\chi}}_{0},$$

   $$\dot{\boldsymbol{F}}_{0}=-\mathbf{A}\dot{\boldsymbol{z}}_{0}+\frac{\operatorname{\textrm{d}\!}}{\operatorname{\textrm{d}\!}{t}}\boldsymbol{f}_{0}(t)\bigg{\rvert}_{t=t_{0}},$$

   $$\boldsymbol{r}_{0}=\mathbf{K}\operatorname{\boldsymbol{y}}_{0}+\mathbf{C}\operatorname{\boldsymbol{v}}_{0}$$

   (The $\mathbf{C}\operatorname{\boldsymbol{v}}_{0}$ is dropped if linear viscous damping is no included.)

2. 2.

   Define

   $$\tilde{\mathbf{M}}=\mathbf{M}+\gamma h\mathbf{C}+(\gamma h)^{2}\,\mathbf{K}.$$

3. 3.

   Define (first stage)

   $$\tilde{\boldsymbol{e}}=h\,{\tilde{\mathbf{M}}}^{-1}\left(\boldsymbol{F}_{0}-\boldsymbol{r}_{0}+h\gamma(\dot{\boldsymbol{F}}_{0}-\mathbf{K}\boldsymbol{v}_{0})\right).$$

   $$\tilde{\boldsymbol{d}}=h(\boldsymbol{v}_{0}+\gamma\tilde{\boldsymbol{e}}).$$

4. 4.

   Define (second stage)

   $$\operatorname{\boldsymbol{F}}_{\operatorname{\frac{1}{2}}}=-\mathbf{A}\left(\boldsymbol{z}_{0}+\frac{h}{2}\dot{\boldsymbol{z}}_{0}\right)+\boldsymbol{f}_{0}\left(t_{0}+\operatorname{\frac{1}{2}}h\right),$$

   $$\boldsymbol{r}_{\operatorname{\frac{1}{2}}}=\mathbf{K}\left(\operatorname{\boldsymbol{y}}_{0}+\operatorname{\frac{1}{2}}\tilde{\boldsymbol{d}}\right)+\mathbf{C}\left(\operatorname{\boldsymbol{v}}_{0}+\operatorname{\frac{1}{2}}\tilde{\boldsymbol{e}}\right),$$

   $$\boldsymbol{e}=h\tilde{\mathbf{M}}^{-1}\left(\operatorname{\boldsymbol{F}}_{\operatorname{\frac{1}{2}}}-\boldsymbol{r}_{\operatorname{\frac{1}{2}}}+(h\gamma)\left(2\gamma-\operatorname{\frac{1}{2}}\right)\mathbf{K}\tilde{\boldsymbol{e}}+\gamma\mathbf{C}\tilde{\boldsymbol{e}}\right),$$

   and

   $$\boldsymbol{d}=h\left(\boldsymbol{v}_{0}+\left(\operatorname{\frac{1}{2}}-\gamma\right)\tilde{\boldsymbol{e}}+\gamma\boldsymbol{e}\right).$$

5. 5.

   Finally

   $$\boldsymbol{y}(t_{0}+h)=\operatorname{\boldsymbol{y}}_{0}+\boldsymbol{d},\qquad\operatorname{\boldsymbol{v}}(t_{0}+h)=\operatorname{\boldsymbol{v}}_{0}+\boldsymbol{e}.$$

6. 6.

   Define

   $$\dot{\boldsymbol{\chi}}_{1}=\dot{\boldsymbol{\chi}}(t_{0}+h)=\mathbf{B}(\operatorname{\boldsymbol{v}}_{0}+\boldsymbol{e})$$

Note that the assignment of $\dot{\boldsymbol{\chi}}_{1}$ to $\dot{\boldsymbol{\chi}}(t_{0}+h)$ above is tentative at this stage; if there is a sign change, we will make a correction, as discussed in the next subsection.

Due to the nonanalyticity of hysteresis models, we integrate the hysteresis part using an explicit step. There are slope changes in the hysteretic response whenever $\dot{\chi}$ at any Gauss point crosses zero in time. We will accommodate the sign change of $\dot{\chi}$ within an explicit time step in a second sub-step, after first taking a time step assuming that there is no sign change. Since each $\dot{z}_{i}$ depends on $\dot{\chi}_{i}$ only, accounting for zero crossings can be done individually for individual Gauss points and after such a preliminary step.

We first check elementwise and identify the entries of $\dot{\boldsymbol{\chi}}_{0}$ and $\dot{\boldsymbol{\chi}}_{1}$ that have same sign, and construct sub-vectors $\dot{\boldsymbol{\chi}}_{0\textrm{u}}$ and $\dot{\boldsymbol{\chi}}_{1\textrm{u}}$ out of those entries. The corresponding elements from the vector $\boldsymbol{z}_{0}$ are used to construct a sub-vector $\boldsymbol{z}_{0\textrm{u}}$. Here the subscript “${\rm u}$” stands for “unchanged”. These elements can be used in a simple predictor-corrector step for improved accuracy as follows.

and

Next, we address the remaining entries of $\dot{\boldsymbol{\chi}}_{0}$ and $\dot{\boldsymbol{\chi}}_{1}$, those that have crossed zero and flipped sign. We construct sub-vectors $\dot{\boldsymbol{\chi}}_{0\textrm{f}}$ and $\dot{\boldsymbol{\chi}}_{1\textrm{f}}$ out of them, where the “${\rm f}$” subscript stands for “flipped”. The corresponding elements from the vector $\boldsymbol{z}_{0}$ are used to construct sub-vector $\boldsymbol{z}_{0\textrm{f}}$. Using linear interpolation to approximate the zero-crossing instants within the time step, we use

We define

and

which is consistent with Eq. (8) because $\dot{\boldsymbol{\chi}}=\boldsymbol{0}$ at the end of the sub-step. Next we complete the step using

and

where $\boldsymbol{h}$ is a vector of the same dimensions as $\boldsymbol{h}_{0}$ and with all elements equal to $h$. Finally, having the incremented values $\boldsymbol{z}_{1}$ at all locations, those with signs unchanged and those with signs flipped, we use

We clarify that the explicit integration of the hysteretic variables, as outlined in this subsection, is a compromise adopted to avoid difficulties due to the nonanalyticity of the hysteresis model. Continuing to integrate the inertia and stiffness parts explicitly will still allow us to use usefully large steps, as will be seen in section 5.

Having the numerical integration algorithm in place, however, we must now turn to the second contribution of this paper, namely model order reduction.

The undamped normal modes and natural frequencies are given by the eigenvalue problem

for $N$ eigenvector-eigenvalue pairs $(\boldsymbol{v}_{i},\omega_{i})$. In a finely meshed FE model, $N$ is large. In such cases, we may compute only the first several such pairs using standard built-in iterative routines in software packages.

The $N$ dimensional displacement vector $\boldsymbol{q}$ is now approximated as a linear combination of $r\ll N$ modes. To this end, we construct a matrix $\mathbf{R}$ with the first $r$ eigenvectors as columns so that

where $\mathbf{R}=[\operatorname{\boldsymbol{v}}_{1}\,\;\operatorname{\boldsymbol{v}}_{2}\,\dots\,\operatorname{\boldsymbol{v}}_{r}]$, and where the elements of $\boldsymbol{\xi}(t)$ are called modal coordinates. Substituting Eq. (10) in Eq. (7a) and pre-multiplying with the transposed matrix $\mathbf{R}^{\top}$, we obtain

obtaining equations of the form

In the above, due to orthogonality of the normal modes, the matrices $\tilde{\mathbf{M}}$ and $\tilde{\mathbf{K}}$ are diagonal. Also, $\boldsymbol{\chi}_{\textrm{a}}$ is an approximation of the original $\boldsymbol{\chi}$ because we have projected onto a few vibration modes, but it still has the same number of elements as $\boldsymbol{\chi}$ and requires numerical integration of the same number of nonsmooth hysteresis equations. A reduction in the number of hysteretic variables is necessary.

The system still has a large number ($n_{\textrm{g}}n_{\textrm{e}}$) of hysteretic damping variables. These are arranged in the elements of $\boldsymbol{z}$ in Eq. (7a). Selecting a smaller set of basis vectors and projecting the dynamics of the evolving $\boldsymbol{z}$ onto those has analytical difficulties. Here we adopt a data-driven approach to select a submatrix of $\boldsymbol{z}$, say

The indices $j_{1},\,j_{2},\,\dots j_{m}$ must now be chosen form the set $\{1,2,\dots,n_{\rm e}n_{\rm g}\}$, along with a matrix ${\bf P}$, such that

Working with that reduced set of hysteretic variables, we will be able use a reduced set of driving local curvatures

a submatrix of $\dot{\boldsymbol{\chi}}_{\textrm{a}}$. In other words, we will work with a subset of the original large number of Gauss points used to compute virtual work integrals for the hysteretic dissipation. In our proposed data-driven approach, selection of the indices $j_{1},\,j_{2},\,\dots j_{m}$ is the only significant issue. The matrix ${\bf P}$ can then be found by a simple matrix calculation.

However, for implementing our approach, we must first generate a sufficiently large amount of data using numerical integration of the full equations with representative initial conditions and/or forcing.

These initial conditions are randomly generated as follows. The hysteretic states are uniformly distributed random numbers in the interval (0, 0.1). For the displacement and rotation degrees of freedom, the initial conditions have nonzero values along only the first three modes. These initial conditions, in all cases, are scaled to make the tip displacement to be 2 cm. The three modal coordinate values are three randomly chosen, positive, independent and identically distributed numbers, subsequently rearranged so that the first mode had the largest displacement and the third mode had the smallest displacement. The initial values of the velocity degrees of freedom are taken to be zero. In applications, users who have other preferred criteria for choosing representative initial conditions can freely use them with no consequences for the rest of the algorithm.

Having selected initial conditions, we must simulate the system by integrating the equations forward in time. To that end, we use the semi-implicit integration method presented in section (3).

The data generation process is relatively slow, but once it is done and the reduced order model is constructed, we can use it for many subsequent quicker calculations. Beyond its academic interest, our approach is practically useful in situations where the reduced order model will be used repeatedly after it is developed.

It is perhaps not strictly necessary to place both of our contributions within one paper, but in our implementation the second part relies heavily on the first part; moreover, the second part is a relatively small and simple application of an efficient method to generate the needed data.

We consider a cantilever beam with Young’s modulus $E=200$ GPa and density $\rho=7850\,\textrm{kg}\,{\rm m}^{-3}$; of length $1\,\textrm{m}$ and square cross section of $2\,{\rm cm}\times 2\,{\rm cm}$. This yields a flexural rigidity $EI=2666.7\,\textrm{N}\textrm{m}^{2}$ and mass per unit length $\bar{m}=3.14\,\textrm{kg}\,\textrm{m}^{-1}$.

For many metals, 0.2% strain is near the border of elastic behavior. For the beam parameters above, if the beam is statically deflected under a single transverse tip load, then 0.2% strain at the fixed end of the beam corresponds to a tip deflection of 6 cm, which is taken as a reasonable upper bound for the vibration amplitudes to be considered in our simulations below.

Next, we consider the hysteretic damping itself. The index $n_{\rm h}$ in the Bouc-Wen model is primarily taken to be 0.5 for reasons explained a little later. A larger value is considered in subsection 5.4, and only in subsection 5.4, to clarify an issue in convergence. The parameters $\alpha$ and $\beta$ are somewhat arbitrary; we have found that the values $\alpha=0.8$ and $\beta=0.5$ yield hysteresis loops of reasonable shape. It remains to choose the Bouc-Wen parameter $\bar{A}$.

It is known that for small amplitude oscillations, the Bouc-Wen dissipation follows a power law. An estimate for the upper limit of forcing amplitude in the Bouc-Wen model, below which the power law should hold, is given in [35] as

Choosing $\chi_{\rm max}$ to correspond to the abovementioned 0.2% strain at the fixed end, in turn corresponding to a static free-end deflection of 6 cm for the above beam, we find from Eq. (18) that

Only one physical model parameter remains to be chosen, namely $\gamma_{\textrm{h}}$. To choose this parameter value, we note that the amplitude decay in free vibrations with Bouc-Wen hysteretic dissipation is not exponential. Thus, the idea of percentage of damping is only approximately applicable. Upon looking at the reduction in amplitude after a certain number of oscillations (say $M$), an equivalent “percentage damping” can be approximated using the formula

Metal structures often have 1-2 % damping [34]. We will choose values of the hysteretic dissipation $\gamma_{\textrm{h}}$ (recall Eq. (2)) to obtain $0.01\leq\zeta_{\textrm{equiv}}\leq 0.02$. Numerical trial and error show that

is a suitable value (see figure 3).

Before we examine the performance of our proposed semi-implicit integration method, we will verify that the FE model indeed displays power law damping at small amplitudes, as predicted by theoretical analyses of the Bouc-Wen model.

The small amplitude dissipation per cycle of the Bouc-Wen model is proportional to amplitude to the power $n_{\textrm{h}}+2$ [35]. Since the energy in the oscillation is proportional to amplitude squared, we should eventually observe small amplitude oscillations in the most weakly damped mode with amplitude $A$ obeying

whence for $n_{\textrm{h}}>0$

For the Bouc-Wen model, letting $n_{\textrm{h}}\rightarrow 0$ produces hysteresis loops of parallelogram-like shape, and so we prefer somewhat larger values of $n_{\textrm{h}}$; however, to have a significant decay rate even at small amplitudes, we prefer somewhat smaller values of $n_{\textrm{h}}$. As a tradeoff, we have chosen $n_{\textrm{h}}=\operatorname{\frac{1}{2}}$. We expect an eventual decay of vibration amplitudes like $1/t^{2}$.

For our beam model, with 10 elements and 3 Gauss points per element, and with our semi-implicit numerical integration algorithm²²2Matlab’s ode15s struggles with such long simulations on a small desktop computer; our algorithm works quickly., the computed tip displacement asymptotically displays the expected power law decay rate, as seen on linear axes in Fig. (4a) and more clearly in the logarithmic plot of Fig. (4b). The frequency of the residual oscillation is close to the first undamped natural frequency of the beam. Higher modes are not seen in the long-term power law decay regime.

Having checked the asymptotic power law decay rate, we next ensure that the semi-implicit algorithm does not produce spurious oscillations in higher modes within the computed solution. This issue is important because, with increasing mesh refinement in the FE model, very high frequencies are unavoidable. While those high frequency modes exist in principle, their response should be small if any external excitation is at low frequencies and initial conditions involve only the lower modes. To examine this aspect, we denote the time period of the highest mode present in the structural model by $T_{\textrm{min}}$. As the number of elements increases, $T_{\textrm{min}}$ decreases, as indicated in Fig. (5a).

We now use our semi-implicit method to simulate an FE model with 100 elements and with nonzero initial conditions along only the first 3 modes. In the computed solution, we expect to see frequency content corresponding to the first three modes only. Figure 5b indeed shows only three peaks. Spurious responses of the higher modes are not excited. Note that this simulation was done with a time step of $h=10^{-4}$, which is more than 275 times larger than the time period of the highest mode for $n_{\rm e}=100$, which is $T_{\min}=3.6\operatorname{\times}10^{-7}\,\textrm{sec}$ (see Fig. (5a)). It is the implicit part of the code that allows stable numerical integration with such a large step size.

Having checked that spurious oscillations in very high modes are not excited, it remains to check that accurate results are obtained with time steps $h$ which are sufficiently small compared to the highest mode of interest. To this end, the convergence of the solution with decreasing $h$ is depicted in Fig. (6). For the beam modeled using 10 elements, the first five natural frequencies are

Numerical simulation results using our semi-implicit algorithm are shown in Fig. (6) for $h=10^{-3}$, $10^{-4}$ and $10^{-5}$. The overall solution is plotted on the left, and a small portion is shown enlarged on the right. Only 10 elements were used for this simulation to allow use of Matlab’s ode15s, which has adaptive step sizing and allows error tolerance to be specified (here we used $10^{-10}$). It is seen in the right subplot that although all three solutions from the semi-implicit method are stable, the one with $h=10^{-3}$ does not do very well in terms of accuracy; the one with $h=10^{-4}$ is reasonably close and may be useful for practical purposes; and the one with $h=10^{-5}$ is indistinguishable at plotting accuracy from the ode15s solution. The match between our semi-implicit integration with $h=10^{-5}$ and Matlab’s ode15s with error tolerances set to $10^{-10}$ indicates that both these solutions are highly accurate.

Using simulations with relatively few elements and over relatively short times, we can compare the results obtained from our semi-implicit method (SIM) with ode15s, and examine convergence as $h$ is made smaller.

Figure 6 shows the tip response for an FE model with 10 elements for different time step sizes and compares them with the ode15s solution. Here we will use two different error measures.

Everywhere in this paper except for this single subsection, we have used $n_{\rm h}=0.5$. In this subsection only, we consider $n_{\rm h}=1.5$. Due to the greater smoothness of the hysteresis loop near the point B of Fig. (8), we expect better convergence for this higher value of $n_{\rm h}$.

Some parameter choices must be made again. Using the yield criteria used in section 5.1, for $n_{\rm h}=1.5$, we find we now require

Subsequently, we find an approximate equivalent damping ratio $\zeta_{\rm equiv}=0.015$ for $\gamma_{\rm h}=0.3$. It is interesting to note that with the change in $n_{\rm h}$ and for the physical behavior regime of interest, $\bar{A}$ and $\gamma_{\rm h}$ have individually changed a lot but their product has varied only slightly (195 Nm for about 1.6% damping in the $n_{\rm h}=0.5$ case, and 183 Nm for about 1.5% damping in the $n_{\rm h}=1.5$ case).

The decay of tip response amplitude (Fig. (10)) for these parameters (with $n_{\rm h}=1.5$) looks similar to the case studied in the rest of this paper ($n_{\rm h}=0.5$, with attention to Fig. (3)).

The main point to be noted in this subsection is that with $n_{\rm h}=1.5$, $\bar{A}=608.9$, $\gamma_{\rm h}=0.3$, and all other parameters the same as before, we do indeed observe superior convergence over a significant range of step sizes: see Fig. (11) for FE models with 10 and 30 elements, and compare with Fig. (7). For the case with 10 elements, the convergence is essentially quadratic. With 30 elements the convergence is slightly slower than quadratic, but much faster than linear. Note that these estimates are from numerics only: analytical estimates are not available.

As mentioned above some of the difficulty with accurate integration of hysteretically damped structures comes from the zero crossings of the hysteretic variable $z$ itself. In this paper, for simplicity, we have avoided special treatment of these zero crossings. However, the results of this subsection show that the effect of these zero crossings is milder if $n_{\rm h}$ is larger.

We now turn to using results obtained from this semi-implicit integration method to develop data-driven lower order models of the hysteretically damped structure.

In our beam system, we have 2 degrees of freedom per node (one displacement and one rotation). Let the total number of differential equations being solved, in first order form, be called $N_{\rm D}$. If there are $n_{\rm e}$ elements, we have $N_{\rm D}=4n_{\rm e}+n_{\rm e}n_{\rm g}$ first order differential equations. In a reduced order model with $r$ modes and $m$ Gauss points, we have $N_{\rm D}=2r+m$ first order differential equations. Although the size of the problem is reduced significantly in this way, the accuracy can be acceptable.

For demonstration, we consider two FE models of a cantilever beam, with hysteresis, as follows:

• (i)

  100 elements with 3 Gauss points each ($N_{\rm D}=700$).

• (ii)

  150 elements with 3 Gauss points each ($N_{\rm D}=1050$).

For each of the two systems above, the datasets used for selecting the subset of hysteretic states (or Gauss points) were generated by solving the full systems 60 times each for the time interval 0 to 1 with random initial conditions and time step $h=10^{-4}$, exciting only the first 3 modes ($r=3$). Data from each of these 60 solutions were retained at 1000 equispaced points in time ($N_{t}=1000$; recall subsection 4.2.1).

All reduced order model (ROM) results below are for $r=3$ retained modes. Further, the number of hysteresis Gauss points is denoted by $m$. Results for $m=n_{\rm e}$ are shown in Fig. (12a,12c).

We now quantitatively assess the accuracy of the ROMs. To this end, we write $y^{(m)}_{\rm tip(ROM)}(t)$ and $y_{\rm tip(FM)}(t)$ for the ROM and full model outputs respectively. We then compute the error measure

for different values of $m$ and for an integration time interval $[0,1]$ (with $h=10^{-4}$).

The variation of $\mathcal{E}_{\rm rms}$ with $m$, for both models, is shown in Fig. (13). This error measure is not normalized. If we wish to normalize it, then we should use the quantity obtained when we set $y^{(m)}_{\rm tip(ROM)}(t)$ identically to zero: that quantity is 0.006 for both models. Thus, reasonable accuracy is obtained for $m\geq 50$, and good accuracy (about 1% error) is obtained only for $m>100$. These numbers refer to a specific case with random initial conditions; however, Fig. (13) is representative for other, similar, initial conditions (details omitted).

Finally, we report on run times needed on a modest laptop computer with the different levels of modeling described above. A representative initial conditions similar to the above was selected. The simulation duration was chosen to be 1 second, and the time step taken was $10^{-4}$ seconds. The following results were obtained after averaging run times for 4 runs each. (i) The full simulation with our algorithm took an average of 27.22 seconds per run. (ii) After projecting on to 3 normal modes but retaining all 450 hysteretic states, each run took an average of 1.3 seconds. (iii) Finally, with 3 normal modes and hysteretic states reduced from 450 to 150, each run took an average of 0.84 seconds. Note that we could save further run time by taking larger time steps after model order reduction and/or retaining even fewer hsyteretic states.