PARALLEL PROCESSORS AND

NONLINEAR STRUCTURAL DYNAMICS ALGORITHMS AND SOFTWARE

Principal Investigator: Ted Belytschko

Department of Civil Engineering

Northwestern University

Evanston, Illinois 60208

Semiannual Status Report

May 1, 1987 through October 31, 1987

NASA Research Grant NAG-1-650, Supplement No. 2

£*€&-1-829 25) OBA1LEL FECCESSOBS IHD 1*88-24167

NCNLIKEA E STEOClOFflL DXHABICS AiGCBITUUS ASD

£Cfl»*B£ Seffliancual Status Eepc-xt, 1 Hay -

31 Oct. 1S67 (Scirthwestern Ciiv.) 48 p Unclas

CSC1 09B G3/61 01U642 9
PREFACE

This research was conducted under the direction of Professor Ted

Belytschko. Participating research assistants were Noreen 0. Gilbertsen and

Mark 0. Neal. The help of Argonne National Laboratory, particularly Or. James

M. Kennedy who provided access to several parallel computing machines, is also

appreciated.

The following papers supported by NASA were submitted for publication:

Mark 0. Neal and Ted Belytschko, "Explicit-Explicit Subcycling with

Non-Integer Time Step Ratios for Structural Dynamic Systems", Computers

and Structures, submitted.

P. Smolinski, T. Belytschko, and M. Neal, "Multi-Time-Step

Integration Using Nodal Partitioning", International Journal for

Numerical Methods in Engineering, 26, 349-359, 1988.
Abstract

An explicit-explicit subcycling procedure for the finite element analysis of structural

dynamics is developed. This procedure has relaxed the usual constraint of requiring integer

time step ratios for adjacent nodal groups. This allows for greater advantage to be taken of

local stability criteria, and thus improves the efficiency of the explicit time integrator.

Example problems are included to demonstrate the accuracy and stability of the method.
I. Introduction

Explicit-explicit subcycling is oftentimes used in dynamic finite element problems

in order to improve computational efficiency when the element size varies over the mesh.

This is accomplished by separating the elements or nodes into groups and assigning a

different time step to each group. The time step of each group depends only on the

frequencies of the elements in that group. This results in great computational savings since

less stiff elements can be assigned to the same group. In this way a time step is used for

each group which is much closer to the critical time steps of the elements in the group.

Mixed time integration was first introduced by Belytschko and Mullen[l] where they

presented an "implicit-explicit" method for the second order equations which arise from the

semidiscretization of the continuum equation. This was a nodal partition where both

implicit and explicit methods were used in the time integration process. The stability of this

type of procedure was studied in [2,3]- Hughes and Liu[4,5] later introduced an "implicitexplicit" method which divided the mesh into groups of elements rather than nodes, known

as an element partition. Both methods were found to substantially improve computational

efficiency. An "explicit-explicit" subcycling procedure was then introduced by Belytschko

and Mullen[6]. In this procedure each group of nodes was integrated explicitly, but

different time steps were used for each group.

Mixed time integration procedures for first order semidiscretizations have also been

proposed. Belytschko et al.[7] present sufficient conditions for stability of their element

partition by using the element eigenvalue inequality directly on the evolution equations.

The stability analysis of nodal partitions is complicated by the fact that their amplification

matrices are unsymmetric. Smolinski et al.[8] avoid solving any unsymmetric systems in

their nodal partition, and are able to develop stability conditions for semidiscretizations of

the diffusion equation.
A major drawback to all these previous methods is that integer time step ratios are

required. That is, if one group with a time step of At shares nodes or elements with a

second group, then the only possibilities for the time step of the second group are [(n)At]

or [At / n] where n is a positive integer. This introduces two difficulties when dividing the

mesh into groups. First, it can severely restrict the division of the mesh into more than just

a few groups, and therefore decrease the efficiency of the mixed time integration procedure

being used. Second, since the time step is restricted to just a few possibilities, a time step

may have to be used for a group that is much smaller than the stable time steps of the

elements in the group.

Recently, Mizukami[9] has relaxed this restriction for first order ordinary

differential equations which result from semidiscretization of the diffusion equation. He

presents an explicit-explicit nodal partition that allows for much greater flexibility when

assigning time steps to the different nodal groups. He then gives a stability proof of the

method, however, this proof only applies to a particular class of elements. Donea and

Laval[10] have recently extended this proof of stability to arbitrary elements.

The purpose of this paper is to present an explicit-explicit subcycling procedure 'for

structural dynamics which does not require integer time step ratios for adjacent groups. It

will be shown that this procedure requires very little additional storage and appears to be

stable for a wide variety of example problems.
n. Nomenclature and definitions

n = time step iteration counter

At M = master time step

Atj = nodal time step

At £ld = nodal time step for the previous step

^

t

 critical _ critical time step of an element

At e = time step used to update an element

^ MLS

 = master level time step

tM = master clock

tj = nodal clock

te = element clock

CMLT = master level time

d = global vector of displacements

v = global vector of velocities

a = global vector of accelerations

d

e

 = element level vector of displacements

v

e

 = element level vector of velocities

CFe = Cauchy stress tensor for an element

CTe = stress rate tensor for an element

E

e

 = strain tensor for an element

£

e

 = strain rate tensor for an element

r>

f = global vector of nodal forces

f exi = global vector of external nodal forces

fint = global vector of internal nodal forces

f fnl = element level vector of internal nodal forces
Wext

B

L

e

M

mi

Q

e

NUMEL

NUMEL

MOD(i j)

= external work done on the system

= kinetic energy of the system

= internal energy of the system

= element level increment in internal energy

= gradient operator

= element assembly operator

= lumped mass matrix

= diagonal component of the mass matrix for node I

= domain of the element

= total number of elements

= total number of nodes

= remainder of the division of i by j
HI. Governing equations

The equations of motion for structural dynamics are given by

Ma + f = 0 (3.1)

where,

M = mass matrix

a = global vector of accelerations

f = global vector of nodal forces.

The nodal forces are comprised of external and internal parts.

f = fint-fext (3-2)

where,

fext = global vector of external nodal forces

f int = global vector of internal nodal forces.

The internal nodal forces for a single element are calculated from the element stresses and

can be written

f

e I T>T *m c A f\ /o o\

int = ne

 B

 «

 d" (3.3)
where,

f

1

?

int

B

<7e

= element internal force vector

= domain of the element

= gradient operator

= Cauchy stress tensor for an element

These internal forces are then assembled to the global vector

f;mt

(3.4)

where,

L

e

 = element assembly operator.

For nonlinear analysis a strain rate formulation will be used to calculate the stress of

an element.

= CEe(3.5)

(3.6)

where,

strain rate tensor for an element, such as the velocity strain or rateof-deformation

frame-invarient stress rate tensor for an element
C = constitutive matrix

At e = time step used for the stress update of an element.

The strain rate tensor is given by the strain displacement relation of the element as follows:

E

e

 = Bve (3.7)

where v

e

 is the element velocity vector.

In this paper, central difference time integration with variable time steps will be

used to to advance the nodal velocities and displacements in time. The equations are given

by the following:

v L n . i / 2 (At I + At°ld)aLn - (3.8)

d

l.n + At I vI,n+l/2

 (3'9)

where,

aLn = acceleration of node I at time (n)At

vI>n+i/2 = velocity of node I at time (n+V2)At

^i.n-t-l

 = displacement of node I at time (n+l)At

At j = time step of node I

At £ld = time step of node I for the previous step.

Since the time step for a node could change for each step, it is necessary to average

the time steps as shown in equation (3.8) for the velocity update. Note that equation (3.6)

has been written so that different time steps may be used for different elements. This
modification of the usual equations has also been made on equations (3.8) and (3.9) so that

different time steps can be used for each node. The use of different time steps for different

nodes and for different elements is common to all subcycling procedures. For the

procedure to be presented here, however, the element update time steps of equation (3.6)

are not necessarily the same as the nodal update time steps of equations (3.8) and (3.9).

This added feature is unique to this new procedure and allows for the use of non integer

time step ratios of adjacent nodes as will be explained in the following sections.
IV. Subcycling Procedure

The basic idea underlying this subcycling procedure is the use of a master clock

which is incremented by the smallest time step which occurs within the mesh; a master

clock was also used by Belytschko and Liu[l 1]. This smallest time step will be referred to

as the master time step. While the master clock is updated by the master time step, the

nodal clocks are incremented by the time step of each node. Whenever a nodal clock time

is equal to the master clock the velocity and displacement of the node are updated with

equations (3.8) and (3.9) using the time step of the node. The time step of each node must

be an integer multiple of the master time step and cannot be greater than the smallest critical

time step of the elements surrounding the node. Although the nodal time steps must be

integer multiples of the master time step, the attractive feature of this procedure is that it is

not necessary for the time steps of two adjacent nodes to form an integer ratio.

The times at which each element is updated are also governed by the nodal clocks of

the nodes of the element Each element is updated every time one of the nodes of the

element will undergo a velocity and displacement update. For the purpose of this element

update, it is necessary to use an element clock for each element. This element clock is set

equal to the master clock each time the element is updated. In this way a record is kept of

when each element was last updated. The time step for the incremental update of the

element stresses, equation (3.7), is then given by the difference between the master clock

and the element clock. It is important to realize that this element update time step does not

necessarily coincide with the time step of any of the element's nodes.

As can be seen, the structure of this procedure is very similar to the subcycling

algorithm described by Belytschko and Liu[l 1] and by Belytschko[12]. The major

difference is that in the previous algorithm the times at which an element was updated was

governed by element clocks while in this new procedure they are governed by the nodal

clocks of the nodes of the element Figure 1 gives the nodal influence diagram for the

10
partition described in references [11] and [12]. In this figure an arrow represents the

influence of one node on the velocity and displacement calculation of another node, a solid

dot represents a time at which the velocity and displacement of a node are updated, and an

open dot represents a time at which the displacement of a node would have to be

interpolated between two known values. In order for a partition to work correctly the

element stresses must be calculated whenever a node of the element is updated. In figure 1

it can be seen that the stresses of element 1 are updated every 6At, while the stresses of

elements 2 and 3 are updated every 2At as is shown by the bold horizontal lines. Because

the time steps of these two nodes form an integer ratio, element stresses are calculated

whenever they are needed for nodal updates. However, if this procedure is used in the

case of non integer time step ratios as shown in figure 2, then the appropriate element

stresses would not be available to update each node's acceleration, velocity, and

displacement. For example, to update the nodal values of node 2 at time t = (n+3) At, it is

necessary to have the element stresses at this time for all elements that are adjacent to node

2. From figure 2 it is seen that the stress values of element 2 are known at times (n) At,

(n+2) At, and (n+4) At, but not at the required time, (n+3) At If nodal clocks are used

instead, as shown in figure 3, the values of the element stresses are calculated whenever

they are needed for a nodal update so that even in the case of non integer time step ratios it

is possible to correctly update the nodal values.

11
V. Implementation of the Subcycling Procedure

A. Nodal Grouping

The first step in the nodal grouping procedure is to assign a time step to

each node. This is accomplished as follows:

1. Assign a very large time step to each node, At j.

2. Loop on the elements, (JE).

Calculate the critical time step of the element, At£ritical(JE).

Loop on the nodes of this element, (I).

At j = Minimum(At x, At |ritical(JE)). ' (5.1)

In this way each node is given the time step of the smallest element time step

of the elements connected to the node. Next, it is necessary to determine a master

time step so that all nodal time steps are an integer multiple of the master time step.

This is most easily accomplished by letting the master time step be given by the

smallest nodal time step.

AtM = Minimum (At z). (5.2)

The remaining nodal time steps are then reduced by the equation

Atj = AttInteger(Atj/At)]. (5.3)

12
By equation (5.3) it is obvious that the master time step must be chosen so

that it is smaller than or equal to all nodal time steps. It can also be seen from

equation (5.3) that all nodal time steps are integer multiples of the master time step.

Although equation (5.2) was used to determine the master time step, it should be

noted that the master time step can be set to any value as long as the following two

restrictions are satisfied.

1. All nodal time steps are integer multiples of the master time step.

2. A nodal time step is not greater than the smallest critical time step of the

elements surrounding the node.

Next, the master level time step must be calculated. This identifies the time

increment that occurs between steps when all nodes are updated and is defined,

therefore, as the least common multiple of all nodal time steps.

^ MLS

 = Least common multiple (At l). (5.4)

Once this value is calculated the master level time can be calculated by

adding the master level time step to the current master time.

'MLT = 'M+^MLT ' (5-

5

>

This master level time identifies the next time that all nodal and element

clocks will be the same value, that is t j^y. This means that all nodes and element

stresses will have been updated to this master time level. It is useful to have all

nodal values and element stresses updated to the same time level at some point in

13
the analysis so that the energy balance of the system may be checkedOnce the

master level time is determined the nodal grouping is complete.

To illustrate the nodal grouping procedure, an example is given in figure 4,

along with the element time steps, nodal time steps, master time step, master level

time step, and master level time.

B. Nodal Force Calculations

After the master time step, nodal time steps, and the master level time have

all been assigned as described in the previous section, the time stepping procedure

can begin. The first step in the procedure is to calculate the nodal forces. The

external nodal forces at each node are calculated as usual, however the calculation

of internal nodal forces for each element must be modified. Before an element is

updated the clocks of all the nodes of the element are checked by the condition

Is (tl = tM) for node I (5.6)

where,

= nodal clock

tM = master clock.

If condition (5.6) is true for any node of a given element, then the nodal

clock is equal to the master clock and the node will be updated during this step.

This nodal update will require the correct element stresses at this time for any

element connected to this node. Therefore, this element must be updated during

this step. This can be seen more easily with the help of an example. If it is

14
assumed that condition (5.6) is true for node I of an element, JE, then the velocity

and displacement of node I will be updated by equations (3.8) and (3.9) which

depend on the nodal acceleration, ar The calculation of al requires current internal

nodal force values at node I, so that any element connected to node I must undergo

a stress update at this time, so that the element stress coincides with the current

time, t M. Therefore, if condition (5.6) is true for any node of an element, the

stress of the element must be updated. For the purpose of this stress update an

element time step must be calculated. This element time step is given by the

equation

At. = t M - t e (5.7)

where,

.•

Ate = element time step

te = element clock.

This will give the time increment that has occurred since the last stress

update of this element. Note that this element time step is not necessarily the same

for an element each time it is updated. This can be seen by examining element 2 of

figure 3. Once the element time step is calculated, the stress update of the element

can take place. This update is given by equations (3.5) through (3.7) with the

element update time step just calculated, At e. In this way the element is updated to

the current master clock time, t M, using the appropriate time increment for each

element. After the stresses have been updated the nodal forces are calculated in the

usual manner by equations (3.3) and (3.4). The final step in the internal nodal

15
force calculation is to increment the clocks of all elements that have been updated

with the equation

te = tM. (5.8)

C. Nodal Updates

The next step in the time stepping procedure is to update the nodal variables.

Again, central difference time integration is used, however, the procedure given by

equations (3.1), (3.8) and (3.9) must be modified slightly. The calculation of the

nodal acceleration at time (n)At is the same as before and is given by

(5-9)

mi

for node I. To update the velocity of node I at time (n)At, the validity of the

following condition is checked:

Is (tj = tM) fornodel? (5.10)

If condition (5.10) is true for node I, then equations (3.8) and (3.9) are

used to update the velocity and displacement of node I. If condition (5.10) is false,

then the nodal values are not changed. Once the new velocities and displacements

have been calculated, the nodal clock is incremented by the nodal time step.

(5.11)

16
VI. Energy Calculations

A major drawback to explicit time integration schemes is the fact that they are

conditionally stable. When performing explicit calculations a limit must be placed on the

time increment being used in order to prevent a loss of stability. In linear calculations this

presents no problem since in this case the stability criteria remain the same for the duration

of the analysis. However, in nonlinear calculations, the stability conditions of the mesh are

constantly changing due to changes in problem geometry or material properties. Therefore,

some means of checking the stability of the problem must be introduced. In linear

problems, when an instability occurs it is characterized by exponential oscillatory growth in

the nodal displacements. A periodic check on displacements can usually guarantee that no

instability has occurred in the problem. However, in nonlinear situations where plastic

deformations may dissipate large amounts of energy, the decreased modulus may enable

the system to regain stability (known as an arrested instability, Belytschko[13]). In these

situations simply checking for large unrealistic results may not always provide assurance

that no instability has occurred, since an instability may be confined to a small region of the

mesh and the problem may even regain stability at a later time. If this happens, the results

obtained for the problem may seem quite reasonable, but the earlier instability in the

problem will render the final results totally inaccurate. An alternative method for checking

the stability of explicit finite element calculations is the use of an energy balance check.

Hughes and Belytschko[14] gives this energy balance check in the form

I W^+W^-Wr I <8llW H (6.1)

where,

WnXt = external work performed on the system up to time (n)At

17
W£" = internal energy of the system at time (n)At

Wn^ = kinetic energy of the system at time (n)At

8 = tolerance

and "W" is some measure of the total energy of the system. The use of any one of the

terms on the left hand side of equation (6.1) alone would not provide a useful measure of

this energy. For impulsively loaded problems when an initial velocity is prescribed, there

would be no external work. For vibration problems the values of the kinetic energy and

internal energy would oscillate and become very small at some point in the problem.

Therefore, a measure of the total energy of the system should include all three terms.

(6.2)

The use of absolute value sign in the left hand side of equation (6.1) is not really

necessary since an instability is usually characterized by a growth in kinetic and internal

energies which are nonnegative except for roundoff errors; therefore, this left hand side

would always be positive in the case of an instability.

In time integration without subcycling equation (6.1) is checked at each time step to

insure that the work done on the system and the energy of the system are nearly equal.

Since a numerical instability in a problem always implies a loss of energy balance, it can be

assumed that if equation (6.1) is true at every step, there is no such instability in the

problem. For the case with subcycling, this energy balance cannot be checked at every

step. This is because the internal energy of each element is not updated at every step. This

can be seen in the equation for the increment of internal energy for an element.

e e T

18

)

 e ^ (6 3)
where,

= increment of internal energy for an element from

time (n-m) At M to (n) At M.

In equation (6.3) the subscript m represents the number of master time steps that

have passed since this element was last updated. After this increment of internal energy for

an element is calculated, it is added to the total internal energy of the system.

(6.4)

Note that for the subcycling procedure, this element internal energy increment is not

calculated for every element at every step. It is only calculated for an element when the

stress of the element has been updated. At other steps the internal energy in the element is

not known. Therefore, the energy balance can only be checked at steps when the stress of

every element is updated. At these steps the contributions of all elements are added to the

internal energy of the system and the correct value of the internal energy of the system is

known. By definition these steps occur only at the master level times described earlier.

The external work done on the system is calculated at every master time step and is

given by the equation

(6>5)

19
The kinetic energy can only be computed at half steps since the nodal velocities are

only calculated at half steps.

(6'5)

The kinetic energy at each step is then found by averaging the half step values.

W

kin

n

20
. Additional Storage Requirements

The addition of subcycling to an existing finite element code requires only four new

airays of storage. The additional storage requirements are as follows:

t e = a clock for each element

tj = a clock for each node

Atj = nodal time step for each node

At °ld = nodal time step for each node for the previous step

At MLS = master level time step.

The necessity of these additional storage requirements and a general overview of the

procedure can be understood with the help of the flowchart given in figure 5.

21
. Numerical Examples

To test the stability and accuracy of the proposed subcycling scheme, two example

problems have been examined. Problem 1 consists of a bar subjected to a suddenly applied

load at one end while the opposite end is fixed to a rigid wall. The bar is modelled with 32

constant strain rod elements which have been divided into four groups. The problem

dimensions and material properties are shown in figure 6, while the nodal grouping data is

given in table 1. In case 1 the external force applied to the bar is small so that the problem

remains completely elastic, while in case 2 plastic deformation takes place due to the larger

external force. The elastic problem is solved with the proposed subcycling procedure.

Results of the elastic problem are given in figures 7 and 8. Figure 7 gives the time history

of the stress at x=5.5 for both the subcycled and non subcycled cases. The exact solution

is also included for reference. Figure 8 gives the energy error for both the subcycled and

non subcycled cases. Figures 9 and 10 give similar results for the plastic problem. These

results indicate that the subcycling procedure performs well for the proposed problem,

since the stresses are similar with and without subcycling, and there is minimal energy

balance loss for the subcycled case.

The second example problem is a rotating ring subjected to a suddenly applied

internal pressure. The problem dimensions and material properties are given in figure 1 1,

and the finite element meshes are given in figures 12 and 14, Key[15]. Figure 13. gives

the time history of the maximum stress of a single element for the mesh shown in figure

12. Figure 15 compares the effective plastic strain of an element as a function of time with

and without subcycling. Figure 16 gives the effective plastic strains of two different

elements in the subcycled case, while the energy error is shown in figure 17. Figures 15

through 17 all refer to the 100 element mesh of figure 14 with the nodal grouping data in

table 2. From the figures it can be seen that the subcycled solution retains the symmetry of

the problem and appears to be stable since there is no appreciable loss of energy balance.

22
Vin. Conclusions

The subcycling, or mixed time integration, procedure described here has great

flexibility in exploiting the increased stability of domains with larger of more flexible

elements. Thus the difficulties associated with the presence of a few stiff elements are

avoided. The procedure is easy to program, requiring on the order of 40 additional

FORTRAN statements. Furthermore, it requires little additional storage:

2 (NUMEL) + 2 (NUMNP)

where,

NUMEL = number of elements in the mesh

NUMNP = number of nodes in the mesh.

Stability proofs are not available for these subcycling procedures for second-order

systems; stability proofs have been developed only for linear first-order systems,

Belytschko, Smolinski, and Liu[7], Mizukami[9]. However, the energy stability checks

clearly indicate that the procedure developed here is stable and we have not encountered any

instabilities.

Because of the great computational savings which can be achieved by this method,

it is an attractive feature for any computer program used for engineering analysis by explicit

time integration.

23
ACKNOWLEDGEMENT

The support of NASA-Langley under Grant No. NAG-1-650 to Northwestern

University is gratefully acknowledged

24
References

[1] T. Belytschko and R. Mullen. Mesh partitions of explicit-implicit time integrators.

In: Formulations and Computational Algorithms in Finite Element Analysis , K. J.

Bathe, J. T. Oden, and W. Wunderlich, eds., 673-690 MIT Press, Cambridge, MA

(1976).

[2] T. Belytschko and R. Mullen. Stability of explicit-implicit mesh partitions in time

integration. Internal. J. Numer. Meths. Engrg. 12, 1575-1586 (1978).

[3] K. C. Park. Partitioned transient analysis procedure for coupled field problems:

Stability analysis. /. Appl. Mech. 47, 370-376 (1980).

[4] T. J. R. Hughes and W. K. Liu. Implicit-explicit finite elements in transient

analysis: Stability theory. /. Appl. Mech. 45, 371-374 (1978).

[5] T. J. R. Hughes and W. K. Liu. Implicit-explicit finite elements in transient

analysis: Implementation and numerical examples. J. Appl. Mech. 45, 375-378

(1978).

[6] T. Belytschko and R. Mullen. Explicit integration of structural problems. In: Finite

Elements in Nonlinear Mechanics 2, P. Bergen et al. eds., 697-720 (1977).

[7] T. Belytschko, P. Smolinski, and W. K. Liu. Stability of multi-time step partitioned

integrators for first-order finite element systems. Comp. Meths. Appl. Mech.Engrg.

49, 281-297 (1985).

25
[8] P. Smolinski, T. Belytschko, and M. Neal. Muld-time-step integration using nodal

partitioning. Internal. J. Numer. Meths. Engrg. 26, 349-359 (1988).

[9] A. Mizukami. Variable explicit finite element methods for unsteady heat conduction

equations. Comp. Meths. Appl. Mech. Engrg. 59, 101-109 (1986).

[10] J. Donea and H. Laval. Nodal partition of explicit finite element methods for

unsteady diffusion problems. To appear in Comp. Meths. Appl. Mech. Engrg.

[11] T. Belytschko and W. K. Liu. Time integration with explicit/explicit partitions in

EPIC-2. Report to Ballistics Research Laboratory, July, 1982.

[12] T. Belytschko. Partitioned and adaptive algorithms for explicit time integration. In:

Nonlinear Finite Element Analysis in Structural Mechanics , E. Wunderlich et al.

eds. 572-584 (1980).

[13] T. Belytschko. Transient Analysis, in: Structural Mechanics Computer Programs,

W. Pilkey et al. eds. 255-276 Univ. Press of Virginia (1974).

[14] T. Belytschko. Overview of semidiscrtization. In: Computational Methods for

Transient Analysis, eds. T. Belytschko and T. J. R. Hughes, 67-155 North

Holland, Amsterdam (1983).

[15] S. Key. Personal communication (1987).

26
List of figures

Table 1. Nodal grouping and time steps for example problem 1.

Table 2. Nodal grouping and time steps for example problem 2.

Figure 1. Belytschko and Liu[l 1] nodal partition for integer time step ratios

Figure 2. Belytschko and Liu[l 1] nodal partition for non integer time step ratios.

Figure 3. New nodal partition for non integer time step ratio case.

Figure 4. An example of the nodal grouping part of the subcycling procedure

Figure 5. Flowchart of the subcycling procedure

Figure 6. Example problem 1. A bar subjected to a sudden axial load

Figure 7. Stresses of the rod for the elastic case

Figure 8. Energy error of the elastic problem

Figure 9. Stresses of the rod for the plastic case

Figure 10. Energy error for the plastic case

Figure 11. Example problem 2, rotating ring with internal pressure

Figure 12. 36 element mesh for example problem 2

Figure 13. Maximum stress of 36 element mesh

Figure 14. 100 element mesh for example problem 2

Figure 15. Effective plastic strains for rotating ring

Figure 16. Effective plastic strains for the rotating ring with subcycling

Figure 17. Energy error of example 2
Table 1. Nodal grouping and time steps for example problem 1.

Group

A

B

C

D

Elements

1-10

11-15

16-25

26-32

Nodes

1-10

11-15

16-26

27-33

Element length

1.0

0.4

0.1

1.0

Subcycling

rime step

lOAt

4At

At

lOAt

No subcycling

rime step

At

At

At

At
Table 2. Nodal grouping and time steps for example problem 2.

Group

1

2

3

4

5

6

7

Elements

1-28

29-37

38-48

49-55

56-65

66-78

79-100

Nodes

3-54

55-72

73-94

95-108

109-128

129-154

155-200

and 1-2

Subcycling

time step

12At

8At

6At

4At "

3At

2At

At

Nosubcycling

rime step

At

At

At

At

At

At

At
1 1 22 3 3 4

Figure 1
1 1 2 2 3 3. 4

Figure 2
l l

2 3

Figure
8

r

i

1

10

Element rime steps

At! = 3.4

At2= 1.1

At3= 1.1

At4 = 2.3

Nodal time steps

x, At8 = 3.4

Ats.At6.At7 =1.1

= 2.3

Master time step

At =1.1

Truncated, nodal time steps

Ati, Atg = 3At = 3.3

At2. At3> At4,

5, At6>At7 =At =1.1

= 2At = 2.2

Master level time step

At j^ = least common multiple of (1.1,2.2,3.3) = 6.6

Master level time

1

 MLT =

 6'6

Figure 4
Initialization

1. Read in data and set initial conditions, tM = t e = t I = 0

2. Determine the values of At r, At M, and t MLS as described in section

VILA.

Tune integration loop

3. Calculate the external nodal forces.

4. Loop on the elements of the mesh (JE)

If any node of element JE satisfies equation (5.5)

THEN for element JE

A. Calculate At M.with equation (5.6)

B. Update the stress by equations (3.4) through (3.6).

C. Calculate the internal nodal forces and assemble with

equation (3.3).

D. Calculate the increment in internal energy for this element and

assemble to the internal energy of the system, equations (6.2) and (6.3).

E. At e=At M

5. Use equation (1) to calculate the acceleration vector, ai.

6. Loop on all the nodes of the mesh (I)

A. If node I satisfies condition (5.9)

THEN the velocity and displacement of node I are updated by equations

(3.7) and (3.8).

B. At?ld = Atj

C. If node I does not satisfy condition (5.9)

THEN the velocity and displacement of node I are not updated.

7. Calculate the external work done of the system with equation (6.4).

8. Calculate the kinetic energy of the system with equations (6.5) through

(6.6)

9. If this is a master level time step, i.e., t M = t MLT

THEN

A. Check the energy balance with equation (6.1).

B. If desired, recalculate the time steps, etc., as described in section VILA.

10. Output the results.

11. If the maximum time of the problem is not exceeded

THEN

A

- 'M

 = tM + AtMB. GOTO step 3.

12. If the maximum time of the problem is exceeded

THEN STOP.

Figure 5
f»-H A

1 f\ r\ , , 1U.U

B

2.0

,

c

1

 D

i

/. nu

1.0

Modulus of elasticity =

Density =

/

'/

'

100.0

1.0

Plastic modulus = 25.0

Yield point (isotropic hardening) =

Case fext

1 0.2

2 2.0

1.0

Figure 6
0.1

0.0

no subcycling

subcycling

exact

2.0 3.0

Time

4.0 5.0

Figure 7
D.U ~

4.0-

| 3.0-

£

to 2.0-

i— 1

<u

§ 1.0-

^

o.o-

-1.0.

«^ M«UA.MI:M<«

_^^^ _ pn i^fwrMiTirr

[/V.V,A_^.

0.0 1.0 2.0 3.0

Time

4.0 5.0

Figure 8
0.5

0.0

-0.5-

C/5

C/3 -i f\

<u -1.0

CO

-1.5-

-2.0-

-2.5

no subcycling

subcycling

exact

' v v

0.0 1.0 2.0 3.0

Time

4.0 5.0

Figure 9
4.0-

p

 3.0 H

no

1.0 2.0 3.0 4.0 5.0

Time

Figure 10
CO

E

'-tangent

P

"*•"

= 1.03 e+7 psi

= 0.333

= 5.17 e+5 psi

= 1 (isotropic)

= 4.15 e+4 psi

P

Ri

R2

CO

P

= 2.508 e-4

= 1.05

= 0.95

= 4000

= 9760 psi

Msec2 /in4

Figure 11
Figure 12
2.10

.60 1.20 1.80 2.40

Time (xlOe+4)

3.00

Figure 13
Figure . 14
§C/3

o

0.3

0.2-

-2

'S«

<u

•i. o.iO

ffl

0.0 no subcycling

subcycling

0.0 0.2 0.4 0.6 0.8 1.0

Time (xlOe+4)

Figure 15
.s

cd

t-4

4-^

C/3

00

"Su

element 10

element 98

0.0 0.2 0.4 0.6 0.8 1.0

Time (xlOe+4)

Figure 16
_1

g

1?

a

ZU.Ufi r\ _ u.u -

-20.0 -

x /*v /"v

-40.0 -

-&c\ n.

;/ ' li

i*

y &

0.0 0.2 0.4 0.6 0.8

Time (xlOe+4)

1.0

Figure 17