Extreme events and their optimal mitigation in nonlinear structural systems excited by stochastic loads: Application to ocean engineering systems

TL;DR

This paper introduces an efficient numerical method to accurately predict extreme response statistics of nonlinear structural systems under stochastic loads, with applications to ocean engineering, enabling optimized mitigation strategies.

Contribution

A novel numerical approach combining statistical linearization and effective measures for fast, accurate extreme response prediction in nonlinear systems under stochastic excitation.

Findings

The method captures heavy-tail response statistics orders of magnitude faster than direct simulations.

Optimized asymmetric nonlinear springs outperform traditional linear dampers and cubic energy sinks.

Application to ocean engineering systems demonstrates effective mitigation of extreme events.

Abstract

We develop an efficient numerical method for the probabilistic quantification of the response statistics of nonlinear multi-degree-of-freedom structural systems under extreme forcing events, emphasizing accurate heavy-tail statistics. The response is decomposed to a statistically stationary part and an intermittent component. The stationary part is quantified using a statistical linearization method while the intermittent part, associated with extreme transient responses, is quantified through i) either a few carefully selected simulations or ii) through the use of effective measures (effective stiffness and damping). The developed approach is able to accurately capture the extreme response statistics orders of magnitude faster compared with direct methods. The scheme is applied to the design and optimization of small attachments that can mitigate and suppress extreme forcing events delivered to a primary structural system. Specifically, we consider the problem of suppression of extreme responses in two prototype ocean engineering systems. First, we consider linear and cubic springs and perform parametric optimization by minimizing the forth-order moments of the response. We then consider a more generic, possibly asymmetric, piecewise linear spring and optimize its nonlinear characteristics. The resulting asymmetric spring design far outperforms the optimal cubic energy sink and the linear tuned mass dampers.

Tables

Table 1. Table 1 : Parameters and relevant statistical quantities for the suspended seat system.

$m_s$	$1$	$m_a$	$0.05$
${\lambda }_s$	$0.01$	${\lambda }_a$	$0.021$
$k_s$	$1$	$k_a$	$0$
—	—	$c_a$	$3.461$
$T_{\alpha }$	$5000$	${\sigma }_{\eta }$	$0.0227$
${\mu }_{\alpha }=7\times {\sigma }_{\dot{h}}$	$0.1$	$q$	$1.582\times {10}^{-4}$
${\sigma }_{\alpha }={\sigma }_{\dot{h}}$	$0.0141$	${\sigma }_h$	$0.0063$
${\mathbb{P}}_r^x$	$0.0214$	${\mathbb{P}}_r^v$	$0.0107$
${\mathbb{P}}_r^{\dot{x}}$	$0.0210$	${\mathbb{P}}_r^{\dot{v}}$	$0.0100$
${\mathbb{P}}_r^{\ddot{x}}$	$0.0212$	${\mathbb{P}}_r^{\ddot{v}}$	$0.0096$

Table 2. Table 2 : Parameters and relevant statistical quantities for the suspended deck-seat system.

$m_h$	$1$	$m_s$	$0.05$	$m_a$	$0.05$
${\lambda }_h$	$0.01$	${\lambda }_s$	$0.1$	${\lambda }_a$	$0.035$
$k_h$	$1$	$k_s$	$1$	$k_a$	$0$
—	—	—	—	$c_a$	$5.860$
$T_{\alpha }$	$5000$	${\mu }_{\alpha }=7\times {\sigma }_{\dot{h}}$	$0.1$	$q$	$1.582\times {10}^{-4}$
${\sigma }_{\eta }$	$0.0232$	${\sigma }_{\alpha }={\sigma }_{\dot{h}}$	$0.0141$	${\sigma }_h$	$0.0063$
${\mathbb{P}}_r^y$	$0.0245$	${\mathbb{P}}_r^x$	$0.0247$	${\mathbb{P}}_r^v$	$0.0162$
${\mathbb{P}}_r^{\dot{y}}$	$0.0234$	${\mathbb{P}}_r^{\dot{x}}$	$0.0202$	${\mathbb{P}}_r^{\dot{v}}$	$0.0161$
${\mathbb{P}}_r^{\ddot{y}}$	$0.0238$	${\mathbb{P}}_r^{\ddot{x}}$	$0.0081$	${\mathbb{P}}_r^{\ddot{v}}$	$0.0146$

Table 3. Table 3 : Suspended seat system parameters.

$m_s$	$1$	$m_a$	$0.05$
${\lambda }_s$	$0.01$	$k_s$	$1$
$T_{\alpha }$	$5000$	$-$	$-$
${\mu }_{\alpha }=7\times {\sigma }_{\dot{h}}$	$0.1$	$q$	$1.582\times {10}^{-4}$
${\sigma }_{\alpha }={\sigma }_{\dot{h}}$	$0.0141$	${\sigma }_h$	$0.0063$

Table 4. Table 4 : Suspended deck-seat system parameters.

$m_h$	$1$	$m_s$	$0.05$
$m_a$	$0.05$	${\lambda }_h$	$0.01$
$k_h$	$1$	${\lambda }_s$	$0.1$
$k_s$	$1$	$T_{\alpha }$	$5000$
${\mu }_{\alpha }=7\times {\sigma }_{\dot{h}}$	$0.1$	$q$	$1.582\times {10}^{-4}$
${\sigma }_{\alpha }={\sigma }_{\dot{h}}$	$0.0141$	${\sigma }_h$	$0.0063$

Equations

$$m_s\ddot{x}$$

$$m_a\ddot{v}$$

$$m_h\ddot{y}$$

$$m_s\ddot{x}$$

$$m_a\ddot{v}$$

$$\ddot{\xi }(t)=\ddot{h}(t)+\sum _{i=1}^{N(t)}{\alpha }_i\delta (t-{\tau }_i),0<t\le T,$$

$$S_{hh}(\omega)=q\frac{1}{\omega^{5}}\exp\biggl{(}-\frac{1}{\omega^{4}}\biggr{)},$$

$$\ddot{x}+\lambda \dot{x}+kx=\ddot{\xi }(t),$$

$$p_x(r)=p_{x_b}(r)(1-{\mathbb{P}}_r)+p_{x_r}(r){\mathbb{P}}_r,$$

$${\ddot{x}}_b+\lambda {\dot{x}}_b+k{x}_b=\ddot{h}(t).$$

$$S_{x_b{x}_b}(\omega )=\frac{{\omega }^4{S}_{hh}(\omega )}{{\left(k-{\omega }^2\right)}^2+(\lambda \omega {)}^2},{\sigma }_{x_b}^2={\int }_0^{\infty }{S}_{x_b{x}_b}(\omega )d\omega .$$

$$p_{x_r}(r)=\int p_{x_r\mid \eta }(r\mid n)p_{\eta }(n)dn,$$

$$\eta \equiv {\dot{x}}_b+\alpha \sim \mathcal{N}({\mu }_{\alpha },{\sigma }_{{\dot{x}}_b}^2+{\sigma }_{\alpha }^2).$$

$$p_{x_{r}\mid\eta}(r\mid n)=\operatorname{\text{Hist}}\bigl{\{}x_{r\mid\eta}(t\mid n)\bigr{\}},\quad t\in[0,\tau_{e}],$$

$$\displaystyle x_{r\mid\eta}(t\mid n)=\frac{n}{2\omega_{o}}\biggl{(}\mathrm{e}^{-(\zeta\omega_{n}-\omega_{o})t}-\mathrm{e}^{-(\zeta\omega_{n}+\omega_{o})t}\biggr{)}.$$

$$x_{r}(\tau_{e})=\rho_{c}\,\max\big{\{}|x_{r}|\big{\}},$$

$${\mathbb{P}}_r={\nu }_{\alpha }{\tau }_e={\tau }_e/T_{\alpha }.$$

$$\displaystyle p_{z}(r)=\frac{1-\nu_{\alpha}\tau_{e}^{z}}{\sigma_{z_{b}}\sqrt{2\pi}}\exp\biggl{(}-\frac{r^{2}}{2\sigma_{z_{b}}^{2}}\biggr{)}+\nu_{\alpha}\tau_{e}^{z}\int^{\infty}_{0}\operatorname{\text{Hist}}\bigl{\{}z_{r\mid\eta}(t\mid n)\bigr{\}}p_{\eta}(n)\,dn,$$

$$m_s\ddot{x}$$

$$m_a\ddot{v}$$

$$m_s{C}_{xx}^{\prime \prime }+{\lambda }_s{C}_{xx}^{\prime }+k_s{C}_{xx}+{\lambda }_a\left(C_{xx}^{\prime }-C_{vx}^{\prime }\right)+k_a\left(C_{xx}-C_{vx}\right)+c_a\overline{{\left(x(t)-v(t)\right)}^{3}x(s)}=-m_s{C}_{hx}^{\prime \prime },$$

$$m_s{C}_{xv}^{\prime \prime }+{\lambda }_s{C}_{xv}^{\prime }+k_s{C}_{xv}+{\lambda }_a\left(C_{xv}^{\prime }-C_{vv}^{\prime }\right)+k_a\left(C_{xv}-C_{vv}\right)+c_a\overline{{\left(x(t)-v(t)\right)}^{3}v(s)}=-m_s{C}_{hv}^{\prime \prime },$$

$$m_s{C}_{xh}^{\prime \prime }+{\lambda }_s{C}_{xh}^{\prime }+k_s{C}_{xh}+{\lambda }_a\left(C_{xh}^{\prime }-C_{vh}^{\prime }\right)+k_a\left(C_{xh}-C_{vh}\right)+c_a\overline{{\left(x(t)-v(t)\right)}^{3}h(s)}=-m_s{C}_{hh}^{\prime \prime },$$

$$m_a{C}_{vx}^{\prime \prime }+{\lambda }_a\left(C_{vx}^{\prime }-C_{xx}^{\prime }\right)+k_a\left(C_{vx}-C_{xx}\right)+c_a\overline{{\left(v(t)-x(t)\right)}^{3}x(s)}=-m_a{C}_{hx}^{\prime \prime },$$

$$m_a{C}_{vv}^{\prime \prime }+{\lambda }_a\left(C_{vv}^{\prime }-C_{xv}^{\prime }\right)+k_a\left(C_{vv}-C_{xv}\right)+c_a\overline{{\left(v(t)-x(t)\right)}^{3}v(s)}=-m_a{C}_{hv}^{\prime \prime },$$

$$m_a{C}_{vh}^{\prime \prime }+{\lambda }_a\left(C_{vh}^{\prime }-C_{xh}^{\prime }\right)+k_a\left(C_{vh}-C_{xh}\right)+c_a\overline{{\left(v(t)-x(t)\right)}^{3}h(s)}=-m_a{C}_{hh}^{\prime \prime }.$$

$$\overline{{\left(x(t)-v(t)\right)}^{3}x(s)}=\left(3{\sigma }_x^2-6{\sigma }_{xv}+3{\sigma }_v^2\right)C_{xx}-\left(3{\sigma }_x^2-6{\sigma }_{xv}+3{\sigma }_v^2\right)C_{vx},$$

$$\overline{{\left(x(t)-v(t)\right)}^{3}v(s)}=\left(3{\sigma }_x^2-6{\sigma }_{xv}+3{\sigma }_v^2\right)C_{xv}-\left(3{\sigma }_x^2-6{\sigma }_{xv}+3{\sigma }_v^2\right)C_{vv},$$

$$\overline{{\left(x(t)-v(t)\right)}^{3}h(s)}=\left(3{\sigma }_x^2-6{\sigma }_{xv}+3{\sigma }_v^2\right)C_{xh}-\left(3{\sigma }_x^2-6{\sigma }_{xv}+3{\sigma }_v^2\right)C_{vh}.$$

$$S_{xx}(\omega ;{\sigma }_x^2,{\sigma }_{xv},{\sigma }_v^2)$$

$$S_{vv}(\omega ;{\sigma }_x^2,{\sigma }_{xv},{\sigma }_v^2)$$

$$S_{xv}(\omega ;{\sigma }_x^2,{\sigma }_{xv},{\sigma }_v^2)$$

$$S_{xh}(\omega ;{\sigma }_x^2,{\sigma }_{xv},{\sigma }_v^2)$$

$$S_{vh}(\omega ;{\sigma }_x^2,{\sigma }_{xv},{\sigma }_v^2)$$

$$\mathcal{A}(\omega ;{\sigma }_x^2,{\sigma }_{xv},{\sigma }_v^2)=$$

$$\mathcal{B}(\omega ;{\sigma }_x^2,{\sigma }_{xv},{\sigma }_v^2)=$$

$$\mathcal{C}(\omega ;{\sigma }_x^2,{\sigma }_{xv},{\sigma }_v^2)=$$

$${\sigma }_x^2={\int }_0^{\infty }{S}_{xx}(\omega ;{\sigma }_x^2,{\sigma }_{xv},{\sigma }_v^2)d\omega ,{\sigma }_{xv}={\int }_0^{\infty }{S}_{xv}(\omega ;{\sigma }_x^2,{\sigma }_{xv},{\sigma }_v^2)d\omega ,{\sigma }_v^2$$

$$p_{z}(r)=\frac{1-\nu_{\alpha}\tau_{e,\text{dis}}^{z}}{\sigma_{z_{b}}\sqrt{2\pi}}\exp\biggl{(}-\frac{r^{2}}{2\sigma_{z_{b}}^{2}}\biggr{)}+\nu_{\alpha}\tau_{e,\text{dis}}^{z}\int^{\infty}_{0}\operatorname{\text{Hist}}\bigl{\{}z_{r\mid\eta}(t\mid n)\bigr{\}}p_{\eta}(n)\,dn,$$

Full text

Keywords

Impact mitigation in nonlinear structural systems; Response under extreme events; Optimization and design under stochastic loads; Nonlinear Energy Sinks; Optimization of suspended seats and decks in high speed craft motion.

1 Introduction

For a plethora of structural systems it is essential to specify their reliability under uncertain environmental loading conditions and most importantly provide design guidelines using knowledge of their response characteristics. This involves accurate estimation of the structural systems probabilistic response. Environmental loads are typically random by nature and are likely to include intermittently occurring components of an extreme magnitude, representing abnormal environmental events or conditions. Although extreme loadings occur with lower probability than typical conditions, their impact is significant and cannot be neglected since these events determine the systems behavior away from the average operating conditions, which are precisely the conditions that are important to quantify for safe assessment and design. Important examples include mechanical and ocean engineering systems. High speed crafts in rough seas [29, 28], wave impacts on fixed or floating offshore platforms and ship capsize events [23, 2, 24, 17, 14], vibrations of buildings or bridge structures due to earthquakes or strong wind excitations [15, 4, 16, 37] are just a few examples where extreme responses occur infrequently but are critical in determining the overall systems reliability.

Numerous research endeavors have been dedicated on the effective suppression and rapid dissipation of the energy associated with extreme impacts on structures. Many of these schemes rely on linear configurations, known as tuned mass damper (TMD) and result in a halving of the resonance frequency. Although the mitigation performance is highly effective when most of the energy is concentrated at the characteristic frequency of the system, their effectiveness drastically drops if there is a mistuning in frequency. Moreover, it is not clear how these configurations perform in the presence of rare impulsive loads. Many of these limitations can be overcome by utilizing small attachments coupled with the primary system through nonlinear springs, also known as nonlinear energy sinks (NES). If carefully chosen these nonlinear attachments can lead to robust, irreversible energy transfer from the primary structure to the attachment and dissipation there [39, 40]. The key mechanism behind the efficient energy dissipation in this case is the targeted energy transfer phenomenon which is an essentially nonlinear mechanism and relies primarily on the energy level of the system, rather then the resonant frequency [41, 11]. Such configurations have been proven to be successful on the mitigation of deterministic impulsive loads on large structures  [33, 34, 18] and their performance has been measured through effective nonlinear measures such as effective damping and stiffness  [32, 26].

Despite their success, nonlinear configurations have been primarily developed for deterministic impulsive loads. To quantify and optimize their performance in the realistic settings mentioned previously it is essential to understand their effects on the statistics of the response and in particular in the heavy tails of the probability distribution function (PDF). However, quantifying the PDF of nonlinear structures under random forcing containing impulsive type extreme events, poses many challenges for traditional methods. Well established approaches for determining the statistics of nonlinear dynamical systems include the Fokker-Planck equation [36, 35], the joint response-excitation method [31, 42, 10, 1], Gaussian closure schemes, moment equation or cumulant closure methods [3, 43], the Polynomial-Chaos approach [44], and stochastic averaging methods [46]. For systems associated with heavy tails, however, these methods either cannot capture the statistics of rare and extreme type events due to inherent limitations [19] or are far too computationally expensive in practice, even for low-dimensional systems [20, 6]. Alternatively, one can study the statistics of the extreme events alone (by ignoring the background ‘non-extreme’ forcing fluctuations) through a Poisson process representation and then analyze the response using the generalized Fokker-Planck or Kolmogorov-Feller equations [35], which governs the evolution of the corresponding PDF, or by applying the path integral formalism [13, 8], or even through special stochastic averaging techniques [45]. While attractive, these ideas lead, in general, to analytical results for a very limited number of special cases. Besides, it is still an important aspect to account for the background random fluctuations in the forcing term in order to fully characterize the systems overall probabilistic properties (e.g. this is important in order to fully determine all the moments of the response). Moreover, even though the background forcing component does not directly correspond to extreme events, the background term may have important consequences for the initiation of intermittent type extreme responses [21].

In this work we consider the problem of nonlinear structural systems under general time-correlated stochastic forcing that includes extreme, impulsive type random events. We address two important challenges related to this problem. The first is the development of a fast and accurate estimation method for the response statistics, expressed through the PDF, with emphasis on the accurate estimation of the tail form (events far away from the mean). The second is the design and parameter optimization of small attachments that can mitigate or suppress the effects of the extreme forcing events on the system response while they also improve the system behavior during the regular regime. The two problems are connected since extreme event suppression is directly reliant upon a fast and accurate estimation method for the response pdf under different designs or parameters. Indeed, without a fast and reliable method to evaluate response statistics, in particular tail statistics, optimization cannot be performed because of the inherent computational cost associated with typical quantification methods such as Monte-Carlo. This aspect highlights the practical utility of the proposed fast PDF estimation scheme. We will illustrate the pdf estimation method and shock mitigation design analysis throughout the manuscript with a practical motivating prototype system related to high speed vehicle motion in rough seas, however we emphasize the proposed method broad applicability.

The probabilistic quantification scheme formulated here is based on the most general probabilistic decomposition-synthesis framework [21, 22], that has recently been applied in linear systems subjected to stochastic forcing containing extreme events [9] and can be used to efficiently estimate the PDF for the response displacement, velocity, and acceleration. We begin by formulating the response pdf quantification method (developed for linear multi-degree-of-freedom (MDOF) systems in [9]) for the case of nonlinear MDOF systems. This is achieved by combining the probabilistic decomposition-synthesis framework [21, 22] with the statistical linearization method  [30]. The scheme circumvents the rare-event problem and enables rapid design and optimization in the presence of extreme events. We emphasize the statistical accuracy of the derived scheme, which we have validated through extensive comparisons with direct Monte-Carlo simulations. Next, we consider two prototype ocean engineering systems and perform a quantitative comparison of the performance of TMD and NES, evaluating their effectiveness at shock suppression under stochastic excitation containing extreme events. Finally, we perform optimization on a very generic, possibly asymmetric family of piecewise linear springs. Previous endeavors in the context of single-sided vibro-impact NES have shown that asymmetries in the NES can improve the shock mitigation properties (see [34]). In agreement with these results, our optimization scheme leads to the derivation of a new asymmetric NES which significantly improves the shock mitigation properties of the system in the realistic setting of stochastic excitation.

The paper is structured as follows. In section 2 we describe the prototype models for high speed craft motion that we utilize throughout the paper as practically relevant example. Next, in section 3 we provide a brief review of the probabilistic decomposition-synthesis (PDS) framework for the response pdf quantification of a linear single-degree-of-freedom system subject to a random forcing term containing extreme impulse type events. Section 4 describes the proposed general semi-analytical PDF estimation method for nonlinear MDOF structures and also includes a section on quantifying the conditionally rare response via the effective stiffness and damping framework. In section 5 we present the mitigation of extreme events analysis on the prototype high speed craft designs for both TMD and cubic NES attachments. Next, in section 6 we propose a new piecewise linear and asymmetric NES design that we optimize for extreme event mitigation. Finally in section 7 we offer concluding remarks.

2 Prototype models for high speed vehicle motion in rough seas

Here we describe the prototype models that we apply the quantification method for extreme event analysis and optimization. Specifically, we model the motion of a high-speed craft in random seas through two prototype systems: one being a two-degree-of-freedom system consisting of a suspended seat attached to the hull and the second being a three-degree-of-freedom system where the seat is attached to a suspended deck, which is attached on the hull; both prototypes contain a small linear or nonlinear energy sink (NES) vibration absorber.

2.1 2DOF Suspended seat system

In figure 1 we illustrate the first model consisting of a linear primary structure under base excitation that is attached to a small oscillator connected through a nonlinear spring (with cubic nonlineariry). This is a prototype system modeling the suspended seat of a high speed craft [25, 5]. The vibration absorber is attached to the seat with the aim to minimize ocean wave impacts on the operator of the vehicle and naturally we require that the attachment mass is much lower than the seat mass (i.e. $m_{a}<0.1m_{s}$). The equation of motion for this two-degree-of-freedom system is given by:

[TABLE]

where $x,v$ are the relative displacements of the seat response and attachment response, respectively, with reference to the base motion $\xi(t)$ (that is, $x=\hat{x}-\xi$ and $v=\hat{v}-\xi$).

2.2 3DOF Suspended deck-seat system

The second prototype system is a suspended deck design for a high speed craft [38, 27, 12] and is illustrated in figure 2. In this case, the vibration absorber is attached to the suspended deck. The attachment mass is comparable to the seat mass and both are considerably smaller than the deck (i.e. $m_{a}\simeq m_{s}<0.1m_{h}$). The governing equations for this three-degree-of-freedom system are given by:

[TABLE]

where, again, $x,y,v$ are the relative displacements of the seat response, the deck response and the attachment response, respectively, with reference to the base motion $\xi(t)$.

In both prototypes the aim of the vibration mitigating attachment is to minimize extreme impacts on the seat attachment as this represents an operator on the vehicle. We first examine the case of tuned-mass damper vibration absorber $k_{a}\neq 0,c_{a}=0$ and the essentially nonlinear energy sink absorber $k_{a}=0,c_{a}\neq 0$, that has been studied extensively in the context of shock mitigation [40]. In the last section we will examine the performance of an asymmetric, piecewise linear, spring.

2.3 The structure of the intermittently extreme stochastic forcing

Motivated by the ocean engineering systems in section 2, we consider base motion of the form,

[TABLE]

In the expression above, $h(t)$ denotes a zero-mean smooth motion characterized by a Pierson-Moskowitz spectrum,

[TABLE]

where $q$ controls the magnitude of the motion. The second term in equation 3 describes rare and extreme impulses in terms of a random impulse train ($\delta(\,\cdot\,)$ is a unit impulse), occurring due to slamming events. For this component, $N(t)$ is a Poisson counting process that represents the number of impulses that arrive in the time interval $0<t\leq T$, $\alpha$ is the impulse magnitude, which we assume is normally distributed with mean $\mu_{\alpha}$ and variance $\sigma_{\alpha}^{2}$, and the constant arrival rate is given by $\nu_{r}$. We take the impulse magnitude as being $\beta$-times larger than the standard deviation of the excitation velocity $\dot{h}(t)$: $\mu_{\alpha}=\beta\ \sigma_{\dot{h}},\text{ with }\beta\ >1.$

3 Review of the probabilistic decomposition-synthesis (PDS) method

We first provide a brief review of the semi-analytical response quantification method for a linear single-degree-of-freedom system [9] subjected to stochastic excitation containing rare events. The purpose of this section is to provide a self-contained review of the core ideas, since the scheme for nonlinear structural systems that is described in the following section depends upon these concepts.

Consider the following linear system

[TABLE]

$k$ is the stiffness, $\lambda$ is the damping, and $\zeta=\lambda/2\sqrt{k}$ is the damping ratio. Despite the simplicity of this system, the structure of the statistical response may be significantly complex and posses heavy-tails.

The framework to estimate the response PDF of equation 5 is the probabilistic decomposition-synthesis (PDS) method [21]. The basic idea is to decouple the rare events regime from the background fluctuations and then quantify the statistics of the two components separately. The results are then synthesized to obtain the full response PDF by using the total probability law:

[TABLE]

where $p_{x_{b}}(r)$ is the conditional PDF due to the smooth motion of the base, $p_{x_{r}}(r)$ is the conditional PDF due to the extreme impacts and $\mathbb{P}_{r}$ is the overall probability that the system operates in the extreme events regime.

3.1 Background response PDF

We first obtain the statistical response of the system under the condition that only the background (smooth) forcing component is acting. We have,

[TABLE]

In this case the analysis is particularly simple since the system is linear and time-invariant and the response PDF, $p_{x_{b}}$, is a zero-mean Gaussian. The spectral density of the response displacement and the variance are given by:

[TABLE]

The computations for the response velocity and acceleration can be similarly obtained.

3.2 Numerical histogram for rare events

The next step is to compute the rare event distribution $p_{x_{r}}$ and the rare event probability $\mathbb{P}_{r}$. Specifically, the rare event distribution can be written as,

[TABLE]

where $p_{\eta}(n)$ is the distribution of the impulse magnitude, and $p_{x_{r}\mid\eta}$ is the conditional PDF of the response for an impact of magnitude $\eta$.

It is important to note that once an impulse of magnitude $\alpha$ hits the system, the momentum of the system right after the impact would be $\dot{x}_{b}+\alpha$, since the momentum of the system right before the impact is $\dot{x}_{b}$. As these two variables are both Gaussian distributed and independent, their sum is also Gaussian distributed and is given by,

[TABLE]

We estimate the conditional PDF $p_{x_{r}\mid\eta}(r\mid n)$ by the numerically computed histogram:

[TABLE]

where $\tau_{e}$ is the typical duration of the rare event (see next subsection) and the conditional response $x_{r\mid\eta}$ is given by,

[TABLE]

The conditionally extreme event distribution for velocity and acceleration are derived in a similar fashion.

3.3 Numerical estimation of the rare event probability

In order to compute the histogram of a rare impulse event, the duration of a rare response needs to be obtained numerically. We define the typical duration of a rare response by

[TABLE]

where $\rho_{c}=0.1$, or in other words, the histogram is taken over the time it takes for the system response to decay to 10% of its maximum value. The absolute value of the maximum of the response needs to be estimated numerically.

Once this rare event duration has been specified, we can also obtain the probability of a rare event by

[TABLE]

Note that the extreme event duration for the displacement $\tau_{e}^{x}$, velocity $\tau_{e}^{\dot{x}}$, and acceleration $\tau_{e}^{\ddot{x}}$ are in generally different.

3.4 Semi-analytical response probability distributions

With the description above, we obtain the response PDF using the total probability law. The resulting response PDF takes the form,

[TABLE]

where the argument $z$ is either $x$, $\dot{x}$, or $\ddot{x}$. The validity of this approximation has been thoroughly verified in [9].

4 PDF quantification method for nonlinear MDOF systems

Here we formulate the probabilistic-decomposition method for multi-degree-of-freedom, nonlinear mechanical systems. There are some important differences with respect to the case of linear systems studied in [9]. Firstly, for the background component the system nonlinearities can be important and to this end we must utilize an appropriate statistical quantification method. Here we employ the statistical linearization approach [30]. Secondly, to characterize the statistics in the rare event regime it is even more crucial to take into account the nonlinear properties of the system, since these control the shock mitigation capabilities of the attachment.

To achieve this we use two alternative approaches. The first one is based on the direct simulation of the system for a range of initial conditions corresponding to all possible impact magnitudes. The second is based on the notion of effective stiffness and damping [32], which are measures that characterize the system response under various excitation magnitudes taking into account the presence of the nonlinear attachment. We provide comparisons with direct Monte-Carlo simulations to demonstrate the accuracy of both approaches. We first present the analysis for the background component

4.1 Quantification of the response pdf for the background component

For the background regime, we must account for nonlinearities and their interaction with the background part of the excitation. We use the statistical linearization method, since we are only interested in resolving the low-order statistics of the background response of the system (the rare events component defines the tails of the PDF).

Consider the response of the suspended seat problem, section 2.1, under the excitation term $\ddot{h}(t)$:

[TABLE]

We first multiply the above two equations by $x(s)$, $v(s)$, $h(s)$ at different time instant $s\neq t$, and take ensemble averages to write the resulting equations in terms of covariance functions.

[TABLE]

Here ′ indicates the partial differentiation with respect to the time difference $\tau=t-s$. We then apply Isserlis’ theorem based on the Gaussian process approximation for response to express the fourth-order moments in terms of second-order moments [7].

[TABLE]

This leads to a set of linear equations in terms of the covariance functions. Thus, the Wiener-Khinchin theorem can be applied to write the equations in terms of the power spectrum, giving

[TABLE]

where,

[TABLE]

At this point $\sigma_{x}^{2}$, $\sigma_{v}^{2}$, and $\sigma_{xv}$ are still unknown, but can be determined by integrating both sides of equations 27, 28, and 29 and forming the following system of equations:

[TABLE]

By solving the above we find $\sigma_{x}^{2},$ $\sigma_{v}^{2},$ and $\sigma_{xv}$. This procedure determines the Gaussian PDF approximation for the background regime response. Further details regarding the special case of a linear attachment and the analysis for the suspended deck-seat problem can be found in appendix A.

4.2 Quantification of the response pdf for the extreme event component

We are going to utilize two alternative methods for the quantification of the statistics in the extreme event regime. The first approach is to obtain the conditional statistics based on direct simulations of the system response. The second method is utilizing effective measures [32] that also characterize the system nonlinear response in the presence of attachments.

4.2.1 Rare response PDF using direct simulations of the system under impulsive excitation

To compute the conditionally extreme distribution $p_{x_{r}}$ and the probability of rare events $\mathbb{P}_{r}$ we follow the steps described in  algorithm 1, which provides a high-level description for a single mode. The procedure is repeated for each degree of freedom of interest (in this case it is more efficient to simply store all the impulse realizations and then run the procedure for each degree of freedom of interest). We emphasize that the numerical simulation of impulse response for nonlinear systems is efficient, since the integrations are necessarily short due the impulsive nature of the forcing and the condition on the rare event end time in equation 13. Moreover, throughout these simulations we do not take into account the background excitation since this is negligible compared with the effect of the initial conditions induced by the impact.

Algorithm 1 Calculation of $\mathbb{P}_{r}$ and $p_{x_{r}}(r)=\int p_{x_{r}\mid\eta}(r\mid n)p_{\eta}(n)\,dn$.

1:discretize $p_{\eta}(n)$

2:for all $n$ values over the discretization $p_{\eta}$ do

3:     solve ODE system for $x^{n}(t)$ under impulse $n$, neglecting $\ddot{h}$

4:     ${\tau_{e}}^{n}\leftarrow\{t_{e}\mid\rho_{c}\max_{t}{\lvert x^{n}(t)\rvert}=x^{n}(t_{e})\}$       // we set $\rho_{c}=0.1$

5:     p_{x_{r}\mid\eta}^{\,n}\leftarrow\operatorname{\text{Hist}}\bigl{\{}x^{n}(t)\;\big{|}\;t\in[0,{\tau_{e}}^{n}]\bigr{\}}

6:end for

7:$p_{x_{r}}\leftarrow\int p_{x_{r}\mid\eta}^{\,n}\,p_{\eta}^{\,n}$

8:$\tau_{e}\leftarrow\int\tau_{e}^{\,n}\,p_{\eta}^{\,n}$

9:$\mathbb{P}_{r}\leftarrow\nu_{\alpha}\tau_{e}$

10:output: $\mathbb{P}_{r},p_{x_{r}}$

Comparison with Monte-Carlo Simulations

The full response PDF is composed using the total probability law,

[TABLE]

where $z$ is either the displacement, velocity or acceleration of the seat/attachment response. We utilize a shifted Pierson-Moskowitz spectrum $S_{hh}(\omega-1)$ for the background forcing term in order to avoid system resonance. Details regarding the Monte-Carlo simulations are provided in appendix B.

In figure 3 we show comparisons for the suspended seat problem with parameters and relevant statistical quantities given in table 1. In figure 4 we also show comparisons for the suspended deck-seat problem with parameters and relevant statistical quantities in table 2. For both cases the adopted quantification scheme is able to compute the distributions for the quantities of interest extremely fast (less than a minute on a laptop), while the corresponding Monte-Carlo simulations take order of hours to complete.

Note that our method is able to capture the complex heavy tail structure many standard deviations away from the mean (dashed vertical line denotes 1 standard deviation). We emphasize that similar accuracy is observed for a variety of system parameters that satisfy the assumptions on the forcing. The close agreement validates that the proposed scheme is applicable and can be accurately used for system optimization and design.

4.3 Rare response PDF using effective measures

Here we describe an alternative technique to quantify the rare event PDF component using the effective stiffness and damping framework described in [32]. These effective measures express any degree-of-freedom of the coupled nonlinear system, for a given initial energy level, as an equivalent linear single-degree-of-freedom system. Specifically, these effective measures correspond to the values of damping and stiffness for a linear system that has (for the same initial conditions) a response that is as close as possible to that of the original system, in the mean square sense.

We focus on the suspended seat problem to illustrate this strategy. It should be pointed out that the accuracy and applicability of this approach has some limitations:

• •

The accurate estimation of the PDF requires the knowledge of the effective measures over a sufficiently large range of initial impulses.

• •

The motion of the system should have an oscillatory character so that it can be captured by effective measures.

• •

The statistics of the attachment motion cannot be obtained directly from the effective measures.

To derive the PDF in the rare event regime we reduce the system to an effective linear system for the degree-of-freedome of interest. Consider the suspended seat system under an impulse,

[TABLE]

with initial conditions, at an arbitrary time say $t_{0}=0$,

[TABLE]

To determine the effective linear system for this system, we follow the strategy in [32] and compute the effective stiffness and damping:

[TABLE]

where $\langle\,\cdot\,\rangle$ denotes spline interpolation of the local maxima of the time series. We can then compute the weighted-average effective stiffness and damping:

[TABLE]

With the weighted-average effective measures we rewrite the original two-degree-of-freedom system during rare events into an equivalent linear single-degree-of-freedom system with coefficients that depend on the initial impact (or the initial energy level of the system):

[TABLE]

Using the effective system in equation 41 we can obtain the conditionally rare PDF using the analysis for the linear system in section 3.1. The damping ratio and natural frequency now become functions of the initial impact, $n$:

[TABLE]

Subsequently, the PDF is obtained by taking a histogram of

[TABLE]

In figure 5 (top) we present the suppression of the probability for large motions of the primary structure due to the presence of the NES (parameters given in table 1). This suppression is fully expressed in terms of the effective damping measure shown in the lower plot. Note that the suppression of the tail begins when the effective damping attains values larger than one.

We emphasize that in the context of effective measures the motion of the system is assumed oscillatory. Motions with radically different characteristics will not be captured accurately from the last representation and the resulted histograms will not lead to an accurate representation of the tail. This problem is, in general, circumvented if we employ the first approach for the computation of the conditional PDF during extreme impacts. On the other hand, the advantage of the second approach is that we can interpret the form of the tail in the various regimes with respect to the properties of the effective measures (figure 5). This link between dynamics (effective measures) and statistics (heavy tail form) is important for the design process of the NES.

Comparison with Monte-Carlo simulations

Here we compare the PDS method combined with the effective measures with direct Monte-Carlo simulations. In figure 6 we show the response PDF for the primary structure for parameters given in table 1. Details regarding the Monte-Carlo computations are provided in section 4.2. We observe that the PDS method utilizing effective measures performs satisfactorily over a wide range similarly with the first general scheme, based on individual trajectories computation.

4.4 Quantification of the absolute response pdf

The developed PDF quantification schemes provide statistical description for relative quantities (with respect to the base), that is $x=\hat{x}-\xi$, $y=\hat{y}-\xi$ and $v=\hat{v}-\xi$. However, for the prototype systems that we consider we are more interested for the suppression of absolute quantities, instead of relative ones. As we illustrate below, the absolute response PDF can be derived from the relative response PDF in a straightforward manner .

Background component

For the background regime, we the absolute motion is expressed as:

[TABLE]

As the relative motion and base motion $h(t)$ are both Gaussian distributed (but not independent), their sum is also Gaussian distributed and it is given by,

[TABLE]

In the previous section we have derived both $\sigma_{x}^{2}$ and $\sigma_{h}^{2}$, and what remains is the covariance term $\sigma_{xh}$ whose spectral density function is given in equation 29. This is given by:

[TABLE]

Extreme event component

For the extreme event component the motion of the motion is assumed very small (compared with the magnitude of the impact), in which case we have:

[TABLE]

The estimation of the conditional PDF for $x_{r}$ has already been described in section 4.2.

Comparison with Monte-Carlo Simulations

The full absolute response PDF is expressed using eq. (36), where $z$ is either relative or absolute displacement, velocity or acceleration of the seat/attachment response. We compare the PDS method with direct Monte-Carlo simulations for the case of absolute motions. In figure 7 we show the absolute response PDF for the primary structure for parameters given in table 1. Details regarding the Monte-Carlo computations are provided in section 4.2.

5 System optimization for extreme event mitigation

We now consider the problem of optimization in the presence of stochastic excitation containing extreme events. The developed method provides a rapid and accurate semi-analytic estimation scheme for the statistical response of the nonlinear structural system. In particular, we can efficiently obtain the response statistics of the primary structure (the seat) for any given shock mitigating attachment and accurately capture the heavy-tailed structure of the distribution. This allows us to explore rare event mitigation performance characteristics of different attachment parameters and perform optimization. Such analysis is not practically feasible via a direct Monte-Carlo approach since a single parameter set takes on the order of hours to compute the resulting response PDF with converged tail statistics.

We consider the prototype systems described in section 2 with the aim to suppress the large energy delivered to the passenger (i.e. the seat). In all cases we optimize the attachment parameters, while the parameters of the primary structure are assumed to be fixed.

5.1 Optimization objective

We adopt the forth-order moment as our measure to reflect the severity of extreme events on the seat:

[TABLE]

where the argument $\hat{z}$ can be either absolute displacement of the seat or absolute velocity depending on the optimization objective. The goal here is to minimize this measure and analyze the performance characteristics of the attachment when its parameters are varied.

We illustrate the results of the optimization using the following normalized measure:

[TABLE]

where $\hat{z}_{a}$ is either $\hat{x}$ or $\hat{\dot{x}}$, and $z_{o}$ is the corresponding quantity without any attachment. Values of this measure which are less than 1 ($\gamma<1$) denote effective extreme event suppression.

5.2 Optimization of NES and TMD parameters

Results are shown for the suspended seat problem with an attachment mass $m_{a}=0.05.$ For a NES attachment ($k_{a}=0$) we optimize over $c_{a}$ and $\lambda_{a}$, while for a TMD ($c_{a}=0$) we vary $k_{a}$ and $\lambda_{a}$ (figure 8). The resulted response PDF that minimize the displacement moments are illustrated in figure 9. The same analysis is performed for the suspended deck-seat problem with the same attachment mass $m_{a}=0.05$ for both systems (figure 10). The resulted response PDF are illustrated in figure 11.

In both cases of systems we observe that the TMD and the optimal cubic NES can improve significantly the behavior of the primary structure in terms of reducing the displacement during impacts, with a reduction of 66-68% of the fourth-order moment. We also observe that the NES design is more robust to variations in the attachment parameters over the TMD design, which requires more stringent attachment parameter values for best performance with respect to $\gamma$. This is in line with the fact that the NES attachment performs better over a broader excitation spectrum than the TMD configuration, which requires carefully tuning. Note that for the case of the deck-seat problem (figure 10) we can achieve much larger mitigation of the absolute velocity at the order of 32-34% compared with the simpler system of the seat attached to the hull directly (figure 8), where the suppression is much smaller, 2-4%.

We performed the grid search for demonstration purposes to illustrate the performance characteristics as the stiffness and damping are varied; clearly, if we are only interested in the optimal attachment the use of an appropriate global optimizer (such a particle swarm optimizer) would be more appropriate. All the results shown where computed using the proposed PDF estimation method. As a further check and validation, we benchmarked the semi-analytical PDF estimates and compare them with Monte-Carlo results for the extremity measure $\gamma$ over a coarse grid of the attachment parameters.

6 Design and optimization of a piecewise linear NES

To further improve the shock mitigation properties of the attachment, we utilize a more generic form of NES consisting of a possibly asymmetric, piecewise linear spring. Similarly with the cubic NES and TMD attachments, we perform parameter optimization on the NES spring restoring characteristics and obtain a new optimal design that outperforms the TMD and cubic NES for the considered problems.

Here, we focus on suppressing large displacements of the seat, although velocity or acceleration would also be appropriate depending on the desired objectives. The general form of the considered spring consists of a linear regime with slope equal to that of the optimal TMD within a range of 4 standard deviations of the expected seat motion (e.g. when the TMD is employed). For motions (displacements) outside this range the spring has also a linear structure but with different slopes, $\alpha_{-1}$ for negative displacements (beyond 4 standard deviations) and $\alpha_{1}$ for positive displacements (beyond 4 standard deviations). Therefore, the optimal linear stiffness operates for small to moderate displacement values and outside this regime, when the response is very large, we allow the stiffness characteristics to vary. The objective is to determine the optimal values for the curve in the extreme motion regime with respect to optimization criterion.

Therefore, the analytical form of the piecewise linear spring is given by:

[TABLE]

where, $\sigma_{\zeta}$ is the standard deviation of the relative displacement $\zeta=x-v$ between the primary structure (the seat) and the attachment for the case of a TMD attachment. The parameters, $\alpha_{1}\geq 0$ and $\alpha_{-1}\geq 0$ define the slopes in the positive and negative extreme response regimes, which we seek to optimize. Moreover, the values for $\beta_{1}$ and $\beta_{-1}$ are obtained by enforcing continuity:

[TABLE]

The value of the stiffness in the center regime, $k_{o}$, is chosen using the optimal TMD attachment.

6.1 Application to the suspended seat and deck-seat problem and comparisons

We illustrate the optimization using the fourth-order moment of the seat response, employing the following measure:

[TABLE]

where $\hat{z}_{t}$ is the system response with the optimal TMD attachment (from the previous parametric grid search optimization) and $\hat{z}_{n}$ is the response of the system with the piecewise linear NES attachment. Parameters corresponding to values less than 1 ($\gamma^{\prime}<1$) denote additional extreme event suppression, compared with the utilization of optimal TMD.

The result of the optimization for minimum fourth-order moment for the displacement, on the suspended seat problem, is shown in figure 12 while the corresponding PDF for the displacement, velocity and acceleration are shown in figure 13. We note the strongly asymmetric character of the derived piecewise linear spring. This is directly related with the asymmetric character of the impulsive excitation, which is in general positive. The performance of the optimized piecewise linear spring is radically improved compared with the optimal cubic NES and TMD as it is shown in the PDF comparisons. Specifically, for rare events (probability of 1%) we observe a reduction of the motion amplitude by 50%, while for the velocity the reduction is smaller. A representative time series illustrating the performance of the optimal design for the suspended seat problem is shown in figure 14. The PDF for the acceleration for this set of parameters is not changing significantly. Our results are in agreement with previous studies involving single-sided vibro-impact NES that have been shown to improve shock mitigation properties in deterministic setups [34].

The result of the optimization for the suspended deck-seat problem is shown in figure 15 and the corresponding PDF are shown in figure 16. Similarly with the previous problem, the optimization in this case as well leads to a strongly asymmetric piecewise linear spring. The reduction on the amplitude of the displacement during extreme events is radical (with an additional reduction of 32%) while the corresponding effects for the velocity and acceleration are negligible. This small improvement for the velocity is attributed to the fact that we have focused on minimizing the fourth-order moments for the displacement.

7 Summary and conclusions

We have formulated a parsimonious and accurate quantification method for the heavy-tailed response statistics of nonlinear multi-degree-of-freedom systems under extreme forcing events. The computational core of our approach is the probabilistic decomposition-synthesis method which is formulated for nonlinear MDOF systems under stochastic excitations containing extreme events. Specifically, the excitation is modeled as a superposition of a Poisson distributed impulse train (with extreme magnitude and large inter-arrival times) and a background (smooth) component, modeled by a correlated stochastic excitation with broadband spectral density. This algorithm takes the form of a semi-analytical formula for the response PDF, allowing us to evaluate response statistics (having complex tail structure) on the order of seconds for the nonlinear dynamical structures considered.

Based on this computational statistical framework, we proceed with the design and optimization of small attachments that can optimally mitigate and suppress the extreme forcing events delivered to the primary system. We performed the suppression of extreme responses on prototype ocean engineering dynamical structures, the suspended seat and the suspended deck-seat of high speed crafts, via optimal TMD and cubic NES attachments through parametric optimization. As an optimization criterion we selected the forth-order moments of the response displacement, which is a measure of the severity of large deviations from the mean. Quantitative comparisons of TMD and cubic NES were presented, evaluating the effectiveness and robustness in terms of extreme event suppression. We then proposed a new piecewise linear NES with asymmetries, for extreme event mitigation. The optimization of the new design led to a strongly asymmetric spring that far outperforms the optimal cubic NES and TMD for the considered problem.

We emphasize the statistical accuracy of the PDF estimation schemes, which we demonstrated through comparisons with direct Monte-Carlo simulations. The presented schemes are generic, easy to implement, and can profitably be applied to a variety of different problems in structural engineering where similar characteristics are present, i.e. structures excited by extreme forcing events represented by impulsive-like terms that emerge from an otherwise random excitation background of moderate magnitude.

Acknowledgments

We would like to acknowledge Dr. Timothy Coats, NSWC Carderock, for discussions on the effects of high-speed planing craft wave slams on human performance. T.P.S. has been supported through the ONR grants N00014-14-1-0520 and N00014-15-1-2381 and the AFOSR grant FA9550-16-1-0231. H.K.J. and M.A.M. have been supported through the first and third grants as graduate students. We are also grateful to the Samsung Scholarship Program for support of H.K.J. as well as the MIT Energy Initiative for support under the grant ‘Nonlinear Energy Harvesting From Broad-Band Vibrational Sources By Mimicking Turbulent Energy Transfer Mechanisms’.

Appendix A Statistical linearization of the background regime

Suspended seat system with a linear attachment

For the special case of a linear attachment, $c_{a}=0,$ the operators for the suspended seat problem, $\mathcal{A},$ $\mathcal{B}$, and $\mathcal{C}$ in equations 32, 33, and 34 reduce to

[TABLE]

In this case, we can directly integrate equations 27, 28, and 29 to obtain the second order response statistics.

Suspended deck-seat system

For the suspended deck-seat design the background response is governed by the following system

[TABLE]

As before we first multiply the above two equations by $y(s)$, $x(s)$, $v(s)$, $h(s)$ at different time instant $s\neq t$, and take ensemble averages to write the resulting equations in terms of covariance functions.

[TABLE]

where $\eta$ can be either $y$, $x$, $v$, or $h$, and ′ indicates the partial differentiation with respect to the time difference $\tau=t-s$. We then apply Isserlis’ theorem based on the Gaussian process approximation for response to express the fourth-order moments in terms of second-order moments [7].

[TABLE]

This leads to a set of linear equations in terms of covariance functions and thus the Wiener-Khinchin theorem can be applied to write the equations in terms of the power spectrum. The spectral equations in this case are given by

[TABLE]

where

[TABLE]

Now $\sigma_{y}^{2}$, $\sigma_{v}^{2}$, and $\sigma_{yv}$ are still unknown, but can be determined by integrating both sides of equations 64, 66, 68, and 70 and forming the following system of equations,

[TABLE]

from which we obtain $\sigma_{y}^{2},\sigma_{yv},\sigma_{v}^{2}$.

Suspended deck-saet system with a linear attachment

If the attachment is linear $c_{a}=0$, $\mathcal{A},$ $\mathcal{B}$, and $\mathcal{C}$ in equations 78, 75, and 74 reduce to

[TABLE]

which can be directly integrated to obtain the second order response statistics.

Appendix B Monte-Carlo simulations

For the Monte-Carlo simulations the excitation time series is generated by superimposing the background and rare event components. The background excitation, described by a stationary stochastic process with a Pierson-Moskowitz spectrum (equation 4), is simulated through a superposition of cosines over a range of frequencies with corresponding amplitudes and uniformly distributed random phases. The intermittent component is the random impulse train, and each impact is introduced as a velocity jump at the point of the impulse. For each of the comparisons performed in this work we generated $10$ realizations of the excitation time series, each with a train of $100$ impulses. Once each ensemble for the excitation is computed, the governing ordinary differential equations are solved using a 4th/5th order Runge-Kutta method (we carefully account for the modifications in the momentum that an impulse imparts by integrating up to each impulse time and modifying the initial conditions that the impulse imparts before integrating the system to the next impulse time). We verified that this number of ensembles and their durations leads to converged response statistics for the displacement, velocity, and acceleration.