Mathematical Modeling and Parameters Estimation of Car Crash Using Eigensystem Realization Algorithm and Curve-Fitting Approaches
An eigensystem realization algorithm (ERA) approach for estimating the structural system matrices is proposed in this paper using the measurements of acceleration data available from the real crash test. A mathematical model that represents the real vehicle frontal crash scenario is presented. The model’s structure is a double-spring-mass-damper system, whereby the front mass represents the vehicle-chassis and the rear mass represents the passenger compartment. The physical parameters of the model are estimated using curve-fitting approach, and the estimated state system matrices are estimated by using the ERA approach. The model is validated by comparing the results from the model with those from the real crash test.
Car crash test is usually performed in order to ensure safe design standards in crashworthiness (the ability of a vehicle to be plastically deformed and yet maintain a sufficient survival space for its occupants during crash scenario). Nowadays, due to advanced research in computer simulation software, simulated crash tests can be performed beforehand the full-scale crash test. Therefore, cost associated with real crash test can be reduced. Vehicle crashworthiness can be evaluated in four distinct modes: frontal, side, rear, and rollover crashes. Several researches have been carried out in this field, which resulted in several novel computational models of vehicle collisions in literature. In , a mathematical model is proposed to estimate the maximum occupant deceleration, which is one of the main tasks in the area of crashworthiness study by a Kelvin model which contains a mass together with spring and damper connected in parallel. An application of physical models composed of springs, dampers, and masses joined together in various arrangements for simulating a real car collision with a rigid pole was presented in .
In , the authors presented an overview of the kinematic and dynamic relationships of a vehicle in a collision, whereby the work was to identify the parameters of the vehicle crash model using experimental dataset. In [4, 5], a lumped parameter modeling in frontal crash was investigated and analyzed in five degrees of freedom and has been used to analyze the response of occupant during the impact. In [6, 7], an optimization procedure to assist multibody vehicle model development and validation was proposed. The authors first devised the topological structure of the multibody system representing the structural vehicle components and described the most relevant mechanisms of deformation. In the work of , the authors proposed an approach to control the seat-belt restraint system force during a frontal crash to reduce thoracic injury.
The main challenge in accident reconstruction is the system identification described as the process of constructing mathematical models of dynamical systems using measured input-output data. In case of vehicle crash, system identification algorithm consists in retrieving the unknown parameters such as the spring stiffness and damping coefficient. A possible approach is to identify these parameters directly from experimental dynamical data.
From literature, system identification algorithms (SIA) have been developed for different applications. Among others, we can state the following: subspace identification, genetic algorithm, eigensystem realization algorithm, and data-based regressive model approaches. Typical examples where these SIA have been used can be found in [9, 10], [4, 11], and [12, 13], respectively.
System identification using ERA has so far received considerable attention, as evidenced by the works of Juang , Chiang and Lin , Ko and Hung , Juang and Pappa , and Yang and Yeh . For instance, in , the author developed the ERA to estimate the natural frequencies and damping ratios of a dynamical system from known Markov parameters, and in , the authors used the ERA to identify the system matrices of a vibrating structure from the displacement-based Markov parameters, which were estimated from measured displacement responses together with the excitation forces.
The main contribution of this paper is the development of a mathematical model for a double-spring-mass-damper system which reconstructs a vehicle frontal crash scenario and estimates structural parameters such as natural frequencies, spring stiffness, and damping coefficients of the system. In this paper, an ERA is used to identify the system matrices of a vehicle impacting a rigid barrier, modeled by a double-spring-mass-damper system. The model represents the inertia of the vehicle chassis and the passenger compartment. The state-space representation of the model is estimated from the acceleration-based Markov parameters which are extracted from the measured acceleration response. To estimate the physical parameters (stiffness and damping coefficient) of the model, a curve-fitting method is used. It is worthy noting that the effectiveness and accuracy of simulation modeling results are verified by the real physical experiments. The novelty in this paper as compared to those referred to is that the physical parameters, stiffness and damping coefficients, were first estimated, and, finally, the model was simulated, and the results were compared with experimental results. For instance, in [2, 4, 11], the authors validated their models from numerical examples with known parameters.
2. Vehicle Crash Experimental Test
The real vehicle crash experiment was conducted on a typical mid-speed vehicle for pole collision. Its elaboration was the initiative of Robbersmyr . A test vehicle was subjected to impact with a vertical, rigid cylinder. The acceleration field was 100 meters long and had two anchored parallel pipelines. The vehicle was steered using those pipelines that were bolted to the concrete runaway. Setup scheme is shown in Figure 1.
During the test, the acceleration was measured in three directions (-longitudinal,-lateral, and-vertical) together with the yaw rate from the center of gravity of the car. Using normal-speed and high-speed video cameras, the behavior of the safety barrier and the test vehicle during the collision was recorded. The initial velocity of the car was 35 km/h, and the mass of the vehicle (together with the measuring equipment and dummy) was 873 kg. The obstruction was constructed with two steel components—a pipe filled with concrete and a baseplate mounted with bolts on a foundation. The car undergoing the deformation is shown in Figure 2. The accelerometer is located at the mass center of gravity of the vehicle in the passenger compartment. Since we are interested in the frontal crash, only the measured acceleration in the longitudinal direction is considered in this study. The acceleration data is imported and processed in MATLAB for analysis. The deformation of the vehicle is obtained by integrating twice the acceleration signal.
3. Mathematical Modeling Theoretical Background
Mathematical models describe the dynamic behavior of a system as a function of time. During frontal crash, the vehicle is subjected to an impulsive force caused by the obstacle. The model for vehicle crash simulates a rigid barrier impact of a vehicle, where and represent the frame rail (chassis) and passenger compartment masses, respectively.
Parameters to be estimated are springs and and dampers and , as shown in Figure 3. When the vehicle impacts on a rigid barrier, the two masses will experience an impulsive force during collision. The method for solving the impact responses of the two masses is adapted from the method used in the free vibration analysis of a two-degrees-of-freedom damped system .
The dynamic equations of the double-mass-spring-damper model are shown as follows: or The solution for , can be represented as with .
After substituting , , , and into (4), we get the follwing: Now, for a nontrivial response, that is, for nonzero values of and , the determinant of their coefficient matrix must vanish. That is, Expansion of (10) leads to a characteristic equation of the system obtained as follows: where
Equation (11) is a fourth-order polynomial in s and is to be solved to get four roots. All of the coefficients of this polynomial are physical parameters of the system shown in Figure 3 and are all positive. For that reason, such a polynomial cannot have positive roots. Three allowable configurations of roots are given as follows : (1)two pairs of complex conjugates, (2)one pair of complex conjugates and two real and negative roots, (3)four real and negative roots.
Case 1 (two pairs of complex conjugates). The system in this case has moderate damping. The rate of decay is defined by , the real part of the root, and the frequency of vibration is specified by , the imaginary part. The two pairs of complex conjugates are the following:(1),,(2),,where , , , and are all positive. and are the first pair of complex conjugates, and and are the second pair.
The two roots and in the first pair will yield the solutions and , where the first subscript refers to the mass index and the second subscript refers to the pair number of the complex conjugate. The displacement components and due to and , respectively, are given by where,,,,, .
The general solution is the following:
Case 2 (one pair of complex conjugates and two real and negative roots). The general displacement solutions are shown as follows: where .
Case 3 (four real and negative roots). The system has a large damping. When it is disturbed, the system will settle to its equilibrium configuration without oscillation. The solutions of the 4th-order polynomial yield four real and negative roots. In , the authors focused just on the first case. In this paper, we will also focus on the third case. The third case is for the system which has a large damping. When it is disturbed, the system will settle to its equilibrium configuration without oscillation. The displacement signal of the real crash is similar to that of a case of an overdamped vibrating system. Hence, the third case would represent the vehicle frontal crash reconstruction. The solution for Case 3 is the following: with , where are the roots of the characteristic equation (11).
3.1. Estimation of Model Parameters by Curve Fitting
When (16) is curve fitted into displacement experimental data, the constants and (, and ) can be easily found resulting in the following system of equations that can be solved for and :
For , (17) can be written in a matrix form as Let Then, (19) can be represented as and, using the pseudoinverse, we can obtain the following: The spring stiffness and damping coefficient are calculated from (25) as follows:
3.2. Eigensystem Realization Algorithm
In this section, a state-space representation (SSR) of the model is derived from the dynamic equation of a vehicle subjected to a frontal crash. The same representation is also retrieved from the system Markov parameters.
3.2.1. Formulation of the SSR from the Model Dynamic Equations
Considering an input force acting on the front mass , (1) can be rewritten as
The dynamic equation (26) can be rewritten in a matrix compact form as with In general, the equation of motion fordegrees of freedom is expressed in a matrix form as or, equivalently, where , , and are the mass, damping, and stiffness matrices, respectively, while is the input matrix.
, , andare vectors of generalized acceleration, velocity, and displacements, respectively, and the vectorof dimension is the input force containing external excitations acting on the systems.
Let us define the new state variables , , , and .
Substituting these variables into (26) and combining in state equations, we get the following: Using the original state variables and interchanging rows 2 and 3, columns 2 and 3 of (31) and we get the following:
The output equation or the measurement vector , which may contain any combination of modal displacements, velocities, and/or accelerations, is given by where , , and are the output influence matrices for position, velocity, and acceleration, respectively. In our experiment, the acceleration is measured. Therefore, the output equation is the acceleration measurement From (32),
Therefore, the continuous-time state-space model of the dynamic system is written as with where is the state matrix, is the input matrix or the state influence matrix, is the output matrix or the measurement influence matrix, and is the feedforward matrix or the direct transmission matrix.
Once the , , , and matrices are known, it is easy to find the transfer function (TF) and impulse response function (IRF) of the system. By using these matrices, the system’s response to any input can be found in time domain or frequency domain. The state-space representation is useful for constructing the mathematical model in MATLAB environment.
The discrete-time state-space representation of a MIMO system is given by the following: where is the integer discrete-time index at time instant , is the state vector at the discrete-time , is the force vector, is the output vector, is the discrete state system matrix, and is the discrete input influence matrix for the state vector . The output matrix and the direct transition matrix during the zero-order-hold operations. Because experimental data are discrete in nature, equations (38) form the basis for the system identification of linear time-invariant, dynamical systems. The state matrix and the influence matrix of the discrete-time model are related to the matrices , of the continuous-time model by the following expression : The continuous-time model is calculated from the discrete-time model by where is a constant interval.
The dimensions of the discrete-time system are equal to those of the continuous system.
3.2.2. ERA from the System Markov Parameters
ERA is a minimum-order realization technique that uses singular value decomposition technique. A flowchart for the ERA is shown in Figure 4.
If the excitations of the dynamic system are measured by the input quantities in the vector , the equations of motions and the set of output equations can both be, respectively, rewritten in terms of the state vector.
ERA begins with the definition of the Markov parameter of a state-space model. The method for deriving the expression for the system matrices is adapted from . Consider a discrete-time state-space model (38).
The state-space model (38) has an impulse response
The discrete-time Markov parameters can be defined in the same way as (41). The term of is called the Markov parameter of the system. By using these parameters, one can define the impulse response of the system as follows:
Consequently, the identification problem is the following: given the values of s, construct the constant matrices to identify the system.
The algorithm begins by constructing an generalized Hankel matrix.
Given a number of input and output measurements and generated by a system of unknown parameters, it is requested to identify the order of the system as well as the discrete state matrices () of the system. Then, the identified continuous state matrices () that have the same size as in the physical model can be estimated from (40), and, finally, extract the respective parameters.
All minimum realizations have the same set of eigenvalues and eigenvectors, which are the modal parameters of the system itself. Assume that the state matrix of order has a complete set of linearly independent eigenvectors with corresponding eigenvalues; then where is a diagonal matrix of the eigenvalues and is the matrix of the eigenvectors. The realization can be transformed into the realization by using the eigenvalues and eigenvectors matrices. The diagonal matrix contains the information of modal damping rates and damped natural frequencies. The matrix defines the initial modal amplitudes, and the matrix denotes the mode shapes at the sensor points.
All of the modal parameters of a dynamic system can, thus, be identified by the triplet .
The real part of , into the continuous-time model, gives the modal damping rates, while the imaginary part gives the damped natural frequencies.
After identifying the combined system and observer gain Markov parameters, the following step consists in forming the generalized Hankel matrix : When , we get the following:
In order to compute a minimum-order realization of the system , it is necessary to construct a shifted Hankel matrix as follows:
The block matrix is the observability matrix, whereas the block matrix is the controllability matrix.
Denote column submatrices of by and row submatrices of by . The ERA data block matrix can be expressed by Assume that there exists a matrix satisfying the relation where is an identity matrix of order. The matrix which is the pseudoinverse of plays a major role in deriving the ERA. It is observed that
The ERA process starts with the factorization of the block data matrix (46) using the singular value decomposition where the columns of matrices and are orthonormal and is a rectangular matrix given as with .
Let and be matrices formed by the firstcolumns of R and S, respectively. Hence, the matrices and become where .
Examining the singular value of the Hankel matrix , it is possible to determine the order of the system.
From (51), the first columns of the observability matrix form the input matrix whereas the first rows of the controllability matrix form the output matrix .
With in (51), we get that
One obvious solution for the state matrix becomes Let be a null matrix of order , let be an identity matrix of order , and let the matrices and be defined as where is the number of outputs and is the number of inputs.
Finally, using (50), (51), (52), (56), and (57), the basic formulation of the minimum-order realization for the ERA/OKID (OKID means Observer/Kalman filter identificationis) is given as follows: Recall that in (42) .
Hence, a minimal realization is given as follows: The continuous state matrix and the input influence matrix are obtained from (40), and the physical parameters-matrices , , and are embedded in the state matrix .
4. Results and Discussion
4.1. Parameters Estimation from Curve-Fitting Approach
From the curve fitting, the values for and () are found to be the following:,,,,, , ,.
The result from the curve fitting is shown in Figure 5.
4.2. Vehicle Crash Experimental Data Analysis
It is observed from Figure 6 that the dynamic crash from the real vehicle crash test is 53.17 cm and occurs at time , when the unfiltered data are used in the analysis. The filtered data result in a dynamic crash of 51.11 cm at time as shown in Figure 7.
The initial velocity for both filtered and unfiltered data is closer to 35 km/h (i.e., 34.99 km/h for the unfiltered data and 35.28 km/h for the filtered data).
4.3. Results from the Model
Four different cases are considered in this section as a sample of results. Let be the mass of the chassis, the mass of passenger compartment, and the total mass of the vehicle.
Case 1 (). One has that, and . From Figure 8, the dynamic crash of is 80 cm which is the displacement of the passenger compartment. Therefore, this model cannot represent the vehicle crash scenario. It is observed that the time for dynamic crash is longer than that for the real crash (i.e., 0.17 s instead of 0.078 s). The dynamic crash of is 42.5 cm and occurs after 0.17 s.
Case 2 (). One has that, and . From Figure 9, the dynamic crash of the passenger compartment is 66.2 cm. The time for dynamic crash increases further up to 0.15 s. The dynamic crash of is 45.5 cm and occurs after 0.13 s. Therefore, for this case, the model cannot represent the vehicle crash scenario.
Case 3 (). One has that, and . From Figure 10, the dynamic crash of the passenger compartment is 69 cm. The time for dynamic crash is 0.15 s. The dynamic crash of is 30.9 cm and occurs after 0.14 s. This also cannot represent the real vehicle crash test.
Case 4 (). One has that, and . From Figure 11, the dynamic crash of the passenger compartment is 49.8 cm, and the time for dynamic crash is 0.11 s. The dynamic crash of is 35.5 cm and occurs after 0.1 s. Therefore, this case can represent the vehicle crash senario because the dynamic crash is much closer to that from the real vehicle crash and the time is relatively small as compared to other cases.
A summary of main results is shown in Table 1. The values for and are the first and second entries of vector in (22). The values for and depend on the value of mass taken into consideration. Values for and are obtained from vector in (25). The stiffness coefficients which result in a closer vehicle crash reconstruction are found to be , and , and the damping coefficients are: and when the mass of the chassis is the total mass of the vehicle , where the dynamic crash of the passenger compartment is equal to 49.8 cm and occurs after 0.11 s (see Subsection 4.3, Case 4, Figure 11).
It is observed that the passenger compartment is the one that reconstructs the vehicle crash. When the mass of the chassis is greater than that of the passenger compartment, the results from the model are closer to the expected values. For example, when (i.e., of the total mass of the vehicle) and (i.e., of the mass of the vehicle), the dynamic crash of the passenger compartment is 49.8 cm which is closer to 51.11 cm (the dynamic crash from the real vehicle crash).
Remark 2. It is worthy noting that optimal values for stiffness and damping coefficients are not fixed. They are dependent on the mass of passenger compartment taken into consideration.
4.4. State-Space Realization of the System by ERA
Consider a 2nd-order system (i.e., ) for a single degree of freedom and , the number of samples to assemble the Hankel matrix. It is observed that the continuous system matrices from the ERA are the following: From the state matrix, the eigenvectors and eigenvalues were found to be and .
The natural frequency and the damping ratio of a single degree of freedom are 533 rad/s and 0.0338, respectively. The Bode diagram for the system is shown in Figure 12.
Considering a 4th-order system (i.e, ) for a two-degree-of-freedom system and the number of samples to assemble the Hankel matrix equal to , it is observed that the continuous state-space realizations from the ERA are found to be the following: with and .
The natural frequencies and damping ratios are given as= 3592.8 rad/s, = 9337.4 rad/s and = 0.481, = 0.112.
For higher sample size (i.e., with order ), the system matrices are given as follows: with and .
The natural frequencies and damping ratios are given as = 617 rad/s, = 575 rad/s and = 0.00712, = 0.006680.
Remark 3. Based on the state-space realization obtained by ERA, we should be able to transfer the state matrix to the original obtained from (36). Hence, extract the mass, stiffness, and damping matrices.
5. Conclusion and Future Work
In this paper, a method to estimate the parameters of a double-spring-mass-damper model of a vehicle frontal crash is presented. It was observed that the model results were much closer to the experimental test results. The overall behavior of the model matches the real vehicle’s crash. Two of the main parameters characterizing the collision are the maximum dynamic crash—which describes the highest car’s deformation—and the time at which it occurs—. They are pertinent to the occupant crashworthiness since they help to assess the maximum intrusion to the passenger’s compartment.
It can be concluded from this study that a double-spring-mass-damper model can represent the real vehicle crash scenario when the mass representing the chassis of the vehicle is less than , the mass representing the passenger compartment. In all cases, the front part of the vehicle undergoes smaller deformation than the passenger compartment. The time at which the maximum chassis displacement occurs is slightly shorter than the time for the passenger compartment because of its additional compression by the rest of the car. In our future work, more investigation on data analysis would be interesting to consider effect of sensor delays in measurement [21–23]; extraction of mass, stiffness, and damping matrices from the identified state-space representation will be investigated; we would like to extend our study to a three-mass-spring-damper model taking into consideration the nonlinearity of the spring and damper: the three masses will be representing the engine, the suspension, and the passenger compartment mass, respectively, interconnected by springs and dampers; injury mechanisms of the occupant such as Head Injury Criterion (HIC) would be considered in the future work.
J. Marzbanrad and M. Pahlavani, “Paremaeter determination of a vehicle 5-DOF model to simulate occupant deceleration in a frontal crash,” World Academy of Science, Engineering and Technology, vol. 79, no. 55, pp. 336–341, 2011.View at: Google Scholar
A. A. Alnaqi and A. S. Yigit, “Dynamic analysis and control of automotive occupant restraint systems,” Jordan Journal of Mechanical and Industrial Engineering, vol. 5, no. 1, pp. 39–46, 2011.View at: Google Scholar
D. Guida, F. Nilvetti, and C. M. Pappalardo, “Parameter identification of a two degrees of freedom mechanical system,” International Journal of Mechanics, vol. 3, no. 2, pp. 23–30, 2009.View at: Google Scholar
J. N. Juang, Applied System Identification, PTR Printice Hall, 6th edition, 1994.
K. G. Robbersmyr, “Calibration test of a standard Ford Fiesta 1.1l, model year 1987, according to NS-EN 12767,” Tech. Rep. 43/2004, Agder Research, Grimstad, Norway, 2004.View at: Google Scholar
M. Huang, Vehicle Crash Mechanics, CRC Press, Boca Raton, Fla, USA, 2002.