Novel constitutive models for capturing anisotropy, tension-compression asymmetry, and anisotropic hardening of titanium

(1)

Research Collection

Doctoral Thesis

Novel constitutive models for capturing anisotropy,

tension-compression asymmetry, and anisotropic hardening of titanium

Author(s): Raemy, Christian Publication Date: 2018 Permanent Link: https://doi.org/10.3929/ethz-b-000253718 Rights / License:

In Copyright - Non-Commercial Use Permitted

This page was generated automatically upon download from the ETH Zurich Research Collection. For more information please consult the Terms of use.

(2)

DISS. ETH NO. 24864

Novel constitutive models for capturing anisotropy,

tension-compression asymmetry, and anisotropic

hardening of titanium

A thesis submitted to attain the degree of DOCTOR OF SCIENCES of ETH ZURICH

(Dr. sc. ETH Zurich)

presented by Christian Raemy

MSc ETH in Mechanical Engineering born on 27 December 1987

citizen of Plaeien FR

accepted on the recommendation of Prof. Dr. Pavel Hora, examiner

Prof. Dr. Ir. Ton van den Boogaard, co-examiner

(3)

(4)

Acknowledgements

This thesis was written during my time as a doctoral student at the Institute of Virtual Manufacturing (IVP) of ETH Zurich. I thank the head of the institute and my doctoral adviser, Prof. Dr. Pavel Hora for the great opportunity to work on this project and for his interest and encouragement. Many thanks also go to Prof. Dr. Ir. Ton van den Boogard for accepting the co-examination and the fruitful discussions we had.

Further, I want to thank my colleagues at IVP for their steady support. Especially the help of Dr. Niko Manopulo was highly appreciated, he was a constant source of valuable hints during the development of the models and their FE-implementation. Also Dr. Bekim Berisha and Dr. Christoph Becker always oered helpful inputs.

Last but not least I thank my family for the many years of support during my studies.

Zurich, March 2018 Christian Raemy

(5)

(6)

Abstract

Numerical simulations are state of the art for process planning, especially in the forming industry. In order to design tools and auxiliary devices such that the rst produced work piece complies with the technical specications, an accurate prediction of the nal shape and the forces is required. Exact simu-lation results, however, rely heavily on suciently accurate material models. Titanium has an increasing scope of application because of its extraordinarily advantageous properties but standard models fall short of describing its com-plex plastic behaviour. Titanium and its alloys do not only show anisotropy, its yield strength also depends on the sign of the applied load (tension-compression asymmetry). Furthermore, the hardening behaviour is greatly aected by twin formation, a mechanism which was found to be strongly direction dependent. As a mathematical description of these phenomena, two novel constitutive laws will be presented. One is a ow function on the basis of a two-dimensional Fourier series which allows for the accurate modelling of anisotropy as well as tension-compression asymmetry. Almost arbitrarily complex ow surfaces can be constructed without altering its formal structure. Through this, numerous experimental data (ow stresses and Lankford parameters) can be captured exactly. It will be shown that the criterion can be also applied to other material classes with ease.

Crystal twins form new grain boundaries, which impose an obstacle for dislocation slip and, therefore, lead to the hardening of the material. The second model is a modied Hall-Petch relation which accounts for this eect. An internal variable is introduced to capture the twin volume fraction and it is updated depending on the deformation history of the material.

Both models are implemented in commercial nite element software. Vali-dation experiments show without increase of calculation time a signicant improvement of result accuracy compared to state of the art models.

(7)

(8)

Zusammenfassung

Numerische Simulationen sind aus der industriellen Produktionsplanung speziell in der Umformtechnik nicht mehr wegzudenken. Nur mit einer genauen Voraussage der fertigen Bauteilgeometrie und der benötigten Kräfte lassen sich die Werkzeuge und Hilfsvorrichtungen so konstruieren, dass sie ohne Nacharbeiten produktiv eingesetzt werden können, wodurch Kosten eingespart werden. Genaue Simulationsergebnisse setzen aber hinreichend genaue Mate-rialmodelle voraus. Gerade bei Titan, das aufgrund seiner günstigen Eigen-schaften zunehmende Anwendung in vielen Industriezweigen ndet, scheitern die verbreiteten Modelle jedoch an der Komplexität der Materialeigenschaften. Titanwerkstoe weisen nämlich nicht nur senkrechte und planare Anisotropie auf, die Festigkeit hängt auch vom Vorzeichen der anliegenden Belastung ab (Zug-Druck-Asymmetrie). Zudem wird das Verfestigungsverhalten stark durch Kristallzwillingsbildung beeinusst, dessen Auftreten jedoch ebenfalls vom Spannungszustand abhängt.

Zur Beschreibung dieser Phänomene werden zwei neuartige Materialmodelle vorgeschlagen. Erstens ist dies eine Fliessfunktion auf der Grundlage einer zweidimensionalen Fourierreihe. Mit ihrer Hilfe lassen sich sowohl Anisotropie wie auch Asymmetrie abbilden. Ohne Veränderung ihrer formalen Struk-tur können praktisch beliebig komplexe Fliessächen konstruiert werden, so-dass auch eine grosse Zahl von Messdaten (Fliessspannungen und Lankford-Parameter) exakt wiedergegeben werden können. Es wird gezeigt, dass dieses Fliesskriterium ohne Weiteres auch auf andere Materialklassen angewendet werden kann.

Beim zweiten Modell handelt es sich um ein modiziertes Hall-Petch-Gesetz, womit der Verfestigung durch Zwillingsbildung Rechnung getragen wird. Die Kornzwillinge bilden nämlich neue Korngrenzen, die die Verformung des Mate-rials erschweren und somit dessen Festigkeit erhöhen. Der Anteil der Zwillinge wird als innere Variable behandelt, die abhängig von der Belastungsgeschichte verändert wird.

Beide Modelle wurden in kommerzielle FE-Software integriert. Validierung-sexperimente zeigten ohne Erhöhung der Rechenzeit eine signikante Verbesserung der Vorhersagegenauigkeit gegenüber den bekannten Modellen.

(9)

(10)

5.2.1 Isothermal model . . . 78 5.2.2 Non-isothermal extension . . . 79 5.2.3 Application to titanium . . . 81 5.3 Conclusions . . . 85 6 Numeric implementation 87 6.1 Stress integration . . . 87 6.1.1 Cutting-plane algorithm . . . 87 6.1.2 Closest-point-projection algorithm . . . 90 6.2 Continuum tangent . . . 93 6.3 Thermo-mechanical coupling . . . 95 7 Validation experiments 97 7.1 Deep drawing experiment . . . 97

7.1.1 Set-up . . . 97

(12)

Contents xi 7.1.3 Results . . . 99 7.2 Die bending . . . 105 7.2.1 Set-up . . . 105 7.2.2 FE-Modell . . . 106 7.2.3 Results . . . 108 7.3 Sheet-bulk forming . . . 110 7.3.1 Set-up . . . 111 7.3.2 FE model . . . 111 7.3.3 Results . . . 112 7.4 Conclusions . . . 112 8 Concluding remarks 115 Bibliography 117 A Specimen geometries 125 B Jacobian between spherical and Cartesian coordinates 127 C Hessian operator in polar coordinates 131 D Gaussian curvature 133 E Hessian of FAY 135 F Calibrated model parameters 139 F.1 Titanium . . . 139

(13)

(14)

List of Figures

1.1 Measured material properties and model predictions. . . 2 1.2 Measured and simulated ow curves for uni-axial compression

in rolling and transverse direction. . . 3 2.1 Bravais lattice and atomic arrangement of an hcp- unit cell. . 9 2.2 Slip systems of the titanium lattice. . . 10 2.3 Schematic view of the crystal lattice after twinning. . . 10 3.1 Experimental iso-thermal ow curves at strain rate 0.5 s−1

(Ti-CP Grade 4). Tensile ow curves measured on deformation dilatometer. . . 25 3.2 Quasi-static ow curves in dierent directions at room

tem-perature (Ti-CP Grade 4). Tensile ow curves measured by standard tensile tests. . . 26 3.3 Compressive ow curves in TD at dierent temperatures and

strain rates (Ti-CP Grade 4). . . 27 3.4 Optical micrographs (polarized light) after uni-axial

compres-sion, |log (ε)| = 15 %. . . 28 4.1 Flow surface after von Mises in space of plane stresses . . . . 34 4.2 Embedding of the canonical not used spherical coordinates

in space of plane stress states. . . 36 4.3 Embedding of the used spherical coordinates in space of plane

stresses . . . 37 4.4 Yield loci after Hosford (1972) for dierent values of the

exponent a. For a = 2 the plain-strain points are depicted. . . 52 4.5 Calibration of the proposed ow criterion . . . 54 4.6 Yield loci predicted by the presented models τxy = 0 (Ti-CP

Grade 4). . . 65 4.7 Uni-axial ow stresses and R values (Ti-CP Grade 4). . . 66 4.8 Predicted yield loci, uni-axial ow stresse and R values after

(15)

xiv List of Figures

4.9 Predicted yield loci, uni-axial ow stresse and R values after Park and Chung (2012) and FAY for AA5042. . . 70 4.10 Predicted yield loci, uni-axial ow stresse and R values after

Stoughton and Yoon (2004) and FAY for AA2008-T4. . . 72 4.11 Predicted yield loci, uni-axial ow stresse and R values after

Stoughton and Yoon (2004) and FAY for AA2090-T3. . . 73 5.1 Geometric interpretation of the angle ω . . . 79 5.2 Qualitative course of the function ky(T )and graphic

interpre-tation of its parameters. . . 80 5.3 Simulated ow curves (FAY combined with the twinning

model) for uni-axial compression at dierent angles to rolling direction (isothermal, strain rate independent model at room temperature, Ti-CP Grade 4) . . . 82 5.4 Uni-axial compressive ow curves, ˙¯ε = 0.5 s−1_{, Ti-CP Grade 4.}

Simulation with FAY and non-isothermal twinning model. . . 83 5.5 Eect of temperature (T = 20◦_{C, 200}◦_{C, 400}◦_C_{and 500}◦_C)

on uni-axial compressive ow curves, ˙¯ε = 0.5 s−1_{, Ti-CP Grade}

4. Simulation with FAY and non-isothermal twinning model. 84 5.6 Eect of temperature (T = 200◦_{C, 400}◦_C_{and 500}◦_{C) on}

uni-axial tensile ow curves, ˙¯ε = 0.5 s−1_{, Ti-CP Grade 4.}

Simula-tion with FAY and non-isothermal twinning model. . . 85 6.1 Visualization of the Cutting-Plane Algorithm. . . 89 6.2 Visualization of the Closest-point-projection algorithm. . . . 90 6.3 Visualization of the outer iteration scheme of the multi-stage

stress integration. . . 93 6.4 Graphic interpretation of the factor w (6.25), which is needed

at the elasctic-plastic transition. . . 95 7.1 Drawn cup. . . 97 7.2 Shape of the nal cup (simulated with FAY) and used mesh. 98 7.3 Optically measured and predicted ange radius r (θ) . . . 100 7.4 Measured and predicted thickness strain at outer ange radius 101 7.5 Measured and predicted punch forces . . . 103 7.6 Predicted volume fraction of twinned grains (simulation with

FAY). . . 104 7.7 Calculation time of the explicit forming step, normalized by

(16)

List of Figures xv

7.8 Set-up of the die bending simulation (bending operation). . . 105

7.9 Meshed specimen of the die bending experiment. Symmetry was considered. . . 105

7.10 Cross section of the punch (die bending) . . . 106

7.11 Set-up of the die bending simulation (trimming operation). . 106

7.12 Predicted yield loci, ows stresses and R values of Ti-CP Grade 4 according to Hill48 . . . 107

7.13 Approximation of the tensile ow curve in RD with a Hockett-Sherby function . . . 108

7.14 Measured and calculated cross-sections of the work piece after springback. . . 109

7.15 Projection to the shear-free plane of the stress state (normal-ized with the corresponding von Mises stress) of all plastically deformed elements in the die bending simulation. Increment at the end of the rst bending. . . 110

7.16 Simulation set-up for the bulk forming experiment . . . 111

7.17 Final work piece (bulk forming). . . 112

7.18 Meshed specimen (bulk forming). . . 112

7.19 Measured and predicted maximal punch forces during the bulk forming experiment. . . 113

(17)

(18)

List of Tables

0.1 Greek letters . . . xix 0.2 Latin letters . . . xx 0.3 Diacritic . . . xxii 3.1 Allowed impurities (alloying elements) in commercially pure

titanium (Ti-CP Grade 4), weight per cent. . . 23 4.1 Flow stresses and Lankford parameters of Ti-CP Grade 4 for

dierent experiments at room temperature. . . 64 5.1 Hockett-Sherby parameters for uni-axial compression at room

temperature (Ti-CP Grade 4). . . 76 5.2 Parameters of the Spittel and Spittel (2009) model for

uni-axial compression in RD (Ti-CP Grade 4). . . 77 5.3 Optimized parameters of the isothermal twinning model

(Ti-CP Grade 4). . . 81 5.4 Optimized parameters of the non-isothermal twinning model

(Ti-CP Grade 4). . . 82 F.1 Coecients for Cazacu et al. 2006, a = 2, Ti-CP Grade 4. . 139 F.2 Coecients for Plunkett et al. 2008, a = 4, Ti-CP Grade 4. 139 F.3 Coecients for FAY, q = 2, Ti-CP Grade 4. . . 140 F.4 Coecients for FAY, q = 6, AA209-T3 (Park and Chung 2012). 141 F.5 New coecients for Yld2000, a = 8, AA5042, associated model

(Park and Chung 2012). . . 142 F.6 Coecients for FAY, q = 6, AA5042 (Park and Chung 2012). 142 F.7 Coecients for FAY, q = 6, AA2008-T4 (Stoughton and

Yoon 2004). . . 143 F.8 Coecients for FAY, q = 6, AA209-T3 (Stoughton and

(19)

(20)

Nomenclature

Table 0.1: Greek letters

Symbol Meaning(s)

α Material parameter (twinning model) γ Shear strain

δij Kronecker-Delta

ε Linear strain tensor

η Material parameter (twinning model) ζ Virtual temperature

Θ Heaviside function

θ Angle between test and rolling direction λ Thermal conductivity

µ Friction coecient ν Poisson's ratio

ξ Virtual displacement ρ Density

Σi Stress in reduced space (FAY)

σ Cauchy stress tensor σ Normal stress σy Flows stress

σTW Material parameter (twinning model)

τ Shear stress

Φ Auxiliary quantity (Cazacu et al. 2006) ϕ Polar angle (FAY)

ψ Azimuth angle (FAY)

ω Angle to uni-axial stress in transverse direction Ω Domain of integration

(21)

xx Nomenclature

Table 0.2: Latin letters

Symbol Meaning(s) A Surface

Material parameter (Hockett-Sherby model) am,n Coecient of Fourier series

a Edge length of a hcp cell Exponent (Cazacu et al. 2006) B Displacement-strain matrix

B Material parameter (Hockett-Sherby model) b Volumetric forces

bm,n Coecient of Fourier series

C Stiness tensor

C Transformation tensor (Cazacu et al. 2006) Cλ Heat capacity matrix

cm,n Coecient of Fourier series

c Height of a hcp cell Heat capacity

D Rate of deformation tensor D Average grain diameter dm,n Coecient of Fourier series

E Material tangent (Closest-Point-Projection) F Nodal force vector

F Material parameter (Hill 48) Auxiliary quantity (FAY)

Auxiliary quantity (stress integration) f Fourier series

G Material parameter (Hill 48)

Derivative of ow stress wrt. to strain rate g Auxiliary quantity (twinning model) H Interpolation matrix

Hesse-Matrix

H Material parameter (Hill 48)

Derivative of ow stress wrt. to strain h Heat-transfer coecient

K Stiness matrix Kλ Thermal stiness

K Gaussian curvature

(22)

Nomenclature xxi

Table 0.2: Latin letters (continued)

Symbol Meaning(s)

k Material parameter (Cazacu et al. 2006) Iteration counter

ky Material parameter (twinning model)

M Mass matrix

m Material parameter (Hockett-Sherby model) mi Material parameter (Spittel-Spittel model)

n Outward normal Eigen vector

N Material parameter (Hill 48) n Time step

Material parameter (Hockett-Sherby model) Qλ Sum of heat introduced to the element

Q Heat

q Source term (heat conduction) Exponent (FAY)

R Residuum Rotation matrix R Lankford parameter

r Radial coordinate (FAY) S Transformed stress deviator

Si Material parameter (twinning model)

s Stress deviator u Displacement T Temperature

T0 Material parameter (twinning model)

t Stress vector t Time V Volume

Volume fraction of twins w Weighting factor

x Position vector

(23)

xxii Nomenclature Table 0.3: Diacritic Symbol Meaning ¯ x Equivalent quantity ˆ

x Nodal quantity (FEM)

˙x Total time derivative of a quantity ˇ

Xi i-th eigenvalue of a tensor

Xe Elastic part of a tensor Xp Plastic part of a tensor

(24)

1 Introduction

1.1 Motivation

Because of their excellent strength-to-density ratio, corrosion resistance and biocompatibility, titanium alloys are used more frequently in aerospace and the medical industry, as an example. Lightweight construction is a necessity to improve fuel eciency, which leads also to more widespread application in the automotive industry. Accurate simulations of manufacturing processes are a prerequisite of an ecient production. However, they rely on constitutive mod-els which are able to capture the material behaviour accurately. For titanium, specically, this behaviour is rather complex and currently available models may sometimes only deliver unsatisfactory approximations. In particular, they must be able to describe anisotropy as well as tension-compression asymmetry (strength-dierential). Simple ow functions do not oer sucient degrees of freedom to model both ow stresses and ow directions accurately, while more sophisticated functions cause unacceptably high computational costs.

Therefore, a completely new framework for the construction of ow surfaces of virtually arbitrary complexity will be presented. The basic idea is the parametrization of the stress tensor by spherical coordinates and the expression of the ow function as a Fourier series of the angular coordinates. Experimental data are used to dene support functions, which span the ow surface. In-between, an interpolation is carried out, respecting physical constraints.

In this theoretical context, both symmetric as well as non-symmetric ow surfaces can be constructed, which although using an associated ow rule are able to reproduce even a high number of experimental data. Furthermore, the computational eciency is at least comparable to state-of-the-art models. Fig. 1.1 shows experimental data and model predictions of a widely used model (Cazacu et al. 2006) and the proposed function. A signicant improvement of accuracy is observable specically in the R value plot.

Further, a physically based model capturing directional hardening of tita-nium will be proposed. It has been recognised that crystal twinning is the main cause of this phenomenon. It is assumed that the twins create new grain boundaries which impose an obstacle to dislocation slip, leading to the

(25)

hard-2 1 Introduction −1.5 −1 −0.5 0 0.5 1 1.5 −1.5 −1 −0.5 0 0.5 1 1.5 σxx/σy σy y /σ y

(a) Yield locus

0 15 30 45 60 75 90 500 600 700 800 θ◦ σy [MP a]

(b) Flow stresses (tension)

0 15 30 45 60 75 90 500 600 700 800 θ◦ σy [MP a]

(c) Flow stresses (compression)

0 15 30 45 60 75 90 1 1.5 2 2.5 θ◦ R (d) R values (tension) 0 15 30 45 60 75 90 1 1.5 2 2.5 θ◦ R

(e) R values (compression)

Experiments Reference model Proposed function

Figure 1.1: Measured material properties and model predictions.

ening of the material. Hence, a Hall-Petch eect takes place which is modelled by a corresponding mathematical law. Because twin formation was found to be dependent on the direction of the applied stress, twin volume fraction is considered as an internal variable, which is updated based on the deformation history.

By that means, it is possible to capture the non-monotonically decreasing hardening rate observed in some experiments (Fig. 1.2).

It will be shown that by the use of both models ow surface and directional hardening the accuracy of forming simulations can be signicantly improved compared to current state of the art models.

(26)

1.2 Outline of the thesis 3 0 0.1 0.2 0.3 600 800 1000 1200 εplog |σ | [MP a] 0◦Model 90◦Model 0◦Exp. 90◦Exp.

Figure 1.2: Measured and simulated ow curves for uni-axial compression in rolling and transverse direction.

1.2 Outline of the thesis

After a general introduction to the basics of material science and the nite element method, the plastic behaviour of the investigated material is charac-terized in chapter 3. It will be shown that classical models cannot capture the occurring eects with sucient accuracy and, therefore, enhanced mod-els are needed. A ow criterion able to map plastic anisotropy and tension-compression asymmetry will be developed in the fourth chapter. Modelling of anisotropic hardening due to directional twinning is discussed in chapter 5. A mechanism-based empirical model is introduced there. The integration of the developed models into commercial nite-element solvers will be out-lined in chapter 6. In the seventh chapter the usability of the models for real applications will be discussed.

(27)

(28)

2 Fundamentals

2.1 State of the art

2.1.1 Flow criteria

If a metallic material is subjected to mechanical strain, it will rst deform elastically and after reaching a certain limit plastically (ow). The predic-tion of the onset of ow can be dicult, especially for complex stress states. Numerous models have been created for this purpose.

Huber (1904), Mises (1913), and Hencky (1924) independently developed a criterion based on strain energy. Later it was recognised that this function cannot describe the behaviour of sheet metals with sucient accuracy, since they have a texture and, therefore, direction dependent properties. By adding empirical coecients to von Mises' model, Hill (1948) presented a rst simple concept to describe anisotropy. Still, this function is in broad use because of its simple calibration and low computational costs.

Hershey (1954), Hosford (1972), and Hill (1979) formulated a gener-alized notion of mechanical energy and proposed non-quadratic models, in-troducing further exibility. Nevertheless, as these criteria were formulated in terms of principal stresses, planar anisotropy could not be modelled. Lo-gan and Hosford (1980) conrmed the applicability of non-quadratic models based on crystal plasticity calculations. Hosford (1985) re-introduced planar anisotropy.

First theoretical considerations about the determination of yield surfaces based on crystal plasticity calculations were already presented by Bishop and Hill (1951a) and Bishop and Hill (1951b). Because of the lack of computa-tional power, a numerical application could only be carried out decades later; Barlat and Richmond (1987) showed good agreement between the Bishop and Hill model and the macroscopic model of Hill (1979). Barlat and Lian (1989) adopted this formulation and introduced linear transformations for an improved modelling of anisotropy. An extension from plane to general stress states was proposed by Barlat et al. (1991b). By adding additional transformations, Barlat et al. (2003) and Banabic et al. (2005) designed

(29)

6 2 Fundamentals

a more exible criterion for plane stress states; Barlat et al. (2007) could prove the equivalence of these two models. The same approach was applied to general stress states by Bron and Besson (2004) and Barlat et al. (2005). A geometric representation of the ow surface was proposed by Vegter and van den Boogaard (2006), who interpolated between experimental points in stress space by means of Bézier splines. Moreover, the spline parameters were dened as functions (Fourier series) of the angle between the rst principle di-rection and the rolling didi-rection. Tong (2006) presented a similar idea, where the parameters of the equivalent models after Hill (1979) and Logan and Hosford (1980) were modelled as Fourier series of the same angle.

Besides the modelling of anisotropy, many authors also tackled the mapping of tension-compression asymmetry. This phenomenon can be observed pre-dominantly in metals with hexagonal crystal lattice such as titanium, magne-sium and zirconium, but also in some aluminium grades and martensitic steels. The eect was already reported by Hirth and Cohen (1970), Pampillo et al. (1972), and Winstone et al. (1973). Its cause, however, remained a matter of debate. Hirth and Cohen (1970) discussed dierent possible reasons, such as micro cracks, residual stresses and internal Bauschinger eects as well as pressure sensitivity. Latter was also proposed byLowden and Hutchinson (1975). Twinning, which is activated depending on the sign of the acting stress, was listed as a possible cause by Hosford and Allen (1970). Further, they explained asymmetry with pressure-sensitive dislocation slip, thus a violation of Schmid's law. The reasoning of Hosford and Allen (1970) (twinning and non-Schmid behaviour) appears to be the most widely accepted in literature. Spitzig and Richmond (1984), however, investigated the dependency of ow stress of steel and aluminium on hydrostatic pressure and could show a linear correlation. Nevertheless, the measured volume dilation was more than one order of magnitude lower than predicted by an associated ow rule.

To model tension-compression asymmetry macroscopically without depen-dence on hydrostatic pressure, Yoon et al. (2000) used a shift of the yield surface in stress space, as it is common for capturing the Bauschinger eect. Cazacu and Barlat (2004) modied Drucker's model involving the third invariant of the stress deviator. By the use of generalized invariants, they could map orthotropy. Cazacu et al. (2006) modelled asymmetry (which was assumed to result solely from twinning) by a sign-dependent shift of the eigen-values of the stress deviator. Anisotropy was considered via a linear transfor-mation. This concept was combined with the idea of multiple transformations by Plunkett et al. (2008). Nixon et al. (2010) also introduced anisotropy

(30)

2.1 State of the art 7

to the isotropic variant of the model after Cazacu and Barlat (2004) by applying a linear transformation.

Stoughton (2002) proposed the use of a non-associated ow rule to de-couple ow criterion and ow direction. Especially, this allows for a pressure dependent model without plastic volume change. Such a function was pre-sented by Stoughton and Yoon (2004). Bai and Wierzbicki (2008) used hydrostatic pressure and the third invariant to model asymmetry. Khan and Yu (2012) superimposed an asymmetric, but isotropic function (involving the third invariant) and a symmetric, but anisotropic function (Hill 48). Yoon et al. (2014) combined the concepts of pressure sensitivity, consideration of the third invariant and multiple linear transformations. Manopulo et al. (2017) distorted symmetric ow surfaces by means of the framework presented by Barlat et al. (2011) which originally was intended to model anisotropic hardening (see below).

2.1.2 Anisotropic hardening

During plastic deformation of metals, dislocations are moved across the crystal lattice. Grain boundaries represent an obstacle for this movement, leading to pile-up of the dislocations. Due to this eect an the continuing creation of new dislocations, the resistance of the material against further deformation rises it hardens. If the material tends to build crystal twins, another mode of deformation is available. However, this mechanism also has an inuence on dislocation slip.

Microscopically, the hardening behaviour is thus strongly dependent on the orientation of the slip and twin systems with respect to the loading direction. If the material is textured (hence anisotropic), this eect is also observable on the macroscopic scale. Moreover, the grains may rotate during deforma-tion to a preferred orientadeforma-tion (texture evoludeforma-tion), which promotes direcdeforma-tional hardening.

Numerous modelling approaches for the mathematical description of aniso-tropic hardening were presented. In some studies, crystal plasticity models were used to predict the plastic properties of a material subjected to non-proportional load paths (e.g. Graff (2007) for magnesium, Peeters et al. (2000), Dancette et al. (2012), Raemy et al. (2015), and Jeong et al. (2016) for steel).

Microscopic simulations are able to predict anisotropic hardening quite ac-curately and give insight to the mechanisms that take place in the material. Nevertheless, for macroscopic simulations, crystal plasticity is not a suitable

(31)

8 2 Fundamentals

method because a real work piece consists of myriads of grains. State-of-the-art computers and even supercomputers are not able to handle and process this vast amount of information. In FEM, macroscopic models are therefore needed.

Many macroscopic model make use of possibly tensor-valued internal variables which are subjected to evolution laws. Chaboche (1986) moved the ow surface in stress space by means of back stresses. The direction and magnitude of those tensors are modied during the plastic deformation. This phenomenological approach was enhanced and improved by many authors.

Feigenbaum and Dafalias (2007) proposed a tensor of fourth order to transform the stresses. Thermodynamic equations were used to modify the parameters of the transformation. Shi and Mosler (2013) and Shi et al. (2015) modied the model and applied it to magnesium.

Barlat et al. (2011) introduced a micro-structure deviator, which repre-sents grain rotation due to plastic ow. This tensor and other internal variables evolve based on the deformation history. A decoupled modelling of dierent phenomena such as Bauschinger eect and permanent softening is made pos-sible by this theoretical framework. An extension thereof was presented by Barlat et al. (2013), who could also capture stress overshoot at an orthogo-nal change of the loading direction.

Even without a change of the loading direction, directional hardening can be observed, since the material hardens in dierent manners when loaded in dierent directions. In order to capture this eect, many authors (e.g. Wang et al. (2009), Nixon et al. (2010), and Peters et al. (2013)) distort the ow surface by dening its parameters in dependency of plastic strain.

2.2 Titanium and titanium alloys

Titanium is the fourth most common metal in earth crust, but only rarely ap-pears in high concentrations (Peters et al. 2002). It is utilized in aerospace, architecture, process engineering, medical industry (because of its biocompat-ibility, bioadhesion and Young's modulus which is similar to that of bones), energy production, o-shore applications as well as sport articles (Peters and Leyens 2002). The major part of the global usage, though, is caused by titanium dioxide, which is a versatile inert white pigment.

Titanium has an extraordinarily high specic strength (Young's modulus/ density ratio) and shows very good resistance against corrosion. In human

(32)

2.2 Titanium and titanium alloys 9

body uid, which has a pH value of around 7.4 and a salt concentration of 0.9 %, it is signicantly more durable than stainless steel (Breme et al. 2002). Pure titanium at low temperatures has a hexagonal crystal lattice (hcp, α titanium), above the transition temperature of 882◦_{C, the lattice changes to}

cubic body-centred (bcc, β titanium). Frequently used alloys in the medical industry are TiAl6V4 und TiAl6Nb7. Vanadium as well as niobium stabilize the β phase, those alloys show an α+β microstructure. Vanadium is suspected to have cytotoxic eects, thus niobium alloys are now preferred (Kobayashi et al. 1998; Niinomi 1998). The advantages of those alloys are the increased strength and corrosion resistance compared with pure titanium.

Commercially pure (CP) titanium is classied into grades. Primarily, they dier in their oxygen content, higher grades having more oxygen. This element occupies interstitial sites, which leads to a distorted lattice and thus increased strength, but diminished ductility. Carbon and iron are impurities due to the manufacturing process and allowed up to a certain content.

2.2.1 Micromechanical behaviour during forming

Like magnesium and zirconium, for instance, titanium has a hexagonal close-packed crystal lattice (hcp, Fig. 2.1). Plastic deformation of metals is enabled by dislocation slip, thus the movement of crystallographic defects. They allow for a relative movement between individual layers of atoms. This process takes place more easily, if the slip plane and direction are densely packed (Gottstein 2007).

a

a c

(33)

10 2 Fundamentals

Further, the number of slip systems has a signicant impact on the deforma-tion resistance. hcp lattices possess less slip systems than face-centred cubic (fcc) and body-centred cubic (bcc) lattices and are therefore comparatively harder to form, although the packing eciency of the hcp lattice is the highest possible. The reason for this is that, if the aspect ratio is ideal, this lattice only has one slip plane the basal plane and three corresponding directions (Fig. 2.2).

(a) Basal slip (b) Primatic slip (c) Pyramidal slip Figure 2.2: Slip systems of the titanium lattice.

For titanium, the c/a ratio (Fig. 2.1) is 1.587 and thus lower than the ideal value of 1.633. Therefore, the packing eciency of the prismatic and pyramidal planes are increased compared to the basal plane. Thereby, slip is rendered possible on the other planes as well. Slip along the vertical axis, however, can only be accommodated on the pyramidal planes. On the other planes, it is only possible in a direction parallel to the base (Fig. 2.2). A change of height is thus hardly feasible. Nevertheless, a part of the deformation can be realised through crystal twinning (Battaini 2008).

Twinning plane Twinning plane

Matrix Twin Matrix

(34)

2.3 The Finite Element Method 11

If twinning occurs, a region of the crystal lattice tilts over into a mirror-inverted position (Fig. 2.3), while the lattice structure is retained. This process takes place spontaneously and propagates with sonic velocity.

2.3 The Finite Element Method

Numerical simulations are nowadays an indispensable tool for the design and optimization of tools and manufacturing processes. E.g. they allow for the prediction of the spring-back of the work piece and for the corresponding cor-rection of the forming tools even before the rst piece is really produced. Enor-mous cost and time saving potential is the consequence. Although numerous software packages for the solution of structural-mechanical problems exist and compete, virtually all of them are based on the nite element method (FEM). In the scope of this thesis, the programs LS-DYNA and simufact.forming, which uses the solver of MSC Marc R_{, were used.}

In the following, the governing equation of the FEM are derived and the in-terfaces to the developed models are illustrated. Viscous damping is neglected for sake of clarity.

2.3.1 Governing equation

The fundamental problem of continuum mechanics consists of a set of coupled partial dierential equations. These are the dynamic equilibrium (momentum balance, Nasdala (2012))

div (σ) + ρb_{− ρ¨}u = 0, (2.1) the kinematic relation between displacements and strains

D := 1 2 h

grad ( ˙u) + grad ( ˙u)Ti (2.2) and the in general non-linear relation between stresses and strains (consti-tutive law)

˙

σ = ˙σ (D) . (2.3) Hereby, u denotes displacement of a material point, ρ the density, b the sum of volumetric forces (e.g. gravity or electromagnetic forces), σ Cauchy stress and D the rate of deformation tensor i.e. the symmetric part of the velocity

(35)

12 2 Fundamentals

gradient (Betten 2001). For suciently small displacements, the latter cor-responds approximately to the linear strain tensor1

D≈ ˙ε. (2.4) The equations are subjected to boundary conditions such as external forces (Neumann boundary conditions) and prescribed displacements (Dirichlet boundary conditions). If it is a dynamic problem, initial conditions are needed as well.

The basic idea of the FEM is the division of the body into small, more or less uniformly shaped domains Ω (nite elements) and to solve the continuum-mechanical problem for those individual elements. It is crucial that the equa-tions will not be fullled at every point of the body, but in an averaged sense over all elements. If the elements are suciently small and therefore numer-ous, an accurate approximation of the solution of the overall problem can be assembled.

To this end, the dynamic equilibrium is transformed from the so-called strong form (2.1) to the weak form by integration over the volume

Z Ω div (σ)dV + Z Ω ρ (b_{− ¨}u) dV = 0. (2.5)

Now, the scalar product with an admissible innitesimal virtual displacement eld δξ is formed. The principle of virtual work (c.f. Dankert and Dankert (2013)) dictates that the result has to vanish identically for all δξ.

Z Ω div (σ)· δξdV + Z Ω ρ (b− ¨u)· δξdV ≡ 0. (2.6)

This can be rewritten by means of the identity div (Ab) = div (A)·b+A·grad b (A denotes a tensor of second order and b a vector)

Z Ω div (σδξ)dV − Z Ω σ· grad (δξ)dV + Z Ω ρ (b− ¨u)· δξdV = 0 (2.7)

1_{In non-linear continuum mechanics, numerous dierent tensors for the measurement of}

stresses and strains are common. LS-DYNA as well as MSC Marc R

use Jaumann rate of Cauchy stress (co-rotated Cauchy stress) and the co-rotated rate of deformation.

(36)

and by the application ofR_Ωdiv (σ)dV =R_∂ΩσndA(Gauss's divergence the-orem) Z Ω σ_{· grad (δξ)dV −} Z Ω ρ (b_{− ¨}u)_{· δξdV −} Z ∂Ω σn_{· δξdA = 0,} (2.8)

where n denotes the unit outward normal on the domain boundary ∂Ω and dAa surface segment (according to Cauchy's fundamental theorem, t = σn is the stress vector at the outer boundary of the volume Ω).

The nite elements are spanned by nodes and the nodal displacements ˆu are treated as variables. By means of the interpolation matrix H (x), the displacement of an arbitrary point within the element can be approximated as u (x) = H (x) ˆu. (2.9) Particularly, for virtual displacements, δξ (x) = H (x) δˆξ holds if δˆξ denotes virtual nodal displacements. Further, a matrix B (x) can be derived which yields the linearised strains from the nodal displacements (Bathe 2002).

To simplify the calculations, and to save memory, Voigt notation is usually used for tensors. By exploiting the symmetry of stress and strain tensors, they can be written as vectors2

σ =σ11 σ22 σ33 σ23 σ31 σ12 T (2.10a) ε =ε11 ε22 ε33 γ23 γ31 γ12 T . (2.10b) It has to be considered, that shear strains εij are replaced by engineering

shear strains γij= 2εij. In that way, the incremental mechanical work can be

calculated by the scalar product

dW = σ_{· dε.} (2.11) Tensors of fourth order, such as the elastic stiness Ce_{, become tensors of}

second order. Thus, Hook's law can be written as σ = Ce

εe, if εe denotes elastic strains.

The (total) strains can hence be written in Voigt notation as

ε (x) = B (x) ˆu. (2.12)

2_{The order of the shear components depends on the software. The used conventions can}

(37)

14 2 Fundamentals

Since Cauchy stress σ is symmetric, the calculation of the inner product in (2.8) only requires the symmetric part of grad (δξ). This corresponds to the virtual strains which result from the virtual nodal displacements.

sym (grad (δξ)) = Bδˆξ. (2.13) As the virtual nodal displacements are independent of the position x, they can be taken out of the integral and

δˆξ·  Z Ω ρHTHdV ˆu +¨ Z Ω BTσdV −  Z Ω ρHTbdV + Z ∂Ω HTσndA     = 0 (2.14) follows. This must hold for all δˆξ and thus

Z Ω ρHTHdV ˆu +¨ Z Ω BTσdV = Z Ω ρHTbdV + Z ∂Ω HTσndA (2.15)

is postulated. The rst integral can be interpreted as element mass matrix M, the second as internal forces (without inertial forces) Fint_{, and the entire}

right-hand-side as the sum of all external forces Fext_{, thus the Neumann boundary}

conditions. The integrals of the left-hand-side are approximated numerically by means of Gaussian integration (weighted sum of the function values at several discrete integration points). Finally, the discrete fundamental equation of FEM reads

M ˆ¨u + Fint= Fext. (2.16) This is valid for the nodal forces and displacements of one element. However, since the elements of a body share several nodes, a global system of equations can be assembled for the problem as a whole.

Numerous reasons can lead to non-linearities in the calculation. For in-stance, boundary conditions can change over time, large displacements and strains can occur, instabilities may arise or the material can show non-linear behaviour. Therefore, (2.16) in general has to be solved incrementally.

2.3.2 Elastic-plastic constitutive law

As mentioned above, metal deforms rstly elastically and after exceeding the yield limit plastically. A total strain increment, therefore, in general consists

(38)

of an elastic and a plastic part

∆ε = ∆εe+ ∆εp. (2.17) The stress increment follows by the linear Hook's law from the elastic strains ∆σ =Ce∆εe. (2.18) To determine whether the material undergoes elastic or plastic deformation for a given stress state, a ow function (chapter 4) is necessary. This function returns a scalar equivalent stress ¯σ from a stress tensor σ. If ¯σ is lower than the current ow stress σy, thus if

Φ := ¯σ (σ)_{− σ}y≤ 0 (2.19)

holds, purely elastic deformation is assumed

∆ε = ∆εe. (2.20) However, if (2.19) is violated, elastic-plastic deformation takes place. In this case, the elastic and plastic part of the strain increment have to be separated (this is treated in section 6.1). To this end, the ow direction has to be calculated by means of the ow rule. In plasticity theory, the existence of a ow potential g in stress space is assumed. The gradient of this eld is coaxial to the plastic strain increment

∆εp= ∆γ∂g

∂σ, (2.21) where ∆γ is the plastic multiplier. The commonly accepted associated ow rule states that the potential coincides with the ow function

g_{≡ ¯σ,} (2.22) in which case the plastic multiplier corresponds to the increment of equivalent plastic strain

∆γ_{≡ ∆¯ε}p (2.23) which follows from the equivalence of plastic work

(39)

16 2 Fundamentals

given that the ow function is homogeneous of degree 1. Bishop and Hill (1951b) support this assumption, but cannot give a general proof. Several authors (e.g. Stoughton 2002), therefore, dispense with this constraint in order to design more exible models.

The hardening of the material is reected by the increase of ow stress during plastic deformation. The modelling of this eect is treated in chapter 5.

For implicit schemes (section 2.3.4), also the derivative of the stress in-crement with respect to the total strain inin-crement Cep

ijkl = ∂∆σij/∂∆εkl is

needed. The determination of this tensor is discussed in section 6.2.

If user-dened material laws shall be considered, they must be provided to the FE solver as subroutines. The major quantities that are passed from the main program to the subroutine are stress at beginning of the increment, and the total strain increment. The expected returned quantities are stress at the end of the increment, and for implicit schemes the continuum tangent. Further, internal variables such as total equivalent plastic strain have to be updated.

2.3.3 Explicit time integration

Explicit programs solve the dierential equation (2.16) through an Euler-forward scheme. Approximating the nodal accelerations by the central dif-ference method, at increment n

ˆ¨un=

ˆ

un+1− 2ˆun+ ˆun−1

(∆t)2 (2.25) follows with a time step ∆t. Inserted to (2.16)

M (ûn+1− 2ûn+ ûn−1) + (∆t)2Fintn = (∆t)2Fextn . (2.26)

This can be solved for the wanted nodal displacements at increment n + 1 ˆ un+1= (∆t) 2 M−1hFextn − F int n i + 2ˆun− ˆun−1. (2.27)

Explicit FEM for quasi-static processes

By means of an explicit scheme, hyperbolic problems such as wave propagation can be solved eciently. However, it is also applicable for quasi-static problems such as forming operations. It must be pointed out that the time step ∆t has to

(40)

be chosen such that a sound wave propagates at most across one (the smallest) element, otherwise the solution becomes unstable. To allow for larger time steps, the density can be articially increased to reduce the speed of sound. However, in this case also inertial forces become more signicant and may aect the accuracy. The mass scaling factor is therefore limited.

If no re-meshing occurs, the mass matrix M is constant and thus only has to be inverted once. Further, the masses can be lumped into one integration point and the matrix can thus be approximated by a diagonal matrix, rendering the inversion trivial. In this case the calculations are even further simplied as the matrix-vector multiplication can be replaced by a element-wise vector-vector multiplication.

The calculation of the term Fint n =

R

Ω

BTσndV requires knowledge of the

stresses at the integration points. Therefore, a generally non-linear constitutive law is needed to determine them based on the strains εn(x) = B (x) ˆun. As

relatively small time steps are used in explicit analyses, the constitutive law has to be evaluated many times for each integration point. The total calculation time is thus strongly dependent on the numerical costs of the law.

For the investigations carried out in the later chapters, the material models were integrated in the explicit solver of LS-DYNA.

2.3.4 Implicit time integration

Implicit solvers are commonly used to solve quasi-static processes. In this case, inertial forces can be neglected and (2.16) is reduced to

Fint= Fext. (2.28) Because of non-linearities arising during (forming) simulations, this equation is discretized in time, and for each increment a linear approximation is solved. At the beginning of increment n, an estimation of the nodal displacements is generated based on the boundary conditions. The internal forces Fint

n =

R

Ω

BTσndV are determined and the dierence to the external forces Rn =

Fintn − F ext

n is calculated. If this residuum exceeds a chosen tolerance3, the

nodal displacements are corrected by means of a Newton-Raphson method.

3_{A norm has to be dened as a measure of the dierence R}

n. In MSC MarcR

, the maximum norm is used, thus the vector component with the highest absolute value.

(41)

18 2 Fundamentals

To this end, the constitutive law (2.3) is linearised; in Voigt notation the tensor operation

∆σn= Cepn ∆εn= Cepn B∆ˆun (2.29)

results, where ∆ˆun denotes the correction of the nodal displacements. The

calculation of the continuum tangent Cep

= ∂∆σ/∂∆εis treated in chapter 6. The correction of the internal forces follows approximately as

∆Fintn =

Z

Ω

BTCnBdV ∆ˆun = Kn∆ˆun (2.30)

with a stiness matrix Kn. Now, if ∆F int

n =−Rn is postulated,

∆ˆun =−K−1n Rn (2.31)

follows. Thus, a system of equations has to be solved in this step. After the update of the nodal displacements, convergence is checked again. This scheme is repeated until the chosen convergence criteria are fullled. Subsequently, the next increment is calculated.

Because of the convergence controls during an implicit analysis, relatively large time steps can be chosen. However, the evaluation of one increment is signicantly slower than for an explicit analysis, because of the system of equation to be solved. Therefore, the complexity of the constitutive law has signicantly less impact on the calculation time.

2.3.5 Thermal problems

FEM is not limited to the solution of structural mechanical problems, it is a general numerical tool. Other important applications are heat conduction and transmission problems. Especially for forming simulations, it is often necessary to solve thermo-mechanically coupled problems, as e.g. plastic deformation generates heat which is transmitted to the environment.

A similar approach as above is used to derive the FE equations for the thermal problem. The fundamental equation of heat conduction (Fourier's law) reads for spacial problems as

cρ∂T

(42)

Here, the specic heat capacity c, density ρ, heat conduction coecient λ and Temperature T are, in general, depending on the location x and time t; the coecients can also depend on T . q denotes inner heat sources or sinks. Boundary conditions such as prescribed temperatures (Dirichlet boundary con-ditions) or heat uxes (Neumann boundary concon-ditions), and initial conditions close the problem. The heat uxes can arise from convection of the ambient medium, body contact, friction, or radiation.

Integrating and multiplying (2.32) with a virtual temperature eld δζ yields Z Ω cρ∂T ∂tδζdV − Z Ω div (λ grad (T ))δζdV ₋ Z Ω qδζdV = 0, (2.33)

which has to be fullled for all admissible δζ. By means of Gauss's divergence theorem, one nds Z Ω cρ∂T ∂tδζdV − Z ∂Ω λ grad (T )· nδζdA + Z Ω λ grad (T )_{· grad (δζ)dV −} Z Ω qδζdV = 0. (2.34)

The term λ grad (T ) · n corresponds to the heat ux qs through the element

boundaries, which is caused by the aforementioned reasons. It can be summa-rized as

qs= (Tu− Ts) h, (2.35)

where Tu denotes ambient temperature, Ts surface temperature and h a

gen-eralized heat transmission coecient.

The nodal temperatures are stored in the vector ˆT . Within the elements, the temperature eld and its gradient can be approximated by means of the interpolation matrices Hλ(x)and Bλ(x)

T (x) = Hλ(x) ˆT

grad (T (x)) = Bλ(x) ˆT .

(2.36)

For the surface temperature, an additional matrix Hs(x)is dened and

(43)

20 2 Fundamentals holds. (2.34) is reformulated δˆζ "Z Ω cρHTλHλdV ˆ˙T + Z ∂Ω hHTsHsdA ˆT− Z ∂Ω hHTsTudA + Z Ω λBTλBλdV ˆT − Z Ω HTλqdV # = 0, (2.38) and thus Z Ω cρHTλHλdV ˆ˙T +  Z Ω λBTλBλdV + Z ∂Ω hHTsHsdA   ˆT = Z Ω HTλqdV + Z ∂Ω hHTsTudA. (2.39)

The rst integral on the left-hand-side can be interpreted as heat capacity matrix Cλ, the sum in parentheses as thermal stiness Kλ. The rst

inte-gral on the right-hand-side is the sum of the heat generated in the element. Together with the second integral, it is summarized in the vector Qλ. The FE

equation for the thermal problem thus reads

Cλˆ˙T + KλT = Qˆ λ. (2.40)

To solve it, again implicit or explicit integration schemes can be used (not shown here).

Possible internal heat sources are electric current or, especially in forming technology, plastic work. The latter is a result of the structural mechanical cal-culation. On the one hand, the constitutive law used therein might depend on temperature. On the other hand, temperature changes induce thermal strains and inhomogeneous temperature elds can cause thermal stresses. Thus, both problems are coupled.

2.3.6 Thermo-mechanical coupling

Usually, the coupled problem is solved by a weakly coupled a staggered scheme. This means, for a specic time step, rstly a purely mechanical analysis is carried out. The plastic work is calculated and passed to the

(44)

ther-2.3 The Finite Element Method 21

mal analysis (see section 6.3), where the temperature eld as well as thermal stresses and strains are updated. This new data enters the mechanical analysis of the next time step, and so on. During the mechanical pass, constant tem-perature is assumed and vice versa, during the thermal pass, all mechanical variables are assumed to constant.

If the two partial problems are actually combined to one system of equa-tions, where e.g. the temperature dependence of the constitutive law has to be considered, a strongly coupled analysis is performed. However, in general, weakly coupled analyses deliver suciently accurate results since the maxi-mal change of temperature during one increment is monitored. If the change exceeds a certain value, the time step is decreased and the increment will be re-evaluated.

(45)

(46)

3 Material characterization

All investigations mentioned in this work concern unless otherwise specied commercially pure titanium sheet (Ti-CP Grade 4). It contains 0.4 % oxygen which acts as an interstitial alloying element. Because of the resulting distor-tion of the lattice, the strength of the alloy is signicantly increased compared to high purity titanium. However, ductility degrades as well. Forming at ele-vated temperatures is, therefore, required to achieve larger strains (Peters et al. 2002). Due to the manufacturing process, further impurities are contained in the alloy. Their maximum permissible mass fractions are summarized in Tab. 3.1.

Table 3.1: Allowed impurities (alloying elements) in commercially pure tita-nium (Ti-CP Grade 4), weight per cent.

O Fe C N H others 0.400 0.500 0.080 0.050 0.015 0.300

The investigated material shows very high corrosion resistance and is, there-fore, preferably used in the chemical and medical industries. The latter also makes use of the excellent biocompatibility and bioadhesion (Breme et al. 2002). It is an α alloy i.e. virtually all grains have an hcp crystal lattice.

In the following, the plastic behaviour of the material is analysed by uni-axial tension and compression tests. This data is later used to calibrate the developed material models which in turn are expected to increase the accuracy of forming simulations.

3.1 Flow behaviour

To quantify the onset of plastic ow and to characterize the hardening be-haviour, uni-axial tension and compression tests were conducted. The param-eters temperature, strain rate, and loading direction (angle to rolling direction) were varied for the dierent experiments.

(47)

24 3 Material characterization

The tests were carried out on a deformation dilatometer Bähr DIL 805. The specimen geometries are listed in appendix A. Additionally, to measure Lank-ford coecients, standard quasi-static tensile tests according to ISO 6892-1 were conducted at room temperature.

The loading directions were 0◦_{, 45}◦_{and 90}◦_{to rolling direction (RD). Normal}

direction (through thickness, ND) was also considered for the compression tests in order to measure biaxial ow stresses1. For the in-plane experiments, sheets of thickness 3 mm was used, compression in normal direction was carried out on a sheet of thickness 5 mm due to physical limitations of the testing device. Strain rate of 0.05 s−1_{, 0.5 s}−1 _{and 5 s}−1 _{were applied at temperatures}

of 20◦_{C, 200}◦_{C, 400}◦_{C, 500}◦_C_{and 600}◦_{C, but not all possible combinations}

were realized. The inuence of the dierent parameters are visualized by ow curves.

3.1.1 Temperature dependence

Fig. 3.1 shows isothermal ow curves at a constant strain rate of 0.5 s−1_{. As}

expected, ow stresses decrease with increasing temperature in all loading directions, as dislocation slip is eased. However, the compression curves in transverse direction (TD) show a particular behaviour. Up to a temperature of around 500◦_{C, the hardening rate (slope of the curve) does not monotonically}

decrease as it is the case in all other directions. These curves show a sigmoidal shape, which is a typical feature of hcp metals. A physical explanation for this is given in the next section.

Further, it is observed that the yield limit is higher in TD than in RD for the tensile tests at all temperatures. This dierence of ow stress vanishes with increasing strain. For the compression tests, the dierence of the yield limits is even more signicant. In this case, however, it is higher for RD than for TD at all temperatures. Because of the eect mentioned above, the situation can be reversed; the ow stress in TD can exceed the one in RD at higher strains.

3.1.2 Direction dependence

The strongly diering hardening behaviour of the dierent loading directions is pointed out in Fig. 3.2, which shows ow curves at room temperature. Diagonal direction (DD) and normal direction are shown here as well. Besides

1_{Uni-axial compression in normal direction and biaxial tension produce the same deviatoric}

stress state. Assuming pressure invariance, one of the experiments can, therefore, be used to replace the other.

(48)

3.1 Flow behaviour 25 RD,T = 20◦C TD,T = 20◦C RD,T = 200◦C TD,T = 200◦C RD,T = 400◦C TD,T = 400◦C RD,T = 500◦C TD,T = 500◦C RD,T = 600◦C TD,T = 600◦C 0 0.02 0.04 0.06 0.08 0.1 0 200 400 600 800 1000 1200 εplog |σ | [MP a]

(a) Uni-axial tension

0 0.1 0.2 0.3 0 200 400 600 800 1000 1200 εplog |σ | [MP a] (b) Uni-axial compression Figure 3.1: Experimental iso-thermal ow curves at strain rate 0.5 s−1

(Ti-CP Grade 4). Tensile ow curves measured on deformation dilatometer.

(49)

26 3 Material characterization Compression, RD Tension, RD Compression, TD Tension, TD Compression, ND Tension, DD 0 0.05 0.1 0.15 0.2 0.25 0.3 600 800 1000 1200 εplog |σ | [MP a]

Figure 3.2: Quasi-static ow curves in dierent directions at room temperature (Ti-CP Grade 4). Tensile ow curves measured by standard tensile tests.

the already observed anisotropy (dierent yield limits in dierent directions), also an asymmetry (dierent yield limits when the sign of stress changes) can be noted. This asymmetry is signicant in RD, where the strength under compression is higher than under tension. Asymmetry is weak in DD and for TD, it even reverts its behaviour. No tensile experiments were possible in ND because of practical reasons. Nevertheless, no sigmoidal shape is recognised in the compressive ow curve for this direction.

3.1.3 Strain rate dependence

The inuence of the strain rate is depicted in Fig. 3.3 exemplarily for com-pression in TD. It is observed that for room temperature, the increase of ow stress at higher rates is relatively small compared to the quasi-static curve the fastest experiment was conducted with a hundredfold higher strain rate than during the slowest. As expected, the strain rate sensitivity is more pro-nounced at higher temperatures. However, after a certain strain, the curves of the fastest experiment show less hardening than the curves of the slower ones. This eect is caused by the produced heat due to plastic work, which cannot be evacuated in a suciently fast manner. During the calibration of the hardening models (chapter 5), this temperature change was considered.

(50)

3.2 Microstructure 27 T = 20◦ T = 200◦ T = 400◦ ˙εlog= 0.05 s−1 T = 20◦ T = 200◦ T = 400◦ ˙εlog= 0.50 s−1 T = 20◦ T = 200◦ T = 400◦ ˙εlog= 5.00 s−1 0 0.05 0.1 0.15 0.2 0.25 0.3 500 1000 εplog |σ | [MP a]

Figure 3.3: Compressive ow curves in TD at dierent temperatures and strain rates (Ti-CP Grade 4).

3.2 Microstructure

As shown in the previous section, signicantly dierent hardening character-istics can be observed for compressive tests in dierent loading directions. Particularly, tests in RD resulted in a ow curve with monotonically decreas-ing hardendecreas-ing rate, as it is known from materials such as steel or aluminium. In TD, however, three phases can be identied: A rst phase with decreasing rate, a second with increasing rate and nally a third phase with decreasing rate again, giving the particular S shape of the ow curve. This eect is pronounced up to a temperature of about 500◦_{C, above this temperature it}

abates.

This three-phased hardening behaviour, which is characteristic of hcp met-als, has been identied in numerous studies. However, usually no direction dependence is reported, provided it was investigated at all (cf. Jain and Ag-new (2007)). Tritschler et al. (2014) predict a direction dependence based on VPSC calculations and nd a three-phased hardening in TD too, but not in the other directions.

The detected dierences suggest that twinning is a directional phenomenon. To substantiate this thesis, micrographs of several compressed specimens were

(51)

28 3 Material characterization

(a) 20◦_{C, Rolling direction} _{(b) 20}◦_{C, Transverse direction}

(c) 200◦

C, Rolling direction (d) 200◦

C, Transverse direction

(e) 400◦

C, Rolling direction (f) 400◦C, Transverse direction Figure 3.4: Optical micrographs (polarized light) after uni-axial compression, |log (ε)| = 15 %.

(52)

3.3 Conclusions 29

taken at a strain of 15 %. Polarized light scatters anisotropically at a hexagonal lattice. Therefore, dierently oriented grains appear at a dierent brightness (Oettel 2011), which makes it possible to clearly identify twins even without the use of a EBSD device. Micrographs of experiments in RD and TD at tem-peratures of 20◦_{C, 200}◦_C_{and 400}◦_C_{and a strain rate of 0.5 s}−1 _{are depicted}

in Fig. 3.4.

All micrographs show crystal twinning. However, for all temperatures, it is much more prominent in the specimens compressed in TD. The dierence is blatant at 200◦_C_{and 400}◦_C_{. At these temperatures, a relatively high increase}

of ow stress can be observed (previous section). Explaining the particular hardening behaviour with directional twinning thus seems reasonable.

3.3 Conclusions

The ow characteristics of commercially pure titanium were investigated under the variation of multiple parameters. Both anisotropy (dierent yield stresses in dierent directions) and asymmetry (dierent yield stresses in tension and compression) were observed.

Strongly dierent hardening characteristics were noted for the individual loading directions. This eect was attributed to directional twinning, micro-graphs supported this thesis.

To capture the material behaviour correctly, a constitutive model capable of modelling anisotropy, asymmetry and anisotropic hardening due to directional twinning is, therefore, needed. It will be developed in the following chapters.

(53)

(54)

4 Flow criteria

The existence of numerous ow functions capturing tension-compression asym-metry was pointed out in the introduction. In the following, two commonly used models will be introduced and a new criterion will be developed there-after. A comparison will emphasise the restrictions of the reference models and clarify the advantageous properties of the proposed model.

4.1 Reference models

In numerous studies, the ow functions developed by Cazacu et al. (2006) and Plunkett et al. (2008) were applied, e.g. in (Jia and Bai 2016; Barros et al. 2016; Steglich et al. 2016; Yoon et al. 2014; He et al. 2014; Tritschler et al. 2014). They are thus used as the reference, to which the proposed model shall be compared.

Cazacu et al. (2006) dened equivalent stress as a homogeneous function

¯ σ =   ˇS1 − k · ˇS1 a + ˇS2 − k · ˇS2 a + ˇS3 − k · ˇS3 a Φ1 + k · Φ1 a +Φ2 + k · Φ2 a +Φ3 + k · Φ3 a   1/a , (4.1)

where ˇSi are eigenvalues of the transformed stress deviator.

S := Cs =         C11 C12 C13 0 0 0 C12 C22 C23 0 0 0 C13 C23 C33 0 0 0 0 0 0 C44 0 0 0 0 0 0 C55 0 0 0 0 0 0 C66                 sxx syy szz syz szx sxy         (4.2)

Anisotropy is introduced through this transformation as it is common with models of the Barlat family. The parameter k describes tension-compression asymmetry. The coecients in the denominator consist of the components of

(55)

32 4 Flow criteria

the transformation tensor. Depending on which ow curve shall be used as reference, another combination can be chosen. Because uni-axial ow stress in RD was selected, the terms were set to

Φ1= 2 3C11− 1 3C12− 1 3C13 Φ2= 2 3C12− 1 3C22− 1 3C23 Φ3= 2 3C13− 1 3C23− 1 3C33. (4.3)

Without the exponent a, the model has 9 free parameters. As out-of-plane shear experiments are dicult to conduct, the coecients C44 and C55 are

often set to the isotropic value 1. Hence, 7 degrees of freedom remain for the calibration.

Plunkett et al. (2008) expanded the model by introducing further lin-ear transformations. Here, the variant using two transformations was used. Equivalent stress is dened as

¯ σ = N D 1/a , (4.4) where N = ˇS10 − k0 · ˇS10 a + ˇS02 − k0 · ˇS20 a + ˇS30 − k0 · ˇS03 a + ˇS100 − k00_{· ˇ}_S00 1 a + ˇS002 − k00_{· ˇ}_S00 2 a + ˇS300 − k00_{· ˇ}_S00 3 a and D =Φ01 + k0 · Φ01 a +Φ02 + k0 · Φ02 a +Φ03 + k0 · Φ03 a +Φ001 + k00 · Φ001 a +Φ002 + k00 · Φ002 a +Φ003 + k00 · Φ003 a . Again, uni-axial ow stress is the reference when those coecients are used. All quantities are dened similarly as above, but with two distinct transformations tensors C0 _{and C}00_{and with}

S0:= C0s, S00:= C00s. (4.5) Thus, the model has 14 free parameters to be adjusted.

(56)

4.2 Proposed ow function 33

4.1.1 Eigenproblems

The necessary eigenvalue calculations were carried out by aid of the hybrid al-gorithm presented by Kopp (2008)1. The algorithm is optimized for symmet-ric 3×3-matsymmet-rices. To calculate the derivatives of the eigenvalues with respect to their corresponding tensor, following method was used:2.

Let A ∈ R3×3 _{be a symmetric tensor of second order, ˇ}_A

i an eigenvalue of

A, and ni∈ R

3_{a corresponding normalized eigenvector. Per denition,}

A_{− ˇ}AiI

ni= 0 (4.6)

holds, where I stands for the unit tensor. The total dierential reads dA_{− d ˇ}AiI

ni+ A− ˇAiI

dni= 0. (4.7)

Pre-multiplying (forming the inner product) with ni yields

ni· dAni− ni· d ˇAiIni+ ni· A − ˇAiI

dni= 0, (4.8)

which can be rewritten to

dA_{· (n}i⊗ ni) + d ˇAini· ni+ A− ˇAiI ni | {z } 0 ·dni= 0. (4.9)

Here, no sum over i is implied. Since ni is normalized according to the

re-quirements,

d ˇAi

dA = ni⊗ ni. (4.10) The needed eigenvectors were also calculated with Kopps algorithm.

4.2 Proposed ow function

The models shown above are based on linear transforms of the stress tensor. It is also possible to think about the ow surface as a hyper surface in stress space which is described geometrically. As an example, the ow surface represented by the von Mises model forms, in the subspace of plane stresses, an ellipsoid (see Fig. 4.1).

1_{Fortran code is available at the URL mentioned in the publication.}

2_{Thanks go to Prof. Hubert J. M. Geijselaers, University of Twente, Enschede NL, for this}

(57)

34 4 Flow criteria σyy σyy τxy σxx σxx

Figure 4.1: Flow surface after von Mises in space of plane stresses (σ2 xx− σxxσyy+ σ 2 yy+ 3τ 2 xy= σ 2

y). Dashed lines correspond to uni-axial stress states.

Hue is reecting the value of the τxy component.

Its only adjustable parameter is a scale factor which correlates linearly with ow stress. For an anisotropic and probably asymmetric material, the points at which plastic ow sets in do not all lie on the ellipsoid any more. Furthermore, to map the R values correctly, the surface normal has to be adapted locally.

Vegter and van den Boogaard (2006) used Bézier splines to interpolate the yield locus between measured points in the space principal plane stresses. Thus, they describe a curve on a plane. By shifting the hinge points of the splines, they managed to adjust the locus normal at the experimental points to capture the correct R values. The coordinates of the hinge points are the free parameters of the model. These coordinates are described as a Fourier series of the angle between RD and the direction of rst principal stress, which allows for the consideration of shear stresses.

A similar albeit not geometrically motivated approach was presented by Tong (2006). He varied the anisotropy coecients of the Logan and Hos-ford (1980) model, by expressing them as a Fourier series of the angle between RD and the direction of rst principal stress as well. Altenbach et al. (2014) showed a function comprising two Fourier series of the stress angle, hydrostatic pressure, and the second invariant of the stress deviator.

An explicit geometric description of the ow surface could be realized through NURBS (non-uniform rational B-splines). But, as it is the case with the model of Vegter and van den Boogaard (2006), there is no di-rect relation between the components of the stress tensor and the surface parametrization. To evaluate the ow function, it would, therefore, be

Novel constitutive models for capturing anisotropy, tension-compression asymmetry, and anisotropic hardening of titanium

Research Collection

Doctoral Thesis

Novel constitutive models for capturing anisotropy,

tension-compression asymmetry, and anisotropic hardening of titanium

Novel constitutive models for capturing anisotropy,

tension-compression asymmetry, and anisotropic

hardening of titanium

Acknowledgements

Abstract

Zusammenfassung

Contents

List of Figures

List of Tables

Nomenclature

1 Introduction

1.1 Motivation

1.2 Outline of the thesis

2 Fundamentals

2.1 State of the art

2.1.1 Flow criteria

2.1.2 Anisotropic hardening

2.2 Titanium and titanium alloys

2.2.1 Micromechanical behaviour during forming

2.3 The Finite Element Method

2.3.1 Governing equation

2.3.2 Elastic-plastic constitutive law

2.3.3 Explicit time integration

2.3.4 Implicit time integration

2.3.5 Thermal problems

2.3.6 Thermo-mechanical coupling

3 Material characterization

3.1 Flow behaviour

3.1.1 Temperature dependence

3.1.2 Direction dependence

3.1.3 Strain rate dependence

3.2 Microstructure

3.3 Conclusions

4 Flow criteria

4.1 Reference models

4.1.1 Eigenproblems

4.2 Proposed ow function

4.2 Proposed ow function