Subcycling: implementation and comparison of two methods

Subcycling

Citation for published version (APA):

Bruijs, W. E. M. (1988). Subcycling: implementation and comparison of two methods. (DCT rapporten; Vol. 1988.020). Technische Universiteit Eindhoven.

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

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.

methods

Zarch 1988

In crash simulation, the finite element mesh is usually non-uniform. If explicit integration is applied - this integration scheme is often used in crash simulation - the small elements of the mesh determine the maximum stable time step. In the larger elements, a larger time step is allotaed. If subcycling is applied, the smaller elements are integrated with a smaller time step than the larger elements, This can yield substantial savings in computing time.

However, a problem occurs on the boundary of the area integrated with small time steps and the area integrated with large time steps. In this report, two alternative treatments to solve this problem are introduced, i.e. the interpolation and the integration method.

If these alternatives are applied to a s h p l e linear example. the maximum stable time step can be determined analytically. For more complex problefns, an analytical technique to determine the wiaximum stable time step is not yet available.

In order to compare the alternatives with respect to numerical stabilityl accuracy and computing time reduction. both alternatives are applied to a linear and a nonlinear example: the impact of an elastic and a plastic bar.

For the linear example the performance o f the alternatives is the same. For the nonlinear example, the performance o f the interpolation method is better than the integration method. with respect to as well numerical stability as accuracy and computing tirae reduction,

ABSTRACT COMVENTIONS 1 INTRODUCTION

2 PROBLEH DEFINITION

3 TWO ALTERNATIVES

3.1 The Information Flow Diagram

3.2 Information Flow with Subcycling

3.3 Interpolation 3.4 Integration 4 IMPL~MENTATION 4.1 Interpolation 4.2 Integration 5 STABILITY ANALYSIS

5.1 Straight Explicit Integration

5.2 Stability of the Interpolation Hethod

5.3 Stability of the Integration Hethod

6 INPACT OF AN ELASTIC BAR

7 IMPACT OF A PLASTIC BAR

8 DISCUSSION

REFERENCES

x : column filled with scalars

,-..

A

: matrix filled with scalars

aT : transpose of a column

N

-

A-' : inverse of a matrix

Ixl

-

: 2-norm o f a column

1 INTRODUCTION

The finite element mesh of a vehicle structure for crash simulation is, like many other engineering applications non-uniform. In such a mesh, the smaller

or stiffer ele-nents establish the maximuìrt stable time step. The larger or softer elements of the mesh are usually integrated with this time step toe. This time step is much smaller than the maximum stable time step allowed in the larger or softer elements. It Bould therefore be useful to integrate the smaller or stiffer elements with a smaller time step than the larger or softer elements. This is called subcycling. Subcycling may lead to savings of the computational costs.

In order to quantify the savings in computational costs, consider a mesh consisting o f n, small elements and n1 large elements, It is assumed that the material properties of the small and the large elements are the same, Let the time step for the small elements be At/m and for the large ele3ents At. Let the process time be T. Let the computing time necessary to calculate one cycle for one element be Tc. The computing time with subcycling: T,. is then given by:

T

Ts = ns

rn

TC + nl

2

TC

and the computing time without subcycling: Tws, by: m

T ,

, = Tc(ns

+

n,) Let the efficiency be defined by:

Eff =

2

Substitution o f (1.11 and

(12)

yields: 1 + -

1i,q

ns Eff = nl (1.1) (1.4) "1

-

O, or? in In fig. 1.1, the efficiency is drawn as a function o f

-.

nl

computational costs are obtained. The savings are larger when

-

"1 is larger? If

-

-

other words! the mesh contains only small elements! no savings in

"S nS

nS so subcycling yields the largest computational savings when the mesh contains many large elements and only few small elements.

Fig. 1.1: The efficiency as a function of the ratio between the number o f

large and small elements.

In the above considerations, the fact that subcycling costs some extra computing time to treat the boundary between the elements that are and are not subcycled, is not taken in consideration. As will be seen later in this report, this can be ignored compared to the savings in computing time.

In crash simulation, the elements must be small in areas vhere large plastic deformations are expected. In these areas explicit integration will be

employed! because the plasticity algorithms make implicit integration procedures unfavourable. In the areas of the large elements implicit or explicit integration procedures can both be used successfully, However, in this report. it is chosen to integrate the area of the small elements arad the area o4 t h e large elements explicitly with the central difference method.

In chapter 2 . the problem occurring on the boundaries between the subcycled and the not subcycled area is discussed. In chapter 3 , two alternative methods to treat the boundary are introduced. using information flow

diagrams. In chapter 4 , some implementation aspects will be discussed- The stability behaviour of the two alternatives is analysed in chapter 5 for a simple linear example. In the chapters 6 and 7, the alternatives are

compared with respect to their numerical stability and accuracy? using a linear one-dimensional problem in chapter 6 and a nonlinear probleiì? in chapter 7 . In chapter 8 f the results of this comparison are discussed.

2 PROBLEH DEFIMITIOM

Consider the following set of differential equations, representing the equations of motion of a structure? after spatial discretization with the finite element method:

..

Ex i- fint = f

N ,ex t

-

( 2 . 1 )

*.

In this equation

g

is a lumped mass matrix, x are the nodal accelerations, f are the internal nodal forces (a function of the nodal displacements, ,int

the nodal velocities and the time):

N

,int f = -int f

(x.X,t)

,

and fext are the external nodal forces, (a function of the nodal displacements and the time):

,ext f = ,ext f (x.t)

( 2 . 2 )

( 2 . 3 3

The set of differential equations (2.1) can be rewritten in the following form:

..

x

= flx.X,t) n N N N where o ( 2 . 4 )

If subcycling is to be applied to the set ( 2 . 4 1 0 this set is divided in two parts a and b. Part a will be integrated using the time step At and part b will be integrated using the time step dt/mp with m a positive integer:

Linearization of ! 2 , 6 ) and (2.7) around the point x1 and

i,

a t time tl yields :

where: df 'c tl -3.

ax

x = Xl'X = Xl't = A = n - C I N

-

and ( 2 . 8 )

Because the equations of motion (2.1) result after spatial discretization with the finite element methodl the matrices

Al

and €i1 will have a band structure. Suppose that the column

x

consists of p elements and the column

, a

xb consists of q elements. Four areas can be distinguished in the matrices -1

A and EZ3.:

'c

The matrix

El

will have a similar structure as the matrix

Al.

If the first p

equations are integrated Eith the time step At and the last q equations with the time step At/m, it can be seen from the structure of the matrices

Al

and

€+, that the first p equations are coupled with the last q equations and vice versa. (The areas 3 and 4 o f equation (2.101 contain nonzero elements.) It is necessary to develop techniques which account for the coupling between the two sets of equations ( 2 . 6 ) and ( 2 . 7 ) .

c

-

k -k -k k -k O -k k -k

-

\ \ \ \ \ \ \ \ \

0

-k k -1; -k k -k -k k i

-

3 TWO ALTERNATIVES

In this chapter two alternatives to account for the coupling between the two equations ( 2 . 6 ) and (2.7) will be discussed. But first the information flow diagram of a one dimensional mesh will be derived. These information flow diagrams will be used to illustrate the difference between the alternatives. 3.1 The Information Flow Diaaram

Consider the equations of motion o f a structure after spatial discretization with the finite element method:

The equations (3.1) will Be integrated using the central difference method:

..

where x are the accelerations at time t, = nAtl

agn+l

are the displacements a t time tn+l = (n91)At etc, Consider (3.1) at time t, and substitute (3,2f:

J

--

Suppose that the mesh consists o f linear rod elements with element stiffness matrix :

(3.41 -e

*

=

i*

-1

The stiffness m t r i x of a one dimensional mesh will have the following form:

(3,59

f - K x

-int

-

--

( 3 . 6 )

Suppose that the mass matrix has a diagonal form. The pth component o f the column

xntl:

xgtl can be calculated using:

Suppose that the displacements of all nodal points are known up to and including time t,. The displacements at time tntl can be calculated using

( 3 . 7 ) for each nodal point. This process is illustrated in fig. 3 . 1 . In this figure, an arrow pointing from point A t o point B means that the

displacement of point A is necessary to calculate the displacement of point

B.

time

T

tn-1

p- 2 p-1 p t l p92 pt3

Fig, 3.1: The inforination flow diagram for a one dimensional mesh with linear elements.

If a nonlinear set of equations is considered, the sase information flow diagram can be applied with the following adaptations:

-

An arrow pointing from point A t o point B means that the displacement and the velocity of point A must be known to calculate the displacement of point B.

-

For the calculation of the displacements at time tnt1 the displacements and velocities at previous times are necessary.

For reasons of surveyability, in the following the arrows pointing from time tnml t o time will be omitted.

3.2 Information Flow with Subcyclina

A s mentioned before, different time steps will be used to integrate the set of differential equations. Suppose that part a will be integrated using the time step 5t and part b will be integrated using the time step At/m. In fig. 3.2, the information flow diagram is given. In this information flow diagram the parameter m has the value four. It can be seen that for node p

information of the displacements at the times t 1, t 2 and t 3 is

necessary. This infomation can be calculated using different nethods. The sections 3 . 3 and 3.4 will discuss two alternatives.

n+H n+z n+q time

i

O O p-2 P-1 B 4 b a p+1 p+2 p-i.3 4 b

The displacement for node p at times tntí., t

linear interpolation between the time t, and tntl (fig. 3 . 3 )

2 and t 3 is calculated using

~ n+q n+T ‘nti t 3

ntH

t 1

n+Z

tn++ time %-I I O I p-2 p-1 P p+1 p+3 4 P 4 b

*

a

3.4 Intearation

The displacements o f node p at times tn+ t

3At

using numerical integration with the time steps T F 7 and -5- (fig. 3 , 4 ) .

2 and t 3 are calculated

n+q At A:+4 t 3

n+T

timie

i

p-2 O O O O P-1 P p+l p+2 p+3

*

b V

In this chapter the implementation of the two alternatives presented in 3 . 3

and 3.4 will be discussed. Consider again the equations of Botion of a structure after spatial discretization with the finite element method:

..

= fext

-

fint (4.1)

Equation (4-1) will be divided in two parts: a and b. Part a will be integrated with the tirne step At and part b with the time step At/E* The subset a of the set o f equations (4.1) is considered at time t, and the subset b of the equations (4.1) is considered at time tn+p where p is the

subcycle number (1 <= p <=

ml:

m

( 4 . 3 )

The central difference forinulas are:

and :

where:

(4.6)

4.1 Internolation

Let the displacements of the boundary nodes

xB

-

boundary nodes are nodes which are a member of partition a, but belong to an element which has nodes belonging to partition b toor e.g. node p in fig. 3.2

-

at the times Between tn and tn+l be calculated using a linear interpolation:

The algorithm to obtain the solution o f (4.1) at the time interval [ O F T ? is given below. (i) begin 1ii) n = O (iii) (iv) p = 1 tv) calculate

:

:

+

d

with ( 4 . 7 ) (vi) (vii) (viii) if p

<

m+1 go to (v) (ix) (XI end

calculate

::+,

with (4.2) and 14.4)

m

calculate

xb

D with (4.3) and (4.5)

p = p

+

1

if nat

<

T then n = n+l and go to (iii)

-n+,

4.2 Intearation

If integration is applied, the displacements of the boundary points at times between t, and tn,l are calculated using:

The algorithm to calculate the solution o f (4.1) at the time interval [O,T] is the same as the algorithm given in section 4.iF except that in step ( v i the displacements o f the boundary points are calculated with (4.8) instead

5 STABILITY ANALYSIS

In this chapter the stability of the alternatives described in the previous chapters will be discussed. Because a general stability analysis cannot be given, the stability is analysed using a simple linear systeli (Eig. 5.11. The method used for the stability analysis is the same as the method used by Belytschko et al. (1979).

Uk m k m

Fig. 5.1: The system used for the stability analysis.

The equations o f motion are given by:

where:

Equation (5.1) can be divided in two equations:

mX1

+

(ui-1) kx’

-

kx2 = O

mi2

+

kx2

-

kx 1 = O

Equation ( 5 . 2 ) is integrated using the time step At12:

and equation ( 5 . 3 ) is integrated with the time step At:

In section 5.1 it is shown that, if straight explicit integration is applied, the method used to determine the stability of the subcycling methods leads to a well known expression for the maximum stable time step. In the sections 5,2 and 5 . 3 , the stability of the two alternative treatments of the boundary nodes will be given.

5.1 Straight ExDlicit Intearation

If straight explicit integration is applied, the equations 15.2) and ( 5 , 3 )

are both integrated with a time step At+ Substitution of (5.4) f o r both degrees of freedom yields:

Assume that the solutions are of the form:

1 2n

x, = A 1 h

xn = -92A 2n

( 5 - 8 )

( 5 . 9 )

In Belytschko et al. (1979) the same assumptions for the solutions are used without further references. More literature research has to be done to prove that the solutions can be written in the above form.

Substitution of (5,8) and ( 5 . 9 ) in (5.6) and ( 5 . 7 ) yields:

with:

(5.11) If the solution x is ctableothe roots o f the above determinant must lie within the unit circle. The determinant is given by:

-

with:

Ecpation (5.12) can be transformed By using the z-transformation:

(5.15) 2 - 1 9 2

- m

The requirement that ! A ! <= 1 is then replaced by the requirement that the real part of z must be smaller than or equal t o zero. The transformed equation is:

Application o f the Routh-Hurwitz criterion (appendix A) yields:

or :

#e are interested in the smallest time step which satisfies (5.18):

The maximu.m eigenfrequency of the system is given by:

Combination of (5.191 and (5.20) yields the well-known expression (Bathe &

idilson. 1976) for the maximuni stable time step o f the explicit method: 2

At <=

-

"max (5.21)

5.2 Stabilitv of the Intermiation Hethod The computation scheme for each cycle is:

1. Consider ( 5 . 3 ) at time t, and substitute (5.5) to calculate x2 ~ + ~ ~

2. Use a linear interpolation between

x,

2 and xaatl. 2

3 . Consider (5.2) at time t, and use (5.4) t o calculate xn-+$.

1

to calculate xn+. 2

4. Consider (5.2) a t tiliie t 1 an expression similar to ( 5 . 0 ) to calculate n+Z

"2 n -

CI

1 1 2 (5.23) x n + i 1

-

2xn 1+ x

"-2

1 1

+

?(a 1

+

llwx,

-

p x n = O

-

2x 1 1

+

xn 1+ 1?(a + l)wx,+~ 1

-

2 + x i ) = O (5.24)

'ntí

n+Z

where w = :At2. Substitution of a solution of the form: x, 1 = A I P xn 2 = A2h2n XA+$ = A3h 2n+l yields:

- p h

1 .L

The characteristic equation is

ai th:

(5.25) (5.26)

h8

+

h1h6

+

h2A4

+

hlh2 i 1 = O _(5.29)

Application of the z-transform yields:

The use of the Routh-Hurwitz criterion yields:

(5.329

15-31]

This is a third degree polynomial in w e Ye are again interested in the

smallest value o f bt (or vi) satisfying (5.331, This value is calculated

numerically and plotted as a function of a in fig. 5.2. The results of the stability analysis of the straight explicit method. are also given in this

1.60 1.40

-

i 1.20 Y :: 1.00 P %

-

, 0.80 0.60 - - - - - - I STABILITY ANALYSIS

=.

Q.. 'y,

I

0.40

2

MAX. TIME STEPS I

Fig. 5.2: The rriaximum stable time step as a function o f alpha. 5.3 Stabilitv of the Intearation I-ethod

The computation scheme f o r each cycle is:

1, Consider (5.3) at time t, and use (5-51 t o calculate x2 ~ + ~ ~

2. Consider (5.2) at tine tn and use (5.4) t o calculate

XI

1.

3 . Consider (5.3) at time t, and use an equation similar to (5.4) to

4 . Consider ( 5 . 2 ) a t time t 1 and use an equation similar t o ( 5 . 4 ) and the n+Z

2 calculate

n+2 1

result of step 3 to calculate x ~ + ~ .

Each time step the following computations have to be performed: - 2xn 2

+

xn-l 2

-

WX; + W X ~ = O %+l x1 1 - 2xn 1

+

x 1-1 t g(a 1

+

l)wxn 1

- p x n

1 2 = O x2 1

-

ax, 2

+

x 1 2

+

;wx;

-

&xi

xn+l 1 - 2x 1 1 t xn 1

+

q(a+l)wx 1 i 1

-

1wx2 i

"+z

n Z

"+z

n-z

n+Z n+Z

P

n+Z (5.34) (5.35) (5.361 (5*37)

k 2

vhere w =

.

Substitute a solution of the form: X I n = Al$n x1 1 = A3h ".+i = A q h x2 _n = 2 n f l ntZ 2 2 n + l This yields: ( 5 . 3 8 ) (5,391 (5.40) (5.41)

The use of the a-transform yields:

The use o f the Routh-Hurwitz criterion yields:

This i s again a third degree polynomial in w. We are again interested in the

smallest value for At satisfying ( 5 . 4 7 ) . This value is calculated numerically for different values of ct and plotted i a i fig. 5.2.

6 IEPACT OF AH ELASTIC BAR

In this chapter the two alternatives are employed to simulate the inpact of

an elastic bar. The alternatives are compared with respect to their

numerical stability, accuracy and computing time consumption. The modeled bar consists o f 19 elements, i-e. 10 small and 9 large elements (fig. 6 * 1 j e

1

4+-\

_P P

10

small e

1

emen t s 9 large elements

Fig. 6.1 The mesh used for the impact simulation.

The sound speed c in this material is equal to 5.00 10 3 [m s-'] I so the

maximm stable time steg in the small elements is:

( 6 - 1 ) 5tma,(sma11 elements)= l = 2.00 is]

and. in the large elements:

M.1

nodes that belong t o a small element (1 up to and including 10) are integrated with the time step Atlrn and nodes that belong to large elements only (11 up t o and including 19) are integrated with the time step Atcit. Because it is not yet possible to predict the inaxiirmm stable time step o f

this problem numerically. this time step is determined by numerical

experiments. The results are plotted as a function o f m in figure 6.2. Por further calculations, the used time step is taken to be 80 percent o f the maximum stable time step, These time steps are given in fig. 6 . 2 also.

In fig. 6 . 3 , the computing time of the two alternatives are plotted as a function of m. For each situation! the same time interval is calculated with the time steps given in Eig. 6.2. Application o f subcycling reduces the necessary computing tiïne with a factor three for this example. The shape o f

IMPACT ELASTIC BAR

I

TIME STEP O . O O ~ . " . ' . , , . I . . . I I . l l . l . . . . , . . . . , . , , . I used tstep 0.0 2.0 ¶.O 6.0 8.0 10.0 m 1-1

I-

interpol. v---+---+ integration Q. ... 0 ... *

Fig. 6 - 2 : The time step as a function o f m.

35.0 30.0 u €i 8 20.0 15.0 10.0

IMPACT ELASTIC BAR

I . . . . , . . . . , , . . . , , . . . , . , , . , . . . . , 0.0 2.0 ¶.O 6.0 8.0 10.0 12.0 14.0 m [-I COMPUTING TIME interpol.

-

_integration r---+---+

used for t h e calculations increases. This roeans that the total number of cycles decreases and the computing tirite decreases, For 9.larger than 10, the

time step used remains the saxte. This means that the number of cycles will remain the salne. Hosever, for increasing roc the nui%ber oi subcycles per cycle will increase. This means that the computing time will increase again.

In the figures 6 - 4 up to 6 . 6 ! the accuracy of the methods are compared. In fig. 6 . 1 the calculated solution is compared to the 'exact' solution.

Actually, in this case the name 'exact' solution is not correct, because an exiict solution is not available, The 'so-called 'exact' solution is the solrition which is calculated by straight explicit integration with a very small time step, i.e. 1 lo-'' Es]. The value which is plotted as a function of xt in fig, 6.4! i s defined by:

max. relative error

In (6.31 xmax is the maximum t, and i for the node numbers.

appearing displacement? E stands for the t i m s

in fig. 6.5. the ritaxiinum relative force error is given as a function o f the time The ì ~ a x i m u ~ ~ relative force error ferror is defined- as:

..

4- f

'

,int' fer,,, = max( f m x I f aa ,int (6.51

where _{,int i} is the n o m of the i~aximus appearing internal force.

in figure 6.QI the maximuro relative energy deviation Eer,,, ia plotted as a function of m. Eerror is defined by:

where Enot is the total energy at time tn and tirrie zero.

is the total energy at

In closing, in fig. 6.7, the computing time of the two alternatives are given as a function of lii? but with equal maximm relative errors. This means

that the time step of %he interpolation method is decreased until the value

of the aiixiaurn relative error was equal to the maxim.urn relatise error of the integration method.

3.50 *io2 3.00 ?. I I E 2.50 e 2 2.00 m 3 X 1.50 1.00

IMPACT ELASTIC BAR r 0 . 0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 m [-I M A X . REL. ERROR interpol.

-

integration t - - - + - - - +

Fig. 6,4: The maximum relative e r r o r as a function of m.

3.00 .lo1 2.50 U u 1.50 3 al 1.00 0.50 0.00

IMPACT ELASTIC BAR

0.0 2.0 4.0 6 . 0 8.0 10.0 12.0 14.0

m [-I

REL. FORCE ERROR

: 4

3.0C .lo2 2.50

-

i 2.00 2 1.50 3 2 1.00 0.50 0.00

IMPACT ELASTIC BAR

I . . , . , . . . . , . . . . , . < . . , . . . . , . , . . ,

0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 m [-I

REL. ENERGY DEV.

Fig. 6.6: The maximum relative energy deviation as a function of m.

IMPACT ELASTIC BAR

0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 m [-I COMPUTING TIME i n t e r p o l .

-

i n t e g r a t i o n t - - - + - - - +

Fig. 6 , 7 : The computing time as a function of m aitol the same maximum

7 IEPACT OF A PLASTIC BAR

The two alternatives are employed to simulate the impact of a plastic bar. The goal is to compare the alternatives with respect to numerical stability, accuracy and computing time consumption. when they are applied to a

nonlinear set of differential equations. The mesh used for the simulation is the same as the rnech used €or the simulation of the impact of an elastic bar

(fig. 6.1). The material properties and some other data are given in fig. 7.1. t a m l t a m 2 p B.. Q = _{E = 6.89 lolo EN m-2]} = 5.38 lo8 [N m-2i = 2.76 IO8 [M m-2] = 2,76 lo3 hkg i-31

a

= 314.16 IO6 Em2] "0 isotropic hardening = 50 Cm s-'1

Fig. 7.1: Data for the plastic impact simulation.

The nodes 1 up to and including 10 are again Integrated with the time step &t/m and the nodes 11 up t o and including 19 are integrated with the time step At,

In fig. 7.2 the maximum stable time step is given as a function of m - These time step are determined by numerical experiments. The maximum stable time step of the integration method decreases if m is larger than 6. No results for the integration method for m larger than 8 are presented, because of this decreasing maximum stable time step, The used time step for the

simulations is again taken to be 80 percent of the maxiinurn stable time step and also plotted in fig. 7.2.

Fig. 7.3 shows the computing time as a function of m.

In figse 7.4! 7.5! 7.6 and 7.7 respectively the maximum relative error of both methodsl the maximum relative error o f the interpolation method, the maximum relative force error and the maximum relative energy deviation are given as a function of m. These quantities are defined by the equations

i n t e r p o l . p_ used i n t e r p . v---+---+ i n t e g r a t i o n 0. ..--... 0- _...__ 0 used i n t e g r . t. ..-+....-+ IMPACT PLASTIC BAR

2.00 *lo6 1.50 10

-

I 4 : 1.00 a Y 0.50 0.00 - . T FT---* ----_- - - - . . ' . I . .8 . 1 . . . . I , . . . , . . , , , _{I . . . . ,}

Fig. 7.2: The time step as a function of an.

IMPACT PLASTIC BAR

5 5 . 0 I 0.0 2.0 4 . 0 6.0 8.0 10.0 12.0 14.0 m [-I - TIME STEP COMPUTING TIME i n t e r p o l .

-

i n t e g r a t i o n ' - - - q . - - - a p

9.00 *lo3 8.00 7.00 .-. I I LI 2 6.00 3 li 1 5.00 E 4.00 3.00 2 . 0 0

IMPACT PLASTIC BAñ

-

: - - - - - - " 1 . ' . , . . I . . . I I . . , . I . . . . I . . . . , . . . . ,

MAX. REL. EPAOR

~ 4.0C *lö2 3.00

-

I LI m I 2 2.00 X E 1 . 0 0 0 . 0 0 1 . " ' I ' " ' " ' . * 1 . , ' . I . . * . ' . . . . ' . * . . I 0.0 2.0 4 . 0 6.0 8 . 0 1 0 . 0 12.0 1 4 . 0 m I-] i n t e r p o l .

-

_{i n t e g r a t i o n} p---+---+

Fig. 7.4: The maximum relative error as a function of m.

MAX. REL. ERROR

Pig. 7.5: The maximum relative error of the interpolation method as a

1 . 2 0 *lo1 1.00 U il 0.60 3 Iu 0.40 0 . 2 0 0.00 IMPACT P L A S T I C BAR 0 . 0 2.0 4 . 0 6 . 0 8 . 0 1 0 . 0 1 2 . 0 1 4 . 0 m 1-1

REL. FORCE ERROR

Fig. 7.6: The maximum relative force error as a function of M.

0.00 IMPACT P L A S T I C BAR 2.50

-

i 2 .0 0

i

1 . 5 0 3 8 1 . 0 0 0 . 5 0 0 . 0 2 . 0 4 . 0 6.0 8 . 0 10.0 12.0 14.0 m 1-1

REL. ENERGY OEV.

8 DISCUSSIOW

If the alternatives are applied to a linear problem - the impact of an elastic bar

-

the performance o f both methods is the same. The computing time as a function of m, shows an optimum for m = 8 . Because the ratio between the length of the large elements and the length of the small elements is equal to 10, it is striking that the optimum is not found for m = 10. For the impact of the elastic bar, the computing time reduction is about a factor three. The maximum relative error increases with a factor three, but, for m = 8 its value is about three percent! which is acceptable. If the alternatives are applied to a nonlinear problem

-

the impact of a plastic bar

-

the performance of the interpolation method is much better

than the performance of the integration method. For the latter, the rnaximum stable time step increases until m = 6 , but decreases for m larger than 6 . Por this method, the optimum of the computing appears for m = 4. The maximum stable time step for the interpolation method shows the same behaviour as in the linear case, except from the fact that the optimum in the computing time appears for m = 10 instead o f m = 8 in the linear case. This can be declared from the fact that in the linear case, the maximum stable time step

increases until m = 8 , and in the nonlinear case until m = 18. For the

integration method the computing time reduction is about a factor three, for the integration method the reduction is about a factor two. If the accuracy of the methods are concerned, the performance of the interpolation method is better too,

For future research! the following things are proposed:

-

It would be interesting to find out why the maximum stable time step

decreases if m increases for the integration method applied to a nonlinear problem.

-

It is necessary to develop a nore general stability criterion for subcycling methods.

-

The subcycling methods should be implemented in a finite element package in order to be able to examine their behaviour for more dimensional problems and for a larger number of degrees of freedom.

-

Hore alternatives should be developed, and compared to the performance of the alternatives in this report, using the simulation of the impact of the elastic and the plastic bar.

Bathe, K.J. & Wilson, E.L. (1976).

Numerical Methods in Finite Element Analvsis. Englewood Cliffs, Mew Yersey: Prentice-Hall, Inc.

Belytschko, T., Yen, H.-J. & Hullen, R. (1979).

Mixed Hethods for Time Integration. Comrmter Hethods in A m l i e d Mechanics and Enaineering, 17/18, 259-275.

Hinorsky, N. (1962).

The theoren of Routh-Hurwitz, which is groved in -Iinorsl;yF ( 1 9 8 2 1 , states that all roots of the polynomial:

tihere ai is real (i = O,.e.pnl, have negative real parts or a real part equal to zero, if:

a0 >= O ( A . 2 )

and all the leading principal ainors o f 0 are not negative, with:

a. a2 a 4 -

- -

-

-

O al a 3 -

- - - -

I]

al a3 a s -

- - - -

O a. a s -

- - - -

I O

1

;

o 0 0 I I i i ! t ! 1

‘J

o

( A . 3 )

For the situations considered in chapter 5, a - is equal to zero if j is odd. This means that all principal 1Einob-s of 0 are zero. The only remaining demand is a0 >= O,