A Dynamic Response Analysis of Vehicle Suspension System

: Automotive engineering is a very important branch of design and development in mechanical engineering devoted to vehicle manufacturing. It demands solid knowledge focusing on the kinematics and dynamics of mechanical systems, either using structural analysis or modelling techniques. In this manuscript, a model based on the application of the fundamentals of the static and dynamic behaviour of the suspension system of a bus vehicle is developed. Therefore, a set of mathematical equations is expanded following the structural dynamic formulations to obtain an analytical solution on natural frequencies and the displacement under an external load. The problem is also solved using the Finite Element Method (FEM) to conﬁrm the analytical results. Equivalent material parameters and boundary conditions were used for the FEM simulation. This work aims to propose a methodology based on modal analysis with a simpliﬁed model to predict the dynamic response of an automotive structure.


Introduction
Automotive suspension systems are designed to control the passenger's comfort and safety when riding in a vehicle.The mechanical design of the suspension system regulates the vehicle's performance under various driving conditions [1].Therefore, analysing the suspension control systems is a key to improving the smoothness of a vehicle's ride without compromising safety.The balance between the comfort and handling stability of a vehicle under traffic conditions must be controlled [2,3].
Conventional suspensions are fabricated and assembled using traditional manufacturing methods.Thus, they encompass some limitations, including strict design constraints, high costs for parts, and modification [4].Currently, there are three types of automotive suspensions commonly used: passive, semi-active, and active [5].
A passive vehicle suspension system is reliable, simple, and inexpensive, being a typical mechanical structure incorporating linear springs and viscous dampers with constant stiffness and damping coefficient, respectively.In such systems, the damper and spring are fastened between the car body and wheel support components.A damper is engaged with hydraulic oil or compressed gas.There is a piston moving through a rod from the exterior [5].Semi-active suspension is similar to the passive design.It differs in the presence of variable damping coefficients but still has a fixed spring constant.Semi-active design contributes to a seamless change between a passive damper and one with a semi-active damping coefficient [5].Compared with the other suspension systems, an active apparatus comprises an actuator that can supply active force regulated by a controlling algorithm.It relies on the information gathered from installed vehicle sensors [5].
Nevertheless, most vehicles on the road do not possess an active suspension.As of now, active suspension is reserved for higher-end and performance vehicles.Though, given Appl.Sci.2023, 13, 2127 2 of 19 the direction of innovation and the addition of more features, active suspension will likely become more commonplace in the near future [6].
Despite the common use of a passive suspension system in recent vehicles, the engineering inherent in its design promotes the inclusion of an efficient stability control system.It will counteract the detrimental effects of sharp external forces that may come from evasive emergency driving manoeuvres [7].Many accidents occur due to an improper suspension design, mechanical malfunction, underinflated tyres, or worn/inadequate dampers [8,9].Therefore, enhancing the quality of transportation services has been recognised as one of the key factors toward sustainable mobility [10,11].

Related Works
The vehicle's behaviour is naturally conditioned by the effect of time-dependent (dynamic) forces on the vehicle body.Furthermore, the tyre/suspension system quality has a direct impact on the vehicle's on-road performance.Miloradović et al. [12] conducted a detailed study on the dynamic interaction between a passenger vehicle's steering and the physical connection of the suspension bodies.They performed experiments to acquire acceleration, deformation, and torque data at various driving speeds on different roads.
To characterise the problem variables and mathematical modelling, there exist a considerable number of mathematical methodologies to study the vehicle's dynamics.Talukdar et al. [13] proposed some ride models to evaluate the parameters affecting the vehicle's behaviour.They studied a half-model of a car, having extended the analysis from a design perspective, e.g., variation of suspension stiffness, damping coefficient, etc.In a particular context of vehicle design for a safe ride, the suspension performance has drawn the attention of mathematicians for applied industrial problems such as suspension optimization for non-linear stiffness in terms of stroke, and it offers a progressive strength to vehicle roll in a curve or a yaw motion [14].
Currently, computational modelling is widely used to solve problems involving fast or transient loading status [15].The high costs associated with developing new vehicle models have also made computer simulations of a vehicle's dynamics increasingly relevant [16].Several studies were performed to characterize the vehicle dynamics behaviour, identifying the suspension and tyre assembly as two of the key components that most influence the ride quality [17,18].Likewise, the advancement of technology in electronics, sensors, and actuators has enabled the application of control systems to enhance the vehicle's ride quality [19].
In the field of automotive engineering tools for design and structural assessment, the Finite Element Method (FEM) is considered a crucial technique for a reliable and accessible approach to studying the structural assessment of a vehicle's dynamics [20].Thus, FEM models contribute to the possibility of simulating the model in any case of static or dynamic loading conditions, thermal effects, fluid-structure interaction, and electromagnetic field effects [19].Its accuracy has been confirmed by reasonable agreement with corresponding experimental tests [21].Furthermore, the FEM was used to study the natural frequencies since it deals with an impressively large set of parameters [22,23].Abdullah et al. [23] sought to explore the level of accuracy of the simulated simplified body using various modelling methodologies.They further validated such models by comparing finite element modal characteristics with experimental modal properties.
A simplified model of a passive suspension system installed on a bus vehicle is proposed here.The vehicle body is assumed to be extremely rigid compared with the tyre/spring suspension group.Therefore, a rectangular, rigid plate supported by four corner springs is established.Then, a set of mathematical formulations is developed relying on the dynamics modal analysis theory.It must be mentioned that the bus body is included in the calculation because of its equivalent mass density.Since the suspension spring stiffness is much smaller than the one assigned to the structural components, only the elastic effects of the vehicle suspension will be taken into account in the calculation.Furthermore, this study intends to model the effect of the road irregularities or the driver's actions via time dependency through transient dynamic formulations integrated with the real-time direct integration analysis of Newmark's method.The normal displacement field of the studied model is established by a 3-DOF model consisting of a concentrated vertical translation of the mass centre superimposed on two rotations about the xand y-axes, respectively.
To confirm the analytical solutions, numerical modelling is performed through FEM software ABAQUS ® .As a result, natural frequencies are computed through a frequency step scheme.Then, the FEM model is solved by a steady-state dynamic, direct step in which two concentrated forces are applied to the front side of the model, which simulates the response of the bump impacts as theoretically considered.

Mathematical Development
Based on the modal analysis theory, the stiffness matrix, together with the displacement and force vectors, is defined to provide the natural frequency solution obtained from the eigenvalues of the suggested geometry.Figure 1 shows a general view of the studied model, in which L b and b represent the bus length and width, respectively.The wheelbase length is identified as L and the track width is defined as b w .Each wheel hub is supported by a spring set [24].Table 1 reports the main characteristics of the studied model.
Furthermore, this study intends to model the effect of the road irregularities or the driver's actions via time dependency through transient dynamic formulations integrated with the real-time direct integration analysis of Newmark's method.The normal displacement field of the studied model is established by a 3-DOF model consisting of a concentrated vertical translation of the mass centre superimposed on two rotations about the x-and yaxes, respectively.
To confirm the analytical solutions, numerical modelling is performed through FEM software ABAQUS ® .As a result, natural frequencies are computed through a frequency step scheme.Then, the FEM model is solved by a steady-state dynamic, direct step in which two concentrated forces are applied to the front side of the model, which simulates the response of the bump impacts as theoretically considered.

Mathematical Development
Based on the modal analysis theory, the stiffness matrix, together with the displacement and force vectors, is defined to provide the natural frequency solution obtained from the eigenvalues of the suggested geometry.Figure 1 shows a general view of the studied model, in which  and  represent the bus length and width, respectively.The wheelbase length is identified as  and the track width is defined as  .Each wheel hub is supported by a spring set [24].Table 1 reports the main characteristics of the studied model.

Matrix Formulation of the Model
This study focuses on the vehicle's dynamic interaction with the road when the vehicle is moving.Therefore, different states, including driving manoeuvres such as acceleration, braking, and curving, are examined according to the suspension flexibility and damping characteristics.The proposed model is defined as a four-node macro-element with a rectangular shape, dimensioned as wheelbase length  and bus width , respectively, in the  − and  − coordinate axes, as depicted in Figure 2.
This rigid plate (the simplified geometry of the bus) materializes the structure of the bus platform.It is assumed to be very rigid compared with the stiffness of the suspension

Matrix Formulation of the Model
This study focuses on the vehicle's dynamic interaction with the road when the vehicle is moving.Therefore, different states, including driving manoeuvres such as acceleration, braking, and curving, are examined according to the suspension flexibility and damping characteristics.The proposed model is defined as a four-node macro-element with a rectangular shape, dimensioned as wheelbase length L and bus width b, respectively, in the x− and y− coordinate axes, as depicted in Figure 2.
This rigid plate (the simplified geometry of the bus) materializes the structure of the bus platform.It is assumed to be very rigid compared with the stiffness of the suspension spring/damper assembly at each wheel.With this characteristic, the displacement field bears on a simple algebraic equation defining a plane.
spring/damper assembly at each wheel.With this characteristic, the displacement field bears on a simple algebraic equation defining a plane.To define the spatial position of the plane points by its transverse displacement  (, ) ( = 1, ..4) at each corner node, some parameters must be assigned in the system of equations as: In which: •  is the transverse displacement related to the plate mass centre, ; •  and  are the plate rotations about  − and  − coordinate axes.
If the plate element is supported by 4 springs at corner nodes with stiffness  ,  ,  , and  , then, for a generic position of the rectangular plate defined by the nodal displacements, referring to Equation (1), the elastic deformation energy stored in the plate is obtained by the total contribution of the deformed springs at each node.Indeed, the total amount of the internal energy (elastic energy) must balance the external work performed by the external forces (or moments) applied to the plate as: where the nodal displacement  is a function of parameters presented in Equation (1).
In relation to the elastic energy  of the deformed plate, the nodal equivalent external force system must be obtained to equilibrate internal forces as follows: Taking the derivatives of these equations, the chain rule for the independent variables  ,  , and  would be acquired, therefore, the linear system of 3 equations can take the following form: To define the spatial position of the plane points by its transverse displacement u i (x, y) (i = 1, . . .4) at each corner node, some parameters must be assigned in the system of equations as: In which: • u G is the transverse displacement related to the plate mass centre, G; • θ x and θ y are the plate rotations about x− and y− coordinate axes.If the plate element is supported by 4 springs at corner nodes with stiffness k 1 , k 2 , k 3 , and k 4 , then, for a generic position of the rectangular plate defined by the nodal displacements, referring to Equation (1), the elastic deformation energy stored in the plate is obtained by the total contribution of the deformed springs at each node.Indeed, the total amount of the internal energy (elastic energy) must balance the external work performed by the external forces (or moments) applied to the plate as: where the nodal displacement u i is a function of parameters presented in Equation (1).In relation to the elastic energy U elastic of the deformed plate, the nodal equivalent external force system must be obtained to equilibrate internal forces as follows: Taking the derivatives of these equations, the chain rule for the independent variables u G , θ x , and θ y would be acquired, therefore, the linear system of 3 equations can take the following form: It must be mentioned that using springs with stiffness, k, at the four corner nodes of the plate contributes to rendering the previous stiffness matrix diagonal, which led to significantly simplifying the problem complexity once a set of uncoupled systems of equations is attained: As an illustration, if a static concentrated transverse force F 1 at node 1 is considered prior to evaluating the RHS (Right Hand Side) of equation system (6), it is important to recover the transverse displacement at node 1 presented in Equation ( 1).Hence, the work W released by a force term of F 1 is: Applying Castigliano's theorem [25], the external force system, performing the work equivalent to the one due to force F 1 , can be formulated as: Therefore, the solution by uncoupled equations is computed as: This is the RHS due to the contribution of force at node 1, which changes the spatial configuration of the chassis plate model.However, in terms of nodal displacement, substituting this vector in Equation ( 1) results in: Figure 3 depicts the aspect of the deformed shape regarding the plane.
If the reactions of each corner due to the mounted springs are taken into account, then the forces are: A dynamic analysis of structures requires not only the definition of their internal reactions due to external force systems but also the computation of inertial forces as a consequence of the acceleration field.To achieve this goal, it is necessary to evaluate the structure mass matrix, a task that can be solved by variational techniques on evaluating the work performed by inertial forces (involving accelerations) on the respective displacement field.If the reactions of each corner due to the mounted springs are taken into account, then the forces are: A dynamic analysis of structures requires not only the definition of their internal reactions due to external force systems but also the computation of inertial forces as a consequence of the acceleration field.To achieve this goal, it is necessary to evaluate the structure mass matrix, a task that can be solved by variational techniques on evaluating the work performed by inertial forces (involving accelerations) on the respective displacement field.
Considering that the bus platform behaves as a plane, any displacement at a point can be interpolated from the 4 nodal displacements using shape or interpolation polynomial functions.Then, in a matrix, the displacement vector of any internal point can take the following form: where  …  and  …  are accordingly the nodal shape functions and the displacement.Indeed, a similar interpolation applies to the acceleration field, which is calculated in a similar mode to the displacements: Table 2 reports the corresponding interpolation functions.Considering that the bus platform behaves as a plane, any displacement at a point can be interpolated from the 4 nodal displacements using shape or interpolation polynomial functions.Then, in a matrix, the displacement vector of any internal point can take the following form:

Shape 𝑵 𝒊
where N 1 . . .N 4 and u 1 . . .u 4 are accordingly the nodal shape functions and the displacement.Indeed, a similar interpolation applies to the acceleration field, which is calculated in a similar mode to the displacements: ..
Table 2 reports the corresponding interpolation functions.
Table 2. Integral form of shape functions products.

Shape N i
Product Shape Figure 4 represents the finite element with four nodes.Assuming that each variable is attributed to the relation of x = L 2 (1 + ξ); y = b 2 (1 + η), the corresponding differential relation is described as dx = L 2 dξ; dy = b 2 dη.Hence, the work performed by inertial forces, assuming a consistent mass distribution, takes the following form: relation is described as  = ;  = .Hence, the work performed by inertial forces, assuming a consistent mass distribution, takes the following form: The consistent 4 × 4 DOF's mass matrix is obtained through solving the integrals with the shape function combination for all four nodes of the plate element, as demonstrated in Equation ( 15) and Equation ( 16).
The recently presented mass matrix is associated with a 4-node element defined by 4 nodal displacements.However, such a displacement field, as shown in Equation ( 1), can be condensed into a 3-DOF as: The scenario is the same regarding the acceleration vector.It is possible to reformulate Equation (8) for previous transformations to achieve a 3 × 3 mass matrix.Previous algebraic operations lead to a useful 3 × 3 mass matrix that is diagonal as: Consequently, the stiffness matrix was formulated as a 3×3 matrix, which contributes to providing a simple and practical solution to the problem of natural frequencies.It is remarkable to note that, in this orthogonalized mass matrix, a mass is associated with a displacement DOF, while moments of inertia are associated with rotation DOFs,  , and  .The consistent 4 × 4 DOF's mass matrix is obtained through solving the integrals with the shape function combination for all four nodes of the plate element, as demonstrated in Equations ( 15) and ( 16).

Modal Analysis
The recently presented mass matrix is associated with a 4-node element defined by 4 nodal displacements.However, such a displacement field, as shown in Equation ( 1), can be condensed into a 3-DOF as: The scenario is the same regarding the acceleration vector.It is possible to reformulate Equation (8) for previous transformations to achieve a 3 × 3 mass matrix.Previous algebraic operations lead to a useful 3 × 3 mass matrix that is diagonal as: Consequently, the stiffness matrix was formulated as a 3 × 3 matrix, which contributes to providing a simple and practical solution to the problem of natural frequencies.It is remarkable to note that, in this orthogonalized mass matrix, a mass is associated with a displacement DOF, while moments of inertia are associated with rotation DOFs, θ x , and θ y .

Modal Analysis
A modal analysis of a rigid, rectangular, solid plate with nodal-assigned elastic elements can be presented as a 4-node macro element with specific spring stiffness.This example is not representative of a design solution for a bus chassis (there are no realistic solution options in adopting a rigid plate as a passenger vehicle chassis base); however, this can work as a dynamically equivalent model regarding the vehicle mass to the centre of mass.In addition, inertial rotation elements can be considered for any of the coordinate x-y axes passing the centre of mass, as shown in Figure 1.Polar moments of inertia, which are the sum of two axial moments over y− and z− axes can take the following form:
Therefore, the eigenvalue equation can be computed as: The first three natural frequencies of the studied model can be derived from Equation (19).Hence, Table 3 reports the obtained natural frequencies from the developed analytical methodology identified as f ANL i .In addition, the deformed shape of the studied model is shown in Figure 5.
this can work as a dynamically equivalent model regarding the vehicle mass to the centre of mass.In addition, inertial rotation elements can be considered for any of the coordinate x-y axes passing the centre of mass, as shown in Figure 1.Polar moments of inertia, which are the sum of two axial moments over  − and  − axes can take the following form: The first three natural frequencies of the studied model can be derived from Equation (19).Hence, Table 3 reports the obtained natural frequencies from the developed analytical methodology identified as  .In addition, the deformed shape of the studied model is shown in Figure 5.
Rotation over X G axis;

Transient-Dynamic Formulation
As mentioned, the vehicle mass and inertial properties are approached by a simple rectangular, solid plate element with spring suspension at each corner.Once it was possible to obtain the solution of the modal analysis of the problem from a set of uncoupled equations, the transient analysis would be obtained in a similar practical mode.The dynamic equilibrium equation for the oscillating plate can be re-written as follows: As previously mentioned, the RHS of this dynamic equilibrium must be transformed to equivalent moments while dealing with nodal forces.As an illustration, having four independent nodal forces acting transversely on the plate plane, the inverse of Equation ( 1) produces the expression of the condensed displacement vector (composed of the mass centre displacements and rotations about coordinate axes), therefore: Then, the corresponding force vector, if composed of four time-dependent components, will be: Considering that there exists a force F u 1 acting at node 1, the equivalent (condensed) force vector is represented in Equation (23).The complete Equation ( 21) can be extended to dynamic situations of transient forces applied at nodes as: A practical and efficient mode to solve this equation for generalized time-dependent loads relies on direct time integration methods.It could be developed from the mode superposition technique, another way to setup the solution for this dynamic problem.However, in the case of complex loads, where their decomposition into elementary modes is necessary, direct time integration is a more practical tool, limited only by the computation time in the case of very large DOF structure modelling.We adopted the Newmark's constant acceleration method [26] for its accuracy, stability, and ease to program for computational purposes.An additional reason for this option is the fact that the solution is structured as an uncoupled system of equations.
The dynamic equilibrium of a structure subjected to external time-dependent loads can be analysed by a direct time integration method such as Newmark's algorithm [26].This algorithm calculates the shape of a structure based on the concept of Constant Acceleration Method, where at each time step the structure's acceleration vector is supposed to remain constant.To develop the algorithm, the structure equilibrium at an instant t + ∆t can be formulated as follows: .
U t+∆t , and U t+∆t are accordingly the acceleration, velocity, and displacement vectors of the structure, while F t+∆t is the external force vector, all at time step t + ∆t.If velocity and displacement vectors are considered from the previous instant value, t, it is feasible to obtain the structural evolution in the next time step under dynamic loads.Thus, the constant acceleration concept and the trapezoidal rule in the integration of a step-function can provide the following equations: .
These equations are indeed applied to the plate rotations instead of translation.The evaluation of .. U t+∆t is the key solution of the algorithm; for that, previous expressions are substituted in the structure equilibrium equation at the time: t + ∆t: On separating terms referred to time t + ∆t and t, respectively, at the left and righthand sides of the equation, it is possible to solve it for the updated structure acceleration vector as: Once solved for ..
U t+∆t , knowing the previous dynamic status at instant t (that is the RHS of the equation), an updated vector for acceleration is attained.Thus, it is leading to new velocity and displacement vectors, as in Equation (25).The equations for each DOF are uncoupled.The iterative incremental solutions can take the following form: .
As an illustration, the analysed plate provided the largest natural frequency f 3 = 2.0 (Hz), where t 3 = 0.5 (s) and thus, ∆t ≤ 0.5 (s).With this value, the solution is warranted for its stability.Hence, the smaller the time step, the better the accuracy, although a value much smaller than the recommended one would be useless since the computational time would increase largely.For instance, if a vehicle accelerates during t = 0.5 (s) with no damping, an impulsive force is due to the load on the front axis during F = 2 × 2.5 (kN) = 5.0 (kN) (hubs 1 and 4).Basic equation for centre of mass vibration with no damping is: Regarding the natural frequency detected by Newmark's method, one cycle = 1.02 (s); therefore, an analytical evaluation can lead to: Referring to Figure 6, the maximum displacement of the mass centre is u Max Regarding the natural frequency detected by Newmark's method, one cycle = 1.02 (s); therefore, an analytical evaluation can lead to: Referring to Figure 6, the maximum displacement of the mass centre is  = 1.23 × 10 () at the transient step.Concerning the rotation about  −  (forces are supposed to act at nodes 1 and 4), 2nd uncoupled equation: dynamic inertial moment over transverse  − .The dimensions of the bus body chassis are listed below.• Equivalent forces: Rotation  after applying  = 2 × 2.5 (kN) suddenly (and sustainably) at nodes 1 and 4. Thus, the equivalent moment is  = −30.0(kN.m).
From the graph plotted in Figure 7, the maximum rotation is obtained as θ y = −0.00929(rad).Hence, analytically, true displacement at front of vehicle model (nodes 1 and 4):

Finite Element Analysis
Computational simulations have been performed to confirm the analytical hypothesis presented in Section 3. In this regard, the FEM stands as an accurate methodology to highlight the most critical aspects where the analytical approach faced difficulties.Regarding the numerical model, the vehicle chassis is mainly manufactured from steel, as its geometry shows in Figure 8a.
Numerical modelling has been performed based on the FEM formulations implemented in ABAQUS© [27].The problem was modelled according to the real conditions, c.f. Figure 8.The vehicle chassis is modelled as a shell supported by four clamped springs simulating the vertical tyre/suspension stiffness as reported in [24,28,29].To define the boundary conditions, eight reference points ( ,  = 1, … ,8) were assigned to the model.Therefore, the springs (stiffness of  = 185 (kN m ⁄ )) were defined through pairing the defined points.Consequently,  , …  are located on the wheel hubs, representing the spring interaction, while the other points signify the restricted displacement at the spring's end, defined by  =  =  = 0.
The non-linear geometrical effects remain deactivated, and the numerical simulation has been carried out through a static analysis.At the rear (nodes 2 and 3), there will be a lowering of the pavement level, thus:

Finite Element Analysis
Computational simulations have been performed to confirm the analytical hypothesis presented in Section 3. In this regard, the FEM stands as an accurate methodology to highlight the most critical aspects where the analytical approach faced difficulties.Regarding the numerical model, the vehicle chassis is mainly manufactured from steel, as its geometry shows in Figure 8a.

Finite Element Analysis
Computational simulations have been performed to confirm the analytical hypothesis presented in Section 3. In this regard, the FEM stands as an accurate methodology to highlight the most critical aspects where the analytical approach faced difficulties.Regarding the numerical model, the vehicle chassis is mainly manufactured from steel, as its geometry shows in Figure 8a.
Numerical modelling has been performed based on the FEM formulations implemented in ABAQUS© [27].The problem was modelled according to the real conditions, c.f. Figure 8.The vehicle chassis is modelled as a shell supported by four clamped springs simulating the vertical tyre/suspension stiffness as reported in [24,28,29].To define the boundary conditions, eight reference points ( ,  = 1, … ,8) were assigned to the model.Therefore, the springs (stiffness of  = 185 (kN m ⁄ )) were defined through pairing the defined points.Consequently,  , …  are located on the wheel hubs, representing the spring interaction, while the other points signify the restricted displacement at the spring's end, defined by  =  =  = 0.
The non-linear geometrical effects remain deactivated, and the numerical simulation has been carried out through a static analysis.Numerical modelling has been performed based on the FEM formulations implemented in ABAQUS© [27].The problem was modelled according to the real conditions, Figure 8.The vehicle chassis is modelled as a shell supported by four clamped springs simulating the vertical tyre/suspension stiffness as reported in [24,28,29].To define the boundary conditions, eight reference points (P i , i = 1, . . ., 8) were assigned to the model.Therefore, the springs (stiffness of k = 185 (kN/m)) were defined through pairing the defined points.Consequently, P 1 , . . .P 4 are located on the wheel hubs, representing the spring interaction, while the other points signify the restricted displacement at the spring's end, defined by u x = u y = u z = 0.

FE Mesh Convergence Study
The non-linear geometrical effects remain deactivated, and the numerical simulation has been carried out through a static analysis.

FE Mesh Convergence Study
This section is devoted to a convergence study carried out on the FE mesh to assess how the mesh would affect the results and particularly to determine the optimal mesh density.The chassis is considered a plate on which shell elements can be employed to analyse the problem.Different shell element types were studied: quadratic elements with/without reduced integration for four and eight nodes identified as S4R, S4, and S8R.
To achieve the optimal mesh, the model presented in Figure 8 was analysed.Regarding the boundary condition, a vertical displacement as V(x, y) = 10 (mm) was imposed at Node 1.The numerical simulations have been performed on different FE mesh types with an element size of 200 (mm), which is reported in Table 4.However, Table 5 presents a summary of the obtained displacement along the z-direction.As a conclusion, different element types do not have a significant impact on the results.Based on the findings of the FE convergence investigation, it can be concluded that the S4R element type would be the best option for ensuring solution correctness and minimizing computing costs.To study the FE mesh size, various element sizes were considered, as listed in Table 6.Maintaining the same boundary conditions mentioned before and considering the different FE meshes reported in Table 6, the problem was numerically solved, and the displacement variation over z-direction was obtained for all nodes as reported in Table 7.
Table 7. Displacement along z-direction for different mesh densities of S4R.Owing to the results obtained from various FE densities, it can be inferred that medium mesh density produced accurate results while minimizing computational costs and increasing convergence.

Static Analysis
An external load impacting the connection ground-vehicle is generated by modelling a ride up of 10 (mm) or more sharp road imperfections, for instance 100 (mm).This vertical displacement is imposed on the node 1 as demonstrated in Figure 8. Figure 9 depicts the displacement field in the direction of the displacement imposition.Referring to Figure 9a, it is notable that the maximum displacement was reached where the displacement was imposed, while the other two neighbours (Node 2 and Node 4) tend to move in the same direction as Node 1. Regarding Node 3, the displacement is moving in the opposite direction of the other nodes.In Figure 9b, it can be inferred that there is a tendency for rotation along a parallel diagonal connecting Node 2 and Node 4.

Static Analysis
An external load impacting the connection ground-vehicle is generated by modelling a ride up of 10 (mm) or more sharp road imperfections, for instance 100 (mm).This vertical displacement is imposed on the node 1 as demonstrated in Figure 8. Figure 9 depicts the displacement field in the direction of the displacement imposition.Referring to Figure 9a, it is notable that the maximum displacement was reached where the displacement was imposed, while the other two neighbours (Node 2 and Node 4) tend to move in the same direction as Node 1. Regarding Node 3, the displacement is moving in the opposite direction of the other nodes.In Figure 9b, it can be inferred that there is a tendency for rotation along a parallel diagonal connecting Node 2 and Node 4. Table 8 reports the obtained displacement on different corners.An external load that acts unexpectedly on the vehicle can cause the excitation of the body to step into resonance.Hence, it is rational to study the natural excitation frequencies of the model, which are discussed in the next section.

Frequency Analysis
This section is devoted to the numerical analysis of the example that was analytically solved in Section 3.2.The problem was resolved using a Frequency step to compute the natural frequencies following Lanczos eigensolver [30].The FEM analysis was carried out excluding the acoustic structural coupling due to the problem conditions.The eigenvectors were normalized by the displacement.The displacement normalization is a technique to represent mode shapes in a modal analysis where the peak amplitude is normalized to a value of 1.This is a common method to signify mode shapes not only using commercial tools but also for general purposes while representing an analytical solution.The eigenvector normalization allows scaling up or down the mode shape amplitude and hence does not affect frequencies.An external load that acts unexpectedly on the vehicle can cause the excitation of the body to step into resonance.Hence, it is rational to study the natural excitation frequencies of the model, which are discussed in the next section.

Frequency Analysis
This section is devoted to the numerical analysis of the example that was analytically solved in Section 3.2.The problem was resolved using a Frequency step to compute the natural frequencies following Lanczos eigensolver [30].The FEM analysis was carried out excluding the acoustic structural coupling due to the problem conditions.The eigenvectors were normalized by the displacement.The displacement normalization is a technique to represent mode shapes in a modal analysis where the peak amplitude is normalized to a value of 1.This is a common method to signify mode shapes not only using commercial tools but also for general purposes while representing an analytical solution.The eigenvector normalization allows scaling up or down the mode shape amplitude and hence does not affect frequencies.
Conforming to the model shown in Figure 8b, the frequency analysis has been performed, requesting the first six eigenvalues.The first three modes are zero due to the representation of rigid body movement.Table 9 summarizes the obtained results on numerical natural frequencies identified as f FEM i .Regarding the second harmonic mode, it reflects a rotation about the x-axis according to the coordinate system presented in Figure 8.The suspension system is established by four springs, which vertically vibrate with a 180-degree phase shift.It must be mentioned that the third harmonic mode accounted for a higher frequency due to a rotation that occurred along the diagonal of the corners.As a conclusion, it is possible to state that the natural frequencies are independent of the displacement imposed.
Figure 10 shows the magnitude displacement profile obtained on the second and third natural frequencies if the imposed displacement is 10 (mm) on Node 1.As shown in Figure 10a, the suspension compresses and extends vertically, while the third frequency experienced a rotation along the diagonal direction as depicted in Figure 10b.This behaviour should be similar for different displacement conditions.Therefore, it is feasible to predict the shape of the response of the model if subjected to such displacement imposition.
Conforming to the model shown in Figure 8b, the frequency analysis has been performed, requesting the first six eigenvalues.The first three modes are zero due to the representation of rigid body movement.Table 9 summarizes the obtained results on numerical natural frequencies identified as  .Regarding the second harmonic mode, it reflects a rotation about the x-axis according to the coordinate system presented in Figure 8.The suspension system is established by four springs, which vertically vibrate with a 180-degree phase shift.It must be mentioned that the third harmonic mode accounted for a higher frequency due to a rotation that occurred along the diagonal of the corners.As a conclusion, it is possible to state that the natural frequencies are independent of the displacement imposed.Figure 10 shows the magnitude displacement profile obtained on the second and third natural frequencies if the imposed displacement is 10 (mm) on Node 1.As shown in Figure 10a, the suspension compresses and extends vertically, while the third frequency experienced a rotation along the diagonal direction as depicted in Figure 10b.This behaviour should be similar for different displacement conditions.Therefore, it is feasible to predict the shape of the structural response of the model if subjected to such displacement imposition.

Steady-State Dynamic Analysis
This section deals with the numerical analysis of the problem that was analytically solved in Section 3.3.ABAQUS/Standard provides a "direct" steady-state linear dynamic analysis procedure in which the equations of steady harmonic motion of the system are solved directly without using the eigenmodes.Thus, the model was simulated following a steady-state dynamic, direct step in which two concentrated forces were applied to the front side of the model (Node 1 and Node 2) as shown in Figure 11a.
The direct solution scheme solves the complete set of nodal DOFs at each excitation frequency, and the frequency excitation is not required since the eigenmodes can be used to identify the excitation frequency.It is rather expensive since a complex solution must be obtained for the full system of equations at each excitation (driving) frequency.Applied forces can be defined as real (in-phase) and imaginary (out-of-phase) loads, and non-zero boundary conditions may be applied to any of the active DOF.
However, due to the nature of the problem, the system is considered an in-phase one.Therefore, the applied force has been defined as a real component with a magnitude of  = 2.5 (kN) on each node, as it was prescribed in the analytical model.The analysis has

Steady-State Dynamic Analysis
This section deals with the numerical analysis of the problem that was analytically solved in Section 3.3.ABAQUS/Standard provides a "direct" steady-state linear dynamic analysis procedure in which the equations of steady harmonic motion of the system are solved directly without using the eigenmodes.Thus, the model was simulated following a steady-state dynamic, direct step in which two concentrated forces were applied to the front side of the model (Node 1 and Node 2) as shown in Figure 11a.
The direct solution scheme solves the complete set of nodal DOFs at each excitation frequency, and the frequency excitation is not required since the eigenmodes can be used to identify the excitation frequency.It is rather expensive since a complex solution must be obtained for the full system of equations at each excitation (driving) frequency.Applied forces can be defined as real (in-phase) and imaginary (out-of-phase) loads, and non-zero boundary conditions may be applied to any of the active DOF.
However, due to the nature of the problem, the system is considered an in-phase one.Therefore, the applied force has been defined as a real component with a magnitude of F z = 2.5 (kN) on each node, as it was prescribed in the analytical model.The analysis has been conducted in the requested frequency range of f 1 = 1.07 to f 3 = 2.17 (Hz).Figure 11b demonstrates the displacement variation in the z-direction at a frequency of 1.07 (Hz).
been conducted in the requested frequency range of  = 1.07 to  = 2.17 (Hz).Figure 11b demonstrates the displacement variation in the z-direction at a frequency of 1.07 (Hz).

Results Discussion
Table 10 reports a comparison between the results acquired from the developed analytical solution and the numerical FEM.A reasonable agreement has been found, allowing the FEM model to confirm the proposed analytical methodology.From the acquired frequency data, it can be stated that the frequency increased for higher modes in the analytical solution since a higher vibrational deformation energy is associated with increasing natural frequencies, as confirmed by the FEM results.Consequently, it can be inferred that a low deviation for  and  has been obtained, and while it is a bit higher for  , it remains below 8%.
Regarding the steady-state dynamic analysis, the displacements are associated with the intensity of the natural vibrational mode.The external force system is superimposed with the internal inertial forces, leading to displacement as a consequence of the dynamic equilibrium regime of the structure.Thus, the maximum and minimum displacement are not equal in intensity and sign due to the superposition of a uniform transverse displacement and a rotation about the x-axis.

Results Discussion
Table 10 reports a comparison between the results acquired from the developed analytical solution and the numerical FEM.A reasonable agreement has been found, allowing the FEM model to confirm the proposed analytical methodology.From the acquired frequency data, it can be stated that the frequency increased for higher modes in the analytical solution since a higher vibrational deformation energy is associated with increasing natural frequencies, as confirmed by the FEM results.Consequently, it can be inferred that a low deviation for f 1 and f 2 has been obtained, and while it is a bit higher for f 3 , it remains below 8%.Regarding the steady-state dynamic analysis, the displacements are associated with the intensity of the natural vibrational mode.The external force system is superimposed with the internal inertial forces, leading to displacement as a consequence of the dynamic equilibrium regime of the structure.Thus, the maximum and minimum displacement are not equal in intensity and sign due to the superposition of a uniform transverse displacement and a rotation about the x-axis.

Conclusions
In this work, a single modular rectangular element was developed to accurately characterise the dynamic behaviour of a bus vehicle model with a passive suspension subjected to generalised transient forces.A mathematical set of equations was developed based on the modal analysis hypothesis to offer an analytical solution for natural frequency calculation.The obtained results infer that the performance of the proposed methodology is accurate in static and steady-state dynamics applications.However, the main conclusions of the present study can be drawn as follows: • The most valuable practical implication is a comprehensive vehicle suspension diagnostic system concept, from which a simplified model helps to produce outputs with sufficient accuracy in which an approximate structural behaviour of the physical model can be predicted; • Only the chassis-induced transverse displacement vibrational modes allow the direct comparison in which the basic model would be cross-checked;

•
In this paper, by means of simulation, the analytical control algorithm is used to evaluate the dynamic equilibrium of the structure, subjected to time-dependent domain loads, at any instant.Thus, the obtained results revealed a deviation of less than 8% about the displacement, allowing the model confirmation.This may provide a useful reference for further experimental investigations.

•
The analytical solution relies on a single finite element demonstrating a coherent dynamic response under transverse forces.It can practically model a less smooth track in addition to the presence of in-plane transient forces, such as as those due to sharp manoeuvres or the effect of side gusts; • A quite low number of DOFs and its accuracy make sound arguments for the applicability of this single modular element to automotive engineering in the fields of driving practice and suspension turning for dynamic safety.
In the future development of this single chassis block-like modulus, the incorporation of rotating inertial elements in the mass matrix and the steady state dynamic analysis integrating damping will result in an improved modal analysis and dynamic vehicle behaviour when subjected to external loads.The presentation of this structural dynamic modelling could be applied to the modal analysis of more elaborate construction vehicles in the future.

Data Availability Statement:
The data presented in this study are available on request from the corresponding authors.The data are not publicly available due to their containing information that could compromise the privacy of research participants.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 1 .
Figure 1.A general perspective of the studied model.

Figure 1 .
Figure 1.A general perspective of the studied model.

Figure 2 .
Figure 2. Platform model as a plane resting on 4 springs at corner nodes.

Figure 2 .
Figure 2. Platform model as a plane resting on 4 springs at corner nodes.

Figure 3 .
Figure 3. Deformed shape of the rigid plate under prescribed force at node 1.

Figure 3 .
Figure 3. Deformed shape of the rigid plate under prescribed force at node 1.

Figure 4 .
Figure 4. Notation and non-dimensional coordinates of rectangular finite element.

Figure 4 .
Figure 4. Notation and non-dimensional coordinates of rectangular finite element.

G = 1 . 23 ×
10 −2 (m) at the transient step.Concerning the rotation about y − axis (forces are supposed to act at nodes 1 and 4), 2nd uncoupled equation: dynamic inertial moment over transverse y − axis.The dimensions of the bus body chassis are listed below.Appl.Sci.2023, 13, x FOR PEER REVIEW 11 of 19

Figure 6 .
Figure 6.Displacement of mass centre, u G .

Figure 8 .
Figure 8. Numerical model; (a) geometrical perspective and (b) FE mesh, interactions and essential boundary condition presentation.An equivalent density ( = 3.7 × 10 tonne mm ⁄ ) was considered to obtain the real mass of the vehicle ( = 12 tonne).

Figure 8 .
Figure 8. Numerical model; (a) geometrical perspective and (b) FE mesh, interactions and essential boundary condition presentation.An equivalent density ( = 3.7 × 10 tonne mm ⁄ ) was considered to obtain the real mass of the vehicle ( = 12 tonne).

Figure 8 .
Figure 8. Numerical model; (a) geometrical perspective and (b) FE mesh, interactions and essential boundary condition presentation.

Figure 11 .
Figure 11.Steady-state dynamic analysis by FEM; (a) natural boundary condition and (b) displacement in z-direction at a frequency of 1.07 (Hz).

Figure 11 .
Figure 11.Steady-state dynamic analysis by FEM; (a) natural boundary condition and (b) displacement in z-direction at a frequency of 1.07 (Hz).

Table 1 .
Characteristics of the studied model.

Table 1 .
Characteristics of the studied model.

Table 2 .
Integral form of shape functions products.

Table 3 .
Analytical natural frequency solutions.

Table 4 .
FE mesh properties for different element types.

Table 5 .
Displacement along z-direction for different FE types.

Table 6 .
FE mesh properties for the element type of S4R.
Table 8 reports the obtained displacement on different corners.

Table 10 .
Comparison amongst the results obtained from FEM and developed analytical approaches.

Table 10 .
Comparison amongst the results obtained from FEM and developed analytical approaches.Result Analytical − Result FEM /Result FEM .