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

Authors: Bernard B. Munyazikwiye; Hamid Reza Karimi; Kjell Gunnar Robbersmyr
Journal: Mathematical Problems in Engineering
 (2013)

---

## 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.


---

## Body

## 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 (x-longitudinal,y-lateral, andz-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, wherem1 and m2 represent the frame rail (chassis) and passenger compartment masses, respectively.Parameters to be estimated are springsk1 and k2 and dampers c1 and c2, 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:(1)m1x¨1+(c1+c2)x˙1+(k1+k2)x1-c2x˙2-k2x2=0,m2x¨2-c2x˙1+c2x˙2+k2x2-k2x1=0,
or
(2)[m100m2][x¨1x¨2]+[k1+k2-k2-k2k2][x1x2]+[c1+c2-c2-c2c2][x˙1x˙2]=[00].
The solution for xi,i=1,2, can be represented as
(3)xi=Cieskt
with k=1,…,4.WhereCi and sk may be complex numbers. Substituting (3) into (1), we get
(4)R1kC1-R2kC2=0,-S1kC1+S2kC2=0,(5)C2C1=R1kR2k=S1kS2k,(6)R1k=R2k*S1kS2k,(7)R1k*S2k-R2k*S1k=0,
where
(8)R1k=m1sk2+(c1+c2)sk+(k1+k2),R2k=c2sk+k2,S1k=c2sk+k2,S2k=m2sk2+c2sk+k2.After substitutingR1k, R2k, S1k, and S2k into (4), we get the follwing:
(9)[m1sk2+(c1+c2)sk+k1+k2-c2sk-k2-c2sk-k2m2sk2+c2sk+k2][C1C2]=[00].
Now, for a nontrivial response, that is, for nonzero values of C1 and C2, the determinant of their coefficient matrix must vanish. That is,
(10)[m1sk2+(c1+c2)sk+k1+k2][m2sk2+c2sk+k2]+(c2sk+k2)2=0.
Expansion of (10) leads to a characteristic equation of the system obtained as follows:
(11)sk4+tsk3+usk2+vsk+w=0,
where
(12)t=m1c2+m2(c1+c2)m1m2,u=m1k2+m2(k1+k2)+c1c2m1m2,v=k1c2+k2c1m1m2,w=k1k2m1m2.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 byp1, the real part of the root, and the frequency of vibration is specified by q1, the imaginary part. The two pairs of complex conjugates are the following:(1)
s1=-p1+iq1,s2=-p1-iq1,(2)
s3=-p2+iq2,s4=-p2-iq2,where p1, p2, q1, and q2 are all positive. s1 and s2 are the first pair of complex conjugates, and s3 and s4 are the second pair.
The two rootss1 and s2 in the first pair will yield the solutions X11 and X21, 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 X11 and X21 due to s1 and s2, respectively, are given by
(13)X11=A11e-p1t×sin(q1t+ϕ11),X21=A21e-p1t×sin(q1t+ϕ21),
whereA112=4C11C12,A21=4C21C22,C11=m2s12+c2s1+k2,C21=c2s1+k2,C12=m2s22+c2s2+k2, C22=c2s2+k2.
The general solution is the following:(14)Xi=∑j=12‍xijXi=Ai1e-p1tsin(q1t+ϕi1)+Ai2e-p2tsin(q2t+ϕi2).Case 2 (one pair of complex conjugates and two real and negative roots).
The general displacement solutions are shown as follows:(15)Xi=Ai1e-ptsin(qt+ϕi1)+Ci3es3t+Ci4es4t,
where i=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:
(16)Xi(t)=Ci1es1t+Ci2es2t+Ci3es3t+Ci4es4t
with i=1,2, where si 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 Cik and sk (i=1,2, and k=1,…,4) can be easily found resulting in the following system of equations that can be solved for k2 and c2:
(17)C1k=m2sk2+c2sk+k2,(18)C2k=c2sk+k2.Fork=1,…,4, (17) can be written in a matrix form as
(19)[1s11s21s31s4][k2c2]=[C11-m2s12C12-m2s22C13-m2s32C14-m2s42].
Let
(20)A=[1s11s21s31s4],B2=[C11-m2s12C12-m2s22C13-m2s32C14-m2s42],V2=[k2c2].
Then, (19) can be represented as
(21)A*V2=B2,
and, using the pseudoinverse, we can obtain the following:
(22)V2=(AT*A)-1*B2.
The spring stiffness k1 and damping coefficient c1 are calculated from (25) as follows:
(23)[1s11s21s31s4][k1c1]=[R11-(m1s12+c2s1+k2)R12-(m1s22+c2s2+k2)R13-(m1s32+c2s3+k2)R14-(m1s42+c2s4+k2)].Remark 1.
From (6), it was shown that R1k=R2k*S1k/S2k. Let
(24)V1=[k1c1],B1=[R11-(m1s12+c2s1+k2)R12-(m1s22+c2s2+k2)R13-(m1s32+c2s3+k2)R14-(m1s42+c2s4+k2)].
Then, we have that
(25)V1=(AT*A)-1*B1.
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 forceu1 acting on the front mass m1, (1) can be rewritten as
(26)x¨1=-k1+k2m1x1+k2m1x2-c1+c2m1x˙1+c2m1x˙2+1m1u1,x¨2=k2m2x1-k2m2x2+c2m2x˙1-c2m2x˙2.The dynamic equation (26) can be rewritten in a matrix compact form as
(27)M[x¨1x¨2]+L[x˙1x˙2]+K[x1x2]=[u10]
with
(28)M=[m100m2],L=[c1+c2-c2-c2c2],K=[k1+k2-k2-k2k2].
In general, the equation of motion forNdegrees of freedom is expressed in a matrix form as
(29)Mx¨+Lx˙+Kx=Bcu,
or, equivalently,
(30)x¨=-M-1Lx˙-M-1Kx+M-1Bcu,
where M∈ℝN×N, L∈ℝN×N, and K∈ℝN×N are the mass, damping, and stiffness matrices, respectively, while Bc∈ℝN×r is the input matrix.x
¨, x˙, andxare vectors of generalized acceleration, velocity, and displacements, respectively, and the vectoruof dimension r×1 is the input force containing r external excitations acting on the systems.Let us define the new state variablesx1=x1, x2=x˙1, x3=x2, and x4=x˙2.Substituting these variables into (26) and combining in state equations, we get the following:
(31)[x˙1x˙2x˙3x˙4]=[0100-k1+k2m1-c1+c2m1k2m1c2m10001k2m2c2m2-k2m2-c2m2][x1x2x3x4]+[01m100]u.
Using the original state variables and interchanging rows 2 and 3, columns 2 and 3 of (31) and we get the following:
(32)[x˙1x˙2x¨1x¨2]=[00100001-k1+k2m1k2m1-c1+c2m1c2m1k2m2-k2m2c2m2-c2m2][x1x2x˙1x˙2]+[001m10]u.The output equation or the measurement vectory(t), which may contain any combination of modal displacements, velocities, and/or accelerations, is given by
(33)y=[CpxCvx˙Cax¨],
where Cp, Cv, and Ca 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)y=Cax¨1.
From (32),
(35)y=x¨1=[-k1+k2m1k2m1-c1+c2m1c2m1][x1x2x˙1x˙2]+[1m1]u.Therefore, the continuous-time state-space model of the dynamic system is written as(36)x˙=Acx+Bcu,y=Ccx+Dcu
with
(37)Ac=[00100001-k1+k2m1k2m1-c1+c2m1c2m1k2m2-k2m2c2m2-c2m2],Bc=[001m10],Cc=[-k1+k2m1k2m1-c1+c2m1c2m1],Dc=[1m1],
where Ac is the state matrix, Bc is the input matrix or the state influence matrix, Cc is the output matrix or the measurement influence matrix, and Dc is the feedforward matrix or the direct transmission matrix.Once theAc, Bc, Cc, and Dc 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:(38)x(k+1)=A^x(k)+B^u(k),y(k)=C^x(k)+D^u(k),
where k is the integer discrete-time index at time instant t=kΔt, x(k) is the state vector at the discrete-time k, u(k) is the force vector, y(k) is the output vector, A^ is the discrete state system matrix, and B^ is the discrete input influence matrix for the state vector x(k). The output matrix C^=Cc and the direct transition matrix D^=Dc 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 A^ and the influence matrix B^ of the discrete-time model are related to the matrices Ac,Bc of the continuous-time model by the following expression [14]:
(39)A^=eAcΔt,B^=∫0Δt‍eAcτdτ·Bc.
The continuous-time model is calculated from the discrete-time model by
(40)Ac=ln(A^)Δt,Bc=Ac(A^-I)-1*B^,
where Δt 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 Figure4.Figure 4
Flowchart for the ERA by Juang [14].If the excitations of the dynamic system are measured by them input quantities in the vector u, 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
(41)y(k)={D^,k=0,C^A^k-1B^,k≥1.The discrete-time Markov parameters can be defined in the same way as (41). The term of C^A^k-1B^ is called the Markov parameter of the system. By using these parameters, one can define the impulse response of the system as follows:
(42)Y[k]=C^A^k-1B^,(43)Y[0]=D^,Y[1]=C^B^,Y[2]=C^A^B^,…,Y[k]=C^A^k-1B^.Consequently, the identification problem is the following: given the values ofY[k]′s, construct the constant matrices to identify the system.The algorithm begins by constructing anr×s generalized Hankel matrix.Given a number of input and output measurementsuk and yk generated by a system of unknown parameters, it is requested to identify the order of the system as well as the discrete state matrices (A^,B^,C^,D^) of the system. Then, the identified continuous state matrices (Ac,Bc,Cc,Dc) 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 matrixA^ of order n has a complete set of linearly independent eigenvectors {Ψ1,Ψ2,…,Ψn} with corresponding eigenvalues{λ1,λ2,…,λn}; then
(44)A^Ψ=ΨΛ,
where Λ is a diagonal matrix of the eigenvalues and Ψ is the matrix of the eigenvectors. The realization {A^,B^,C^} can be transformed into the realization {Λ,Ψ-1B^,C^Ψ} by using the eigenvalues and eigenvectors matrices. The diagonal matrix Λ contains the information of modal damping rates and damped natural frequencies. The matrix Ψ-1B^ 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{Λ,Ψ-1B^,C^Ψ}.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 matrixH(k-1):
(45)H(k-1)=[YkY(k+1)⋯Y(k+β+1)Yk+1Y(k+2)⋯Y(k+β)⋯⋯⋯⋯Y(k+α-1)Y(k+α)⋯Y(k+α+β-2)].
When k=1, we get the following:
(46)H(0)=[Y1Y2⋯YβY2Y3⋯Y(β+1)⋯⋯⋯⋯YαYα+1⋯Yα+β-1].In order to compute a minimum-order realization of the system{A^,B^,C^}, it is necessary to construct a shifted Hankel matrix H(1) as follows:
(47)H(1)=[Y2Y3⋯Yβ+1Y3Y4⋯Yβ+2⋯⋯⋯⋯Y(α+1)Y(α+2)⋯Y(α+β)].Substituting the Markov parameters from (43) into (45) and decomposing H(k-1) into three matrices yield the following:
(48)H(k-1)=ObA^k-1Ct,
where Ob and Ct are given as
(49)Ob=[C^C^A^C^A^2⋮C^A^α-1],Ct=[B^A^B^A^2B^⋯A^β-1B^].The block matrixOb is the observability matrix, whereas the block matrix Ct is the controllability matrix.Denote column submatrices ofB^ by B^i and row submatrices of C^ by C^j. The ERA data block matrix can be expressed by
(50)H(k-1)=[Ysi+k+tj],Ysi+k+tj=CjA^si+k-1+tjBi,(51)H(k)=ObA^kCt.
Assume that there exists a matrix H† satisfying the relation
(52)ObH†Ct=In,
where In is an identity matrix of ordern. The matrix H† which is the pseudoinverse of H(0) plays a major role in deriving the ERA. It is observed that
(53)H(0)H†H(0)=ObCtH†ObCt=ObCt=H(0).The ERA process starts with the factorization of the block data matrix (46) using the singular value decomposition
(54)H(0)=R∑‍ST,
where the  columns of matrices R and S are orthonormal and ∑‍ is a rectangular matrix given as(55)∑‍=[∑n000]
with ∑n‍=diag{σ1,σ2,…,σi,σi+1,…σn}.LetRn and ∑n‍Sn be matrices formed by the firstncolumns of R and S, respectively. Hence, the matricesH(0) and H† become
(56)H(0)=Rn∑n‍SnT,(57)H†=Sn∑n-1‍RnT,
where RnTRn=In=SnTSn.Examining the singular value∑n‍ of the Hankel matrix H(0), it is possible to determine the order of the system.Comparing (56) and (51) with k=0, we get that Ob=Rn∑n1/2‍ and Ct=∑n1/2‍SnT.From (51), the first r columns of the observability matrix Ob form the input matrix B,^ whereas the first m rows of the controllability matrix Ct form the output matrix C^.Withk=1 in (51), we get that
(58)H(1)=ObA^Ct=Rn∑n1/2‍A^∑n1/2‍SnT.One obvious solution for the state matrixA^ becomes
(59)A^=∑n-1/2‍RnTH(1)Sn∑n-1/2‍.
Let Oi be a null matrix of order i, let Ii be an identity matrix of order i, and let the matrices EmT and ErT be defined as
(60)EmT=[Im0m⋯0m],ErT=[Ir0r⋯0r],
where m is the number of outputs and r 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:
(61)Yk=EmTRn∑n-1/2‍[∑n-1/2‍RnTH(1)Sn∑n-1/2‍]k-1∑n-1/2‍SnTEr.
Recall that in (42) Y[k]=C^A^k-1B^.Hence, a minimal realization is given as follows:(62)A^=∑n-1/2RnTH(1)Sn∑n-1/2,B^=∑n-1/2SnTEr,C^=EmTRn∑n-1/2,D^=Y(0).
The continuous state matrix Ac and the input influence matrix Bc are obtained from (40), and the physical parameters-matrices M, L, and K are embedded in the state matrix Ac.


## 3.1. Estimation of Model Parameters by Curve Fitting

When (16) is curve fitted into displacement experimental data, the constants Cik and sk (i=1,2, and k=1,…,4) can be easily found resulting in the following system of equations that can be solved for k2 and c2:
(17)C1k=m2sk2+c2sk+k2,(18)C2k=c2sk+k2.Fork=1,…,4, (17) can be written in a matrix form as
(19)[1s11s21s31s4][k2c2]=[C11-m2s12C12-m2s22C13-m2s32C14-m2s42].
Let
(20)A=[1s11s21s31s4],B2=[C11-m2s12C12-m2s22C13-m2s32C14-m2s42],V2=[k2c2].
Then, (19) can be represented as
(21)A*V2=B2,
and, using the pseudoinverse, we can obtain the following:
(22)V2=(AT*A)-1*B2.
The spring stiffness k1 and damping coefficient c1 are calculated from (25) as follows:
(23)[1s11s21s31s4][k1c1]=[R11-(m1s12+c2s1+k2)R12-(m1s22+c2s2+k2)R13-(m1s32+c2s3+k2)R14-(m1s42+c2s4+k2)].Remark 1.
From (6), it was shown that R1k=R2k*S1k/S2k. Let
(24)V1=[k1c1],B1=[R11-(m1s12+c2s1+k2)R12-(m1s22+c2s2+k2)R13-(m1s32+c2s3+k2)R14-(m1s42+c2s4+k2)].
Then, we have that
(25)V1=(AT*A)-1*B1.
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 forceu1 acting on the front mass m1, (1) can be rewritten as
(26)x¨1=-k1+k2m1x1+k2m1x2-c1+c2m1x˙1+c2m1x˙2+1m1u1,x¨2=k2m2x1-k2m2x2+c2m2x˙1-c2m2x˙2.The dynamic equation (26) can be rewritten in a matrix compact form as
(27)M[x¨1x¨2]+L[x˙1x˙2]+K[x1x2]=[u10]
with
(28)M=[m100m2],L=[c1+c2-c2-c2c2],K=[k1+k2-k2-k2k2].
In general, the equation of motion forNdegrees of freedom is expressed in a matrix form as
(29)Mx¨+Lx˙+Kx=Bcu,
or, equivalently,
(30)x¨=-M-1Lx˙-M-1Kx+M-1Bcu,
where M∈ℝN×N, L∈ℝN×N, and K∈ℝN×N are the mass, damping, and stiffness matrices, respectively, while Bc∈ℝN×r is the input matrix.x
¨, x˙, andxare vectors of generalized acceleration, velocity, and displacements, respectively, and the vectoruof dimension r×1 is the input force containing r external excitations acting on the systems.Let us define the new state variablesx1=x1, x2=x˙1, x3=x2, and x4=x˙2.Substituting these variables into (26) and combining in state equations, we get the following:
(31)[x˙1x˙2x˙3x˙4]=[0100-k1+k2m1-c1+c2m1k2m1c2m10001k2m2c2m2-k2m2-c2m2][x1x2x3x4]+[01m100]u.
Using the original state variables and interchanging rows 2 and 3, columns 2 and 3 of (31) and we get the following:
(32)[x˙1x˙2x¨1x¨2]=[00100001-k1+k2m1k2m1-c1+c2m1c2m1k2m2-k2m2c2m2-c2m2][x1x2x˙1x˙2]+[001m10]u.The output equation or the measurement vectory(t), which may contain any combination of modal displacements, velocities, and/or accelerations, is given by
(33)y=[CpxCvx˙Cax¨],
where Cp, Cv, and Ca 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)y=Cax¨1.
From (32),
(35)y=x¨1=[-k1+k2m1k2m1-c1+c2m1c2m1][x1x2x˙1x˙2]+[1m1]u.Therefore, the continuous-time state-space model of the dynamic system is written as(36)x˙=Acx+Bcu,y=Ccx+Dcu
with
(37)Ac=[00100001-k1+k2m1k2m1-c1+c2m1c2m1k2m2-k2m2c2m2-c2m2],Bc=[001m10],Cc=[-k1+k2m1k2m1-c1+c2m1c2m1],Dc=[1m1],
where Ac is the state matrix, Bc is the input matrix or the state influence matrix, Cc is the output matrix or the measurement influence matrix, and Dc is the feedforward matrix or the direct transmission matrix.Once theAc, Bc, Cc, and Dc 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:(38)x(k+1)=A^x(k)+B^u(k),y(k)=C^x(k)+D^u(k),
where k is the integer discrete-time index at time instant t=kΔt, x(k) is the state vector at the discrete-time k, u(k) is the force vector, y(k) is the output vector, A^ is the discrete state system matrix, and B^ is the discrete input influence matrix for the state vector x(k). The output matrix C^=Cc and the direct transition matrix D^=Dc 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 A^ and the influence matrix B^ of the discrete-time model are related to the matrices Ac,Bc of the continuous-time model by the following expression [14]:
(39)A^=eAcΔt,B^=∫0Δt‍eAcτdτ·Bc.
The continuous-time model is calculated from the discrete-time model by
(40)Ac=ln(A^)Δt,Bc=Ac(A^-I)-1*B^,
where Δt 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 Figure4.Figure 4
Flowchart for the ERA by Juang [14].If the excitations of the dynamic system are measured by them input quantities in the vector u, 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
(41)y(k)={D^,k=0,C^A^k-1B^,k≥1.The discrete-time Markov parameters can be defined in the same way as (41). The term of C^A^k-1B^ is called the Markov parameter of the system. By using these parameters, one can define the impulse response of the system as follows:
(42)Y[k]=C^A^k-1B^,(43)Y[0]=D^,Y[1]=C^B^,Y[2]=C^A^B^,…,Y[k]=C^A^k-1B^.Consequently, the identification problem is the following: given the values ofY[k]′s, construct the constant matrices to identify the system.The algorithm begins by constructing anr×s generalized Hankel matrix.Given a number of input and output measurementsuk and yk generated by a system of unknown parameters, it is requested to identify the order of the system as well as the discrete state matrices (A^,B^,C^,D^) of the system. Then, the identified continuous state matrices (Ac,Bc,Cc,Dc) 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 matrixA^ of order n has a complete set of linearly independent eigenvectors {Ψ1,Ψ2,…,Ψn} with corresponding eigenvalues{λ1,λ2,…,λn}; then
(44)A^Ψ=ΨΛ,
where Λ is a diagonal matrix of the eigenvalues and Ψ is the matrix of the eigenvectors. The realization {A^,B^,C^} can be transformed into the realization {Λ,Ψ-1B^,C^Ψ} by using the eigenvalues and eigenvectors matrices. The diagonal matrix Λ contains the information of modal damping rates and damped natural frequencies. The matrix Ψ-1B^ 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{Λ,Ψ-1B^,C^Ψ}.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 matrixH(k-1):
(45)H(k-1)=[YkY(k+1)⋯Y(k+β+1)Yk+1Y(k+2)⋯Y(k+β)⋯⋯⋯⋯Y(k+α-1)Y(k+α)⋯Y(k+α+β-2)].
When k=1, we get the following:
(46)H(0)=[Y1Y2⋯YβY2Y3⋯Y(β+1)⋯⋯⋯⋯YαYα+1⋯Yα+β-1].In order to compute a minimum-order realization of the system{A^,B^,C^}, it is necessary to construct a shifted Hankel matrix H(1) as follows:
(47)H(1)=[Y2Y3⋯Yβ+1Y3Y4⋯Yβ+2⋯⋯⋯⋯Y(α+1)Y(α+2)⋯Y(α+β)].Substituting the Markov parameters from (43) into (45) and decomposing H(k-1) into three matrices yield the following:
(48)H(k-1)=ObA^k-1Ct,
where Ob and Ct are given as
(49)Ob=[C^C^A^C^A^2⋮C^A^α-1],Ct=[B^A^B^A^2B^⋯A^β-1B^].The block matrixOb is the observability matrix, whereas the block matrix Ct is the controllability matrix.Denote column submatrices ofB^ by B^i and row submatrices of C^ by C^j. The ERA data block matrix can be expressed by
(50)H(k-1)=[Ysi+k+tj],Ysi+k+tj=CjA^si+k-1+tjBi,(51)H(k)=ObA^kCt.
Assume that there exists a matrix H† satisfying the relation
(52)ObH†Ct=In,
where In is an identity matrix of ordern. The matrix H† which is the pseudoinverse of H(0) plays a major role in deriving the ERA. It is observed that
(53)H(0)H†H(0)=ObCtH†ObCt=ObCt=H(0).The ERA process starts with the factorization of the block data matrix (46) using the singular value decomposition
(54)H(0)=R∑‍ST,
where the  columns of matrices R and S are orthonormal and ∑‍ is a rectangular matrix given as(55)∑‍=[∑n000]
with ∑n‍=diag{σ1,σ2,…,σi,σi+1,…σn}.LetRn and ∑n‍Sn be matrices formed by the firstncolumns of R and S, respectively. Hence, the matricesH(0) and H† become
(56)H(0)=Rn∑n‍SnT,(57)H†=Sn∑n-1‍RnT,
where RnTRn=In=SnTSn.Examining the singular value∑n‍ of the Hankel matrix H(0), it is possible to determine the order of the system.Comparing (56) and (51) with k=0, we get that Ob=Rn∑n1/2‍ and Ct=∑n1/2‍SnT.From (51), the first r columns of the observability matrix Ob form the input matrix B,^ whereas the first m rows of the controllability matrix Ct form the output matrix C^.Withk=1 in (51), we get that
(58)H(1)=ObA^Ct=Rn∑n1/2‍A^∑n1/2‍SnT.One obvious solution for the state matrixA^ becomes
(59)A^=∑n-1/2‍RnTH(1)Sn∑n-1/2‍.
Let Oi be a null matrix of order i, let Ii be an identity matrix of order i, and let the matrices EmT and ErT be defined as
(60)EmT=[Im0m⋯0m],ErT=[Ir0r⋯0r],
where m is the number of outputs and r 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:
(61)Yk=EmTRn∑n-1/2‍[∑n-1/2‍RnTH(1)Sn∑n-1/2‍]k-1∑n-1/2‍SnTEr.
Recall that in (42) Y[k]=C^A^k-1B^.Hence, a minimal realization is given as follows:(62)A^=∑n-1/2RnTH(1)Sn∑n-1/2,B^=∑n-1/2SnTEr,C^=EmTRn∑n-1/2,D^=Y(0).
The continuous state matrix Ac and the input influence matrix Bc are obtained from (40), and the physical parameters-matrices M, L, and K are embedded in the state matrix Ac.


## 3.2.1. Formulation of the SSR from the Model Dynamic Equations

Considering an input forceu1 acting on the front mass m1, (1) can be rewritten as
(26)x¨1=-k1+k2m1x1+k2m1x2-c1+c2m1x˙1+c2m1x˙2+1m1u1,x¨2=k2m2x1-k2m2x2+c2m2x˙1-c2m2x˙2.The dynamic equation (26) can be rewritten in a matrix compact form as
(27)M[x¨1x¨2]+L[x˙1x˙2]+K[x1x2]=[u10]
with
(28)M=[m100m2],L=[c1+c2-c2-c2c2],K=[k1+k2-k2-k2k2].
In general, the equation of motion forNdegrees of freedom is expressed in a matrix form as
(29)Mx¨+Lx˙+Kx=Bcu,
or, equivalently,
(30)x¨=-M-1Lx˙-M-1Kx+M-1Bcu,
where M∈ℝN×N, L∈ℝN×N, and K∈ℝN×N are the mass, damping, and stiffness matrices, respectively, while Bc∈ℝN×r is the input matrix.x
¨, x˙, andxare vectors of generalized acceleration, velocity, and displacements, respectively, and the vectoruof dimension r×1 is the input force containing r external excitations acting on the systems.Let us define the new state variablesx1=x1, x2=x˙1, x3=x2, and x4=x˙2.Substituting these variables into (26) and combining in state equations, we get the following:
(31)[x˙1x˙2x˙3x˙4]=[0100-k1+k2m1-c1+c2m1k2m1c2m10001k2m2c2m2-k2m2-c2m2][x1x2x3x4]+[01m100]u.
Using the original state variables and interchanging rows 2 and 3, columns 2 and 3 of (31) and we get the following:
(32)[x˙1x˙2x¨1x¨2]=[00100001-k1+k2m1k2m1-c1+c2m1c2m1k2m2-k2m2c2m2-c2m2][x1x2x˙1x˙2]+[001m10]u.The output equation or the measurement vectory(t), which may contain any combination of modal displacements, velocities, and/or accelerations, is given by
(33)y=[CpxCvx˙Cax¨],
where Cp, Cv, and Ca 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)y=Cax¨1.
From (32),
(35)y=x¨1=[-k1+k2m1k2m1-c1+c2m1c2m1][x1x2x˙1x˙2]+[1m1]u.Therefore, the continuous-time state-space model of the dynamic system is written as(36)x˙=Acx+Bcu,y=Ccx+Dcu
with
(37)Ac=[00100001-k1+k2m1k2m1-c1+c2m1c2m1k2m2-k2m2c2m2-c2m2],Bc=[001m10],Cc=[-k1+k2m1k2m1-c1+c2m1c2m1],Dc=[1m1],
where Ac is the state matrix, Bc is the input matrix or the state influence matrix, Cc is the output matrix or the measurement influence matrix, and Dc is the feedforward matrix or the direct transmission matrix.Once theAc, Bc, Cc, and Dc 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:(38)x(k+1)=A^x(k)+B^u(k),y(k)=C^x(k)+D^u(k),
where k is the integer discrete-time index at time instant t=kΔt, x(k) is the state vector at the discrete-time k, u(k) is the force vector, y(k) is the output vector, A^ is the discrete state system matrix, and B^ is the discrete input influence matrix for the state vector x(k). The output matrix C^=Cc and the direct transition matrix D^=Dc 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 A^ and the influence matrix B^ of the discrete-time model are related to the matrices Ac,Bc of the continuous-time model by the following expression [14]:
(39)A^=eAcΔt,B^=∫0Δt‍eAcτdτ·Bc.
The continuous-time model is calculated from the discrete-time model by
(40)Ac=ln(A^)Δt,Bc=Ac(A^-I)-1*B^,
where Δt 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 Figure4.Figure 4
Flowchart for the ERA by Juang [14].If the excitations of the dynamic system are measured by them input quantities in the vector u, 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
(41)y(k)={D^,k=0,C^A^k-1B^,k≥1.The discrete-time Markov parameters can be defined in the same way as (41). The term of C^A^k-1B^ is called the Markov parameter of the system. By using these parameters, one can define the impulse response of the system as follows:
(42)Y[k]=C^A^k-1B^,(43)Y[0]=D^,Y[1]=C^B^,Y[2]=C^A^B^,…,Y[k]=C^A^k-1B^.Consequently, the identification problem is the following: given the values ofY[k]′s, construct the constant matrices to identify the system.The algorithm begins by constructing anr×s generalized Hankel matrix.Given a number of input and output measurementsuk and yk generated by a system of unknown parameters, it is requested to identify the order of the system as well as the discrete state matrices (A^,B^,C^,D^) of the system. Then, the identified continuous state matrices (Ac,Bc,Cc,Dc) 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 matrixA^ of order n has a complete set of linearly independent eigenvectors {Ψ1,Ψ2,…,Ψn} with corresponding eigenvalues{λ1,λ2,…,λn}; then
(44)A^Ψ=ΨΛ,
where Λ is a diagonal matrix of the eigenvalues and Ψ is the matrix of the eigenvectors. The realization {A^,B^,C^} can be transformed into the realization {Λ,Ψ-1B^,C^Ψ} by using the eigenvalues and eigenvectors matrices. The diagonal matrix Λ contains the information of modal damping rates and damped natural frequencies. The matrix Ψ-1B^ 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{Λ,Ψ-1B^,C^Ψ}.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 matrixH(k-1):
(45)H(k-1)=[YkY(k+1)⋯Y(k+β+1)Yk+1Y(k+2)⋯Y(k+β)⋯⋯⋯⋯Y(k+α-1)Y(k+α)⋯Y(k+α+β-2)].
When k=1, we get the following:
(46)H(0)=[Y1Y2⋯YβY2Y3⋯Y(β+1)⋯⋯⋯⋯YαYα+1⋯Yα+β-1].In order to compute a minimum-order realization of the system{A^,B^,C^}, it is necessary to construct a shifted Hankel matrix H(1) as follows:
(47)H(1)=[Y2Y3⋯Yβ+1Y3Y4⋯Yβ+2⋯⋯⋯⋯Y(α+1)Y(α+2)⋯Y(α+β)].Substituting the Markov parameters from (43) into (45) and decomposing H(k-1) into three matrices yield the following:
(48)H(k-1)=ObA^k-1Ct,
where Ob and Ct are given as
(49)Ob=[C^C^A^C^A^2⋮C^A^α-1],Ct=[B^A^B^A^2B^⋯A^β-1B^].The block matrixOb is the observability matrix, whereas the block matrix Ct is the controllability matrix.Denote column submatrices ofB^ by B^i and row submatrices of C^ by C^j. The ERA data block matrix can be expressed by
(50)H(k-1)=[Ysi+k+tj],Ysi+k+tj=CjA^si+k-1+tjBi,(51)H(k)=ObA^kCt.
Assume that there exists a matrix H† satisfying the relation
(52)ObH†Ct=In,
where In is an identity matrix of ordern. The matrix H† which is the pseudoinverse of H(0) plays a major role in deriving the ERA. It is observed that
(53)H(0)H†H(0)=ObCtH†ObCt=ObCt=H(0).The ERA process starts with the factorization of the block data matrix (46) using the singular value decomposition
(54)H(0)=R∑‍ST,
where the  columns of matrices R and S are orthonormal and ∑‍ is a rectangular matrix given as(55)∑‍=[∑n000]
with ∑n‍=diag{σ1,σ2,…,σi,σi+1,…σn}.LetRn and ∑n‍Sn be matrices formed by the firstncolumns of R and S, respectively. Hence, the matricesH(0) and H† become
(56)H(0)=Rn∑n‍SnT,(57)H†=Sn∑n-1‍RnT,
where RnTRn=In=SnTSn.Examining the singular value∑n‍ of the Hankel matrix H(0), it is possible to determine the order of the system.Comparing (56) and (51) with k=0, we get that Ob=Rn∑n1/2‍ and Ct=∑n1/2‍SnT.From (51), the first r columns of the observability matrix Ob form the input matrix B,^ whereas the first m rows of the controllability matrix Ct form the output matrix C^.Withk=1 in (51), we get that
(58)H(1)=ObA^Ct=Rn∑n1/2‍A^∑n1/2‍SnT.One obvious solution for the state matrixA^ becomes
(59)A^=∑n-1/2‍RnTH(1)Sn∑n-1/2‍.
Let Oi be a null matrix of order i, let Ii be an identity matrix of order i, and let the matrices EmT and ErT be defined as
(60)EmT=[Im0m⋯0m],ErT=[Ir0r⋯0r],
where m is the number of outputs and r 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:
(61)Yk=EmTRn∑n-1/2‍[∑n-1/2‍RnTH(1)Sn∑n-1/2‍]k-1∑n-1/2‍SnTEr.
Recall that in (42) Y[k]=C^A^k-1B^.Hence, a minimal realization is given as follows:(62)A^=∑n-1/2RnTH(1)Sn∑n-1/2,B^=∑n-1/2SnTEr,C^=EmTRn∑n-1/2,D^=Y(0).
The continuous state matrix Ac and the input influence matrix Bc are obtained from (40), and the physical parameters-matrices M, L, and K are embedded in the state matrix Ac.


## 4. Results and Discussion

### 4.1. Parameters Estimation from Curve-Fitting Approach

From the curve fitting, the values forsi and C1j (j=1,…,4) are found to be the following:C11=-253.6,C12=83.75,C13=333,C14=-163.2,s1=-7.7540, s2=-6.238, s3=-9.603,s4=-10.87.The result from the curve fitting is shown in Figure5.Figure 5
Result from curve fitting.


### 4.2. Vehicle Crash Experimental Data Analysis

It is observed from Figure6 that the dynamic crash from the real vehicle crash test is 53.17 cm and occurs at time tc=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 tc=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. Letm1 be the mass of the chassis, m2 the mass of passenger compartment, and mt=873kg the total mass of the vehicle.Case 1 (m1<m2).
One has thatm1=(1/3)mt, and  m2=(2/3)mt. From Figure 8, the dynamic crash of m2 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 m1 is 42.5 cm and occurs after 0.17 s.Figure 8
Comparative analysis between vehicle crash test and model results form1=(1/3)mt.Case 2 (m1>m2).
One has thatm1=(2/3)mt, and m2=(1/3)mt. From Figure 9, the dynamic crash of the passenger compartment m2 is 66.2 cm. The time for dynamic crash increases further up to 0.15 s. The dynamic crash of m1 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 form2=(1/3)mt.Case 3 (m1<m2).
One has thatm1=(1/4)mt, and m2=(3/4)mt. From Figure 10, the dynamic crash of the passenger compartment m2 is 69 cm. The time for dynamic crash is 0.15 s. The dynamic crash of m1 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 form1=(1/4)mt.Case 4 (m1>m2).
One has thatm1=(3/4)mt, and m2=(1/4)mt. From Figure 11, the dynamic crash of the passenger compartment m2 is 49.8 cm, and the time for dynamic crash is 0.11 s. The dynamic crash of m1 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 form2=(1/4)mt.A summary of main results is shown in Table1. The values for k2 and c2 are the first and second entries of vector V2 in (22). The values for k2 and c2 depend on the value of mass m2 taken into consideration. Values for k1 and c1 are obtained from vector V1 in (25). The stiffness coefficients which result in a closer vehicle crash reconstruction are found to be k1=74681N/m, and k2=45821N/m, and the damping coefficients are: c1=18176Ns/m and c2=11196Ns/m when the mass of the chassis is 3/4 the total mass of the vehicle (mt), 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.
Parameters
Case1
Case2
Case3
Case4
m
1
=
(
1
/
3
)
m
t
m
1
=
(
2
/
3
)
m
t
m
1
=
(
1
/
4
)
m
t
m
1
=
(
3
/
4
)
m
t
Dynamic crash (cm) = 80
Dynamic crash (cm) = 66
Dynamic crash (cm) = 69
Dynamic crash (cm) = 49.8
Time of crash (s) = 0.17
Time of crash (s) = 0.15
Time of crash (s) = 0.14
Time of crash (s) = 0.11
k
1 [N/m]
45929
74681
k
2 [N/m]
40731
45821
c
1 [Ns/m]
13687
18176
c
2 [Ns/m]
9952
11196It is observed that the passenger compartmentm2 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 m1=582kg (i.e., 3/4 of the total mass of the vehicle) and m2=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.,n=2) for a single degree of freedom and N=100, the number of samples to assemble the Hankel matrix. It is observed that the continuous system matrices from the ERA are the following:
(63)Ac=[95.5776-545.1246545.12-131.6485],Bc=[-1474.85798.1],C=[0.16390.5717],D=0.1964.
From the state matrixAc, the eigenvectors Ψ and eigenvalues Λ were found to be
(64)Ψ,Ψ*=[0.1474±0.6916i0.7071],
and Λ=diag{-18.04±533.15i}.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 Figure12.The Bode diagram forn=2, N=100, and N=50.
(a)
(b)Considering a 4th-order system (i.e,n=4) for a two-degree-of-freedom system and the number of samples to assemble the Hankel matrix equal to N=60, it is observed that the continuous state-space realizations from the ERA are found to be the following:
(65)Ac=[-22.5-526.356.2141.5526.3-243.8631.5279.956.2-631.5-171.2-1363.2-141.5279.9136.2-445.9],Bc=[-1990.85436.83302.4-4233.2],C=[0.21490.5366-0.3368-0.3826],D=0.1964
with
(66)Ψ1,Ψ1*=[0.68640.2979∓0.5573i-0.1126±0.1269i-0.0715±0.3092i],Ψ2,Ψ2*=[0.0537∓0.0828i-0.1269∓0.3363i0.70020.0587∓0.6061i],
and Λ=diag{-274.9±501.4i,-166.8±1476.7i}.The natural frequencies and damping ratios are given asωn1= 3592.8  rad/s, ωn2 = 9337.4 rad/s and ζ1 = 0.481, ζ2 = 0.112.One can see the stability of the model from pole-zero map as shown in Figure13 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 forn=2 and N=100.Figure 14
Pole-zero map forn=2 and N=50.Figure 15
Responses of a single-degree-of-Freedom model to excitation force.For higher sample size (i.e.,N=2501 with order n=4), the system matrices are given as follows:
(67)Ac=[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],Bc=[1280.1-5366.11636.4-946.8],C=[-0.144-0.5322-0.1656-0.0897],D=0.1964
with
(68)Ψ1,Ψ1*=[0.0053±0.6800i0.6840-0.0869∓0.1733i-0.1561±0.0880i],Ψ2,Ψ2*=[0.0878±0.1732i0.1572∓0.0850i-0.0027±0.6799i0.6841],
and Λ=diag{-4.39±617i,-3.85±575i}.The natural frequencies and damping ratios are given asωn1 = 617 rad/s, ωn2 = 575 rad/s and ζ1 = 0.00712, ζ2 = 0.006680.Remark 3.
Based on the state-space realization obtained by ERA, we should be able to transfer the state matrixAc to the original obtained from (36). Hence, extract the mass, stiffness, and damping matrices.


## 4.1. Parameters Estimation from Curve-Fitting Approach

From the curve fitting, the values forsi and C1j (j=1,…,4) are found to be the following:C11=-253.6,C12=83.75,C13=333,C14=-163.2,s1=-7.7540, s2=-6.238, s3=-9.603,s4=-10.87.The result from the curve fitting is shown in Figure5.Figure 5
Result from curve fitting.


## 4.2. Vehicle Crash Experimental Data Analysis

It is observed from Figure6 that the dynamic crash from the real vehicle crash test is 53.17 cm and occurs at time tc=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 tc=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. Letm1 be the mass of the chassis, m2 the mass of passenger compartment, and mt=873kg the total mass of the vehicle.Case 1 (m1<m2).
One has thatm1=(1/3)mt, and  m2=(2/3)mt. From Figure 8, the dynamic crash of m2 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 m1 is 42.5 cm and occurs after 0.17 s.Figure 8
Comparative analysis between vehicle crash test and model results form1=(1/3)mt.Case 2 (m1>m2).
One has thatm1=(2/3)mt, and m2=(1/3)mt. From Figure 9, the dynamic crash of the passenger compartment m2 is 66.2 cm. The time for dynamic crash increases further up to 0.15 s. The dynamic crash of m1 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 form2=(1/3)mt.Case 3 (m1<m2).
One has thatm1=(1/4)mt, and m2=(3/4)mt. From Figure 10, the dynamic crash of the passenger compartment m2 is 69 cm. The time for dynamic crash is 0.15 s. The dynamic crash of m1 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 form1=(1/4)mt.Case 4 (m1>m2).
One has thatm1=(3/4)mt, and m2=(1/4)mt. From Figure 11, the dynamic crash of the passenger compartment m2 is 49.8 cm, and the time for dynamic crash is 0.11 s. The dynamic crash of m1 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 form2=(1/4)mt.A summary of main results is shown in Table1. The values for k2 and c2 are the first and second entries of vector V2 in (22). The values for k2 and c2 depend on the value of mass m2 taken into consideration. Values for k1 and c1 are obtained from vector V1 in (25). The stiffness coefficients which result in a closer vehicle crash reconstruction are found to be k1=74681N/m, and k2=45821N/m, and the damping coefficients are: c1=18176Ns/m and c2=11196Ns/m when the mass of the chassis is 3/4 the total mass of the vehicle (mt), 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.
Parameters
Case1
Case2
Case3
Case4
m
1
=
(
1
/
3
)
m
t
m
1
=
(
2
/
3
)
m
t
m
1
=
(
1
/
4
)
m
t
m
1
=
(
3
/
4
)
m
t
Dynamic crash (cm) = 80
Dynamic crash (cm) = 66
Dynamic crash (cm) = 69
Dynamic crash (cm) = 49.8
Time of crash (s) = 0.17
Time of crash (s) = 0.15
Time of crash (s) = 0.14
Time of crash (s) = 0.11
k
1 [N/m]
45929
74681
k
2 [N/m]
40731
45821
c
1 [Ns/m]
13687
18176
c
2 [Ns/m]
9952
11196It is observed that the passenger compartmentm2 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 m1=582kg (i.e., 3/4 of the total mass of the vehicle) and m2=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.,n=2) for a single degree of freedom and N=100, the number of samples to assemble the Hankel matrix. It is observed that the continuous system matrices from the ERA are the following:
(63)Ac=[95.5776-545.1246545.12-131.6485],Bc=[-1474.85798.1],C=[0.16390.5717],D=0.1964.
From the state matrixAc, the eigenvectors Ψ and eigenvalues Λ were found to be
(64)Ψ,Ψ*=[0.1474±0.6916i0.7071],
and Λ=diag{-18.04±533.15i}.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 Figure12.The Bode diagram forn=2, N=100, and N=50.
(a)
(b)Considering a 4th-order system (i.e,n=4) for a two-degree-of-freedom system and the number of samples to assemble the Hankel matrix equal to N=60, it is observed that the continuous state-space realizations from the ERA are found to be the following:
(65)Ac=[-22.5-526.356.2141.5526.3-243.8631.5279.956.2-631.5-171.2-1363.2-141.5279.9136.2-445.9],Bc=[-1990.85436.83302.4-4233.2],C=[0.21490.5366-0.3368-0.3826],D=0.1964
with
(66)Ψ1,Ψ1*=[0.68640.2979∓0.5573i-0.1126±0.1269i-0.0715±0.3092i],Ψ2,Ψ2*=[0.0537∓0.0828i-0.1269∓0.3363i0.70020.0587∓0.6061i],
and Λ=diag{-274.9±501.4i,-166.8±1476.7i}.The natural frequencies and damping ratios are given asωn1= 3592.8  rad/s, ωn2 = 9337.4 rad/s and ζ1 = 0.481, ζ2 = 0.112.One can see the stability of the model from pole-zero map as shown in Figure13 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 forn=2 and N=100.Figure 14
Pole-zero map forn=2 and N=50.Figure 15
Responses of a single-degree-of-Freedom model to excitation force.For higher sample size (i.e.,N=2501 with order n=4), the system matrices are given as follows:
(67)Ac=[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],Bc=[1280.1-5366.11636.4-946.8],C=[-0.144-0.5322-0.1656-0.0897],D=0.1964
with
(68)Ψ1,Ψ1*=[0.0053±0.6800i0.6840-0.0869∓0.1733i-0.1561±0.0880i],Ψ2,Ψ2*=[0.0878±0.1732i0.1572∓0.0850i-0.0027±0.6799i0.6841],
and Λ=diag{-4.39±617i,-3.85±575i}.The natural frequencies and damping ratios are given asωn1 = 617 rad/s, ωn2 = 575 rad/s and ζ1 = 0.00712, ζ2 = 0.006680.Remark 3.
Based on the state-space realization obtained by ERA, we should be able to transfer the state matrixAc 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—tm. 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 massm1 representing the chassis of the vehicle is less than m2, 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.



---

*Source: 262196-2013-11-14.xml*