Constitutive models based on dislocation density: formulation and implementation into finite element codes

DOCTORA L T H E S I S DOCTORA L T H E S I S

Luleå University of Technology

Department of Applied Physics and Mechanical Engineering Division of Computer Aided Design

2005:35|: 02-5|: - -- 05⁄35 -- 

2005:35

Constitutive models based on dislocation density

Formulation and implementation into finite element codes

Konstantin Domkin

DOCTORAL THESIS

Constitutive models based on dislocation density

Formulation and

implementation into finite element codes

Konstantin Domkin

September 2005

ii

iii

Abstract

Correct description of the material behaviour is an extra challenge in simulation of the materials processing and manufacturing processes such as metal forming. Material models must account for varying strain, strain rate and temperature, and changing microstructure. Physically based models that can catch the essential phenomena dominating the deformation are expected to have a larger range of usability and validity. Also the models should still be tractable for large scale finite element simulations. This work focuses on the dislocation density models of metal plasticity, their numerical implementation and parameter identification.

The basic concepts of dislocation density modelling are introduced, including the effects of static and dynamic recovery, influence of strain path (PAPER “A”) and modelling of the back-stress (PAPER “B”). Possible mechanisms controlling athermal and thermally activated processes involving dislocations, vacancies and solute atoms are also discussed (PAPER “D”).

Mobile and immobile dislocation densities, vacancy concentrations and other variables are treated as internal state variables. The dislocation models are incorporated in a classical continuum plasticity or viscoplasticity framework by means of the evolution equations for these internal variables which effectively control the hardening behaviour.

Numerical implementation of the models into finite element codes is straightforward. With numerically accurate and computationally efficient standard plasticity stress-update algorithms it is possible to apply the dislocation models in large-scale simulations. In the PAPER “E” a special extended version of a return-map stress-update algorithm and consistent tangent is derived to accommodate the complex coupling effects in the material model. Some of the dislocation models are implemented in user material subroutines for ABAQUS and MSC.MARC and used in simulations of metal forming processes (PAPERS “A” and “C”). Simulations of the cup forming tests shows the capability of the model to capture the main trends in the shape of the forming limit diagrams. Also the models are implemented in a custom toolbox for parameter optimisation in Matlab.

Numerical difficulties of parameter optimisation such as non-uniqueness of the solution, high sensitivity to the starting guess-value and to the choice of the error function, appear to be a common problem with advanced material models (PAPER “B”). Simultaneous curve-fitting of multiple experimental curves of different mechanical testing types is advised to achieve more robust optimisation results. Parameters of dislocation density models usually have clear physical interpretation, and it is possible to obtain values of some of them from sources other than mechanical testing.

The accuracy of physically based models is totally dependent on finding the adequate equations to describe the physical process(es) dominating the material behaviour during deformation.

These equations may be more or less accurate than standard engineering models. However, the use of physically significant parameters connected to the microstructure features such as grain size etc gives a natural way to couple them to models for microstructure evolution which is important in simulations of material processing and manufacturing processes.

Keywords: physically based model, dislocation density, finite element method (FEM),

return-map algorithm, simulation of metal forming, parameter optimisation

iv

v

Preface. Acknowledgements

This thesis summarises my doctoral studies carried out at the Dalarna University (Högskolan Dalarna – HDa), Borlänge, during the years 2000-2005. During this period I was a registered graduate student at the Luleå University of Technology (Ltu), division of Computer Aided Design. The work was supervised by Professor Lars-Erik Lindgren, both at HDa and Ltu.

The financial support for this research project from the Knowledge Foundation (KK-Stiftelsen) is gratefully acknowledged, as well as support from the industrial partners in this project:

Sandvik Materials Technology, Outokumpu Stainless, and SSAB Tunnplåt.

I wish to express my sincere gratitude to my examiner and supervisor, Professor Lars-Erik Lindgren for his guidance throughout my studies, for his support, advice, and critical suggestions to my work, and most importantly for his great patience with me over all these years.

Thank you, Lars-Erik!

The research subject was initially thought of as a revival and continuation of the past work of Professor Yngve Bergström, HDa – the work which I now consider outstanding (and even more so in comparison with my own achievements). I very much appreciate the memorable advice he gave me during the spin-off discussions in the first weeks of my studies.

I would like to thank Dr. Lars Troive, HDa, and Dr. Göran Engberg, MIKRAB, who didn’t hesitate to share their expertise and ideas in our discussions of the research problems, and who also contributed to producing the results of some of the appended papers.

I would also like to send warm thanks to all former and present fellow graduate students and staff at Högskolan Dalarna for making it a very pleasant working environment.

Finally, I’d like to thank my friends and my family, especially my parents for their love, support and encouragement.

Konstantin Domkin

Borlänge, Sweden. September 2005

vi

vii

Contents

1 I

...1

1.1 Research questions and approach... 1

1.2 Limitations and scope of this study ... 2

1.3 Empirical and physically based models ... 3

1.4 Dislocations in crystals ... 4

1.5 Direct and indirect modelling of the crystal microstructure... 4

1.5.1 Discrete dislocation dynamics ... 4

1.5.2 Crystal plasticity... 5

1.5.3 Dislocation density based models... 6

2 C

-

...6

2.1 Yield surface and flow rule... 7

2.2 Internal variables and hardening... 7

2.3 Rate-dependent plasticity (viscoplasticity)... 8

3 S

-

FEM ...8

3.1 Large deformation aspects of constitutive routines ... 9

3.2 Stress-update algorithms for plasticity – general remarks ... 9

3.3 Return-map algorithm. Special cases... 10

4 D

. B

...12

4.1 Deformation mechanisms. Dislocation density... 12

4.2 Mobile and immobile dislocations. Evolution equations ... 13

4.3 Further advances in dislocation density modelling ... 13

5 S

-

...14

5.1 Simplified and semi-empirical models ... 15

5.2 Dislocation models with vacancy concentration... 16

6 O

...17

6.1 Strain-path dependent hardening... 17

6.2 Back-stress and kinematic hardening ... 17

6.3 Anisotropy and texture evolution ... 18

7 P

...18

8 S

’

...19

9 C

...22

10 R

...23

A

viii

1

1 Introduction

Computer simulation of materials processing and manufacturing processes such as metal forming is an important part of product design and development in the industry. Simulations using the finite element method (FEM) are commonplace in design of components. There exist commercial finite element codes that have quite advanced features. Also a number of constitutive models of the plastic behaviour of metals are available for finite element modelling.

Nonetheless, modelling the material behaviour poses an extra challenge in analysis of manufacturing processes. Correct description of material behaviour is crucial for simulations of material deformation processes. Material models must account for varying strain, strain rate and temperature, and sometimes also changing microstructure. Therefore, it is preferable to use physically based models that can catch the essential phenomena dominating the deformation based on the underlying physics of the deformation coupled to microstructure evolution.

However, the models must still be tractable for large scale computations.

The desired features of such models also include the possibilities to obtain model parameters from conventional test data and to implement the model in a finite element code via user-defined material subroutines. Additional advantage of using physically based models is the possibility to rely on information from sources other than mechanical testing to determine the model parameters or verify that their values are in accord with their physical interpretation.

1.1 Research questions and approach The material modelling process includes these steps:

i) obtaining test data,

ii) choosing constitutive model,

iii) numerical implementation of stress calculation algorithm, and iv) finding material parameters.

This work focuses on the physically based material models of metal plasticity, their numerical implementation and parameter identification. The mechanical testing procedures are not in the focus of the study.

Material models based on consideration of the underlying physical processes are expected to have a larger range of usability and validity than more commonly used empirical models. As can be seen in the literature, there has been some progress made to model the evolution of the microstructure for some limited processes but it is not expected that it will be possible to find one universal model covering all aspects of the behaviour of a given material (McDowell 2000).

The research question(s) can be formulated as:

What physically based models for plastic deformations exist? Do they have any advantages compared with more common empirical models?

The approach used in the work consists of the following parts:

x Study of the current state and advances in the field of physically based material models for metal plasticity, especially (but not limited to) dislocation density based models.

x Implementation of some of the dislocation models.

x Parameter identification.

x Study of the performance (applicability, accuracy) of a particular dislocation model by

Bergström (1983), and research the possibilities to improve the model.

2

Accordingly, the thesis content is arranged as follows. This section, Chapter 1, is an extended introduction to physically based material modelling. In Chapter 2, a framework to constitutive modelling of rate-dependent and rate-independent plasticity is presented as well as some common types of hardening. Next, Chapter 3, a numerical implementation of the models of plasticity with internal variables in the context of finite element solution is discussed. In Chapter 4, the basic concepts of dislocation density modelling are introduced and illustrated with examples from the literature. Further aspects of dislocation models are discussed in subsequent chapters: temperature- and strain-rate dependence in Chapter 5, influence of strain path and modelling of back-stress in Chapter 6. The parameter identification procedures are briefly reviewed in Chapter 7. Finally, the overview of appended papers and discussion conclude the thesis.

1.2 Limitations and scope of this study

The physical basis of inelastic deformations in metals is well described by Stouffer & Dame (1996). The book by Reed-Hill & Abbaschian (1992) gives a more fundamental description of the possible dislocation mechanisms for plastic flow. The physical mechanism involved is mainly the movement of dislocations and interaction between dislocations and crystal structure.

This is discussed in more detail in the Chapters 4-6 of this thesis.

The material models intend to capture real material behaviour within a continuum based representation. Thus we deal with different length scales. There are four important length-scales of material modelling, Figure 1-1.

FE analysis

STRUCTURE

Material

properties Dislocation model

TEST SPECIMEN GRAINS DISLOCATIONS

FE analysis

STRUCTURE

Material

properties Dislocation model

TEST SPECIMEN GRAINS DISLOCATIONS

Figure 1-1. Length-scales in material modelling of metals.

(1). Structure (Macro) level: e.g. forming process etc, which are typically analysed using the finite element method (FEM). Our aim here is to provide a material model in terms of macroscopic stresses and strains.

(2). Material test level: e.g. uniaxial tensile test, torsion, etc. Material model should be fitted to this experimental data.

(3). Microstructure level: grain structure of polycrystalline metals, etc. It is usually treated only

in an average sense in the macroscopic models, for example, introducing metallurgical variables

3 into the formulation. To certain extent, this is based on the experimental results of the studies of crystallographic texture, Transmission Electron Microscopy (TEM) and other studies of dislocation cells and patterns formation due to straining are applied at this level.

On the other hand, microstructure problems can also be solved by FEM, using e.g. poly-crystal plasticity theory. This approach will not be used in this study as the focus is on the macro-level material models.

(4). Micro scale (and nano-scale). TEM studies of dislocation systems have provided the experimental evidence to the classical theoretical discussions of physical principles of plastic deformation.

The discrete dislocation dynamics can be applied to directly simulate material behaviour at the micro scale.

To summarise the scope and the limitations of the present work:

The evaluated material models should account for the basic physics of plastic deformation, but they should still be applicable in finite element analysis on the structure/macro level for the simulations of the metal forming processes.

1.3 Empirical and physically based models

Engineering or empirical models are determined by means of fitting model equations and parameters with experimental data without considering the physical processes causing the observed behaviour. These empirical models are also named engineering models as they are more common in engineering applications than the physically based material models.

Physically based material models, on the other hand, are models where knowledge about the underlying physical process, dislocation processes etc, is used to formulate the constitutive equations. Naturally, the division between these kinds of models is somewhat arbitrarily. Both types of models can be considered “engineering”.

There exist a large number of papers dealing with the short-term properties (yield strength, ultimate strength, ductility etc) and long-term properties (creep, fatigue etc). Many references are found in the book by Marshall (1984). In the majority of FE codes the hardening behaviour of metals is represented either by direct interpolation (e.g. piecewise-linear) of the tabulated stress- strain data, or by empirical relationships, e.g. power laws with strain and strain rate. Such equations for the flow stress and plastic strain are not based on any reasoning about the nature of the deformation. Therefore, most of the equation parameters are lacking physical interpretation.

Besides, despite their sometimes remarkably good fit to the measured stress-strain curves within certain range of strains, strain rates and temperatures, empirical relations have no predictive power beyond that range of deformation conditions and material microstructure. Hence, their usage is limited only to the range of deformation conditions at which they were curve-fitted, and the accuracy outside that range is often not satisfactory. The tabulated data approach has the same limitations as interpolation is not possible outside the range of the data.

Two different types of physically based models exist. One option is to explicitly include parameters representing something of the physics of the deformation process. This are computed from evolution equations and are coupled to the constitutive model.

The other possibility is to do it somewhat implicitly: determine the form of the constitutive

equation based on knowledge about the physical process causing the deformation. Although

these equations are laid out as phenomenological, they are based to some extent on physical

considerations and experimental observation, which mainly concerns the strain rate and

4

temperature terms. This approach is a so-called “model-based-phenomenology” (Frost and Ashby 1982).

The same hardening model can be obtained from empirical models and models of dislocation mechanisms as in the case of power law creep.

1.4 Dislocations in crystals

In crystalline materials, the initiation of plastic flow and subsequent permanent deformation under applied stress is related to the presence of certain crystallographic defects – the dislocations (Figure 1-2). The idea of a dislocation was introduced in the continuum mechanics in early 1900s by Volterra and others, and it was purely geometrical. The physical nature of dislocations and their crucial importance in plastic deformation was demonstrated theoretically in 1930s. The presence of dislocations in real crystals was confirmed for the first time in the mid 1950s by direct observations using transmission electron microscopy.

Edge- dislocation

Burger’s vector b

Edge- dislocation

Burger’s vector b

Figure 1-2. A classical sketch of an edge dislocation in a crystal.

There are four major geometric types of crystal defects: point, line, surface, and volume defects.

Point defects represent vacancies and interstitial atoms. Line defects – dislocations and also disclinations – correspond to certain imperfections in the structure of a real crystal. Surface defects are e.g. grain boundaries, twins, internal phase boundaries, stacking faults. Voids and cavities are volume defects.

Geometrically, a dislocation is an isolated singular line in a medium, such that displacement of a point is a multi-valued function on any closed circuit that is threaded by this line. In a metal crystal, the dislocation can be constructed by cutting the crystal along a plane, with a dislocation line being an edge of such plane, and then shifting the adjusting parts of the crystal through the distance equal to inter-atomic spacing. The resulting crystal structure has not changed anywhere except in the vicinity of the dislocation line. The direction of the shift is called the Burgers vector (Figure 1-2).

1.5 Direct and indirect modelling of the crystal microstructure

Physically based models may follow different approaches how to describe the changes of the microstructure and the processes which occur in the crystal.

1.5.1 Discrete dislocation dynamics

One of the direct approaches to modelling is the discrete dislocation dynamics: motion and

interaction of individual dislocation lines is considered, and the stress-strain response of the

material is a result of direct simulation of a huge assemble of dislocations. As pointed out by

5 Estrin (1999), “recent developments in computer modelling make it possible, at least in principle, to simulate the dislocation dynamics on the basis of atomistic calculations, while macroscopic deformation can be modelled using mesoscopic dislocation dynamic simulations.

This scale-bridging approach is still in its infancy. Yet, one can notice a steady progress and some promising results.” Some good examples can be found in the work of Devincre & Kubin (1994), Rhee et al (1998), Zbib et al (2000). However, the algorithms are extremely demanding computationally, and are not readily available for implementation in standard engineering software.

1.5.2 Crystal plasticity

Direct modelling of the crystal structure and deformation of the crystal is the basis for crystal plasticity and its generalised version - the poly-crystal plasticity. It is a very effective approach to deal with material microstructure. It was proposed in the classical work by Asaro & Rice (1977), and since then it has become very popular and proved successful. The material behaviour is studied at the grain level, each grain being a single crystal, and the overall polycrystalline properties are obtained by averaging procedures, e.g. using classical Taylor-Bishop-Hill method.

Recently, the direct FEM modelling of the grain structure has been addressed, but truly realistic models are difficult to create and to solve, and often geometrically simplified models are used.

The crystal plasticity formulation relies on the multiplicative decomposition of the total deformation gradient into inelastic component, associated with pure slip while lattice remains undistorted, and an elastic component, which accounts for elastic stretches and rigid-body rotation.

Figure 1-3. Multiplicative decomposition of the total strain gradient in the crystal plasticity theory (after Asaro and Rice, 1977).

Thermal deformation gradient can also be included as discussed in Meissonnier et al. (2001).

From the kinematics of dislocation motion, the rate of the inelastic deformation gradient tensor is assumed as

p

a

a a a

p

m n F

F ¿ ¾ ½

¯ ®

­ ¦ ^{J } …



where J  represents the shear strain rate on the slip system a, m

and n

are the slip system unit

vectors, defining the slip direction and the slip plane normal.

6

To complete the formulation, the flow and evolutionary equations that describe the behaviour of each individual slip system are required. Various phenomenological models are typically employed to describe self- and latent-hardening of the slip system.

1.5.3 Dislocation density based models

Alternatively it is possible to model the microstructure indirectly, so that the effects of the micro level processes are somehow accounted for, perhaps in an average way, on the macro level. Such type of approach, using the notion of dislocation density, is the subject of the study in this thesis.

Opposite to the discrete dislocation dynamics approach, the dislocation density based models are formulated at the macro level, i.e. all the quantities (dislocation densities, flow stress etc) are calculated for a representative material volume that can be considered homogeneous.

During the 1970´s, in a series of investigations by Swedish researchers, several basic dislocation material models were developed (Bergström 1969, Roberts & Bergström 1973, Bergström 1983).

Similar ideas were independently exploited by Kocks (1966, 1976), and later developed by Estrin & Mecking (1984) and Estrin (1998). Many other researchers incorporated dislocation density logic either implicitly (Bammann 1990) or explicitly (Busso 1998) into material models.

With the above indirect approach, dislocation density models provide a bridge between the micro-level phenomena and macro-level continuum quantities, such as stress and strain. The next chapter introduces the constitutive framework of plasticity, in which the dislocation models should be incorporated.

2 Constitutive framework - plasticity

Within the scope of this study, material models are targeted for application in a finite element analysis of manufacturing or material processes such as sheet metal forming or extrusion. Thus a material model should describe thermo-mechanical deformation in terms of the macroscopic stresses and strains. Several phenomenological classes of constitutive models are well-known in the classical continuum mechanics: elasticity, plasticity, creep, viscosity etc. This study is concerned with elasto-plasticity and visco-plasticity only, as these phenomena are most relevant to the deformation of metals (this will be further discussed in the Chapter 4).

From the continuum mechanics point of view, the constitutive equations should satisfy the second principle of thermo-dynamics (dissipation inequality), and also the objectivity principle (frame indifference). The latter is also important in the context of a numerical solution, which is touched upon in the next chapter. The dissipation inequality imposes additional constraints on material parameters. More in-depth information on this subject can be found e.g. in the book by Lemaitre & Chaboche 1990.

Classical formulations of plasticity are well established since the mile-stone work of Hill (1950).

Originally developed for small strains case, plasticity theory was also generalised to the case of large deformations, for isotropic as well as for anisotropic materials. The formulations can be found in numerous text-books (e.g. Crisfield 1991, Crisfield 1997, Belytschko et al 2000) The rate-independent plasticity theory describes the time-independent irreversible (inelastic) deformations. The domain of validity of pure plasticity is restricted by the following limitations:

moderate temperature usage and non-damaging loads (Lemaitre & Chaboche 1990).

7 The rate-dependent viscoplasticity theory can describe more complex behaviour, with varying strain rates and higher temperatures.

2.1 Yield surface and flow rule

The yield criterion and yield stress define a surface in the stress space that separates the elastic domain inside the surface and inelastic domain outside the surface:

0 ) , (

f ı V (2-1)

In most of the appended papers of this thesis, von Mises yield criterion is used (for isotropic materials). In the PAPER “A” a formulation with an anisotropic yield surface is used.

The flow rule establishes relation between the stress and the plastic strain increments (or rates).

The associated plasticity uses the yield surface function f ( ı ) as a plastic dissipation potential.

This has the advantage of wide applicability to different metals, and can be effectively implemented in computations. The associated flow rule

İ ı w w f d

d

O ^(2-2)

defines that increments of plastic strain are normal to the yield surface. The scalar d is usually O referred to as the ‘plastic multiplier’

In the rate-independent plasticity theory, the plastic multiplier is determined from the solution of the consistency condition:

0

f (2-3)

The above condition together with the yield criterion (2-1) imply that during continuous plastic deformation the stress stays on the yield surface (yield function is zero) and the increments of the yield function are zero too. In the numerical solution the latter property can be used in addition to the consistency condition.

2.2 Internal variables and hardening

If the yield stress V

and other yield criterion parameters depend on plastic strain, one can describe the hardening of material. The kinematic hardening translates the yield surface whereas the isotropic hardening changes the size of the yield surface. Models that include the effect of developing anisotropy will also change the shape of the surface.

Modelling of the isotropic hardening is the main focus of this study, and a review of existing models is presented in the subsequent chapters of the thesis. In the appended PAPER “B”, also kinematic hardening model is adopted for the one-dimensional stress-strain tests. In the literature there exist various formulations to generalise such models to three dimensions.

The plastic multiplier can be treated as an internal variable controlling plastic deformation. It is simply equal to the equivalent plastic strain in the case of von Mises yield criterion.

Using the concept of internal variables, it is also possible to describe more complex hardening

behaviour, when the amount of hardening can not be directly determined by the equivalent

plastic strain alone. For example, the microstructure characteristics of the material can be defined

as internal variables which influence the plastic deformation. In this study, a dislocation density

8

and related parameters are treated as internal variables, thus providing the bridge between the micro-level phenomena (discussed in the Chapters 4, 5 and 6) and macro-level continuum theory.

2.3 Rate-dependent plasticity (viscoplasticity)

The basic idea in rate-dependent plasticity is similar to rate-independent theory. A yield surface is used as a reference. The stress state can be outside this surface, and the plastic strain rate is a function of the distance to the yield surface, the overstress. Creep can be considered as a special case of viscoplasticity without elastic domain.

The magnitude of the effective plastic strain rate is obtained from the flow strength equation:

f

p

g

H ^ (2-4)

where g f is a general notation of the increase in rate when the stress state is outside the yield surface and its value is zero when the stresss state is inside the yield surface. There exists a large range of the choices of functions in different viscoplasticity models. Similar expressions are also used in so called unified models (Miller 1987).

3 Stress-update algorithms in FEM

In this chapter the implementation of material models in a finite element program is discussed.

The three main tasks of a constitutive subroutine for plasticity in an FEM-code:

x Yield criterion calculation (elastic-to-plastic and plastic-to-elastic transition).

x Updating the stress state and internal variables from given total strain increment.

x Calculation of a consistent constitutive matrix.

If the full implicit Newton-Raphson procedure is executed for equilibrium iterations, the consistent constitutive matrix is required in the computation of the tangent stiffness matrix to achieve overall 2nd order rate of convergence. Otherwise various approximations of the tangent matrix can be used. Furthermore, the consistency of the constitutive matrix is not important in the explicit FE solvers where no Newton-Raphson iterations are used.

The choice of either rate-independent plasticity or rate-dependent viscoplasticity framework dictates the choice of a corresponding stress-update algorithm: with the consistency condition for the former model, or with a flow-strength equation for the latter model (e.g. Belytschko et al.

2000, and see PAPER “B” for more references).

In the case of rate-dependent models, a unified-viscoplasticity form of stress-update algorithm could be used (Bammann 1984, 1990). However, it is also possible to formulate the rate- dependent and rate-independent plasticity in parallel, and in the actual implementation to use just one generalised formulation where the former model is a sub-case of the latter (Ponthot 2002).

Similar approach is used in the present study.

Essentially, the rate-independent plasticity framework is chosen, and the rate-effects are

introduced by means of dependency of the hardening parameters on the plastic strain rate

(PAPERS “B” to “D”). Thus we may have a yield limit that is rate-dependent but applied in the

context of rate-independent plasticity.

9 3.1 Large deformation aspects of constitutive routines

Simulation of typical manufacturing processes, e.g. sheet metal forming or extrusion, requires a large deformation analysis. The stress-strain relations can be either of hyperelastic type or hypoelastic type. Hypoelastic approach is based on the choice of an objective stress rate measure defining the increment in stress. The Cauchy stress and the rate of deformation tensor (velocity strain) are commonly used as a conjugate stress-strain pair.

In the current study, implicit solvers of the two commercial finite element codes are employed:

ABAQUS and MSC.MARC. Both of these FE programs provide a set of “user subroutines” for implementation of the user-defined constitutive stress-update algorithms. The internal workings of the above-mentioned software are somewhat different. However they share certain common theoretical background for the solution of the large-deformation problems.

An additive decomposition of strain rate is assumed. The Green-Naghdi objective stress rate is proportional, via Hooke´s law, to the elastic strain rate tensor (Simo & Pister 1984).

Alternatively, Jaumann objective stress rate can be used depending on the FE software, choice of element type and other program options. The hypo-elastic relations are the basis for applying the constitutive model, and the stress computation is performed in a given reference configuration.

The stress computation algorithm will then be form-identical to the small strain case (Hughes 1983). PAPER “E” contains a detailed review of this approach and references to the original articles with rigorous proofs on this matter and for more in-depth explanations.

Accordingly, in the constitutive routines of these implicit FE programs the incremental form of constitutive equations is used with an objective co-rotational stress rate. The numerical algorithm ensures objectivity (frame-indifference) of the constitutive relation by means of a co-rotational transformation:

T n

n

ı R ı R

ı

'

 ' ˜ ˜ ' (3-1)

where 'R is an incremental rotation matrix, appropriately computed.

Hence, the large-rotation related issues of the constitutive routine are taken care of by the ABAQUS or MARC solvers. The user-defined material subroutine will only perform calculations which are form-identical to the small strain case: update stresses, plastic strains and internal variables, and compute consistent tangent matrix:

*

*

į

į ' ı C

˜ ' İ (3-2)

These calculations are made in a reference configuration and the rotation to the final configuration at the end of an increment is made outside the stress-strain algorithm. Therefore, in the subsequent section, the discussion of specific algorithms will use stress and strain relations without distinction of large or small strains case.

3.2 Stress-update algorithms for plasticity – general remarks

The choice of a yield criterion and a flow rule usually demands a certain stress updating algorithm, which in turn defines the derivation of a tangent moduli matrix consistent with the solution procedure (Crisfield 1991, 1997). There exist a number of algorithms and techniques to update the stresses and plastic strains, such as operator splitting, the standard predictor,

“predictor + corrector”, sub-incrementation, pure incremental or forward Euler scheme,

generalised trapezoidal or mid-point algorithms, a backward Euler (return-map), etc.

10

In this work the implicit, backward Euler algorithm was chosen. The method has some desirable properties and has a good performance for large strain increments. Enforcing the consistency condition at the end of the time step will lead to a symmetric consistent constitutive matrix (Simo & Taylor 1986) for many material models. The radial-return algorithm is a special case of backward Euler return-map, and it is especially effective for von Mises criterion: it is indeed a simple ‘radial return’ in deviatoric stress space. It is used in the PAPERS “B” and “C”.

User-defined material subroutine (ABAQUS - UMAT)

Update internal variables Calculate yield stress Calculate hardening modulus

¿ ¾

½

¯ ®

­ 'İ

İ ı

n

U

¿ ¾

½

¯ ®

­

'









İ İ ı

ı į

į

1 1 1

n n n p

n

U

1 n

V

1

U

p Y

n d

H d

H

V Return-map algorithm

INPUT OUTPUT

User-defined material subroutine (ABAQUS - UMAT)

Update internal variables Calculate yield stress Calculate hardening modulus

¿ ¾

½

¯ ®

­ 'İ

İ ı

n

U

¿ ¾

½

¯ ®

­

'









İ İ ı

ı į

į

1 1 1

n n n p

n

U

1 n

V

1

U

p Y

n d

H d

H

V Return-map algorithm

INPUT OUTPUT

Figure 3-1 Dataflow and tasks of a user-defined material subroutine.

The dataflow and the task of a user-defined subroutine is schematically shown on the Figure 3-1 (for the case of isotropic hardening, rate-independent plasticity). The generalisation of the return- map algorithm to the case of anisotropic yield criteria is straightforward and is used in the PAPER “A”.

The computation of the so-called consistent tangent moduli matrix (which corresponds to the implicit iterative equations of full Newton–Raphson method) instead of standard elasto-plastic moduli matrix (which corresponds to the original stress-strain relation in the rate form) is favourable, because the full Newton-Raphson’s solution procedure for the non-linear system of coupled equations in finite element formulation exhibits second order convergence only if the used stiffness matrix is the consistent tangent matrix (Simo 1988; Simo & Hughes 1997).

The importance of the efficiency of constitutive update algorithms within FE analysis workflow is apparent, considering the number of times the routine is called during the solution: every cycle of external load increments, every cycle of equilibrium iterations, in every finite element, in every integration point of an element.

3.3 Return-map algorithm. Special cases

The return-map stress-update algorithm is constructed in two steps. First, the so called trial state is defined as the state of stress when the given total strain increment is taken as elastic and zero plastic strain increment is assumed (elastic predictor step). Second step is the plastic corrector step: the trial stress is relaxed back onto the yield stress in the direction of the closest projection, which is also the direction of the plastic strain increment in the radial-return method.

According to the implicit backward Euler approach the plastic strain increment is computed in

the updated configuration, which depends on the updated stresses to be determined.

11 In the heart of the return-map algorithm is the calculation of the equivalent plastic strain increment (plastic multiplier) using the condition that the stress remains on the yield surface:

) (

)

1

(

1

p p n Y p n

f

V

' H  V H  ' H (3-3)

In the case of von Mises plasticity with radial-return, the expression of effective stress is a linear function of the plastic strain increment, which is used e.g. in the PAPER “E”. Otherwise it may be a non-linear function, e.g. for anisotropic yield criteria used in the PAPER “A”.

In either case, the expression of the yield stress is a non-linear function of plastic strain (the hardening curve). Hence, the whole condition, Eq. (3-3), becomes a non-linear equation, which should be solved for the plastic strain increment. In other words, in the case of non-linear hardening, the solution of the above equation requires scalar iterations.

Thus one may consider possible simplifications or complications of the algorithm depending on the particular hardening law of the material, especially when additional internal variables are used.

(1) The most trivial case, when an analytical expression is available for the yield stress as a function of the plastic strain. If the yield stress is expressed via internal variables,. then the evolution equations for internal variables w.r.t. plastic strain must have an analytical solution over the increment of strain. This kind of model is used in the PAPER “A” to create the user- defined material subroutine UMAT for the ABAQUS implicit solver. The same formulation is used in the custom optimisation toolbox for Matlab in the PAPER “B” to carry out the parameter identification.

(2a) A less trivial case, when there is no analytical expression of the yield stress available, and the evolution equations for the internal variables have to be integrated numerically. In such case, one might still choose the first simplified approach and perform numerical integration for the yield stress in every iteration during the solution of the Eq. (3-3). Such approach was used in the PAPERS “C” and “D”. The advantage is that for uniaxial and proportional loadings it is possible to compute the hardening curve with higher accuracy e.g. using sub-incrementation. However, such method is obviously not the most efficient computationally .

(2b) A numerically more efficient approach to the previous case is simultaneous solution of both the consistency condition and the evolution equations. The effect of the used numerical integration method for the internal variables on the hardening modulus must then be accounted for. The hardening curve is then computed with a reduced accuracy when a lower order update formulas are applied to the evolution equations.

(3) In the PAPER “E”, an even more complicated situation is encountered. The hardening equations are of the form “generalised hardening law”, i.e. a combination of differential and algebraic equations, involving several internal variables and implicitly coupled with the current stress tensor. In addition, elastic moduli and thermal dilatation are coupled with plastic strain.

To achieve better numerical efficiency, an extended radial return method had to be used. The formulation is derived similar to the one found in the book by Belytschko et al. (2000) with certain modifications. The details are given in the PAPER “E”.

As it was already pointed out in the Chapter 2, the material models used in the appended

PAPERS “A” to “D” are physically based models, involving dislocation densities and related

parameters. In the logic of the return-map algorithm these dislocation related variables are

introduced as internal variables to define the hardening behaviour and yield stress. The

fundamentals and derivations of the dislocation models are discussed in the subsequent chapters.

12

4 Dislocation density models. Background

In this chapter, the three mainstream dislocation density based models (Bergström 1983, Gottstein et al 1987, and Estrin et al 1984) are discussed. They all have similar basic assumptions and ideas about dislocation density evolution processes, which basically consist of storage and recovery. But the particular forms of the equations are often different. In the subsequent chapter it will be shown also that such phenomena as dislocation cell substructure, strain-rate and temperature effects, cyclic loading conditions can be accounted for by this kind of dislocation models. The accuracy of one or another model should be validated by further experimental observations.

There are also numerous papers in which similar semi-empirical and dislocation models are discussed, and some suggestions made concerning temperature and strain-rate dependence of the flow stress, and some other micro-structural aspects of the deformation. Those papers are not reviewed here, as all the key issues are already addressed in the presented models.

4.1 Deformation mechanisms. Dislocation density

Dislocations provide the mechanism of plastic deformation of metals. Despite the fact that there exist many elementary geometrical types of dislocations in crystals (e.g. jogs, loops, dipoles etc), it is possible to establish a fundamental relation between macroscopic plasticity effects and the properties of dislocations. The plastic strain is directly related to the motion of dislocations, and hardening/softening of metals is attributed to the interaction of dislocations with each other and with surrounding crystal microstructure.

The microstructure is evolving during the deformation. This is also affected by the strain rate and temperature. The active mechanisms depend on the current structure of the material and applied stress or strain rate at the given temperature. An overview of the dominating deformation mechanisms can be outlined in a deformation map (Frost. & Ashby 1982). The models presented in this thesis concern the following deformation mechanisms: Low-temperature plasticity by dislocation glide (limited by a lattice resistance or Peierl’s stress, by discrete obstacles etc) and Power-law creep by dislocation glide or glide plus climb (limited by glide processes; by lattice- diffusion controlled climb; etc). PAPER “D” contains a brief review of these mechanisms.

Furthermore, dislocations in the crystal can form loops, can pile up on the grain boundaries and precipitate particles, and arrange themselves in various types of cells or substructures called dislocation networks or low energy dislocation structures. These arrangements act as obstacles to the motion of other dislocations, thus providing an important mechanism of hardening. The point, surface and volume defects interact with dislocations and also play important role in hardening mechanisms.

The dislocation density concept links the macroscopic stresses and strains to the underlying micro-structural processes of plastic deformation. A generally accepted assumption for the influence of dislocation density on hardening behaviour is

U V

V

 aGb (4-1)

where U is the so called “forest” dislocation density in cm of dislocations per unit volume, G is

the shear modulus, b is the magnitude of the Burger’s vector, and a is a material constant related

to the crystal and grain structure. The first term, V

is the extrapolated stress to zero dislocation

density (strain-independent “friction stress”), corresponds to the short-range interactions. It is a

friction stress needed to move dislocations through the lattice and to pass short-range obstacles.

13 Thermal vibrations can assist the stress to overcome these obstacles. The second term is due to long-range interactions with the dislocation substructure. It is an athermal stress contribution.

This model does not account for different elementary types of dislocations and dislocation arrangement structures.

4.2 Mobile and immobile dislocations. Evolution equations

The essential property of a dislocation is its mobility, and the simplest dislocation model should distinguish at least two types of dislocations: mobile and immobile. Motion of mobile dislocation carries the plastic strain, and immobile dislocations contribute to the plastic hardening.

The direct simulation of the process, in the spirit of discrete dislocation dynamics, is costly to compute. Hence, an average treatment of dislocation processes is favourable, and the concept of dislocation density is found useful (Estrin 1999).

The dislocation processes include generation of new dislocations, annihilation, immobilisation and re-mobilisation (recovery) of dislocations. It is necessary to have a model for strain hardening of a material due to these processes. Most commonly a kind of evolution equation is derived for each type of dislocation density, which may generally be expressed in the form of

“hardening term – recovery term”:

) ( ) (

 U

U

U ^{d d}

d (4-2)

The increments in the above equation are w.r.t. either time or plastic strain.

Often the following assumptions are also accepted (Bergström 1969, 1983):

(a) mobile dislocation density is strain independent and much smaller than the immobile one.

(b) mobile dislocations move, on average, a distance (mean free path) before they are immobilised or annihilated.

The basis of this modelling is related to the statistical theory of work-hardening by Kocks (1966).

4.3 Further advances in dislocation density modelling

The dislocation density based models can account for static and dynamic strain-aging effects, if extra internal variables are introduced: partially immobilized and locked dislocation densities (Roberts & Bergström 1973). Their evolution laws are similar to presented above.

The dislocation model by Estrin & Mecking (1984) is based on the developments of Kocks (1966, 1976), and is very similar to the Bergström’s model. A one-internal-variable model is usually sufficient for describing monotonic deformation, without abrupt changes of the deformation rate or the deformation path. More general two-internal-variables models were developed by (Estrin 1998) that can account for transients associated with such changes. Then a second, ‘rapid response’ internal variable is added. The model distinguishes between the mobile dislocation density and the density of less mobile (forest) dislocations. The one-internal-variable variant of the model is obtained as a special case of the two-internal-variable one.

The need for a second internal variable also arises when a complex strain hardening behaviour at large strains is to be accounted for. The one-internal- variable model predicts a linear decrease of the strain-hardening coefficient with stress, the stress asymptotically approaching saturation.

This is referred to as stage III hardening. In reality, a new stage of hardening intervenes, stage IV, when strain-hardening rate stays at a nearly constant low level until it drops again in stage V.

The emergence of a dislocation cell structure may be responsible for these hardening features. A

14

single-phased material then effectively becomes a two-phase one, with soft cell interior regions separated by hard cell walls with a much higher dislocation density.

This situation is also assumed in a more sophisticated three-internal-variable model described in Roters, Raabe & Gottstein 2000, which is en extension of a previous basic concept by Gottstein

& Argon 1987. The details of the proposed dislocation model are also discussed in Aretz et al.

2000, and Luce et al. 2001. The model distinguishes three dislocation populations (Figure 4-4):

x mobile dislocations travelling through cell structure, U

x the immobile dislocation density inside cells, U

x the density of immobile dislocations in the cell walls, U

.

Figure 4-4. Schematic drawing of the arrangement of the three dislocation classes (after Roters, Raabe & Gottstein 2000).

The mean free path of mobile dislocations is determined by the effective grain size and three obstacle spacing: the forest dislocation spacing in the cell walls, in the cell interior, and the spacing of the precipitates. The effective shear stresses are calculated separately for the inside of the cell and inside of the cell walls. The macroscopic flow stress is an average of the effective stress in cells and cell walls.

5 Strain-rate and temperature effects on hardening

The current study (appended PAPERS “A” to “E”) uses the concept of rate-independent plasticity and a yield surface. It is an approximation where the plastic strains are assumed to develop instantaneously in such a way that the stress state stays on the yield surface during a plastic process. The applied stress is then equal to the yield stress during a plastic deformation.

On the other hand, the physical processes that determine this yield surface are rate-dependent.

Therefore, it leads to a class of models where a rate-dependent yield limit is introduced in the context of rate-independent plasticity.

The thermal effects are introduced by means of the temperature-dependence of the yield limit

and other hardening parameters. (In one of the appended papers, PAPER “E”, elastic moduli and

thermal dilatation are temperature-dependent properties too.) Furthermore, it is assumed that the

deformation process does not affect the temperature and also the thermal conductivity

phenomena are left outside the scope of the material models. This will be the case e.g. in a

mechanical step of a staggered solution of coupled thermo-mechanical problems, with a thermal

15 step followed by a mechanical step. Then the temperature is simply prescribed in the mechanical stage. However, the models themselves can be generalised to include any form of coupling between temperature and deformation, e.g. by means of appropriate equation for heat generated due to plastic dissipation.

The physical basis of the strain-rate and temperature dependence in the dislocation models is discussed in this chapter. There exists a fundamental relation between plastic strain rate and average dislocation velocity, the Orowan equation:

 ^{ı,ı ,ȡ , İ ,...} 

v v , ȡ v

İ 

v



(5-1)

where v is a function of stresses, dislocation densities, strain rate, temperature etc. Most

dislocation models accept this relation to introduce rate-dependency. However, there is no agreement in particular choices of the velocity function.

5.1 Simplified and semi-empirical models

First, an extension of the basic dislocation model to account for rate-dependent yield limit (Bergström & Hallen 1982) is described. See also the PAPERS “A” and “B” for further explanations of the symbols and parameters. The dislocation multiplication process is assumed to be rate-independent. The recovery process is strongly dependent on temperature and strain rate.

The rocevery term of the equaton (4-2) is assumed to be proportional to the dislocation density, ȍ U H

U 

^ (5-2)

The expression for recovery parameter, : , consists of an athermal component :

and a thermal component:

1

0

( )

,



¸ ¸

¹

·

¨ ¨

©

 §

ref p

ref

p

T ȍ ȍ T

ȍ H

H H



  (5-3)

The formula is derived on the basis of the following assumption: the re-mobilisation of immobile dislocations takes place predominantly in the cell walls by means of diffusion controlled vacancy climb.

Dislocations in motion experience resistance in their glide plane, and certain stress should be applied in order for a dislocation to start moving. The expression for the friction stress relating these effects is assumed to be

) (

,

T m

ref p

ref p

f

T T

¸ ¸

¹

·

¨ ¨

©

 § H V H V H

V 

  , V

  T V

 V

 V

(5-4)

where V

, V

, V

are the solution, precipitation and grain-size hardening components (Bergström 1983). As temperature rises, the strain rate dependent part of the friction stress quickly drops to zero.

The above two equations are used in the material model in the PAPER “B”.

The relation in Eq. (5-4) is equivalent to the kinetic equation of viscoplasticity (Hasegawa et al.

2000):

n i y m

p

v v ¸¸

¹

·

¨¨ ©

§  v v

0

) ˆ (

, V

U V U V

H ^{ (5-5)}

which could be used in a unified-viscoplastic form of stress-update algorithm.

16

Estrin (1998) uses unified viscoplasticity approach, but with a different formula:

p m m

p

1

0

0

¸ ¸

¹

·

¨ ¨

©

˜ §

¸ œ

¹

¨ ·

©

§

H V H V V

H V

H 

 

 (5-6)

Thus, a multiplicative split of strain and strain-rate terms is assumed, which is often replaced by an additive split in other dislocation models in order to fit experimental data better, and based on some physical reasoning.

According to Gottstein (1987), their dislocation model (see section 4.3) can be further enhanced to include rate- and temperature-dependency. The effective shear stresses in cells and walls are derived from the kinetic equation of state for this model – the Orowan equation J   H M U

bv , where v is the glide velocity according to the following formula,

¸¸ ¹

·

¨¨ ©

§

¸¸ ¹

·

¨¨ ©

§ 

T k

V T

k f Q v

B eff

B

O W

~ sinh

exp (5-7)

where O is the jump width, i.e. the mean spacing of obstacles (the immobile forest dislocation in this case), f the attack frequency (assumed to be constant), Q activation energy for dislocation glide (or “forest cutting”), and V is the activation volume.

5.2 Dislocation models with vacancy concentration

A different and somewhat more systematic approach is used in the PAPER “D” to introduce various effects, including strain-rate and temperature. The ideas are collected from the literature devoted to physical modeling of the plastic deformation, and equations are organised in sub- models responsible for particular effect or process. An abridged and modified version of this model is also used in the PAPER “C”.

First, in the Orowan equation, Eq. (5-1), the velocity is related to the time it takes for a dislocation to pass an obstacle, the waiting time t

, as the time of flight to the next obstacle is negligible. The velocity is written as

kT G a kT G

a

e e

b

v E X

/ X

^(5-7)

where / E b is the mean free path between two successful events. X

is the attempt frequency and depends on the obstacles. G ' is the activation energy, k is the Boltzmann constant and T is the temperature in Kelvin. Considering possible mechanisms controlling the dislocation glide in a crystal structure and using a general formula for the activation energy, the expression of the friction stress as a function of the effective plastic strain rate is obtained,

q p

p ref

F kT

1 1

ln 1

ˆ

*

¸ ¸

¸

¹

·

¨ ¨

¨

©

§

¸ ¸

¹

·

¨ ¨

©

§

¸ ¸

¹

·

¨ ¨

©

§

 '

H W H

V ^ _ ^(5-8)

where ' is the free energy required to overcome the lattice resistance or obstacles without aid F from external stress. The quantity W ˆ is the athermal flow strength that must be exceeded in order to move the dislocation across the barrier without aid of thermal energy. It reflects not only the strength but also the density and arrangement of the obstacles.

Thus, the explicit rate- and temperature-dependence of the yield stress is introduced in this short-

range stress component.

17 The “hidden” rate- and temperature-dependence of the yield stress is also introduced: in the long-range stress component via dislocation density evolution equation. Following the logic of Militzer et al. (1994) with some modifications, an equation for the dislocation recovery is obtained,

²

v v v

i

kT

Gb c D c

c U U

U 

 (5-9)

where c

is the vacancy concentration, and the meaning of other parameters is explained in the PAPER “D”. Furthermore, an evolution equation for the vacancy concentration is used, coupled with dislocation density and other parameters. Together, they render the process of dislocation recovery as being controlled by temperature and strain-rate. See the PAPER “D” for detailed derivation and an overview of the literature on the subject.

6 Other aspects of material modelling

In this chapter modelling of certain aspects of material behaviour is addressed, which typically may be encountered in the process of sheet metal forming such as deep drawing, and in cyclic tests.

6.1 Strain-path dependent hardening

The formability of sheet metals (characterised by forming limit diagrams, FLD) is observed to depend strongly on the strain path during the deformation. Besides, experiments with rolled steel sheets show (Laukonis & Ghosh 1978) that the effective stress–equivalent strain curves are different in uniaxial and in balanced biaxial tension. Another noteworthy feature of the FLD is that the decrease of failure strains after a pre-strain is material dependent. All this suggests that the material model should account for such behaviour.

A plausible explanation was suggested by Bergström & Ölund (1982): the hardening of the material is dependent on the strain path. Their theoretical treatment is adopted in the PAPER “A”

to compute forming limit diagrams. However, the physical interpretation of some of the resulting equation parameters is not straightforward.

In the numerical implementation, the strain path is defined by the ratio of principal strains, and all hardening parameters are assumed to depend on this ratio. It is quite similar to the logic of the rate-dependent hardening discussed in the Chapter 3. Effectively, material model has isotropic hardening, but the hardening parameters are evaluated based on a strain path parameter.

Thus, a kind of strain path dependent isotropic hardening is obtained.

6.2 Back-stress and kinematic hardening

The most prominent factor in the cyclic tension-compression tests is the presence of kinematic hardening associated with the back-stress. In the appended PAPER “B” an advanced dislocation model by Estrin et al. (1996) is used to accommodate the evolution of the back-stress internal variable. The back-stress is attributed to the formation of channel-like dislocation structures upon cyclic deformation, with the immobile dislocation density in channel walls consisting of

‘trapped’ dislocations and of ‘recoverable’ dislocations. The changing curvature of the dislocation segments in the channels is believed to be the driving force of the process.

Accordingly, an evolution equation for this quantity is chosen in the following form: