Sage Journals: Discover world-class research

Abstract

A dual-chamber hydro-pneumatic suspension, a modified version of the traditional hydro-pneumatic suspension, can provide enhanced isolation performance. In this article, we use two very different theories, nonlinear vibration theory and fractional calculus theory, to analyze the characteristics of a dual-chamber hydro-pneumatic suspension. We find that both the storage stiffness and damping coefficient show a strong dependency on the excitation frequency. The effect of various suspension parameters on storage stiffness and damping coefficient is discussed in detail using numerical simulation. We also performed experiments that confirm the validity of the mathematical expressions for both the storage stiffness and the damping coefficient.

Keywords

Dual-chamber harmonic balance method fractional calculus equivalent stiffness equivalent damping coefficient

Introduction

For several decades, due to their nonlinear stiffness, pneumatic vibration isolators have been widely used for vibration isolation in many engineering structures such as vehicle suspensions, high-precision manufacturing equipment, and seat suspensions, which were gradually improved.^1–16

A dual-chamber hydro-pneumatic (DCHP) suspension system was reported in El-Tawwab¹ that provides improved suspension performance while avoiding a major increase in complexity and cost. The structure was shown schematically, and a simplified linear mathematical model was employed to describe the suspension dynamics. Simulation results in the time domain show that the twin accumulator suspension leads to significant improvements over the conventional passive suspension system. Erin et al.² point out that a dual-chamber pneumatic isolator would benefit from an accurate mathematical model to improve the design of vibration isolation tables because significant deviations were found between predicted and observed behavior in several studies. The article introduced a modified model for a pneumatic isolator taking into account the effect of the diaphragm. As a result, the accuracy of the predictions using the modified model was better. Sun et al.³ describe the structure of a pneumatic suspension with an auxiliary reservoir, and the test results show the characteristic of varying damping and stiffness. A compact variable stiffness and damping shock absorber was designed, prototyped, tested, and modeled in his later research.⁴ The intergeneration projection genetic algorithm (IP-GA) method was used to identify model parameters. Lee and Kim⁵ adopt a pneumatic spring suspension with a similar structure to the one described in Erin et al.,² which is to be used for vibration isolation for precision manufacturing instruments. A mathematical model capable of high precision considers thermodynamic behavior of air in chambers, fluid-dynamics of the air flow in capillary tubes, and the dynamic characteristic of a diaphragm. It acts as an effective tool to describe the behavior of the DCHP suspension. The complex stiffness derived from a set of motion equations shows a dependency on both the excitaion frequency and the excitation amplitude. The validity of the expression for the complex stiffness was confirmed through experiments. The two terms, storage stiffness and loss factor, are the parameters to be determined for designing a dual-chamber pneumatic isolator. In Lee and Kim’s⁶ later research, the transmissibility from base excitation to the displacement of the spring mass was derived. To minimize the transmissibility at the resonance frequency and to retain good isolation performance for higher frequencies, the effect of the parameters on transmissibility was discussed in detail. Similar studies were described in Porumamilla et al.⁷ and Holtz and van Niekerk⁸ with the aim to reduce the seat effective amplitude transmissibility (SEAT) value used for evaluating the isolation performance of the seat suspension. Other contributions from Nieto and colleages,^9–11 Liu and Lee,¹² Pu et al.,¹³ Li et al.,¹⁴ and Chai et al.¹⁵ are helpful to understand the frequency and amplitude characteristics of a DCHP suspension.

The article consists of six sections. The structure and components are presented in section “Structure of the DCHP suspension.” Nonlinear isolators have attracted increasing attention recently because all isolators are inherently nonlinear. Furthermore, the use of nonlinear stiffness and damping enables superior isolation performance. The harmonic balance method (HBM) is one of the most commonly used tools for nonlinear vibration analysis. Preliminary studies indicate that both the storage stiffness and the coefficient damping are frequency-dependent. Modeling and analysis of the DCHP suspension were carried out using the HBM described in section “Modeling and analysis of the DCHP suspension using nonlinear oscillation theory.” In section “Modeling and analysis of the DCHP suspension based on fractional order calculus theory,” a model of the DCHP suspension was developed based on fractional calculus. We then compared the two modeling methods. In section “Experiment results and discussion,” experiments for varying parameters are described. The measurements confirm the analytical expression for the equivalent storage stiffness and the damping coefficient.

Structure of the DCHP suspension

A schematic of the DCHP suspension is shown in Figure 1. The suspension system consists mainly of five components. A cylinder that supports the spring mass is located between the vehicle body and the front or rear wheel. A piston rod can move axially inside the suspension cylinder. The piston divides the cylinder into two chambers: the compression chamber and the extension chamber. When the piston travels up and down, the oil flows through the valve block. As a result, the main damper can provide a damping force to reduce the oscillations of the vibration system. Two accumulators, which act as elastic elements, can be found in the DCHP suspension. These two accumulators are linked via a variable orifice, that is, an adjustable throttle valve. The valve opening of the adjustable throttle valve can be regulated to adjust the damping coefficient between the accumulators. The effect of additional damping on both storage stiffness and damping coefficient of the isolation system will be discussed in the following sections.

Figure 1.

Schematic of the dual-chamber hydro-pneumatic suspension.

Modeling and analysis of the DCHP suspension using nonlinear oscillation theory

Modeling of the pressure of the rebound chamber

Figure 1 shows the flow of oil in the cylinder. The geometrical change of volume of the rebound chamber due to the moving piston is balanced by the flow of oil through the piston valve and the compression or expansion of oil in this chamber. The damping force is typically calculated using fluid mechanics theory assuming that the used oil is incompressible. Lang¹⁷ considered the compressibility of oil in a mathematic model for the shock absorber. Both simulation and experiment indicate that the compressibility of oil is one of the major causes of the damper’s nonlinear dynamic behavior, such as hysteresis loops that occur at higher excitation frequencies.

The compressibility $β$ of the fluid is defined as¹⁷

β = \frac{d ρ}{ρ dp} = - \frac{dV}{Vdp}

(1)

where V and $ρ$ represent the volume and density of the fluid, respectively, and p is the chamber pressure. Considering the law of conservation, the relationship between mass flow and volumetric flow can be written as¹⁸

Q = \frac{dV}{dt} + β V \frac{dp}{dt}

(2)

The motion of the main piston, x, is enforced by the test rig actuator. As a result, a certain amount of oil, $V_{13}$ , flows into the rebound chamber from the compression chamber due to the pressure difference between the two chambers. Hence, the relationship between the rate of flow and the pressure derivative in the rebound chamber of the DCHP spring can be formulated as

Q_{13} = (A_{piston} - A) \overset{\cdot}{x} + β [V_{30} + (A_{piston} - A) x] \frac{d p_{3}}{dt}

(3)

where $A_{pistion}$ and A represent the section areas of the main piston and the piston rod, respectively. $V_{30}$ is the volume of the rebound chamber when the main piston remains in its initial position. The compressibility parameter $β$ takes into account the oil’s compressibility and the mechanical compliance of the container and seals.

Modeling of the pressure in the accumulators

The compression chamber is connected to the main accumulator (see Figure 1). As a result, the pressure in the compression chamber and the main accumulator is the same. The gas volume change in the main accumulator is caused by the geometrical change of the volume of the compression chamber during the motion of the piston. Apart from the oil flowing into the rebound chamber, there is another volume $V_{12}$ of oil flowing into the additional accumulator. The latter is due to the pressure difference between the two accumulators. The pressure equations for the main and additional accumulators can be derived from the polytropic equation of state, $p V^{r} = constant$

\begin{matrix} p_{1} & = p_{10} {(\frac{V_{10}}{V_{10} - A_{piston} \cdot x + V_{13} + V_{12}})}^{r} \\ p_{2} & = p_{20} {(\frac{V_{20}}{V_{20} - V_{12}})}^{r} \end{matrix}

(4)

Here, $p_{10}$ $(p_{20})$ and $V_{10}$ $(V_{20})$ are the initial pressures and initial volumes of the main accumulator (additional accumulator), and $p_{1}$ $(p_{2})$ is the pressure of the main accumulator (additional accumulator). The adiabatic exponent r, which represents the ratio of the specific heat capacity at constant volume for a particular gas, is about 1.4 for diatomic gases (e.g. nitrogen). However, due to the working condition for the hydro-pneumatic suspension, there is the possibility of a marginal heat exchange between the gas and the environment or surrounding components. As a result, there cannot be a truly polytropic process. The value for the adiabatic exponent is typically between 1 and 1.4. The accurate number can be obtained through experiments.

The oil-compressibility factor for the rebound chamber is often considered during the modeling process. However, for the compression chamber, similar to the mono-tube damper,¹⁹ the oil compressibility is negligible compared with the compressibility of the nitrogen gas enclosed in diaphragm accumulators. As a result, we can simplify as follows

- A_{piston} \cdot x + V_{13} = - A \cdot x

(5)

In other words, the pressure model of the compression chamber can be modified using the simplification assuming that the oil is incompressible and both accumulators have the same initial pressure, $p_{10} = p_{20} = p_{0}$ . Therefore, equation (4) can be simplified as

\begin{matrix} p_{1} & = p_{0} {(\frac{V_{10}}{V_{10} - A \cdot x + V_{12}})}^{r} \\ p_{2} & = p_{0} {(\frac{V_{20}}{V_{20} - V_{12}})}^{r} \end{matrix}

(6)

Modeling of the damping force of the throttle valve

The volumetric flow rate Q for a liquid through an orifice is generally described by the equation $Q = c_{d} A_{d} \sqrt{2 (p_{1} - p_{2}) / ρ}$ , where A_d refers to the cross-section area of the orifice. The discharge coefficient, $c_{d}$ , is often assumed constant. In the absence of any special conditions, $c_{d} = 0.6 ~ 0.7$ can be used.

The liquid damping force in the adjustable throttle can be expressed as the quadratic term of the velocity, plus a linear term

{\begin{matrix} \begin{matrix} p_{1} = p_{2} + \frac{ρ}{2 c_{d}^{2} A_{d}^{2}} {(Q_{12})}^{2} + \frac{c_{L}}{A_{d}^{2}} Q_{12} & Q_{12} \geq 0 \end{matrix} \\ \begin{matrix} p_{1} = p_{2} - \frac{ρ}{2 c_{d}^{2} A_{d}^{2}} {(Q_{12})}^{2} + \frac{c_{L}}{A_{d}^{2}} Q_{12} & Q_{12} < 0 \end{matrix} \end{matrix}

(7)

Simplification of the dynamic characteristic of the DCHP suspension

According to the above-discussed modeling of the components, the transmitted force of the DCHP suspension for a base excitation x is given by

{\begin{matrix} F = p_{0} {(\frac{V_{10}}{V_{10} - Ax + V_{12}})}^{r} A + F_{d} \\ p_{0} {(\frac{V_{10}}{V_{10} - Ax + V_{12}})}^{r} A = \frac{ρ}{2 c_{d}^{2} A_{d}^{2}} | Q_{12} | Q_{12} A + \frac{c_{L}}{A_{d}^{2}} Q_{12} A + p_{0} {(\frac{V_{20}}{V_{20} - V_{12}})}^{r} A \end{matrix}

(8)

where $F_{d}$ is the damping force provided by the main damper. According to the general expression for the Taylor expansion, the pressure $p_{1}$ and $p_{2}$ of the main and additional accumulators can be expanded at $x - z$ and z, respectively. Therefore, equation (8) related to the transmitted force can be rewritten as

{\begin{matrix} F = p_{0} A + \sum_{i = 0}^{n} k_{1 i} {(x - z)}^{i} + F_{d} \\ \sum_{i = 0}^{n} k_{1 i} {(x - z)}^{i} = \sum_{j = 0}^{n} k_{2 j} z^{j} + c_{2 Q} | \overset{\cdot}{z} | \overset{\cdot}{z} + c_{2 L} \cdot \overset{\cdot}{z} \end{matrix}

(9)

Here, $k_{1 i} = p_{0} A^{i + 1} Π_{m = 1}^{i} (r + m - 1) / V_{20}^{i}$ and $k_{2 j} = p_{0} A^{j + 1} Π_{n = 1}^{j} (r + n - 1) / V_{10}^{j}$ $(i = 1, 2, 3, \dots)$ represent the stiffness of the spring force provided by the main and additional accumulator. $z = V_{12} / A,$ $\overset{\cdot}{z} = Q_{12} / A$ , and $δ = x - z .$

For simplicity, the storage stiffness of the accumulator fits the linear model in earlier research. However, the linear model is unable to describe the effect of the nonlinearity of a restoring force on the equivalent stiffness, damping coefficient, and isolation performance for a DCHP suspension. Because it is very difficult to solve equation (9) analytically for high orders of $x - z$ and $δ$ , only the first three items of the restoring force are used in this article, the equation can be rewritten as

{\begin{matrix} F = p_{0} A + k_{11} δ + k_{12} δ^{2} + k_{13} δ^{3} + F_{d} \\ k_{11} δ + k_{12} δ^{2} + k_{13} δ^{3} = k_{21} z + k_{22} z^{2} + k_{23} z^{3} + c_{2 Q} | \overset{\cdot}{z} | \overset{\cdot}{z} + c_{2 L} \cdot \overset{\cdot}{z} \end{matrix}

(10)

The quartic term in the simplified storing force model means that the DCHP suspension is a hardening isolation system in the direction of the compression, but a softening isolation system in the direction of the extension. Overall, it is a hardening system in a complete isolation cycle, which is one of the advanced properties of a hydro-pneumatic suspension.

Solution of the equations for motion for harmonic excitation

Mixed-parity nonlinear oscillation problems, especially isolators with a nonlinear restoring force with an additional quartic term, have been studied previously. Lu et al.²⁰ transform the isolation system with nonlinear restoring and damping forces into a standard duffing isolation system. They derive the displacement transmissibility using the HBM, which is considered suitable for both weak and strong nonlinear vibration problems. The simulation results show that the quartic term for the restoring force affects the static equilibrium position of the oscillation system.

The first-order HBM is used to solve equation (10) using the assumption of the approximate solutions

{\begin{matrix} z = Z_{0} + Z_{1} \cos ω t \\ δ = Δ_{0} + Δ_{1} \cos (ω t + θ) \end{matrix}

(11)

The bias terms $Z_{0}$ and $Δ_{0}$ refer to the static equilibrium positions, and there is a phase difference $θ$ between the displacement solution z and the dynamic deformation $δ .$ Substituting the approximate solutions into equation (10) and equating the coefficients of $\sin (ω t)$ and $\cos (ω t)$ yields the following equations

\begin{matrix} Δ_{0} k_{11} + Δ_{0}^{2} k_{12} + Δ_{0}^{3} k_{13} + \frac{1}{2} Δ_{1}^{2} k_{12} + \frac{3}{2} Δ_{0} Δ_{1}^{2} k_{13} \\ = Z_{0} k_{21} + Z_{0}^{2} k_{22} + Z_{0}^{3} k_{23} + \frac{1}{2} Z_{1}^{2} k_{22} + \frac{3}{2} Z_{0} Z_{1}^{2} k_{23} \\ (Δ_{1} k_{11} + 2 Δ_{0} Δ_{1} k_{12} + 3 Δ_{0}^{2} Δ_{1} k_{13} + \frac{3}{4} Δ_{1}^{3} k_{13}) \cos θ \\ = Z_{1} k_{21} + 2 Z_{0} Z_{1} k_{22} + 3 Z_{0}^{2} Z_{1} k_{23} + \frac{3}{4} Z_{1}^{3} k_{23} \\ (Δ_{1} k_{11} + 2 Δ_{0} Δ_{1} k_{12} + 3 Δ_{0}^{2} Δ_{1} k_{13} + \frac{3}{4} Δ_{1}^{3} k_{13}) \sin θ \\ = \frac{ω 8 c_{2 Q} Z_{1}^{2}}{3 π} + c_{2 L} \end{matrix}

(12)

Squaring and adding the last two equations of equation (12) yield

\begin{matrix} Δ_{0} k_{11} + Δ_{0}^{2} k_{12} + Δ_{0}^{3} k_{13} + \frac{1}{2} Δ_{1}^{2} k_{12} + \frac{3}{2} Δ_{0} Δ_{1}^{2} k_{13} \\ = Z_{0} k_{21} + Z_{0}^{2} k_{22} + Z_{0}^{3} k_{23} + \frac{1}{2} Z_{1}^{2} k_{22} + \frac{3}{2} Z_{0} Z_{1}^{2} k_{23} \\ {(Δ_{1} k_{11} + 2 Δ_{0} Δ_{1} k_{12} + 3 Δ_{0}^{2} Δ_{1} k_{13} + \frac{3}{4} Δ_{1}^{3} k_{13})}^{2} \\ = {(Z_{1} k_{21} + 2 Z_{0} Z_{1} k_{22} + 3 Z_{0}^{2} Z_{1} k_{23} + \frac{3}{4} Z_{1}^{3} k_{23})}^{2} \\ + {(\frac{ω 8 c_{2 Q} Z_{1}^{2}}{3 π} + c_{2 L})}^{2} \\ Δ_{1} = (X_{1} - Z_{1}) \\ Δ_{0} = - Z_{0} \end{matrix}

(13)

where $X_{1}$ is the amplitude of the harmonic excitation, and the bias terms $Z_{0}$ and $Δ_{0}$ have equal values but opposite signs.

For a given sectional area of the piston rod A, the initial pressure $p_{0}$ and the volume $V_{0}$ of the accumulators, the parameters of $k_{11}, k_{12}, k_{13}, k_{21}, k_{22}, k_{23}$ can be calculated. Furthermore, for a given excitation angular frequency $ω$ , the harmonic term $Z_{1}$ and the bias $Z_{0}$ can be solved.

Equivalent stiffness and damping coefficient

According to equation (10), the linear model for the DCHP suspension under base excitation can be given by

{\begin{matrix} F = k_{1} (x - z) + c_{1} \overset{\cdot}{x} \\ k_{1} (x - z) - k_{2} z = c_{2} \overset{\cdot}{z} \end{matrix}

(14)

and the transmitted force using a Laplace transformation of equation (14) is

\begin{array}{l} F = (k_{1} - \frac{{k_{1}}^{2} (k_{1} + k_{2})}{{(k_{1} + k_{2})}^{2} + c_{2}^{2} ω^{2}}) X \\ + j (c_{1} + \frac{{k_{1}}^{2} c_{2}}{{(k_{1} + k_{2})}^{2} + c_{2}^{2} ω^{2}}) ω X \end{array}

(15)

The real part of the transmitted force represents the storage stiffness. The imaginary part is the loss stiffness,¹³ which actually denotes the damping characteristic of the DCHP suspension. The dynamic stiffness and the damping coefficient can be written in the following form

{\begin{matrix} k_{d} = k_{1} - \frac{{k_{1}}^{2} (k_{1} + k_{2})}{{(k_{1} + k_{2})}^{2} + c_{2}^{2} ω^{2}} \\ c_{d} = c_{1} + \frac{{k_{1}}^{2} c_{2}}{{(k_{1} + k_{2})}^{2} + c_{2}^{2} ω^{2}} \end{matrix}

(16)

Although the spring and the damping component are linear, equation (16) shows that both the dynamic stiffness and the damping coefficient have a frequency-dependent characteristic. The dynamic stiffness $k_{d}$ increases with the increase in the excitation frequency $ω$ , while $c_{d}$ shows the opposite trend.

Using the concepts of equivalent stiffness and damping coefficient,^21,22 the equivalent stiffnesses of the main and additional accumulators are

k_{eq 1} = \frac{k_{11} + 2 Δ_{0} k_{12} + 3 Δ_{0}^{2} k_{13} + 3 Δ_{1}^{2} k_{13}}{4}

(17)

k_{eq 2} = \frac{k_{21} + 2 Z_{0} k_{22} + 3 Z_{0}^{2} k_{23} + 3 Z_{1}^{2} k_{23}}{4}

(18)

The adjustable throttle valve with mixed damping has an equivalent damping coefficient

c_{eq 2} = \frac{ω 8 c_{2 Q} Z_{1}}{3 π} + c_{2 L}

(19)

Here, $c_{2 Q}$ and $c_{2 L}$ represent quadratic and linear parts for the mixed damping coefficient, respectively. As seen in equations (17)–(19), both amplitude and frequency dependence represent the essential differences between linear and nonlinear stiffness and damping.

The equivalent stiffness and damping coefficient for the DCHP suspension are given by

\begin{matrix} k_{eq} = k_{11} + 2 Δ_{0} k_{12} + 3 Δ_{0}^{2} k_{13} + \frac{3}{4} Δ_{1}^{2} k_{13} \\ - \frac{{(k_{11} + 2 Δ_{0} k_{12} + 3 Δ_{0}^{2} k_{13} + \frac{3}{4} Δ_{1}^{2} k_{13})}^{2} (k_{11} + 2 Δ_{0} k_{12} + 3 Δ_{0}^{2} k_{13} + \frac{3}{4} Δ_{1}^{2} k_{13} + k_{21} + 2 Z_{0} k_{22} + 3 Z_{0}^{2} k_{23} + \frac{3}{4} Z_{1}^{2} k_{23})}{{(k_{11} + 2 Δ_{0} k_{12} + 3 Δ_{0}^{2} k_{13} + \frac{3}{4} Δ_{1}^{2} k_{13} + k_{21} + 2 Z_{0} k_{22} + 3 Z_{0}^{2} k_{23} + \frac{3}{4} Z_{1}^{2} k_{23})}^{2} + {(\frac{ω 8 c_{2 Q} Z_{1}}{3 π} + c_{2 L})}^{2} ω^{2}} \\ c_{eq} = c_{eq 1} + \frac{{(k_{11} + 2 Δ_{0} k_{12} + 3 Δ_{0}^{2} k_{13} + \frac{3}{4} Δ_{1}^{2} k_{13})}^{2} (\frac{ω 8 c_{2} Z_{1}}{3 π} + c_{2 L})}{{(k_{11} + 2 Δ_{0} k_{12} + 3 Δ_{0}^{2} k_{13} + \frac{3}{4} Δ_{1}^{2} k_{13} + k_{21} + 2 Z_{0} k_{22} + 3 Z_{0}^{2} k_{23} + \frac{3}{4} Z_{1}^{2} k_{23})}^{2} + {(\frac{ω 8 c_{2 Q} Z_{1}}{3 π} + c_{2 L})}^{2} ω^{2}} \end{matrix}

(20)

Unlike the dynamic stiffness and damping coefficient for the DCHP suspension with linear elastic and damping components, both $k_{eq}$ and $c_{eq}$ depend on the excitation frequency $ω,$ the response oscillation amplitude $Δ_{1}$ and $Z_{1}$ , but also on the bias terms $Δ_{0}$ and $Z_{0}$ , because of the mixed-parity of the nonlinear restoring force of the hydro-pneumatic spring. The equivalent damping coefficient, $c_{eq 1},$ provided by the main damper has no effect on the dynamic and equivalent damping coefficient apart from the initial damping value. Therefore, it will not be the focus of this study, and to keep the complexity low, the explicit expression for the equivalent main damping coefficient $c_{eq 1}$ is not provided in equation (20).

The characteristics of the dynamic stiffness and damping coefficient for the linear and nonlinear DCHP suspension in the frequency domain are shown in Figure 2.

Figure 2.

Computational simulation of the (a) storage stiffness and (b) equivalent damping coefficient for different types of damping behavior.

According to Figure 2(a), the storage stiffness curve follows an “S” shape in the frequency domain. The stiffness remains a small constant but increases with increasing frequency to a certain level. When the curve becomes flat, the stiffness reaches its own upper limit. The minimum and maximum storage stiffness of the nonlinear DCHP suspension are deduced by setting $ω = 0$ and $ω = \infty$ , respectively

k_{eqMIN} = \frac{k_{eq 1} \cdot k_{eq 2}}{k_{eq 1} + k_{eq 2}}, k_{eqMAX} = k_{eq 1}

Based on the expressions above, $k_{eqMIN}$ is equal to the stiffness of two elastic elements in series, and $k_{eqMAX}$ is equal to the main spring only. This means that only the main accumulator of the DCHP suspension provides the elastic force for high frequencies.

For the sake of studying the effect of additional nonlinear damping, the maximum as well as the minimum storage stiffness for different types of nonlinear additional damping are set to be equal. Figure 2(a) shows storage stiffness curves indicating similar behavior but not exactly identical. The stiffness curve is much steeper under the influence of the additional damper with hardening flow resection behavior.

The damping curve of the linear model shows a reverse-“S” shape in the frequency domain (see Figure 2(b)). However, affected by the frequency-dependent characteristic, the curve for quadratic damping follows a bell-shape. It is the hardening damping that causes the transcendence for the damping coefficient in the frequency range that $f_{1}$ refers to. The lower limit for the equivalent damping coefficient is $c_{eq 1}$ when the excitation frequency $ω = \infty$ , which means that only the main damper is active.

The additional damping determined by the degree of opening the throttle valve is essential for the isolation performance of the DCHP suspension. Adjusting the throttle valve mainly affects the secondary damping coefficient. Therefore, some simulation results are shown for different values of $c_{2 Q}$ in Figure 3.

Figure 3.

Computational simulation of the (a) storage stiffness and (b) equivalent damping coefficient for different values of $c_{2 Q}$ .

We find that the storage stiffness curves shift along the frequency axis for different $c_{2 Q}$ values, and the shape remains almost invariable. The peak value of the equivalent damping coefficient decreases with decreasing additional damping $(c_{2 Q} > c'_{2 Q} > c_{2 Q}^{″})$ . The equivalent damping coefficient with higher $c_{2 Q}$ has higher values in the high-frequency range, but lower values in the low-frequency range. The exchange points are located at the frequency segment marked by $f'_{1} .$

Modeling and analysis of the DCHP suspension based on fractional order calculus theory

Modeling of the DCHP suspension

This section describes a very different modeling method for the dynamic behavior of the DCHP suspension.

Some mechanical isolators,²³ such as engine mounts and viscoelastic dampers, are made of viscoelastic materials that have the same frequency-dependent characteristic. The properties of viscoelastic materials have been extensively studied in their own material-specific field. Both a theoretical analysis and experiments of mechanical isolators were carried out to study viscoelastic materials. Modeling is done assuming that a certain number of springs and dampers are combined into a mechanical network with the specific series-parallel logic, which is a common method to describe and analyze the behavior of the viscoelastic material.²⁴ The Kelvin–Voigt model consists of a spring and a damper connected in parallel, while the Maxwell model uses them in series. The Zener model consists of a spring connected with a Kelvin model in series, or a spring connected with a Maxwell model in parallel. It is commonly used to describe the dynamic behavior of viscoelastic isolators. Two types of Zener models are shown schematically in Figure 4(a) and (b).

Figure 4.

The schematic diagram of Zener models: (a) Kelvin form and (b) Maxwell form.

Adopting the linear Zener model (Kelvin form) to describe the dual-chamber pneumatic suspension, VM Ryaboy¹⁶ derived the expressions for both the dynamic stiffness and the loss stiffness. Moreover, the transmissibility for a single degree of freedom vibration system installed with the dual-chamber pneumatic suspension from base excitation to spring mass displacement could be established mathematically, together with an analysis of the effect of parameters on the transmissibility. To improve the accuracy of the frequency-dependent behavior described by conventional models, fractional calculus was used. The fractional models have been used as modified versions of these traditional linear models. They were shown to be powerful tools.^24–26 Fractional models have the ability to model the viscoelastic isolators mathematically using only a small number of parameters. More specifically, they use only a single equation instead of the set of differential equations used in traditional models.

T Pritz²⁵ summarized prior studies and put forward a fractional derivative model with only four parameters to describe the variations of the dynamic properties of the viscoelastic materials for a wide frequency range. In later studies,²⁶ a new version of the fractional derivative model characterized by five parameters was introduced. It describes even more accurately the asymmetrical loss factor peak as well as the dynamic behavior for high frequencies.

For viscoelastic dampers, Makris and Constantinou²⁷ proposed a fractional derivative Maxwell model, which was validated by comparing simulation results for dynamic stiffness with experiment observations.

The general form of the standard linear viscoelastic model is given as²⁸

F (t) + \sum_{m = 1}^{M} b_{m} \frac{d^{m} F (t)}{d t^{m}} = E_{0} x (t) + \sum_{n = 1}^{N} E_{n} \frac{d^{n} x (t)}{d t^{n}}

(21)

where $F (t)$ and $x (t)$ have the same definition as above. The symbol E in material research represents the stiffness of the elastic component. The introduction of fractional calculus theory was implemented replacing the integer order time-derivatives with fractional order time-derivatives. The general form of the fractional derivative viscoelastic model can be deduced

F (t) + \sum_{m = 1}^{M} b_{m} \frac{d^{β_{m}} F (t)}{d t^{β_{m}}} = E_{0} x (t) + \sum_{n = 1}^{N} E_{n} \frac{d^{α_{n}} x (t)}{d t^{α_{n}}}

(22)

By keeping only the first-order fractional derivative of $x (t)$ , and ignoring all $β$ terms, equation (22) turns into the fractional Kelvin-Voigt model. By keeping only the first-order fractional derivative of $x (t)$ and $F (t)$ , and ignoring the term $E_{0} x (t),$ equation (22) represents the fractional Maxwell model. The fractional derivative Zener model, which is considered a compromise between complexity and effectiveness, is regarded as the most efficient model to predict dynamic behavior of viscoelastic materials in the frequency domain. It has been verified through experimental observations in previous papers. After simplification of equation (22), the fractional derivative Zener model with five parameters, $b,$ $α,$ $E_{0}$ , $E_{1}$ , and β, is given by

F (t) + b \frac{d^{β}}{d t^{β}} F (t) = E_{0} x (t) + E_{1} \frac{d^{α} x (t)}{d t^{α}}

(23)

The $α th$ order fractional derivation of $x (t)$ can be defined using the gamma function²⁶

\frac{d^{α}}{d t^{α}} x (t) = \frac{1}{Γ (1 - α)} \frac{d}{dt} \int_{0}^{t} \frac{x (t)}{{(t - τ)}^{α}} d τ

(24)

Storage stiffness and equivalent damping coefficient

According to the fractional Fourier transform theory, the Fourier transform of the fractional derivative of order α has the same form as the one of the integer derivative order

F {\frac{d^{α}}{d t^{α}} x (t)} = (i ω)^{α} F {x (t)}

(25)

Here, $F$ denotes the Fourier transform. The transmitted force in the frequency domain is given by

\hat{F} (i ω) = \frac{E_{0} + E_{1} {(i ω)}^{α}}{1 + b {(i ω)}^{β}} \hat{x} (i ω)

(26)

where $\hat{F} (i ω)$ and $\hat{x} (i ω)$ denote the Fourier transforms of $F (t)$ and $x (t) .$ Substituting Euler’s formula into equation (26), and rearranging terms yields

\hat{F} (j ω) = (K_{S} (ω) + i ω \cdot C_{e} (ω)) \hat{x} (i ω)

(27)

where the real part, $K_{S} (ω),$ of the complex stiffness denotes the storage stiffness. Similarly, the imaginary part $(ω \cdot C_{e} (ω))$ , often named “loss stiffness,” refers to the damping characteristics of the DCHP suspension. The expressions for $K_{s} (ω)$ and $C_{e} (ω)$ can be expressed as follows

\begin{matrix} K_{S} (ω) & = \frac{(E_{0} + E_{1} ω^{α} \cos \frac{π}{2} α) ((1 + b ω^{β} \cos \frac{π}{2} β) + E_{1} ω^{α + β} b \sin \frac{π}{2} α \sin \frac{π}{2} β)}{1 + 2 b ω^{β} \cos \frac{π}{2} β + b^{2} ω^{2 β}} \\ c_{e} (ω) & = \frac{(E_{1} ω^{α} \sin \frac{π}{2} α) (1 + b ω^{β} \cos \frac{π}{2} β) - (E_{0} + E_{1} ω^{α} \cos \frac{π}{2} α) (b ω^{β} \sin \frac{π}{2} β)}{(1 + 2 b ω^{β} \cos \frac{π}{2} β + b^{2} ω^{2 β}) ω} \end{matrix}

(28)

Assuming the fractional factor $α = β$ , the number of parameters in the fractional Zener model can be reduced without perceivable reduction in accuracy.^24–26,28 The expressions for $K_{S} (ω)$ and $C_{e} (ω)$ can now be rewritten as simplified versions

\begin{matrix} {K'}_{S} (ω) = \frac{E_{0} + (E_{1} + E_{0} b) ω^{α} \cos \frac{π}{2} α + E_{1} b ω^{2 α}}{1 + 2 b ω^{α} \cos \frac{π}{2} α + b^{2} ω^{2 α}} \\ {C'}_{e} (ω) = \frac{(E_{1} - E_{0} b) ω^{α} \sin \frac{π}{2} α}{(1 + 2 b ω^{α} \cos \frac{π}{2} α + b^{2} ω^{2 α}) ω} \end{matrix}

(29)

Comparison of the simulation results of two different models

Equations (20) and (29) are the general expressions for the equivalent stiffness and damping coefficient of the nonlinear suspension system. However, the storage stiffness $(K'_{S} (ω))$ and the dynamic stiffness $(k_{d})$ have the same expression, and the equivalent damping $(C'_{e} (ω))$ is equal to the dynamic damping $(c_{d}) ω$ when assuming

E_{0} = \frac{k_{1} k_{2}}{k_{1} + k_{2}}, E_{1} = \frac{k_{1} c_{2}}{k_{1} + k_{2}}, b = \frac{c_{2}}{k_{1} + k_{2}}, α = 1

According to the simulation results for the different $α$ values shown in Figure 5, it is concluded that both the storage and the equivalent damping curves, $K'_{S} (ω)$ and $C'_{e} (ω)$ , display similar tendencies. However, they are not exactly the same as the behavior of the DCHP suspension shown in Figure 2. From Figure 5(a), it can be seen that the growth of the storage stiffness curves is small in the low- and high-frequency range, but it is high in the mid-frequency range. The intersection points of different curves are located within a narrow frequency band. When $α \leq 1,$ the stiffness curves rise monotonically. The slope of the curve in the mid-frequency range depends on the fractional derivative order $α$ , which plays the same role as the softening or hardening behavior for the additional damper.

Figure 5.

Computational simulation of the (a) storage stiffness and (b) equivalent damping coefficient for different values of the fractional derivative order α.

Both minimum and maximum storage stiffness can be calculated using equation (29) (a) at $ω = 0$ and $ω = \infty$

K'_{SMIN} (ω) = E_{0}, K'_{SMAX} (ω) = \frac{E_{1}}{b}

A set of equivalent curves obtained by varying the fractional derivative $α$ from 0.9 to 1 is shown in Figure 5(b). Equivalent damping curves with a lower $α$ are bigger in the low- and high-frequency range and the two intersection points located in frequency segments marked $ω_{1}$ and $ω_{2}$ . After comparing these curves in Figure 5(b), it is clear that an increase in the fractional derivative order $α$ changes the shape of the curve. When $α = 1$ , the curves for equivalent stiffness and damping coefficient have the same shape as the corresponding curves for linear damping in Figure 2. However, the decrease in the fractional derivative order $α$ enables softer damping in the suspension system. As a result, the corresponding stiffness curve $(α < 1)$ becomes more gentle in Figure 5(a).

Experiment results and discussion

Set-up of experiments

Two completely different theories are applied to describe the dynamic behavior of the DCHP suspension. The validity of equations (20) and (29), which are compatible with nonlinear vibration theory and fractional calculus, respectively, will be verified through the following experiments. To reduce complexity, the valve disk was removed from the piston, so that the damping force produced by main damper can be ignored. The adjustment of the additional damping coefficient will be realized via the throttle valve.

The DCHP suspension is installed in the test setup, which consists of a computer-controlled electro-hydraulic shaker platform (Japan Saginomiya), frame, S-type load cell, and a linear variable differential transformer (LVDT) displacement sensor. For the transmitted force measurements, the bottom of the DCHP spring is mounted onto the shaker head and the top is linked with the frame. The excitation frequency varies between 0.1 and 9 Hz. The displacement input signal, which is a sine wave with fixed amplitude (10 mm), is measured with an LVDT displacement sensor. The signal acquisition of the force and displacement data is accomplished using a data acquisition instrument (COINV INV3060V). Photos of the experimental device are shown in Figure 6.

Figure 6.

Photos of the experiment device.

Experimental results

In the experiment, the frequency of the excitation sine wave and the opening of the throttle valve are the parameters that were adjusted. All parameters for the DCHP spring used in simulations and experiments are listed in Table 1.

Table 1.

Parameters of the DCHP suspension.

Symbol	Value	Unit	Meaning
$p_{0}$	2.00	$MPa$	The initial pressure in both accumulators
A	706.86	$m m^{2}$	The sectional area of the rod
$V_{10}$	$0.16$	$L$	The initial volume of gas in the main accumulator
$V_{20}$	$0.16$	$L$	The initial volume of gas in the additional accumulator

The force (F) provided by the DCHP suspension consists of frictional force (F_f), damping force $(c \cdot X_{1} \cos θ)$ and spring force $(k \cdot X_{1} \sin θ)$

F = F_{f} + c \cdot X_{1} \cos θ + k \cdot X_{1} \sin θ

(30)

Here, $X_{1}$ is the amplitude of the harmonic excitation; k and c represent the linear stiffness and damping coefficient, which are parameters to be identified using the least-square method with different excitation frequencies.

Computational simulations in the frequency domain were carried out, and the simulation results were compared with the experiment (see Figures 7 and 8).

Figure 7.

Comparison of storage stiffness and equivalent damping coefficient between experiment and simulation for throttle opening level 1: (a) storage stiffness based on nonlinear vibration theory, (b) equivalent damping based on nonlinear vibration theory, (c) storage stiffness based on fractional calculus, and (d) equivalent damping based on fractional calculus.

Figure 8.

Comparison of storage stiffness and equivalent damping coefficient between experiment and simulation for throttle opening level 2: (a) storage stiffness based on nonlinear vibration theory, (b) equivalent damping based on nonlinear vibration theory, (c) storage stiffness based on fractional calculus, and (d) equivalent damping based on fractional calculus.

The behavior of the DCHP suspension strongly depends on the excitation frequency. Although the predictions for equivalent damping show some discrepancies with the experimental measurements for the low-frequency ranges, the behavior of the DCHP suspension is accurately predicted for the whole frequency range using the models based on both nonlinear vibration theory and fractional calculus theory.

Conclusion

This article describes a traditional hydro-pneumatic suspension connected with an additional accumulator and an adjustable throttle valve. Nonlinear vibration theory and fractional calculus theory were used to model the DCHP suspension. The expression for the storage stiffness and the equivalent coefficient was deduced and comprehensive simulations and experiments carried out. The frequency-dependent characteristic can be summarized as follows:

In the low-frequency range, both accumulators act as elastic elements, and the storage stiffness is numerically equal to the stiffness values of two springs in series. With increasing excitation frequency, the additional accumulator gradually stops functioning, resulting in only the main accumulator providing the elastic force.

An additional damping coefficient that can be controlled by a throttle valve is the most important parameter that causes the storage stiffness curve to shift along the frequency axis and affect the peak value of the equivalent damping.

The validity of two complete different models describing the behavior of the DCHP suspension was confirmed. Good agreements between analytical equations and experiment results were obtained.

We conclude that the performance of vibration systems can be improved significantly and optimized through suitable selection of the parameters in the DCHP suspension.

Footnotes

Academic Editor: Rahmi Guclu

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

El-Tawwab

Twin-accumulator suspension system. SAE technical paper 970384, 1997.

Erin

Wilson

Zapfe

An improved model of a pneumatic vibration isolator: theory and experiment. J Sound Vib 1998; 218: 81–101.

Sun

Deng

WH.

Variable stiffness and damping suspension system for train. Proc SPIE 2014; 9057: 90570P.

Sun

Deng

. A compact variable stiffness and damping shock absorber for vehicle suspension. IEEE-ASME T Mech 2015; 20: 2621–2629.

Lee

Kim

KJ.

Modeling of nonlinear complex stiffness of dual-chamber pneumatic spring for precision vibration isolations. J Sound Vib 2007; 301: 909–926.

Lee

Kim

KJ.

A method of transmissibility design for dual-chamber pneumatic vibration isolator. J Sound Vib 2009; 323: 67–92.

Porumamilla

Kelkar

Vogel

JM.

Modeling and verification of an innovative active pneumatic vibration isolation system. J Dyn Syst-T ASME 2008; 130: 031001.

Holtz

van Niekerk

JL.

Modelling and design of a novel air-spring for a suspension seat. J Sound Vib 2010; 329: 4354–4366.

Nieto

Morales

González

. An analytical model of pneumatic suspensions based on an experimental characterization. J Sound Vib 2008; 313: 290–307.

10.

Nieto

Moralesa

Trapero

. An adaptive pneumatic suspension based on the estimation of the excitation frequency. J Sound Vib 2011; 330: 1891–1903.

11.

Morales

Nieto

Chicharro

. An adaptive pneumatic system for the attenuation of random vibrations. J Vib Control 2015; 21: 907–918.

12.

Liu

Lee

JC.

Model development and experimental research on an air spring with auxiliary reservoir. Int J Automot Techn 2011; 12: 839–847.

13.

Luo

Chen

XD.

Modeling and analysis of dual-chamber pneumatic spring with adjustable damping for precision vibration isolation. J Sound Vib 2011; 330: 3578–3590.

14.

Chen

XD.

Generic vibration criteria-based dual-chamber pneumatic spring vibration isolation table design. Proc IMechE, Part B: J Engineering Manufacture 2014; 228: 1621–1629.

15.

Chai

Dreyer

Singh

Nonlinear dynamic properties of hydraulic suspension bushing with emphasis on the flow passage characteristics. Proc IMechE, Part D: J Automobile Engineering 2015; 229: 1327–1344.

16.

Ryaboy

VM.

Static and dynamic stability of pneumatic vibration isolators and systems of isolators. J Sound Vib 2014; 333: 31–51.

17.

Lang

HH.

A study of the characteristics of automotive hydraulic dampers at high stroking frequencies. PhD Dissertation, University of Michigan, Ann Arbor, MI, 1977.

18.

Cossalter

Doria

Pegoraro

. On the non-linear behaviour of motorcycle shock absorbers. Proc IMechE, Part D: J Automobile Engineering 2010; 224: 15–27.

19.

Duym

SWR

. Simulation tools, modelling and identification, for an automotive shock absorber in the context of vehicle dynamics. Vehicle Syst Dyn 2000; 33: 261–285.

20.

Yang

Brennan

. On the performance of a two-stage vibration isolation system which has geometrically nonlinear stiffness. J Vib Acoust 2014; 136: 064501.

21.

Tang

Brennan

MJ.

A comparison of two nonlinear damping mechanisms in a vibration isolator. J Sound Vib 2013; 332: 510–520.

22.

Worden

Tomlinson

GR.

Nonlinearity in structural dynamics: detection, identification and modelling. Boca Raton, FL: CRC Press, 2000.

23.

Park

SW.

Analytical modeling of viscoelastic dampers for structural and vibration control. Int J Solids Struct 2001; 38: 8065–8092.

24.

Lewandowski

Chorazyczewski

Identification of the parameters of the Kelvin-Voigt and the Maxwell fractional models, used to modeling of viscoelastic dampers. Comput Struct 2010; 88: 1–17.

25.

Pritz

Analysis of four-parameter fractional derivative model of real solid materials. J Sound Vib 1996; 195: 103–115.

26.

Pritz

Five-parameter fractional derivative model for polymeric damping materials. J Sound Vib 2003; 265: 935–952.

27.

Makris

Constantinou

MC.

Fractional-derivative Maxwell model for viscous dampers. J Struct Eng-ASCE 1991; 117: 2708–2724.

28.

Bagley

Torvik

PJ.

Fractional calculus: a different approach to the analysis of viscoelastically damped structures. AIAA J 1983; 21: 741–748.

Numerical analysis of a dual-chamber hydro-pneumatic suspension using nonlinear vibration theory and fractional calculus

Abstract

Keywords

Introduction

Structure of the DCHP suspension

Modeling and analysis of the DCHP suspension using nonlinear oscillation theory

Modeling of the pressure of the rebound chamber

Modeling of the pressure in the accumulators

Modeling of the damping force of the throttle valve

Simplification of the dynamic characteristic of the DCHP suspension

Solution of the equations for motion for harmonic excitation

Equivalent stiffness and damping coefficient

Modeling and analysis of the DCHP suspension based on fractional order calculus theory

Modeling of the DCHP suspension

Storage stiffness and equivalent damping coefficient

Comparison of the simulation results of two different models

Experiment results and discussion

Set-up of experiments

Experimental results

Conclusion

Footnotes

Declaration of conflicting interests

Funding

References