Journal of Sound and Vibration (1992) 159( 1), 57-83
IDENTIFICATION
OF A CRASH
MODEL
J. P. MIZZI INRETS (Znstitut National de Recherche sur Ies Transports et leur Stkuritt!), 109 avenue Salvador Allende, Case 24, 69675 Bran Ckdex, France AND
L.
JEZEQUEL
ECL (Ecole Centrale de Lyon), 36 avenue Guy de Collongue, 69130 Ecully, France (Received 13 July 1990, and injinalform
22 July 1991)
Our knowledge of the behaviour of passenger cars and road safety systems during a crash trial is based on experimental studies. A survey was carried out on the modelling of the front compartment of a passenger car: the model will make it possible to enlarge the conclusions drawn from a test by extending the results to different situations. We have designed a mathematical spring-masses model which simulates the behaviour of a passenger car during various frontal crash configurations. However, the main difficulty is to know perfectly the laws of behaviour of the springs. That is why an identification methodology was envisaged from the configuration of the experimental results. To know the vehicle’s real behaviour during a crash trial, it is necessary to have experimental devices which make it possible to rebuild the space kinematics of the components. We thus designed, in each case, suitable acquisition and processing software. Different non-parametric and parametric identification methods were then tested on simple and then complex models. The results have permitted us to determine which is the best adapted to solve our problem.
1.
INTRODUCTION
The mathematical modelling of a vehicle when crashing into a fixed or deformable obstacle must allow both the deceleration field inside the car and the structural and kinematic behaviour of the vehicle to be reconstituted. Although a detailed view of the behaviour can be obtained by using finite element modelling, the more global spring-mass modelling has the advantage of requiring considerably lower computing costs, and therefore of allowing the behaviour for many types of impact to be analyzed quickly. The model is adapted to three test configurations which are representative of what happens in real accidents. Despite how variable the entry data may be, its processing is remarkably flexible and rapid. However, to begin with, we are restricted to using a simplified working diagram of the springs: the ends are not clamped. The mechanical equations were obtained by using the orthogonal group of a transformation as a Lie group. The main difficult in executing the procedure is to determine the laws governing the behaviour of the connections. These can first be estimated either by finite element calculations, or by static or dynamic tests. This approach can be improved by using an identification methodology, provided that the behaviour in a crash trial is well known The displacement at a point on the outside of the vehicle is reconstituted on film from data supplied by at least two video cameras set at different angles. The absolute acceleration of reference marks within the vehicle in the laboratory frame is recorded by accelerometry. 57 0022460X/92/220057+ 27 .WS.OO/O
@ 1992Academic
Press Limited
58
J
P. MIZZI
AND
L. JEZEQUEL
As there are a priori so many visible exterior points, most of the information available for identification is in the form of displacements. It is useful to have information on speeds and accelerations. We have therefore tested and used a method based on spline functions to filter and derive the experimental noisy data reconstituted from the film. A mathematical description of a system can be specified from an identification, by using a priori and/or experimental information. The difficulty of the inverse problem derives from its intrinsic unstable nature. That is why a methodology which is best adapted to the resolution of the problem with which we are dealing must be determined. It would appear to be imposible to completely determine the stresses in the connections from filmed data. However, if we take the hypothesis that for each connection the force depends linearly on a certain number of parameters, the estimation can be made from the trajectory. A nonparametric method has been considered, based on the breaking down of the response into a double series of orthogonal polynomials. Conversely, if an initial estimate of the set of parameters is available, a parametric identification can be preferred. Several have been tested. Once their respective stability and strength had been determined on one-dimensional test models with one and two degrees of freedom, predictability and reliability trials in the presence of noisy data were made.
2. MODEL DEFINITION
2.1. FINITE ELEMENT MODEL The new digital methods, such as the finite element method, backed up by increasing high-performance computing systems, mean that problems which could not even be tackled until a few years ago can now be resolved. Thus studies on a whole vehicle have been made [ 11. However, an impact concerns both non-linear behaviour and dynamic stress, so the meshing which underpins the computations must be based on a structure plan which is closed in its final form, and the resulting computation, even when the equipment and algorithms used are sophisticated, requires lengthy processing. MODEL 2.2. SPRING-MASS During impact tests, some constitutive elements are considerably deformed while others remain practically undefonnable. The vehicle can therefore be represented by a system of rigid bodies with deformable joints between each. The ensuing model gives kinematic information at discrete points of the structure only. Several authors [2-121 in the U.S.A. have already considered the problem of modelling the front compartment as bodies and joints. However, they considered only collisions with symmetrical orthogonal obstacles, and American vehicles are designed very differently from European ones. In France, the RNUR [ 1l] has used this method to study the behaviour of some of its vehicles. Very few studies have been made on two- or three-dimensional models [ 131.
2.3. BALANCE EQUATIONS Once we know the geometry of the front compartment and the type of obstacle, the vehicle can be modelled (see Figure 1). Each body, being a rigid solid, is characterized by its geometry and its initial kinematics. Each joint, undergoing as it does only normal stress, is defined by the knowledge of its points of connection to the bodies, and the law governing its mechanical behaviour, which is the result of the assembly of several real parts. It accounts for possible clearances, buckling phenomena and loading/unloading cycles. For our first approximation, we assumed that its static representation could be a succession of
IDENTIFlCATION
OF A CRASH
MODEL
Figure 1. Lumped mass/non-linear spring model.
affine straight lines (see Figure 2). The effective dynamic force is obtained by applying a correction a : F&n = F,,,( 1 + ah
(1)
Here ci is the joint deformation speed. The methods usually used to form equations which govern the movement of a complex system require the use either of general theorems, or of analytical mechanics techniques (Lagrange or Hamilton equations). Our idea is to find a method which allows these equations to be written more synthetically, so as to leave the computer to carry out the t Stress
Deformotion
(b)
Figure 2. Laws of behaviour of the non-linear springs. (a) Only compression; (b) traction and compression.
60
J. P. MIZZI
AND
L. JEZEQUEL
greatest proportion of the developments. This method [14-161 is based on two essential concepts: the organization of the joints must be described with the help of a multi-graph ; the dynamics of each solid of the system must be expressed in a very synthetic form by using the differential geometry techniques of Lie groups. That is to say, for a system of n rigid bodies subjected to joints, the positions of the system at every moment are described by elements (s,, . . , s,,)ES~. The bodies of the system which defines the assembly I= (122, ” . 9n) are interlinked by joints, the existence of which is described by a graph structure on I. The description of the interactions transmitted by the joints to one or other extrema, are opposed by virtue of the action/reaction law, which justifies the use of an oriented graph. As we wish to distinguish various joints linking the same two bodies, we shall use a U multi-graph. Obstacle W is then considered either as a body or as a set of bodies which is peculiar in that if it is fixed, the links with W impose conditions on the movement of the bodies, and if it is movable the links impose conditions on the relative movements. If a reference configuration I’[for each body of the system is chosen, the movement of the system is described by a family Di 1[4= {orthogonal transformation group in an Euclidian space} and si = Diri )
iEI.
(2)
The Di)s are then the unknowns of the problem, which must be determined in time. The hexadimensional inertia force which acts on body i is expressed in Lagrangian terms by [I41 .z’=-Hj(v:)-[vf,
Hj(v;)],
(3)
where d denotes the Lie algebra associated to D, HiEcT(d) are the hexadimensional inertia operators of body i in reference configuration ri, vf is the Lagrange speed, and [. , . ] denotes a Lie bracket. A law governing the interaction force is associated with each arc a of U, as a function of the positions and speeds of the bodies at the extremities of a. Fajk is then defined as the wrench of the forces exerted on bodies j by the intermediary of a. The objectivity principle [ 16-191 shows that Fnjk must be of the form Fajk=Di Yy,j,JDF’Dk,
D,"(Uk-
vi)),
(4)
where E is Euclidian affine space, EqE is the set of equiprojective field wrenches, Fajk denotes an arbitrary application of D x EqE + EqE, D,r’Dk is the relative displacement of k to j, D,T’(vk- Vj) is the relative speed of bodies k and j in the reference configuration of body j, and Vi- Vk is the relative Euclidian speed field of bodies j and k. During the compression strains develop, the resultant of which shall be termed Rj for body j: Rj= T,+ N/.
(5)
Nj is an unknown reaction, which does not build up any power for the admissible kinematic speeds. The balance equations are written as Hj(UT) + [I$‘, Hj(Vj)]=F,+
Faj+ Rj,
(6)
IDENTIFICATION
OF A CRASH
61
h&DEL
where Fj is the given stress (in our case it results from the weight of the body) and Foj are the stresses in the a springs in contact withj. The equations which determine the N/stresses must be added to those above. However, in some cases, these equations overlap with the movement equations. Before impact, the vehicle is assumed to have a uniform movement. In the case of a one-dimensional link a, and a plane movement deformation d, one has d=L--
IMa, i)--4o,_i)ll,
(7)
where 1, is the unloaded length of the link, 11.II denotes the norm in the vector space E associated to 8, ao(a, i) is the position of attachment of a on body i and ~(a, i) = DiuO(cx,i), ao(aJ) is the position of attachment of a on bodyj and a(a,~)=Q~,(a,~), and, as Di is an isometry, d=l,-
IMa,
O-Df’&da,j)ll,
(8)
where J={(DY’(bi_bj)[4a,
i)l/(ao(a,
i)-Df’Dj4a,j)ll,
i)-D;‘Djao(a,i)))/llao(a,
(9)
and ( * , - ) denote a scalar product in E. Vi-4 is applied to ~,,(a, i). The resulting stress for a spring a is thus written as @==F~,,{a(o,
~)-4~,.O)llb(~,
G--4a,i)ll,
~~=@a~4~,jlg(~),
(l&11)
where g(i) is the position of body i’s centre of gravity. The results show that this modelling, thanks to the existence of computing tools which are both user friendly and modular, and to the use of efficient integration algorithms which can be adapted to the model types and the impact configurations, does not require model conception and equation calculation times which are excessively long. In the case of simulations against rigid orthogonal obstacles, the results correlate well with those of the literature (see Figure 3).
3. ANALYSIS OF TRIAL TESTS
3.1. METHODS OF OBTAINING DATA The various mechanical parameters of the general springs can be determined with greater accuracy by experimentation, which also enables critical comparisons to be made with simulation results. However, this does require the right experimental equipment to gather the information and analyze it. We can have at our disposal high frequency cinematography and accelerometry. With cinematography, the three-dimensional movement of a vehicle can be fixed onto flat images. The initial trajectory of the points can be reconstructed by reversing the image formation process using a mathematical model. Since the mid-1960’s, many papers have been published on the use of non-metric cameras for photogrammetric applications. However, the spatial kinematics of a deformable object cannot be reconstructed from a single shot. The object must be filmed by at least two cameras set at different angles. Synchronizing the two views is a problem that has several possible solutions. We have chosen a timemarker of the images at a given frequency, using photo-emissive diodes. The data can be processed either analytically or numerically. We preferred the former, the final precision obtained is greater. These methods have not often been used to reconstruct an object’s movement in space. The technique usually used is that of the “multiplier”, which can measure only plane or almost plane movements. We can choose either iterative or direct analytical methods. The former require knowledge of the optical properties of
62
J. P. MIZZI
AND
L. JEZEQUEL
160 100 z E z 5 6
60 0 -60
-160
60
I
I
70,
1
I
60
I
I
At++
t
+++++
I
I ++++++++++++++
I
(c)
i
50
6 ._ z L 2 x
40 30
s
20 ‘0 i _.lsm
O8 l.c
20 I
60 1
40I
00 1
I
100
Time (ms)
Figure 3. Simulation results; plane model with four bodies and 14 springs; displacements (a), speeds (b), accelerations (c), us. time (in ms). 0, Body 1; x, body 2; +, body 3; *, body 4.
the cameras, which is difficult to obtain for non-metric cameras. That is why our final choice was to use a direct analytical method, as it studies each perspective beam from a camera separately. The only constraint is that the control points should be well spaced out to reverse the perspective projection. The efficiency of the technique can be limited by errors arising from the degree of precision with which the control points are located, the quantity of the data recording and processing measurements, and the spatial configuration of the cameras and their synchronization. To obtain the absolute acceleration of a point moving through space from a fixed frame, it is not enough to place a triaxial accelerometer at that point, for the acceleration data obtained are absolute only in the frame linked to the sensor. Given the type of information collected during an impact test, piezoresistive acceleration sensors must be used. Various methods [8-12,201 allow the angular speed of a rigid solid to be determined from data supplied by accelerometers, as well as the acceleration at its centre. However,
IDENTIFICATION
OF A CRASH
MODEL
63
the use of any piece of equipment implies experimental uncertainties. In particular, when an accelerometer is excited in a direction, the output signal should be a cosine function of the angle between the excitation input axis and the sensor’s sensitive axis, but this is not the case. The output/input ratio, expressed as a percentage of sensitivity following the measurement axis is defined as the transverse sensitivity coefficient. Other kinds of errors, due to linearity, hysterisis, resolution, etc. coexist. Fortunately, most of these can be reduced to a minimum. In addition to errors attributable to the sensor itself, there are errors due to its position in space, as it is not reduced to a material point. 3.2. EQUATIONS OF OBSERVATIONS BASED ON THE FILM The principle of this method has been desribed by Fligor [20] for the plane case, and by Abdel Aziz and Karara and Walton [21] for the three-dimensional case. Let M(X, Y, Z) be a point of the object space, and M(X, y) its image. The basic relation of the direct analytical method is x=
AX+BY+CZ+D
IX+JY+KZ+L
EX+FY+GZ+l’
(12)
‘=&X+FY+GZ+l’
where A, B, . . . , L are coefficients which are functions of the internal elements and of the spatial orientation of the camera. An object point A4iof known absolute co-ordinates (Xi, YipZi) gives two equations, so that with six control points the down projective coefficients can be determined with a slight over-determination. The control points should be well spaced out, and not share the same plane. It is best to have these points as far apart as possible within the depth of field limits of the camera. Adding other points to the six minimum ones should be done cautiously, for it will not necessarily result in a better estimation of the projective coefficients. It may lead to mutual dependence. Even when the image plane of the measuring system is not parallel to that of the film expression (12) is still valid. On the other hand, the lenses used may give rise to different types of error: either a deterioration in the image’s quality, or a change in the image’s geometry, depending on the case. The distortion of the second type is of the greatest importance in photogrammetry. There are two parts to it: optical distortion (a property of the lens) and off-centre distortion (a property of the image, due to poor optical element alignment). We must add in the deformation of the film during processing and its positioning in the measuring system. The transformation is then written as
X=
A’X+B’Y+
Y=
C’Z+ D’+M’Au+N’Av+~,~~
E’X+F’Y+GZ+l I’X+J’Y+K’Z+L’+P’AU+Q’AU+~~~. E’X+F’Y+
,
(13)
G’z+ 1
Note that it is no longer possible to determine the coefficients directly, as Au, Au, Ax and Ay intervene in the expression, and they depend on the image co-ordinates of the principal point of the camera. To use this method, the projective coefficients must first be determined. We have designed a mobil device (see Figure 4)-mobile in that first it does not remain on or around an impact zone, and second in that the same reference points are used whatever the test configuration used. The object and image frames can be any, provided that some conditions
64
J. P. MIZZI
AND
L. JEZEQUEL
Figure 4. Spatial cinematography.
of definition are laid down. Expression (12) applied to all the control points supplies the matrix system, [d]C=x
(14)
where C is the coefficient vector, [A] is a matrix dependent on xi, yi, Xi, Yi and Zi and x is a matrix dependent on xi and yi. As we have over-determined, an estimator of C is calculated by least squares. JIx- [n]C 11is minimized by an orthogonalization method. The matrices and vectors are transformed by a Pi sequence of Householder’s matrix. The coefficients of expression (14) are calculated in two stages. First, the above 11 coefficients are calculated. With them, an order of magnitude of the co-ordinates of the image at the camera’s principal point can be determined. Then the 16 coefficients of expression (13) can be calculated by using an iterative method. Once the calibration has been made for each of the views, the trajectory of a point of the space can be reconstructed from the measurement on the film, X‘ Y Z_ where the matrix p and the vector I- are functions of the projective coefficients of the image co-ordinates of the point on the various views used. Measurement of the views can be made manually by using a digitizing table. As there was so much information to process, we preferred to collect it semi-automatically, by digitizing the images and recognizing shapes from test patterns fixed on the vehicle which defined the required points.
IDENTIFICATION
3.3.
EQUATIONS
FOR
OBSERVATIONS
OF A CRASH
WITH
THE
MODEL
65
ACCELEROMETERS
The calculation of the absolute acceleration of a point M of a rigid solid can be written in the following form if, at points Pi, expression ( 16) is applied : fj=f(J+
[Pj]ui + [O][W]pi.
(16)
Here fi is the absolute acceleration of Pi, i$ is the absolute acceleration of the centre of the frame linked to the solid, [Pi] is an antisymmetric matrix depending on the position of point Pi, pi is the position vector of point Pi, [w] is an antisymmetric matrix depending on angular speeds and ci, is the angular acceleration vector. All these expressions are expressed in the moving frame fixed to the solid. If this equation is applied at three points, on each of which is mounted a triaxial accelerometer sensor [21], an equation of the following form is obtained, d,+[P]ti+u=O.
(17)
where R. is a nine-line column vector, [P] is a nine-line, three-column matrix and u is a nine-line column vector. However, due to measurement inaccuracies, equation ( 17) is not proven. An error vector can be introduced, leading to &+[P]w+u=&.
(18)
To minimize the quadratic error E~CE,R,, and w must be found, C being the covariance marix which is assumed to be equal to the unit matrix. Taking POto be the barycentre of points PI, PZ and PS affected with identical weights, and POto be the centre of the moving frame, one obtains [P]‘& = 0,
(19)
whence c3 = -([P]‘[P]-‘)[P]‘u.
(20)
Relation (20) gives us the absolute accelerations at one point of the solid in the moving frame. To have them in a fixed frame, the transformation between the two frames must be determined for all time. This can be done classically by using Euler’s angles. However, a different approach has been used in navigation [4]. Taking Z, to be the components of a point M in a moving frame, and Zj its components in a fixed frame, one has Zi= DgZj
where Do is orthogonal.
(21)
Upon supposing
.n,=-D&
(22)
so d,=-D$j,
aI= DboiDg,
L$= DgUie
(23-25)
The orientation of the moving reference frame compared to the fixed frame can then be calculated by integrating (25). The problem of the non-commutativity of the finite rotations leads to an incorrect D, matrix.
66
J
P. MIZZI
AND
12. JEZEQUEL
The suggested method [4] isolates the vector which depends on the non-commutativity. Taking R to be the moving frame, it turns with angle v, around a unitary axis u with time. Upon supposing
where cp is the direction vector of u, one has
(27) where [q.J is an antisymmetric matrix associated with q. The vector cpsatisfies the differential equation
fj=o+&o+-!- l-pcl+cosq) 2 sin ij9 v2i
CpAcpAW,
1
for a v, which is different from kn/2 and k’~, k, k’eZ, with Z the set of relative integers. The connection existing between the preceding rotations and w is represented by the last two terms of equation (28). A device was designed (see Figure 5) in order to reconstruct the absolute acceleration in a fixed frame. As the device is carried on cars, it must stand up to the mechanical strains brought about by the impacts without any notable deterioration or deformation which might lead to important measurement errors. The device must be held in the chosen area by an interface which will not significantly modify either the resistance or the inertia. The car/device connection is optimized so that no string vibrations which would lead to measurement errors are introduced.
Figure
5. Accelerometers
inside a vehicle.
IDENTtFtCATtON
OF A CRASH MODEL
67
Comparisons and cross-checking of results of the two cinematographic and aczelerametric methods were made. The results show that there is good correlation, both as regards the shape of the curves and the values obtained {see Figure 6). “..
50 40 30 20 10 0
-101 0
f 50
I
100
,
I
I
150
200
250
I
Figure 6. Cross-checking cinematography-acceterometry ; displacements us. time; (a) and (c) points on the vehicle sideways (cinematography); (b) vehicle centre of gravity {acceleromerry).
3.4. ESTIMATION OF THE DERIVATtVES To estimate a differential of a function from noisy sampfed data is an ill-posed problem. Differentiating is the same as solving a particular Fredholm integral equation of the first kind. Tihonov fl] has studied how to pose the problem correctly by regularizing the equation. In fact, this boils down to searching for a family with a well-posed problem optimization parameter, so that the corresponding solution family estimates the solution of the integral equation as well as possible. The final object is to determine an optimal low filtering which will suppress all the frequency components added by noise. If the power spectrum of the signal and noise were known, a Wiener filter would be suitable for an optimal least squares separation. As we never have this information, additional constraints must be included. The movement is presumed to be smooth enough not to condition high frequency components; otherwise that would mean that there are large inertia forces. The sampling frequency must be high
68
J. P. MIZZI
AND
L. JEZEQUEL
enough to ensure that the signal and noise have negligible components outside the Nyquist band. Two kinds of method coexist: frequency methods. the aim of which is to determine the position and shape of the transition curve between the signal and the cut-off frequency; and time methods, which select a smooth function by estimating its value at the sampling points by at least squares criterion-the functions determined from splines are very efficient. Anderssen and Bloomfield [9] have perfected a differential procedure using Cullum’s regularization process [5] which was itself based on Tihonov’s paper [22]. It takes as an hypothesis that the variables are stationary and equally spaced and uses a spectral analysis procedure. It generates the no differential directly from sampled noisy data, and so minimizes the increase of error effects. The recording times are generally short, so the data cannot be seen as stationary or weakly stationary. Transformations can then be envisaged to render them almost stationary without increasing error sensitivity or resulting in the derivatives not being obtained. A function ge W: with g(0) =g( 1) is given. We have sought the functions fE WY which minimize C(a, f ),
cl..f)=,,nia~:t[~~'f(x)dll;io~f,,:. few:,
(2%
where a E [0, 1] fixed. Wf is the space of twice differentiable functions with second derivatives square summable. I H(t-s)f(s) ds=Af, (30) g(t) = s0 where H(t -3) is the Heaviside unit function. By using the likelihood function, we have obtained the parameter a. A filter window best adapted to the data may be defined. The time method is based on the determination of a spline function which minimizes a function of the form
'p(f)= ' (f"'(t))* dt+i pi(zi-f(ti))*v sa i-1
(31)
where zi are the sampled data at points ti and f is the function being sought. The solution of this minimization problem is an adjustable smooth natural spline of order q which depends on a parameter p = (pi, . . . , p,)‘. The problem is to obtain a good value for this parameter. Several attempts have been made to determine it automatically according to the data. Some methods give only the theoretical behaviour of p as a function of generally unknown quantities. Others allow a practical determination. Various studies have been made. The most feasible is that of Craven and Wahba [23], improved by Utreras [24], who considerably reduced the computation time when the amount of data is large. The smooth values at the sampling points are linked to the noisy data by the relation x=A(p)z, where A is a (n x n) matrix. To find the spline parametered by the optimal p involves minimizing 2 v(p)=,$,
k(A(p)-z)i
1 -it’d
ii
*
(32)
The denominator takes the longest time to evaluate. Utreras [24] has shown that determining the eigenvalues of A means solving for the eigenvalues of a variational equation. We
IDENTIFICATION
OF A CRASH MODEL
69
applied this formulation to cubic (q = 2) and quintic (q = 3) splines for unequally sampled data. Determining the derivative is automatic, as the parametered spline sought is defined from divided differences. Here again, the accumulation of error effects is avoided. A comparison of the two methods from noisy test data shows that the time method determines a better filter than the frequency method with the same boundary conditions of the interval studied (see Figure 7). The primary derivatives obtained are smoothed. To obtain better secondary derivatives, the order of the spline function must be increased. The results obtained are of satisfactory quality. This methodology was also applied to real data from cinematography. Comparisons with parallel external measurements also showed very good matching of the results (see Figure 8).
Figure 7. Smoothing and differentiation. (a) Noisy initial function; (b), (c) smoothing and differentiation by frequential method; (d), (e) smoothing and differentiation by spline method.
70
J
P. MIZZI
AND
I_. JEZEQUEI
(b)
0
Figure
8. Smoothing
I
I
50
100
and differentiation
4.1.
DEFINITION
OF THE
INVERSE
150
200
on cinematographic
4. IDENTIFICATION
i
-
OF
250
data.
(a) Displacements
THE
MODEL
vs. time; (b) speeds.
PROBLEM
When we have a physical system to be studied, the first step is to parameterize the system : that is, to define a minimum series of parameters of the model, the values of which characterize it completely. The second step is to discover some physical laws with which, for given parameter values, predictions of the results of observable parameter measurements can be made. This is the direct modelling stage. If the physical laws are known precisely, the predictions should match the measurements, unless there is a modelling error. If the laws are not precisely known, they can be more accurately determined by using the actual results of certain observable parameter measurements in order to deduce the real values of the model parameters. This is reciprocal modelling. An inverse problem can be solved by using various techniques. They may be analytical : all the possible values of the parameters may be scanned systematically or randomly by Monte-Carlo methods, which have the advantage of not using any linear approximation. However, normal problems do not generally have few parameters, and determining the probability density is expensive in computing time. Then the only approach is to define a strategy to attain the point which maximizes the probability density quickly. The most used criterion is a Gaussian one, for most sources of error can be modelled by using Gauss functions. The resolution methods then are based on the classical optimization theory. As regards our models, the linkage strains cannot be completely determined at each instant, as there are not enough dynamic equations, unless one presumes, as we did, that the generalized force depends linearly on a certain number of internal parameters. The differential equation giving the evolution can always be expressed as 2==f(x)u+g(X)
for
tE[ti,
&+I],
(33)
where X denotes the set of state variables of the system, CJdenotes the set of parameters which enter into the definition of the springs and f(X) and g(X) denote non-linear functions of X. The impact lasts for as long as the obstacle exerts a force. The inertial speed, as can easily be imagined, is of prime importance. When it decreases, the maximum global crushing of the whole vehicle decreases. This implies that, from one test to another, the spring effort
IDENTIFICATION
OF A CRASH
MODEL
71
paths will be different. On the other hand, the duration of the impact and the moment of maximum deformation have little to do with the initial speed, for they are characteristics functions of the system. The various parameters which characterize a behaviour law have preferential zones of influence, and this is where they will be most easily identifiable. This is a problem of dynamics, and the influence of the deformation speed is not negligible in the system’s response. There will therefore be damping throughout the collision, which will be superimposed on each area of influence of the other parameters. That is why we consider that we know its value by linking the static force to the dynamic force by relation (1) defined in section 1. It can be noted that the spring parameters influence the duration of the impact and the moment of maximum deformation. Another important factor in the resolution of inverse problems is the steps of processing of the various equation resolutions. It should be fine enough to be able to follow the various influence zones of the parameters. However, the optimum step depends not only on the springs themselves, but also on the initial speed. The presence of parameters which enter non-linearly into the movement equation make the identification problem more difficult. Because the impact phenomenon is so brief, we cannot obtain the equivalent of several periods of the system to study its response domain. Moreover, the excitation can be controlled only by varying the initial speed. Finally, as each spring has a fragmented affine behaviour law, and therefore its own slope change sequence, the multiplication of the springs in the model requires that a common sequence be determined, so that between two consecutive elements the set of linkage parameters is constant. 4.2. NON-PARAMETRIC IDENTIFICATION To avoid the drawback of having to determine complex spring behaviour laws, we considered non-parametric identification. This consists of first seeking to determine a function h(X), where X represents the state of the system; the search is made in a function rather than a parameter space. It is not necessary to make any hypothesis about the structure of the springs: the number of degrees of freedom of the system only needs to be known for the answer to be estimated by a double series of orthogonal polynomials. We choose Chebyshev T, polynomials for the following reasons: they satisfy orthogonality relations; they satisfy T,,(cos 0) = cos n&--this is useful for the integrations; they facilitate the identification of non-linear behaviour. The method used is derived from recent research [25]. To test its feasibility, we used a mode1 made up of a spring and a body (see Figure 10 of section 4.3) against a rigid orthogonal obstacle. The state variable of model X is written
[I x
x= ; . x
Its dynamics are governed by the equation h(x, a) = -5
(34)
The function h is known at a number of discrete points which correspond to a speed and displacement measurement during a test with various initial speed configurations. The coefficients of the series curve are calculated from the interpolation or extrapolation of the measurements. Non-parametric identification is difficult because of this determination.
72
J. P. MlZZl
AND
L. JEZEQUEL
One has h(a, 6)&z’,
b’)= c” 2 Ck,Tk(U’)T,(b’), k=O
(35)
I=0
where the Tj are Chebyshev’s N, or Nb degree polynomials
4 ck’=(~+60,)(*+~O&)(~~-!)(~~-l)i_;,
j-1
The P, are the discretization points M,Mb in number, and the formula for the c&l)s is given by approximating h^in the least mean squares sense. To determine the C&;s, h must be known at points P,,which is not a priori the case unless the measurement points coincide with the discretization points. Two cases must be considered depending on whether Pji is close to a measurement point (that is, there is a triangle which is not too flat, made up of three experimental points around it), or whether it is far away from it. The results from measurements points derived from digital simulations we made were then exploited, for four different initial speeds. The precision of the non-parametric identification is of course conditioned by the convergence of the double series of polynomials, which is a function of the number of MJ4, points where It is sampled, of the number of NM measurement points present, and of the maximum degree of polynomials used (N,, N/,). The method was applied for two sets of data: ( 1) N, = Nh = 8, M, = Mh = 15, NM = 20; (2) N, = Nb = 2, A4, = Mb = 15, NM= 20. The results (see Figure 9) cannot be said to give satisfactory identification. Only the general aspect is correct. The quality of the approximation depends on the choice of the initial speed, because it is a function of the measurement points from which the Chebyshev coefficients were calculated. It takes a very long time to calculate these, for the whole phase plan must be explored in order to carry out the interpolations. Moreover, this methodology requires many test results to compensate for the information contained in the spring definition. As it costs so much to submit a vehicle to these tests, the fact that the initial estimation of the parameters is avoided does not justify its adoption. 4.3. PARAMETRIC IDENTIFICATION The parametric problem is expressed as follows. Minimize a functional
s T
J=
UX, U) dt,
(38)
0
where T is the final impact time, with strains X(0) = vo f
Jf=f(Jw+g(X),
U(t)=CJ
if
tE[tj-1,
tj],jE{l,.
. . ,tZ>, (39)
theunknownsbeingt,,...,t,_,,V,,...,U,withto=Oandt,=T. A one-dimensional test model with one spring only was defined (see Figure 10). It is represented in theoretical form. The spring behaviour law is also defined by a succession of fragmented affine straight lines (see Figure 2). This model is a special case of the general model which interests us.
IDENTIFICATION 0.0
73
OF A CRASH MODEL
1 ‘8:
,/
-1.2.
0.0,
#,’ ,/’
- 0 2 :,
: : >’
‘:3
-1.75.
:.. 0.0
05
1.0
1.5
0.0
2.0
0.5
/’
1.o
1.5
2.0
Figure 9. Non-parametric identification; one-degree-of-freedom model. Accelerations us. time (dotted line, experimental results) (continuous line, identification results). Set I : crash speed 1 for (a) and 1.5 for (b); set 2: crash speed 1 for (c) and 1.5 for (d).
dv &&----
l-
8 I-&-
f*
t
Figure 10. One-degree-of-freedom
model.
The various parameters are expressed in non-dimensional form by choosing T,= m as the time unit, and a nominal value VUas the speed unit. The values to be identified are a = k/k, ,
P =Wh
,
y=e/L,
6 =_L/h&,
(40)
with x, = x/x’,
(41)
where X’is the non-dimensional value of X (the displacement of the body’s centre of gravity). The reference law will be obtained from a set of parameters and an initial speed, by integrating the movement equation. We shall then use a set of parameters which will
.I. P.
74
MIZZI
AND
I... JEZEQUEI
represent an initial estimation. The various methods of parametric identification will be applied with these data as a basis, in order to find the reference set of parameters, or at least to obtain an heuristic law which is as close as possible to the reference one. 4.3.1.
Review of tested methods
The classical gradient method is applied where for the test model
s T
J(U, T)=;
--Jk,)‘Q(f(WJ+g(~)
tf(X)U+gV)
-.&x,4 dt.
(42)
0
In all the following, the final time of impact is presumed to be known. J then depends only on U. Moreover, as X depends on U, its gradient can be expressed only as a function of X and U, but a CX/alJ must be introduced. In the procedure used the time interval [0, T] was divided into segments. This results in J being replaced by T, L(X. U) df. J’(U, Ti)= (43) s0 By minimizing this new functional, a set of parameters Vi will be obtained, which will be used as the initial estimation for the minimization of J’( U, Ti+ ,). Its minimum will be sought by seeking the minimum fii of J’( Vi+ si(aJ’/aU)i):
To sum up, this method could be efficient, but its main drawback, even with a simple model, is its very long computation time. Another method based on the theories of dynamic programming, and described by DiStefano and Rath [lo], was tested. The spring parameters are integrated to the state of the system. We always end up with a first order differential system: with A=(X,
dA/dt=G(A)
U),
A(O)=Ao,
G(A) = [-F(X,
U), 01.
(44)
The initial state A0 is made up of the usual conditions to which an estimation of the parameters U is added. The problem is therefore to find the optimal state law of the system. A is such that k = G(A) and minimizes I(A)=
T(W-NA)‘(W-NA)dt+(A(0)-A,)‘P(A(O)-Ao), s0
(45)
with W= NA + E being the vector of experimental observations with error E; N is a matrix. The numerical integration of the evolution equation gives a family of solutions as an heuristic law of state A form. The solution which minimizes I is made more precise by determining A at every instant. The optimal value A(T) is called the optimal least squares filter of A. Where only one pass of duration T does not conclusively improve the value of the parameters, several passes are made. The results are not very satisfactory. The problems encountered are due to the fact that the interaction time is limited, and the derivations from the parameters are based on instant values. A new procedure was worked out to avoid, as far as possible, the drawbacks of the two above methods; namely, the fact that the gradient method is long and cumbersome, and
IDENTIFICATION
OF A CRASH
75
MODEL
that the impact is too brief to determine the optimal filter. The effects of the parameter variations should be studied on cost functions added up over the total impact duration. To increase the speed of convergence as many cost functions as parameters should be considered. These costs will favour precise intervals corresponding to the parameters’ fields of influence. A number m of basic independent functions B=&(t), i= 1, . . . , m, are constructed, and with them m cost functions can be determined in the form T Ji(X,
Cl)
Yi(X, U, t) dt
=
where K(X, U, t) = L(X, U)Bi(t).
s0
(46)
The basic functions chosen were of three types: power functions, Bi(t) = I’- l/T’-’ ; sinusoidal functions, Bi( t) = sin (nit/ST), so that they are positive throughout the impact ; and square Chebyshev polynomials, to prevent different sign divergences cancelling out. The power functions with a fine step give good results, but the sinusoidal functions are the most efficient. If in L(X, U) the speed and displacement, weighted compared to the acceleration, are considered, the results show that an appreciable improvement is made in the final estimation of the parameters. However, the method is still very dependent on the initial estimation. If we suppose that the threshold values of the behaviour law are known, then the determination of the parameters can be seen to be perfect, and practically independent of the initial estimation (see Figure 11 and Table 1). 4.3.2. Sensitivity method The preceding method of section 4.3.1 gives good results when the initial values of the set of parameters are not too far from the reference set (error of about 20%). However, the essential problem in the preceding procedures is how to determine the derivative of a function, in particular of the gradient of the functional to be minimized. The evolution equation is now written as X = F(X, U) with X a function of U, whence 2=4(U), and g
= (8~18~ )(8X/a u) +
aF/a u,
(47)
whence
1axau
d ax aFax aF ---_--+dt
[au
(48)
au’
This equation is the sensitivity equation of X for a variation of U. The criterion to be optimized is
s T
J(U)=
L(X, U)d~,~=joT((;~~+~)di.
(49950)
0
The series of parameters (Uk) converges towards a value which minimizes J and is defined by the following relation (obtained by using a Newton-Raphson method): -’ dJ uk+,
In our case, F(X, U) =f(X) in the form
=
uk-p
U+g(X)
(51)
,,(uk)*
on intervals [I~, ti+ ,I. The function L(X, U) is chosen
L(X, U)=f(X-X~~~)‘Q(X-X,,,)+~(U-
UW,)‘R(U-
Uinir),
(52)
J. P. MIZZI
76
AND
I.. JEZEQUEL
0,40.3-
-.~,._...h.-.-._.-.-.-.-.~.-.-.-. 0.2 -.’ 0.1 0.0
0
(b) 5
/ 10
I 15
I 20
_L.-.-,
I 25
I 30
I 35
I 40
1 45
1 50
-1.2-(c) 00
0.5
1.0
1.5
2.0
Figure 1 I. Integral filter method, crash speed= 1. (a) Alpha and beta parameters us. iteration number; (b) gamma and delta parameters vs. iteration number; (c) crash deceleration vs. time. Acceleration us. time: , dotted line, result with initial parameters; - - -, experimental result; --, result with identified parameters.
where X,, is the measured data of the state variable X and Uiniris the initial estimation of the set of parameters. We are thus led to solve the following system (53): d ax -=-aFax IaF dt [1 au ax au au’
(53)
77
IDENTIFICATIONOF A CRASH MODEL TABLE 1 Initial speed = 1. Weighting coeficients:
displacements = 2; speeds = 1
Thresholds Free
Clamped
Clamped
Initial ;
0.9
0.9
0.4 1.5
:
0.14 0.35
0.14 0.35
0.1 0.85
;
1.200 1.202
1.200
I.200 1.200
0.24 0.600
0.240 0.600
0.240 0.600
0.49 x 10-l 044 x 1o-6
0.49 x 10-I 0.19 x 1o-9
0.65 0.58 x 1O-9
Final
; Initial cost Final cost
However, the a priori instants (ti) are unknowns of the problem. The J’s minimum for each ti must be calculated. A variation of ti will bring about a variation of X. By setting Xi=X(ti)
and
GXi=X(ti+
6tf) -X(ti),
(54)
then, to first order, where U= Vi-1 if at,>0 and U= Vi if 6ti<0.
sxi=(f(xi)+g(Xi))Sti+6X(ti)
(55)
First we presumed that the threshold values were known. The differential system was then integrated on each interval for the one-dimensional model. The results are very conclusive (see Figure 12). The method converges very rapidly and very poor initial estimations can be made. As we have already noted, the initial impact speed strongly influences the identification of the parameters. At high speeds, the zones of influence of some parameters become extremely brief and their indentification is harder to obtain.
0.5
1.0
1.5
2.0
Figure 12. Crash deceleration, crash speed = 1; sensitivity method. Acceleration us. time: * . . . * , dotted line, result with initial parameters; - - -, experimental result; -, result with identified parameters.
J. P. MlZZl
78 4.3.3.
AND
I,. JEZEQUEL
Optimal control method
Another method consists in solving the following non-linear programming problem (P) : /u=f(x)u+g(X),
U=lJ,
.w)=~o,
if lE[ti, tp 8 I].
I- I,
. II _ 1; (56)
minimize with respect to U f,t I J(U)=
t,+l).
(57)
MX, U, A, t) = L(X, U, 0 + k’(.f‘(X) +g(X)).
(58)
L(X, U, t) dt+ ‘f’(X(t,+,),
s I, The Hamiltonian
is now defined as the scalar function
It is presumed from now on that Hand L do not explicitly depend on time. d is any vector. To minimize J under the constraints defined by problem (P) means calculating
We have already seen that X depends on U. To avoid solving the sensitivity equation, we choose the vector d such that it satisfies
naIv/ax=o for
;i+aH/ax=o with
Therefore we suppose that at instant ti (dX),,=O written as
t=ti+I.
(60)
because X, is given. J’s gradient is then
fZ+l aH ’ dJ=
au
sof,
(61)
dUdt.
Now
aH aL z=(l+f(XYI
oi=[;+‘((~j+L~(X))dUdt.
whence
(62)
As above, the links of the interval are also unknowns of the problem. Here too we consider that the threshold value is known. As there are so many operations to be made, the optimization of too great a number of parameters simultaneously should be avoided. Upon taking L(X, U, t) in the form L(X, U)=~(X-x~.~~)‘Q(x-x~,)+f(u-
Uinit)‘R(U-
Uinir),
(63)
we must then solve the following system (S;): u+
adw ax
-
Ia 1 ’
tit s I
~=ftx)u+g(n
X(O)=x0
9
minimize dJ/d U =
((U-
Ui”,,)‘R’+;Itf(X))
dt.
0
(64)
In system (S;) the most difficult term to calculate is {af(X)l U. We must also simultaneously solve a differential equation with an initial condition and one with a final condition. However, in no case doesf, which is a sparse matrix, need to be explicit. The
IDENTIFICATION
OF A CRASH
MODEL
79
Figure 13. Crash deceleration, crash speed= 1; optimal control method. Acceleration US.time: . . . ., dotted line, result with initial parameters; - - -, experimental result; -, result with identified parameters.
term f (U- Uhi,)‘R( U- U,,) strongly conditions the estimation of the parameters U. That is why we shall first suppose that R = 0, but the a*H/aU* matrix is singular. The results show, as with the preceding method, that convergence is rapid and the values of the parameters obtained are of excellent quality (see Figure 13). 4.3.4. Comparison of the methods on a two-dimensional model The model with one degree of freedom is relatively simple. Yet we already have information on the qualities and faults of the methods used. We shall now consider two systems, the design of which is identical to the preceding system, put end to end (see Figure 14). It is made up of 10 parameters for which one has the following expressions, once the equations have been rendered non-dimensional : al =Wh,
PI =Wh,
YI =e/xu,
61 =f,lk,Xu,
(65)
a{=k;/kl,
p i =G/ki9
Yi = e’/X,,
Si =fI/WU,
(66)
jt=
m/m',
(67)
~=Wk~)l(m’/m).
For the gradient formula, the results (see Table 2) confirm and greatly amplify the poor quality and speed of convergence. The procedure is used with a 0.08 step and segments of 100 steps for a test at nominal speed. The optimal filter method was no longer considered, as it already gave the poorest results. For the integral filter method, the best suited basic function had first to be determined. The more efficient sinusoidal functions were adopted. We divided the impact time into 20 segments, cut up at most 20 times. The related results reduce the cost by about lo4 compared to the initial estimation. However, it is very dependent on this estimation (see
Figure 14. Two-degree-of-freedom
model.
80
J. P. MIZZI AND L. JEZEQUEL TABLE 2
Step of 0.008, slices of 100 steps; crash speed= I Reference narameters
YI
6, a2 P2
I.2 1.2 0.6 0.24 1.2 1.2 0.6 0.24 2.4 0.7
Initial parameters
From 1 to 100 steps
From 1 to 200 steps
From I to 300 steps
From I to 400 steus
I 1 0.5 0.2 1 1 0.5 0.2 2 0.5
I.05 0.984 0.686 0.223 1.08 1 0.5 0.2 2.05 0.785
I.05 0.978 0.641 0.22 I.07 1 0.5 0.205 2.05 0.761
I.05 0.969 0.646 0,192 1.08 0.995 0.534 0.203 2.05 0.745
I.05 0.969 0.64 0.195 I.08 0.997 0.527 0.197 2.05 0.742
Table 3). Moreover, its predictability is rather poor. We then considered a cost function where speed and displacement as well as acceleration were introduced. The results then show almost perfect predictability. However, the initial estimation is the delicate part of
this method. When it is poor, the method rapidly diverges, and the parameter values become eccentric. If a penalizing coefficient is introduced into the parameters the results improve (see Figure 15), but not to the point of allowing any initial estimation if good predictability correlations are desired. If the thresholds are fixed, in the case of a good initial estimation, the results are satisfactory. On the other hand, when it is poor only a slight improvement is obtained (see Table 4). For the sensitivity and optimal control methods, the results are identical. They show that, whatever the initial estimations, convergence is rapid and of good quality. However, to avoid having to estimate too many parameters simultaneously, and so to improve the speed and quality of the convergence, it was useful to break down the initial model into a set of substructures having a single body and several links. The various parameters are then identified by using the kinematic data of measurements corresponding to a body and to the linkage points. On the preceding model, the results (see Figure 16) showed an increase in the convergence speed and quality. This quality is such that the predictability tests give identical results to the reference results. TABLE 3 Integral filter
method: dependence on initial estimation Final parameters
al
PI YI
6, a2 P2
Y2 s2 rl P
From overestimated initial parameters of 20%
From underestimated initial parameters of 30%
0.967 1.170 0.530 0.195 1.047 1.023 0.616 0.205 0.769 2.4
0.693 0.840 0.516 0.230 1.089 1.903 0.709 0.517 0.41 1.981
IDENTIFICATIONOF A CRASH MODEL
81
Figure 15. Integral filter method, crash deceleration; poor initial estimation. - - -, Result with initial parameters; . . ., dotted line, experimental result; -, result with identified parameters. (a), (b) body I and 2 without penalty; (c), (d) body 1 and 2 with penalty.
TABLE 4 Integral filter method: dependence on initial estimations and thresholds Thresholds
Free
Clamped
Clamped
Parameters
Initial
Final
Initial
Final
Initial
Final
at
0.9 o-9 0.4 0.17 0.9 0.9 0.4 0.17 0.5 1.80
0.460 1.070
0.9 0.9 0.4 0.17 0.9 0.9 0.4 0.17 0.5 1.80
1.181 1.204 0.600 0.236 I.150 1.101 0.600 0.230 o-713 O-239
0.4 1.5 0.850 0.1 1.7 0.6 0.450 0.350 0.5 1.8
0400 1.500 0.600 0.080 1.700 0.600 0.600 0.340 0.5 I.8
PI YI
61 a2 P2
Yz 62 ? P
0.608 0.275 2.244 1.282 0.299 0.377 0.385 2.055
5. CONCLUSIONS
A system of rigid bodies subjected to joints has been proposed to describe the crash behaviour of passenger cars. The equilibrium equations of the discrete model are written with the help of the differential geometry techniques of Lie groups. The hysteretic behaviour of the joints takes into account clearances and buckling phenomena of structures. With the intention of obtaining the best information from the crash tests, a method of observations by camera based on the introduction of control points is used. Moreover, a numerical procedure giving the absolute acceleration in a fixed frame from triaxial accelerometer sensor data is proposed. A non-parametric identification method using the kinematic data has been used to describe the dynamic behaviour of vehicles. However, this needs too much experimental information, and it appears that it is necessary to determine the parameters associated with the joint behaviour laws.
J. P. MIZZI
82
-2.5-
AND
L. JEZEQUEL
-1.50-
-0.75
-
- 1.00 -
-2.5-
-1.25
-
-1.50
-
0
1
2
3
4
Figure 16. Crash deceleration vs. time, crash speed=2. ‘, Dotted line, result with initial parameters; --experimental result; -, result with identified parameters. Sensitivity method (a), (b) body 1 and 2; optkal control method (c), (d) body 1 and 2.
Several methods of identification have been proposed via the filtering approach or optimal control theory. It appears that it is possible to describe the crash behaviour of passenger cars with low-degree-of-freedom systems identified by tests. The model predictability has been analyzed by introducing structural perturbations. Such a strategy of modelling can integrate seat and passenger behaviour and provide a powerful tool to control the efficiency of the design of several road safety systems.
REFERENCES 1. P. GUELLEC 1981SIA 81-18, L’bvolution de la conception de la structure et des #nents de carrosserie. Modelisation de la caisse-methodologie et applications. 2. M. M. tiMAL 1970 SAE Paper 700414. Analysis and simulation of vehicle to barrier impact. 3. K. J. ASTROM and P. EYKHOFF 1971 Automatica 7, 123-162. System identification, a survey. 4. J. E. BORTZ 1971IEEE Transactions on Aerospace and Electronic Systems. AES 7 (1). A new mathematical formulation for strapdown inertial navigation. 5. J. CULLUM 1971 SIAM Journal of Numerical Analysb 8(21). Numerical differentiation and regularization. 6. Y. I. ABDEL-AZIZ and H. M. KARARA 1971Proceedings of the Annual Convention of the American Society of Photogrammetry, 26-29 January 1971, Falls Church, Virginia. Direct linear transformation from comparator-coordinates into object space coordinates in close range photogrammetry. 7. G. PICHON 1973 Groupes de Lie. Paris: Hermann.
IDENTIFICATION
OF A CRASH MODEL
8.A.I. KING, A. J. PADGAONKARand K. W. KRIEGER 1974 Second 9. 10. 11.
12.
13.
14. 15. 16. 17.
83 Annual Session, 6 December
1974, Highway Safety Research Institute, Ann Arbor, Michigan. Measurement of angular acceleration of a rigid body using linear accelerometers. R. S. ANDERSSEN and P. BL~~MFIELD 1974 Numerical Mathematics 22, 157-182. Numerical differentiation procedures for non-exact data. N. DISTEFANO and A. RATH 1975 Computer Methods in Applied Mechanics and Engineering 5, 353-372. System identification in nonlinear structural seismic dynamics. P. VENTRE, J. C. RULLIER and J. P. VEROLLET 1976 6th International Conference on Experimental Safety Vehicles, 12-15 October 1976, Washington, D.C. The behaviour of vehicles and occupants in asymmetrical frontal collisions. J. E. GREENE 1977 International Automotive Engineering Congress and Exposition Cob01 Hall, Detroit, 28 February-4 March 1977, SAE Paper 770015. Computer simulation of car to car collisions. M. BRENNAN, M. MACAULAY and A. WYNN-RUFFHEAD 1984 International Conference on Vehicle Structures, 16-18 July, 1984. Modelling the collapse of cars in asymmetrical barrier impact tests. D. CHEVALLIER 1984 Groupes de Lie et Mecanique des Systemes de Corps Rigides-Seminaire “Modelisation Mathematique”, Kassel, 1984. MacGraw-Hill. D. CHEVALLIER 1986 Cahier du CERMA no. 6, decembre 1986. Le principe des puissances virtuelles en mecanique analytique, ses relations avec la mecanique de Newton Euler. D. CHEVALLIER 1988 Cahier du CERMA no. 9, juin 1988. Principe d’objectiviti: et thtorie Newtonienne de l’inertie du solide rigide. D. CHEVALLIER and J. P. LEBACQUE 1987 Rapport de contrat INRETS-CERMA. Mise au point d ‘un modele theorique de deformation de I ‘avant d ‘un vthicule ltger en configuration plane lors d’un choc frontal.
18. J. P. MIZZI 1987 Presentation d’un Modele Plan a Liaisons non Encastrhes, Rapport intermediaire
NNS 8703, decembre 1987.
19. P. MICHELBERGER, J. BOKOR, A. KEREZTES and P. VARLAKI 1987 International Journal of Vehicle Design 8( 1). Identification of a multivariable linear model for road vehicle (bus) dynamics from test data. 20. P. D. FLICXIR 1967 Proceedings of the Annual Convention of the American Society of Photogrammetry, Washington, D.C., March 1967. Resection without a camera of stationary parameters. 21. J. S. WALTON 1978 Proceedings of the International Congress of Sport Sciences, Edmonton, Canada, July 1978. Close-range tine photogrammetry: another approach to motion analysis. 22. A. N. TIHONOV 1963 Soviet Math. Doklady 4, 1035-1038. Solution of incorrectly formulated problems and the regularization method. 23. P. CRAVEN and G. WAHBA 1979 Numerical Mathematics 31, 377403. Smoothing noisy data with spline functions. 24. F. UTRERAS 1981 SIAM Journal of ScientiJic Statistics and Computation 2(3). Optimal smoothing of noisy data using spline functions. 25. P. ARGOUL and L. JEZEQUEL 1989 Journal of Applied Mechanics 56,697-703. Improvement of a non-parametric identification procedure used in non-linear dynamics.