. .. .. Lhw .q ,
554.32.... ...,.

:mvmfiuw

.3 ., .

 

.1.

in
S...

.l...
.«x1

o

 

w. _. tm 7‘.
$1....“
1H4“...
his“.
, \.

 

3o... . 03.10. vx .. Id“: .o.vn4.v8uv¢!._. U u,....ll z. ,Fvflwdw: . . .t 1‘. a in fivth . .Jl.v4u.|?..4~.o.lan.‘. l):
f. , guru“. =§$fi§§m Ream? A.» . . Amiga? . gamma £4 Egg

 

1 . .

 

 

LIBRARIES
MICHIGAN STATE UNIVERSITY
EAST LANSING, MICH 48824-1048

This is to certify that the
thesis entitled

NUMERICAL INTEGRATION OF CONSTITUTIVE MODELS
APPLICABLE TO METAL FOAMS

presented by

Abhishek Modi

has been accepted towards fulfillment
of the requirements for the

 

 

 

MS. degree in Deparment of Mechanical
Ermineering
W 4 F
Major P{c}és'sor’s Signature
Ea C . / g I Z 0 0 4

 

Date

MSU is an Affinnative Adlai/Equal Opportunity Institution

 

PLACE IN RETURN BOX to remove this Checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 c:/ClRC/DateDue.p65-p.15

 

 

NUMERICAL INTEGRATION OF CONSTITUTIVE MODELS
APPLICABLE TO METAL FOAMS
By

Abhi shek Modi

A Thesis
Submitted to
Michigan State University

in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE
Department of Mechanical Engineering

2004

ABSTRACT
NUMERICAL INTEGRATION OF CONSTITUTIVE MODELS
APPLICABLE TO METAL FOAMS
By

Abhishek Modi

Metal foams have unique combinations of physical and mechanical properties and are
becoming increasingly popular among researchers and designers. Four elasto-plastic
material models for metal foams are studied in detail. All four material models
considered are rate independent elastoplastic models with isotropic hardening rules. A
return mapping algorithm using the Newton-Raphson closest point projection is discussed
for integrating the constitutive equations of the models. A way to do the return mapping
along the spectral directions of the stress tensor is described so as to reduce the
dimensionality of the problem. An algorithm applicable to all metal foam models, both
with associative and non-associative hardening, is presented. A user subroutine is written
in FORTRAN for implementing the algorithm. Numerical simulation of the diagonal
crushing of metal foam cube is done in FEM program ABAQUS to verify the material
model algorithm. The result of the numerical simulation using user defined subroutine, is
compared with the result using the inbuilt metal foam model in ABAQUS. The energy
absorption capacity of metal foams is demonstrated by comparing the crushing process of

a metal foam cube with a solid aluminum cube.

TABLE OF CONTENTS

LIST OF FIGURES .................................................................................. v
CHAPTER 1

Introduction to metal foams and their properties ................................................. l

1.1 Introduction .............................................................................. l

1.2 Applications of metal foam ............................................................. 2

1.3 Classification of metal foam ........................................................... 4

1.4 IssUes in manufacturing and application of metal foams ........................... 6
CHAPTER 2

Experimental studies on metal foam ................................................................ 9

2.1 Introduction .............................................................................. 9

2.2 Scaled equations for strength and modulus of foamed metals ..................... 9

2.3 Effect of micro-structural irregularities on strength and modulus ............... 12

2.4 Uniaxial stress strain response of metal foams ..................................... 13

2.5 Initial yield surface of foams ......................................................... 13
CHAPTER 3

Elastoplastic material models of metal foams ................................................... 16

3.1 Introduction ............................................................................. 16

3.2 General equations of elastoplastic theory as applicable to metal foams ........ 16

3.3 Description of selected material models ............................................ 21

3.3.1 LS-DYNA model ........................................................... 22

3.3.2 ABAQUS model ............................................................ 24

3.3.3 Deshpande-Fleck model ......................... : ........................... 26

3.3.4 Miller model ................................................................. 28

3.4 Comparison of Yield surfaces ........................................................ 31
CHAPTER 4

Integration of the constitutive equations .......................................................... 33

4.1 Introduction ............................................................................. 33

4.2 Constitutive equations in continuous form ......................................... 34

4.3 Time discreet form of the mathematical problem .................................. 37

4.3.1 Closest point projection algorithm ........................................ 37

4.3.2 Predictor-Corrector algorithm ............................................. 38

4.4 Stability and complexity of the numerical algorithms ............................ 40

4.5 Algorithm for two invariant plasticity models ..................................... 42

4.6 Implementation of algorithm for selected material models ....................... 48

iii

CHAPTER 5

Implementation of the algorithm ................................................................ 49

5.1 Introduction ............................................................................ 49

5.2 Input data for the model ............................................................. 49

5.3 Numerical simulation ................................................................ 51

5.4 Comparison of results ................................................................ 54

5.5 Conclusion ............................................................................. 56
Appendix I

Derivatives of the yield function and flow potential for all four models .................. 57
Appendix H

Fortran Subroutines for implementing algorithm in ABAQUS ............................. 66

Bibliography ....................................................................................... 76

iv

LIST OF FIGURES

Fig 3.1 metal foam yield surface in deviaton'c and meridian plane ........................... 20
Fig 3.2 evolution of yield surface and flow potential for LS-DYNA model ................ 23
Fig 3.3 evolution of yield surface and flow potential for ABAQUS model ................. 25
Fig 3.4 evolution of yield surface and flow potential for Deshpande-Fleck model ........ 28
Fig 3.5 evolution of yield surface and flow potential for Miller model ...................... 30
Fig 3.6 comparison of initial yield surface in p-q plane ....................................... 32
Fig 5.1 typical uniaxial compressive response of closed cell foam ........................... 50
Fig 5.2 model used for simulation of diagonal crushing ....................................... 51

Fig 5.3 diagonal crushing of cube with user defined material model

(Deshpande—Fleck model) ................................................................. 52
Fig 5.4 diagonal crushing of cube with metal foam model

included with ABAQUS .................................................................. 52
Fig 5.5 diagonal crushing of a solid aluminum cube ........................................... 53

Fig 5.6 Composition of internal energy in metal foam cube
during stepl and 2 .......................................................................... 55

Fig 5.7 composition of internal energy in solid aluminum cube
during stepl and 2 .......................................................................... 55

Chapter 1

Introduction to Metal Foams and their properties

1.1 Introduction

Foamed metals have become increasingly important to both researchers and designers in
the past two decades owing to the unique and interesting combinations of physical and
mechanical properties of these materials. Increased interest in these materials in recent
times is also due to improved production methodologies [6, 10] which result in easier
availability of foamed metals of better and consistent properties at affordable prices.
Metal foams have unique combinations of desired properties like impact energy
absorption capacity, air and liquid permeability, good acoustic properties as well as low
thermal conductivity and high electrical conductivity. One reason why cellular metals
have not been more popular in engineering design till recently is the poor manufacturing
processes available in the past. New manufacturing techniques developed in the last
decade have now ensured that foamed metals have a greater range and control of physical
properties. Increased process control over the manufacturing parameters has improved
the quality of the metal foams produced. The micro mechanical structure of these new
metal foams is more consistent leading to better mechanical properties which can be
predicted with more accuracy. In this chapter we have listed some of the applications for
metallic foams, classification of cellular metals and the commonly available metal foams

with their market names.

1.2 Applications of metal foams

Metallic Foams are important materials in engineering because of the interesting
combinations of physical and mechanical properties. The existing applications for metal
foams cover a wide range of fields [10]. The applications for metal foams can be broadly

divided into three categories based upon the structure and properties of these foams.

Foams with predominantly closed cells are useful for shock and impact absorption
applications. In energy absorption applications their structure allows these closed cell
foams to absorb large amounts of mechanical energy when they are deformed. Closed
cell foams have a high compressive strength and can undergo large plastic deformations
because of their porous structure. The automotive industry is looking at metal foams for
use in front bumpers to use their high impact absorption capacity. Some other suggested
uses for these materials are as protective envelope for air borne safety equipment like
aircraft black boxes. Disposable landing feet for spacecrafts and disposable clamping

fixtures are other places these materials can be used.

Open cell foams, because of their air and fluid permeability, can be used for functional
applications like filtration, damping etc. Metal foams have good thermal and acoustic
properties which can be used to dampen vibrations and absorb sounds in certain
conditions. Dust and fluid filters, engine exhaust mufflers are some areas these properties
are useful. Inert metal foams (ex. gold, silver) are used for the filtration of corrosive
liquids. Other useful properties like low thermal conductivity in conjunction with high

structural strength is desirable for applications like high temperature gaskets, flame

arrestors etc. Metal foams are good electric conductors which make them good
candidates for application as porous electrodes. Foamed metals have also been used

advantageously in heaters and heat exchangers.

Foamed aluminum is being used as core material in sandwich panels which are used for
structural, load bearing purposes. When used as the filler in sandwich panels they take up
a lot of volume without considerably increasing the overall weight of the structure. These
three dimensional panels which are structurally very stiff can be used effectively in the
aircraft industry as a low cost replacement of honeycomb panels. The automotive
industry is also looking at these panels to use in the body of the car. The high structural
strength of these panels coupled with the thermal properties of the porous metals has
found application as wall and floor heating tiles. Open cell metal foams when used as the
core of reinforcing structural struts allow the flow of liquids through them making them

useful for heating or cooling purposes.

The above list of applications for cellular metals is not a complete list of all the places
these remarkable materials with a combination of useful properties have found
application. As the manufacturing processes for metallic foams improve leading to better

quality foams at lower prices, designers will use them more widely.

1.3 Classification of metal foams
Cellular metals are defined as metals with solid and gaseous phases in their morphology.
To properly identify metallic foams which are a subset of the larger class of cellular

metals, proper classification has to be made for all kinds of cellular metals [6].

Cellular metals- This is the most general classification and refers to a metallic body in
which any kind of gaseous voids are dispersed.

Porous Metals- These are a special kind of cellular metal. The voids are pores which are
round and isolated from each other.

Metal Foams- These are a kind of cellular metal with restricted morphology. The cells
are closed and are separated from each other by thin films.

Metal Sponges- These are open cell cellular metals with interconnected voids.

The subject of this study will be limited to Metal Foams and Metal Sponges.

New and improved manufacturing techniques have resulted in the availability of better
quality metal foams at lower prices. Banhart [6] and Davies [10] have listed the new
manufacturing techniques being used by leading manufacturers of metal foams in the
industry. Alcan, Alporas and Alulight are the trade names of the most commonly
available closed cell aluminum foams. The manufacturing processes of these common
types of aluminum foams are described below along with the range of their physical

properties.

Trade name: ALCAN

Manufacturers: Hydro Aluminum in Norway & Cymat Aluminum Corporation in Canada
Manufacturing process: This is a foaming method for obtaining metal foam from
aluminum melt. The foaming of aluminum melt is done by gas injection. Silicon-carbide,
aluminum-oxide or magnesium-oxide particles are used to increase the viscosity of the
melt. This helps in the stability of the foam when it is being formed. This step requires
sophisticated mixing techniques to ensure a uniform distribution of particles. The melt is
formed in the second step by injecting gas. The foaming is done in a continuous process
in which the mixture of gases and solid metal floats to the top while the liquid metal
drains out. Because of the gravitationally induced drainage, a gradient in density, pore
size and pore elongation is observed in the resulting foam. By this process the foams

obtained have the following average properties:

Average densities 0.069 - 0.54 _g_
cm3

Average pore size 3mm-25mm

Trade name: ALPORAS

Manufacturers: Shinko Wire Company in Japan

Manufacturing process: Metal/Aluminum melts can also be foamed by directly adding a
blowing agent in the melt instead of injecting gas into it. First some 1.5% wt Ca metal is
added to the melt for increasing the viscosity by the formation of calcium oxide and
calcium-aluminum oxide. After the viscosity has reached the desired level, titanium
hydride is added (1.6% wt), serving as a blowing agent by releasing hydrogen gas in the

hot viscous liquid. After cooling the vessel below the melting point of the alloy, the

liquid foam turns into solid aluminum foam. By this process the foams obtained have the
following average properties:

Average densities 0.18 — 0.24 i—

cm3

Average pore size 2mm-10mm

Trade name: ALUIJGHT

Manufacturers: F raunhafer Institute in Germany

Manufacturing process; Foamed aluminum can also be prepared from metal powder.
Aluminum powder or alloy powder are mixed with the right amount of blowing agent,
after which the mix is compacted to yield a dense product. Next, this semi-finished
product is heat treated to temperatures close to the melting point of the matrix material.
This leads the blowing agent to decompose and release gas, which in turn forces the
precursor material to expand and foam. An advantage of this process is that the expansion
during foaming can be controlled inside a mold which results in the formation of near net

shapes of the foam material.

1.4 Issues in manufacturing and applications of metal foam

There has been relatively limited use of metal foams in engineering design despite their
impressive combinations of properties. Notwithstanding the recent developments in
manufacturing of metal foams, there is a lack of sophistication in manufacturing
processes when compared to polymeric foams. Some of the important issues hindering in

the manufacturing and application of metal foams are listed below.

Stability of foams: foams are unstable systems because their large surface area causes
energy to be far from a minimum value. Foams can at best be meta-stable, constantly
decaying at a certain rate. One way to increase the stability of foams is to introduce solid
particles in the melt to increase the viscosity. The action of foam stabilization by the
addition of particle additives is not entirely clear. Some recent work has been done to

study the stability of metal foams [9].

Consistent quality of foam in production: availability of consistent foams in a wide range
of densities is important for wide scale industrial applications. There is a lack of
sufficient knowledge of the mechanisms of foam formation to be able to produce foams

of consistent quality with pro-defined parameters.

Expense of metal foams: There have been some recent improvements in the field of
manufacturing, but there is still need for setting up efficient industrial processes to reduce
the cost and improve the quality of foams. Metal foams are nevertheless less expensive
than honeycomb panels that are popular in the design of weight bearing structures in the

aircraft industry.

Physical properties of foam: In recent years there has been considerable research interest
in experimental studies of the physical and mechanical properties of available metal
foams [2, 12, 17, 18]. There is still a need to establish a wide database for the physical

properties like density, average cell size, cell wall irregularities, and cell wall thickness

and mechanical properties like compressive strength, tensile strength, hydrostatic

strength, elastic modulus and energy absorption capacity.

Modeling of the mechanical response of foams- It is increasingly important for design
engineers to be able to predict the large strain response of materials they use in design.
Efficient and accurate material modeling of metal foams is required for use in FEM
packages. There has been some research which treats metal foams as solid continuum and

applies the equations of elastoplasticity theory [11, 19].

The focus of this study is to provide an algorithm for the numerical integration of the

elastoplastic constitutive equations of metal foams.

Chapter 2

Experimental studies on metal foams

2.1 Introduction

There is a sizable and growing amount of literature on the experimental studies done on
metal foams to record their physical and mechanical properties. Several researchers have
done experimental studies on the uniaxial and biaxial response of foamed metals [17, 18].
A majority of these experimental works have focused on the uniaxial stress strain
response of metal foams both in compression and in tension. Some work has been done to
summarize the results of these tests in simple theoretical equations [2]. Closed cell foams
exhibit remarkably lower values of strength in experiments than the theoretical
predictions. The micro structure of foams has been studied in some detail to explain the
apparent drop in the physical properties of closed cell metal foams. Apart from uniaxial
experiments, some biaxial loading studies [11, 12] have been performed to measure the
performance of foamed metals under various stress situations. Such experiments are
useful in determining the initial yield surface of foams. A presentation of results from the

experimental literature, most pertinent to the present study, is given in this chapter.

2.2 Scaled equations for strength and modulus of foamed metals
Andrews et al [2] have attempted to use simple dimensional arguments to model the yield
strength and elastic stiffness of closed cell and open cell foams in terms of simple

equations. The most important physical property which determines the mechanical

properties of foams is the density of the foam relative to the density of the cell wall
material. Dimensional arguments can be used to relate the Young’s Modulus and yield

strength of the foam to the relative density with respect to the cell wall material.

The mechanical properties of foams can be modeled by considering the mechanisms by
which the cells deform and fail. Open cell foams and closed cell foams fail by different
mechanisms. Open cell foams under uniaxial stress deform by the formation of plastic
hinges in the cell walls. In Closed cell foams on the other hand, bending of cell walls is

accompanied by stretching of cell faces.

A structural mechanical analysis of a low density, open cell, Kelvin foam with a
tetrakaidecahedral structure leads to the following equation for the modulus and yield

strength of open cell foams [2].

a: a- 2
§—=0.98[p—] (2.1)
E ,0
3
at: at: 3
EC. = 0,3 :0— (23)
cc ,0
where

E* - elastic modulus of open cell foam

0': - uniaxial compressive yield stress of open cell foam.

E - elastic modulus of cell wall metal (aluminum)

ac - uniaxial compressive yield stress of cell wall metal (aluminum).

10

*

£— - relative density of the foamed metal with respect to the base metal.

,0

A similar analysis for ideal tetrakaidecahedral closed cell structure leads to the following

set of equations [2].

* * 2 *
E—=o.32 3— +o.32 9-- (2.3)
E .0 P

. .. 2 ..
51:0.” ”— +0.44 £- (2.4)
0c p p
where now

E*- elastic modulus of closed cell foam

0': - uniaxial compressive yield stress of closed cell foam.

Comparison of the predicted values of modulus and strengths from these equations to the
experimental results from studies on the stiffness and strengths of the open cell foam
ERG and the closed cell foams Alcan, Alporas & Alulight show interesting results [2].
The theoretical equations for the open cell foam (Eq 2.1 and 2.2) are quite close to the
experimental results for ERG, but the theoretical equation for closed cell foam (Eq 2.3
and 2.4) over-predicts the experimental results from Alcan, Alporas and Alulight quite
appreciably. The reason for the observed difference between the model and experimental
values is the presence of micro-mechanical imperfections in the morphology of closed

cell foams.

ll

2.3 Effect of micro-structural irregularities on strength and modulus
Micro-structural observations of the closed cell foams (Alcan, Alporas and Alulight)

reveal a number of defects in the structure which is not taken into account in the

theoretical model [15]. Some of these defects are listed below

Cell wall curvature: The walls of some of the cells are not plane but have some initial
curvature which is not taken into account by the model

Cell wall corrugations: The walls of some cells have corrugations and chinks.

Density gradient: The density of the foam varies throughout the thickness of the foam due
to the processing conditions. This is pronounced in Alcan which is cast under the
influence of gravity.

Defects in cell walls: some of the cell walls are missing and others have voids.

The effect of cell wall curvature and corrugations has been modeled using FEM and
compared to the strength of a periodic tetrakaidecahedral cell [15]. The observed cell wall
imperfections can account for up to a 70% drop in the mechanical properties of foams.
The presence of these micro structural defects is a direct result of the production methods
for foams. One of the important challenges for the manufacturers of foamed metals is to
develop production methodologies that increase the uniformity of the microstructure of

foams.

12

2.4 Uniaxial compressive stress strain response of metal foams

Experimental studies on the uniaxial compressive response of various types of closed and
open cell metal foams reveal certain similarities [7]. All foams irrespective of the cell
wall material exhibit the schematic response under uniaxial loading. Initial loading results
in an increase in stress magnitude which appears to be proportional to strain. But the
unloading modulus is greater that the loading modulus suggesting that local yielding
takes place almost immediately on loading. In a macro sense the material can be
considered elastic up to a certain strain. After the stress reaches the yield point, there is a
regime of strong plastic deformation with very little strain hardening. Sometimes a
separate upper and lower yield point is observed. After an extended plateau the stress
strain curve enters a regime of rapid strain hardening. This part of the stress strain curve
is characterized by densification of cell walls as they collapse on one another,
accompanied by a rapid increase in stress. The strain at which the stress strain regime
enters densification is called the densification strain. Fig. 5.1 shows a typical uniaxial

compressive response of closed cell aluminum foam.

2.5 Initial yield surface of foams

The initial yield surface for metal foams is important in the mathematical model of the
elastoplastic response of foams. Some phenomenological yield criteria have been
suggested by Deshpande-Fleck [11], Miller [19] and Doyoyo [12] among others, for
foamed metals. Experimental studies for biaxial and multiaxial loading of foams are

limited when compared to the large experimental literature available for uniaxial loading.

13

Nevertheless attempts have been made in this regard by Deshpande-Fleck [11] and

Doyoyo et al [12] among others.

Doyoyo et al [12] have investigated the yield surface of closed cell aluminum foam using
the modified Arcan apparatus. Butterfly shaped test specimens of the foam were used in
the Arcan apparatus which ensures the yielding of the foam is at the center of the
specimen.

The Yield surface was found to be pressure dependant and quadratic in the meridian
plane. A phenomenological yield surface has been proposed in the principal stress space

which depends on the first invariant of principal stresses (II) and second invariant of

deviatoric stresses (12 ).

II = oi, (2.5)
12 £5,ij (2.6)
Here

oi]. = %ou 6,]. + Si). (2,7)

aij - SITCSS tensor

Si]- - deviatoric stress tensor

Deshpande & Fleck [11] have performed biaxial loading experiments for the yield
surface and have come up with a similar pressure dependant yield surface in the principal

stress space. They used triaxial system consisting of a pressure cell and a screw driven

l4

piston for axial loading. The probing of the yield surface suggests a quadratic shape in the

stress space of mean stress (0,") and Von Mises effective stress (Ge).

where
l
Um = 30-h: (27)
3
03 = ESijSij (2.8)

15

Chapter 3

Elastoplastic material models of metal foam

3.1 Introduction

Metal foams are heterogeneous because of the gaseous voids present in the morphology.
Metal foams are therefore not continuous materials in the strict sense of the definition.
The elastoplastic theory applicable to continuous materials however, can be applied to
metal foam specimens which are considerably larger in dimension than the average cell
dimension. The rule of thumb followed in some references [2, 12] is that the smallest
specimen dimension should be at least 10 times the average cell diameter for continuum
theory to be applicable. Several attempts have been made to model the mechanical
response of metal foams under the broad theory of elastoplastic materials. Commercially
available FEM programs LS-DYNA [13] and ABAQUS [1] have material subroutines for
metallic foams incorporated in the program. Two other material models in literature by
Deshpande-Fleck [11] and Miller [19] are also studied in detail. In this chapter we
present the general equations of elastoplastic theory as applicable to metal foams and also

the constitutive equations of the four models.

3.2 General equations of elastoplastic theory as applicable to metal foams

Uniaxial experiments [18] performed on metal foams have shown that the unloading
elastic modulus of foams is higher than the loading modulus even for small strains. This
is explained by the fact that metal foams undergo local plastic yielding almost

immediately on loading. The local yielding is due to the presence of irregularities at the

16

cellular level. The overall specimen can still be treated as elastic at the macroscopic level
up to the yield stress. All the four models considered assume isotropic elastic behavior in

the elastic region with a constant elasticity tensor.

0:7 = Cans/f1 (3.1)
where

Cijkl - elasticity tensor

£5- - elastic strain

since isotropic elastic behavior is assumed, Cu“ is completely defined by any two elastic

'1

parameters like modulus of elasticity (E) and elastic Poisson’s ratio ( v ).

Elastoplastic theories in general assume the strain rates to decompose additively. So we

can write

é”- : e;- + a; (3.2)
where

815-, - plastic strain

Isotropic elastic behavior is assumed until yielding is reached. Once yielding is reached,
the inelastic response is assumed to be strain rate independent by all four material

models. A yield surface in stress space is used to indicate yielding.
f(a,.j,e,.5.’) = o (3.3)

The presence of 8,5? in the yield rule is to take into account the strain hardening of the

material. As explained later, the yield surface of metal foams can be conveniently

17

represented in the stress space by the Von-Mises equivalent stress (q) and hydrostatic

pressure stress (p). where p and q are defined as follows

1

p = -3017 (3.4)
3

q = 55.75:;- (3.5)

The hardening of the yield surface depends on the amount of inelastic strain in the
material. Three of the models (ABAQUS, Deshpande-Fleck, and Miller) assume an
isotropic hardening rule for strain hardening. LS-DYNA uses a hardening rule that allows
for change in the shape of the yield surface. The hardening for all four models can be
given in terms of two scalar quantities that are used to track the growth of plastic strain.
The following strain-like scalar internal variables are introduced

3” - equivalent plastic strain

85’ - volumetric plastic strain

The equivalent plastic strain is defined later (Eq 4.7). The volumetric plastic strain can be
written as

6,? = eff (3-6)

Now we can write the yield surface in terms of four scalars in the following manner

f (p.q.£"’ .83." ) = 0 (3.7)

Elastoplastic theories in general define the rate of inelastic strain in terms of a flow

potential function
. a
£5 = r—a: (3.8)

1.l

18

where

7- plastic multiplier

¢ - flow potential function

The flow potential function (¢) is defined in the stress space and can evolve with

accumulation of plastic strain just like the yield potential function (f). In fact a lot of
elastoplastic models assume the flow potential to be the same as the yield potential. If this
is the case then the flow rule is said to be associative, otherwise it is said to be non-
associative. The flow rule is associative for two of the models (Deshpande-Fleck &
Miller) while the other two models (LS-DYNA & ABAQUS) use a non-associative flow
rule.

The strain hardening of the yield function and flow potential is incorporated in the
models through the input of experimental data. The four models considered use the data

from uniaxial compressive test and/or hydrostatic compressive test as input. The

equivalent plastic strain (ép ) is used as a measure of the uniaxial plastic strain for all

models except the Miller model (see Sec 3.3.4)

Pc = Fc (gvp ) (3'9)
ac = fic(§p) (3.10)
where

pc - hydrostatic compressive yield stress
ac - uniaxial compressive yield stress

The yield surfaces of elastoplastic materials are most conveniently represented in the

principal stress space. The complete representation of the yield surface in principal stress

19

space requires the description of cross-sectional shape in the deviatoric plane and its trace
on the meridian (p-q) plane. All four models considered here assume the material to
remain isotropic in the inelastic region leading to circular cross-sections in the deviatoric
plane. Therefore the models can be compared in the meridian plane generated by Von
Mises equivalent stress (q) and hydrostatic pressure stress (p). The yield surface for all
the four material models studied here are pressure dependant and quadratic in the p-q

plane.

0'1

 

 

 

03 o ‘ 0
pt pc p
Deviatoric plane Meridian (p—q) plane

Fig 3.1 metal foam yield surface in deviatoric and meridian plane
Experiments done on probing the initial yield surface of metal foams reveal that all metal
foams share some important characteristics [16]. All metal foams have a pressure

dependant yield surface that is quadratic in the p-q plane. The yield functions of all four

material models share this property. The other characteristics of the yield surface of foam

20

can be described in terms of some experimental parameters. We define these parameters

below

0'? - initial compressive uniaxial stress
p2 - initial compressive hydrostatic stress
p? - initial tensile hydrostatic stress

VP - plastic Poisson’s ratio

0

a . . . . . . . .
= —g - ratio of the initial compressrve and tensrle unraxral stress
at

The four material models studied, have yield surface functions that can incorporate one or
more of these additional experimental parameters. The complexity of the yield surface
depends on how many parameters it can accommodate. Miller model [19] and ABAQUS
[l] are both three parameter models while DYNA [l3] and Deshpande-Fleck model [11]

are two parameter models.

3.3 Description of selected constitutive models
Four elastoplastic material models for metal foam are studied in detail and the description

of their yield surfaces, hardening rule and flow rule are given in this section.

21

3.3.1 LS DYNA model (material model # 75) [13]

LS-DYNA has an inbuilt material model for metal foams as material model # 75.

 

 

 

 

The yield function is given by
f (p.q.£"’ .6.” ) = 0 (Eq. 3.7)
where
r 1 2
P'EIPc" p.) q 2
f= + (-) -1 (3.11)
V a b
\
Pc
= _ 3.12
Pt 10 ( )
1
a= 3(1). +pc) (3.13)
1 1 2 .
M = Gel (Peptgadpc aid-30'. (3.14)
b = Ma (3.15)

As can be seen from the above equations, the initial calibration of the yield surface
requires the knowledge of 0'? and p2 . The material model can therefore be described as

a two parameter model.

The strain hardening uses two sets of experimental data as input. The values of ac and
pc in the above equations are obtained from the results of uniaxial compressive test (Eq.

3.10) and hydrostatic compressive stress (Eq. 3.9) respectively.
The ratio of the half axes of the ellipse (M) is not a constant and so a and b can evolve

independently from each other. Hence the shape of the yield surface can change during

22

loading. The yield surface remains an ellipse in p-q plane and no comer development

takes place (Fig 3.2)

The flow rule is not associative and is given by

. 8¢
P = _ .8

where

(Penna/q2 + %p2

 

 

cl
hardened YS 1
\ 0".
initial YS
\\
flow rule 0':
/ I \ \
I] \\
I \
I \
l I
I, ‘\
f I
J '.
I l
I l o 1
pc pc

(3.16)

Fig 3.2 evolution of yield surface and flow potential for LS-DYNA model # 75

23

3.3.2 ABAQUS metal foam model (metal foam plasticity with volumetric hardening) [1]
ABAQUS has two inbuilt material models: metal foam plasticity with volumetric
hardening and metal foam plasticity with isotropic hardening. We will consider metal
foam plasticin with volumetric hardening.

The yield function is the same as the one in DYNA and is given by

 

 

 

f(p.q.é".e.f’) =0 (Eq 3.7)
where
1 2
P§(pc- pt) q 2
f = + [—) -1 (3.17)
a b
1 o
a = 5(1). +13.) (3.18)
1 1

M = 03/J(p3p9303(p3 nib-3032 (3.19)
b = Ma (3.20)

Tensile hydrostatic yield stress (pto) is a constant and is independent of the compressive

hydrostatic yield stress ( Pc) for this model. The initial calibration of the yield surface

requires the knowledge of 0'? , p2 and p?. The material model can therefore be

described as a three parameter model.

The hardening rule for this model is more simplistic than LS-DYNA and is governed

completely by the hydrostatic compressive yield stress - volumetric strain curve.

p. = We?) (Eq 3.9)

The ratio of the half axes, M is a constant during the plastic flow which conserves the

shape of the ellipse. Thus the hardening in this case can be considered isotropic (Fig 3.3)

The flow rule is not associative and is the same as the one in LS-DYNA.

. EM
.1.) = —
a", y 30,-]-

where

¢(p.q)=‘/qz + 3132

hardened YS
initial Y

flow rule

to

 

(3.21)

Pc

Fig 3.3 evolution of yield surface and flow potential for ABAQUS metal foam model

25

3.3.3 Deshpande & Fleck model [11]

Unlike metals, metal foams undergo volume reduction under hydrostatic compressive

pressure. The plastic Poisson’s ratio(v" ) can be taken as a measure of the

compressibility of materials. Solid metals which do not undergo volume reduction under
pure hydrostatic pressure have a constant plastic Poisson’s ratio (VP = é). Experimental

studies done by Deshpande & Fleck indicate that in metal foams the plastic Poisson’s

ratio is an important parameter that can be directly related to the yield surface ellipticity.

They introduced a shape factor a that is related to plastic Poisson’s ratio (VP ) by

a2 =3(l—-3VL) (3.22)

2 (1+VP)
They further suggested a two parameter yield surface that is quadratic in p-q plane.
Initial calibration of the yield surface requires knowledge of shape factor a and initial

uniaxial compressive stress Y °. The yield function is given by

 

f(p.q.€”.£.f’)=0 (Eq3.7)
where
f = __l_(q2 + a2p2) —0'c (3.23)

 

PEI

Note that v” = 0.5 gives elliptic shape factor a = Owhich reduces the Deshpande-Fleck
yield criterion to the Von Mises yield criterion for metals

q—a'c =0 (3.24)

26

Deshpande-Fleck have suggested two different hardening models in their paper.
Differential hardening model allows the yield surface to change shape whereas self-
similar hardening preserves the shape of the yield surface during plastic deformation. The

differential hardening model has a yield function of the form

2 2
Asia-1

Q and P are the yield strengths under deviatoric and hydrostatic stress respectively. Since
these can evolve independently of each other, the shape of the ellipse can change while
always remaining symmetric in tension and compression. Different yield strengths under
tension and compression are therefore not incorporated in either of the Deshpande-Fleck
models.

We have considered only the more simple case of self-similar hardening that is sufficient
for most purposes and is governed by the uniaxial compressive stress and equivalent

plastic strain curve.
ac = 6‘12”) . (Eq 3.10)

The flow rule is associative

 

. a
8.5? = r 3:. (Eq 3.8)
l]
where
. l

¢(p.q.£p,£f) = ————i-(q2 + azpz) ~00 (3.26)

1+ 9'-

3

27

hardened YS \

initial Y

 

 

 

Fig 3.4 evolution of yield surface and flow potential for Deshpande-Fleck model

3.3.4 Miller model [19]

Ron Miller has proposed a yield surface to describe the plastic flow of metal foams by
starting with the Drucker-Prager material model originally proposed for frictional
materials like soil.

The Drucker-Prager yield surface is

q-i’p-Uc =0 (3.27)
By adding a quadratic pressure term to the Drucker-Prager model one can include several

important features of metal foam plasticity. The new model requires three parameters for

28

the initial calibration which can be obtained by the results of uniaxial compression and

tension experiments. The yield surface is now given by

f(p,q,£‘p,€f)=0 (Eq3.7)
where

a

2
- 0' 3.28
(100's P Cl0 c ( )

 

f = q - 7p +
The model has a three parameter yield surface. Apart from the initial compressive
uniaxial yield stress (0'? ), the model incorporates the ratio of the uniaxial compressive
and uniaxial tensile stresses ( ,3 ) and the plastic Poisson’s ratio (vp ).

The coefficients y,a and do are the functions of ,6 and vp

[6fl2-12fl+6+9((f::))]

2(/3+1)2

 

r: (3.29)

_ 45+24y—4y2 +4vp(2+vp)(6y—yz -9)
l6(1+vp)2

do=l[1_Z+J(1_Z) +52] (3.31)
2 3 3 9

a

 

(3.30)

 

The model incorporates hardening through a function that is explicitly entered as input.
The uniaxial stress strain curve is used for the input data. A novelty in this model is the
multiplicative decomposition of the uniaxial stress-strain curve into hardening caused by

strain hardening of base material and hardening caused by cell compaction.

29

ac = "c(£‘p )17(£,f’ ) (3.32)
The advantage of this method of modeling hardening is that the large densification

hardening due to cell collapse can be reached only through large values of volumetric
plastic strain (85’ ). This is unlike the other models like Deshpande-Fleck where

densification stress levels can be reached even in a case where volumetric plastic strain

(6,? ) is low but equivalent plastic strain (1?” ) is high.

The flow rule is associative

 

35: 3.9.. (Eq3.8)
where
,,sP,sP = - + a 2—do 3.33
¢(Pq v) q 7P doacp O c ( )
q

 

 

 

P

Fig 3.5 evolution of yield surface and flow potential for Miller model

30

3.4 Comparison of Yield Surfaces

The complexity of the yield surface depends on the number of experimental parameters
that are used for calibration. Two of the four models (ABAQUS and Miller) are three
parameter models while the other two (DYNA and Deshpande-Fleck) are two parameter

models. Moreover the experimental parameters that these models take as input are

different. For example, ABAQUS uses the values of compressive uniaxial stress (0'? ),

compressive hydrostatic stress ( p2 ) and tensile hydrostatic stress (p? ). On the other

hand, Miller model which is also a three parameter model uses the values of compressive
0

unraxral stress (0'? ), ratro of compressrve to tensrle un1ax1al stress ( ,8 = —‘6) and plastic
0':

Poisson’s ratio (VP ).

In order to compare the shape of the yield surface of these four models we have

calibrated all of the yield surfaces for the same parameters. The values we choose are

02:1180kPa
o _
pc —690kPa
o _
p, —200kPa

ABAQUS and Miller, being three parameter models, pass through all three points while

DYNA and Deshpande-Fleck can accommodate only the first two.

31

q(kPa)

1500

02:11:30

 

L

 

-1000 1000

 

500 p2 = 690
pUcPa)

Fig 3.6 comparison of initial yield surface in p-q plane

32

Chapter 4

Integration of the constitutive equations

4.1 Introduction

Integration of the constitutive equations of the inelastic material models is an important
requirement in the finite element analysis of mechanical structures. All four material
models presented in the previous chapter fall in the category of rate independent
elastoplastic models. Typically rate independent elastoplasticity models define the rate of
inelastic strain in terms of stress and a set of internal state variables (Eq. 3.8). The
solution of the nonlinear differential constitutive equations is typically performed
incrementally and numerical integration becomes necessary. In a finite element method
for solving boundary value problems the integration of the constitutive equations is
performed at each quadrature point for every time increment in a strain driven
framework. The mathematical problem then reduces to a local one in which the updated
stresses and state variables are sought for a given strain increment and initial values of
the internal variables. In problems of finite deformations, where strain increments are two
or three times the yield strains, the stability of the numerical integration scheme becomes
important. An unconditionally stable algorithm specific to the material models of metal

foams is described in this chapter.

33

4.2 Constitutive equations in continuous form

As seen in the previous chapter, the yield potentials (f) and flow potentials (¢) of metal

foams are conveniently represented in terms of the deviatoric stress (q) and hydrostatic
stress (p). It is therefore instructive to write the stress (031-) in terms of p and q.

oi}- = -p6,-j + Si]- (4.1)
where

Sij - deviatoric stress tensor

we can then write

2
0,, = —p5,.,. +?‘1n (4.2)

if
where nij is the unit tensor that gives the direction of the stress point in the deviatoric

plane

(4.3)

The track of the plastic strain (65 ) in the material is kept by the scalar quantities 5“" and

8,? as mentioned in chapter 3. These quantities along with p and q are used to define the

yield rule

f(p.q,ép.el’)=0 (Eq3.7)
and flow potential

“p.23” .65) (Eq 3.8)

34

The complete description of the state of a material point therefore requires the knowledge

of oij, a” and 8”.

The flow rule gives the rate of plastic strain (Eq 3.8)

3¢

80,-]-

épzy—

using Eq. (4.2) we can break it down along p and q directions

ép_y[a¢_ ap a_¢ ]_a_q

”’ mam, an
_ 13¢ 8;»
35 _ y[—§3-6,~ +— aq nij] (4.4)

Now we will define the equivalent plastic strain (5” ) in terms of the plastic multiplier

(7). Consider a complementary equivalent stress (6') corresponding to the equivalent

plastic strain (3” ) . Then we get

A 2p: 8p
08 011-8 ’1
59¢
: yaij 8——
”it

using Eq. (4.2) and Eq. (4.4) we get

~31) = y[pa_¢.+ 431’) (4.5)

35

d¢ 29¢

the term p5—+qa— has the dimensions of stress and can be considered as a measure of
P q

the equivalent stress. Lets define

.. 3;» w
= __ _ 4.6
0 17 3p +q aq ( )
then we get
3” = y (4.7)

Eq (3.8), Eq (4.7) and Eq (3.6) give the rate of plastic strain in terms of the plastic
multiplier (y).

Kuhn-Tucker condition is used to maintain the consistency of the yield function in the

elastic region and the plastic region. The condition is given by

7f = 0 (4.8)
such that

f .<_ 0 (Elastic region)

y 2 0 (Plastic region)

When the loading is in the elastic region, the stress is within the yield envelope (f S 0)
and there is no plastic growth (7 = 0). In the plastic region, r 2 0 and the stress point is

on the yield envelope ( f = 0)

These equations can be taken together to state the mathematical problem. In a strain
driven framework the mathematical problem of integrating the constitutive equations can

be stated in continuous form as follows

36

Given £17

find

aij’ 5P and 8vp
subject to

. d

a; = 78—0?- (Eq. 3.8)

U

5;” = y (Eq. 4.7)
.93," = .95 (Eq. 3.6)
7f = 0 (Eq- 4.8)
such That

0,-1- = €in (8“ —£,5) (Eq. 3.1)
y 2 O

f S 0

4.3 Time discreet form of the mathematical problem

4.3.1 Closest point projection algorithm

The continuous differential equations in sec 4.2 cannot be solved in a closed form and
they need to be integrated numerically in time. For ensuring the stability of the
algorithms in dealing with large strains, we choose backward Euler time integration. In
the following equations all the quantities, unless explicitly mentioned are calculated at

time t+At.

37

The mathematical problem for the time discreet form now takes the following form
Given

Ag"

find

0'-- 6"” and £5)

1] ’
subject to

.3 = e£| EVE“; (4.9)
I

U I} 30--
I]

3" = éPL +Ay (4.10)

85 = 65 (4.11)

Ayf =0 (4.12)

such that

=Cijk,(£k,—£5) (Eq3.l)

4.3.2 Predictor-Corrector algorithms
The numerical solution of the mathematical problem defined above is usually

accomplished following a predictor/corrector strategy. A common consideration for the

predictor is the elastic trial state defined as follows

38

 

 

€5,trial = £15) (4.13)
t

gem“! = g!" (4.14)
r

E‘erial = Evp (4.15)
t

Then we have

034a: = W (8” _ gamut) (4.16)

ftrial ___ f(o.:jrial ’éerial ’E‘pdrial) (4.17)

The advantage of this elastic trial state is that we can determine if the strain increment is
elastic or plastic simply by considering the trial state.

Elastic loading

If f’”“’ s 0 then A7 = o

In this case the strain increment is elastic and the trial state is the final state

a; = 55”“ (4.13)
5p ___ éerial (4.14)
(if = Efrrial (4.15)
Plastic loading

If fl“ 20 then A720

In this case the increment has a plastic component and the plastic corrector is required to

restore the consistency of the yield surface. The mathematical problem is now of the

39

same form as in section 4.3.1 except that there is no need of the Kuhn-Tucker condition.
Using Eq. (4.9) -— Eq. (4.11), the problem is restated as

Given

A817
find

o-- 3” and 8,,”

11’

subject to

at = 53' +Ay— (Eq. 4.9)
t

5P =3P| +Ay (Eq. 4.10)

55’ :39 (Eq. 4.11)

such that

adj:cljkl(£k1_£kpl) (Eq3.1)

4.4 Stability and complexity of the numerical algorithms

Return mapping algorithms employing the predictor/corrector operator split strategy are
the most common methods for the integration of the elastoplastic constitutive equations.
Ortiz & Simo [21] have studied the accuracy and stability of the generalized trapezoidal
and generalized midpoint rules for time integration of the constitutive equations. They

have found that for finite deformations, where the strain increments are the order of a few

40

percent, backward Euler time integration provides unconditional stability. If backward
Euler method is used to discreetize the constitutive equations in time, the algorithm is
called closest point projection algorithm. The mathematical problem statement in the
previous section (4.3.2) is the general set of constitutive equations employing backward
Euler method and the predictor/corrector method. This set of implicit equations can be
solved using Newton-Raphson iteration at each time step. If Newton-Raphson iterations
are employed for solving the backward Euler time integrated set of constitutive
equations, the algorithm is referred to as Newton-Raphson closest point projection
algorithm.

Newton-Raphson closest point projection is the preferred method to solve constitutive
equations becauseof unconditional stability of the algorithm. However a complex yield
surface with corners is extremely unwieldy to implement using the Newton-Raphson
iterations. Yield surfaces are most conveniently represented in principal stress space. If
the yield surface is a function of all the three stress invariants, they have an irregular
cross-sectional shape in deviatoric plane and a general trace in the meridian plane when
represented in principal stress space. R. Borja et al [8] have proposed a return mapping
algorithm in the principal stress directions appropriate for three invariant plasticity
models. This reduces the dimension of the stress space for the return mapping from six to
three. However in describing the evolution of the plastic response the rotation of the
principal stress directions must be accounted for. This introduces numerical complexity

in addition to the complexity of the yield surface and flow potential.

41

Since metal foam models are only two invariant elastoplastic models, (the cross—sections
in deviatoric plane are circular) the return mapping can be done in the p-q plane. An

elegant method to do this is described in the next section.

4.5 Algorithm for two invariant plasticity models

Metal foams have a pressure dependant yield surface with circular cross-sections in the

deviatoric plane. The yield surface is not dependant on the third stress invariant(l3) and

in this sense it can be considered a two invariant yield surface. As seen in chapter 3,

metal foam yield surfaces ( f ) and flow potentials(¢) are conveniently represented in

terms of hydrostatic pressure stress (p) and Von-Mises deviatoric stress (q). Aravas [3]
has proposed an algorithm for return mapping in the p-q plane which we shall implement
for metal foams.

We can decompose the plastic strain rate in the directions corresponding to p and q.

Using Eq. 4.4

after time discreetization this equation becomes

1 a.» w
Asp: -A7(_38—6‘7 +— aq nil-J (4.16)
Let us define
-A7— a¢= (4.17)
3p
A-y-al-A = (4.18)
34

42

Then we have

A5,? = §Ae,,§,,- + Aeqnij

(4.19)

Thus the plastic strain increment is decomposed in the spectral directions of the stress

tensor(5,,,n,.). This is possible because the flow rule (¢) is independent of the third

invariant which in turn is a consequence of circular cross-sections of the yield surface in

deviatoric plane.

Now consider the trial stress (Eq. 4.4)

rial _ ,trial
03 tjkl (Eu “'85 )

The updated stress (0', 03,-) can be found from the trial stress and plastic strain increment

0.1.1.: oilflal_ Cij “A85

use the definition of the isotropic elastic tensor

CI!” = 205,16.” —(K_§G)5ij 6k!

where
G- elastic shear modulus
K- elastic bulk modulus

Using Eq (4.21) and Eq (4.19) in Eq (4.20) we get

_ rial _ . _ ..
6,,- _ 0,4,. “3,6,, ZGAsqn,,

use Eq 4.2 to decompose 0,,- and 0'3”“,

-p5,- +- -qn,-j= —pmaI5,-j+ -3—q qmal nt-"al —KA£,, 5,-- —ZGA£q n,--

43

(4.20)

(4.21)

(4.22)

(4.23)

comparing the coefficients in the spectral directions (5,, ,ng.) we get

r.
nf-‘al =n

,, (4.24)

This means that the stress direction in the deviatoric plane does not change during return
mapping. So the return mapping can be done in the p-q plane.

using Eq. (4.23), the updated stress components (p, q) can then be written as

p = p‘”"’ + K A8,, (4.25)

q = gm“ — 3GA£q (4.26)
Eq 4.25 and 4.26 can replace Eq 3.1 to describe the evolution of stress in terms of two
scalars Asp and A8,,

Moreover Eq 4.17 — Eq 4.19 can be used to describe the evolution of the plastic strain in
terms of three scalars(A£,,,A£q,A}/)

The mathematical problem statement in sec 4.3.2 can therefore be now written in terms of
these three scalar quantities (Asp,A£q,Ay). We have further reduced the dimension of

the problem by eliminating A}! from Eq (4.17) and Eq (4.18). The constitutive equations
can now be written as

Find

Pr+Ath+At

subject to

23¢ 8¢
A8,, 35+ A8,, g = (4.27)

f (134.3” .83? ) =0 (Eq 3.7)

such that

p = lam"I + KAep (Eq 4.25)
q = q‘”'“1 — 30Aeq (Eq 4.26)
. As . A8
,. p = ap,trral _ P = Ap,trral q 4 2
a 8 ___—SE a +—§£ ( . 8)
3p 361
35’ = 85"“, + A8,, (4.29)

The dimension of the problem has been effectively reduced from seven to two. This

greatly simplifies the numerical implementation of Newton-Raphson to solve the explicit

equations.

Implementation of the N-R iterations for the above set of equations is done as follows

 

 

Let
X = (4.30)
A3,,
_ A5,, 92 + A8,, 2Q
F(X)= 3‘1 31’ (4.31)
“(10,23,853 _
let
2"“ =1?" +2" (4.32)

th

where 22"“, I?" are the (k +1)” and (k) iterations and E is defined as the correction

VCCIOI'.

E=[c”] (4.33)
C

45

E is the correction for)? vector. By applying Newton-Raphson method to solve

F(}-(')=0 weget

 

3F(fk) _ _k
85(- ~c=-F(X)

wecanwrite

[A11 A12][CP:|=[bl]
A21 A22 cq b2

where

 

[Air 412])“in
A21 A22 —

l2i=[::;’;]

The derivatives of the flow rule and the yield function can be expressed as

3;» 32¢
=— A
A“ aq + 5” aque,

 

8): 82¢
=_. A _—
A” 3p+ 8q8p3A£q+

1421“”(1'i 8f 81?,

82¢

 

E
q apaAep

32¢

 

8
p dquEq

+8f

- 3P 66"” 348p 385’

A”: 8., 331" 0A5,

_3G_a.-.f_+i_a_£_p_

46

(4.34)

(4.35)

(4.36)

(4.37)

(4.38)

(4.39)

(4.40)

(4.41)

where

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

( 2 "1
32¢ 338; M" 32¢ 62¢
= 1+ ‘1 2 K + (4.42)
aque, [ 3,») aqap aqag‘f’
\ a‘1 J
K 2 )_lf 2 \
82¢ 883:” A8” 82¢ 323:? 82¢
= 1— p K — + (4.43)
apaAep [3,)2 apap 92 apes;
i 8A" J K ap /
f 2 “1
32¢ 888:” A8” 02¢
= 1+ p —3G— (4.44)
deA 6,, [ 3 ¢ )2 aqap
I 31" 1
—1
/ 32¢ \ ( 32¢ \
32¢ aqaa-P " 32¢ ages?
= 1- —3G 4.45
aqus,, [ 3(1)]2 aqaq a_¢ ( )
\ a‘1 ) K aq )
33p =—As 1 32¢ (4 46)
BAEP q 22 2 aque,,
34
.p 2
as A 1 3 ¢ (4.47)

3A8}, = 8” _Qfl 2 apaAeq
3p

Eq (4.38) to Eq (4.47) can be calculated for the four material models discussed in

chapter3 to implement the algorithm for these models.

47

4.6 implementation of the algorithm for selected material models
Implementing the above scheme for a particular material model requires the following

derivatives.

Flow rule ¢(p,q,£‘p,8,f’)

EB? 55!
3p ’aq

32¢ 32¢ 82¢
apap' apaq' 3434

 

32¢ 32¢ 62¢ 62¢
apagP ’ aqaep’ apes: ’ aqaa:

 

Yield function f (p.q, 2” . 8,?)

.31 3f.
ap’aq

1 PL
36"" ’ as,”

The derivatives of the four material models described in chapter 3 have been calculated

and are presented in appendix I.

48

Chapter 5

Implementation of the algorithm

5.1 Introduction

The Newton-Raphson closest point projection algorithm described in the last chapter can
be implemented for all the material models of metal foam. FEM software package
ABAQUS has the functionality of a user defined subroutine built in. A subroutine was
written in FORTRAN for implementing the algorithm (Appendix II). The subroutine can
be used for all the four material models by the input of the appropriate derivatives. A
numerical simulation was done in ABAQUS for diagonal crushing of a cube. The
subroutine was implemented for the Deshpande-Fleck model. In this chapter we present

the results of the numerical simulation.

5.2 Input data for the model

The Deshpande-Fleck model requires the experimental data of uniaxial stress and
uniaxial strain as input for incorporating hardening rule. Due to lack of any experimental
data an input curve has been manufactured for closed cell aluminum foam. To ensure that
the manufactured data is close to the actual data, elastic modulus and the initial uniaxial
yield stress were calculated using the scaling equations introduced in chapter 2. (Eq 2.3
and Eq 2.4) The cell wall material properties were taken to be that of Aluminum alloy

2014.
p = 2800 kg/m3

E = 70x109 Mn2

49

0', = 260 MPa

We will consider closed cell aluminum foam with 5% relative density for the simulation.

1: 0.05
p

In order to account for the reduction of observed values of the modulus and yield stress
due to the presence of irregularities in the microstructure (Sec. 2.3), 30% of the equation
values were taken.

We get

5‘ =0.35 GPa
0: =1.18 MPa

For the uniaxial stress strain data we have assumed shape similar to that typical of metal
foams. The plateau stress is taken as 1.25Mpa. Densification strain of 200% is taken

where the material undergoes rapid strain hardening due to the collapse of cells.

 

 

 

 

3000000..
25000004 ,
it"
vzoooooo-
g .
E 1500000 0
m A A A
3 v v v v
1110000001
>.
500000~
0 V I r l
0 05 1 1.5 2 25 3
EOPLSTRAIN

 

 

 

Fig 5.1 typical compressive uniaxial response of closed cell foam

50

5.3 Numerical simulation

Numerical simulation of the diagonal crushing of a cube has been done in ABAQUS. The
motivation for selecting the diagonal crushing process is to provide biaxial stress loading
as against uniaxial loading while at the same time keeping the load condition simple. At
the same time the diagonal crushing process demonstrates the energy absorption of metal
foam due to the reduction in volume.

Solid elements with reduced integration were used to model both the plate and the cube.
Perfectly elastic steel is: used for the material in the plate. Material models using
Deshpande-Fleck and ABAQUS metal foam were used for defining the foam material of
the cube. Simulation was also done for the case of solid aluminum cube in order to
compare the energy absorption capacity of the metal foam.

Two of the bottom faces of the cube are rigidly constrained. Displacement boundary
conditions are applied on the bottom face of the steel plate to come down and again go up
in two steps. The displacement boundary conditions are applied in a ramp fashion to

ensure that there is no discontinuity in the displacements at any time.

L Hm]
I

 

 

Fig 5.2 model used for simulation of diagonal crushing

51

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

<5 3%
387 W

Fig 5.3 diagonal crushing of cube with user defined material model
(Deshpande-Fleck model)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Fig 5.4 diagonal crushing of cube with metal foam model
included with ABAQUS

52

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

    
 

,_ dfilhkz'-

 

 

 

   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

   

Fig 5.5 diagonal crushing of a solid aluminum cube

53

5.4 Comparison of results

The results of the simulation for both the user defined model (Deshpande-Fleck) and the
metal foam model (ABAQUS) are quite similar to each other (Fig 5.3 and Fig 5.4). The
collapse of the elements as the plate comes down can be seen clearly. This demonstrates
the cell wall collapse in actual metal foams and also shows the volume reduction capacity
of the metal foam. The plastic Poissons ratio is close to zero as can be seen by the
negligible flow of material under the plate. When the plate is taken up in step 2 there is

very little elastic springback. This is due to the small elastic region of the metal foams.

It is instuctive to compare the crushing process of foam cube with an aluminum cube (Fig
5.5). The lateral expansion of the cube in this case is a result of the fact the the volume of
aluminum must remain constant during plastic flow. In metal foam most of the internal
strain energy (Fig 5.6) is in the form of plastic strain energy. The elastic strain energy
stored in step 1 is negligible. This means that most of the energy absorbed during step 1
is used up in the plastic volume reduction and very little is stored in the form of elastic
energy. On the other hand, solid aluminum stores more than 50% of its internal energy as
elastic strain energy and the rest is dissipated as plastic strain energy (Fig 5.7). In step 2

all the stored elastic energy is released resulting in springback action.

Due to the different amounts of external work done in the two crushing processes, the
total energy absorbed by the two materials cannot be compared quantitatively. This type
of comparison is useful for making a statement about the the capacity of metal foams to

absorb plastic energy during a large range of deformation.

54

Energy (N-m)

 

0'6 +—+ Elastic Energy
9—9 Internal Energy

1E 5‘ Plastic Energy

0.4

0.2

 

 

 

L L [JELLLJWILLIL LlLLIAJ

' y r v v v v v v v v vvvvvvvv

Time (sec)

Fig 5.6 composition of internal energy in metal foam cube during stepl and 2

Energy (N -m)

 

60.0

40.0

20.0

 

 

 

 

0-0 0.5 1.0 1.5 2.0

Time (see)

Fig 5.7 composition of internal energy in solid aluminum cube during stepl and 2

55

5.5 Conlusion

An elegant method of integrating the constitutive equations of two invariant constitutive
models specific to metal foams is presented. Backward Euler time integration is chosen
for deriving the constitutive equations of the closest point projection algorithm. This
ensures the convergence of the solution even for large strain increments. Newton-
Raphson iterations are used for solving the closest point projection algorithm because of

their global convergence property which makes the algorithm stable.

The other advantage of using Newton-Raphson iterations is that the tangent modulus of
the algorithm can be derived [3, 21]. The algorithmic tangent modulus is essential in
ensuring the optimal rate of global convergence of the FEM equations. Since only quasi-
static loading was used in this work, we have used explicit solution of FEM equations

and so calculation of the tangent modulus is not necessary.

The present algorithm can be extended to include the derivation of the algorithmic

tangent modulus. In tension, metal foams are susceptible to fracture. A further suggestion

for additional work would therefore be to incorporate failure of material in the model.

56

APPENDIX I

Derivatives of the yield function and flow potential for all four models

LS-DYNA material model # 75

 

 

 

 

 

 

 

 

 

¢= qz+—P2
92-2 P
8p 2 9

q2+_2_p2
92: q
361 2 9 2

+—

4 2p
_aifl._fl_l__+§_1 p2
apap 2 2 9 2 2 2

+
4 2p [q2+2p2)2
2
2 2
38¢: l + q 3
44 2 2
+
J4 2P (q2+2p2)2
2
32¢ 9 q 81 p20
_=_ +—
paq 2 3 2 E
2 22 2 22
+ +

[q 2.») [q 2,.)
82¢ _
deé‘p

57

 

 

 

 

 

 

 

 

 

 

aqasP
32¢ _
apasf
82¢ _
aqaef
_l
1 2 2
__ _ 2 __ _
i: p 2(1). p.) {1) p 2(1). p.)
8p a b 02
_l
2 2
-— - 2
8q a b b2
_l
1 2 2
i=1 P—§(Pc-Pt) +(EJZ :2_q_ab
8.?” 2 a b b3 3,312
_l 2
2
1 2 _ -1 _ -1
8f =1 p-2(pc—p') {if 9(” ZOPCJBpC+ 2(" 20”) aa +—2q2 ab
85,? 2 a b 1002 3.9,? a3 38,? b3 38,?
where

58

8b = llpc 2p3+3acpc 00c

 

 

 

 

 

 

38‘: 44/15 2 10 2 ‘2' 38"
(p8 +30'cpc-‘3‘0'c )
\ /
Ba _fldpc
88,? 20 885’
20 2
36 =110'c ”CFC-“5% 3pc
88,? 4% 10 g 88,?
\[Pg +3O'cpc '"9—002)

 

 

/

59

ABAQUS volumetric hardening metal foam model

 

 

 

 

 

 

 

 

¢= qz+-p2
12-2 2
3p 2 2 9 2
+—
2P
it: 4
34 2 9 2
+—
q 217
apaP-Z 2 9 2 2 .3.
+
q 2P [q2+2p2)2
2
340: 2 9 2+ 3
+__ _
\lq 2P [2+2p2J2
2
2
3(1) 9 q 81 pq
—_=_ +—
pdq 2 2 2 9 2
2 22 2 22
+ +
(q 21») (q 2p]
2
89¢ =0
apes?

ABAQUS volumetric hardening metal foam model

 

 

 

 

 

 

 

 

 

¢= qz+-p2
_3_¢ 2 P
3p 2 9
q2+§p2
92: 4
3q 9
42+-p2
pap 2 2 9 2 2 2
+—
4 2p (q2+2p2)2
2
__qu 2 9 2+ 3
q
+—
Jq 2p (2+3p2)2
2
82¢ 9 q 81 p2,,
_=_ +—
pdq 2 2 9 2
2 22 2 22
+ +
[4 2p) [q ,p]
2
3.) =0
apaép

 

 

 

 

 

 

1 2 2 1
3f ram-p.) q 2 127008-12.)
_= + _
3p a [b] 02
_l

2 2
8i: P-§(Pc-Pt) +(_q.)2 fl—
3q 0 b b2
3:...
38”

 

02 a3 b3 38,?

 

61

Deshpande & Fleck Model [10]

¢=—-:——‘/q21+a2p2—0'c

 

 

 

 

 

 

 

.82: 1 02p
6P \ll+(£)2 q2_,_a,2p2
3
it: 1 q
aq 1+(£)2 #2,ng
3
32 1 4 2 1
_T¢= 2a2q2+a2p2 zap 2 2 2
pp “(2) \(q +2212 \lq +6221)
3
fl: -1 2’24
paq a 2 2 a2 2 2
”(3) (q + 2)
32¢ 1 2 2 q 1
——= q +02p
qaq “(2) [ ([q2+a2p2 q2+a2p
3

62

 

 

 

 

 

 

52L: 1 q

aq 1 (ET Jq2+a2p2
3

1:31;

asp 33p

af =0

885’

63

Miller model [18]

 

 

 

 

32¢ __ Zap fl

apes? doacz asp

 

 

 

 

2
3 ¢ =0
aqasf
3f __ 2a
3p - 7+ doO’c

64

 

 

 

 

 

65

APPENDIX II

Fortran Subroutines for implementing algorithm in ABAQUS

Subroutine for user material in FORTRAN for Deshpande-Fleck model

subroutine vumat(
c Read only (unmodifiable)variables -

l nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,

2 stepTime, totalTime, dt, cmname, coorde, charLength,
3 props, density, strainInc, relSpinInc,

4 tempOld, stretchOld, defgradOld, fieldOld,

S stressOld, stateOld, enerIntemOld, enerInelasOld,

6 tempNew, stretchNeW, defgradNew, fieldNew,

c Write only (modifiable) variables -

C

C

7 stressNew, stateNew, enerIntemNew, enerInelasNew )

include 'vaba_param.inc'
character* 80 cmname

dimension props(nprops), density(nblock), coorde(nblock,*),
1 charLength(nblock), strainInc(nblock,ndir+nshr),

2 relSpinInc(nblock,nshr), tempOld(nblock),

3 suetchOld(nblock,ndir+nshr),

4 defgradOld(nblock,ndir+nshr+nshr),

5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),

7 enerInelasOld(nblock), tempNew(nblock),

8 stretchNew(nblock,ndir+nshr),

8 defgradNew(nblock,ndir+nshr+nshr),

9 fieldNew(nblock,nfieldv),

l stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
2 enerIntemNew(nblock), enerInelasNew(nblock)

c-------Local Variables

dimension

1 xx(50),yy(50),

2 C(6,6),

3 strainI(6,l),stressInc(6,l),

4 stresstrial(6,l),devstress(6,l),xntrial(6,l),
5 Amat(2,2),cor(2,l)

real

66

1 lamda,mu,

2 a,fracturestress,fractureenergy,

3 kmod,gmod,

4 alphatrial,xptrial,xqtrial,

5 deltaep,deltaeq,

6 xp,xq,alpha,

7 Rl,RZ,norml,norm2,

8 Y,H,deltaalpha,

9 dfdp,dqu,dfdpdp,dqudq,dfdpdq,dadep,dadeq

integer
1 npoints,mpoints,kstart,
2 i,j,l(
3 nloop
c
c OPEN(UNIT=l7,FIIEz'legr/research/helmholtz/abhishek/debug.dat')
c --------- input properties
i=1
do 5 j=9,nprops,8
xx<i>=props<i>
YY(i)=Pr0PS(i+1)
i=i+1
5 continue
npoints =i-1
mpoints = 4
kstart = l

lamda = props(l)

mu = props(2)

a = Prop-3(3)
fracturestress = props(4)
fractureenergy = props(5)

do 11 i = 1,6
do 10 j = 1,6

10 C(i,j) = 0

11 continue
do 16 i = 1,3
do 15 j = 1,3

15 C(iJ) = lamda

16 continue
C(1,1) = C(l,l) + 2*mu
C(2,2) = C(2,2) + 2*mu
C(3,3) = C(3,3) + 2*mu
C(4,4) = C(4,4) + 2*mu
C(5,5) = C(5,5) + 2*mu

67

C

c..-

C

C

----- loop over nblock material points

C(6,6) = C(6,6) + 2*mu

kmod = (3*lamda + 2*mu)/3
gmod = mu

 

do 100 i = l,nblock

alphatiial = stateOld(i,l)
do 20 j=l,6
strain1(j,1)=strainlnc(ij)

20 continue

call matmul(C,strainI,stressInc,6, 1 ,6)
do 25 j=l,6
stresstiialG,1)=stress0]d(i,j)+stresslnc(j,1)

25 continue

if (stepTime.eq.0.0) then
stateNew(i,l) = stateOld(i,l)
do 27 j = 1,6

stressNew(i,j) = stresstrial(j,l)

27 continue

go to 100
endif

xptiial = -(1 .0/3.O)*(st1esst1ial(l , 1)+stresst1ial(2,1)+
lstresstiial(3,l))

devstress(l,1) = stresstiial(l,l) + xptrial
devstress(2,1) = stresstrial(2,l) + xptrial
devstress(3,1) = stresstiial(3,l) + xptrial
devstress(4,1) = stresstiial(4,l)

devstress(5,l) = stressttial(5,1)

devstness(6,l) = stresstrial(6,l)

xquial = sqrt(l.5)*sqrt(devstress(1,l)**2+devstress(2,l)**2
1+ devstress(3,l)**2+2*devstress(4,l)**2

1+ 2*devstress(5,1)**2 + 2*devstress(6,1)**2)

if (xqt1ial.eq.0.0) then

do 29 j =1,6

xntrialGJ) = 0.0

29 continue

else
do 30 j =1,6
xntrialGJ) = 3.0/(2.O*xqt1ial)*devstress(j,l)

30continue

endif

68

deltaep = O
deltaeq = 0
xp = xptrial
xq = xqtrial
alpha = alphatrial

R1 = 0
if (alpha.eq.0.0) then
=yy(1)
H=(Y)’(2)-yy(l))/XX(2)
else
call startint(xx,yy,npoints,mpoints,
lkstart,alpha)
call polint(xx(kstart),yy(kstart),
lmpoints,alpha,Y,I-I)
endif
R2 = sqrt(l/(l+a**2/9))*sqrt(xq**2+a**2*xp**2) - Y
norml = sqrt(Rl**2)
norm2 = sqrt(R2**2)

if (R2.LE.0.0) go to 35

nloop = 0

do while (norml.GT.l.OR.norm2.GT. 100)

nloop = nloop + 1

dfdp = sqrt(1/(l+a**2/9))*(a**2*xp)/sqrt(xq**2+a**2*xp**2)
dqu = sqrt(1/(1+a**2l9))*(xq)/sqrt(xq**2+a**2*xp**2)
dfdpdp = sqrt(1/(1+a**2/9))*(a**2*sqrt(xq**2+a**2*xp**2)—
1a**4*xp**2/sqrt(xq**2+a**2*xp**2))*(l/(xq**2+a**2*xp**2))
dqudq = sqrt(l/(1+a**2/9))*(sqrt(xq**2+a**2*xp**2) -
lxq**2/sqrt(xq**2+a**2*xp**2))*(l/(xq**2+a**2*xp**2))
dfqu = -sq11(l/(1+a**2/9))*(a**2*xp*qu
lsqrt(xq**2+a**2*xp**2))*(l/(xq**2+a**2*xp**2))

dadep = -l/dfdp + deltaep*(1/dfdp**2)*kmod*dfdpdp

dadeq = —deltaep*(1/dfdp**2)*3*gmod*dfdpdq

Amat(l,l) = dqu+deltaep*kmod*dfdpdq+deltaeq*kmod*dfdpdp
Amat(l,2) = dfdp-deltaeq*3*gmod*dfdpdq-deltaeq*3*gmod*dqudq
Amat(2,1) = kmod*dfdp - H*dadep

Amat(2,2) = -3*gmod*dqu - H*dadeq

cor(1,l) = -R1

cor(2,l) = -R2

call gaussj(Amat,2,2,cor,l,l)

deltaep = deltaep + cor(l,1)
deltaeq = deltaeq + cor(2,l)

69

C

C
C

xp = xpttial + kmod*deltaep

xq = xqtrial - 3*gmod*deltaeq

dfdp = sqrt(1/(1+a**2/9))*(a**2*xp)/sqrt(xq**2+a**2*xp**2)
dqu = sqrt(1/(1+a**2/9))*(xq)/sqrt(xq**2+a**2*xp**2)
alpha = alphatiial - deltaep/dfdp

R1 = deltaep*dqu + deltaeq*dfdp

if (alpha.LT.alphatrial) then

alpha = alphanial

deltaalpha = 0

else

deltaalpha = -deltaep/dfdp

endif
R2=sqrt(l/(1+a**2/9))*sqrt(xq**2+a**2*xp**2)-(Y+H*deltaalpha)
norml = sqrt(R1**2)

norm2 = sqrt(RZ**2)

if (nloop.GT.lOO) pause 'divergent N-R'
enddo

35 stateNew(i,1) = alpha

do 40 j = 1,3
stressNew(i,j) = -xp + (2.0/3.0)*xq*xnt1ial(j,l)

40 continue

do 45 j = 4,6
stressNew(i,j) = (2.0/3.0)*xq*xntrial(j,1)

45 continue

100 continue

101 return

end

70

Subroutine for multiplying matrices

subroutine matmul(A, B, res, p, q, r)
integer p,q,r
dimension A(p.r),B(r,q).res(p.q)

c
do 105 i = 1,p
do 110j = l,q
res(i,j)=0
do 115 k: l,r

res(i,j) = res(i,j) + A(i,k)*B(k,j)
115 continue
110 continue
105 continue
return
end

71

Subroutine for inverting a matrix

subroutine gaussj(a,n,np,b,m,mp)
integer m,mp,n,np,NMAX

real a(np.np),b(np.mp)
parameter (NMAX = 50)

integer i,j,k,l,ll,irow,icol,indxc(NMAX),indxr(NMAX),ipiv(NMAX)

real big,dum,pivinv

do 211 j = l,n
ipiv(j) = O
211 continue
c
do 222 i = 1,n
c
Mg=0
do 213 j = l,n
if (ipiv(j).ne.1) then
do 212 k = 1,n
if (ipiv(k).eq.0) then
if (abs(a(i,k)).ge.big) then
big = abs(a(j,k))
irow = j
icol = k
endif
endif
212 continue
endif
213 continue
ipiv(icol) = ipiv(icol) + 1

if (irow.ne.icol) then
do 214 l =l,n
dum = a(irow,l)
a(irow,l) = a(icol,l)
a(icol,l) = dum

214 continue
do 215 l = 1,m
dum = b(irow,l)
b(irow,l) = b(icol,l)
b(icol,l) = dum

215 continue
endif

72

 

C

C

C

C

indxr(i) = irow
indxc(i) = icol

if (a(icol,icol).eq.0) pause 'singular matrix'
pivinv = 1/a(icol,icol)

a(icol,icol) = 1

do 216 l = 1,n

a(icol,l) = a(icol,l)*pivinv

216 continue

do 217 l = 1,m
b(icol,l) = b(icol,l)*pivinv

217 continue

do 221 ll = 1,n

if (ll.ne.icol) then

dum = a(ll,icol)

a(ll,icol) = 0

do 218 l = 1,n

a(ll,l) = a(ll,l) - a(icol,l)*dum

218 continue

do 2191: 1,m
b(ll,l) = b(ll,l) - b(icol,l)*dum

219 continue

endif

221 continue

222 continue

do 2241: n,1,-1

if (indxr(l).ne.indxc(l)) then
do 223 k = 1,n
dum = a(k,indxr(l))
a(k,indxr(l))= dum

223 continue

endif

224 continue

return
end

73

Subroutine for interpolating the hardness data

subroutine polint(xa,ya,n,x,y,ydash)
integer n,NMAX
real x,y,xa(n),ya(n),ydash
parameter (NMAX=15)
integer i,m,ns
real dy,den,dif,dift,ho,hp,w,c(NMAX),d(NMAX)
real p(NMAX),pdash(NMAX)
ns=l
dif=abs(x-xa(l))
do 231 i=1,n
dift=abs(x-xa(i))
if (dift.lt.dif) then
ns=i
dif=dift
endif
C(i)=ya(i)
d(i)=ya(i)
P(i)=ya(i)
pdash(i)=0
231 continue
y=ya(ns)
ns=ns-l
do 233 m=1,n-1
do 232 i=l,n-m
ho=xa(i)-x
hp=xa(i+m)-x
w=c(i+1)-d(i)
den=ho-hp
if(den.eq.0)pause 'failure in polint'
den=w/den
d(i)=hp*den
c(i)=ho*den
pdash(i)=((x-xa(i+m))*pdash(i)+(xa(i)-x)*pdash(i+l))
1 /(xa(i)-xa(i+m))<I-(p(i)-p(i+l))/(xa(i)-xa(i+m))
pa)=((x-xa(i+m»*p(i)+<xa(i)—x)*p(i+1)>
l I(xa(i)-xa(i+m))
232 continue
if (2*ns.lt.n-m)then
dy=c(ns+1)
else
dy=d(ns)
ns=ns-1
endif

74

y=y+dy

233 continue
ydash=pdash( l )
return
end

 

subroutine for finding start of table
subroutine startint(xa,ya,n,m,kstart,x)
integer n,m,kstart
real x,xa(n),ya(n)
integer i,ns
real dif,dift
ns=l
dif=abs(x-xa(1))
do 241 i=l,n
dift=abs(x-xa(i))
if (dift.lt.dif) then
ns=i
dif=dift
endif
241 continue
kstart=min(max(ns-(m-1)/2,l),n+l-m)
return
end

75

BIBLIOGRAPHY

1. ABAQUS Theory Manual version 6.4, models for crushable foam, section 4.4.6

2. Andrews E, Sanders W and Ginson L, Compressive and tensile behavior of
aluminum foams, Material science and Engineering, vol. A270 113-124 (1999)

3. Aravas N, On the numerical integration of a class of pressure dependant plasticity
models, International Journal for Numerical Methods in Engineering, vol. 24
1395-1416 (1987)

4. Armero F and Perez-Fouget A, On the formulation of closest point projection
algorithm in elastoplasticity part I: the variational structure, International Journal
for Numerical Methods in Engineering, vol. 53 297-329 (2001)

5. Armero F and Perez-Fouget A, On the formulation of closest point projection
algorithm in elastoplasticity part II: Globally convergent schemes, International
Journal for Numerical Methods in Engineering, vol. 53 331-374 (2001)

6. Banhart J, Manufacturing routes for metallic foams, Journal of Manufacturing,
vol. 52(12) 22-27 (2000)

7. Banhart J and Baumeister J, Deformation characteristics of metal foams, Journal
of Material Science, vol. 33 1431-1440 (1998)

8. Borja R, Sama K and Sanz P, On the numerical integration of three invariant
elastoplastic constitutive models, Computer methods in Applied Mechanics and
Engineering, vol. 192 1227-1258 (2003)

9. Clyne T and Markaki A, The effect of cell wall microstructure on the deformation
and fracture of aluminium-based foams, Acta Mettalurgica, Vol. 49(8) 1677-1686
(2001)

10. Davies G and Zhen S, Metallic foams: their production, properties and
applications, Journal of Material Science, vol. 18 1899-1911 (1983)

11. Deshpande V and Fleck N, Isotropic constitutive models for metallic foams,
Journal of Mechanics and Physics of Solids, vol.48 1253-1283 (2000)

12. Doyoyo M and Wierzbicki T, Experimental studies on the yield behavior of
ductile and brittle aluminum foams, International Journal of Plasticity, vol.19
1195-1214 (2003)

13. LS-DYNA theoretical manual, Bilkhu/Dubois foam model, section 16.9

76

14. Gibson L and Ashby M, Cellular solids: structure and properties, Cambrigde
university press (1997)

15. Gradinger R and Rammerstorfer F, On the influence of meso-inhomogeneities on
the crush worthiness of metal foams, Acta Mettalurgica, vol. 47(1) 143-148
(1998)

16. Hanssen A, Hopperstad O, Langseth M and Ilstad H, Validation of constitutive
models applicable to aluminum foams, International Journal of Mechanical
Sciences, vol. 44 359-406 (2000)

17. Kristz B, Fourougi B and Faure K, Behavior of aluminum foam under uniaxial
compression, Material Science and Technology, vol. 16 792-797 (2000)

18. McCullogh K, Fleck N, Ashby M, Uniaxial stress-strain behavior of aluminum
alloy foams, Acta Mettalurgica, vol. 47(8) 2323-2330, (1999)

19. Miller R, A continuum plasticity model for the constitutive and indentation
behavior of foamed metals, International Journal of Mechanical Sciences, vol. 42
729-754 (2000)

20. Numerical Recipes in Fortran 77: The art of scientific computing, Cambridge
University Press, (1992)

21. Ortiz M and Simo J, Analysis of a new class of integration algorithms for
elastoplastic constitutive equations, International Journal of Numerical methods in
Engineering, vol. 21 353-366 (1986)

22. Reyes A, Hopperstad O, Berstad T and Lansgeth M, Constitutive modeling of
aluminum foam including fracture and statistical variation of density, European
Journal of Mechanics, Vol. 22 815-835 (2003)

23. Simo J and Hughes T, Computational Inelasticity, Springer (1998)

24. Sugimura Y, Meyer J and Evans A, On the mechanical performance of closed cell
al alloy foams, Acta Mettalurgica, Vol. 45(12) 5245-5259 (1997)

77

     

    

  

L llllll

““““““ flfljfllflflllw 111111

111]“ 11111 1111 0,7,

1 293