Mathematical Modeling and Parameters Estimation of Car Crash Using Eigensystem Realization Algorithm and Curve-Fitting Approaches

Abstract

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.

1. Introduction

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 [1], 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 [2].

In [3], 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 [8], 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 [14], Chiang and Lin [15], Ko and Hung [16], Juang and Pappa [17], and Yang and Yeh [18]. For instance, in [17], the author developed the ERA to estimate the natural frequencies and damping ratios of a dynamical system from known Markov parameters, and in [18], 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 [19]. 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.

Figure 1 

Vehicle crash experimental setup [19].

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.

Figure 2 

Vehicle undergoing deformation [19].

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 [20].

Figure 3 

A double-spring-mass-damper model.

The dynamic equations of the double-mass-spring-damper model are shown as follows: 𝑚₁𝑥̈₁+𝑐₁+𝑐₂𝑥̇₁+𝑘₁+𝑘₂𝑥₁−𝑐₂𝑥̇₂−𝑘₂𝑥₂𝑚=0,₂𝑥̈₂−𝑐₂𝑥̇₁+𝑐₂𝑥̇₂+𝑘₂𝑥₂−𝑘₂𝑥₁=0, (1) or 𝑚₁00𝑚₂𝑥̈₁𝑥̈₂+𝑘₁+𝑘₂−𝑘₂−𝑘₂𝑘₂𝑥₁𝑥₂+𝑐₁+𝑐₂−𝑐₂−𝑐₂𝑐₂𝑥̇₁𝑥̇₂=00. (2) The solution for 𝑥_𝑖,𝑖=1,2, can be represented as 𝑥_𝑖=𝐶_𝑖𝑒^𝑠𝑘𝑡 (3) with 𝑘=1,…,4.

Where 𝐶_𝑖 and 𝑠_𝑘 may be complex numbers. Substituting (3) into (1), we get 𝑅_1𝑘𝐶₁−𝑅_2𝑘𝐶₂=0,−𝑆_1𝑘𝐶₁+𝑆_2𝑘𝐶₂𝐶=0,₂𝐶₁=𝑅_1𝑘𝑅_2𝑘=𝑆_1𝑘𝑆_2𝑘,𝑅_1𝑘=𝑅_2𝑘∗𝑆_1𝑘𝑆_2𝑘,𝑅_1𝑘∗𝑆_2𝑘−𝑅_2𝑘∗𝑆_1𝑘=0, (4)(5)(6)(7) where 𝑅_1𝑘=𝑚₁𝑠²_𝑘+𝑐₁+𝑐₂𝑠_𝑘+𝑘₁+𝑘₂,𝑅_2𝑘=𝑐₂𝑠_𝑘+𝑘₂,𝑆_1𝑘=𝑐₂𝑠_𝑘+𝑘₂,𝑆_2𝑘=𝑚₂𝑠²_𝑘+𝑐₂𝑠_𝑘+𝑘₂. (8)

After substituting 𝑅_1𝑘, 𝑅_2𝑘, 𝑆_1𝑘, and 𝑆_2𝑘 into (4), we get the follwing: 𝑚₁𝑠²_𝑘+𝑐₁+𝑐₂𝑠_𝑘+𝑘₁+𝑘₂−𝑐₂𝑠_𝑘−𝑘₂−𝑐₂𝑠_𝑘−𝑘₂𝑚₂𝑠²_𝑘+𝑐₂𝑠_𝑘+𝑘₂𝐶₁𝐶₂=00. (9) Now, for a nontrivial response, that is, for nonzero values of 𝐶₁ and 𝐶₂, the determinant of their coefficient matrix must vanish. That is, 𝑚₁𝑠²_𝑘+𝑐₁+𝑐₂𝑠_𝑘+𝑘₁+𝑘₂𝑚₂𝑠²_𝑘+𝑐₂𝑠_𝑘+𝑘₂+𝑐₂𝑠_𝑘+𝑘₂²=0. (10) Expansion of (10) leads to a characteristic equation of the system obtained as follows: 𝑠⁴_𝑘+𝑡𝑠³_𝑘+𝑢𝑠²_𝑘+𝑣𝑠_𝑘+𝑤=0, (11) where 𝑚𝑡=¹𝑐₂+𝑚₂𝑐₁+𝑐₂𝑚₁𝑚₂,𝑚𝑢=¹𝑘₂+𝑚₂𝑘₁+𝑘₂+𝑐₁𝑐₂𝑚₁𝑚₂,𝑘𝑣=¹𝑐₂+𝑘₂𝑐₁𝑚₁𝑚₂𝑘,𝑤=¹𝑘₂𝑚₁𝑚₂. (12)

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 [20]: (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 𝑋₁₁=𝐴₁₁𝑒^−𝑝1𝑡×sin𝑞₁𝑡+𝜙₁₁,𝑋₂₁=𝐴₂₁𝑒^−𝑝1𝑡×sin𝑞₁𝑡+𝜙₂₁, (13) where𝐴²₁₁=4𝐶₁₁𝐶₁₂,𝐴₂₁=4𝐶₂₁𝐶₂₂,𝐶₁₁=𝑚₂𝑠²₁+𝑐₂𝑠₁+𝑘₂,𝐶₂₁=𝑐₂𝑠₁+𝑘₂,𝐶₁₂=𝑚₂𝑠²₂+𝑐₂𝑠₂+𝑘₂, 𝐶₂₂=𝑐₂𝑠₂+𝑘₂.
The general solution is the following: 𝑋_𝑖=²_𝑗=1𝑥_𝑖𝑗=𝐴_𝑖1𝑒^−𝑝1𝑡sin𝑞₁𝑡+𝜙_𝑖1+𝐴_𝑖2𝑒^−𝑝2𝑡sin𝑞₂𝑡+𝜙_𝑖2. (14)

Case 2 (one pair of complex conjugates and two real and negative roots). The general displacement solutions are shown as follows: 𝑋_𝑖=𝐴_𝑖1𝑒^−𝑝𝑡sin𝑞𝑡+𝜙_𝑖1+𝐶_𝑖3𝑒^𝑠3𝑡+𝐶_𝑖4𝑒^𝑠4𝑡, (15) where 𝑖=1,2.

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 [2], 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: 𝑋_𝑖(𝑡)=𝐶_𝑖1𝑒^𝑠1𝑡+𝐶_𝑖2𝑒^𝑠2𝑡+𝐶_𝑖3𝑒^𝑠3𝑡+𝐶_𝑖4𝑒^𝑠4𝑡 (16) with 𝑖=1,2, 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 𝑠_𝑘 (𝑖=1,2, and 𝑘=1,…,4) can be easily found resulting in the following system of equations that can be solved for 𝑘₂ and 𝑐₂: 𝐶_1𝑘=𝑚₂𝑠²_𝑘+𝑐₂𝑠_𝑘+𝑘₂,𝐶_2𝑘=𝑐₂𝑠_𝑘+𝑘₂. (17)(18)

For 𝑘=1,…,4, (17) can be written in a matrix form as 1𝑠₁1𝑠₂1𝑠₃1𝑠₄𝑘₂𝑐₂=𝐶₁₁−𝑚₂𝑠²₁𝐶₁₂−𝑚₂𝑠²₂𝐶₁₃−𝑚₂𝑠²₃𝐶₁₄−𝑚₂𝑠²₄. (19) Let 𝐀=1𝑠₁1𝑠₂1𝑠₃1𝑠₄,𝐁₂=𝐶₁₁−𝑚₂𝑠²₁𝐶₁₂−𝑚₂𝑠²₂𝐶₁₃−𝑚₂𝑠²₃𝐶₁₄−𝑚₂𝑠²₄,𝐕₂=𝑘₂𝑐₂. (20) Then, (19) can be represented as 𝐀∗𝐕₂=𝐁₂, (21) and, using the pseudoinverse, we can obtain the following: 𝐕₂=𝐀^𝑇∗𝐀⁻¹∗𝐁₂. (22) The spring stiffness 𝑘₁ and damping coefficient 𝑐₁ are calculated from (25) as follows: 1𝑠₁1𝑠₂1𝑠₃1𝑠₄𝑘₁𝑐₁=𝑅₁₁−𝑚₁𝑠²₁+𝑐₂𝑠₁+𝑘₂𝑅₁₂−𝑚₁𝑠²₂+𝑐₂𝑠₂+𝑘₂𝑅₁₃−𝑚₁𝑠²₃+𝑐₂𝑠₃+𝑘₂𝑅₁₄−𝑚₁𝑠²₄+𝑐₂𝑠₄+𝑘₂. (23)

Remark 1. From (6), it was shown that 𝑅_1𝑘=𝑅_2𝑘∗𝑆_1𝑘/𝑆_2𝑘. Let 𝐕₁=𝑘₁𝑐₁,𝐁₁=𝑅₁₁−𝑚₁𝑠²₁+𝑐₂𝑠₁+𝑘₂𝑅₁₂−𝑚₁𝑠²₂+𝑐₂𝑠₂+𝑘₂𝑅₁₃−𝑚₁𝑠²₃+𝑐₂𝑠₃+𝑘₂𝑅₁₄−𝑚₁𝑠²₄+𝑐₂𝑠₄+𝑘₂. (24) Then, we have that 𝐕₁=𝐀^𝑇∗𝐀⁻¹∗𝐁₁. (25) Therefore, the estimated parameters are obtained from (22) and (25).

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 𝑥̈₁𝑘=−¹+𝑘₂𝑚₁𝑥₁+𝑘₂𝑚₁𝑥₂−𝑐₁+𝑐₂𝑚₁𝑥̇₁+𝑐₂𝑚₁𝑥̇₂+1𝑚₁𝐮₁,𝑥̈₂=𝑘₂𝑚₂𝑥₁−𝑘₂𝑚₂𝑥₂+𝑐₂𝑚₂𝑥̇₁−𝑐₂𝑚₂𝑥̇₂. (26)

The dynamic equation (26) can be rewritten in a matrix compact form as 𝐌𝑥̈₁𝑥̈₂𝑥̇+𝐋¹𝑥̇₂𝑥+𝐊¹𝑥₂=𝐮₁0 (27) with 𝑚𝐌=¹00𝑚₂𝑐,𝐋=¹+𝑐₂−𝑐₂−𝑐₂𝑐₂,𝑘𝐊=¹+𝑘₂−𝑘₂−𝑘₂𝑘₂. (28) In general, the equation of motion for𝑁degrees of freedom is expressed in a matrix form as 𝐌𝐱̈+𝐋𝐱̇+𝐊𝐱=𝐁_𝑐𝐮, (29) or, equivalently, 𝐱̈=−𝐌⁻¹𝐋𝐱̇−𝐌⁻¹𝐊𝐱+𝐌⁻¹𝐁_𝑐𝐮, (30) where 𝐌∈ℝ^𝑁×𝑁, 𝐋∈ℝ^𝑁×𝑁, and 𝐊∈ℝ^𝑁×𝑁 are the mass, damping, and stiffness matrices, respectively, while 𝐁_𝑐∈ℝ^𝑁×𝑟 is the input matrix.

𝐱̈, 𝐱̇, and𝐱are vectors of generalized acceleration, velocity, and displacements, respectively, and the vector𝐮of dimension 𝑟×1 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: 𝐱̇₁𝐱̇₂𝐱̇₃𝐱̇₄=−𝑘0100₁+𝑘₂𝑚₁−𝑐₁+𝑐₂𝑚₁𝑘₂𝑚₁𝑐₂𝑚₁𝑘0001₂𝑚₂𝑐₂𝑚₂−𝑘₂𝑚₂−𝑐₂𝑚₂𝐱₁𝐱₂𝐱₃𝐱₄+01𝑚₁00𝐮. (31) Using the original state variables and interchanging rows 2 and 3, columns 2 and 3 of (31) and we get the following: 𝑥̇₁𝑥̇₂𝑥̈₁𝑥̈₂=−𝑘00100001₁+𝑘₂𝑚₁𝑘₂𝑚₁−𝑐₁+𝑐₂𝑚₁𝑐₂𝑚₁𝑘₂𝑚₂−𝑘₂𝑚₂𝑐₂𝑚₂−𝑐₂𝑚₂𝑥₁𝑥₂𝑥̇₁𝑥̇₂+001𝑚₁0𝐮. (32)

The output equation or the measurement vector 𝐲(𝑡), which may contain any combination of modal displacements, velocities, and/or accelerations, is given by 𝐂𝐲=^𝑝𝑥𝐂_𝑣𝑥̇𝐂_𝑎𝑥̈, (33) 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 𝐲=𝐂_𝑎𝑥̈₁. (34) From (32), 𝐲=𝑥̈₁=−𝑘₁+𝑘₂𝑚₁𝑘₂𝑚₁−𝑐₁+𝑐₂𝑚₁𝑐₂𝑚₁𝑥₁𝑥₂𝑥̇₁𝑥̇₂+1𝑚₁𝐮. (35)

Therefore, the continuous-time state-space model of the dynamic system is written as 𝐱̇=𝐀_𝑐𝐱+𝐁_𝑐𝐮,𝐲=𝐂_𝑐𝐱+𝐃_𝑐𝐮 (36) with 𝐀_𝑐=−𝑘00100001₁+𝑘₂𝑚₁𝑘₂𝑚₁−𝑐₁+𝑐₂𝑚₁𝑐₂𝑚₁𝑘₂𝑚₂−𝑘₂𝑚₂𝑐₂𝑚₂−𝑐₂𝑚₂,𝐁_𝑐=001𝑚₁0,𝐂_𝑐=−𝑘₁+𝑘₂𝑚₁𝑘₂𝑚₁−𝑐₁+𝑐₂𝑚₁𝑐₂𝑚₁,𝐃_𝑐=1𝑚₁, (37) 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: 𝐱()=(𝑘)+(𝑘),𝐲(𝑘)=(𝑘)+(𝑘),𝑘+1𝐀𝐱𝐁𝐮𝐂𝐱𝐃𝐮 (38) 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 [14]: 𝐀=𝑒^{𝐀𝑐Δ𝑡},𝐁=₀^Δ𝑡𝑒^𝐀𝑐𝜏𝑑𝜏⋅𝐁_𝑐. (39) The continuous-time model is calculated from the discrete-time model by 𝐀_𝑐=ln𝐀,𝐁Δ𝑡_𝑐=𝐀_𝑐𝐀−𝐼⁻¹∗𝐁, (40) 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.

Figure 4 

Flowchart for the ERA by Juang [14].

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 [14]. Consider a discrete-time state-space model (38).

The state-space model (38) has an impulse response 𝐲(𝑘)=𝐶𝐀𝐷,𝑘=0,_𝑘−1𝐁,𝑘≥1. (41)

The discrete-time Markov parameters can be defined in the same way as (41). The term of 𝐶𝐀^𝑘−1𝐁 is called the Markov parameter of the system. By using these parameters, one can define the impulse response of the system as follows: 𝑌[𝑘]=𝐶𝐀^𝑘−1𝑌[0]=[1]=𝐶𝑌[2]=𝐶𝐀[𝑘]=𝐶𝐀𝐁,𝐷,𝑌𝐵,𝐵,…,𝑌^𝑘−1𝐁. (42)(43)

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 𝐀Ψ=ΨΛ, (44) 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 𝐻_(𝑘−1): 𝐻_(𝑘−1)=𝑌_𝑘𝑌_(𝑘+1)⋯𝑌_{(𝑘+𝛽+1)}𝑌_𝑘+1𝑌_(𝑘+2)⋯𝑌_(𝑘+𝛽)𝑌⋯⋯⋯⋯_{(𝑘+𝛼−1)}𝑌_(𝑘+𝛼)⋯𝑌_{(𝑘+𝛼+𝛽−2)}. (45) When 𝑘=1, we get the following: 𝐻₍₀₎=𝑌₁𝑌₂⋯𝑌_𝛽𝑌₂𝑌₃⋯𝑌_(𝛽+1)𝑌⋯⋯⋯⋯_𝛼𝑌_𝛼+1⋯𝑌_{𝛼+𝛽−1}. (46)

In order to compute a minimum-order realization of the system {𝐀,𝐁,𝐂}, it is necessary to construct a shifted Hankel matrix 𝐻₍₁₎ as follows: 𝐻₍₁₎=𝑌₂𝑌₃⋯𝑌_𝛽+1𝑌₃𝑌₄⋯𝑌_𝛽+2𝑌⋯⋯⋯⋯_(𝛼+1)𝑌_(𝛼+2)⋯𝑌_(𝛼+𝛽). (47)

Substituting the Markov parameters from (43) into (45) and decomposing 𝐻(𝑘−1) into three matrices yield the following: 𝐻(𝑘−1)=𝐎_𝑏𝐴^𝑘−1𝐂_𝑡, (48) where 𝐎_𝑏 and 𝐂_𝑡 are given as 𝐎_𝑏=𝐶𝐶𝐴𝐶𝐴²⋮𝐶𝐴^𝛼−1,𝐂_𝑡=𝐵𝐴𝐵𝐴²𝐵⋯𝐴^𝛽−1𝐵. (49)

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 𝐻_(𝑘−1)=𝑌_{𝑠𝑖+𝑘+𝑡𝑗},𝑌_{𝑠𝑖+𝑘+𝑡𝑗}=𝐶_𝑗𝐴^{𝑠𝑖+𝑘−1+𝑡𝑗}𝐵_𝑖,𝐻(𝑘)=𝐎_𝑏𝐴^𝑘𝐂_𝑡. (50)(51) Assume that there exists a matrix 𝐻^† satisfying the relation 𝐎_𝑏𝐻^†𝐂_𝑡=𝐼_𝑛, (52) where 𝐼_𝑛 is an identity matrix of order𝑛. The matrix 𝐻^† which is the pseudoinverse of 𝐻(0) plays a major role in deriving the ERA. It is observed that 𝐻(0)𝐻^†𝐻(0)=𝐎_𝑏𝐂_𝑡𝐻^†𝐎_𝑏𝐂_𝑡=𝐎_𝑏𝐂_𝑡=𝐻(0). (53)

The ERA process starts with the factorization of the block data matrix (46) using the singular value decomposition 𝐻₍₀₎=𝑅𝑆^𝑇, (54) where the columns of matrices 𝑅 and 𝑆 are orthonormal and ∑ is a rectangular matrix given as∑=^𝑛000 (55) with ∑_𝑛=diag{𝜎₁,𝜎₂,…,𝜎_𝑖,𝜎_𝑖+1,…𝜎_𝑛}.

Let 𝑅_𝑛 and ∑_𝑛𝑆_𝑛 be matrices formed by the first𝑛columns of R and S, respectively. Hence, the matrices𝐻(0) and 𝐻^† become 𝐻₍₀₎=𝑅_𝑛𝑛𝑆^𝑇_𝑛,𝐻^†=𝑆_𝑛⁻¹_𝑛𝑅^𝑇_𝑛, (56)(57) where 𝑅^𝑇_𝑛𝑅_𝑛=𝐼_𝑛=𝑆^𝑇_𝑛𝑆_𝑛.

Examining the singular value ∑_𝑛 of the Hankel matrix 𝐻₍₀₎, it is possible to determine the order of the system.

Comparing (56) and (51) with 𝑘=0, we get that 𝐎_𝑏=𝑅_𝑛∑_𝑛^1/2 and 𝐂_𝑡=∑_𝑛^1/2𝑆^𝑇_𝑛.

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 𝑘=1 in (51), we get that 𝐻₍₁₎=𝐎_𝑏𝐴𝐂_𝑡=𝑅_𝑛^1/2_𝑛𝐴^1/2_𝑛𝑆^𝑇_𝑛. (58)

One obvious solution for the state matrix 𝐴 becomes 𝐴=^−1/2_𝑛𝑅^𝑇_𝑛𝐻(1)𝑆_𝑛^−1/2_𝑛. (59) Let 𝐎_𝑖 be a null matrix of order 𝑖, let 𝐼_𝑖 be an identity matrix of order 𝑖, and let the matrices 𝐄^𝑇_𝑚 and 𝐄^𝑇_𝑟 be defined as 𝐄^𝑇_𝑚=𝐈_𝑚𝟎_𝑚⋯𝟎_𝑚,𝐄^𝑇_𝑟=𝐈_𝑟𝟎_𝑟⋯𝟎_𝑟, (60) 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: 𝑌_𝑘=𝐸^𝑇_𝑚𝑅_𝑛^−1/2_𝑛^−1/2_𝑛𝑅^𝑇_𝑛𝐻(1)𝑆_𝑛^−1/2_𝑛^{𝑘−1−1/2}_𝑛𝑆^𝑇_𝑛𝐸_𝑟. (61) Recall that in (42) 𝑌[𝑘]=𝐶𝐀^𝑘−1𝐁.

Hence, a minimal realization is given as follows: 𝐀=^−1/2_𝑛𝐑^𝑇_𝑛𝐇₍₁₎𝑆_𝑛^−1/2_𝑛,𝐁=^−1/2_𝑛𝐒^𝑇_𝑛𝐄_𝑟,𝐂=𝐄^𝑇_𝑚𝐑_𝑛^−1/2_𝑛,𝐃=𝐘₍₀₎. (62) 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 𝐶_1𝑗 (𝑗=1,…,4) are found to be the following:𝐶₁₁=−253.6,𝐶₁₂=83.75,𝐶₁₃=333,𝐶₁₄=−163.2,𝑠₁=−7.7540, 𝑠₂=−6.238, 𝑠₃=−9.603,𝑠₄=−10.87.

The result from the curve fitting is shown in Figure 5.

Figure 5 

Result from curve fitting.

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 𝑡_𝑐=0.078ms, when the unfiltered data are used in the analysis. The filtered data result in a dynamic crash of 51.11 cm at time 𝑡_𝑐=74.5ms as shown in Figure 7.

Figure 6 

Unfiltered data plot from the real vehicle crash test.

Figure 7 

Filtered data plot from the real vehicle crash test.

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 𝑚_𝑡=873kg the total mass of the vehicle.

Case 1 (𝑚₁<𝑚₂). One has that𝑚₁=(1/3)𝑚_𝑡, and  𝑚₂=(2/3)𝑚_𝑡. 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.

Figure 8 

Comparative analysis between vehicle crash test and model results for 𝑚1=(1/3)𝑚𝑡.

Case 2 (𝑚₁>𝑚₂). One has that𝑚₁=(2/3)𝑚_𝑡, and 𝑚₂=(1/3)𝑚_𝑡. 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.

Figure 9 

Comparative analysis between vehicle crash test and model results for 𝑚2=(1/3)𝑚𝑡.

Case 3 (𝑚₁<𝑚₂). One has that𝑚₁=(1/4)𝑚_𝑡, and 𝑚₂=(3/4)𝑚_𝑡. 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.

Figure 10 

Comparative analysis between Vehicle crash test and model results for 𝑚1=(1/4)𝑚𝑡.

Case 4 (𝑚₁>𝑚₂). One has that𝑚₁=(3/4)𝑚_𝑡, and 𝑚₂=(1/4)𝑚_𝑡. 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.

Figure 11 

Comparative analysis between vehicle crash test and model results for 𝑚2=(1/4)𝑚𝑡.

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 𝑘₁=74681N/m, and 𝑘₂=45821N/m, and the damping coefficients are: 𝑐₁=18176Ns/m and 𝑐₂=11196Ns/m when the mass of the chassis is 3/4 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).

Table 1 

Parameters estimation.

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 𝑚₁=582kg (i.e., 3/4 of the total mass of the vehicle) and 𝑚₂=291kg (i.e., 1/4 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., 𝑛=2) for a single degree of freedom and 𝑁=100, the number of samples to assemble the Hankel matrix. It is observed that the continuous system matrices from the ERA are the following: 𝐀_𝑐=95.5776−545.1246545.12−131.6485,𝐁_𝑐=,−1474.85798.1𝐂=0.16390.5717,𝐷=0.1964. (63) From the state matrix𝐀_𝑐, the eigenvectors Ψ and eigenvalues Λ were found to be Ψ,Ψ^∗=0.1474±0.6916𝑖0.7071, (64) and Λ=diag{−18.04±533.15𝑖}.

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.

(a)

(b)

Figure 12 

The Bode diagram for 𝑛=2, 𝑁=100, and 𝑁=50.

Considering a 4th-order system (i.e, 𝑛=4) for a two-degree-of-freedom system and the number of samples to assemble the Hankel matrix equal to 𝑁=60, it is observed that the continuous state-space realizations from the ERA are found to be the following: 𝐀_𝑐=,𝐁−22.5−526.356.2141.5526.3−243.8631.5279.956.2−631.5−171.2−1363.2−141.5279.9136.2−445.9_𝑐=,,−1990.85436.83302.4−4233.2𝐂=0.21490.5366−0.3368−0.3826𝐷=0.1964 (65) with Ψ₁,Ψ^∗₁=,Ψ0.68640.2979∓0.5573𝑖−0.1126±0.1269𝑖−0.0715±0.3092𝑖₂,Ψ^∗₂=,0.0537∓0.0828𝑖−0.1269∓0.3363𝑖0.70020.0587∓0.6061𝑖 (66) and Λ=diag{−274.9±501.4𝑖,−166.8±1476.7𝑖}.

The natural frequencies and damping ratios are given as𝜔_𝑛1= 3592.8  rad/s, 𝜔_𝑛2 = 9337.4 rad/s and 𝜁₁ = 0.481, 𝜁₂ = 0.112.

One can see the stability of the model from pole-zero map as shown in Figure 13 and Figure 14. It is worth noted that for small sample size the system becomes unstable as evidenced by Figure 15.

Figure 13 

Pole-zero map for 𝑛=2 and 𝑁=100.

Figure 14 

Pole-zero map for 𝑛=2 and 𝑁=50.

Figure 15 

Responses of a single-degree-of-Freedom model to excitation force.

For higher sample size (i.e., 𝑁=2501 with order 𝑛=4), the system matrices are given as follows: 𝐀_𝑐=,𝐁0.1−615.1−3.6−4.9615.1−13.1−24.9−5.3−3.624.9−2.1−577.0−4.9−5.3336577.0−0.9876_𝑐=,,1280.1−5366.11636.4−946.8𝐂=−0.144−0.5322−0.1656−0.0897𝐷=0.1964 (67) with Ψ₁,Ψ^∗₁=,Ψ0.0053±0.6800𝑖0.6840−0.0869∓0.1733𝑖−0.1561±0.0880𝑖₂,Ψ^∗₂=,0.0878±0.1732𝑖0.1572∓0.0850𝑖−0.0027±0.6799𝑖0.6841 (68) and Λ=diag{−4.39±617𝑖,−3.85±575𝑖}.

The natural frequencies and damping ratios are given as 𝜔_𝑛1 = 617 rad/s, 𝜔_𝑛2 = 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, (1) more investigation on data analysis would be interesting to consider effect of sensor delays in measurement [21–23]; (2) extraction of mass, stiffness, and damping matrices from the identified state-space representation will be investigated; (3) 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; (4) injury mechanisms of the occupant such as Head Injury Criterion (HIC) would be considered in the future work.