PAGE 1
MODELING AND COMPARISON OF DIFFERENT CONSTITUTIVE EQUATIONS FOR ELASTICPLASTIC MATERIALS by Ali, A.A Adam B.A., Omar AlMukhtar University 2003 A thesis submitted to the University of Colorado Denver in partial fulfillment of the requirements for the degree of Master of Science Mechanical Engineering 2011 f\
PAGE 2
This thesis for the Master of Science degree by Ali, A.A Adam has been approved by Date
PAGE 3
Adam, Ali, A.A (M.S., Mechanical Engineering) Modeling and Comparison of Different Constitutive Equations for ElasticPlastic Materials Thesis directed by Professor John Trapp ABSTRACT Detailed attention will be given to the objectivity and symmetry of the constitutive relations of isotropic materials at finite elasticplastic deformation. Another approach for the derivation of the special equation for the free energy and stress response under finite elasticplastic deformation of metals which first was developed by Naghdi and Trapp is introduced. Then, a general function for a yield criterion under the restrictions of objectivity and symmetry is obtained. A matlab code is developed to integrate the constitutive equations. Finally, an example simulation with loading and reverse loading will be used to check the validity of the formulas and the numerical integration scheme.
PAGE 4
This abstract accurately represents the content of the candidate's thesis. I recommend its publication. Signed
PAGE 5
CONTENTS Figures 1. Introduction 1.1 Background 1.2 Objective 2. FiniteDeformation Continuum Mechanics 2.1 2.2 Configurations Deformation gradient 2.2.1 Multiplicative Decomposition of Deformation Gradient 2.3 Measures of Finite Deformations 2.4 Yield Criteria . 2.5 Bauschinger Effect 3. A constitutive Equation for ElasticPlastic :rviaterial 3.1 Objectivity and Symmetry 3.1.1 Principle of Rigid body invariance 3.1.1.1 Deformation Gradient 3.1.1.2 Scalar Function . . 3.1.1.3 Internal Energy Function 3.1.1.4 Cauchy Stress Tensor 3 .1.1. 5 FiolaKirchhoff Stress Tensor 3.1.2 Principle of Material Symmetry v Vlll 1 1 2 3 3 4 5 6 8 8 11 11 12 12 14 14 15 17 18
PAGE 6
3.1.2.1 Rotational Symmetry 3.1.2.2 Reflectional Symmetry 3.1.3 3.1.4 3.2 Principle of Intermediate Invariance Principle of Elastic Symmetry . . The Free Energy Function and The Stress Response Rigid Body Invariance Elastic Symmetry 3.2.1 3.2.2 3.2.3 3.2.4 3.3 Material Symmetry Intermediate Frame Invariance Modified Sima's Yield Function 4. Work hardening, Numerical scheme and results 4.0.1 Measures of Work Hardening and State Equations 4.0.1.1 Measures of Work Hardening 4.0.1.2 State Equation 4.0.2 Numerical Scheme 4.0.3 Results 4.1 Conclusion Appendix A. Solving Linear Equations Using Some Functions A.l Properties Function A.2 Stress Function A.3 Yield Function A.4 Differentiation of The First Equation A.5 Differentiation of The Second Equation vi 19 20 22 22 23 25 26 27 28 31 34 34 34 36 37 38 40 46 57 59 60 61 63
PAGE 7
B. Yield Function Plotting References . . . . . vii 69 71
PAGE 8
FIGURES Figure 2.1 2.2 2.3 2.4 3.1 3.2 3.3 3.4 4.1 4.2 4.3 Body and configuration. Multiplicative decomposition of the deformation gradient . The Bauschinger effect model. . . . . . . . The Kinematic hardening model. Observer transformation. Deformation of a material body observed by two observers. Body and configuration. Body and configuration. The stressstrain curve. Beta and the right CauchyGreen tensor curve. The residual stress and the right CauchyGreen tensor curves. 4.4 Input Data. 4.5 Plots of some outputs 4.6 StressStrain Curve. 4. 7 The initial and subsequent yield surface Vlll 4 5 8 10 11 13 20 21 41 42 42 43 43 44 45
PAGE 9
1. Introduction 1.1 Background First, a rough and intuitive definition of elasticity and plasticity theory is introduced with a simple experiment. If someone holds tightly a thin strip of a metal from one side and applies a load to the free end, then the free end will deflect. By removing the load, if the free end jumps back to its original position, then elastic deformation has taken place. On the other hand, if the load is large enough to not allow the free end to spring back, then the deformation is called plastic. Continuum mechanics is a broad realm that includes both the theory of elasticity and plasticity [1]. Numerous everyday objects are constructed of material which have mostly undergone a plastic forming process. Material forming can be defined as the process of plastically shaping the raw material into a product form. Generally, metal forming is categorized into two groups: bulk metal forming and sheet metal forming. Examples of bulk metal forming are rolling, extrusion, forging etc, and examples of sheet metal forming are stretch forming, bending, spinning etc [2]. The product form, that specifies which forming process to be used, has applications in buildings, cars, airplanes, furniture, tools and many other objects. Over the past two decades, the theory of plasticity has received increasing attention and is an active area of research. Reasonable amount of effort has been given to mathematically describe the constitutive relations, yielding criteria and 1
PAGE 10
hardening rules. Also, many experiments have been made to examine some spe cific relations. A numerical approach is a significant way to simulate and test the validity of plasticity constitutive relations. For the purpose of simplicity, some assumptions, simplifications and approximations must be thoughtfully applied to achieve satisfactory analysis. 1.2 Objective Essentially, this thesis intends to define the principles of objectivity and symmetry in detail. Then, a derivation of constitutive equation for elastic plastic metals that satisfy those principles is included. Then, a numerical scheme and matlab code are included to integrate the constitutive equations and a yield equation. Applying the objectivity and some symmetry principles, a stress equation and yield function will be used as an example to test the integration scheme. Finally, some plots will be used to evaluate some results. 2
PAGE 11
2. FiniteDeformation Continuum Mechanics 2.1 Configurations The fundamental constitutive equation of any plasticity theory is to develop a relationship between the stress and strain. Some concepts such as configura tions should be clearly understood to visualize those relations and distinguish between the different stress and strain measures. The situation is fairly simple in the case of infinitesimal strain in an elastic media, since using any linear stressstrain relation will lead to approximately to the same result. However, the case is entirely different for finite deformations [3]. In the case of largedeformation, the difference between the undeformed and deformed coordinates can no longer be ignored. Therefore, constitutive relations have to be described in the undeformed or deformed configuration. Either the Lagrangian description (Ra) or the Eulerian description (R) could be used to describe the finite deformation of a continuous body. In the Lagrangian descrip tion, the reference or initial coordinates are used to describe the deformation of each particle. In the Eulerian description, the spatial or current coordinates are used to describe the deformation of each particle. figure (2.1) shows both descriptions [ 1]. A motion of a particle in the reference configuration (Ro) that has a position vector X into another location in the current configuration (R) with a position vector x can be described mathematically as x = x(X, t) (2.1) 3
PAGE 12
Figure 2.1: Body and configuration. It should be pointed out that uppercase X refers to the reference coordinates and the lowercase x refers to the spatial coordinates. 2.2 Deformation gradient Material body elements such as points, lines, areas and volumes in the refer ence configuration will be transformed into the corresponding ones in the current configuration due to the deformation. Determining the relationships between the geometrical elements in (Ra) and (R) is the goal of kinematic analysis. The deformation gradient tensor (F) is a linear operator that appears in the transformation equations between elements of lines, areas and volumes in the reference and current configuration. Previously, the relationship between points is de scribed in equation (2.1). The formula of the deformation gradient tensor (F) 4
PAGE 13
in the cartesian coordinate systems is[3] F= ox ax (2.2) It is necessary to mention that the components of the deformation gradient are evaluated with respect to both bases, the reference and current configuration. 2.2.1 Multiplicative Decomposition of Deformation Gradient 0 x,X,p Figure 2.2: Multiplicative decomposition of the deformation gradient. The deformation gradient (F) can be expressed as a multiplication of the elastic deformation gradient (Fe) and the plastic deformation gradient ( FP). This can be explained by imagining an intermediate configuration described by the intermediate coordinate system (P). The intermediate coordinate system 5
PAGE 14
( P) is usually defined by unloading the continuum body from the spatial con figuration to a zero stress state. Figure (2.2) shows both deformation gradients (Fe, FP) and the intermediate configuration. Mathematically, the wellknown multiplicative decomposition of the deformation gradient can be written as [3] (2.3) 2.3 Measures of Finite Deformations Essentially, the aim of any plasticity theory is to mathematically relate the stress to the strain. Formulating a relationship between the stress and the strain and also measuring them is complicated in the case of finite strain. Some of those measures are briefly introduced in this section. The right and left CauchyGreen tensors are wellknown and widely used to derive or formulate constitutive relations. C=FTF B=FFT Where C is the right CauchyGreen tensor. B is the left CauchyGreen tensor. (2.4) (2.5) Using the multiplicative decomposition of the deformation gradient, the following can be obtained 6
PAGE 15
The AlmansiHe mel strain tensor is defined as 1 E =(CI) 2 (2.6) The wellknown Cauchy stress tensor u is described in the current configuration (R). However, it is worthwhile to define a stress tensor that is described in the reference configuration (R0).Some significant stress measures that are defined in the reference configuration are the nonsymmetric first PiolaKirchhoff stress tensor (P) and the symmetric Second PielaKirchhoff stress tensor (S). Those stresses are mathematically given by [3] (2.7) (2.8) Also in the literature of finite plasticity, the Kirchhoff stress tensor T is defined as T = det(F)u (2.9) 2.4 Yield Criteria 7
PAGE 16
In a simple tension test, a yield point is the point that separates the elastic region from the plastic region on a stressstrain curve. If many stresses are acting at a point in different directions then yielding will be represented by a surface. The determination of which combination of stresses will cause yielding to occur is called a yield criterion. For solids, in particular, many yielding criteria have been suggested and proposed [3]. Some restrictions for an isotropic yielding criterion will be discussed and if possible a yielding criterion will be formulated. 2.5 Bauschinger Effect Figure 2.3: The Bauschinger effect model. 8
PAGE 17
To illustrate this effect, consider loading a specimen of a material in tension into the plastic region. Now, assume the initial load is removed then reload in the reverse direction until it yields again. The result is that, the yield stress in the reversed direction is less than the yield stress in the loading direction [3]. This is called the Bauschinger effect as show in figure (2.3) The Bauschinger effect is caused by the anisotropy of the dislocation fields due to the former loading. The existence of the Bauschinger effect complicates the modeling of plastic deformation. In order to simplify this effect, the kinematic hardening model is used. In this model, stress increment beyond the yield stress in the original loading direction is equal to the reduction of the yield stress in the reloading direction[3]. This model is illustrated in figure (2.4) 9
PAGE 18
Figure 2.4: The Kinematic hardening model. 10
PAGE 19
3. A constitutive Equation for ElasticPlastic Material 3.1 Objectivity and Symmetry Formulating constitutive equations in finite theory of plasticity is greatly simplified if general principles for all constitutive equations are first applied. Some restrictions are applicable for all materials and others are only applicable on specific materials that meets some requirements. Detailed attention is given to those general principles that apply to all materials. A review of the theory of observer transformation will be introduced first. 0 o ,e Figure 3.1: Observer transformation. The theory of Observer transformation is a wellknown subject in continuum mechanics. To describe the change of a timespace reference frame, two frames 11
PAGE 20
are considered, a reference and movmg frame. Each frame is considered as an observer and all constitutive equations have to be applicable everywhere such observers are located such as on earth, moon and space. Also, it has to be invariant to the changes of observers. Considering two frames by (0, e) and ( 0*, e*) that represent two observers 0, 0*. 0 and 0* are their respective origins and their respective base frames are eb e 2 e 3 and e;, e; [4]. A point x in the (0, e) is related to a point x* in the (0*, e*) frame by x* = Qx + c* where c (t) is a vectorvalued function oft. Q (t) is a proper orthogonal tensorvalued function oft. z.e QQT = QTQ =I and detQ = 1. 3.1.1 Principle of Rigid body invariance (3.1) This principle states that any constitutive equation should be invariant under changes of observer frames. Since constitutive equations have scalars, vectors and tensors as variables, understanding how they are related under changes of frames is basically significant. Now, let us consider two observers watching a material body deforms under an application of load as shown in figure (3.2). 3.1.1.1 Deformation Gradient According to the observers 0 and 0*, the deformation of line will have the following formula 12
PAGE 21
F Figure 3.2: Deformation of a material body observed by two observers. dx = FdX dx* = F*dX Also, we know dx* = Qdx Which leads to F*dX = QFdX Resulting in a useful formula describing a relationship between the deformation gradients of a deformed material being viewed by both observers 0 and 0* F* = QF (3.2) 3.1.1.2 Scalar Function 13
PAGE 22
To illustrate how each observer will view a scalar function under this de formation, let us consider a scalar function such as a temperature described as () = e(F, t) (3.3) Observers have to use the same temperature function with different vari ables. It is obvious that both observers will get the same temperature. Scalars do not change under a motion of a body being observed by different observers. {) (F, T) = {) (F*, t) {) (F, T) = {) (QF, t) Observers 0 and 0* will have their temperatures related by () = ()* 3.1.1.3 Internal Energy Function As another example, let us consider the internal energy function. The in ternal energy function in a simple material can be described by E = (F) According to the principle of rigid body invariance, both observers have to use the same internal energy function but with the corresponding deformation gradient for each observer. Thus, each observer will have E = (F) 14
PAGE 23
c* = E(F*) But we have F*=QF so c* = E(QF) since E is a scalar, then E = c* As a result, the internal energy function of the deformation gradient must satisfy E(F) = E(QF) for any deformation gradient F and for all orthogonal tensors Q. 3.1.1.4 Cauchy Stress Tensor The most important tensors to consider in this theory are stresses and strains. A second order stress tensor function will be used to find out how it will be views by observers 0 and 0* and also how these views of the tensor function is related to each other. Consider a Cauchy stress function as u = iT (F. T) where 15
PAGE 24
T is a scalar. Again, both observers have to use the same stress function to evaluate the stress tensor. Each observer will get a stress tensor calculated as u = ir (F, T) u = ir ( F*, T) As known, the traction vector is related to the Cauchy stress tensor by the following formula where n is the normal to the plane. Each observer will have t = un t = un t* = u*n* Since t and t* are vectors, then they are related by t* = Qt This leads to the wellknown formula 16 (3.4)
PAGE 25
Even though the deformation gradient and Cauchy stress tensor are both second order tensers, equation (3.2) which represents how the deformation gra dient is viewed by 0 and 0* is different from equation (3.4). This difference is due to the fact that the deformation gradient components are evaluated with respect to both bases. 3.1.1.5 PiolaKirchhoff Stress Tensor In addition, the first and second PiolaKirchhoff are introduced for the sake of clarity. The First PiolaKirchhoff Stress Tensor The first PiolaKirchhoff stress tensor is determined from PFT = udet (F) Since Q is an orthogonal tensorvalued function oft. det F* = det (QF) = det Qdet F = det F So P*(F*)T = u* det(F) P*(QFf = QuQT det(F) 17
PAGE 26
P*=QP The Second PiolaKirchhoff Stress Tensor The second PiolaKirchhoff stress tensor is determined from S* = (QF)1 QP S* = plp = S 3.1.2 Principle of Material Symmetry It is more clear to define this principle with an example. Consider a body, composed of elastic material, being deformed relative to a specific reference configuration. If an alternative reference configuration has been chosen and the body is subjected to the same deformation, the resultant stress will generally be different than the first one. However, many elastic materials have material symmetry that allow specific changes of the reference configuration but keeps the stress produced from an arbitrary deformation unchanged. The degree of 18
PAGE 27
the symmetry that is possessed by the material depends on the collection of such transformation [5]. Generally, most materials have a material symmetry in some form. The familiar case is isotropy. An Isotropic material owns no specific direction and it has the same properties in all direction. There exists many real materials that are isotropic. Some common isotropic materials are fluids such as air and water, metals in the polycrystalline form, concrete and others. On the other hand, transversely isotropic materials possess only one preferred direction of symme try. A composite material which is composed of a matrix supported by fibres arranged in parallel lines presents preferred directions. Material symmetries can be categorized into two types; rotation and reflectional [6]. 3.1.2.1 Rotational Symmetry Suppose a body undergoes a homogenous deformation. This deformation can be described by the equation x =(X)= FX (3.5) Now let us consider this element to have a similar deformation to the previ ous one with the exception that the whole deformation field is initially rotated about some axis. If H is an orthogonal tensor in the symmetry group, then the deformation can be then described by the equation x=FHX (3.6) The deformations described by equation (3.5) and equation (3.6) are dis tinct. If there is no material symmetry, the stress responses will be different. 19
PAGE 28
a b c Figure 3.3: Body and configuration. However, for some rotations defined by H, rotating the deformation field will result the same stress filed. This means that if the deformation described by equation (3.5) gives rise to a stress tensor u, then the deformation described by equation (3.6) gives rise to the same stress tensor u. In this case, the material has rotational symmetry characterized by H [6]. 3.1.2.2 Reflectional Symmetry Briefly, to define this symmetry, let us consider the mirror image of the deformation figure (3.3) about some plane such as X 1 = 0. This deformation can be described as x=FH'X (3.7) where H' is a reflection in the (X2 X3 ) plane and has the following com ponents 1 0 0 H' = 0 1 0 0 0 1 20
PAGE 29
Without having material symmetry, the deformations equation (3.5) and equation (3. 7) will cause the two stress to be unrelated. But, if the reflection of the deformation filed is to change the sign. then the material has reflectional symmetry[6]. a b c Figure 3.4: Body and configuration. In brief, tensors, such as the rotation tensor H and the reflection tensors H', are part of the symmetry group of the material if they define the symmetry properties of a material. The material is said to be isotropic, if its group of symmetry transformations includes all rotations about any axes and reflections about any plane. If a material has a symmetry group less than the isotropic material, then this material is called an anisotropic material. Materials that have a symmetry group that includes all rotations about only a specified axes are known as a transversely isotropic materials about that axis. A material is called orthotropic if it has only a reflectional symmetry group with respect to each one of three orthogonal planes. Wood is a wellknown example of this kind of material [ 6]. 21
PAGE 30
3.1.3 Principle of Intermediate Invariance As mentioned in the previous chapter, the deformation gradient in our theory is decomposed into elastic and plastic parts as described in equation (2.3). Let an arbitrary rigid rotation be denoted by N. The decomposition of the deformation gradient is not unique, because N can be inserted so that (3.8) To simply illustrate the concept of this principle, imagine a plastic deforma tion of a material body. At the point between the elastic and plastic deformation, which is described by the intermediate configuration, rotate the material by N and NT then impose the same elastic deformation. In other words, the mate rial is deformed plastically then rotated and rotated back before imposing the original elastic part. Since the internal energy function depends on the deformation gradient, then it can be written in terms of its elastic and plastic parts in the form The principle of the intermediate invariance requires that f(Fe,FP) = e(FeNT,NFP) for all orthogonal N 3.1.4 Principle of Elastic Symmetry 22 (3.9)
PAGE 31
According to the principle of elastic symmetry and before the existence of any plastic deformation, which means FP = I and F = pe, if the elastic part of the deformation gradients deforms as there is no change in the strain response. Where K is an orthogonal tensor in the symmetry group. Applying this principle on the internal energy function results in This means that the internal energy function possess an elastic material symmetry. Now, suppose that the existence of the plastic deformation does not affect the elastic material symmetry. This assumption leads to[7] (3.10) for all K in the elastic symmetry group. 3.2 The Free Energy Function and The Stress Response N aghdi and Trapp have developed special forms for the free energy and the stress response in finite plasticity motivated by the purely mechanical behavior of ductile metals. In this section, a brief review of this paper is introduced with a differen approach for deriving the special forms for the free energy and stress response.The stress response function is used in a matlab code to obtain a stressstrain curve. 23
PAGE 32
In the theory of elasticplastic materials, the constitutive relations consist of two sets: one relates to constitutive equation for the rate of plastic strain and the rate of workhardening. The other set of constitutive equations is for the free energy and the stress, but not of their rates. Only, the second set of constitutive equations for finite deformation of elasticplastic material are discussed here. Basically, attention is given to a special class of response functions for those constitutive equations motivated by the mechanical behavior of ductile metals at finite deformation. Attention is limited to the isothermal theory [7]. The stress response can be obtained from the Helmholtz free energy response function and its partial derivatives. An invariant representation of the Helmholtz free energy function is obtained in terms of kinematic measures. Then, the free energy response is assumed to be a quadratic function of a certain kinematic measure. The obtained constitutive equations are confined to a material that is isotropic in the reference state. The Helmholtz free energy function per unit mass '!jJ and the symmetric PiolaKirchhoff stress tensor S in the purely mechanical theory are given by (3.11) (3.12) Where k is a measure of work hardening. Equations (3.11) and (3.12) hold during loading, unloading and also neutral loading. Assuming the free energy function to depend on pe, FP and k, then the free energy function can be written as 24
PAGE 33
1/J FP, k) (3.13) Now, principles of symmetry and objectivity are applied to make further reductions on the response function 1/J. 3.2.1 Rigid Body Invariance The kinematic quantities under the superposed rigid body motions transform according to FP*+ FP k* = k According to this principle, the Helmholtz free energy function 1/J must be invariant under superposed rigid body motions. This leads to a reduction on this function and it has to satisfy FP, k) = (QFe, FP, k) (3.14) To benefit from this reduction, let us consider the use of the Cauchy theorem of the polar decomposition. According to this theorem, a nonsingular second order tensor A can be decomposed uniquely into the following products [3] A=RU=VR where R is an orthogonal tensor U and V are symmetric tensors. Since Fe is a nonsingular secondorder tensor, then it can be uniquely de composed into 25
PAGE 34
where R is an orthogonal tensor ue is a positivedefinite symmetric tensors. choosing Q = (Re)T in equation (3.14), this leads to the reduction of the free energy function and it can be written as (3.15) 3.2.2 Elastic Symmetry According to this principle, the kinematic quantities in the argument of transform as k* = k Thus, the free energy function must satisfy (3.16) By applying a theorem for the representation of a scalarvalued isotropic function of a secondorder symmetric tensor[8]. the free energy function can be described in terms of the following invariants 26
PAGE 35
thus gives the representation ( KT ee K' FP' k) = (tree' tr ( ee) 2 tr ( ee) 3 FP. k) (3.17) 3.2.3 Material Symmetry Based on this principle, the free energy function has to be invariant under the transformation F*+ FH FP*+ FPH Thus, the free energy function must satisfy ;J; (tree, tr(ee)2 tr(ee)3 FP, k) ='(tree, tr(ee)2 tr(ee)3 FP H, k) (3.18) where H is an orthogonal tensor in the symmetry group. Since FP*+ FPH then (FP*f+ HT(FP)T. To make it easier to apply some mathematical theory, the first, second and third columns of (FPf are called vectors ( v1 v2 v3). The mathematical theory of invariants states that, A scalarvalued function of vectors is isotropic if it can be expressed as a function of the inner products of vi.vi [8]. Since the whole inner product of the three vectors (v11 v2 v3 ) of the tensor (FPf can be express as Then the above restriction on the free energy function can be rewritten as 27
PAGE 36
;fJ (tree, tr(ee)2 tr(ee)3 FP (FP)T, k) =;[;(tree, tr(ee)2 tr(ee)3 BP, k) (3.19) 3.2.4 Intermediate Frame lnvariance The free energy function must be invariant under the transformation where N is an orthogonal tensor. The free energy function must satisfy ;fJ (tree,tr(ee)2,tr(ee)3,BP,k) = ;fJ (tree,tr(ee)2,tr(ee)3,NBPNT,k) (3.20) Applying the mathematical theory results [8] ;fJ (tr(ee), tr(ee)2 tr(ee)3 N BP NT, k) = ;fJ (tr(ee),tr(ee)2,tr(ee)3 (tr(BP),tr(BP)2,tr(BP)3,k) (3.21) Those symmetry implies that ;fJ (tr(ee), tr(ee)2 tr(ee)3 N BP NT, k) = (3.22) (tr(ee),tr(ee)2,tr(ee)3 (tr(BP).tr(BP)2,tr(BP)3.k) 28
PAGE 37
The free energy function can be expressed in terms of C and CP and k. so Also, it follows that Thus equation (3.22) can be written as (tr(Ce),tr(Ce)2,tr(Ce)3 (tr(IJP),tr(1JP)2,tr(1JP)3,k) == tr((CP)1C)2 tr((CP)1C)3 tr(CP), tr(CP)2 tr(CP)3 k) (3.23) Also, subtracting I which a constant and calling M == (CCP)1CI), equation (3.23) can be written as ;j; (tr(Ce), tr(Ce)2 tr(Ce)3 (tr(CP), tr(CP)2 tr(CP)3 k) == 29 (3.24)
PAGE 38
Since the free energy function is assumed to be quadratic in ( CP)1C I, it can be written in the form Po;jJ =I+ (31tr ((CP)1 CI)+ o:1tr ((CP)1 CIr (3.25) + o:2 [tr ( (CP)1 CI) r From equation (3.24) and (3.11), the following equation of the stress can be obtained or Where s = 2(3l(CP)1 + 8o:l(CP)1[EEP](CP)1 (3.26) + 8o:2tr[(CP)1 (EEP)J (CP)1 s = SP8o:l(CP)1[EEP](CP)1 +8o:2tr[(CP)1 (EEP)J (CP)1 (3.27) This result was based on the assumption that the response function for the free energy is elastically isotropic. An alternative constitutive equation for the stress can be obtained by assuming that the response for the difference S SP is elastically isotropic but the residual stress SP will not possess this property [7]. This leads to 30
PAGE 39
SSP= 8o:1(CP)1[EEP](CP)1 (3.28) + 8o:2tr[(CP)1 (EEP)] (CP)1 Assuming that the coefficients o:1 and o:2 are unaffected by plastic deformation and considering the situation in which Ee = E, CP = I, then (3.28) becomes SSP= 8o:1E + 8o:2 (trE) I In this special case it is possible to make Where J.L, and A are constants. 1 0:2 =A 8 The constitutive equation (3.28) can be written as SSP= A tr [ (CP)1 (EEP) J (CP)1 + 2J.L (CP)1[EEP](CP)1 3.3 Modified Simo's Yield Function (3.29) (3.30) Simo derived a formula based on the Von Mises yield criterion. Von Mises criterion states that yielding will take place when the second invariant of the stress deviator tensor reaches a critical value which is the value of the second invariant at yield in simple tension[1]. The classical Von Mises yield condition in stress space ignoring the hardening effects can be written as 31
PAGE 40
j(T) = 1 tr(11)3(tr(1))2R Where 1 is the Kirchhoff stress tensor defined in equation (2.9). R is the radius of the Von Mises Sphere. (3.31) Using Von l\1ises yield condition, Simo has obtained the following formula for the yield function[9] j(S, C)= Jtr(SCSC)R (3.32) Now using the intermediate configuration P instead of the current configuration, the Kirchhoff stress tensor can be written as Applying the Von Mises condition based on the intermediate configuration leads to J(S, FP) = The above formula can be rewritten as f(S, CP) = Jtr(SCPSCP)(tr(SCP))2 R This formula is obtained based on the Von Mises criterion for isotropic hardening. Isotropic hardening means that the yield surface will not move and since R is constant then also it will not expand. On the other hand, kinematic 32
PAGE 41
hardening means the yield surface will move due to the Bauschinger effect [1]. In our case, kinematic hardening will be included. This means that the residual stress Sp is subtracted from S which leads to the following formula f((SSP), FP) = (3.33) Jtr((SSP)CP(SSP))CP)(tr((SSP))CP))2 R 33
PAGE 42
4. Work hardening, Numerical scheme and results 4.0.1 Measures of Work Hardening and State Equations 4.0.1.1 Measures of Work Hardening The plastic work is an important concept in the theory of plasticity. The work done per unit volume on an element during straining is given by[1] dW = ude ( 4.1) = u(dee + deP) = dWe+dWP Since the elastic deformation is reversible then the elastic energy ( dWe = udee) is recoverable. On the other hand the plastic energy can not be recovered since the plastic deformation is permanent. Thus the rest of equation ( 4.1) is called the plastic work per unit volume and defined by (4.2) or dWP = SdEP Two measures of work hardening have been proposed. The first assumes that the amount of hardening depends on the total plastic work only. This implies that any resistance to further yielding depends on the amount of work that has been done on the material. This amount of work is measured by the yield criterion which can be written as 34
PAGE 43
f(u) = k ( 4.3) Where K is known as the work hardening parameter. For work hardening, K keeps changing, but for isotropic hardening f(u) remains constant. The second hypothesis assumes the plastic strain as a measure of work hardening. In this case, the yield criterion is written as (4.4) Drucker work hardening assumption can be expressed in the following mathematical formula (4.5) The general stressstrain relations can be obtained based on the following two assumptions 1The existence of a loading function at each stage of the plastic deformation. This means that any further plastic deformation can take place only if f(u) > k. 2The relation between infinitesimal stress and plastic strain is linear which means[l) d(cP) = Cd(u) Where C is a fourth order tensor. 35
PAGE 44
This linear relation tells us that even though it is a rate equation, it IS independent of time. 4.0.1.2 State Equation Naghdi and Trapp have obtained some results that governing the nonlinear behavior of elasticplastic contained in the work of Green and Naghdi. They assume the existence of plastic strain specified by a secondorder tensor EP at each point of the continuum and a measure of workhardening specified by the scalar function k. Based on work hardening on strain space, an inequality was obtained. This inequality holds for any smooth closed strain cycle beginning and ending at E0 inside the loading surface in the strain space[lO]. The inequality can be written as The above work inequality has an important consequence which leads to the following formula (4.6) The workhardening is given the form (4.7) Where H IS a secondorder tensor function depends on the variables S, Ep, k. A general matlab code has been developed to solve the differential forms of those equations. This code starts with a matrix of inputs to the increment 36
PAGE 45
of the right CauchyGreen tensor and the symmetric PiolaKirchhoff stress tensor. Based on those input increments, the yield function is calculated to check whether the deformation is in the elastic region or the plastic region. No differential equations are computed for the elastic region. However, there are a set of differential equations to be solved in the plastic region. Those equations are rate equations but independent of time since each term of them has time increment. These equations are 1) Rate form of equation (3.12) as as as. ac c + acp CP + ak k s = 0 2) Equation ( 4.6) as as. aJ aCPCP+ akk+roC =O 3) Rate form of yield function 4) A postulated rate equation for hardening parameter k 4.0.2 Numerical Scheme 1 ktr(CP) = 0 2 (4.8) (4.9) (4.10) ( 4.11) Given the increment in C (or some combination of CandS), a set of differential equation are solved for the increment of C, S, CP, k and I Those differential equations are summarized below. 37
PAGE 46
(of )n(cn+l _en)+ ( of )n((CP)n+l (CP)n) +(of t(kn+l kn) oc oCP ok = f(Cn' (CP)n, kn) (4.14) (4.15) This code shows the results in some plots of the stresses and strains. The code is used also to test and show the validity of the derived stressstrain equation (3.30) and the modified Simo's yield function (3.33). 4.0.3 Results It is necessary to find a formula for the residual stress SPin equation (3.30) to obtain some results. Obtaining a valid formula has not been done yet. Preliminary constitutive equations have been used that were guided by the experimental results in [11]. Those experimental results are used to obtain some values of j3 and fit them into a formula. Figure ( 4.1) shows a rough curve for the residual stress Sp and the plastic strain cP. In the following procedure, the residual stress Sp is thought to represent the center of the yield surface. For this reason a middle point on the elastic stressstrain curve is used to calculate the residual stress and it is given the name Spc In figure (4.1), E is the Young's modulus. Obtaining the residual 38
PAGE 47
stress SP and the strain EP from the experimental results, the strain tensor at the middle point of the elastic stressstrain curve is calculated from ( 4.16) Then CP(1, 1) is obtained from the following formula CP(1, 1) = 2EP(1, 1) + 1 ( 4.17) /3 is calculated from the equation The calculated values of /3 and CP(1, 1) are plotted as shown below The following formula for /3 is obtained by using a polynomial curve fitting j3 = 108(1.1858 + 2.3527CP(1, 1)1.1669(CP(1, 1))2 ) (4.18) The experimental results of the residual stress at the middle point and the calculated residual stress from equation (4.18) are plotted verses CP(1, 1) and the results are shown below The first six plots in figure ( 4.4), show the data that will be input to the code. In this test. the material is loaded in a simple tension test until point b. Then, the load is removed and the material is reloaded in the reverse direction until point d. Figure ( 4.6) shows the stressstrain curve for this test. In figure 39
PAGE 48
( 4. 7), the initial and subsequent yield surfaces are plotted. It can be seen how the yield surface translates in the direction of loading. 4.1 Conclusion A specific general form, equation (3.30), of the stressstrain law in large deformation has been derived using the general invariant principles and elastic symmetry arguments. A general matlab code has been developed to integrate the complete set of nonlinear elasticplastic constitutive equation with large de formation. Simo's yield criterion has been modified to include kinematic hard ening so the yield surface can translate. The results appear reasonable and the general integration scheme can now be used to test any general experimentally developed constitutive equations. 40
PAGE 49
s E Figure 4.1: The stressstrain curve. 41
PAGE 50
curve fits + data linear til quadratic fit cubiclit 10000 __ __:.::::.;c:_c;.;__J 8000 2000 +' 1.002 1.004 1.006 1.008 1.01 1.012 1.014 Cp(1.1) data Figure 4.2: Beta and the right CauchyGreen tensor curve. curve fits 80001 I 1 I =:== I 1 I 6000 l / / r /' 5000 :::4000 r /li I I/ 2000 1 _/ 1000 // ! 1 1.002 1.004 1.006 1.008 1.01 1.012 1.014 Cp(1,1) Figure 4.3: The residual stress and the right CauchyGreen tensor curves. 42
PAGE 51
1.031 11 "'/_/_/ __ _Jj 0 100 200 300 400 0 100 Nstep 200 Nstep 300 0 1 00 200 300 400 0 1 00 200 300 Nstep Nstep
PAGE 52
i I Ul 2ll C; "' I I c: .!! II "' 1 "' = ,, 0 II .r:; .r:; 1:! 'i I I .. 0 I a: qI I I 3 0 a SlressSirain curve d I I i 1 I 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 AlmansiHemel slrain tensor E(1,1) Figure 4.6: StressStrain Curve. 44
PAGE 53
x 11!1tial yield surface at point"a" x 1 o"field surface at point"b" 6 6 4 4 2 C) 2 C) N 0 N u; r.n 0 2 2 4 4 6 6 5 0 5 5 0 5 S11 X 104 S11 X 104 x 1 o"field surface at point"c" x 10"field surface at point"d" 6 6 4 4 2 C) 2 C) N 0 N u; u; 0 2 2 4 4 6 6 5 0 5 5 0 5 S11 X 104 S11 X 104 Figure 4. 7: The initial and subsequent yield surface 45
PAGE 54
APPENDIX A. Solving Linear Equations Using Some Functions 0001 function Cpold = plastic2; 0002 clear all;close all;clc 0003 % the purpose of this code is to compare & simulate 0004 % different constitutive equations for elasticplastic 0005% materials. 0006 %0007 % properties 1s a function that has properties of 0008% d1fferent materials. 0009% typical properties of some metals are used from:0010 % (Engineering materials science, No 615917). 0011 %these properties are:0012% ME =the modulus of elasticity cr Young's modulus. 0013% MEW =the shear modulus or the modulus of rig1d1ty. 0014% NEW =Poisson's ratio. 0015 % SEGMA =yield strength. 0016 % LAlvJDA lame modulus. 0017 %0018 % STRESS !JO = #, is the number of the st:cess fun::t1on 0019 % that will be used. 0020 ,, /,. Stress h:nct1ons are ln a funct1on 1f i le called 0021 % stressfunc. 46
PAGE 55
0022% YIELDF NO = #, is the number of the yield function that 0023% will be used. 0024% Yield functions are in a function Mfile called y:eldfunc. 0025 %0026 format compact 0027 format long e 0028 STRESS_NO 1; 0029 YIELDF ...NO 1; 0030 properties 0031 global ME MEW NEW SEGMA LAMDA 0032 global STRESS...NO YIELDF...NO 0033 global NSTEP i M N 0034 %0035 % F = the deformation gradient which is given & increases 0036 % by delta. 0037 % c = the right CauchyGreen tensor, Cp the plastic 0038 % part of it. 0039 ,, E = the AlmansiHemel strain tensor, Ep the plastic /, 0040 % part of it. 0041 % s = the second PiolaKnchhoff stress tensor,Sp the 0042 % plastic part of it. 0043 %. P1ola kirrchoff = the firs PiolaKirchhoff stress tensor. 0044 % Kappc.. the work ha:cden1ng parameter. 0045 '/, alpha 47
PAGE 56
0046 % vield f = vield criterion (Von ses yield condi:ion). 0047 % yield f < 0 ,for elastic deformatJ.cn domain. 0048 % yield f c ,fer plaStlC deformatio:1 domain. 0049 % signa cauchy = the Ca":Jchy stress tensor. 0050 %0051 Cdot = zeros(3,3); 0052 Cpdot = zeros(3,3); 0053 Sdot = zeros(3,3); 0054 7, H is an orthogonal tensor ir: the symmetry group; 0055 H zeros(3,3); 0056 dt = 1. 0; 0057 Kappaold = SEGMA; 0058 Cold = eye(3,3); 0059 Cpold = eye(3,3); 0060 Sold = zeros(3,3); 0061 I= eye(3,3); 0062 NSTEP = 0.0; 0063 for i = 1:150 0064 %0065 NSTEP = NSTEP + 1; 0066 if i <= 100 0067 0068 0069 delta 0.0001; taw 0.0; else if i <= 250 48
PAGE 57
0070 delta = 0.001; 0071 taw = 0.0; 0072 else 0073 delta = 0.001; 0074 taw = 0.0; 0075 end 0076 %0077 W1 = eye(9,9); 0078 W1(5,5) = 0.0; 0079 W1(9,9) = 0.0; 0080 W2 = zeros(9,9); 0081 W2(5,5) = 1.0; 0082 W2(9,9) = 1.0; 0083 DSDC = numdifftwo(Cold, Cpold, Kappaold); 0084 W4 = eye(9,9); 0085 W = [ W1 W2 ; DSDC W4 ] ; 0086 Z = [delta 0 0 0 0 0 0 0 0 0 0 0 0 taw ... 0087 0088 0089 0090 0091 0092 0093 unknowns Cdot(1,1) Cdot(1,2) Cdot(1 ,3) Cdot(2,1) Cdot(2,2) 0 = 0 0 taw] ; inv(W)*Z'; = unknowns(1); Sdot (1, 1) = unknowns(10); = unknowns(2); Sdot (1, 2) = unknowns(11); = unknowns(3); Sdot(1,3) = unknowns (12); = unknowns(4); Sdot(2,1) = unknowns (13); = unknowns(5); Sdot(2,2) = unknowns (14) ; 49
PAGE 58
0094 Cdot(2,3) = unknowns(6); Sdot(2,3) = unknowns (15); 0095 Cdot(3,1) = unknowns(?); Sdot(3,1) = unknowns(16); 0096 Cdot(3,2) = unknowns(8); Sdot(3,2) = unknowns (17); 0097 Cdot(3,3) = unknowns(9); Sdot(3,3) = unknowns (18); 0098 Cnew = Cold + Cdot*dt; Snew = Sold + Sdot*dt; 0099 %0100 % fcllov;ing is a funct1cn called yieldfunc that lS ... 0101 % calculating the yield function. 0102 yield_f = yieldfunc(Cnew, Cpold, Kappaold); 0103 % follov:ing is If Statement to check whether the yield f 0104 % is in the elastic domain or in the plastic domain. 0105 if yield_f < 0.0 % elastic domain, 0106 0107 'inside' ytest (i) = 1. 0; 0108 elseif yield_f >= 0.0 % plastic doma1n, 0109 0110 'at or outside yield surface ytest(i) = 2.0; 0111 %0112 % numdiffone lS a function that v:ill calculate 0113 % DyDS = the der1vative of the peld function 0114 % "'. r. t the secc,nd PlolaKirchhcff stress tenser, 0115 % DyDC = the der:J. vati '.'e of the peld functlon ''. r. t 0116 % the right CauchvGreen tensor, 0117 % DyDCp = the der1vative of the yield fur.cticn ''. r. t 50
PAGE 59
0118 % 0119 % 0120 % 0121 % 0122 % 0123 % the plastic part of the right CauchyGreen tensor, DyDKappa = the derivative of the yield function kappa, those partial are used in the following formula:(dyield f/dC)*(dC/dt) + Cdyield f/dCp)*(dCp/dt) + (dyield f/dKappa)*(dKappa/dt) = 0; 0124 %0125 [DyDC, DyDCp, DyDKappa] = numdiffone(Cold, Cpold, Kappaold); 0126 %0127 % 0128 % 0129 % 0130 % 0131 % 0132 % 0133 % 0134 % 0135 % 0136 % 0137 % 0138 % 0139 % 0140 % 0141 % numdifftwo is a function that will calculate DSDC = the derivative of the second PielaKirchhoff stress tensor w.r.t the right CauchyGreen tenser, DSDCp = the derivative of the second PielaKirchhoff stress tensor w.r.t the plastic part of the right CauchyGreen tensor, DSDKappa = the derivative of the second PielaKirchhoff stress tensor w.r.t kappa, those partial derivatives are used in the following formulas:(dS/dCp)*(dCp/dt) + (dS/dKappa)*(dKappa/dr) + gamma*(dyield f/dC) G.O; (dS/dt) (dS/dC)*(dC/dt; (dS/dCp)(dCp/dt) (dS/dKappa)*(dKappa/dt) 0.0 51
PAGE 60
0142 %0143 [DSDC, DSDCp, DSDKappa] = numdifftwo(Cold, Cpold, Kappaold); 0144 %0145 % 0146 % 0147 % 0148 % 0149 % 0150 % 0151 % 0152 % 0153 % 0154 % set the follow1ng equations :(dyield f/dCp)*(dCp/dt) + (dyield f/dKappa)* ... (dKappa/dt) =(dyieldf/dC)*(dC/dt) ; ..... (1); fdS/dCp)*(dCp/dt) + (dS/dKappa)*(dKappa/dt) + gamma*(dyield f/dC) = 0.0; .... (210); dkappa/dt(1/2)trace[dCp/dt] = 0.0; ...... (11); to the following mitrx form:[A] [d/dt and gamma] = [rhs] ; where [ A ] = 1s a matrix of the coefficients; 0155 %0156 W3 = zeros(9,11); 0157 W5 = zeros(1,9); 0158 W6 = zeros(9,1); 0159 W7 = zeros(9,9); 0160 DSDK = DSDKappa(1:9); 0161 DyDC1 = DyDC(1:9); 0162 DyDCp1 = DyDCp(1:9); 0163 A = [ W1 W2 W3; DSDC W4 DSDCp DSDK' W6; W7 W7 DSDCp_ .. _. 0164 DSDK' DyDC1'; DyDC1 W5 DyDCp1 DyDKappa 0; W5 ... 0165 W5 1/2 0 0 0 1/2 0 0 0 1/2 1 0 ] ; 52
PAGE 61
0166 yield_f = yieldfunc(Cold, Cpold, Kappaold); 0167 S = stressfunc(Cold, Cpold, Kappaold); 0168 diffS = S Sold; 0169 diffS = diffS'; 0170 diffS1 = diffS(1:9); 0171 Z1 = Z(1:9); 0172 B = [ Z1 diffS1 W5 yield_f 0 ] ; 0173 unknowns = inv(A)*B'; 0174 %0175 0176 0177 0178 0179 0180 0181 0182 0183 0184 0185 0186 0187 0188 0189 Cdot (1, 1) Cdot (1, 2) Cdot (1 ,3) Cdot(2,1) Cdot(2,2) Cdot(2,3) Cdot(3,1) Cdot(3,2) Cdot(3,3) Cpdot (1, 1) Cpdot(1,2) Cpdot(1,3) Cpdot(2,1) Cpdot(2,2) Cpdot(2,3) = unknowns(1); = unknowns(2); = unknowns(3); = unknowns(4); = unknowns(5); unknowns(6); unknowns(?); = unknowns(8); = unknowns(9); = unknowns (19); = unknowns(20); = unknowns(21); = unknowns(22); = unknowns(23); = unknowns(24); 53 Sdot (1, 1) = unknowns(10); Sdot(1, 2) = unknowns (11) ; Sdot(1 ,3) = unknowns (12) ; Sdot(2,1) = unknowns (13); Sdot(2,2) = unknowns (14); Sdot(2,3) = unknowns (15); Sdot(3,1) = unknowns (16); Sdot(3,2) = unknowns (17); Sdot(3,3) = unknowns (18) ;
PAGE 62
0190 0191 0192 0193 0194 0195 0196 % 0197 0198 0199 0200 0201 0202 Cpdot(3,1) unknowns(25); Cpdot(3,2) unknowns(26); Cpdot(3,3) = unknowns(27); Kappadot = unknowns(28); gamma= unknowns(29); Cnew = Cold + Cdot*dt; Cpdot = (1/2)*(Cpdot + Cpdat'); Cpnew = Cpold + Cpdot*dt; Snew = Sold + Sdot*dt; Kappanew = Kappaold + Kappadot*dt; gamma; Cpold = Cpnew; Kappaold = Kappanew; 0203 end 0204 %0205 Cpp(:, :,i)= Cpold; 0206 %0207 Cold = Cnew; 0208 Sold = Snew; 0209 %0210 [V,D] = eig(Cold); 0211 Q = sqrt(D); 0212 U = V*Q*inv(V); 0213 S = Snew; 54
PAGE 63
0214 % Pidak S*F'; 0215 % s1gnac (1/det(F))*F*S*F'; 0216 E = (1/2)*(CnewI); 0217 % Defining variables tc be plotted 0218 E11(i) = E(1,1); E12(i) =E(1,2); E13(i) = E(1,3); 0219 E21(i) = E(2,1); E22(i) = E(2,2); E23(i) = E(2,3); 0220 E31(i) = E(3,1); E32(i) = E(3,2); E33(i) = E(3,3); 0221 CP11(i) = Cpold(1,1); CP12(i) = Cpold(1,2); 0222 CP13(i) = Cpold(1,3); 0223 CP21(i) = Cpold(2,1); CP22(i) = Cpold(2,2); 0224 CP23(i) = Cpold(2,3); 0225 CP31(i) = Cpold(3,1); CP32(i) = Cpold(3,2); 0226 CP33(i) = Cpold(3,3); 0227 811(i) = 8(1,1); 812(i) = s (1, 2); 813(i) =S(1,3); 0228 821(i) =8(2,1); 822(i) = 8(2,2); 823(i) = 8(2,3); 0229 S31(i) = S(3,1); 832(i) = S(3,2); S33(i) = S(3,3); 0230 U11 (i) =U(1,1); U12(i) =U(1,2); U13(i) =U(1,3); 0231 U21(i) = U(2,1); U22(i) = U(2,2); U23(i) = U(2,3); 0232 U31(i) = U(3,1); U32(i) = U(3,2); U33(i) = U(2,3); 0233 %0234 subplot(441), plot(E11), grid 0235 xlabel 'NSTEP' 0236 ylabel 'E11' 0237 subplot(442), plot(E22), grid 55
PAGE 64
0238 xlabel 'N8TEP' 0239 ylabel 'E22' 0240 subplot(443), plot (CP11), grid 0241 xlabel 'N8TEP' 0242 ylabel 'CP11' 0243 subplot(444), plot(CP22), grid 0244 xlabel 'NSTEP' 0245 ylabel 'CP22' 0246 subplot(445), plot(U11), grid 0247 xlabel 'N8TEP' 0248 ylabel 'U11' 0249 subplot(446), plot(U22), grid 0250 xlabel 'N8TEP' 0251 ylabel 'U22' 0252 subplot(447), plot (S11), grid 0253 xlabel 'N8TEP' 0254 ylabel '811' 0255 subplot(448), plot(S22), grid 0256 xlabel 'NSTEP' 0257 ylabel '822' 0258 subplot(449), plot(E11,S11), grid 0259 xlabel 'Ell' 0260 ylabel ) 811) 0261 subplot(4,4,10), plot(E12,S12), grid 56
PAGE 65
0262 0263 0264 0265 0266 0267 0268 0269 0270 0271 0272 0273 % 0274 % 0275 % 0276 % 0277 % 0278 % 0279 end 0280 pause xlabel 'E12' ylabel 'S12' subplot(4,4,11), plot(E22,S22), grid xlabel 'E22' ylabel 'S22' subplot(4,4,12), plot(ytest), grid xlabel 'NSTEP' ylabel 'ytest' hold on subplot(4,4,12), plot([12],[1 2],'.') hold off subplot(441), plot(PKll), grid xlabel 'NSTEF' ylabel 'PK11' subplct(442), plot(PK22), grid xlabel 'HSTEP' ylabel 'PK22' 0281 DONE= plotting(Cpp); A.l Properties Function 0001 function properties 0002% typical properties of s0me metals used from: 57
PAGE 66
0003 %CEng1neering materials science, JJo 615917). 0004 %these properties are:0005% ME =the modulus of elasticity or Young's modulus. 0006 %MEW =the shear modulus or the modulus of rigidity. 0007% NEW =Poisson's ratio. 0008 % SEGMA =yield strength. 0009 % LAlvJDA lame modulus. 0010% Units:1psi 6.895 kPa. .. and 1ksl 1000 psi. 0011 global ME MEW NEW SEGMA LAMDA 0012 %0013% steel 0.2%C, hotrolled,(Epsi), (mewpsi). (segnayksi): 0014 ME = 30*106*6.895; 0015 MEW 12*106*6.895; 0016 NEW 0.27; 0017 SEGMA = 40*1000*6.895; 0018 LAMDA = (2*MEW*NEW)/(1 2*NEW)*6.895; 0019 %0020 % steel 0.2%C, coldrolled, (Epsi), (mewpsi), (segna yksi): 0021 % ME = 30*106; 0022 % !'1E'I'f 12*106; 0023 % !JEw 0.27; 0024 I. SEGIA 6.5*1000; 0025 % LAI'JDA (2*11EW*IJE\0 I C 1 2*NE'vv'),.6.895; 0026 %58
PAGE 67
0027 % ah:.rninum alley, 202414, (Epsl), Cmev.1psi), (segna yksi) : 0028 % HE = 10.6*106; 0029 % 4*1C6; 0030 % 0.33; 0031 % SEGMA 44*1000; 0032 % LA!1DA (2*HEW*IJEW) I ( 1 2*NEW)*6.895; 0033 %0034 % copper annealed,(Epsl), (mevcpsi), ( segna yksi) : 0035 % HE = 17*106; 0036 % HEW 6.4*106; 0037 % !JEW 0.33; 0038 % SEG!1A 10*1000; 0039 % LAHDA (2*MEW*NEW)/(: 2*NEW)*6.895; A.2 Stress Function 0001 % th1s is a function called stressfunc that is calculating .. 0002 % the stress tensor, Pick which formula that w1ll be used .. 0003 % to find the stress tensor; 0004 function Stressf = stressfunc(C, Cp, Kappa); 0005 global ME MEW NEW SEGMA LAMDA 0006 global STRESS_NO YIELDF_NO 0007 Sp = [0 0 0 ; 0 0 0 ; 0 0 0] ; 0008 if STRESS_NO == 1.0 0009 % the f1rst stress funcion is from a paper called ... 59
PAGE 68
0010 0011 0012 0013 0014 0015 0016 0017 % 0018 else 0019 0020 0021 end % (On Finite ElasticPlastic Deformation of Metals BY ... % P. M. Naghadi & J. A. March, 1974, Pg 259, ... % Equation .. (63) which has the form % Sress(Sp,C,Cp,lame,mew) = Sp + lame*tdinv(Cp) ... % *(CCp)]*inv(Cp)+mew*lnv(Cp)*(CCp)*inv(Cp); Stressf = Sp + 0.5*LAMDA*trace(inv(Cp)*(CCp))* inv(Cp)+ MEW*inv(Cp)*(CCp)*inv(Cp); Stress = Stressf/MEW; stress = 'No stress function has been chosen' pause A.3 Yield Function 0001 % this 1s a function called yieldfunc that 1s calculating .. 0002 % the yield function, Pick wh1ch formula that will be used .. 0003% to find yield f; 0004 function yieldf = yieldfunc(C, Cp, Kappa); 0005 global STRESS_NO YIELDF_NO 0006 global NSTEP 0007 ternS = stressfunc(C, Cp, Kappa); 0008 if YIELDF_NO == 1 0009 % the first yield func1on 1s frcm a paper called (A 0010 % FRAMWDRK FOR FIIJITE STRAIN ELASTOPLASTICITY BASED ON. 60
PAGE 69
0011 % 11AXI11IUJ1 PLASTIC DISSIPATION A!iD THE MULTILICATIVE .. 0012 % DECOMPOSITION: PART I. CONTINUUM FORMULATION BY J.C ... 0013 % SIMO ,26 NOVEMBER 1986,, Page 214, Equation .. (2.14b) .. 0014 % has the form 0015 % yield f(S,C) = sqrt(S(i,J)*S(k,l)*C(l,k)*C(J ,1) 0016 % 1/3*(S(k,l)*C(k,ll)2)R 0017 yieldf = sqrt(trace(temS*C*temS*C) (1/3)_ .. _. 0018 *(trace(temS*C))2) Kappa; 0019 % y1eldf = sqrt(trace(temS*Cp*temS*Cp) 0020 % (1/3)*(trace(temS*Cp))2)Kappa; 0021 else 0022 yieldfunc = 'No yield function has been chosen' 0023 pause 0024 end A.4 Differentiation of The First Equation 0001 % this is a function called numdiffone that is calculating .. 0002% the part1al derivative of the v1eld function C, .. 0003 % Cp, Kappa) numerically are needed to be used in 0004% the following fJrmula :0005% (dyield f/dC)*C + (dyleld f/dCpJ*Cp + (dyield f/dKappa)* ... 0006 % Kappa = 0; 0007 function [DyDC, DyDCp, DyDKappa] numdiffone(Cold, Cpold,_ .. _. 0008 Kappaold) ; 61
PAGE 70
0009 global NSTEP 0010 y = yieldfunc(Cold, Cpold, Kappaold); 0011 %0012 % first determining DyDC = the derivat1ve of the yield 0013% w.r.t the right CauchyGreen tensor, 0014 DyDC = zeros(3,3); 0015 for 1 = 1:3 0016 fork= 1:3 0017 diff = zeros(3,3); 0018 if Cold(l,k) = 0.0 0019 diff(l,k) = 1.0e06*Cold(l,k); 0020 else diff(l,k) 1.0e06; end 0021 0022 0023 0024 ydiff = yieldfunc(Cold + diff, Cpold, Kappaold); DyDC(l,k) = ( ydiffy )/diff(l,k); 0025 end 0026 end 0027 %0028 % s8cond determining DyDCp = the derivative of the vield 0029 % function w.r.t the plastic part ofthe right 0030 % CauchyGreen tensor, 0031 DyDCp = zeros(3,3); 0032 for 1 = 1:3 62
PAGE 71
0033 fork= 1:3 0034 diff = zeros(3,3); 0035 if Cpold(l,k) = 0.0 0036 diff(l,k) = 1.0e06*Cpold(l,k); 0037 else 0038 diff(l,k) = 1.0e06; 0039 end 0040 ydiff = yieldfunc(Cold, Cpold + diff, Kappaold); 0041 DyDCp(l,k) = ( ydiffy )/diff(l,k); 0042 end 0043 end 0044 %0045 % finally determining DyDKappa = the der1vative of the 0046 % y1eld function Kappa, 0047 diff = 1.0e06*Kappaold 0048 0049 ydiff = yieldfunc(Cold, Cpold, Kappaold + diff); DyDKappa = ( ydiff y )/diff; A.5 Differentiation of The Second Equation 0001 % this 1s a function called that is calculating ... 0002 % the partial derivative of the second PielaKirchhoff 0003% stress tensor w.r.t (C. Cp, numericallv. 0004 % which are needed to be used 1n the fcllow1ng formulas 0005% (dS/dCp)* dCp/dt) + tdS/dKappa:*ldVappa/dt) ... 63
PAGE 72
0006% Cdy1eld f/dC) = 0.0; 0007% (dS/dt)"dS/dC)*CdC/dt)(dS/dCp)*(dCp/dt)(dS/dKappa)* ... 0008 % (dKarpa/dt) = 0.0: 0009 function [DSDC, DSDCp, DSDKappa] = numdifftwo(Cold, Cpold,_ .. _. 0010 Kappaold); 0011 SS = stressfunc(Cold, Cpold, Kappaold); 0012 % First determining DSDC = the derivative of the second 0013 % PlolaKirchhoff stress tensor the right Cauchy... 0014 %Green tensor, 0015 DSDC = zeros(3,3); 0016 for 1 = 1:3 0017 for k = 1:3 0018 diff = zeros(3,3); 0019 if Cold(l,k) = 0.0 0020 diff(l,k) = 1.0e06*Cold(l,k); 0021 else 0022 diff(l,k) 1.0e06; 0023 end 0024 SSdiff = stressfunc(Cold + diff' Cpold, Kappaold); 0025 DSDC(: ,: ,l,k) = ( SSdiff ss )/diff(l,k); 0026 end 0027 end 0028 DSDC= [DSDC(1,1,1,1) DSDC(1,1,1,2) DSDC(1,1,1,3) 0029 DSDC(1,1,2,1) DSDC(1,1,2,2) DSDC(1,1,2,3) 64
PAGE 73
0030 DSDC(1,1,3,1) DSDC(1,1,3,2) DSDC(1,1,3,3); 0031 DSDC(1,2,1,1) DSDC(1,2,1,2) DSDC(1,2,1,3) 0032 DSDC(1,2,2,1) DSDC (1, 2, 2, 2) DSDC(1,2,2,3) 0033 DSDC(1,2,3,1) DSDC(1,2,3,2) DSDC (1 2 3 3) ; 0034 DSDC(1,3,1,1) DSDC(1,3,1,2) DSDC(1,3,1,3) 0035 DSDC(1,3,2,1) DSDC(1,3,2,2) DSDC(1,3,2,3) 0036 DSDC(1,3,3,1) DSDC(1,3,3,2) DSDC (1 3 3 3) ; 0037 DSDC(2,1,1,1) DSDC(2,1,1,2) DSDC(2,1,1,3) 0038 DSDC(2,1,2,1) DSDC(2,1,2,2) DSDC(2,1,2,3) 0039 DSDC(2,1,3,1) DSDC(2,1,3,2) DSDC(2,1,3,3); 0040 DSDC(2,2,1,1) DSDC(2,2,1,2) DSDC(2,2,1,3) 0041 DSDC(2,2,2,1) DSDC(2,2,2,2) DSDC(2,2,2,3) 0042 DSDC(2,2,3,1) DSDC(2,2,3,2) DSDC(2,2,3,3); 0043 DSDC(2,3,1,1) DSDC(2,3,1,2) DSDC(2,3,1,3) 0044 DSDC(2,3,2,1) DSDC(2,3,2,2) DSDC(2,3,2,3) 0045 DSDC(2,3,3,1) DSDC(2,3,3,2) DSDC(2,3,3,3); 0046 DSDC(3,1,1,1) DSDC(3,1,1,2) DSDC(3,1,1,3) 0047 DSDC(3,1,2,1) DSDC(3,1,2,2) DSDC(3,1,2,3) 0048 DSDC(3,1,3,1) DSDC(3,1,3,2) DSDC(3,1,3,3); 0049 DSDC(3,2,1,1) DSDC(3,2,1,2) DSDC(3,2,1,3) 0050 DSDC(3,2,2,1) DSDC(3,2,2,2) DSDC(3,2,2,3) 0051 DSDC(3,2,3,1) DSDC(3,2,3,2) DSDC(3,2,3,3); 0052 DSDC(3,3,1,1) DSDC(3,3,1,2) DSDC(3,3,1,3) 0053 DSDC(3,3,2,1) DSDC(3,3,2,2) DSDC(3,3,2,3) 65
PAGE 74
0054 DSDC(3,3,3,1) DSDC(3,3,3,2) DSDC(3,3,3,3)]; 0055 %0056 % second determining DSDCp = the derivative of the 0057 % second PiolaKirchhoff stress tensor w.r.t the 0058 %plastic part of the r1ght CauchyGreen tensor, 0059 DSDCp = zeros(3,3); 0060 for l = 1:3 0061 for k = 1:3 0062 diff = zeros(3,3); 0063 if Cpold(l,k) = 0.0 0064 diff(l,k) = 1.0e06*Cpold(l,k); 0065 else diff(l,k) 1.0e06; end 0066 0067 0068 0069 0070 SSdiff stressfunc(Cold, Cpold + diff, Kappaold); DSDCp(:,: ,l,k) = ( SSdiffSS )/diff(l,k); 0071 end 0072 0073 0074 0075 0076 0077 end DSDCp=[DSDCp(1,1,1,1) DSDCp(1,1,2,1) DSDCp(1,1,3,1) DSDCp(1,2,1,1) DSDCp(1,2,2,1) DSDCp(1,1,1,2) DSDCp(1,1,1,3) DSDCp(1,1,2,2) DSDCp(1,1,2,3) DSDCp(1,1,3,2) DSDCp (1 1 3 3) ; DSDCp (1 2 1 2) DSDCp(1,2,1,3) DSDCp(1,2,2,2) DSDCp(1,2,2,3) 66
PAGE 75
0078 DSDCp(1,2,3,1) DSDCp(1,2,3,2) DSDCp (1 2 3 3) ; 0079 DSDCp(1,3,1,1) DSDCp(1,3,1,2) DSDCp (1, 3, 1, 3) 0080 DSDCp(1,3,2,1) DSDCp(1,3,2,2) DSDCp(1,3,2,3) 0081 DSDCp(1,3,3,1) DSDCp(1,3,3,2) DSDCp(1,3,3,3); 0082 DSDCp(2,1,1,1) DSDCp(2,1,1,2) DSDCp(2,1,1,3) 0083 DSDCp(2,1,2,1) DSDCp(2,1,2,2) DSDCp(2,1,2,3) 0084 DSDCp(2,1,3,1) DSDCp(2,1,3,2) DSDCp(2,1,3,3); 0085 DSDCp(2,2,1,1) DSDCp(2,2,1,2) DSDCp(2,2,1,3) 0086 DSDCp(2,2,2,1) DSDCp(2,2,2,2) DSDCp(2,2,2,3) 0087 DSDCp(2,2,3,1) DSDCp(2,2,3,2) DSDCp(2,2,3,3); 0088 DSDCp(2,3,1,1) DSDCp(2,3,1,2) DSDCp(2,3,1,3) 0089 DSDCp(2,3,2,1) DSDCp(2,3,2,2) DSDCp(2,3,2,3) 0090 DSDCp(2,3,3,1) DSDCp(2,3,3,2) DSDCp(2,3,3,3); 0091 DSDCp(3,1,1,1) DSDCp(3,1,1,2) DSDCp(3,1,1,3) 0092 DSDCp(3,1,2,1) DSDCp(3,1,2,2) DSDCp(3,1,2,3) 0093 DSDCp(3,1,3,1) DSDCp(3,1,3,2) DSDCp(3,1,3,3); 0094 DSDCp(3,2,1,1) DSDCp(3,2,1,2) DSDCp(3,2,1,3) 0095 DSDCp(3,2,2,1) DSDCp(3,2,2,2) DSDCp(3,2,2,3) 0096 DSDCp(3,2,3,1) DSDCp(3,2,3,2) DSDCp(3,2,3,3); 0097 DSDCp(3,3,1,1) DSDCp(3,3,1,2) DSDCp(3,3,1,3) 0098 DSDCp(3,3,2,1) DSDCp(3,3,2,2) DSDCp(3,3,2,3) 0099 DSDCp(3,3,3,1) DSDCp(3,3,3,2) DSDCp(3,3,3,3)]; 0100 lc0101 lc fina:ly determining DSCKappa = the derivative of the 67
PAGE 76
0102 %second PielaKirchhoff Kappa is zero, 0103 0104 0105 diff = 1.0e06*Kappaold; SSdiff = stressfunc(Cold, Cpold, Kappaold + diff); DSDKappa = ( SSdiff SS )/diff; 68
PAGE 77
APPENDIX B. Yield Function Plotting 0001 % this is a function called yieldfunc that is calculating .. 0002 % the yield function, Pick which formula that be used .. 0003% to find yield f; 0004 function yieldf = yieldfunc(C, Cp, Kappa); 0005 global STRESSJNO YIELDF_NO 0006 global NSTEP 0007 temS = stressfunc(C, Cp, Kappa); 0008 if YIELDFJNO == 1 0009 % the f1rst yield funcion is from a paper called (A 0010 % FRAMWORK FOR FINITE STRAIN ELASTOPLASTICITY BASED ON .. 0011 % MAXIMIUM PLASTIC DISSIPATION AND THE MULTILICATIVE .. 0012 % PART I. CONTIIHJUI"i BY J. C ... 0013 % SIMO ,26 NOVEMBER 1986,, Page 214, Equation .. (2.14b) .. 0014 % which has the form 0015 % yield f(S,C) = sgrt(S(i,J)*SCk,l)*CCI,k)*C(J,l) 0016 % 1/3*(S(k,l)*C(k,l))2)R 0017 yieldf = sqrt(trace(temS*C*temS*C) (1/3)_ . _. 0018 *(trace(temS*C))2) Kappa; 0019 % yield = sqrt(trace(temS*Cp*tem8*Cp) 0020 % (1/3;*rtrace(temS*Cp))2)Kappa; 0021 else 69
PAGE 78
0022 yieldfunc = 'No yield function has been chosen' 0023 pause 0024 end 70
PAGE 79
REFERENCES [1) MENDELSON, ALEXANDER .. Plasticity; Theory and Application. New York: Macmillan, 1968. [2) DIXIT, PRAKASH MAHADEO, AND UDAY S. DIXIT .. Modeling of Metal Forming and Machining Processes by Finite Element and Soft Computing Methods. London: Springer, 2008. [3) KHAN, AKHTAR S., AND SuJIAN HUANG .. Continuum Theory of Plasticit New York: Wiley, 1995. [4] ] JAUNZEMIS, WALTER. Continuum Mechanics. New York: Macmillan, 1967. [5] CHADWICK, PETER Continuum Mechanics: Concise Theory and Prob lems. Mineola, NY: Dover Publications, 1999. [6] SPENCER, A. J. M .. Continuum Mechanics. NMineola, NY: Dover Pub lications, 2004. [7) NAGHDI, P.M., AND J. A. TRAPP. On Finite ElasticPlastic Deformation of Metals, Journal of Applied Mechanics 41.1, 1974, 254. [8) ) TRUESDELL, C., AND W. NOLL. The Nonlinear Field Theories of Me chanics Berlin: SpringerVerlag, 1992. [9) SIMO, J .c. "A Framework for Finite Strain Elastoplasticity Based on Maximum Plastic Dissipation and the Multiplicative Decomposition: Part I. Continuum Formulation.", Computer Methods in Applied Mechanics and Engineering 66.2, 1988, 199219. [10] NAGHDI, P. M., AND J. A. TRAPP. Restrictions On Constitutive Equa tions Of Finitely Deformed ElasticPlastic Materials, The Quarterly Journal of Mechanics and Applied Mathematics 28.1 (1975): 2546. [11] PHILLIPS, A. The Effect of Loading Path on the Yield Surface at Elevated Temperatures, International Journal of Solids and Structures 8.4 (1972): 46374. 71
