Dynamic Finite Element Analysis of Flexible Double Wishbone Suspension Systems with Different Damping Mechanisms

Alaa A. Abdelrahman*, Ayman E. Nabawy, Soliman S. Alieldin and Ayman M. M. Abdelhaleem

Mechanical Design and Production Dept., Faculty of Engineering, Zagazig University, Zagazig, Egypt

Email: alaaabouahmed@gmail.com

* Corresponding Author

Received 28 June 2019; Accepted 05 February 2020; Publication 09 March 2020

Abstract

Suspension systems in running vehicles keep the occupants comfortable and isolated from road noise, disturbances, and vibrations and consequently prevent the vehicle from damage and wearing. To attain comfortable and vibration isolation conditions, both material flexibility and damping should be considered in the considered suspension model. This paper presents an incremental finite element model to study and analyze the dynamic behavior of double wishbone suspension systems considering both material flexibility and damping effects. The flexibility of the suspension links are modeled with plane frame element based on Timoshenko beam hypothesis (TBH). On the other hand, the flexibility of joints connecting the suspension links together and with the vehicle chassis is modeled with the revolute joint element. To incorporate the damping effect, viscoelastic, viscous and proportional damping are considered. An incremental viscoelastic constitutive relations, suitable for finite element implementation, are developed. The developed finite element equations of motion are solved using the Newmark technique. The developed procedure is verified by comparing the obtained results with that obtained by the developed analytical solution and an excellent agreement is found. The applicability and effectiveness of the developed procedure are demonstrated by conducting parametric studies to show the effects of the road irregularities profiles, the vehicle speed, and the material damping on the transverse deflection and the resultant stresses of suspension system. Results obtained are supportive in the mechanical design, manufacturing processes of such type of structural systems.

Keywords: Viscoelastic, incremental finite element, double wishbone suspension, links flexibility, revolute joint, road irregularities, Timoshenko beam hypothesis (TBH).

1 Introduction

Through the study and analysis of automotive engineering it is most important to study how disturbances from the engine are led into the coupe. Also, one of the most important interest to investigate how disturbances from the roadway propagate through tire, and wheel suspension into the framing and car body. These disturbances may cause undesired noise levels and vibrations, Jonsson [1]. The main objectives of the vehicle suspension systems are listed as follows: Carrying the vehicles and its weight. Maintaining the correct alignment of the wheels. Controlling the travel direction of the vehicles. Finally, minimizing the chock forces effect. Selection of a certain type of suspension system in automotive industry depends on important criteria such as costs, space requirements, kinematic properties, and compliance attributes Rill [2].

Suspension system in its simplest form may be thought of as a set of flexible links connected to each other and to the vehicle chassis through ball joints to allow the wheel to move relative to the body in addition to some elastic element to support loads while allowing that motion, Kamal and Rahman [3]. One of the most important parts in the vehicle suspension systems is the lower control arm which is called as A-type control arm by which the wheel trace is controlled and the load exerting on the wheel by the road is transmitted to other parts of the car, Zhao et al. [4]. In addition, the lower control arm is considered as the most critically loaded member as the load is transferred from the upper arm to the lower control arm causing impact load due to the road irregularities (road bumps), Patil et al. [5]. To avoid failure of the lower control arm, it is essential to focus on the stress analysis study of lower wishbone arm to improve and modify the existing design, these stresses are caused by the dynamic forces which are caused by the dynamic deflections, Dehkordi [6].

Due to the difficulties associated with the boundary conditions and geometrical and/ or material nonlinearities, it is difficult to obtain an analytical closed form solution for such type of problems, the Finite Element Method (FEM) is one of the most effective numerical techniques which can give better conception of the dynamical behavior of this kind of analysis for these systems.

Many considerable research efforts have been dedicated to the design and analysis of the suspension systems.

Based on interval analysis, the kinematic analysis of the mechanisms involved in suspension of vehicles is performed using some exact computational methods, Papegay et al. [7]. Based on gradient determination using exact differentiation, a multi-objective dimensional synthesis technique is used to present the kinematic analysis of double wishbone suspension systems in vehicles, Sancibrian et al. [8]. Using Rodrigue’s parameters and computation of Gröbner basis, the double wishbone and the MacPherson, Reddy et al. [9] developed a complete spatial model for independent suspension systems to investigate the kinematics analysis of these systems. The variable-geometry suspension systems have been synthesized and analyzed in nonlinear manner using the nonlinear polynomial Sum-of-Squares (SOS) programming method Németh and Gáspár [10]. Using variable length arms in an adaptive manner to improve the handling characteristics of a vehicle, camber and toe angles of a double wishbone suspension system are controlled, Kavitha et al. [11]. Based on the point joint coordinate formulation, the dynamic analysis of the double wishbone motor vehicle suspension system is discussed, Attia [12]. Girder and Hos-sack suspension systems, which have the common double wishbone design have been accommodated by modifying the front end geometry based on an existing mathematical model of high-fidelity motorcycle, Ramírez et al. [13].

Many research articles discussed the dynamic behavior of the critically loaded suspension lower control arm. FEM-based suspension arm model and the modal frequencies are established to determine the loads affecting an arm operating in a complete geometric model of a car suspension, Burdzik, and Łazarz [14]. Using the FEM, the fatigue life, von Misses stress, factor of safety and stability of the suspension system are analyzed. Moreover, the weight and size of the suspension system which affect the fuel economy can be reduced, Ijagbemi et al. [15]. By incorporating the frame flexibility and the contacts of the wheel ground, suspension, unsprung mass, and compliance parameters are incorporated to investigate the dynamic behavior using a linear quarter car model, Kanchwala and Chatterjee [16].

To illustrate, the different performance of elastic and composite lower control arms, a comparative study of composite lower wishbone arm and steel lower wishbone arm is performed to investigate the stress-strain analysis using the FEM, Ramesh et al. [17]. The upper control arm has been optimized using Kriging model and the response surface model (RSM) in which the design objective is the control arm weight and the constraint objective is the maximum allowable von Mises stress, Song et al. [18]. Traditional static load optimization techniques and multi-body dynamics are combined using equivalent static load to perform the dynamic optimization of the lower control arm, Zhao et al. [4]. The lower arm of the vehicle suspension system was found to be the most important design variable to achieve noise reduction using sensitivity analysis performed by the FEM determining sensitive design factors, these factors are obtained also from the experimental vehicle tests which compared with those obtained from a robust and efficient sequential approximate optimization method, Kim et al. [19]. Baranowski et al. [20] developed numerical and experimental models for investigating the dynamic response of vehicle tyres under impulse loading conditions. Considering the nonlinear suspension components behaviors, Sayyaadi and Shokouhi [21] developed complete four axle rail vehicle model for investigating vehicle dynamics.

Many research articles present the dynamic behavior of beams based on Timoshenko theory, all Timoshenko beam forced and natural boundary conditions are satisfied using the finite element model, Thomas and Abbas [22]. Analytical expressions for the terms of the coupled bending-torsional dynamic stiffness matrix are derived for a Timoshenko beam with a uniform axial load in an exact sense, Williams [23]. The rail vertical vibration reac-ceptances are calculated by developing a double Timoshenko beam model using both continuously and discretely supported models, Thompson [24]. A cracked Timoshenko beam is presented using a new spectral finite element for the modal and elastic wave propagation analysis, Krawczuk et al. [25]. The dynamic response of the Timoshenko cantilever beam is evaluated with a homogeneous non-prismatic pulse loaded using the Lagrange-d’Alembert principle, Navadeh et al. [26].

According to the best of authors’ knowledge, there is no attempts to develop a finite element model capable of studying and analyzing the dynamic behavior of vehicle suspension systems considering both joint and links flexibilities as well as the internal material damping effects. Thus, the present article aims to fill this gap. Based on Timoshenko beam theory, the plane frame element is adopted to model the flexibility of links while the revolute joint element is used for joint flexibility. The internal viscoelastic as well as the external viscous damping models are considered to investigate the damped response. The developed numerical model is validated and checked by comparing the obtained results with the developed analytical solution and good agreement is obtained.

2 Damping Mechanisms

In spite of large amount of research articles that addressed damping mechanisms in solids and structures, understanding of damping mechanisms is quite basic compared to the other aspects of modelling. This because, in contrast to stiffness and inertia forces, the state variables relevant in determining the damping force are not generally clear, Inman [27]. The most common damping mechanisms in solids and structures are the viscous and the viscoelastic damping mechanisms. In viscously damped systems the damping force is proportional to velocity. The proportionality constant is the damping coefficient. Identifying the material viscous damping coefficient is a major problem. One of the most common idealization to identify the viscous damping coefficient is the Rayleigh or proportional damping model. In this model the viscous damping coefficient is presented by a linear combination of the mass and stiffness of the system, Nakamura [28].

$$\mathrm{<mo\; stretchy="false">[</mo>}C\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>}\alpha \mathrm{<mo\; stretchy="false">[</mo>}M\mathrm{<mo\; stretchy="false">]</mo>}+\beta \mathrm{<mo\; stretchy="false">[</mo>}K\mathrm{<mo\; stretchy="false">]</mo>}$$

(1a)

where [C], [M], [K] are the damping, mass and stiffness matrices; respectively. While α and β are the mass and stiffness matrices multipliers contestants; respectively. These constants are obtained through modal analysis. The damping ratio for the nth mode of such a system is, Cook [29], Chopra [30]:

$${\xi }_n=\frac{1}{2}\left[\frac{\alpha }{{\omega }_n}+\beta {\omega }_n\right]$$

(1b)

Equation (1b) can be expressed for the two modes i and j as

$$\left\{\begin{array}{c}{\xi }_i\\ {\xi }_j\end{array}\right\}=\frac{1}{2}\left[\begin{array}{cc}\frac{1}{{\omega }_i}& {\omega }_i\\ \frac{1}{{\omega }_j}& {\omega }_j\end{array}\right]\left\{\begin{array}{c}\alpha \\ \beta \end{array}\right\}$$

(1c)

Assuming constant damping factor, ξ for two modes, i and j, the mass and stiffness multipliers can be expressed as

$$\alpha =\frac{2\xi {\omega }_i{\omega }_j}{{\omega }_i+{\omega }_j}\mathrm{,}\beta =\frac{2\xi }{{\omega }_i+{\omega }_j}$$

(1d)

Thus proportional damping systems preserve the simplicity of the real normal modes as in undamped systems.

Considering the viscoelastic damping mechanism, according to the superposition principle, the linear viscoelastic constitutive relations can be represented by the convolution integral, Lakes and Lakes [31], Shames [32]

$$\begin{array}{l}{\sigma }_{ij}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\int }_{-\infty }^t}{R}_{ijkl}\left(t-s\right)\frac{\partial {\epsilon }_{kl}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}}{\partial s}ds\hfill \\ =R_{ijkl}\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\epsilon }_{kl}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ -{\displaystyle {\int }_0^t}\frac{d{R}_{ijkl}\left(t-s\right)}{ds}{\epsilon }_{kl}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}ds\hfill \end{array}$$

(2)

where σ_ij(x,t) and ε_kl(x, t) are stress and strain tensors; respectively, R_ijkl(t − s) is the relaxation modulus, and s is a time variable which defines the time at which the excitation input is applied. The relaxation modulus depends on the viscoelastic mechanical model used to simulate the viscoelastic constitutive response. The viscoelastic behavior is characterized by displacement creep under constant load and stress relaxation under constant prescribed displacement. Modeling the viscoelastic constitutive response using simple Kelvin Voigt or Maxwell mechanical models leads to inaccurate investigation of the real viscoelastic material response. The Maxwell model is good in predicting the stress relaxation response while it can poorly predict the creep response. On the other hand, Kelvin-Voigt is efficient in predicting the creep response while it is not effective for investigating the stress relaxation response. To overcome the drawbacks associated with the simple viscoelastic mechanical models, generalized viscoelastic mechanical models are used. Theses generalized models are generally constructed by using chains of either Maxwell models in parallel combinations or kelvin-Voigt elements in series combination. Throughout this study, the generalized Wiechert model, shown in Figure 1 is used to simulate the linear viscoelastic constitutive response, Abdelrahman and Mahmoud [33]. According to the generalized Wiechert model, the relaxation modulus can be repressed by

$$R_{ijkl}\mathrm{<mo\; stretchy="false">(</mo>}t-s\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{E}_{ijkl\infty }+{\displaystyle \sum _{n=1}^N}{E}_{jjkln}\left(1-e^{\frac{-\mathrm{<mo\; stretchy="false">(</mo>}t-s\mathrm{<mo\; stretchy="false">)</mo>}}{{\rho }_{ijkln}}}\right)\mathrm{,}{\rho }_{nijkl}=\frac{{\eta }_{jjkln}}{E_{jjkln}}$$

(3)

Figure 1 Wiechert viscoelastic material model.

where E_ijkl∞ and E_ijkln are the fully relaxed and transient modulus of the viscoelastic material, ρ_n is the material relaxation time and N is the number of Maxwell chains in the generalized Wiechert model. These constants are determined experimentally through either creep or relaxation tests.

3 Mathematical Formulation of Viscoelastically Damped Suspension System

Consider the quarter car model with a planar double wishbone suspension model as shown in Figure 2. According to the proposed planar suspension model, both lower and upper control arms are made of deformable materials (flexible). Moreover, flexible bushings or revolute joints used to connect these arms with vehicle chassis and wheel spindle are modeled by spring elements. On the other hand, suspension strut is modeled by a combination of translational spring and dashpot while tire is simulated by translational spring. To investigate the material damping effect on the dynamic behavior of the system, both viscous and viscoelastic damping models are considered for flexible suspension links, Bahrampour et al. [34]. The tire motion over the road bump, y(t) provides a vertical input which excites the body of the vehicle. Based on the displacement fit conditions and kinematic relations of Timoshenko beam theory (TBT), Reddy [35, 36], Applying the Hamilton’s principle, the dynamic equations of motion of viscoelastically damped systems can be written as, Martin [37] and Hilton [38].

$${\displaystyle {\int }_{t_1}^{t_2}}\left(\delta T-\delta U+\delta W_d+\delta W_E\right)dt=0$$

(4)

Figure 2 Quarter-car model of a double wishbone suspension.

where t₁ and t₂ are two arbitrary times, δU and δT are the variations of the strain energy and kinetic energy, respectively which can be expressed as

$$\begin{array}{l}\delta U={\displaystyle {\int }_0^L}{\displaystyle {\iint }_A}\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_{xx}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\delta {\epsilon }_{xx}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}+{\tau }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\delta {\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}dAdx\mathrm{,}\hfill \\ \delta T=\delta \left(\frac{1}{2}{\displaystyle {\int }_0^L}{\displaystyle {\iint }_A}\rho \left({\dot{U}}_x^2+{\dot{U}}_y^2\right)dAdx\right)\hfill \end{array}$$

(5)

δW_d and δW_E are variations of the work done by damping and external forces, respectively.

$$\delta W_d=-{\displaystyle {\int }_0^L}{\displaystyle {\iint }_A}\left({\sigma }_{xx}^v\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\delta {\epsilon }_{xx}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}+{\tau }_{xy}^v\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\delta {\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\right)dAdx$$

(6)

$$\delta W_E={\displaystyle {\int }_0^L}\left(f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\delta u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}+q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\delta w\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}+M\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\delta \phi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\right)dx$$

(7)

where, σ_xx(x, t), τ_xy(x, t) and ${\sigma }_{xx}^v\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$, ${\tau }_{xy}^v\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ are the elastic and viscous components of normal and shear stress; respectively. Based on the hereditary integral, Findley et al. [39], Vallala et al. [40], both elastic and viscous components of normal and shear stress can be expressed according to the Wiechert viscoelastic mechanical model as, Abdelrahman and Mahmoud [33]

$$\begin{array}{l}{\sigma }_{xx}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}E\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\epsilon }_{xx}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\text{and}\hfill \\ {\sigma }_{xx}^v\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{n=1}^N}{\displaystyle {\int }_0^t}\frac{E_n}{{\rho }_n}\left(e^{\frac{t-s}{{\rho }_n}}\right){\epsilon }_{xx}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}ds\hfill \end{array}$$

(8)

$$\begin{array}{l}{\tau }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}G\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\text{and}\hfill \\ {\tau }_{xy}^v\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{n=1}^N}{\displaystyle {\int }_0^t}\frac{G_n}{{\rho }_n}\left(e^{\frac{t-s}{{\rho }_n}}\right){\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}ds\hfill \end{array}$$

(9)

Substituting from Equations (5–9) and evaluating the variations in Equation (4), the dynamic equation of motion for viscoelastically damped Timoshenko beams can be expressed as

$$\begin{array}{l}E\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}A\frac{{\partial }^{2}u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x^2}-A{\displaystyle \sum _{n=1}^N}\left({\displaystyle {\int }_0^t}\frac{E_n}{{\rho }_n}\left(e^{\frac{t-s}{{\rho }_n}}\right)\frac{{\partial }^{2}u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x^2}ds\right)+f\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ =\rho A\frac{{\partial }^{2}u\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}}{\partial t^2}\hfill \end{array}$$

(10)

$$\begin{array}{l}\rho A\frac{{\partial }^{2}w\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}}{\partial t^2}-k_{s}AG\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\left(-\frac{\partial \varphi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x}+\frac{{\partial }^{2}w\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x^2}\right)\hfill \\ +k_{s}A{\displaystyle \sum _{n=1}^N}\left({\displaystyle {\int }_0^t}\frac{G_n}{{\rho }_n}\left(e^{\frac{t-s}{{\rho }_n}}\right)\left(-\frac{\partial \varphi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x}+\frac{{\partial }^{2}w\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x^2}\right)ds\right)\hfill \\ =q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \end{array}$$

(11)

$$\begin{array}{l}\rho I\frac{{\partial }^2\varphi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}}{\partial t^2}-IE\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\frac{{\partial }^2\varphi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x^2}\hfill \\ -k_{s}AG\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}\left(-\varphi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}+\frac{\partial w\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x}\right)\hfill \\ +I{\displaystyle \sum _{n=1}^N}\left({\displaystyle {\int }_0^t}\frac{E_n}{{\rho }_n}\left(e^{\frac{t-s}{{\rho }_n}}\right)\frac{{\partial }^2\varphi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x^2}ds\right)\hfill \\ +k_{s}A{\displaystyle \sum _{n=1}^N}\left({\displaystyle {\int }_0^t}\frac{G_n}{{\rho }_n}\left(e^{\frac{t-s}{{\rho }_n}}\right)\left(-\varphi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}+\frac{\partial w\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x}\right)ds\right)=0\hfill \end{array}$$

(12)

4 Numerical Evaluation of the Viscoelastic Constitutive Equations

To develop a computational finite element model capable of investigating the dynamic response of viscoelastically damped suspension system, the hereditary integral should be expressed in an incremental recursive form suitable for the finite element implementation, Yi and Hilton [41]. According to Wiechert viscoelastic model, the normal and shear stress components for linear isotropic materials can be expressed as

$$\begin{array}{l}{\sigma }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}E\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ -{\displaystyle \sum _{n=1}^N}{\displaystyle {\int }_0^t}\frac{E_n}{{\rho }_n}\left(e^{\frac{t-s}{{\rho }_n}}\right){\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}ds\hfill \end{array}$$

(13)

$$\begin{array}{l}{\tau }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}G\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ -{\displaystyle \sum _{n=1}^N}{\displaystyle {\int }_0^t}\frac{G_n}{{\rho }_n}\left(e^{\frac{t-s}{{\rho }_n}}\right){\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}s\mathrm{<mo\; stretchy="false">)</mo>}ds\hfill \end{array}$$

(14)

Equations (13, 14) can be written as

$$\begin{array}{l}{\sigma }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}E\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}-{\displaystyle \sum _{n=1}^N}(\frac{E_n}{{\rho }_n}{\displaystyle {\int }_0^{\Delta t}}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-s\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-s}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}ds\hfill \\ +\frac{E_n}{{\rho }_n}{\displaystyle {\int }_{\Delta t}^t}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-s\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-s}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}ds)\hfill \\ =E\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}-{\displaystyle \sum _{n=1}^N}{I}_{1n}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \end{array}$$

(15)

$$\begin{array}{l}{\tau }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}G\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}-{\displaystyle \sum _{n=1}^N}(\frac{G_n}{{\rho }_n}{\displaystyle {\int }_0^{\Delta t}}{\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-s\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-s}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}ds\hfill \\ +\frac{G_n}{{\rho }_n}{\displaystyle {\int }_{\Delta t}^t}{\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-s\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-s}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}ds)\hfill \\ =G\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}-{\displaystyle \sum _{n=1}^N}{I}_{2n}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \end{array}$$

(16)

Applying the trapezoidal rule, Equations (15, 16) can be numerically evaluated and the viscoelastic constitutive relations an incrementally expressed as

$${\sigma }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\sigma }_x^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}+{\sigma }_x^h\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\text{and}{\tau }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\tau }_{xy}^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}+{\tau }_{xy}^h\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$$

(17)

where ${\sigma }_x^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ and ${\tau }_{xy}^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ are the normal and shear stress components at the current time instant, t which can be expressed as

$${\sigma }_x^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\left(E_{\infty }+{\displaystyle \sum _{n=1}^n}{E}_n\left[\frac{{\rho }_n}{\Delta t}\left(1-e^{\frac{-\Delta t}{{\rho }_n}}\right)\right]\right){\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{E}_{ve}\epsilon \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$$

(18)

$$\begin{array}{l}{\tau }_{xy}^c\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\left(G_{\infty }+{\displaystyle \sum _{n=1}^n}{G}_n\left[\frac{{\rho }_n}{\Delta t}\left(1-e^{\frac{-\Delta t}{{\rho }_n}}\right)\right]\right){\kappa }_s{\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ ={\kappa }_s{G}_{ve}\gamma \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \end{array}$$

(19)

while ${\sigma }_x^h\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ and ${\tau }_{xy}^h\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}$ are the normal and shear stress history components which can be expressed as

$$\begin{array}{l}{\sigma }_x^h\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{n=1}^N}{E}_n(-e^{\frac{-\Delta t}{{\rho }_n}}+\frac{{\rho }_n}{\Delta t}\mathrm{<mo\; stretchy="false">(</mo>1}-e^{\frac{-\Delta t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-\Delta t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ +e^{\frac{-\Delta t}{{\rho }_n}}{h}_{\sigma n}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-\Delta t\mathrm{<mo\; stretchy="false">)</mo>})\hfill \end{array}$$

(20)

$$\begin{array}{l}{\tau }_{xy}^h\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{n=1}^N}{\kappa }_s{G}_n(-e^{\frac{-\Delta t}{{\rho }_n}}+\frac{{\rho }_n}{\Delta t}\mathrm{<mo\; stretchy="false">(</mo>1}-e^{\frac{-\Delta t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}{\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-\Delta t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ +e^{\frac{-\Delta t}{{\rho }_n}}{h}_{\tau n}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-\Delta t\mathrm{<mo\; stretchy="false">)</mo>})\hfill \end{array}$$

(21)

with h_nσ(x, t) and h_nτ(x, t) are given by

$$\begin{array}{l}{h}_{n\sigma }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{E}_n\left(\begin{array}{l}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\left[1+\frac{{\rho }_n}{\Delta t}\left(1-e^{\frac{-\Delta t}{{\rho }_n}}\right)\right]\hfill \\ +{\epsilon }_x\left(x\mathrm{,}t-\Delta t\right)\left[-e^{\frac{-\Delta t}{{\rho }_n}}+\frac{{\rho }_n}{\Delta t}\left(1-e^{\frac{-\Delta t}{{\rho }_n}}\right)\right]\hfill \end{array}\right)\hfill \\ +e^{\frac{-\Delta t}{{\rho }_n}}{h}_{\sigma n}\left(x\mathrm{,}t-\Delta t\right)\hfill \end{array}$$

(22)

$$\begin{array}{l}{h}_{n\tau }\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\kappa }_s{G}_n\left(\begin{array}{l}{\gamma }_{xy}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\left[1+\frac{{\rho }_n}{\Delta t}\left(1-e^{\frac{-\Delta t}{{\rho }_n}}\right)\right]\hfill \\ +{\gamma }_{xy}\left(x\mathrm{,}t-\Delta t\right)\left[-e^{\frac{-\Delta t}{{\rho }_n}}+\frac{{\rho }_n}{\Delta t}\left(1-e^{\frac{-\Delta t}{{\rho }_n}}\right)\right]\hfill \end{array}\right)\hfill \\ +e^{\frac{-\Delta t}{{\rho }_n}}{h}_{\tau n}\left(x\mathrm{,}t-\Delta t\right)\hfill \end{array}$$

(23)

More details for the numerical evaluation of the convolution integrals in the viscoelastic constitutive form is presented in Appendix A.

5 Finite Element Model

To derive the incremental finite element equations of motion, the double wishbone suspension links are divided into finite elements. Both two nodded and three nodded isoparametric Timoshenko beam elements are considered, Capsoni et al. [42]. On the other hand, the revolute joint element developed in ANSYS Kohnke [43] is adopted to model the joint flexibility of the suspension system. Based on TBT, the dynamic finite element equations of motion is established by applying the virtual work principle, taking into account d’Alembert’s principle and considering viscoelastic damping forces as

$$\begin{array}{l}{\displaystyle {\int }_v}\delta {\epsilon }_x^T{\sigma }_{x}dv+{\displaystyle {\int }_v}\delta {\gamma }_{xz}^T{\tau }_{xz}dv\hfill \\ =-{\displaystyle {\int }_v}\rho \delta U^T\ddot{U}dv+{\displaystyle {\int }_0^L}\left[\delta u_0^T{P}_x+\delta w^T{P}_z+\delta {\varphi }^{T}m\right]dx\hfill \end{array}$$

(24)

where P_x, P_z, and m refer to the distributed normal, shear and bending moment while ρ is the mass density. Equation (24) can be rewritten as

$$\begin{array}{l}{\displaystyle {\int }_0^L}{\displaystyle {\int }_A}\left[\delta {\epsilon }_x^T{\sigma }_x+\delta {\gamma }_{xz}^T{\tau }_{xz}\right]dAdx\hfill \\ =-{\displaystyle {\int }_0^L}{\displaystyle {\int }_A}\rho \delta U^T\ddot{U}dAdx\hfill \\ +{\displaystyle {\int }_0^L}\left[\delta u_0^T{P}_x+\delta w^T{P}_z+\delta {\varphi }^{T}m\right]dx\hfill \end{array}$$

(25)

The displacement, strain displacement relation and the constitutive relations for linear viscoelastic materials can be expressed in matrix form as, Oñate [44]

$$\begin{array}{l}{u}^e=\left[N_1{N}_1{N}_1{N}_2{N}_2{N}_2{N}_3{N}_3{N}_3\right]{\left[u_i{w}_i{\varphi }_i{u}_j{w}_j{\varphi }_j{u}_k{w}_k{\varphi }_k\right]}^T\hfill \\ =N^e{q}^e\hfill \end{array}$$

(26a)

$$\begin{array}{l}\left\{\begin{array}{c}{\epsilon }_x\\ {\kappa }_x\\ {\gamma }_{xy}\end{array}\right\}=\left(\frac{\partial \zeta }{\partial x}\right)\hfill \\ \times \left[\begin{array}{ccccccccc}\frac{\partial N_1}{\partial \zeta }& 0& 0& \frac{\partial N_2}{\partial \zeta }& 0& 0& \frac{\partial N_3}{\partial \zeta }& 0& 0\\ 0& \frac{\partial N_1}{\partial \zeta }& 0& 0& \frac{\partial N_2}{\partial \zeta }& 0& 0& \frac{\partial N_3}{\partial \zeta }& 0\\ 0& \frac{\partial N_1}{\partial \zeta }& N_1& 0& \frac{\partial N_2}{\partial \zeta }& N_2& 0& \frac{\partial N_3}{\partial \zeta }& N_3\end{array}\right]\hfill \\ \times \mathrm{<mo\; stretchy="false">[</mo>}{q}^e\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>}{B}^e{q}^e\hfill \end{array}$$

(26b)

$$N_1=\frac{\zeta \mathrm{<mo\; stretchy="false">(</mo>}\zeta -\mathrm{1<mo\; stretchy="false">)</mo>}}{2}\mathrm{,}{N}_2=\mathrm{<mo\; stretchy="false">(</mo>1}-{\zeta }^2\mathrm{<mo\; stretchy="false">)</mo>,}{N}_3=\frac{\zeta \mathrm{<mo\; stretchy="false">(</mo>}\zeta +\mathrm{1<mo\; stretchy="false">)</mo>}}{2}$$

(26c)

where ε_x, κ_x, γ_xy are the normal strain, the curvature and the shear strain components; respectively, N^e is the element shape function matrix and B^e is the strain displacement relation matrix for the three noded Timoshenko beam (TB) element, shown in Figure 3.

$$\mathit{\sigma }={\mathit{\sigma }}^{\mathbf{c}}+{\mathit{\sigma }}^{\mathbf{h}}\equiv {\mathbf{E}}_{\mathbf{v}}\mathit{\epsilon }+{\mathit{\sigma }}^{\mathbf{h}}$$

(27a)

Figure 3 Three noded isoparametric TB element.

$$\begin{array}{l}{\mathbf{E}}_{\mathbf{v}}=\left(E_{\infty }+{\displaystyle \sum _{n=1}^n}{E}_n\left[\frac{{\rho }_n}{\Delta t}\left(1-e^{\frac{-\Delta t}{{\rho }_n}}\right)\right]\right)\hfill \\ \times \left[\begin{array}{ccc}A& 0& 0\\ 0& I& 0\\ 0& 0& \frac{k_{s}A}{2\left(1+v\right)}\end{array}\right]\text{}\mathrm{,}\hfill \\ \epsilon =\left\{\begin{array}{c}\frac{\partial u}{\partial x}\\ \kappa \\ \gamma \end{array}\right\}\hfill \end{array}$$

(27b)

where A, k_s, ν and I are the beam element cross sectional area, the TB shear correction factor, Poisson’s ratio, and the beam second moment of area; respectively. Substituting Equations (26), (27) into Equation (25), the dynamic finite element equation of motion for viscoelastic system can be written in global coordinates as

$${\mathbf{R}}^T{\mathbf{M}}_{\mathbf{l}}^{\text{e}}\mathrm{<mi>\; </mi>}\mathbf{R}\ddot{\mathit{q}}+{\mathbf{R}}^T\mathrm{<mi>\; </mi>}{\mathbf{K}}_{\mathbf{l}}^{\text{e}}\mathrm{<mi>\; </mi>}\mathbf{R}\mathit{q}={\mathbf{R}}^T\mathrm{<mo\; stretchy="false">(</mo>}{\mathit{F}}_{ext}^e-{\mathit{F}}_{ve}^e\mathrm{<mo\; stretchy="false">)</mo>}$$

(28)

By the same way, the dynamic finite element equation of motion for viscously damped TBT can be written in global coordinates as

$${\mathbf{R}}^T{\mathbf{M}}_{\mathbf{l}}^{\text{e}}\mathrm{<mi>\; </mi>}\mathbf{R}\ddot{\mathit{q}}+{\mathbf{R}}^T\mathrm{<mi>\; </mi>}{\mathbf{C}}_{\mathbf{l}}^{\text{e}}\mathbf{R}\dot{\mathit{q}}+{\mathbf{R}}^T\mathrm{<mi>\; </mi>}{\mathbf{K}}_{\mathbf{l}}^{\text{e}}\mathbf{R}\mathit{q}={\mathbf{R}}^T{\mathit{F}}_{ext}^e$$

(29)

where ${\mathbf{M}}_{\mathbf{l}}^{\mathbf{e}}$, ${\mathbf{K}}_{\mathbf{l}}^{\mathbf{e}}$, ${\mathbf{C}}_{\mathbf{l}}^{\mathbf{e}}$, R and ${\mathit{F}}_{\mathit{ve}}^e$ are the element mass, stiffness, and viscous damping, transformation matrices and the force vector due to viscoelastic material history; respectively.

$$\begin{array}{l}{\mathbf{K}}_{\mathbf{l}}^{\mathbf{e}}={\displaystyle {\int }_{-1}^1}{\mathbf{B}}^T{\mathbf{E}}_{\mathbf{v}}\mathbf{B}\left|\mathbf{J}\right|d\xi ,{\mathit{F}}_{ve}^e={\displaystyle {\int }_v}{\mathbf{B}}^T{\mathit{\sigma }}^{\text{h}}dv,\hfill \\ {\mathbf{C}}_{\mathbf{l}}^{\mathbf{e}}={\displaystyle {\int }_{-1}^1}{C}_v{\mathbf{N}}^{\mathbf{T}}\mathbf{N}\left|\mathbf{J}\right|d\xi \text{and}{\mathbf{M}}_{\mathbf{l}}^{\mathbf{e}}={\displaystyle {\int }_{-1}^1}\rho {\mathbf{N}}^{\mathbf{T}}\mathbf{N}\left|\mathbf{J}\right|d\xi \hfill \end{array}$$

(30)

Considering the revolute joint element, COMBIN7 as expressed in ANSYS, the dynamic equation of motion can be expressed in matrix form as

$$\begin{array}{l}\left[\begin{array}{cccccc}\frac{M}{2}& 0& 0& 0& 0& 0\\ 0& \frac{M}{2}& 0& 0& 0& 0\\ 0& 0& \frac{I}{2}& 0& 0& 0\\ 0& 0& 0& \frac{M}{2}& 0& 0\\ 0& 0& 0& 0& \frac{M}{2}& 0\\ 0& 0& 0& 0& 0& \frac{I}{2}\end{array}\right]\{{\ddot{q}}^e\}\hfill \\ +\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& c_v& 0& 0& -c_v\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& -c_v& 0& 0& c_v\end{array}\right]\{{\dot{q}}^e\}\hfill \\ +\left[\begin{array}{cccccc}{k}_x& 0& 0& -k_x& 0& 0\\ 0& k_y& 0& 0& -k_y& 0\\ 0& 0& k_t& 0& 0& -k_t\\ -k_x& 0& 0& k_x& 0& 0\\ 0& -k_y& 0& 0& k_y& 0\\ 0& 0& -k_t& 0& 0& k_t\end{array}\right]\{q^e\}\hfill \end{array}$$

(31)

where k_x and k_y are dependent on the geometry and material characteristics of the revolute joint while k_t is dependent on the generated frictional torque during the relative motion between the two mating surfaces. While c_v is the viscous damping coefficient of the revolute joint which depends on the type and the extent of the contact area between the two mating surfaces in the joint. Finally, m and I are the total mass and mass moment of inertia of the revolute joint; respectively.

6 Validation of the Developed Model

Within this section the capability of the developed model to analyze the dynamic behavior of elastic and/or viscoelastic structures under different excitation patterns is verified. Consider a simply supported beam under uniformly distributed load with intensity of q₀, as shown in Figure 4. To check the validity of the developed model for joint flexibility, both hinged and roller supports can be represented by flexible translational springs, as shown in Figure 4(b). Froio and Rizzi [45] and Da Costa Azevêdo [46]. Higher values of these spring stiffness leads to fully rigid support consequently the same response should be detected for both rigid and flexible supports with higher stiffness values. Moreover, both EBBT and TBT are considered during the verification procedure. The parameters used for numerical analysis are listed as follows: The load intensity q₀ = 100 N/m. The beam has a length of L = 4 m, width b = 0.2 m, and height, h = 0.4 m. The modulus of elasticity of elastic material, mass density, shear correction factor, and Poisson’s ratio are taken to be E = 9.8 × 10⁷ Pa, ρ = 500 kg/m³, k_s = 5/6, and v = 0.3, respectively. The viscoelastic material parameters for a single repeated model are, the fully relaxed modulus, E₈ = 1.96 × 10⁷ Pa, the transient modulus E₁ = 7.84 × 10⁷ Pa, and the material relaxation time ρ_v = 2.24 s while the instantaneous creep compliance D₀ = 1.2041 × 10⁻⁸ Pa⁻¹, the transient creep compliance D₁ = 4.0816 × 10⁻⁸ Pa⁻¹, and the viscoelastic retardation time τ = 11.2 s. Beams are discretized into beam elements. The elastodynamic and viscoelastodynamic transient responses are detected by applying the Newmark’s numerical integration technique.

The analytical quasi-static viscoelastic solution based on the elastic viscoelastic correspondence principle can be expressed for both TBT and EBBT as, Vallala et al. [40]

$$W_{max}=\{\begin{array}{ll}\left(\frac{5q{L}^4}{384I}+\frac{5qL\left(1+v\right)}{8{\kappa }_{s}A}\right)D\left(t\right)\hfill & TBT\hfill \\ \frac{5q{L}^4}{384I}D\left(t\right)\hfill & EBBT\hfill \end{array}$$

(32)

Figure 4 Geometry of the simply-supported and elastically supported beams.

where D(t) is the viscoelastic creep compliance function which can be defined as

$$D\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{D}_0+{\displaystyle \sum _{n=1}^N}{D}_n\left(1-exp\left(\frac{-t}{{\tau }_n}\right)\right)$$

(33)

with D₀ is the equivalent elastic compliance of the viscoelastic material and D_n is the transient component of the viscoelastic material. On the other hand, the semi analytical elastodynamic solution for both TBT and EBBT can be written as, Abdelrahman et al. [47]

$$W_n\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{j=1}^n}{\beta }_j\sin \left(\frac{j\pi x}{L}\right)\mathrm{<mo\; stretchy="false">(</mo>1}-cos\mathrm{<mo\; stretchy="false">(</mo>}{\alpha }_{j}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$$

(34)

where

$$\begin{array}{l}{\beta }_j=\left\{\frac{\left(\frac{2q}{j\pi }\right)}{{\left(\frac{j\pi }{l}\right)}^{4}EI}\right\}\mathrm{<mo\; stretchy="false">(</mo>1}-cos\mathrm{<mo\; stretchy="false">(</mo>}j\pi \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}\text{and}\hfill \\ {\alpha }_j=\{\begin{array}{cc}{\left(\frac{{\left(\frac{j\pi }{l}\right)}^4\mathrm{<mo\; stretchy="false">(</mo>}EI\mathrm{<mo\; stretchy="false">)</mo>}}{\rho A\left[1+\left(\frac{EI}{k_{s}GA}\right){\left(\frac{j\pi }{l}\right)}^2\right]}\right)}^{\mathrm{1<mo>/</mo>2}}& \mathrm{<mo\; stretchy="false">(</mo>}TBT\mathrm{<mo\; stretchy="false">)</mo>}\\ {\left(\frac{{\left(\frac{j\pi }{l}\right)}^4\mathrm{<mo\; stretchy="false">(</mo>}EI\mathrm{<mo\; stretchy="false">)</mo>}}{\left[\rho A+\rho I{\left(\frac{j\pi }{l}\right)}^2\right]}\right)}^{\mathrm{1<mo>/</mo>2}}& \mathrm{<mo\; stretchy="false">(</mo>}EBBT\mathrm{<mo\; stretchy="false">)</mo>}\end{array}\hfill \end{array}$$

(35)

Moreover the semi analytical viscoelastodynamic deflection can be obtained using the mixed Galerkin’s Laplace techniques for EBBT as, Abdelrahman et al. [47]; Martin [37]

$$W_n\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{j=1}^n}\sqrt{\frac{2}{L}}\sin \left(\frac{j\pi x}{L}\right)T_j\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$$

(36)

where

$$\begin{array}{l}{T}_j\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{c}_{1j}+c_{2j}\exp \mathrm{<mo\; stretchy="false">(</mo>}-0.0893t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ +c_{3j}\exp \mathrm{<mo\; stretchy="false">(</mo>}-0.1786t\mathrm{<mo\; stretchy="false">)</mo>}\left\{\cos \mathrm{<mo\; stretchy="false">(</mo>}{\omega }_{j}t\mathrm{<mo\; stretchy="false">)</mo>}-\frac{0.1786}{{\omega }_j}\sin \mathrm{<mo\; stretchy="false">(</mo>}{\omega }_{j}t\mathrm{<mo\; stretchy="false">)</mo>}\right\}\hfill \\ +\frac{c_{4j}}{{\omega }_j}\exp \mathrm{<mo\; stretchy="false">(</mo>}-0.1786t\mathrm{<mo\; stretchy="false">)</mo>}\sin \mathrm{<mo\; stretchy="false">(</mo>}{\omega }_{j}t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \end{array}$$

(37)

where

$$\begin{array}{l}{\omega }_j=\sqrt{994.3845\times {\mathrm{<mo\; stretchy="false">(</mo>}j\mathrm{<mo\; stretchy="false">)</mo>}}^4-{\mathrm{<mo\; stretchy="false">(</mo>0.1786<mo\; stretchy="false">)</mo>}}^2}\mathrm{,}{c}_{1j}=\frac{P_j\times 0.44643}{88.798533\times {\mathrm{<mo\; stretchy="false">(</mo>}j\mathrm{<mo\; stretchy="false">)</mo>}}^4}\mathrm{,}\hfill \\ P_j=\frac{2.2508}{j}\mathrm{<mo\; stretchy="false">[</mo><mo\; stretchy="false">(</mo>}-{\mathrm{1<mo\; stretchy="false">)</mo>}}^{j+1}+\mathrm{1<mo\; stretchy="false">]</mo>}\hfill \end{array}$$

(38)

$$\begin{array}{l}{c}_{2j}=\frac{P_j\times 0.35713}{0.002136-88.78433\times {\left(j\right)}^4}\mathrm{,}{c}_{3j}=-\mathrm{<mo\; stretchy="false">(</mo>}{c}_{1j}+c_{2j}\mathrm{<mo\; stretchy="false">)</mo>,}\hfill \\ c_{4j}=-\{0.44643\times c_{1j}+0.35713\times c_{2j}+0.0893\times c_{3j}\}\hfill \end{array}$$

(39)

The maximum transverse deflection is investigated for both numerical and analytical analyses over a time interval [0, 40] s for both EBBT and TBT. The dependency of the transverse deflection on time for both viscoelastic and elastic materials behavior for quasistatic and dynamic analyses for both simply supported rigid support (S-S) and restrained supported beam with large values of the spring stiffness (flex. support) are illustrated in Figure 5. It is noticed that an excellent agreement is obtained between both numerical and analytical solutions. Moreover, at very high values of the support spring stiffness, both S-S rigid support and flexible support analyses result in the same response for both elastic and viscoelastic behavior.

Figure 5 Variation of the maximum deflection with time for different materials behavior for both EBBT and TBT.

7 Computational Results

The applicability of the developed model to study and analyze the dynamic behavior of the double wishbone suspension system under trapezoidal road bump profile and at different values of travelling velocity is demonstrated in this section. Moreover, to study the material damping effect on the transient response, both undamped, proportional, and viscoelastic damped responses are investigated. The maximum transverse deflection and the resultant stresses at the critically loaded point on the lower control arm of the double wishbone suspension system are investigated. Exploiting the symmetry of the vehicle, only a quarter car model is simulated. The wheel center is subjected to the different excitations based on the road bump profile. The comprehensive quarter car model incorporating double wishbone linkage shown in Figure 2. is presented for the discussion including the tire and strut as flexible translational springs while the chassis and the tire axle are considered rigid compared to the lower and upper control arms. The model geometrical parameters and material characteristics used to obtain the computational results are summarized in Table 1.

The road disturbance as the system excitation results in vertical vibrations in the whole system. Trapezoidal road disturbance profile is considered to investigate the effect of road irregularity (bump) profile on the dynamic response of the suspension systems. Variation of the trapezoidal road bump profile with time at different vehicle velocities is illustrated in Figure 6, Hassaan [46].

7.1 Effect of Vehicle Velocity on the Transient Response of the Lower Control Arm

The vehicle velocity during passing over the road bump greatly affects the transient response of the vehicle suspension system. The transient response of the lower control arm is investigated for the different vehicle velocities; 10, 15, 25, 35, and 45 km/hr for different material behaviors. Variations of the maximum normalized transverse deflection (W/L) at the critically loaded point on the lower control arm with the normalized time (T/τ), (τ = ρ_v) at different travelling velocities for elastic, proportional viscously damped and viscoelastic damped behaviors are depicted in Figure 7. It is noticed that, the time needed to reach the peak value of the transient response and to drop to the minimum amplitude greatly affected by the vehicle velocity. It is also noticed that for the maximum transverse deflection, although the peak amplitude of vibration is detected at the maximum value of the vehicle velocity (v = 45 km/hr) the smallest amplitude is obtained at a vehicle velocity of 25 km/kr. Thus the smallest vibration amplitude is not obtained at the smallest value of the vehicle velocity. On contrary to the transient elastic response which proceeds with a constant amplitude for all vehicle velocities, the transient proportional viscously damped response is decaying rapidly due to the constant damping ratio while the viscoelastic transient response is continuously decaying gradually with time due to the algorithmic viscoelastic damping. The decaying rate and consequently the time required to reach the steady state oscillations for the viscoelastic transient response depends on the viscoelastic relaxation time.

Table 1 The physical properties of the suspension elastic model

Parameter	Value
E (GPa), (Elastic suspension links)	210
E∞ (GPa), (VE suspension links)	42
E1 (GPa), (VE suspension links)	168
ρv (s), Relaxation time (VE suspension links)	1.0
E (GPa), (chassis and tire axle)	1050
ρ, kg/m3, (chassis and suspension links)	7830
Ks (N/mm), (Strut stiffness)	340
Ktire (N/mm), (Tire stiffness)	235
Kx (N/mm), (flex. Joint)	17.544 × 107
Ky (N/mm), (flex. Joint)	17.544 × 107
Kt (N/mm), (flex. Joint)	19.737 × 104
A m2	6 × 10−4
L (Lower and upper arm length)	0.4 m
I m4	12 × 10−8
Cv (N.s/m)	4 × 104
ξ(—)	75 × 10−4

Figure 6 Variation of road profile with time at different vehicle velocities.

Figure 7 Variation of the maximum transverse deflection with time for different material behaviors at different vehicle velocities.

The dependency of the maximum normalized normal force (N/τ) on the normalized time (T/τ) for different material behaviors is illustrated in Figure 8. It is noticed that, as detected for the maximum normalized transverse deflection, the preferable transient response for trapezoidal road bump profile is detected at vehicle velocity of 25 km/hr. It is noticed that, higher and constant steady state oscillation amplitude is detected for the transient elastic response compared with both proportional viscously damped and viscoelastically damped response. Also, the proportional viscously damped response needs short time interval to decay to smaller oscillation amplitude due to the constant damping ratio. On the other hand, depending on the viscoelastic relaxation time, long time interval is needed for the viscoelastic transient response to reach the corresponding smaller oscillation amplitude.

The maximum normalized shear force, (V/EA) versus the normalized time (T/τ) for elastic, proportional viscously damped, and viscoelastically damped responses is depicted in Figure 9. It is noticed that, smaller values of the instantaneous and steady state oscillation amplitudes of the maximum normalized shear force are detected compared with the obtained maximum normal force for all material behaviors. It is also noticed that preferable response is detected at vehicle velocity of 25 km/hr.

Figure 8 Variation of the maximum normal force with time for different material behaviors at different vehicle velocities.

Figure 9 Variation of the maximum shear force with time for different material behaviors at different vehicle velocities.

Figure 10 Variation of the maximum bending moment with time for different material behaviors at different vehicle velocities.

Variation of the maximum normalized bending moment with the normalized time for different material responses is depicted in Figure 10. It is noticed that, due to material damping, both proportional viscously and viscoelastically damped response possesses stress relaxation leading to decayed oscillation amplitudes compared with the corresponding transient elastic response. The decaying rate depends on the both the viscous damping ratio and the viscoelastic material relaxation time.

7.2 Effect of the Viscoelastic Relaxation Time on the Transient Response of the Lower Control Arm

The relaxation time of the viscoelastic control arm greatly affects the decaying rate of the transient response of the double wishbone suspension system. The transient response of the lower control arm is investigated for different values of the relaxation time (ρ_v); 0.25 ρ_v, 0.5 ρ_v, ρ_v, 1.5 ρ_v, and 2 ρ_v. The dependency of the maximum normalized transverse deflection, the normal force, the shear force, and the bending moment on the normalized time at the critically loaded point of the viscoelastic lower control arm at different values of the relaxation time is illustrated in Figures 11 and 12, respectively. It is noticeable that the same instantaneous response is detected for all values of the relaxation time. On the other hand the decaying rate of the transient response oscillation amplitudes are increased with decreasing the viscoelastic relaxation time and consequently, the time required to reach the steady state amplitude is decreased with decreasing the relaxation time.

Figure 11 Variation of the maximum transverse deflection and the maximum normal force for different values of the relaxation time.

Figure 12 Variation of the shear force and the maximum bending moment at for different values of the relaxation time.

8 Conclusion

An incremental finite element model capable of investigating the dynamic behavior of double wishbone vehicle suspension system is proposed considering both links and joints flexibilities as well as the internal vis-coelastic material damping effects. Besides the internal viscoelastic damping proportional damping mechanisms is considered to show the effect of different damping techniques on the dynamic behavior of the double wishbone suspension systems. The TBT is adopted to model the control arm flexibility. On the other hand, the revolute joint element developed in ANSYS is adopted to model the flexibility of bushings. An incremental viscoelastic constitutive relation suitable for finite element procedure is developed. The proposed procedure is verified by comparing the obtained results with the available analytical forms and an excellent agreement is detected. Parametric studies are conducted to study the effects of the travelling velocity and the viscoelastic relaxation time on the transient response of the suspension lower control arm. The obtained results shows that

1. The vehicle travelling velocity over the road bump affects both the instantaneous and the steady state response of the suspension system. Smaller oscillation amplitude of the transient response does not correspond to the smallest vehicle travelling velocity over the road bump. The most preferable travelling velocity over the road bump depends on the bump profile and the material characteristics.
2. The transient elastic response possesses higher values of the steady state amplitudes compared with the corresponding proportional viscously and viscoelastically damped responses due to the absence of damping.
3. The proportional viscously damped response exhibits rapidly decaying oscillation amplitudes for both the transverse deflection and the resultant stress due to the constant damping factor which is determined based on the real normal modal analysis of undamped system. Unfortunately, most real physical structures possess complex modes instead of real normal modes. Thus proportional viscous damping model does not simulate the actual behavior of the real physical materials.
4. Depending on the viscoelastic relaxation time, the viscoelastically damped response produces slowly decaying oscillation amplitude compared with the proportional viscously damping mechanism. Thus longer time interval is needed to reach a steady state oscillation amplitudes.
5. The viscoelastic relaxation time has a great effect on the transient response of the vehicle suspension lower control arm. Increasing the viscoelastic relaxation time results in decreasing the value of the steady state amplitude and the time required to reach it.
6. The viscoelastic damping model is more suitable for modelling the damping mechanism in solids and structures. This is because, all the viscoelastic material damping parameters are determined experimentally through either creep or relaxation tests not on modal analysis bases. Thus simulating the actual behavior of real life solids and structures.

Appendix A
A Numerical Approximation for the Convolution Integral

$$\begin{array}{l}{\displaystyle \sum _{n=1}^N}{I}_{1n}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \sum _{n=1}^N}(\frac{E_n}{{\rho }_n}{\displaystyle {\int }_0^{\Delta t}}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-s\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-s}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}ds\hfill \\ +\frac{E_n}{{\rho }_n}{\displaystyle {\int }_{\Delta t}^t}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-s\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-s}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}ds)\hfill \end{array}$$

(A1)

The first term in Equation (A1) can be simplified using the trapezoidal rule as

$$\begin{array}{l}\frac{E_n}{{\rho }_n}{\displaystyle {\int }_0^{\Delta t}}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-s\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-\mathrm{<mo\; stretchy="false">(</mo>}s\mathrm{<mo\; stretchy="false">)</mo>}}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}ds=-\frac{E_n}{{\rho }_n}{\displaystyle {\int }_t^{t-\Delta t}}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\eta \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{\eta -t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}d\eta \hfill \\ =\frac{E_n}{{\rho }_n}\mathrm{<mo\; stretchy="false">(</mo>}{e}^{\frac{-t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle {\int }_{t-\Delta t}^t}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}\eta \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{\eta }{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}d\eta \hfill \end{array}$$

(A2)

Evaluating the integral in Equation (A2) yields;

$$\frac{E_n}{{\rho }_n}\left(e^{\frac{-t}{{\rho }_n}}\right)\left[\left({\epsilon }_x\left(x\mathrm{,}\eta \right){\rho }_n\left(e^{\frac{\eta }{{\rho }_n}}\right)\right)|\begin{array}{c}t\\ t-\Delta t\end{array}-{\displaystyle {\int }_{t-\Delta t}^t}\frac{\partial {\epsilon }_x\left(x\mathrm{,}\eta \right)}{\partial \eta }{\rho }_n\left(e^{\frac{\eta }{{\rho }_n}}\right)d\eta \right]$$

(A3)

Simplifying Equation (A3) yields;

$$\begin{array}{l}\frac{E_n}{{\rho }_n}\mathrm{<mo\; stretchy="false">(</mo>}{e}^{\frac{-t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}({\rho }_n\left[\left({\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}\right)-{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}t-\Delta t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{t-\Delta t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}\right]\hfill \\ -{\displaystyle {\int }_{t-\Delta t}^t}\frac{\partial {\epsilon }_x\left(x\mathrm{,}\eta \right)}{\partial \eta }{\rho }_n\mathrm{<mo\; stretchy="false">(</mo>}{e}^{\frac{\eta }{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}d\eta )\hfill \end{array}$$

(A4)

Evaluating the integral and simplifying Equation (A4) yields

$$\begin{array}{l}\frac{E_n}{{\rho }_n}({\rho }_n\mathrm{<mo\; stretchy="false">[</mo>}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}-{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-\Delta t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-\Delta t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}\hfill \\ -{\dot{\epsilon }}_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle {\int }_{t-\Delta t}^t}{\rho }_n\mathrm{<mo\; stretchy="false">(</mo>}{e}^{\frac{\eta }{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}d\eta )\hfill \\ \cong \frac{E_n}{{\rho }_n}({\rho }_n\mathrm{<mo\; stretchy="false">[</mo>}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}-{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-\Delta t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-\Delta t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}\hfill \\ -{\rho }_n^2{\dot{\epsilon }}_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{\eta }{{\rho }_n}}{\mathrm{<mo\; stretchy="false">)</mo><mo>|</mo>}}_{t-\Delta t}^t)\hfill \\ \cong \frac{E_n}{{\rho }_n}({\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\left[{\rho }_n-\frac{{\rho }_n^2}{\Delta t}\right]\mathrm{<mo\; stretchy="false">(</mo>1}-e^{\frac{-\Delta t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ +{\epsilon }_x\left(x\mathrm{,}t-\Delta t\right)\left[-{\rho }_n{e}^{\frac{-\Delta t}{{\rho }_n}}+\frac{{\rho }_n^2}{\Delta t}\mathrm{<mo\; stretchy="false">(</mo>1}-e^{\frac{-\Delta t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}\right]\text{})\hfill \end{array}$$

(A5)

The second integral in Equation (A1) can be simplified as

$$\begin{array}{l}\frac{E_n}{{\rho }_n}{\displaystyle {\int }_{\Delta t}^t}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-s\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-s}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}ds\hfill \\ =\frac{E_n}{{\rho }_n}{\displaystyle {\int }_0^{t-\Delta t}}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-\Delta t-u\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>}{e}^{\frac{-u-\Delta t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}du\hfill \\ =\mathrm{<mo\; stretchy="false">(</mo>}{e}^{\frac{-\Delta t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}{h}_{n\mathit{\sigma }}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-\Delta t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \end{array}$$

(A6)

Consequently, the viscoelastic constitutive relation can be expressed as

$$\begin{array}{l}{\mathit{\sigma }}_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}E\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ -{\displaystyle \sum _{n=1}^n}[\frac{E_n}{{\rho }_n}\left(\begin{array}{l}{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}\left[{\rho }_n+\frac{{\rho }_n^2}{\Delta t}\mathrm{<mo\; stretchy="false">(</mo>1}-e^{\frac{-\Delta t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}\right]\hfill \\ +{\epsilon }_x\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{,}t-\Delta t\mathrm{<mo\; stretchy="false">)</mo>}\left[-{\rho }_n{e}^{\frac{-\Delta t}{{\rho }_n}}-\frac{{\rho }_n^2}{\Delta t}\mathrm{<mo\; stretchy="false">(</mo>1}-e^{\frac{-\Delta t}{{\rho }_n}}\mathrm{<mo\; stretchy="false">)</mo>}\right]\hfill \end{array}\right)\hfill \\ +e^{\frac{-\Delta t}{{\rho }_n}}{I}_{1n}\left(x\mathrm{,}t-\Delta t\right)]\hfill \end{array}$$

(A7)

References

[1] M. Jonsson, Simulation of dynamical behaviour of a front wheel suspension. Vehicle Syst Dyn, 20(5), (1991) 269–281.

[2] G. Rill, Road vehicle dynamics: fundamentals and modeling. Crc Press, 2011.

[3] M. Kamal and M. M. Rahman, Study on dynamic behaviour of wishbone suspension system. In IOP Conference Series: Materials Science and Engineering, 36(1), (2012), 012019, IOP Publishing.

[4] L. Zhao, S. Zheng, J. Feng and Q. Hong, Dynamic Structure Optimization Design of Lower Control Arm Based on ESL. Research J App Sci, Eng Tech, 4(22), (2012) 4871–4878.

[5] A. M. Patil, A. S. Todkar, R. S. Mithari and V. V. Patil, Experimental and Finite Element Analysis of Left Side Lower Wishbone Arm of Independent Suspension System. IOSR J Mech Civil Eng, 2013 (2013).

[6] H. D. Banitalebi, Vibration and force analysis of lower arm of suspension system, Doctoral dissertation, Université du Québec á Chicoutimi, 2014.

[7] Y. A. Papegay, J. P. Merlet and D. Daney, Exact kinematics analysis of Cars suspension mechanisms using symbolic computation and interval analysis. Mechanism and Machine Theory, 40(2005) 395–413.

[8] R. Sancibrian, P. Garcia, F. Viadero, A. Fernandez and A. De-Juan, Kinematic design of double-wishbone suspension systems using a multi-objective optimization approach. Vehicle System Dynamics: International Journal of Vehicle Mechanics and Mobility, 48(7), (2010) 793–813.

[9] K. V. Reddy, M. Kodati, K. Chatra and S. Bandyopadhyay, A comprehensive kinematic analysis of the double wishbone and MacPherson strut suspension systems. Mechanism and Machine Theory, 105 (2016) 441–470.

[10] B. Németh and P. Gáspár, Nonlinear analysis and control of a variable-geometry suspension system. Control Engineering Practice, 61 (2017) 279–291.

[11] C. Kavitha, S. A. Shankar, K. Karthika, B. Ashok and S. D. Ashok, Active camber and toe control strategy for the double wishbone suspension system. Journal of King Saud University-Engineering Sciences, (2018).

[12] A. H. Attia, Dynamic modelling of the double wishbone motor-vehicle suspension system. Eur J Mech A/Solids, 21 (2002) 167–174.

[13] C. M. Ramírez, M. Tomás-Rodríguez and S. A. Evangelou. Dynamic analysis of double wishbone front suspension systems on sport motorcycles. Nonlinear Dynamics, 91(4) (2018), 2347–2368.

[14] R. Burdzik and B. Łazarz, Analysis of properties of automotive vehicle suspension arm depending on different materials used in the MSC. Adams environment. Archives of Materials Science and Engineering, 58(2), (2012) 171–176.

[15] C. O. Ijagbemi, B. I. Oladapo, H. M. Campbell and C. O. Ijagbemi, Design and simulation of fatigue analysis for a vehicle suspension system (VSS) and its effect on global warming. Procedia Engineering, 159, (2016) 124–132.

[16] H. Kanchwala and A. Chatterjee. A generalized quarter car modelling approach with frame flexibility and other nonlocal effects. Sâdhanâ, 42(7), (2017), 1175–1192.

[17] U. A. Ramesh, D. Siddhartha, R. A. Sethu, N. Ramanan and V. Karthi, Modeling and analysis of lower wishbone for independent suspension system for commercial vehicles. Journal of Chemical and Pharmaceutical Sciences, 2 (2017), 166–168.

[18] X. G. Song, J. H. Jung, H. J. Son, J. H. Park, K. H. Lee and Y. C. Park, (2010). Metamodel-based optimization of a control arm considering strength and durability performance. Computers & Mathematics with Applications, 60(4), 976–980.

[19] Y. D. Kim, J. E. Jeong, J. S. Park, I. H. Yang, T. S. Park, P. B. Muhamad and J. E. Oh, Optimization of the lower arm of a vehicle suspension system for road noise reduction by sensitivity analysis. Mechanism and Machine Theory, 69(2013) 278–302.

[20] P. Baranowski, J. Malachowski and L. Mazurkiewicz, Numerical and experimental testing of vehicle tyre under impulse loading conditions. Int J Mech Sci, 106 (2016) 346–356.

[21] H. Sayyaadi and N. Shokouhi, A new model in rail–vehicles dynamics considering nonlinear suspension components behavior. International Journal of Mechanical Sciences, 51(3), (2009). 222–232.

[22] J. Thomas and B. A. H. Abbas. Finite element model for dynamic analysis of Timoshenko beam. Journal of Sound and Vibration 41, no. 3 (1975): 291–299.

[23] J. R. Banerjee and F. W. Williams. Coupled bending-torsional dynamic stiffness matrix of an axially loaded Timoshenko beam element. International Journal of Solids and Structures 31, no. 6 (1994): 749–762.

[24] T. X. Wu and D. J. Thompson. A double Timoshenko beam model for vertical vibration analysis of railway track at high frequencies. Journal of Sound and Vibration 224, no. 2 (1999): 329–348.

[25] M. Krawczuk, M. Palacz and W. Ostachowicz. The dynamic analysis of a cracked Timoshenko beam by the spectral element method. Journal of Sound and Vibration 264, no. 5 (2003): 1139–1153.

[26] N. Navadeh, R. W. Hewson and A. S. Fallah. Dynamics of transversally vibrating non-prismatic Timoshenko cantilever beams. Engineering Structures 166 (2018): 511–525.

[27] D. J. Inman, Engineering vibration. Pearson: Pearson Education; (4th ed.), (2014).

[28] N. Nakamura, Time history response analysis using extended Rayleigh damping model. Procedia Engineering, 199, 1472–1477, (2017).

[29] R. D. Cook, D. S. Malkus and M. E. Plesha, “Stress Stiffening and Buckling”, in Concepts and Applications of Finite Element Analysis, (John Wiley & Sons, New York), pp. 429–448 (1989).

[30] A. K. Chopra, “Damping in structures”, in Dynamics of Structures, (Prentice-Hall Inc., Upper Saddle River, New Jersey), pp. 409–429, (1995).

[31] R. Lakes and R. S. Lakes, Viscoelastic Materials. Cambridge university press, (2009).

[32] I. H. Shames, Elastic and Inelastic Stress Analysis. CRC Press, (1997).

[33] A. A. Abdel Rahman and F. F. Mahmoud, Analysis of nanocontact problems of layered viscoelastic solids with surface energy effects under different loading patterns. Acta Mech. 227(2), (2016) 527–548.

[34] M. Bahrampour, S. J. Hamzeh and S. Shojaee, New Insight into Vis-coelastic Finite Element Modeling of Time-Dependent Material Creep Problems Using Spherical Hankel Element Framework. International Journal of Applied Mechanics, 1850085, 2018.

[35] J. N. Reddy, Theory and analysis of elastic plates and shells. 2nd ed. Boca Raton, FL: CRC Press; 2007.

[36] J. N. Reddy, Microstructure-dependent couple stress theories of functionally graded beams. J Mech Phys Solids, 59(2018) 2382–99.

[37] O. Martin, A modified variational iteration method for the analysis of viscoelastic beams. Applied Mathematical Modelling, 40(17–18), (2016) 7988–7995.

[38] H. H. Hilton, Viscoelastic Timoshenko beam theory, Mech of Time-Depend Mat, 13(1), (2009) 1–10.

[39] W. N. Findley and F. A. Davis, Creep and relaxation of nonlinear viscoelastic materials. Courier Corporation, 2013.

[40] V. P. Vallala, G. S. Payette and J. N. Reddy, A spectral/hp nonlinear finite element analysis of higher-order beam theory with viscoelasticity. Int J Appl Mech, 4(1), 1250010, (2012).

[41] S. Yi and H. H. Hilton, Dynamic finite element analysis of viscoelastic composite plates in the time domain. International Journal for Numerical Methods in Engineering, 37(23), 4081–4096, (1994).

[42] A. Capsoni, G. M. Viganò and K. Bani-Hani, On damping effects in Timoshenko beams. Int J mech sci, 73 (2013) 27–39.

[43] P. Kohnke, ANSYS Theory Reference, Release 5.5 ANSYS. Inc., Canonsburg, (1997).

[44] E. Oñate, Structural analysis with the finite element method. Linear statics: volume 2: beams, plates and shells. Springer Science & Business Media (2013).

[45] D. Froio and E. Rizzi, Analytical solution for the elastic bending of beams lying on a variable Winkler support. Acta Mechanica, 227(4), (2016) 1157–1179.

[46] A. S. Da Costa Azevêdo, A. C. A. Vasconcelos and S. dos Santos Hoefel, Dynamic analysis of elastically supported Timoshenko beam. Revista Interdisciplinar de Pesquisa em Engenharia-RIPE, 2(34), (2016) 71–85.

[47] A. A. Abdelrahman, M. A. Eltaher, A. M. Kabeel, A. M. Abdraboh and A. A. Hendy, Free and forced analysis of perforated beams. Steel and Composite Structures, 31(5), (2019), 489–502.

[48] G. A. Hassaan, Car Dynamics Using Quarter Model and Passive Suspension, Part IV: Destructive Miniature Humps (Bumps). Global Journal of Advanced Research, 2(2), (2015) 451–463.

Biographies

Alaa A. Abdelrahman Associate Professor in mechanical design and production Engineering Department at Zagazig University, Faculty of Engineering since July 2019. His research focuses on Computational Mechanics of Deformable solids and structures.

Ayman E. Nabawy received his M.Sc. degree in Mechanical Engineering at the Zagazig University in 2019. He is currently a PhD student at the department of Mechanical Design and Production Engineering, faculty of Engineering at Zagazig University. His research focuses on Modeling and analysis of the dynamic behavior of deformed solids and structures.

Soliman S. Alieldin, Professor in mechanical design and production Engineering Department at Zagazig University, Faculty of Engineering. His research focuses on Theoretical, Experimental, and Computational Mechanics of materials and structures.

Ayman M. M. Abdelhaleem, Associate Professor in mechanical design and production Engineering Department at Zagazig University, Faculty of Engineering. His research focuses on Analysis, Design, control of motor vehicle systems.

European Journal of Computational Mechanics, Vol. 28-6, 573–604.
doi: 10.13052/ejcm2642-2085.2862
© 2020 River Publishers

Abstract

1 Introduction

2 Damping Mechanisms

3 Mathematical Formulation of Viscoelastically Damped Suspension System

4 Numerical Evaluation of the Viscoelastic Constitutive Equations

5 Finite Element Model

6 Validation of the Developed Model

7 Computational Results

7.1 Effect of Vehicle Velocity on the Transient Response of the Lower Control Arm

7.2 Effect of the Viscoelastic Relaxation Time on the Transient Response of the Lower Control Arm

8 Conclusion

Appendix A
A Numerical Approximation for the Convolution Integral

References

Biographies