Abstract
Fractional derivative viscoelastic model is used in the analysis of a frontal impact of a vehicle against a rigid obstacle. The frontal part of the vehicle is first modeled as a viscoelastic fractional rod and then it is modeled as two different viscoelastic fractional rods with a different length. In the second model also the friction is taken into account. A motion is analyzed during several phases because of both different lengths of the rods and the presence of a dry friction force in the later model. Governing systems of differential equations together with the corresponding initial conditions are derived. Parameter identification is done on the basis of the existing experimental results using the solution of a posed impact problem. What makes the problem more complex, regarding the second model, is the fact that it belongs to the class of nonsmooth fractional order systems, which require special treatment when dealing with deformation history during different motion phases.
1. Introduction
The most of the fatal injuries from vehicle accidents happen in frontal crashes. Nowadays, various crash tests are performed in order to improve safety of vehicle occupants. Both offset and full width barrier tests are mainly conducted in the analysis of frontal impacts. Also, collision of a vehicle against a tree is not rare, so car to pole frontal impact experiments are performed for determination of deformations and velocities during the impact phase and for risk estimation. National Highway Traffic Safety Administration (NHTSA), US Department of Transportation, performed twelve staged collisions, as well as their computational simulations [1]. Parameter identification of a vehicle and occupant model in the case of a side impact crash, on the basis of experimental data from NHTSA, was done by Sun et al.; see [2]. Impact of different vehicles on a road safety barrier was studied in [3]. Rigid body collision modeling, using experimental results from NHTSA, which requires estimation of restitution coefficients in advance is presented in [4]. Mechanical models of a viscoelastic body, such as Kelvin-Voigt, Maxwell, and Zener, are widely used in the analysis of crash tests and accident reconstructions; see [5]. Experimental results of a car to pole frontal impact test are presented in the papers of Pawlus et al. [6, 7]. Numerical simulations obtained by using one viscoelastic body described by some of the standard viscoelastic models, analyzed in the two last mentioned papers, did not show satisfactory agreement with the experimental data. Improved agreement was achieved by the use of a double-spring-mass-damper system, presented in [8]. The aim of this paper is to find a model which could predict the behavior of a body during a crash test with acceptable accuracy. Due to the complexity of impact phenomena, better agreement with experimental results could be expected by using a more complex viscoelastic model, such as the Zener model of a standard linear viscoelastic body, or the fractional Zener model, which comprises fractional derivatives of stress and strain; see [9]. It is well known that the fractional Zener model is widely used to describe the behavior of real materials which exhibit viscoelastic properties. Reasons for that are both good agreement with experimental results and small number of model parameters; see [10, 11]. Mainardi and Spada analyzed different generalized viscoelastic models by means of fractional calculus in [12]. Some applications of fractional derivative viscoelastic models in the problem of impact of a body against a rigid wall are studied in [13, 14]. In this study, in order to achieve better agreement between a model prediction and a crash test experiment than the agreement presented in [6], the classical and the fractional Zener model of a viscoelastic body will be applied to model a deformable rod. Another direction for improvement of a crash test model is to reconsider the mechanism of vehicle deformation occurrence during a crash and to pose a question: “Is there a better way to model the behavior of a vehicle during a crash test than taking a complex viscoelastic model?” Namely, during the frontal impact of a car against an obstacle, different frontal parts of the vehicle start to deform in different time instants. Probably, this effect cannot be described by making a viscoelastic model more complex. In this paper, we also propose the mechanical model consisting of two parallel viscoelastic rods of different lengths, which will not start to deform at the same time. Furthermore, a dry friction effect is introduced to simulate the frontal wheels blockade. In that way, the model becomes nonlinear, firstly because of different lengths of the parallel viscoelastic rods and, secondly, due to the dry friction force modeled by the Coulomb friction law. On the basis of experimental data presented in [6], parameter identification procedure is performed firstly for the model with one viscoelastic rod and secondly for the model consisting of two rods. Viscoelastic bodies are described using the standard Zener model as well as the fractional Zener model.
2. Mechanical Model
In the present analysis, a vehicle in a crash test is modeled by a rigid block and a viscoelastic material connected to it, like in [13, 15], according to Hertz impact theory. The viscoelastic material represents frontal parts of the car which are exposed to deformation during the collision (deformable layer), while other parts are presented by the rigid body. Deformable layer is taken in the form of a rod of constant cross sectional area and its deformation simulates deformation of frontal parts of the vehicle. However, those parts are of different material, size, position, and shape, and not all of them start to deform at the same time during the impact process. In order to simulate these effects, the second model is considered as well, in which the deformable layer is modeled by two parallel viscoelastic rods of different lengths. In that way, the most frontal parts which deform first are modeled by the longer rod, while the parts which start to deform after certain time (it is measured in milliseconds) are modeled by the shorter rod. Both the standard and the fractional Zener models are used.
2.1. System with One Viscoelastic Rod
The system under consideration consists of a rigid block of mass and a viscoelastic rod attached to it, as it is shown in Figure 1, which moves translatory without friction along the horizontal surface. At certain moment, say , when the velocity of the block is , the viscoelastic rod contacts a rigid obstacle and impact occurs. This problem is thoroughly studied in [15].
Governing equations describing this problem with the initial conditions readwhere , and if derivative order is , the expressionstands for the Riemann-Liouville fractional derivative of a function ; see [16]. The Euler Gamma function is denoted by symbol . In the constitutive equation (2), derivative order is bounded , stands for the force in the viscoelastic rod, and , , and are modulus of elasticity, cross-sectional area, and length of the viscoelastic rod. Constants , , and must satisfy conditionsaccording to the second law of thermodynamics; see [9]. Equation (2) represents the fractional Zener model for , while in the case when , it represents the standard Zener model with classical integer order derivatives. On the basis of experimental data of a crash test presented in [6], the task is to determine the parameters of the fractional Zener model for the rod, which is a generalization of the standard Zener model. Introducing dimensionless quantitiesand by omitting the bar for the sake of simplicity, the system of equations readswith dimensionless restrictions in the form of (5)2,3. Derivatives in dimensionless equations denote derivatives with respect to dimensionless time. Parameters , , , and of the linear viscoelastic model should be determined by the least squares method using the given system of nondimensional equations and experimental data; see, for example, [17].
2.2. System with Two Parallel Viscoelastic Rods
Now, a different mechanical system is analyzed in order to take into account both different time instants of the beginning of deformation process for different frontal parts of a vehicle and frontal wheels blockade. Two viscoelastic rods of different properties are attached to the block of mass , as in Figure 2, where the initial conditions for the impact are the same as in the previous case. Because of different lengths of the rods, the motion during collision is divided into several phases. In the first (approach 1) phase () only upper rod deforms, while in the second (approach 2) phase () both rods deform until the block stops at , which represents the end of the approach phase. It is reasonable to suppose that frontal wheels may be blocked in the end of the approach phase, so additional (dry friction) force should be considered as well. From , there are two possibilities: the block will either stay captured or it will change the direction of motion and start moving backward. More about possible scripts of a block with viscoelastic rod impacting against a rigid wall or against another block is presented in [13, 18]. If the block with two viscoelastic rods of different lengths is not captured in the approach phase, then the third (rebound) phase () occurs during which the block moves backward until the impact is finished in . The possible scripts during the rebound phase are capture in the rebound phase with rods in a contact with the wall, capture in the rebound phase with one rod in a contact with the wall, or separation from the wall. The friction force will be taken into account for Contrary to the model with one viscoelastic rod, this model is nonlinear for two reasons: firstly due to the presence of the gap of the shorter rod and secondly for the presence of dry friction.
The system of equations describing the first part of the approach phase () for the model from Figure 2 is presented by (1) and (2) together with initial conditions (3) and restrictions (5). This phase ends at , when the shorter rod, with properties , , , , , and , gets in touch with the obstacle. Note that
During the second phase of motion () governing equations are presented by the following equations:together with the constitutive equation for the longer rod (2). The derivative order is in a range for the fractional Zener model, while for the standard Zener model. The corresponding initial conditions readThe elongation and the force in the shorter viscoelastic rod are denoted by and . Restrictions to the parameters of that rod are analogous to the ones described by (5) in which is replaced by . This phase ends when the block stops, that is, at the time instant .
Fundamental axiom of dynamics for the third phase () readswhere represents the dry friction force in a set-valued form; see [19]Contact forces and satisfy constitutive equations (2) and (11) with corresponding restrictions. In (14), and stand for the dry friction coefficient and the normal contact force between the block and the ground surface. During this phase, the friction force is taken as . The initial conditions for the third phase read
Governing equations in dimensionless form for each of the three phases of impact, with dimensionless quantitiestogether with (6) read as follows.
Equations (7)–(9), with dimensionless restrictions in the form of (5)2,3, describe phase 1.
Constitutive equation (8), together withwhere the dimensionless initial conditions are of the same form as (12), describes phase 2.
Constitutive equations (8) and (18) and fundamental axiom of dynamicswith initial conditions in nondimensional form, which are the same as in (15), describe phase 3.
Restrictions to the model parameters of the fractional Zener model for the second and the third phases are the same and they are of the form (5)2,3 for both the longer rod and the shorter rod, where is replaced by in the latter case.
Parameters , , , , , , , and should be determined on the basis on experimental data, using the corresponding system of dimensionless equations for all three phases.
3. Solution
The two posed impact problems are similar to each other but procedures for their solving differ significantly. They both include system of equations of integer and fractional order, but unlike the first problem, the system with two viscoelastic rods belongs to a class of multiphase dynamical systems because different external loading corresponds to different motion phases. The Grünwald-Letnikov definition of fractional derivative can be used instead of the Riemann-Liouville definition (4) in the case of the presented real physical problem; see [20]. In this section, we present numerical algorithms, including Grünwald-Letnikov representation of a fractional derivative, for obtaining the solutions of both impact problems.
3.1. Solution to the Case with One Viscoelastic Rod
In order to solve systems (7) and (8) with initial conditions (9), we introduce a new dimensionless variable , instead of , in order to remove nonhomogeneous initial condition (9)2, according to requirements which arise when the Riemann-Liouville fractional derivative is used. The governing system of equations with the variable and with homogeneous initial conditions readsDiscretization of the time domain is performed , , where is a time step. The first- and the second-order derivative of a function over , at time instant , can be presented aswhile the Grünwald-Letnikov fractional derivative of order , , readswhere coefficients are given asin the form of the recurrence relationship. By solving systems (20) and (21), using (23)2, (24), and (25), the numerical solution is obtained and it is presented by the following expressions:for , where according to the initial conditions (22) and to (23)1. The position of the block can now be easily calculated by the equationFor calculation of unknown parameters , , , and of the fractional Zener model, the position of the block should be defined as the function of time and those parameters By forming the functionwhere stands for the number of experimental points taken into account and represents the experimental values of the position at time instants , , the model parameters can be calculated by the least squares method. Note that in the case of the standard Zener model.
3.2. Solution to the Case with Two Viscoelastic Rods
The motion of the system with two rods is divided into three phases, which need to be treated separately.
In the first phase of impact () only the longer rod is in contact with the wall and dynamics of the system during that phase is equivalent with dynamics of the model with one viscoelastic rod. Therefore, solutions of the governing equations for that phase are presented by (26) and (27).
The second phase () starts at when the shorter rod contacts the wall with the velocity and starts to deform. Hence, this inhomogeneous initial condition should be removed by introducing a new variable , instead of Note that although the initial conditions given by (12)1,2,3 regarding the longer rod are inhomogeneous for , that rod started to deform at , which was before , and initial conditions are homogeneous when is replaced by because the history of deformation of that rod is taken from the beginning of deformation process at . The governing system of equations with the variables and includes (21) and the following equations:together with the initial conditions for the second phase:It should be noted that due to the relation between variables and which reads Numerical algorithmtogether with expression (26)1 represents the solution to the governing system of equations for the second phase, for , where According to the initial conditions (31) the initial values of the variable and the contact force for the second phase are Coefficients are presented by (25) in which replaces
The position of the block during this phase, which ends at when the velocity of the block becomes zero , readswhere stands for the step number in which the block stops.
The third phase () represents the process of moving the block backward. Governing equations, with respect to variables and , for this phase include constitutive equations (21), (30) and fundamental axiom of dynamics in the formwhile the initial conditions are
Solution to systems (21), (30), and (34) with initial conditions (35) is given in the form of the numerical algorithm presented by expressions (26)1, (32)1, andfor , where is the step number in which the impact process is finished.
Note that the values of discretized functions at the beginning of the current phase (the second or the third phase) correspond to their values at the end of the previous phase.
Parameters , , , , , , , and , of the fractional Zener model for both viscoelastic rods and (the dry friction coefficient), can be obtained by finding the minimum of the following function:where stands for the position of the block as the function of time and model parameters. Note that in the case of the standard Zener model.
Parameter identification results for the model with one viscoelastic rod and for the model with two rods are shown in the next section.
4. Results
Experimental data from [6], where a car-to-pole impact experiment was performed, are used for parameter identification in this paper. Data such as initial velocity of the vehicle, duration of the approach phase, and duration of the crash test are known from the experiment and they read , , ms, and ms, respectively. Values , of the block position are chosen at ms time intervals of the impact duration, where for the case with one viscoelastic rod; see Figure 1. In order to apply the fitting procedure (28) together with (26) and (27), experimental values are converted into dimensionless form and the following values of model parameters are obtained: , , and , where the parameter is taken to be equal, , for simplicity reasons, without loss of generality, because the cross sectional area, the modulus of elasticity, and the length of the rod, which form the dimensionless group (see (6)6), are not explicitly included into the nondimensional numerical algorithm. All results are presented in dimensionless form. Numerical and experimental results are shown in Figure 3, from which it is seen that agreement is not satisfactory. One should note that the obtained derivative order is very close to , which is characteristic for the standard Zener model of a viscoelastic body.
In the authors’ opinion, the reason for disagreement could be due to the simplicity of the model, that is, the presented mechanical model with one viscoelastic rod does not take into account all physical effects of the very complex impact process analyzed here. Some effects, such as different material, size, position, and shape of frontal parts of the vehicle and different time instants of the beginning of their deformation, as well as frontal wheels blockade during the collision, might be taken into account by the more complex mechanical model, which is presented in Figure 2.
Since there are nine unknown parameters within the model with two fractional viscoelastic rods (seven for the classical Zener rods), the process of their determination by the fitting procedure (37) would be too demanding. Hence, it would be useful to determine some of them before applying (37). Parameters and are taken to be equal, , for the same reason as in the first case. Further, from Figure 3, it can be noted that is almost straight line from the beginning to certain time instant, say , when its slope starts to decrease significantly. That could mean that at other frontal parts of the car start to deform. If we take this assumption and determine time instant dimensionless time units from experimental curve, then it is possible to calculate model parameters of the longer rod using (26) to (28) within time interval because only one rod deforms in that case. The remaining set of parameters is to be determined by finding the minimum of the function (see (37)), which now depends on four instead of nine unknown parameters. After applying the fitting procedure with described numerical algorithm for the first phase only, with five experimental points presented in Table 1, the model parameters of the longer rod modeled by the fractional Zener model are obtained: , , and . Agreement between experimental and predicted values of during the first phase is shown in Figure 4.
Unlike the nonsatisfactory agreement for the model with one viscoelastic rod presented in Figure 3, excellent agreement is shown in Figure 4 for the first phase. In the case of the standard Zener model, where , the model parameters for the longer rod read , Agreement with experimental results for this model during the first phase is excellent, like in Figure 4.
Parameters , , , and of the fractional Zener model are calculated by the least squares method with ten experimental points given in Table 2, using (26)1, (32), and (33) for the second phase and (26)1, (32)1, and (36) for the third phase. The obtained results for the fractional Zener model are , , , and , while for the classical Zener model they are , , and . Position of the block , as a result of the numerical simulation during the whole impact process, is presented in Figure 5(a) for the fractional Zener model and in Figure 5(b) in the case of standard Zener model. In these figures, very good agreement with experimental results is achieved for both viscoelastic models. Normalized root mean square error is about 1% for both cases. In Figure 6, velocity of the block and forces in the fractional viscoelastic rods are shown for .
(a)
(b)
(a)
(b)
During the first phase, when only the longer rod deforms (most frontal parts of the vehicle), the slope of the curve is steep and almost constant, Figure 5(a), and the velocity slightly decreases while the force in the longer viscoelastic rod increases, Figure 6(a). During the second phase, which starts at , the slope of apparently decreases. That is, in the phase during which the shorter rod also deforms (parts of the vehicle behind the most frontal ones), contact force in that rod impulsively increases, while the velocity of the block abruptly decreases. After , the third (rebound) phase begins, where decreases. Both contact forces decrease during the rebound phase and at the force in the shorter rod becomes , which means that the rod is separated from the wall, while the longer one is still being in the contact. At the end of the experiment, at , the longer rod still does not separate from the wall, , and the velocity of the block is not zero, , which means that impact process does not finish at (it finishes at when either contact forces or the velocity of the block becomes zero). The time instant when the impact process finishes can be determined by the presented mechanical model with two rods, but it will be omitted due to the lack of experimental results for The same qualitative analysis holds for the classical Zener model of viscoelastic body. One should note that, after the separation from the obstacle, the rods will take the original configuration if tends to infinity. This is due to the fact that the presented model with two viscoelastic rods is piecewise linear. During the approach phase, the proposed model is a bilinear one because of the difference in the lengths of the parallel viscoelastic rods, which start to deform at different time instants. Some other bilinear models are presented in [21, 22]. However, we analyzed the behavior of the vehicle only during its contact with the obstacle, which is a short time duration process, after which, for , the deformations of the rods are nonzero.
5. Conclusions
In the analysis of the frontal impact of a car against a rigid obstacle, the car is modeled as a rigid block with a viscoelastic layer. Two different mechanical models are considered in which the deformable layer is presented either by one or by two viscoelastic rods. In the latter model, the rods of different lengths are parallel to each other. Different properties of the viscoelastic rods are due to different material, size, shape, and position of frontal parts of the vehicle. Both the classical and the fractional Zener model are used for modeling the viscoelastic properties of the rods and the results are compared. Dynamics and deformations during the impact process are described by systems of differential equations, which includes fractional derivatives of both contact forces and deformations of viscoelastic rods in the case of the fractional viscoelastic model. For the classical Zener model standard procedures are used for obtaining the solution, while in the case of the fractional Zener model the numerical algorithm, based on the Grünwald-Letnikov definition of fractional derivatives, is set up for solving governing systems of equations. Values of the model parameters are calculated by parameter identification using the least squares method, on the basis of experimental data from [6]. On the one hand, the results for the model with one viscoelastic rod do not show satisfactory agreement with experimental data, although more complex viscoelastic models than the models used in [6] are applied. On the other hand, the results accomplished by the model with two viscoelastic rods show good agreement with the experiment, much better then results obtained using standard viscoelastic models with integer order derivatives presented in the paper [6]. It is known that fractional derivative models are better than integer order ones for modeling of viscoelastic properties of real materials. However, a very good agreement with experimental results is achieved by both the classical and the fractional Zener model. According to the results of this paper, it can be concluded that taking a more complex viscoelastic model does not necessarily lead to significant improvement of the agreement between theoretical and experimental results for the crash test. Much better results are obtained by changing the mechanical model, that is, when two parallel viscoelastic rods of different lengths are used instead of one, regardless of which viscoelastic model is applied, the classical or the fractional one. As it is commented above, different impact scripts are possible, even the one without the third phase, that is, the one in which the block is captured in the end of the approach phase. It would be interesting to investigate and numerically recover all possible scripts for the system with two viscoelastic rods, as it was done for the model with one rod; see [13, 14, 18].
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
Funding for this work was provided by the Serbian Ministry of Education, Science and Technological Development (Project no. 174016).