Investigation of the Energy Dissipation Characteristics of a Nonlinear Suspension System

1. Introduction

To meet the growing needs of modern society, innovative modes of transportation are being developed to improve ride quality, safety and comfort. These improvements are also expected to boost proficiency and ensure cost-effectiveness for all modes of transportation. All ground vehicles experience different levels of vibration, due to irregularities of the road surface. These vibrations are transmitted to the body of the vehicle (including the passengers), through the tyres of the vehicle and the suspension systems (shock absorbers) of the vehicle. Such vibrations can significantly reduce ride comfort, quality and also affect ride stability [1,2,3,4]. In [2], the control of roll and vertical vibrations, of a seat suspension (a two-layer multi-DOF active seat suspension), triggered by the irregular road levels, underneath both sides of the tyres, was investigated. Two controllers, a non-singular terminal sliding mode controller, and an H_∞ controller, were employed. The effectiveness of both the seat suspension system and the control techniques employed, were validated using both simulation and experimental analyses. The authors in [3], investigated the impact of asymmetrical damping and geometrical nonlinearities, in the suspension system of a vehicle model, towards ride comfort. It was demonstrated, both numerically and experimentally, that a combination of asymmetric and geometric nonlinearities, improved the ride comfort of passengers, as well as the performance of the suspension system. In [4], an in-wheel vibration absorber, comprising a spring, an annular rubber bushing, and an adjustable damper, whose parameters were determined using an improved particle swarm optimisation (IPSO) algorithm, was proposed for attenuating vibrations generated by the motor. Two control algorithms, a linear quadratic regulator (LQR) and a fuzzy proportional-integral-derivative (PID) approach, were employed to facilitate the control of the suspension damper and the in-wheel damper, respectively. Both methods displayed excellent performance.

Vehicle suspension systems are designed to curb the effect of adverse ground vibrations and typically comprises a spring, damper and frames/linkages [5]. However, the mass constitutes the body of the vehicle and passengers [5]. A major significance of an isolation spring is its ability to absorb and store energy, from motion bumps and acceleration, with ease. However, isolation springs such as rubber, coil and leaf springs are not capable of dissipating the absorbed energy easily. This can result in unacceptable spring performance before eventual failure, which makes the need for an adequately designed damping system paramount [6]. The damper plays a vital role in the isolation of vibrations from the vehicle’s mass, as it damps out vertical motions by dissipating the transmitted vibration energy as heat energy [6,7]. This ensures the life cycle of the spring is extended and the structural integrity of the system of interest is maintained. The design of a good isolation system is important in vehicle suspension systems to ensure a decent ride comfort [7,8,9]. Lu and his co-authors, in [10], classified nonlinear dampers into three groups namely- nonlinear stiffness dampers, nonlinear-stiffness nonlinear-damping dampers, and nonlinear damping dampers. Based on this classification, the authors reviewed three dampers that are commonly employed in real-world engineering - nonlinear energy sink (NES), particle impact damper (PID), and nonlinear viscous damper (NVD), respectively [10].

Figure 1. Applications of the suspension system.

Figure 1. Applications of the suspension system.

Due to the inherent nonlinearity present in all physical shock and vibration isolation systems, more interests have been shown in nonlinear vibration isolation systems. This has led to nonlinear isolation designs of engine mounts being widely reported in literature [11,12,13,14,15].

The authors in [11] compared two models of ride comfort for a multibody dynamic system, based on different damping characteristics - nonlinear damping and its equivalent linear damping. The responses from the driver’s seat were compared using these models and it was revealed that the model based on the nonlinear damping provided a more accurate result compared to that provided by the equivalent damping system. However, in [12], the authors developed a novel nonlinear quasi-zero-stiffness (QZS) vibration isolator, whereby the suspension system was made nonlinear from the geometric arrangement of the inclined springs. The QZS isolator is characterised by a high-static and low-dynamic stiffness. Results obtained showed that the QZS vibration isolator, with its active controller, demonstrated a good isolation performance for large amplitude excitations. The authors in [13] conducted a computational study on a similar model as in [12], which is the QZS isolator. However, for their study, they considered various road profiles as per ISO 8608. The study revealed that the proposed system, reduced about 92% and 88% of external excitations, on an average, in road class B-D and E-H, respectively. The impact of the velocity of a vehicle, with a nonlinear suspension system, on the dynamic behaviour of a Bernoulli-Euler bridge, was investigated in [14]. It was established that the speed of the vehicle, as well as its suspension stiffness and damping characteristics, are the only significant parameters that can influence the dynamic behaviour of a bridge subjected to a moving vehicle. In [15], the influence of asymmetrical damping and geometrical nonlinearities, in the performance of a vehicle’s suspension system, was explored. Numerical and experimental results verified that geometrical nonlinearity significantly impacts the behaviour of the damping and stiffness characteristics of a suspension system, thus, facilitating its ability to minimise vibrations.

Several analytical methods exist for the study of nonlinear suspension systems, such as the harmonic balance methods [16,17], averaging methods [18,19] and describing function methods [20,21]. However, these methods do not clearly reveal the relationship between the vehicle suspension performance metrics and the system’s nonlinear parameters. The frequency-based method, the Output Frequency Response Function (OFRF) method [22,23,24,26], is a systematic frequency domain approach that enables the derivation of the system’s output spectrum, in terms of the nonlinear parameters of interest. This method is employed in the present study, as it facilitates the analysis and design of nonlinear systems, in the frequency domain.

A model for a single degree of freedom (SDOF) passive vibration isolation system, with cubic stiffness and damping characteristics, was studied in [25]. The method of averaging perturbation was used, to observe the effects of nonlinearities, on the dynamic behaviour of the system. The same nonlinear passive mount was also considered by Peng and Lang in [26]. However, the OFRF method was employed in the theoretical analysis of the output frequency response behaviour, with respect to the system nonlinearities and subsequently verified by simulation studies. A similar system was investigated in [27], where the effects of the system nonlinear parameters, on the relative and absolute displacement transmissibility, were explored.

In this study, the vehicle suspension system, shown in Figure 2, and similar to that in the articles [25,26,27], is investigated. However, an analysis and design of the suspension system is conducted, using the OFRF method. Optimal values of the design parameters are obtained for a desired design specification, and the designed parameters are validated using numerical simulations. Furthermore, for the first time, the energy dissipation performance, of a vehicle suspension system, is investigated using the OFRF method. The effects of the nonlinearities of the designed suspension system, on the energy dissipated by the isolation system, is explored. An investigation is then performed using the OFRF method, where the nonlinear system parameters are designed, for a desired energy dissipation level, of the suspension system.

The rest of the paper are organised as follows: Section 2 provides the model formulation of the suspension system to be investigated; Section 3 introduces the OFRF method, while Section 4 describes the determination of the OFRF representation, for the system output spectra; Section 5 presents the system design and optimisation process for a desired design specification, while Section 6 demonstrates, with numerical simulations, the effects of the system parameters; Section 7 analyses the energy dissipation level of the designed suspension system, as well as the effects of system nonlinearities on the dissipation level; Conclusions are reached in Section 8.

2. Model Formulation

A schematic diagram of the nonlinear Single degree of freedom (SDOF) suspension system, considered in this study, is shown in Figure 2. This includes a suspended oscillating sprung mass, $m$, representing the quarter mass of a vehicle body, with an absolute displacement of $x(t)$, due to a base motion of y(t). The suspension system has two nonlinear isolation components comprising nonlinear damping (c₁, c₂) and nonlinear stiffness $(k_1,k_2)$ elements.

Figure 2. Schematic of the nonlinear suspension system under a base excitation.

Figure 2. Schematic of the nonlinear suspension system under a base excitation.

The equation of motion for the suspension system is given by

$$m\ddot{z}+c\dot{z}+kz=-m\ddot{y}$$

(1)

where $z=x-y$ is the relative displacement between the base and the sprung mass, while $c$ and $k$ are the damping and stiffness functions, respectively. However, the damping and stiffness functions are described as

$$\{\begin{array}{l}k=k_1+k_2{z}^2\\ c=c_1+c_2{z}^2\end{array}$$

(2)

Assuming the base excitation is harmonic, that is

$$y=Y\sin (\omega t)$$

(3)

then system (1) becomes

$$\{\begin{array}{l}m\ddot{z}+(c_1+c_2{z}^2)\dot{z}+(k_1+k_2{z}^2)z=m{\omega }^{2}Y\sin (\omega t)\\ f_{\mathrm{trans}}=(c_1+c_2{z}^2)\dot{z}+(k_1+k_2{z}^2)z\end{array}$$

(4)

where $\omega$ and $Y$ are the excitation frequency and amplitude of the base displacement, respectively, while $f_{\mathrm{trans}}$ is the force transmitted through the vibration isolation system (a parallel combination of the spring and damper elements).

To ensure ride comfort, the minimisation of the transmitted force is required. The transmitted force output spectrum is maximum at the resonant frequency,${\omega }_r$. Consequently, ${\omega }_r$ becomes the frequency of interest and the objective will be to minimise $f_{\mathrm{trans}}$, at the resonant frequency. The procedure outlined by Zhu and Lang in [28], for the analysis and design of a SIMO system, using the OFRF method, will be employed here.

In the next section, the OFRF method will be introduced, considering a general Volterra system, described by a nonlinear differential equation.

3. The Output Frequency Response Function (OFRF)

The differential equation presented in Equation (5) is used to describe a class of Volterra systems

$${\displaystyle \sum _{\overline{q}=1}^Q{\displaystyle \sum _{v=0}^{\overline{q}}{\displaystyle \sum _{r_1,\dots r_{\overline{q}}=0}^R{c}_{v,\overline{q}-v}}}}(r_1,\cdots r_{\overline{q}}){\displaystyle \prod _{i=1}^v\frac{d^{r_i}y(t)}{d{t}^{r_i}}}{\displaystyle \prod _{i=v+1}^{\overline{q}}\frac{d^{r_i}u(t)}{d{t}^{r_i}}}=0$$

(5)

where, $Q$ is the maximum degree of nonlinearity, with respect to the input $u(t)$ and output $y(t)$ parameters of the system, and $R$ is the order of the derivative. Conforming to the OFRF method [22,23,24,26], the output frequency response, of the system described in Equation (5), can be characterised by a polynomial function, in terms of the system nonlinear parameters as

$$Y(\mathrm{j}\omega )={\displaystyle \sum _{{\gamma }_1=0}^{z_1}\cdots {\displaystyle \sum _{{\gamma }_{S_N}=0}^{z_{S_N}}{\mathsf{\Upsilon }}_{({\gamma }_1,\dots ,{\gamma }_{S_N})}}}(\mathrm{j}\omega ){\eta }_1^{{\gamma }_1}\dots {\eta }_{S_N}^{{\gamma }_{S_N}}$$

(6)

where, $z_i$ are the maximum order of ${\eta }_i$, for $i=1,\dots ,S_N$, in the output spectrum $Y(\mathrm{j}\omega )$, of the system described by the polynomial in Equation (6). The system of symbols, ${\mathsf{\Upsilon }}_{({\gamma }_1,\dots ,{\gamma }_{S_N})}(\mathrm{j}\omega )$ are frequency functions, which are also referred to as the OFRF coefficients. These functions are complex values and are dependent on the system’s input and linear characteristic parameters, where ${\mathsf{\gamma }}_i=0,\dots ,z_i$ and $i=1,\dots ,S_N$. Furthermore, ${\eta }_1^{{\gamma }_1}\dots {\eta }_{S_N}^{{\gamma }_{S_N}}$ is a set of monomials and also known as the OFRF structure. They are derived in terms of the nonlinear characteristic parameters of the system. It should be noted that ${\mathsf{\Upsilon }}_{({\gamma }_1,\dots ,{\gamma }_{S_N})}(\mathrm{j}\omega )$ represents the coefficients of ${\eta }_1^{{\gamma }_1}\dots {\eta }_{S_N}^{{\gamma }_{S_N}}$, with the same dimension, and needs to be evaluated in order to obtain the OFRF representation of system (5). However, the monomials need to be determined, first, using a recursive algorithm proposed in Peng et al. [26]. For a Single-Input-Multiple-Output (SIMO) system, an elaborate procedure has been described by Zhu and Lang in [28]. Let the set of monomials, in the OFRF representation of the nth-order output spectrum, be denoted as $\mathfrak{M}$ and the frequency function vector, be denoted as $\mathsf{\Phi }(\mathrm{j}\omega )$, then the OFRF can be described as

$$Y(\mathrm{j}\omega )=\mathfrak{M}\cdot \mathsf{\Phi }{(\mathrm{j}\omega )}^{\mathrm{T}}$$

(7)

Where

$$\mathfrak{M}={\displaystyle \underset{n=1}{\overset{N}{\cup }}{M}_n}$$

(8)

Here, the monomials, $M_n$ can be determined using [26]

$$\begin{array}{l}{M}_n=\left[{\displaystyle \underset{r_1,\dots ,r_{\mathrm{n}}=0}{\overset{R}{\cup }}\left[c_{0,n}(r_1,\dots ,r_{\mathrm{n}})\right]}\right]\hspace{0.17em}\cup \left[{\displaystyle \underset{\overline{q}-v=1}{\overset{\mathrm{n}-1}{\cup }}{\displaystyle \underset{v=1}{\overset{\mathrm{n}-(\overline{q}-v)}{\cup }}{\displaystyle \underset{r_1,\dots ,r_{\mathrm{n}}=0}{\overset{R}{\cup }}\left(\left[c_{v,(\overline{q}-v)}(r_1,\dots ,r_{\overline{q}})\right]\otimes M_n{\hspace{0.17em}}_{-(\overline{q}-v),v}\right)}}}\right]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\cup \left[{\displaystyle \underset{v=2}{\overset{\mathrm{n}}{\cup }}{\displaystyle \underset{r_1,\dots ,r_{\mathrm{n}}=0}{\overset{R}{\cup }}\left(\left[c_{v,0}(r_1,\dots ,r_{\overline{q}})\right]\otimes M_n{\hspace{0.17em}}_{,v}\right)}}\right]\end{array}$$

(9)

Where

$${\mathrm{M}}_{\mathrm{n},v}\hspace{0.17em}\hspace{0.17em}={\displaystyle \underset{i=1}{\overset{\mathrm{n}-v+1}{\cup }}\left({\mathrm{M}}_i\otimes {\mathrm{M}}_{\mathrm{n}-i,v-1}\right),}{\mathrm{M}}_{\mathrm{n},1}={\mathrm{M}}_{\mathrm{n}},\hspace{0.17em}\hspace{0.17em}{\mathrm{M}}_1=1$$

(10)

The symbol ‘$\otimes$’ is the Kronecker product.

Then the set of monomials can be obtained as $\mathfrak{M}={\displaystyle \underset{n=1}{\overset{N}{\cup }}{M}_n}$

In Section 4, the OFRF method will be used to derive an analytical relationship, between the output spectra of system (4) and its nonlinear parameters of interest. Thereafter, a system analysis and design will be conducted, in order to obtain the optimal nonlinear parameters for the desired performance metrics, based on design specifications. The effects of the nonlinear parameters, $c_2$and$k_2$, on the system output spectra, are also explored using numerical simulations.

4. Determination of the OFRF Representation

In this section, the OFRF representation of the output spectra of system (4) is obtained in terms of the nonlinear parameters of interest, $c_2$and$k_2$. As indicated in Equation (6), the OFRF of the nonlinear system (5), is a polynomial function of the nonlinear characteristic parameters of the system. The coefficients of the OFRF are the functions of the output frequency of the system, determined by the specific structure of Equation (5), the system input spectrum, and the linear characteristic parameters of the system. To enable the use of the OFRF representation of system (4), for system analysis and design, its coefficients and monomials have to be determined first. The following analyses have been performed, using the system parameters presented in Table 1.

It can be inferred that the system in Equation (4), is the instance of system (5), with $c_{0,1}(0)=-m{\omega }^{2}Y$, $c_{1,0}(2)=m$, $c_{1,0}(0)=k_1$, $c_{1,0}(1)=c_1$, $c_{3,0}(0,0,1)=c_2$, $c_{3,0}(0,0,0)=k_2$, else $c_{p,q}(·)=0$.

4.1. Determination of the OFRF Monomials

To determine the OFRF representation, the algorithm described in [26], as well as in [28], for SIMO systems, is employed. The recursive algorithm therein, is applied to system (4), to determine the set of monomials, $\mathfrak{M}$, up to the N = 7th order, which yields

$$\mathfrak{M}={\displaystyle \underset{\mathrm{n}=1}{\overset{7}{\cup }}{\mathrm{M}}_{\mathrm{n}}=\left[1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}_2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{k}_2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}_2^2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}_2{k}_2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{k}_2^2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}_2^3\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}_2^2{k}_2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}_2{k}_2^2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{k}_2^3\right]}$$

(11)

It should be noted, that the same set of monomials is obtained for both outputs of system (4). This implies that the OFRF of the outputs of system (4), can be expressed as

$$\{\begin{array}{l}Z(\mathrm{j}\omega )=\mathfrak{M}\cdot \underset{\_}{\mathsf{\Phi }}(\mathrm{j}\omega )\\ F_{\mathrm{trans}}(\mathrm{j}\omega )=\mathfrak{M}\cdot \underset{˜}{\mathsf{\Phi }}(\mathrm{j}\omega )\end{array}$$

(12)

where $\underset{\_}{\mathsf{\Phi }}(\mathrm{j}\omega )$and $\underset{˜}{\mathsf{\Phi }}(\mathrm{j}\omega )$are the OFRF coefficients of the output spectra of system (4), $Z(\mathrm{j}\omega )$and $F_{\mathrm{trans}}(\mathrm{j}\omega )$ respectively. Next, the OFRF coefficients are determined.

4.2. Determination of the OFRF Coefficients

To determine the coefficients of the OFRFs for the system outputs, $Z(\mathrm{j}\omega )$ and $F_{\mathrm{trans}}(\mathrm{j}\omega )$, the steps outlined in [22,23,24,26,28] are followed. The least square approach is used to determine $\underset{\_}{\mathsf{\Phi }}(\mathrm{j}\omega )$ and $\underset{˜}{\mathsf{\Phi }}(\mathrm{j}\omega )$ as

$$\{\begin{array}{l}\underset{\_}{\mathsf{\Phi }}{(\mathrm{j}\omega )}^{{\rm T}}={\left({\psi }^{\mathrm{T}}\psi \right)}^{-1}{\psi }^{\mathrm{T}}\cdot \tilde{Z}(\mathrm{j}\omega )\\ \underset{˜}{\mathsf{\Phi }}{(\mathrm{j}\omega )}^{{\rm T}}={\left({\psi }^{\mathrm{T}}\psi \right)}^{-1}{\psi }^{\mathrm{T}}\cdot {\tilde{F}}_{\mathrm{trans}}(\mathrm{j}\omega )\end{array}$$

(13)

where,

$$\begin{array}{l}\underset{\_}{\mathsf{\Phi }}{(\mathrm{j}\omega )}^{{\rm T}}=\left[\begin{array}{c}{\underset{\_}{\varphi }}_{0,0}(\mathrm{j}\omega )\\ {\underset{\_}{\varphi }}_{1,0}(\mathrm{j}\omega )\\ ⋮\\ {\underset{\_}{\varphi }}_{0,3}(\mathrm{j}\omega )\end{array}\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\underset{˜}{\mathsf{\Phi }}{(\mathrm{j}\omega )}^{{\rm T}}=\left[\begin{array}{c}{\underset{˜}{\varphi }}_{0,0}(\mathrm{j}\omega )\\ {\underset{˜}{\varphi }}_{1,0}(\mathrm{j}\omega )\\ ⋮\\ {\underset{˜}{\varphi }}_{0,3}(\mathrm{j}\omega )\end{array}\right],\hspace{0.17em}\hspace{0.17em}\\ \hspace{0.17em}\\ \tilde{Z}(\mathrm{j}\omega )=\left[\begin{array}{c}Z{(\mathrm{j}\omega )}_{(1)}\\ Z{(\mathrm{j}\omega )}_{(2)}\\ ⋮\\ Z{(\mathrm{j}\omega )}_{(P)}\end{array}\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\tilde{F}}_{\mathrm{trans}}(\mathrm{j}\omega )=\left[\begin{array}{c}{F}_{\mathrm{trans}}{(\mathrm{j}\omega )}_{(1)}\\ F_{\mathrm{trans}}{(\mathrm{j}\omega )}_{(2)}\\ ⋮\\ F_{\mathrm{trans}}{(\mathrm{j}\omega )}_{(P)}\end{array}\right]\hspace{0.17em}\hspace{0.17em}\end{array}$$

(14)

and

$$\psi =\left[\begin{array}{cccccccccc}1& c_{2(1)}& k_{2(1)}& c_{2(1)}^2\hspace{0.17em}& c_{2(1)}{k}_{2(1)}& k_{2(1)}^2& c_{2(1)}^3& c_{2(1)}^2{k}_{2(1)}& c_{2(1)}{k}_{2(1)}^2& k_{2(1)}^3\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& ⋮& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ 1& c_{2(P)}& k_{2(P)}& c_{2(P)}^2& c_{2(P)}{k}_{2(P)}& k_{2(P)}^2& k_{2(P)}^2& c_{2(P)}^2{k}_{2(P)}& c_{2(P)}{k}_{2(P)}^2& k_{2(P)}^3\end{array}\right]$$

$Z(\mathrm{j}\omega )$ and $F_{\mathrm{trans}}(\mathrm{j}\omega )$ represents the output spectra of system (4), when $c_{2(i)}=c_2(i)$ and $k_{2(i)}=k_2(i)$, for $i=1,\hspace{0.17em}2,\cdots ,\hspace{0.17em}P$, where $P=25$ for a varied set of $c_2\in [0,\hspace{0.17em}\hspace{0.17em}28]\hspace{0.17em}$ N.s.m⁻³ and $k_2\in [0,\hspace{0.17em}\hspace{0.17em}2]\times {10}^2$ N.m⁻³. Therefore, considering N=7, the OFRF representations of the output spectra of system (4) can be expressed as

$$\{\begin{array}{c}\begin{array}{l}Z(\mathrm{j}\omega )={\displaystyle \sum _{\overline{n}=0}^3{\displaystyle \sum _{\overline{m}=0}^{\overline{n}}{\underset{\_}{\varphi }}_{\overline{m},\overline{n}-\overline{m}}(\mathrm{j}\omega )c_2^{\overline{m}}{k}_2^{\overline{n}-\overline{m}}}}\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}{\underset{\_}{\varphi }}_{0,0}(\mathrm{j}\omega )+\hspace{0.17em}{\underset{\_}{\varphi }}_{1,0}(\mathrm{j}\omega )c_2+{\underset{\_}{\varphi }}_{0,1}(\mathrm{j}\omega )k_2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}{\underset{\_}{\varphi }}_{2,0}(\mathrm{j}\omega )c_2^2\hspace{0.17em}+{\underset{\_}{\varphi }}_{1,1}(\mathrm{j}\omega )c_2{k}_2\hspace{0.17em}+\hspace{0.17em}{\underset{\_}{\varphi }}_{0,2}(\mathrm{j}\omega )k_2^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\underset{\_}{\varphi }}_{3,0}(\mathrm{j}\omega )c_2^3\hspace{0.17em}+{\underset{\_}{\varphi }}_{2,1}(\mathrm{j}\omega )c_2^2{k}_2+{\underset{\_}{\varphi }}_{1,2}(\mathrm{j}\omega )c_2{k}_2^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\underset{\_}{\varphi }}_{0,3}(\mathrm{j}\omega )k_2^3\end{array}\hfill \\ \mathrm{and}\hfill \\ \begin{array}{l}{F}_{\mathrm{trans}}(\mathrm{j}\omega )={\displaystyle \sum _{\overline{n}=0}^3{\displaystyle \sum _{\overline{m}=0}^{\overline{n}}{\underset{˜}{\varphi }}_{\overline{m},\overline{n}-\overline{m}}(\mathrm{j}\omega )c_2^{\overline{m}}{k}_2^{\overline{n}-\overline{m}}}}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\underset{˜}{\varphi }}_{0,0}(\mathrm{j}\omega )+\hspace{0.17em}{\underset{˜}{\varphi }}_{1,0}(\mathrm{j}\omega )c_2+{\underset{˜}{\varphi }}_{0,1}(\mathrm{j}\omega )k_2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}{\underset{˜}{\varphi }}_{2,0}(\mathrm{j}\omega )c_2^2\hspace{0.17em}+{\underset{˜}{\varphi }}_{1,1}(\mathrm{j}\omega )c_2{k}_2\hspace{0.17em}+{\underset{˜}{\varphi }}_{0,2}(\mathrm{j}\omega )k_2^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\underset{˜}{\varphi }}_{3,0}(\mathrm{j}\omega )c_2^3\hspace{0.17em}+{\underset{˜}{\varphi }}_{2,1}(\mathrm{j}\omega )c_2^2{k}_2\hspace{0.17em}+{\underset{˜}{\varphi }}_{1,2}(\mathrm{j}\omega )c_2{k}_2^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\underset{˜}{\varphi }}_{0,3}(\mathrm{j}\omega )k_2^3\end{array}\hfill \end{array}$$

(15)

Given the output spectra, in terms of the nonlinear parameters of interest, as in Equation (15), parameter optimisation can be easily conducted. For example, for the resonant frequency, which is the frequency of interest, and given as ${\omega }_r=50\hspace{0.17em}\mathrm{rad}{.\mathrm{s}}^{-1}$, the OFRF representations of the system output spectra, at the frequency of interest, are determined as

$$\begin{array}{l}Z(\mathrm{j}{\omega }_r)=(-0.1585-0.2344i)+\hspace{0.17em}(0.0022+0.0040i)c_2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+(2.6195\times {10}^{-4}-1.6383\times {10}^{-4}i)k_2+\hspace{0.17em}(-3.7058\times {10}^{-5}-1.2109\times {10}^{-4}i)c_2^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+(-1.2745\times {10}^{-5}+5.3623\times {10}^{-6}i)c_2{k}_2+(7.1624\times {10}^{-8}+2.6437\times {10}^{-7}i)k_2^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+(1.1426\times {10}^{-7}+1.6910\times {10}^{-6}i)c_2^3+(2.3307\times {10}^{-7}-4.9742\times {10}^{-8}i)c_2^2{k}_2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+(3.8650\times {10}^{-10}-9.5687\times {10}^{-9}i)c_2{k}_2^2+(-1.9915\times {10}^{-10}+1.2264\times {10}^{-11}i)k_2^3\end{array}$$

(16)

and

$$\begin{array}{l}{F}_{\mathrm{trans}}(\mathrm{j}{\omega }_r)=(-6.1703\times {10}^1-1.2908\times {10}^{2}i)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}(1.1198+2.0110i)c_2+(0.1309-0.0820i)k_2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}(-0.0186-0.0606i)c_2^2\hspace{0.17em}+(-0.0064+0.0027i)c_2{k}_2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+(3.5977\times {10}^{-5}+1.3211\times {10}^{-4}i)k_2^2\hspace{0.17em}+(5.7607\times {10}^{-5}+8.4602\times {10}^{-4}i)c_2^3\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+(1.1653\times {10}^{-4}-2.4955\times {10}^{-5}i)c_2^2{k}_2\hspace{0.17em}+(1.9106\times {10}^{-7}-4.7826\times {10}^{-6}i)c_2{k}_2^2\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+(-9.9691\times {10}^{-8}+6.3415\times {10}^{-9}i)k_2^3\end{array}$$

(17)

respectively.

In this case, the magnitude of the output spectra, $\left|Z(\mathrm{j}{\omega }_r)\right|$ and $\left|F_{\mathrm{trans}}(\mathrm{j}{\omega }_r)\right|$, of system (4), are plotted for a varied set of $c_2\in [0,\hspace{0.17em}\hspace{0.17em}28]\hspace{0.17em}$ N.s.m⁻³ while keeping $k_2=240\mathrm{N}{.\mathrm{m}}^{-3}$ fixed, as presented in Figure 3.

The results, as presented in Figure 3, are obtained using the OFRF models, of the output spectra, as given in Equations (16) and (17) respectively, and compared with results obtained from numerical simulations. Note that some of these results have been achieved for parameters beyond those used to determine the OFRF representations. The outcome shows a good match, between the OFRF-based results and results obtained from numerical simulations.

This indicates the usefulness of the OFRF method and, as revealed, it also shows that it provides a good representation of the actual system output spectra. The implication of this is that the OFRF representations, of the output spectra of system (4), can be used for the analysis and design of the system output spectra. This is because, in using the OFRF, an explicit association of the system performance metrics, with the system design parameters, can be established.

It should be noted that a poor OFRF representation of the actual system output spectra, will provide poor analytical results. Such poor representations can be as a result of large numerical errors generated during the determination of the OFRF coefficients. This typically occurs when the parameters are chosen to be large values, with large numerical difference among them, and also when the set of parameters selected, cover a wide range. This will cause the values of the elements in the matrix $\mathfrak{M}$, to be extremely large. Computational errors may result when the inverse of matrix $\mathfrak{M}$ is computed in MATLAB. Such a problem can be solved by following the steps described in [23,28].

5. System Optimisation

The objective now is to minimise the unwanted force transmitted at the resonant frequency, $\left|F_{\mathrm{trans}}(\mathrm{j}{\omega }_r)\right|$ through the isolation elements (spring-damper system), to an acceptable limit. Let $\left|{\mathbb{F}}_{\mathrm{trans}}(\mathrm{j}{\omega }_r)\right|$ be the acceptable threshold of the force transmitted through the vibration isolator, then the optimisation problem can be formulated as

$$\begin{array}{l}\underset{\left[c_2,k_2\right]}{\min }\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left|F_{\mathrm{trans}}(\mathrm{j}{\omega }_r)-{\mathbb{F}}_{\mathrm{trans}}(\mathrm{j}{\omega }_r)\right|\\ \mathrm{subject}\mathrm{to}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\{\begin{array}{l}{c}_2-35\le 0\\ k_2-2.4\times {10}^2\le 0\end{array}\end{array}$$

(18)

Using the OFRF of Equation (17), the relationship between the design parameters, $c_2$ and $k_2$, and the output spectrum, $F_{\mathrm{trans}}(\mathrm{j}{\omega }_r)$, is established and provided in Figure 4 and Figure 5. Looking at Figure 5, it is apparent that the contour map can be used as a design guide. The results of two design specifications are presented in Table 2. To determine the actual output response, $F_{\mathrm{trans}}(\mathrm{j}{\omega }_r)$, the designed parameters are substituted into the system model (4), and simulated, numerically, using the MATLAB Runge-Kutta 4 algorithm.

It should be noted that the second design specification was achieved, for parameter values beyond the range of parameters used to determine the OFRF, i.e., $c_2\in [0,\hspace{0.17em}\hspace{0.17em}28]$ N.s.m⁻³ and $k_2\in [0,\hspace{0.17em}\hspace{0.17em}2]\times {10}^2$ N.m⁻³. In Table 2, it is observed that the actual output responses, obtained using the designed parameters, $c_2$ and $k_2$, clearly matches the specified output responses, with insignificant percentage errors. This reaffirms the effectiveness of the OFRF-based optimal design.

Looking at the time histories of the force transmitted at the resonant frequency in Figure 6, it is obvious that the designed nonlinear vibration isolation system, performs better compared to its corresponding linear counterpart (where the nonlinear characteristic parameters, $c_2=k_2=0$). Moreover, a close observation of the time-history of the transmitted force, due to the designed nonlinear isolation system, reveals that the desired transmitted force, at the resonant frequency, is achieved. It should be noted that minimising the force transmitted along the isolation elements, to an acceptable threshold,  $\left|F_{\mathrm{trans}}(\mathrm{j}{\omega }_r)\right|$, results in a decrease in the system output spectrum, $\left|Z(\mathrm{j}{\omega }_r)\right|$.

Furthermore, though the system design has been performed at the resonant frequency, it is imperative to verify that the designed nonlinear optimal parameters, maintain the same performance on the system output spectrum, across the entire frequency range. This is verified in Figure 7 which compares the performance of the nonlinear design, to its corresponding linear counterpart (i.e., when $c_2=k_2=0$), on the system output spectrum. It is apparent that not only do the designed optimal parameters achieve a specified output spectrum, at the resonant frequency, they are also adequate at suppressing the vibration levels over other frequencies.

The above system analysis and design have been performed using the OFRF method. It is evident that this method facilitates a systematic analysis and design of nonlinear Volterra systems, compared to other methods. The next section shows the individual effect of the system parameters, on the output spectra of the system, considered in this study. Each parameter is varied while fixing others. These have been performed using numerical simulations.

6. Numerical Studies

To demonstrate the effects of $c_1$, $c_2$, and $k_2$, on the output spectra of system (4), numerical studies were performed, and the corresponding results are presented in Figure 8, Figure 9 and Figure 10, respectively. Figure 8 shows the well-known influence of linear damping characteristic parameter, on the output spectra of system (4), where a good isolation performance is witnessed over the resonant region. However, its detrimental effect, over the region where isolation is required ($\mathsf{\Omega }\gg 1$), when high damping level is employed, suggests the need for the integration of a nonlinear damping element. In Figure 9, the effect of $k_2$ is evident around the resonant region. It is apparent that as $k_2$ increases, the resonant frequency, ${\omega }_r$, is increased beyond the natural frequency; however, its harmful effect (instability) is controlled by an increased damping characteristic. Besides, the effect of increasing the nonlinear viscous damping characteristic parameter, $c_2$, is observed in Figure 10, to suppress the vibration at the resonant region. Meanwhile, the vibration response, at the non-resonant regions (i.e., at low and high frequency ranges), remain unaffected.

7. Energy Dissipation Analysis

For the SDOF system described in Figure 2, the vibration isolation elements (spring and damper), experience energy loss due to their dissipative characteristics. The power dissipated at any instant is $f_{\mathrm{trans}}\dot{z}$, and the energy lost over one complete cycle, $E_{\mathrm{d}}$, is given as

$$\begin{array}{l}{E}_{\mathrm{d}}={\displaystyle {\int }_0^{\frac{2\pi }{\omega }}{f}_{\mathrm{trans}}\dot{z}dt}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}={\displaystyle {\int }_0^{\frac{2\pi }{\omega }}\left[(c_1\dot{z}+c_2{z}^2\dot{z}+k_{1}z+k_2{z}^3)\cdot \dot{z}\right]dt}\end{array}$$

(19)

Let the relative displacement, $z=Z\sin (\omega t)$, then $\dot{z}=\omega Z\cos (\omega t)$. Equation (19) becomes

$$E_{\mathrm{d}}={\displaystyle {\int }_0^{\frac{2\pi }{\omega }}\left[\left(\begin{array}{l}{c}_1\omega Z\cos (\omega t)+c_2\omega Z^3{\sin }^2(\omega t)\cos (\omega t)\\ +k_{1}Z\sin (\omega t)+k_2{Z}^3{\sin }^3(\omega t)\end{array}\right)\cdot \omega Z\cos (\omega t)\right]dt}$$

(20)

which yields

$$\begin{array}{l}{E}_{\mathrm{d}}=\pi c_1\omega Z^2+{\displaystyle \frac{1}{4}}\pi c_2\omega Z^4\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\pi \omega Z^2(c_1+{\displaystyle \frac{1}{4}}{c}_2{Z}^2)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\pi \omega Z^2(c_1+c_e)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\pi c_{eq}\omega Z^2\end{array}$$

(21)

where $c_{eq}=c_1+c_e$ and $c_e={\displaystyle \frac{1}{4}}{c}_2{Z}^2$ It is observed, from Equation (21), that the energy dissipated by the suspension system is, mainly, due to its damping characteristics. However, it should be recalled, as seen in Equation (16), that $Z$ is a function of both nonlinear damping, $c_2$ and nonlinear stiffness, $k_2$. The OFRF of the energy lost by the suspension system, over one complete cycle, can be derived by substituting Equation (15a) into Equation (21), which yields

$$E_{\mathrm{d}}(\omega ,c_2,k_2)=\pi c_1\omega {\left|Z(\mathrm{j}\omega )\right|}^2+{\displaystyle \frac{1}{4}}\pi c_2\omega {\left|Z(\mathrm{j}\omega )\right|}^4$$

(22)

Equation (22) indicates $E_{\mathrm{d}}$ is a function of $c_1$, $\omega$, $Z$, and the design parameters, $c_2$ and $k_2$. This indicates that the OFRF, of the energy dissipated, can be derived in terms of the design parameters.

Just as the OFRF representation of the output spectrum, $Z(\mathrm{j}\omega )$, of system (4), was determined and presented in Equation (15a), the OFRFs of the squared and quartic magnitudes of $Z(\mathrm{j}\omega )$, are also derived as

$$\{\begin{array}{l}{\left|Z(\mathrm{j}\omega )\right|}^2={\displaystyle \sum _{\overline{n}=0}^{\lfloor (N-1)/2\rfloor }{\displaystyle \sum _{\overline{m}=0}^{\overline{n}}{P}_{\overline{m},\overline{n}-\overline{m}}(\omega )c_2^{\overline{m}}{k}_2^{\overline{n}-\overline{m}}}}\hspace{0.17em},\\ \mathrm{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ {\left|Z(\mathrm{j}\omega )\right|}^4={\displaystyle \sum _{\overline{n}=0}^{\lfloor (N-1)/2\rfloor }{\displaystyle \sum _{\overline{m}=0}^{\overline{n}}{Q}_{\overline{m},\overline{n}-\overline{m}}(\omega )c_2^{\overline{m}}{k}_2^{\overline{n}-\overline{m}}}}\end{array}$$

(23)

where N = 7, while $P_{\overline{m},\overline{n}-\overline{m}}$ and $Q_{\overline{m},\overline{n}-\overline{m}}$ are functions of frequency, and represent the OFRF coefficients of ${\left|Z(\mathrm{j}\omega )\right|}^2$ and ${\left|Z(\mathrm{j}\omega )\right|}^4$, respectively. Substituting Equation (23) into Equation (22) yields

$$\begin{array}{l}{E}_{\mathrm{d}}(\omega ,c_2,k_2)=\left(\pi c_1\omega \right)\cdot {\displaystyle \sum _{\overline{n}=0}^{\lfloor (N-1)/2\rfloor }{\displaystyle \sum _{\overline{m}=0}^{\overline{n}}{P}_{\overline{m},\overline{n}-\overline{m}}(\omega )c_2^{\overline{m}}{k}_2^{\overline{n}-\overline{m}}}}\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\left({\displaystyle \frac{1}{4}}\pi \omega \right)\cdot {\displaystyle \sum _{\overline{n}=0}^{\lfloor (N-1)/2\rfloor }{\displaystyle \sum _{\overline{m}=0}^{\overline{n}}{Q}_{\overline{m},\overline{n}-\overline{m}}(\omega )c_2^{\overline{m}+1}{k}_2^{\overline{n}-\overline{m}}}}\end{array}$$

(24)

Equation (24) indicates that an estimate of the energy dissipated, per cycle, by the suspension system, can be determined for a set of design parameters. The relationship between the energy dissipation and the design parameters, at the resonant frequency, $\mathsf{\Omega }=1$, is provided in Figure 11. It is seen, in Figure 12, that a large amount of energy is dissipated when a linear vibration isolator is used (when $c_2=0,k_2=0$). However, the energy dissipated reduces with the integration of nonlinear damping and stiffness characteristics. This behaviour is due to the impact of the nonlinear components, on the relative displacement of the suspension system. The relative displacement decreases as the nonlinearities increase, which, consequently, reduces the force transmitted and energy dissipated, by the suspension system. It is also evident that the energy dissipation level is very sensitive to the nonlinear damping characteristic, and less sensitive to the nonlinear hardening stiffness characteristic. Surface plots between the energy dissipated by the suspension system and the design parameters, beyond the resonant region, that is, at $\mathsf{\Omega }=0.4$ and $\mathsf{\Omega }=2$, are presented in Figure 13 and Figure 14, respectively.

Besides, Equation (24), which is the OFRF of the energy dissipation of the suspension system, can be used to estimate the energy dissipated, for a set of optimal design parameters, $c_2$ and $k_2$. An estimation guide can be developed using an OFRF-generated contour map which can facilitate the selection of design parameters for a specified energy dissipation level, as illustrated in Figure 15.

Choosing $k_2=200$ N.m⁻³, along the contour line of 20J, gives a corresponding nonlinear damping value of $c_2=27.5$ N.s.m⁻³. Substituting $c_2$ and $k_2$ into system (4) and solving for the energy dissipation level, numerically, yields a value of 21.08J, which matches the design requirement well, with a percentage error of -5.1%.

The results from numerical studies, to show the individual effects of the system parameters, $c_1$, $c_2$ and $k_2$, on the energy dissipation level of the suspension system, over the entire spectrum, are presented in Figure 16, Figure 17 and Figure 18, respectively. This was done by varying either one of the parameters, $c_1$, $c_2$ or $k_2$, while keeping the other two fixed. In Figure 16, the unwanted impact of linear damping, over the high frequency region, experienced by the output spectra of system (4), is also apparent for the energy dissipation spectra. It is seen that while an increase in linear damping reduces the energy dissipation, over the resonant region, it increases the energy dissipation, significantly, across the non-resonant regions. Likewise, increasing the nonlinear damping, minimises the energy dissipation, at the resonant region. However, unlike the linear damping, the corresponding increase in the energy dissipation, in the non-resonant regions, is far less significant, as illustrated in Figure 17.

The influence of hardening stiffness, as shown in Figure 18, is similar to the effect it has on the output spectra of system (4). In this case, as nonlinear stiffness is increased, the resonant peak of the energy dissipation level, is reduced at $\mathsf{\Omega }=1$. This is due to the resonant-shift effect of $k_2$, on the resonant energy, beyond the natural frequency.

Using the OFRF representations derived, the variations of both the energy dissipated and force transmitted by the suspension system, as a function of the nonlinear stiffness parameter, $k_2$, are shown in Figure 19. It is observed that while the force transmitted by the suspension system, is less sensitive to $k_2$, the sensitivity of the energy dissipated, to $k_2$, is insignificant. This reaffirms the deduction made from Figure 12. This deduction is largely because the spring element stores potential energy when compressed/stretched, and dissipates insignificant amounts of energy. It is also seen that for a fixed $c_2$, increasing $k_2$, causes an insignificant decrease in the energy dissipated per cycle. However, the vibration force transmitted is increased. Similarly, the results illustrated in Figure 20, are the variations of both the energy dissipated and force transmitted, by the suspension system, as a function of the nonlinear damping parameter, $c_2$. In this case, it is apparent that both performance metrics are significantly sensitive to $c_2$. Both performance metrics also, progressively, decrease as $c_2$ increases. This is because the relative displacement of the vehicle suspension system, decreases as $c_2$ increases, and this causes a decrease in both the force transmitted, and energy dissipated, by the suspension system.

Illustrated in Figure 21, is the variation of the suspension damper force, $F_{\mathrm{d}}\hspace{0.17em}(\mathrm{N})$, and the energy dissipated by the suspension system, $E_{\mathrm{d}}\hspace{0.17em}(\mathrm{J})$, with increase in excitation frequency. It is observed that both quantities attain maximum values at the resonant frequency $(\mathsf{\Omega }=1)$. This is expected as a high damping force is required at resonance (i.e., region of maximum displacement), which corresponds to a high energy dissipation level, by the suspension system. However, the energy dissipation, beyond resonance, is significantly low, with respect to the level of damping force, compared to other frequency regions. This is due to the low relative displacement experienced at high frequencies.

Four hysteretic curves (A, B, C, and D), are presented in Figure 22, from different combinations of the designed nonlinear parameters of system (4). For a vibration isolation system, the area of the Force-displacement curve ($F_{\mathrm{trans}}-Z$), gives an indication of the energy dissipated per cycle, by the isolation system, at any frequency of interest (typically the resonant frequency). It is observed that the hysteretic curves C and D, show a smaller displacement span, compared to A and B. This is evidently due to the contribution of the nonlinear damping characteristic, integrated in the suspension systems of C and D. Those of A and B possess a linear damping characteristic only. Nonlinear dampers exhibit a superior vibration isolation performance compared to the corresponding linear case (when $c_2=0$). The presence of a nonlinear damping characteristics, results in the decrease in the magnitude of the relative displacement, $Z$.

It should be noted that, compared to linear dampers, nonlinear dampers dissipate, significantly, more amount of heat energy, for an equal amount of excitation [29]. However, such a comparison has not been considered here. Furthermore, while for a linear suspension system configuration, the hysteretic curve assumes an elliptic shape, this is not the case for a nonlinear suspension system.

As observed in Figure 23, the effect of $c_2$, on curves C and D, with intercept $c_1{Z}_{C,D}{\omega }_r=18.75\hspace{0.17em}\mathrm{N}$, on the vertical axis, is apparent, compared to the curves A and B (with intercept, $c_1{Z}_{A,B}{\omega }_r=24.96\mathrm{N}$), without the effect of $c_2$. This is seen to be as a result of the influence of $c_2$, on the system output response, $Z$, which is a function of the nonlinear damping characteristics,$c_2$.

8. Conclusions

In this study, the vibration isolation performance, of a nonlinear vehicle suspension system, modelled as a SDOF vibration isolation system, was investigated. The Output Frequency Response Function (OFRF) method was the analytical tool employed, for the analysis and design of the suspension system, based on the design criteria. This method was used because of its ability to explicitly express the system performance metrics, in terms of the design parameters of interest. The results obtained, from the analytical and simulation studies, demonstrate the effectiveness of this method, in the analysis and design of vibration isolation systems. To demonstrate the extensive effects of the designed nonlinear system parameters, a comparative study was conducted, between the nonlinear system model and the corresponding linear system (where the nonlinear components are fixed to zero). The outcome reveals the influence of the system nonlinearities, especially at the resonant region, as the resonant peak is largely reduced without a detrimental effect at higher frequencies. This demonstrates the superior vibration isolation capability, of nonlinear vibration isolation systems, compared to their linear equivalent.

Furthermore, the energy dissipated per cycle, by the designed suspension system, was also investigated, using the OFRF method. By employing the OFRF method, a relationship between the Energy dissipation level, of the nonlinear suspension system, and the design parameters, was derived. This polynomial relationship facilitated the estimation of an energy dissipation level, for a set of system nonlinear parameters. The results obtained confirmed the advantage of nonlinear suspension systems, over its linear counterparts. Numerical studies were conducted to demonstrate the effects of nonlinear system parameters, on the energy dissipation level of the suspension system.

Finally, the energy dissipation level of the suspension system, was investigated, using hysteretic curves. The effect of the system nonlinear parameters, on the shape of the hysteresis loop, was explored. The sensitivity of the energy dissipation level, to the nonlinear damping characteristic, compared to the nonlinear stiffness characteristic, was evident. The results in this study reaffirm the vibration suppression advantage of nonlinear suspension systems, over the corresponding linear counterpart.