*1.3. Energy Harvesting Hydraulic Interconnected Suspension*

To combine the merits of the hydraulic interconnected suspension and the energy harvesting shock absorber, various types of the energy harvesting hydraulic interconnected suspension came into individuals' sight, with the aim to simultaneously enhance the vehicle dynamics performance and harvest energy.

Chen et al. [37] integrated the energy harvesting shock absorbers into the hydraulic interconnected suspension system and indicated that the system performed better than the traditional suspension in terms of rolling dynamics and could harvest 421 W energy at 4 Hz and 40 mm (peak) excitation. Guo et al. [38] proposed a hydraulic interconnected suspension system based on hydraulic electromagnetic shock absorbers, which only adopted one set of hydraulic motor-generator system and greatly reduced the cost while improving the energy recovery efficiency. The hydraulic integrated interconnected regenerative suspension (HIIRS) [39], which is studied in this paper, consists of a two-way hydraulic cylinder installed between the wheel and the body, an oil pipe connecting the hydraulic cylinders, a high-pressure accumulator, a low-pressure accumulator, two hydraulic rectifiers and a hydraulic motor-generator unit. Its structure and working principle are shown in Figure 1.

**Figure 1.** The structure and working principle of the HIIRS. 1 hydraulic cylinder; 2 high-pressure accumulator; 3 check valve; 4 hydraulic rectifier; 5 hydraulic motor; 6 generator; 7 low-pressure accumulator.

The hydraulic rectifier, composed of four check valves, ensures the one-way flow of fluid to drive the hydraulic motor. There is a high-pressure accumulator at the inlet of the hydraulic motor, which can stabilize the fluid flow through the hydraulic motor. This ensures the hydraulic motor to maintain a stable speed and the generator to generate electricity efficiently. The low-pressure accumulator at the outlet can compensate for the variation of the fluid volume in the HIIRS system. The high-pressure fluid flows through

the port of the high-pressure accumulator, and thus, when the high-pressure accumulator is working, it provides an extra rigidity to the suspension.

Recent studies on HIIRS mainly focused on design, modeling and experiments. Their results proved that the energy regenerative suspension had certain energy harvesting capability while ensuring comfort. However, the existing research did not model the HIIRS in the frequency domain. The frequency domain modeling method can easily allow for obtaining the evaluation index of vehicle ride comfort and the energy harvesting power at different road surfaces and vehicle speeds. Compared with the time domain analysis method, the frequency domain analysis method has several distinct advantages when applied to HIIRS. (1) The solution is simpler and more convenient; (2) the natural frequency can be calculated, which could guide the future design of the HIIRS; (3) the frequency domain analysis method can demonstrate the responses of the HIIRS under various excitations in an expedient way. In this paper, a half car model coupled with an HIIRS system is developed in the frequency domain. Based on the model, the vibration isolation characteristics and energy harvesting power of HIIRS are studied.

The rest of the paper is organized as follows. Section 2 develops the model of the HIIRS-equipped half vehicle in the frequency domain based on the block modeling and hydraulic impedance method. Free vibration analysis and forced vibration analysis are performed in Section 3. Energy harvesting power is calculated and estimated in Section 4. Finally, Section 5 concludes this paper.

#### **2. Modeling**

In this chapter, the idea of modular modeling is adopted. As the HIIRS is introduced into the vehicle, the whole system is divided into two parts: the mechanical system and the hydraulic system. The two parts and their boundary coupling conditions are discussed separately, and finally, the coupled dynamic equation is obtained.

#### *2.1. Mechanical System*

Considering the simplicity of modeling, while still accounting for fluid interconnections between the wheel stations and the lumped mass, a four-DOF half-car model, as shown in Figure 2, is used in this investigation.

**Figure 2.** Schematic of a half-car with an HIIRS system.

In this section, the mechanical system is what we concerned about, hence the force of the hydraulic system is regarded as an external force. According to Newton's second law, the kinematics equation of the half-car model is written by

$$\mathbf{M}\ddot{\mathbf{y}} + \mathbf{C}\dot{\mathbf{y}} + \mathbf{K}y = f(t) \tag{1}$$

where y is the displacement vector, y = [y*wl*,y*wr*,y*<sup>b</sup>* ,θ]; *f*(*t*) is the resultant force applied to the vehicle, which can be written as *f*(*t*) = *D*1*Ap*(*t*) + *fx*(*t*). In the equation of *f*(*t*), *D*1*Ap*(*t*) is the hydraulic cylinder force, and *f<sup>x</sup>* is other external forces vector. We define the pressure vector as *p*(*t*) = [*P*<sup>1</sup> *P*<sup>2</sup> *P*<sup>3</sup> *P*4] *T* and the area matrix as *A = diag*(*A*1, *A*2, *A*3, *A*4). *D*<sup>1</sup> 

is a linear transformation matrix *D*<sup>1</sup> = −1 1 0 0 0 0 −1 1 1 −1 1 −1 −*b<sup>l</sup> b<sup>l</sup> b<sup>r</sup>* −*b<sup>r</sup>* . *A<sup>i</sup>* is the cross-sectional

area of the hydraulic cylinder chamber in Figure 1, *b<sup>l</sup>* and *b<sup>r</sup>* are, respectively, the distances from the plane of gravity center of the vehicle to the left and right suspension.

Equation (1) can now be rewritten as

$$M\ddot{y} + \mathbf{C}\dot{y} + \mathbf{K}y = D\_1 A p(t) + f\_\mathbf{x}(t) \tag{2}$$

$$\begin{aligned} \text{where } M &= \begin{bmatrix} m\_l & 0 & 0 & 0 \\ 0 & m\_r & 0 & 0 \\ 0 & 0 & M & 0 \\ 0 & 0 & 0 & I \end{bmatrix}, C = \begin{bmatrix} c\_{tl} & 0 & 0 & 0 \\ 0 & c\_{tr} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}, \\ K &= \begin{bmatrix} k\_{sl} + k\_{tl} & 0 & -k\_{sl} & b\_{l}k\_{sl} \\ 0 & k\_{sr} + k\_{sl} & -k\_{sr} & -b\_{r}k\_{sr} \\ -k\_{sl} & -k\_{sr} & k\_{sl} + k\_{sr} & -b\_{l}k\_{sl} + b\_{r}k\_{sr} \\ b\_{l}k\_{sl} & -b\_{r}k\_{sr} & -b\_{l}k\_{sl} + b\_{r}k\_{sr} & b\_{l}^{2}k\_{sl} + b\_{r}^{2}k\_{sr} \end{bmatrix}. \end{aligned}$$

## *2.2. Mechanical–Fluid System Boundary Conditions*

In a double-acting hydraulic actuator cylinder, any piston movement will cause liquid to flow into or out of the chamber of the cylinder. Assuming the pressure difference between the upper and lower chambers produces linear leakage at the piston [40], then the linear leakage of the left and right actuators, *q<sup>l</sup>* and *q<sup>r</sup>* , are given by

$$\begin{cases} \begin{array}{c} q\_l = \frac{p\_{v1} - p\_{v2}}{R\_l} \\ q\_r = \frac{p\_{v3} - p\_{v4}}{R\_r} \end{array} \end{cases} \tag{3}$$

where *R<sup>l</sup>* and *R<sup>r</sup>* are linearized loss coefficients.

The mechanical–fluid boundary condition, shown in Figure 3, are given by

$$\begin{cases} q\_{\rm v1} = q\_{A1} - q\_l \\ q\_{\rm v2} = q\_{A2} - q\_l \\ q\_{\rm v3} = q\_{A3} - q\_r \\ q\_{\rm v4} = q\_{A4} - q\_r \end{cases} \tag{4}$$

where *qAi*(*i* = 1 2 3 4) is the fluid volume caused by piston motion, which equals the product of the relative speed of the piston *v* and the piston area *A*.

**Figure 3.** The mechanical–fluid system boundary condition.

For a small roll angle, the relative speed of left and right pistons is

$$\begin{cases} \ v\_l = y\_{wl}\dot{-y\_b} + b\_l \dot{\theta} \\\ v\_r = y\_{wr} - \dot{y\_b} - b\_r \dot{\theta} \end{cases} \tag{5}$$

Combining Equations (3)–(5), there are

$$q(t) = A D\_2 \dot{y} + R p(t) \tag{6}$$

where *p*(*t*) is the flow vector,*p*(*t*) = [*P*<sup>1</sup> *P*<sup>2</sup> *P*<sup>3</sup> *P*4] *T* .The matrix *D*<sup>2</sup> and *R* are defined as

$$D\_2 = \begin{bmatrix} 1 & 0 & -1 & b\_l \\ 1 & 0 & -1 & b\_l \\ 0 & 1 & -1 & -b\_r \\ 0 & 1 & -1 & -b\_r \end{bmatrix}, \\ R = \begin{bmatrix} -1/R\_l & 1/R\_l & 0 & 0 \\ -1/R\_l & 1/R\_l & 0 & 0 \\ 0 & 0 & -1/R\_r & 1/R\_r \\ 0 & 0 & -1/R\_r & 1/R\_r \end{bmatrix}.$$

#### *2.3. Fluid System*

In order to solve Equations (2) and (6), which relate to the mechanical system and mechanical–fluid boundary, the fluid system equation in the form *q* = *f*(*p*) must be obtained. *f*(*p*) depends on the modelling approach used. For the sake of computational efficiency and analytical advantages, only the linear function between *p* and *q* is considered. In particular, the focus of this study is frequency domain modeling. Therefore, the target is to seek the linear relationship between the flow rate *Q(s)* and the pressure *P*(*s*) in the frequency domain.

According to the definition of hydraulic impedance, the relationship between the flow rate and the pressure of a fluid system can be expressed as

$$Q(s) = Z(s)^{-1}P(s) \tag{7}$$

where *Z* is the impedance matrix.

Equation (6) describes the fluid system of the HIIRS, while Equation (7) describes the mechanical system of the HIIRS. Next, we must combine the two systems, to solve the HIIRS. By comparing the Laplace form of Equations (6) and (7), we can obtain the relationship between the pressure *P*(*s*) and the exciting displacement *Y*(*s*), as

$$P(s) = sE(s)AD\_2Y(s)\tag{8}$$

where *E(s)* = (*Z*(*s*) <sup>−</sup><sup>1</sup> <sup>−</sup> *<sup>R</sup>*) −1 . Without considering leakage, *E*(*s*) = *Z*(*s*).

The Laplace transform form of Equation (2) is the system differential equation in the frequency domain

$$s\text{Ms}^2Y(\text{s}) + s\text{C}Y(\text{s}) + K\text{Y}(\text{s}) = s\text{D}\_1AE(\text{s})AD\_2\text{Y}(\text{s}) + F\_\text{x}(\text{s})\tag{9}$$

$$[s^2M + s\overline{C}(s) + \mathcal{K}]Y(s) = F\_{\mathcal{X}}(s) \tag{10}$$

where *C*(*s*) = *C* − *D*1*AE*(*s*)*AD*2.

Defining the state vector *x* = [*y*, *sy*] *T* , Equation (10) can be rewritten as

$$sX(s) = \hat{A}(s)X(s) + \hat{B}\mathcal{U}(s) \tag{11}$$

where *A*ˆ(*s*) = 0 *I* <sup>−</sup>*M*−1*<sup>K</sup>* <sup>−</sup>*M*−1*C*(*s*) , *B*ˆ = 0 *M*−<sup>1</sup> , *U*(*s*) = *Fx*(*S*).

It should be noted that it is not easy to determine the impedance correlation matrix *A*ˆ(*s*), which depends on the fluid circuit and the arrangement of the various components. The circuit layout in this study adopts anti-oppositional interconnection, as shown in Figure 4.

**Figure 4.** Schematic of a general anti-oppositional half-car HIIRS arrangement.

Using the original boundary flow definitions, *q<sup>i</sup>* , the nodal state vectors for the mechanical–fluid interfaces in this arrangement are related by

$$
\begin{bmatrix} P\_4 \\ Q\_4 \end{bmatrix} = \begin{bmatrix} T\_{11}^a & T\_{12}^a \\ T\_{21}^a & T\_{22}^a \end{bmatrix} \begin{bmatrix} P\_1 \\ Q\_1 \end{bmatrix} \begin{bmatrix} P\_2 \\ Q\_2 \end{bmatrix} = \begin{bmatrix} T\_{11}^b & T\_{12}^b \\ T\_{21}^b & T\_{22}^b \end{bmatrix} \begin{bmatrix} P\_3 \\ Q\_3 \end{bmatrix} \tag{12}
$$

Combining Equation (12) with Equation (7), *Z*(*s*) −1 can be written as

$$Z(s)^{-1} = \begin{bmatrix} -\frac{T\_{22}^{a}}{T\_{21}^{b}} & 0 & 0 & \frac{1}{T\_{21}^{a}}\\ 0 & \frac{T\_{11}^{b}}{T\_{21}^{b}} & -\frac{T\_{11}^{b}T\_{22}^{b}}{T\_{21}^{b}} & 0\\ 0 & \frac{1}{T\_{21}^{b}} & \frac{-T\_{22}^{b}}{T\_{21}^{b}} & 0\\ T\_{12}^{a} - \frac{T\_{11}^{a}T\_{22}^{a}}{T\_{21}^{a}} & 0 & 0 & \frac{T\_{11}^{a}}{T\_{21}^{a}} \end{bmatrix} \tag{13}$$

By substituting Equation (13) into Equation (12), it yields the complete system equations. For Equation (13), the values of the elements in the matrix *T <sup>a</sup>* and *T b* can be obtained by the hydraulic impedance method and the transfer matrix method. In this method, the state vector of the adjacent state node is related to the transfer matrix *T*. If the state vector is defined as fluid pressure *P* and flow rate *Q*, then

$$
\begin{bmatrix} P \\ Q \end{bmatrix}\_o = \begin{bmatrix} T\_{11} & T\_{12} \\ T\_{21} & T\_{22} \end{bmatrix} \begin{bmatrix} P \\ Q \end{bmatrix}\_i \tag{14}
$$

where the subscript *o* represents the fluid output node, and *i* represents the fluid input node. As long as the previous output node is regarded as the next input node, according to the fluid flow direction in the circuit, the start node and the end node in the circuit can be connected.

As the layout of the hydraulic circuit is shown in Figure 5, the transfer matrix *T a* can be expressed as

*T a* (*s*) = T*X*11→*X*12*TX*10→*X*11*TX*9→*X*10*TX*8→*X*9T*X*7→*X*8*TX*6→*X*7T*X*5→*X*6*TX*4→*X*5T*X*3→*X*4T*X*2→*X*3T*X*1→*X*<sup>2</sup> (15)

**Figure 5.** Schematic of a typical half-car with the HIIRS system. v1—check valves; la—low-pressure accumulator; ha—high-pressure accumulator; v2—accumulator valve; x, x0—state nodes.

The transfer matrix of each segment in the circuit depends on the characteristics of the component, which are modeled as follows.

The two-dimensional viscous compressible flow model is applied to model the pipeline, according to [41–43]. The basic fluid equations of the model can be summarized into the equation of state, continuity equation and momentum equation:

Equation (16) of state:

$$\frac{dp}{d\rho} = a\_0^2 \tag{16}$$

Continuity Equation (17):

$$\frac{\partial p}{\partial t} + \overline{\rho} \left( \frac{\partial v\_x}{\partial x} + \frac{\partial v\_r}{\partial r} \frac{v\_r}{r} \right) = 0 \tag{17}$$

Momentum Equation (18):

$$
\overline{\rho}\frac{\partial v\_x}{\partial t} = \frac{-\partial p}{\partial x} + \overline{\mu}\left(\frac{\partial^2 v\_x}{\partial r^2} + \frac{1}{r}\frac{\partial v\_x}{\partial r}\right) \tag{18}
$$

where *v<sup>x</sup>* and *v<sup>r</sup>* are the axial and radial velocity of fluid; *a*<sup>0</sup> = q *β*/*ρ* is the fluid sonic velocity; *β*, *µ* and *ρ* are the mean value of fluid bulk modulus, viscosity and density. The field transfer matrix of the pipelines is

$$T\_p = \begin{bmatrix} \cosh \Gamma(s) & -Z\_{\mathbb{C}}(s) \sinh \Gamma(s) \\ -\frac{\sinh \Gamma(s)}{Z\_{\mathbb{C}}(s)} & \cosh \Gamma(s) \end{bmatrix} \tag{19}$$

$$\text{where } Z\_{\mathbb{C}}(\mathbf{s}) = \frac{d\overline{\rho}}{d} \left( l - \frac{2l\_{1} \left( \text{ir} \sqrt{s/\overline{\upsilon}} \right)}{\text{ir} \sqrt{s/\overline{\upsilon}} l\_{0} \left( \text{ir} \sqrt{s/\overline{\upsilon}} \right)} \right)^{-1/2}, \Gamma(\mathbf{s}) = \frac{l\overline{s}}{d} \left( l - \frac{2l\_{1} \left( \text{ir} \sqrt{s/\overline{\upsilon}} \right)}{\text{ir} \sqrt{s/\overline{\upsilon}} l\_{0} \left( \text{ir} \sqrt{s/\overline{\upsilon}} \right)} \right)^{-1/2}, \text{ and }$$

*l* is the length of the line element, *A* is the pipeline internal cross-sectional area, *r* is the pipeline internal radius, *ρ* is the mean fluid density, *v* is the mean fluid kinematic viscosity, and *J*<sup>0</sup> and *J*<sup>1</sup> are Bessel functions of the first kind with orders zero and one, respectively. Considering the compressibility of fluid and pipeline, the propagation velocity of transmission line is

$$\mathbf{a} = \sqrt{\frac{a\_0^2}{1 + \left(2ra\_0^2 \overline{\rho} / t\_p E\right)}}\tag{20}$$

where *t<sup>p</sup>* is the thickness of the pipe wall, and *E* is the Young's modulus of the pipe material.

The strict model of the flow dynamics through the damping valves usually involves complicated geometrical parameters. The simplified model used here assumes that the damping valve has a negligible fluid volume. The transfer matrix of the damping valves is then expressed as

$$
\Omega\_V = \begin{bmatrix} 1 & -Z\_v \\ 0 & 1 \end{bmatrix} \tag{21}
$$

where *Z<sup>v</sup>* = *Rv*, and *R<sup>v</sup>* is the constant linear pressure loss coefficient.

For the model of the accumulator, the following assumptions have been made. The compressibility of the liquid is far lower than that of the gas in the accumulator; the elasticity of the diaphragm is neglected; there is no heat exchange between the gas in the accumulator and the outside world. The axle moves quickly relative to the small amplitude of the car body, resulting in a rapid reduction or expansion in the volume of the accumulator, and satisfies the gas adiabatic balance equation. At this time, the accumulator can be regarded as a linear system [16], and the impedance of the accumulator is

$$Z\_{\rm A} = -\frac{\gamma \overline{P}^2}{sp\_p v\_p} \tag{22}$$

The three-way junctions and the accumulator are connected, and hence, the transfer matrix of the three-way junctions and the accumulator are combined here. According to impedance definition and accumulator impedance

$$P\_{\mathbf{X}\_{\ $,2}} = Q\_{\mathbf{X}\_{\$ ,2}} Z\_{\mathbf{A}} \tag{23}$$

Set the forward flow direction to be outside the accumulator

$$P\_{X5.1} = P\_{X5.2} - Z\_{V2}Q\_{X5.2} \text{ and } Q\_{X5.2} = Q\_{X5.1} \tag{24}$$

Combining Equations (23) and (24), we obtain

$$\frac{P\_{X5.1}}{Q\_{X5.1}} = Z\_A - Z\_{V2} = Z\_{X5.1} \tag{25}$$

Now, applying the fluid continuity equation at the tee-junction yields

$$Q\_{X5} = Q\_{X4} + \frac{P\_{X5.1}}{Z\_{X5.1}}\tag{26}$$

Ignoring the pressure difference between the three nodes, the transfer matrix can be written as

$$
\Omega\_{\bar{I}} = \begin{bmatrix} 1 & 0 \\ \frac{1}{\mathbb{Z}\_A - \mathbb{Z}\_{V2}} & 1 \end{bmatrix} \tag{27}
$$

The flow rate at the inlet of the hydraulic motor is defined as *QM*. One part of *QM*, which drives the rotation of the hydraulic motor is *ηVQM*, where *η<sup>v</sup>* is the volumetric efficiency of the hydraulic motor. The other part of the flow is the leakage flow ∆*Q<sup>M</sup>* from

the high-pressure cavity to the low-pressure cavity. The flow at the outlet of the hydraulic motor is *ηVQM*.

When the oil flows through the motor, the relationship between the flow rate *Q<sup>m</sup>* of the hydraulic motor and the speed *n<sup>m</sup>* (rev/s) of the hydraulic motor is

$$m\_{\rm m} = \frac{Q\_M \eta\_{\rm v}}{Q\_{\rm m}} \tag{28}$$

Assuming that the equivalent moment of inertia of the motor-generator coupling system is *Jm*, the speed of the motor is *nm*, *nm*= *ω* 2*π* , and ω is the angular velocity. Then, there is a torque balance .

$$T\_{m2} = T\_\varepsilon + f\_m \dot{\omega} \tag{29}$$

With reference to the principle of hydraulic transmission, the input torque *Tm*<sup>1</sup> and output torque *Tm*<sup>2</sup> of the hydraulic motor can be calculated according to the following Equation (30)

$$\begin{cases} \begin{array}{c} T\_{m1} = \frac{\Delta P\_m Q\_m}{2\pi} \\ T\_{m2} = T\_{m1} \eta\_m \end{array} \end{cases} \tag{30}$$

where ∆*p<sup>m</sup>* is the pressure drop, *η<sup>m</sup>* is the total efficiency of the hydraulic motor, *η<sup>m</sup>* = *ηvη<sup>p</sup>* is the product of the volumetric efficiency and pressure utilization efficiency.

Combining Equations (26)–(30), we obtain

$$\frac{\Delta P\_m Q\_m}{2\pi} \eta\_m = T\_\varepsilon + 2\pi f\_m n\_m^\cdot \tag{31}$$

The electromagnetic induction torque *T<sup>e</sup>* of the generator is determined by the armature current *I* and the torque constant *k<sup>t</sup>* of the motor

$$T\_{\mathcal{E}} = k\_l I \tag{32}$$

The winding armature current *I* is related to the design of the energy-regenerative circuit, and the equation can be obtained from Kirchhoff's voltage law, as

$$
\mathcal{U}\_{emf} = k\_\varepsilon \omega \tag{33}
$$

Combining Equations (31)–(33) to obtain

$$
\Delta p\_m = \frac{4\pi^2 k\_e k\_l \eta\_v}{(R\_1 + R\_2) Q\_m^2 \eta\_m} Q\_M + \frac{4J\_m \pi^2 \eta\_v}{Q\_m^2 \eta\_m} \dot{Q}\_M \tag{34}
$$

Laplace transform of the above Equation (34) yields

$$
\Delta P\_m = \frac{4\pi^2 k\_e k\_l \eta\_v}{(R\_e + R\_{\rm in}) Q\_m^2 \eta\_m} Q\_M + \frac{4J\_m \pi^2 \eta\_v}{Q\_m^2 \eta\_m} s Q\_M \tag{35}
$$

where *R<sup>e</sup>* is the external resistance, *Rin* is the circuit resistance, *k<sup>e</sup>* is the motor speed constant, *Q<sup>M</sup>* is the motor inlet flow, *Q<sup>m</sup>* is the motor displacement, *η<sup>v</sup>* is the volumetric efficiency, and *J<sup>m</sup>* is the motor-generator rotational inertia.

Therefore, the impedance matrix of the motor-generator unit is

$$Z\_M = \frac{\Delta P\_M}{Q\_M} = \frac{4\pi^2 k\_e k\_l \eta\_v}{(R\_e + R\_{\rm in}) Q\_m^2 \eta\_m} + \frac{4J\_m \pi^2 \eta\_v}{Q\_m^2 \eta\_m} s. \tag{36}$$

The transfer matrix of the combined model of hydraulic motor and generator can be written as

$$
\begin{bmatrix} P \\ Q \end{bmatrix}\_{out} = \begin{bmatrix} 1 & -Z\_M \\ 0 & \eta\_v \end{bmatrix} \begin{bmatrix} P \\ Q \end{bmatrix}\_{in} \tag{37}
$$

Now, the hydraulic impedance and transfer matrix of the pipeline, check valve, accumulator and motor-generator are obtained. Substituting Equations (19), (21), (27) and (36) into Equation (15), *T <sup>a</sup>* and *T b* can be clear. By substituting the relative elements of *T a* and *T b* into Equation (13), the matrix *Z* −1 (*s*) is obtained with definite elements. As a result, the governing equation of the HIIRS system is determined. Then, Equations (10) and (11) can be applied to solve the vibration problem of the HIIRS-equipped half-vehicle.

#### **3. Vibration Analysis**

After the above work, the system Equation (11) can be used to analyze the vibration characteristics of the half-car model equipped with HIIRS. The values of each parameter in this study are shown in Table 1.


**Table 1.** Nomenclature.

#### *3.1. Free Vibration Analysis*

When the external input is zero input, the half-vehicle system vibrates freely, and the corresponding mathematical description is the homogeneous form of Equation (11) as

$$\mathbf{s}\mathbf{X}(\mathbf{s}) = \hat{A}(\mathbf{s})X(\mathbf{s})\tag{38}$$

The solution of the free vibration of the system can be obtained by solving

$$\det(\hat{A}(\mathbf{s}) - \mathbf{s}\mathbf{I}) = \mathbf{0} \tag{39}$$

To solve Equation (38), we must solve the eigenvalue of *A*ˆ(*s*). Several elements of *A*ˆ(*s*) are functions of *s*, unless the frequency is known, *A*ˆ(*s*) cannot be completely determined, and therefore, the solution of Equation (38) cannot be obtained by conventional methods. The method used here transforms the process of finding the root of the characteristic equation into the process of finding the local optimal solution.

Assuming that the minimum value of the characteristic equation is the objective function, that is, min J(s) = *<sup>A</sup>*ˆ(*s*) <sup>−</sup> *SI* is the objective function, the fminsearch function in MATLAB is used to find the local minimum of the objective function J(s). The specific method is mainly divided into two processes. In the first process, the relatively rough root finding is performed. The value of the objective function is calculated after initializing the Laplacian operator. If the local extremum is not found, the value of the Laplacian operator is reinitialized, and the above process is performed again. If the local extremum is found, then the following second process is performed. The second process is a test process, and the Laplacian found in the first process is substituted into the characteristic matrix to obtain a fixed eigen matrix. The eigenvalue of the fixed value matrix is then solved, and the value is compared with the local extremum found in the first process. If the two values are not equal, skip the second process and return to the first process to find the local eigenvalue again. If the two values are the same, the root finding process ends.

With the determined eigenvalues *λ<sup>i</sup>* and eigenvectors *α<sup>i</sup>* from Equation (38), the natural frequency and damping ratio for each mode are given by

$$f\_{n\_i} = \frac{\sqrt{\left(real(\lambda\_i)\right)^2 + \left(imag(\lambda\_i)\right)^2}}{2\pi} \text{ and } \mathfrak{f}\_i = \frac{|\left(real(\lambda\_i)\right)|}{\sqrt{\left(real(\lambda\_i)\right)^2 + \left(imag(\lambda\_i)\right)^2}}\tag{40}$$

Based on the analysis of the above complex modal vibration theory, as long as the eigenvalues of the characteristic matrix in Equation (37) and the corresponding eigenvectors are solved, the natural frequency, damping ratio and main vibration mode of the system can be obtained. In order to find the roots conveniently, the three-dimensional graph is obtained, whose horizontal and vertical coordinates are the real and imaginary parts of the Laplacian operator, and the vertical coordinates are the objective function to initially determine the number of roots and the range of the roots according to the relevant parameters of the system. The three-dimensional image of the half-vehicle model obtained is shown in Figure 6.

**Figure 6.** Three-dimensional plot showing the four approximate roots of the characteristic equation of the HIIRS-equipped vehicle.

It is not difficult to see from Figure 6 that the system eigen matrix has four eigenvalues. The system eigenvalues obtained by the local optimization method after initially determining the range of the eigenvalues are shown in Table 2.


**Table 2.** The approximate eigenvalue solutions to the system matrix of the HIIRS-equipped vehicle.

The relatively rough roots obtained in the first process need to be tested in the test process. Substituting the eigenvalues in Table 2 into matrix Aˆ (*s*), the system matrix can be completely determined. The traditional methods of finding eigenvalues and eigenvectors can be applied. If the approximate root in process one is equal to the result obtained by the traditional method, it can be considered that the approximate root is the eigenvalue of matrix Aˆ (*s*). The natural frequency, damping ratio and mode shape corresponding to the eigenvalues that passed the inspection are shown in the Table 3.

**Table 3.** Vibration modes of the HIIRS-equipped vehicle.


The results show that the half-vehicle roll model has four main vibration modes. The first-order vibration mode is dominated by the vertical vibration of the vehicle body and the corresponding natural frequency is 1.49 Hz and the damping ratio is 0.40. Usually, this natural frequency of off-road vehicles is 1.3~2 Hz, and the damping ratio is 0.2~0.4 [44]. The second-order mode is dominated by the body roll vibration. The third-order vibration is dominated by wheel vibration, and the two wheels move in the same direction. The fourth-order vibration mode is dominated by wheel vibration, and the wheels on the left and right sides do reverse vibration. At this time, the body does not vibrate vertically, but has a slight roll vibration.
