Numerical crash simulation of vehicle structures: a literature survey

Citation for published version (APA):

Bruijs, W. E. M. (1987). Numerical crash simulation of vehicle structures: a literature survey. (EUT report. WFW, vakgr. Fundamentele Werktuigbouwkunde; Vol. WFW-87.051), (DCT rapporten; Vol. 1987.051). Technische Universiteit Eindhoven.

Document status and date: Published: 01/01/1987

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

a literature survey

Report number: WF\V 87.051 Faculty of Mechanical Engineering

Author: \V.E.:M. Bruijs

Numerical crash simulation of vehicle structures: a literature survey / author: W.E.M. Bruijs. -Eindhoven: University of Technology, Department of Mechanical Engineering. - Ill. - (Eindhoven University of Technology research reports / Department of Mechanical Engineering, ISSN 0167-9708 ;

WFW 87.051)

Met lit. opg.

ISBN 90-6808-013-X

SISO 657 UOC 629.018(048.8)

This is the report of an investigation of literature dealing with the numerical simulation of the behaviour of vehicles under impact load. This study has been carried out on the basis of three codes that, among others, are used for numerical crash simulation:

- PAM-GRASH; -CRASHMAS; - DYCAST.

The finite element method for nonlinear dynamics is discussed firstly. After that, the results of the three codes mentioned above are presented. These results are promising, but the cpu times are unacceptably large for application in the early phases of the design process.

It should be noted that the codes used for crash simulation, are based on general purpose finite element codes. No special features for crash simulation are present in the three codes mentioned above.

ABSTRACT

CONVENTIONS

1 INTRODUCTION

2 FINITE ELEMENT CODES FOR CRASH SIMULATION 2.1 Introduction

2.2 PAM-CRASH

2.3 CRASHMAS

2.4 DYCAST

3 THE FI~ITE ELEMENT METHOD

3.1 Introduction

3.2 Derivation of the iterative form 3.3 Finite element discretisation 3.4 Introduction of the inertia term

4 NUMERICAL TIME INTEGRATION 4.1 Introduction

4.2 Explicit methods 4.3 Implicit methods 4.4 Discussion

5.2 Results achieved with PAM-CRASH 5.3 Results achieved with CRASHMAS 5.4 Results achieved with DYCAST 5.5 Discussion

6 CONCLUSIONS

-; a

a

N -;

a

N A A

A

C AT

II

a

II

\I

A

\I

A:B : vector

: column filled with scalars : column filled with vectors : 2nd order tensor

: nth order tensor

: matrix filled with scalars : matrix filled with tensors : conjugate of a tensor : transpose of a matrix : maximum norm of a vector : maximum norm of a tensor : dot product of two vectors

: dot product of a tensor and a vector : double dot product of two tensors : cross product of two vectors : dyadic product of two vectors

1 INTRODUCTION

Within the frame work of a cooperation between the TNO Road-Vehicles Research Institute and the Faculty of Mechanical Engineering of the Eindhoven University of Technology, an investigation of the numerical simulation of vehicles under impact load is being carried out. In order to make an inventory of the 'state of the art' of the existing codes for crash simulation, a literature study has been done. This is reported in this literature survey.

Only the three most commonly used codes for crash simulation will be discussed in this report. These codes are:

- PAM-CRASH (Engineering System International); - CRASHMAS (Industrieanlagen Betriebsgesellschaft); - DYCAST (Grumman Aerospace Corporation).

In chapter 2, these programs will be discussed. In chapter 3, the finite element equations for nonlinear dynamics

will

be derived. Various numerical techniques to integrate the differential equations obtained in chapter 3 are presented in chapter 4. In chapter .5, results of crash simulations using P AM-CRASH, CRASHMAS and DYCAST are presented. In closing some concluding remarks are made in chapter. 6.

2 FINITE ELEMENT CODES FOR CRASH SIMULATION

2.1 Introduction

In this chapter some finite element codes, that are used for crash simulation, will be discussed.

2.2 P AM-CRASH

PAM-CRASH is developed at Engineering furstem International (ESI), Paris, and is an explicit Lagrangian finite element code, specialised for crash simulation. P A:t-.1-CRASH is a simplified version of the EFHYD3D finite element code for the analysis of large

displacement dynamic response of coupled hydrodynamic, viscous fluids and elastoplastic continua, shell and beam structures.

The equations of motion are integrated using the central difference method. The structures can be discretised by:

- solid elements;

- plate or shell elements; - beam elements.

A contact-impact algorithm, which is based on the impact logic of Hallquist (1976), permits gaps and sliding along material interfaces.

The yield criterion used, is the von Mises yield criterion. The isotropic or kinematic hardening theory can be applied.

For more information it is referred to Haug, Arneaudeau, Dubois, de Rouvray & Chedmail,

(1983).

2.3 CRASHMAS

CRASHMAS is developed at lndustrie~nlagen Betriebs~esellschaft (IABG), Ottobrun and is an explicit finite element code, specialised for crash simulation. The equations of motion are integrated using the central difference method. The finite element library contains solid elements, shell elements, beam elements and nonlinear springs.

A contact processor of the master-slave concept (Hallquist, 1976) makes it possible to model impacts. Nonlinear material models can be used. For more information it is referred to Bretz, Jarzab

&

Raasch, (1986).

2.4 DYCAST

DYCAST (llinamic Crash Analysis of Structures) is developed at the Grumman Aerospace Corporation (GAC), Bethpage. DYCAST is an updated Lagrangian finite element program and an outgrowth of the PLANS system of finite element programs for static nonlinear structural analysis. The equations of motion can be integrated using explicit - modified Adams (Garnet & Armen, 1975) - methods or implicit - Newmark-,B or Wilson-0 (see section 4.3) - integration methods. The finite element library contains membrane elements. beam elements and nonlinear springs.

Contacts between surfaces can be described using gap elements. A maximum strain failure criterion can be applied.

3 THE FINITE ELEME:\,T METHOD

3.1 Introduction

In this chapter the finite element method, applied to nonlinear dynamics, is discussed. In section 3.2, the statical equilibrium equation of continuum mechanics is transformed into an iterative form. The finite element discretisation is introduced in 3.3 and applied to the iterative form of section 3.2. The result is an equation that can be used to calculate the unknown nodal position vectors by iteration. In section 3.4 the inertia term is introduced

and added to the iterative finite element equation. This leads to a differential equation. This equation can be solved by numerical integration. Some numerical integration methods are discussed in chapter 4.

3.2 Derivation of the iterative form

Consider a three dimensional body. The material points of this body are characterised by the material coordinates

~T

=

[(1 (2 (3]' The set of all material points is denoted by B. The position of the material pOints will depend on the time t. The position vector ~l of point ~ at time tl can be written as:

(3.1 )

The symmetric Cauchy stress tensor q at time tl will be a function of the position in the

body;

The fields xl and 0"1 are unknown. For the calculation of the unknown fields, several equations are available. The first equation is the equilibrium equation that has to be satisfied at any place in the body at time tl (Veldpaus, 1982):

(3.3)

where

PJ

is the mass density at time tl and <II is the force per unit mass at time t_l. The operator ~ is the gradient operator defined by: dx.~O = O(x

+

dx) - O(x), with 0 a certain plac€-{}ependent quantity and

dX

an arbitrary infinitesimal change of the position. Another equation that is needed to obtain the unknowns, is the relation between the deformation and the stresses, the so-called constitutive equation, that can generally be expressed by:

(3.4 )

where F(t) is the deformation tensor of the configuration at time t, with respect to a known reference configuration. Furthermore, the body can be subject to dynamic or kinematic boundary conditions: O"lenl is prescribed along a part of the boundary of the body and

Xl

is prescribed on the remaining part of this boundary (n_l is the unit outward normal vector on the surface of the body at time t

l).

Let VI be the volume of the body at time tr The requirement that (3.3) has to be satisfied at any place in the body is equivalent to the requirement that:

J

W.(~eO"I

+

_{Pl<Il)dV I} =0 V

w

(3.5)

VI

(3.6)

where

PI

=

0'1-°1

and Al is the outer surface of the body at time t_l. Equation (3.6) is called the

weak weighted residuals formulation.

Suppose an estimation of the configuration is

*

known at time t_l. The estimation of the position vector is xl' of the mass density

p,

of the

*

*

stress field

t,

of the loads q 1 and

P

l' The deviations of the estimations are defined by:

*

ql

=

Cll

+

oq

*

PI

=

PI

+

oj)

Substitution in (3.6) yields:

J

(vw)c:(t I

+ oO')dV

1

=

J

_{we(a l}

+

bq)(Pl

+ bp)dVI

+

VI VI

+

J

we(~l

+

bp)d,A₁ Al (3.7) (3.8 ) (3.9) (3.10) (3.11) (3.12)

The calculation of

bx

with this equation is as complicated as the calculation of xl using (3.6). If

bx

is relatively small, an approximation of

Ox

can be found by linearising (3.12). A difficulty is that the integrations must be carried out over the unknown volume VI and

*

*

surface AI' A transformation to the known ~olume VIand the known surface Al is

necessary. The integrands of (3.12) are then no longer considered to be a function of Xl but

*

a function of Xl' A transformation for the gradient operator will be necessary too, because

v

is measured in the real configuration at t = t l' Let (I

+

bB) be the deformation tensor of the real configuration, measured in the estimated configuration. The transformations of the integrals for a function cp, and for the gradient operator are given by (Breke1mans, 198.5):

J

c.pdV 1

J

¢det( 1

+

bB)dV

*

_I (3.13) V

_*

1 VI

*

*

with 1/{x 1 ) cp(x_I (xl)) ,

J

c.pdA 1

J

7,bdet(1

+

bB) 11(1

*

*

- ( ; - l

+

bH) en111dAl Al

_*

Al (3.14 )

*

*

with ill the known unit outward normal vector of the surface Al in the estimated configuration, and:

Substitution of (3.13), (3.14) and (3.15) in (3.12) yields:

f

{(~l~·)Ce(1

+

bHf1}:(d-_l

+

ba)det(1

+

bH)d~l

=

*

VI

=

f

_{we(a l}

+

bq)(Pl

+

bp)det(1

+

bH)d~l

+

*

VI

*

*

*

+

f

we(Pl

+

bp)det(1

+

bH) 11(1

+

bH)-cenlll dAl (3.16)

*

Al

If the estimated configuration differs only little from the real configuration, it must hold that IlbH11

< <

1. The following approximations are allowed if IlbH11

< <

1:

(1

+

bH)-l

~

1 - bH (:3.17)

det(1

+

bH) ~ 1

+

tr(bH) (3.18)

(3.19)

Substitution of (3.17), (3.18) and (3.19) and using the law of conservation of mass:

J

*

* * *

+

~·.(j)l

+

op)(l

+

_{tr(bH))(l - bH:n1nl)dA1}

(3.20)

*

Al

If

ox

is small, bH will be small. Furthermore,

bu,

<Xi and

op

are assumed linear in

Ox.

Linearisation of (3.20) yields:

*

V I

(3.21 )

The terms

oH, ou, oq

and

oj)

can be approximated by linear expressions (Brekelmans, 1985):

*

(-t -t)c

(3.23 ) * * * -l -l 3 -l-l)

oq

= Q-Ox

+

Q:(Vlox (3.24 ) * * * -l -l 3 (-l ~-l)

Op

=

P-Ox

+

P: Vlux (3.25) * * * *

The tensor 4M is defined completely with the constitutive equation. The tensors Q, 3Q, P

*

and 3p are dependent of the volume forces and the surface forces. They are not discussed in * 3* * 3*

more detail here. It should be noted that the tensors Q, Q, P and P are dependent on the

*

estimated configuration. The tensor 4M is dependent on the estimated configuration. and in general, on the deformation history of the material. Substitution of (3.22) to (3.2.5) and the use of:

(3.26)

4*

where Dl is defined by:

(3.27)

and:

(3.28)

f

-l * -l 3* ~ ---+ *

-

w.{P.8x

+

S:('V 18 x) }dA

_l

=

*

Al

(3.29)

Usually, the terms that are related to the external loads in the left hand side of (3.29) are omitteq. This leads to:

= -

f

(~I~')C:dId~I

+

f

w·~IPId~l

+

f

w·g_{1d1 I} (3.30)

*

*

*

VI VI Al

3.3 Finite element discretisation

The body B is divided in ne finite elements of a relatively simple shape. In this chapter isoparametric elements will be used. For each element a separate set of material

coordinates ,e is defined within a definition volume Be. Each element has a number of nodes, the position vectors of which are stored in the column

2

e( t). Let

2

be the column

containing the position vectors of the nodes of the construction and let the nodes of the elements coincide wi th the nodes of the construct ion. In that case

x

_N e can be seen as a subset of

x.

_N For the position vector of a material point of element e at time t can be written:

(3.31)

where p€ are interpolation functions. These interpolation functions have to satisfy certain conditions that are not discussed here. If the functions pe(~) of element e are specified and the nodal point positions

2

e are known, the form of element e is completely determined. If

*

*

the estimated volume of element e is V~, and the surface of element e is A~, (3.30) can be rewritten as:

(3.32)

It should be noted that only element boundaries that are part of the construction

*

boundaries, belong to A ~. Similar to (3.31) is written:

*

*

-+ _ IT,eT(re)-+e

(3.34 )

The functions ~ are chosen to be of the form:

(3.3.5)

The choice for w (in which we is the column that contains the elemental nodal values of w) N

is known as the method of Galerkin. The functions w have to satisfy certain conditions, that are not discussed here. The quantities in the integrand of (3.32) can be considered as

*

functions of the position vector}tl of the material points. Because of the discretisations

*

(3.33) to (3.35), the positions xl are known as a function of the material coordinates £. The quantities in (3.32) can therefore be regarded as being a function of £e. This means that the integrations can also be carried out over the definition volume Be of £e and the surface of this volume Be of element e. It is noted that Be is not the set of all material boundary points of element e, but only the subset belonging to the outer surface of the body. The

*

integral transformations and the transformations of ~ are:

*

3

*

b'!/l

*

*

3

*

b'!/l

... L'" - ...

t:leT... . h ,-Ie

L'"

-Vw = _{' 1} '1' _I

'F"C

_{. N} w = °1 w, WIt °1 _{N N N .} = ₁ '1' _I

'F"C

_. 1= I 1= 1 (3.36 ) (3.37) (3.38)

*

3

*

b'lil

*

*

3

*

b1/l

-+-; ~ -I .:. -I reT-I . h re ~ -I

-\lx = _{· 1}

'I·

_I

or::

. N x = N

°1

N x, WIt

°1

N = ·1

'I·

I

or::

.

I= I 1= 1

(3.39)

*

* .

*

For

sf,

!~r and ~li (i=1,2,3) it is referred to Brekelmans (1985). Substitution of (3.33) to (3.39) in (3.32) yields: ne

*

~

-I eT Ke ,,-Ie W • .vX = l tv - N e= (3.40) with: (3.41 ) (3.42) (3.43)

*

r

e

=

J

ifp)1

dBe

(3.44 )

Be

*

The elemental quantities x~, c5x~ and c5w~ can be regarded as elements of columns

*

containing the nodal quantities of the construction xl ' c5x_I and c5w 1. The relationship

matrix pe for each element:

*

*

-Ie €-I xl N = -P xl N (3.45) -Ie e -I

bx

= P

bx

N - N (3.46) (3.47) Substitution of (3.4.5) to (3.47) in (3.40) yields: Let: (:3.49) (3 . .50 ) (3 . .51) (3.52)

(3.53)

This must hold for all

w,

so:

N

*

* * *

K.b~

=

-i

+

E

+

r.

(3.54)

This is the iterative form that we wanted to derive. The integrals (3.41) to (3.44) can be calculated numerically.

*

In closing the iteration process will be described. Let ~ be the approximation after k iterations:

*

.... ....k x=x N N (3 .. 55 ) .... k

Let bx _N be defined as:

.... k+ 1 .... k + J; .... xk+ 1

x =x U

N N N (3.56)

With equation (3.54) it is possible to calculate

b~k+1:

(3.57)

With the known

b~k+

1 and (3.56) a better approximation

~k+

1 of the exact configuration is known. This process can be repeated until a convergence criterion is satisfied. For the exact solution ~1 within the class of interpolation functions it must hold that: _N

o

N (3 . .58 )

The components of the col umn

g

can be interpreted as the internal nodal forces. The components of the columns

k

_N and

r

_N can be interpreted as the nodal forces due to the external loads. If the desired solution is found, equation (3.58) must be satisfied.

*

If the calculation of Kk is carried out for each iteration, the iteration process is called full

*

Newto71- Raphson.

If Kk is calculated only once, equation (3.57) changes to:

(3 .. 59)

The iteration procedure based on equation (3.59) is called

modified lVewton-Raphson.

3.4 Introduction of the inertia term

If the inertia effects are included, (3.3) changes in the equation of motion:

(3.60)

In this case (3.6) must be changed in:

Consider the term

J

~·~lPl

dV

l· \Vith the use of (3.6) and (3.14) this integral can be VI

*

transformed to the integral over the estimated volume VI:

For the second equation in (3.62) the law of conservation of mass is used:

*

P = (PI

*

+

bp)det(I

+

bH

)

If the construction is divided in ne elements, it must hold that:

Analogous to (3.31),

x

is discretised as:

(3.62)

(3.63)

(3.64 )

(3.65)

where

2f

are the nodal accelerations. Substitution of (3.65) and (3.35) in (3.64) yields, using (3.36):

n n

ef ..

_{-t -t}

*

*

_e

e

_{-t e} T _{e -te} ..

L

wexlPldV I =

L

~

eM

e~

e=l

e=l

*

(3.66) VI with: (3.67)

Because it must hold that:

(3.68)

the mass matrix is independent of the configuration. Addition of the inertia term to (3 . .54) yields:

(3.69)

with:

(3.70)

The accelerations can be integrated numerically. Different numerical integration techniques are available. Some of these techniques are discussed in the following chapter.

4 NUMERICAL TIME INTEGRATION

4.1 Introduction

The differential equation (3.69), obtained in subsection 3.4 can be solved by numerical integration. In this section some numerical integration techniques will be discussed. Suppose that the configuration is known at times 0, At, 2At, ... , t. The purpose is to calculate an approximation of the configuration at time t+At, by direct numerical

integration. There are basically two methods for direct numerical integration: explicit and implicit methods. Using explicit methods, equation (3.68) is considered at time t, and used to calculate the acceleration at time t. Numerical integration of the accelerations yields the position vectors at time t+At: X(t+At). Using implicit methods, equation (3.68) is

considered at time t+At. In that case an iteration has to be performed to obtain a

satisfying estimation of the configuration at time t+At. The explicit and implicit methods are discussed more extensively in the sections 4.2 and 4.3 respectively. In section 4.4 the solution methods are discussed and some possible adaptations are given.

For reasons of simplicity of notation, some conventions are introduced. A right subscript denotes the time at which the quantity is considered. This means that gt+At are the accelerations at time t+At. A right superscript denotes the iteration in which the quantity is considered. So

~~

are the position vectors at time t after iteration k. The quantity

b~k

is defined by:

-+k -+k-l -+k

4.2 Explicit methods

Equation (3.69) is considered at time t, which is a known configuration, so

bx

=

o.

N N

Substitution in (3.70) yields:

(4.2)

The most common explicit formula for numerical integration of ~t is the

central difference

formula:

(4.3)

where

£

(~t)

is the error being made. It can be proved that

11£ II

is proportional to

~

t 2 (Veltkamp & Geurts, 1980/1981). Substitution of (4.3) in (4.2) and omission of

£

yields:

(4.4 )

If the matrix M is diagonal, time integration using the central difference formula does not

involve the solution of any equations.

4.3 Implicit methods

For implicit methods, the equation of motion is considered at time t+~t:

.=.tk k-1 -!k -!k-1 ,.+k-1 rk-1

(4.6)

If the modified Newton-Raphson scheme is applied, Kk-1 must be replaced by KO. One of

the implicit integration methods is the

Newmark-;J method.

The relationship for this integrator is:

Combination of (4.6) and (4.7) yields:

Substitution in (4.5) yields:

.... k-1 ~k-1 ""k-1 - gt+~t + -':t+~t + .tt+~t

Equation (4.9) can be used for an iterative method to obtain an estimation of the configuration at time t+~t. The velocities can be calculated from:

The parameter

a

is free to be chosen. A commonly used value for Q is 1/2.

(4.7)

(4.8)

( 4.9)

Other implicit integration schemes are for instance the

Wilson-O

and the

Houbolt

integration scheme.

Applications of these schemes lead to equations that can be used for iterative methods. These equations are similar to (4.9). For more detail about these schemes it is referred to Bathe (1982).

Consider (4.9) and choose

f3

=

~.

This corresponds to the assumption that the acceleration is linear in the interval (t,t+~t). Substitution yields:

(4.11)

If go and

~O

are given, Eo can be calculated using:

( 4.12)

Let the first estimation of

Elt

be given by the formula:

-! 1 -! ..; 1 2=+

~~t

=

~o + ~t~o + 2~t ~o ( 4.13)

This would have been the solution if an explicit method was applied. Application of (4.11) yields:

Let

~t

be small. Then

IIKll1 « - 4 IIMII.

This yields: ~t

After one iteration, the estimation for the position vectors is:

(4.1.5)

( 4.16)

F _{or t IS equatIOn It as een use t at} h· " h b d' h _{_.~} M ~ 1 _{= -g~t} -d

₊

r' _;!S~t 1

₊

r

_"'~t" 1 R _{eca t at} 11 h -d . 1 _{~~t IS t le}

solution when the finite difference equations are applied. The term

~ ~

2

(~1

-

~o)

is thus the first correction after one iteration. Substitution of 21t in the iteration formula, and calculation of

bE

3 with the assumption

IIKll1 «-4IIMII

yields:

~t

(4.17)

This process can be repeated until the convergence criteria are satisfied. The solution at t=~t is used to calculate a first estimation of the configuration at t=2~t:

( 4.18)

4.4 Discussion

The appealing features of explicit time integration are its simplicity and efficiency. Complex nonlinear models and phenomena. can be modelled easily. In exchange for its simplicity and efficiency, the explicit methods are only conditionally stable. This means that the stable time step is bounded. For a guideline to calculate stable time steps, it is referred to Belytschko (1982).

The implicit time scheme has the disadvanta.ge that the updating of the displacements involves an iteration process. An advantage of these methods is that unconditional stability may be achieved. This means that the time step is not limited by stability conditions introduced by the numerical integration algorithm, and can be chosen based on other criteria.

In wave propagation problems, like crash problems, the explicit stable time step is bounded by the time the wave needs to pass through the smallest element. This means that explicit

integration with a constant time step is often inefficient. Therefore, Belytschko, Yen

&

Mullen (1979) introduced a method which integrates different partitions of the mesh with different time steps.

Another possibility is introducing an explicit-implicit mesh partition, so that part of the

mesh is integrated implkitly and another part is integrated explicitly (Belytschko

&

Mullen, 1977). Belytschko and Mullen (1978) have shown that the stability limit of the time step is determined strktly by the highest frequency of the explicit partition.

5 RESULTS OF CRASH SIMULATIONS

5.1 Introduction

In this chapter, the results of the numerical simulation of a vehicle crash, using PAM-CRASH, CRASHMAS and DYCAST will be discussed. The results of crash simulations with these programs are included in this chapter.

5.2 Results achieved with P AM-CRASH

The application of PAM-CRASH to a crash simulation is shown in the frontal impact of a VVv-Polo on a rigid barrier (Haug, Scharnhorst & Dubois, 1986). The finite element model

(fig. 5.1) exists of 5661 elements (5555 elements and 106 beam elements). The interaction between the rigid barrier and the coach work is described using the master-slave algorithm of Hallquist (Hallquist, 1976). A variable contact between different parts of the

coach-work at places where contacts are expected, is made possible by defining master-slave surfaces at these places.

The time step was equal to 1.5 10-6 [s] and the impact time was 60 10-3 [s]. The cpu time was four hours on a CRA Y 1 super computer. For the results of the calculations, it is referred to Haug et al. (1986). The deformations are given in fig. 5.2.

In order to verify the results of the calculations, a crash of a VW-Polo was carried out. After the results have been filtered, the results of the crash and the crash simulation were compared. The finite element model turned out to be about 12% to stiff. For further details it is referred to Haug et al. (1986).

Fig. 5.1: P AM-CRASH, the mesh. Source: Haug et al. (1986)

SHELL-ELEPlEHTS 106 BEAPI-ELEPlEHTS 5661 ELEPlEHTS 5100 HOOES

Fig. 5.2: The deformations. Source: Haug et al.

Om. 18m' 30 m. 36m. 40ml

PAM-CRASH (lSI)

Fig. 5.2: PAM-CRASH, the deformations (continued). Source: Haug et al. (1986)

Oms "ms

ao

ms 36"" 40ms

'AM-ClASH (ESI) 'OLO 90

5.3 Results achieved with CRASHMAS

In

order to verify the program, CRASHMAS was used to simulate a frontal crash of a car against a rigid barrier (Bretz, Jarzab

&

Raasch, 1986). The mesh is given in fig. 5.3. It is remarked that construction parts which do not contribute, or contribute only little to the absorption of the kinetic energy, are not taken in account. The material behaviour was dependent of the strain rate. The mesh contained 2604 nodes and 2800 elements. The smallest element boundary was equal to 18 [mm]. The car, with a mass of 1140 [kg] had a speed of 12.5 [m s-l] when it hit the rigid barrier.

Fig. 5.3: CRASHMAS, the mesh. Source: Bretz et al. (1986)

The calculation was stopped when the main part of the kinetic energy was absorbed and when no more interesting data were expected from the absorption of the rest of the kinetic energy. The cpu time was 89 hours on a vector computer. It is remarked that the program

In fig. 5.4 the deformation of the mesh is given. For more information concerning the results of the crash, it is referred to Bretz et al. (1986).

Fig. 5.4: CRASHMAS, the deformations. Source: Bretz et al. (1986)

In order to verify the results of the numerical simulation, a crash test has been carried out. The results of this test were in reasonable accordance with the results of the calculations (Bretz et al., 1986).

5.4 Results achieved with DYCAST

DYCAST was used to simulate the 13.5 [m s-l] crash of the steel frame of a Chevrolet Corvette to a rigid barrier (Winter, Crouzet-Pascal

&

Pifko, 1984). The finite element model (fig. 5.5) included the frame plus the other structure, except for the fiberglass body. The model used 157 nodes, 220 elements and 597 degrees of freedom. The contact between

the structure and the barrier was modelled using gap elements. The model included nonlinear springs, the stiffness of which was taken from test data. For the steel, bilinear stress-strain curves were used up to the failure strain. The static values of the yield- and ultimate stresses were increased and the failure strains were decreased to account for strain rate effects. No changes were made to the properties of the aluminum, because of its lack of sensitivity to the expected strain rates.

The time step was equal to 50 10--6 [s] and the impact time was 0.1 [s]. The simulation consumed 200 minutes of cpu time on an IBM 3033 computer system.

The deformed model is displayed in fig. 7.6 at some time intervals. The results are not validated by a comparison with experimental data. For more information it is referred to \\Tinter et al. (1984).

5.5 Discussion

As can be seen from the results achieved with CRASHMAS and PAM-CRASH, the results of the simulations are in good agreement with experimental data. The cpu times are, even on a super computer, at the current state of the art, too large to apply numerical crash simulation as a design tool.

\\

,

'; \ I • t-~;,.""---"'-_

'.

. v' ~ .. I.d; : t _ '.-.

Fig. 5.5: DYCAST, the mesh. Source: Winter et al. (1984)

Fifo. 2 5.0. WI.'" vI "1'''1 • •• emenr lTI<;.Oei

-.'.:.."., .. ,..

...

"-"'-,

:. t

(.-- , , ' , - - - -. --...--~~ ----14 ;oJ "ft'.. - _ . I~ ::.. roo ... J -"IM;>;, ·r~ ___ --,",-~ / ..

-

... \- ... _-- - -r-=-'--':'-V- J -ll.4MI'li -14.1111i>~

1 -2il"'~lC III - 08i1~ .... to

>

-'. ;,.

..

"/ . / . /

---,

!

~ -.

-

.... ---.--~.---",. -- l ... , ~ --I~~"'~

----

J -".J - - b ':. ~t~ _{- - li.l:..} .. .. I. ~.... ~ . • • '" ~.;:. ",,' ~Il. _ _ .. _ _ _ . --.,. -~--

...

----

.. - - • I ' ... ., -~---,.-..

-\':':"

.

- w.l t.(;;i T • lOU "':.\ C. Ii.t

Fig. 5.6: DYCAST: the deformations. Source: Winter et al. (1984)

1

i

---1·,;. ,...

6 CONCLUSIONS

The results achieved with the codes mentioned in chapter 5 are promising. The necessary cpu time for a full scale impact of a vehicle front structure is unacceptably large for an application in the early stages of the design process.

In crash simulation, an explicit integration method based on the central difference scheme, is most commonly used. This means that the maximum stable time step is bounded by a stability criterion. The elements used in crash simulation are shell elements, beam

elements, solid elements and rigid elements. The contacts are described using gap elements or using methods, based on the impact logic of Hallquist (1976). It is not clear which

method is most commonly suited for crash simulation.

Finally, it is important to realise that the codes used for crash simulation are based on general purpose finite element codes. No methods, especially suited for crash simulation are used in the crash simulation codes.

Bathe, K.J. (1982).

Finite element Procedures in En~ineerin~ Analysis. Englevwod Cliffs, New Yersey: Prentice Hall. ISBN 0-13-317305-4.

Belytschko, T.B. & Mullen, R. (1977).

Mesh Partitions of Explicit-Implicit Time Integration. In J. Bathe (Ed.), Formulations and Computational Al~orithms in Finite Element Analysis. (pp. 673-{j90). MIT Press.

Belytschko, T.B.

&

Mullen, R. (1978).

Stability of Explicit-Implicit Mesh Partitions in Time Integration. International Journal of Numerical Methods in En~ineering, 12, 157.5-1.586.

Belytschko, T.B., Yen, H.J. & Mullen, R. (1979).

Mixed Methods for Time Integration. Computer Methods in Applied Mechanics and Engineering, 17, 259-275.

Belytschko, T.B. (1982).

Computer Methods for the Nonlinear Dynamics of Structures. In S. Holmes

&

K. Saczalski (Eds.), Crash Analysis Methods for Vehicle Structures. (pp. 101-131). Rome: ICTS.

Brekelmans W.A.M. (1985).

Niet Lineaire Mechanica: Numerieke Aspekten. Technische Universiteit Eindhoven, collegedictaat 4583.

Tagung: Berechnung im Automobilbau, 613, Wuerzburg.

Garnet, H. & Armen, H. (1975).

A Variable Time Step Method for Determining Plastic Stress Reflections from Boundaries. AIAA J., 13,532-534.

Hallquist, J.O. (1976).

A Procedure for the Solution of Finite-Deformation Contact-Impact Problems bv the Finite Element Method. Lawrence Livermore Laboratory, University of California, UCRL-52066.

Haug E., Arneaudeau, F., Dubois, J., de Rouvray, A.

&

Chedmail, J.F. (1983). Static and Dynamic Finite Element Analysis of Structural Crashworthiness in the Automotive and Aerospace Industries. In W. Jones & T. Wierzbicky (Eds.), Structural

Crashworthiness. (pp. 175-217) Butterworths, London.

Haug, E., Scharnhorst, T.

&

Dubois, P. (1986).

FEM Crash, Berechnung eines Fahrzeugfrontalaufpralls. VDI Tagung: Berechnung im Automobilbau, 613, Wuerzburg.

Veldpaus, F.E. (1982).

Niet Lineaire Mechanica: Basis. Technische Universiteit Eindhoven, collegedictaat 4582.

Veltkamp, G.W.

&

Geurts, A.J. (1980/1981).