A constitutive model for polymers: Formulation and Integration algorithm

A constitutive model for polymers: Formulation and Integration algorithm Negar Bahramsari Mestrado em Engenharia Matemática Departamento de Matemática 2014 Orientador Francisco Manuel Andrade Pires, Professor Auxiliar, FEUP Coorientador Sílvio Marques de Almeida Gama, Professor Associado, FCUP Todas as correções determinadas pelo júri, e só essas, foram efetuadas. O Presidente do Júri, Porto, ______/______/_________ Dedicated to my husband ... Acknowledgements I would like to gratefully acknowledge Prof. Pires for accepting to be my thesis supervisor and letting me integrate is his research group. I would also like to thank Prof. Gama for being my thesis co-supervisor and his support during accomplishment of this work. Eng. Mirkhalaf, a PhD student from research group of Prof. Pires, helped me to get acquainted with concepts of continuum mechanics and Finite Element Method and also constitutive modelling of polymers. I am grateful for his support and his comments while I was involved in my thesis. Being far from my country has been difficult and my friends company has been helping me to keep moving. I thank them all although I can not name all of them. At last, not the least, I am going to deeply thank my husband and my family for their undeniable support without which my studies could be much harder to finish. Abstract The objective of this thesis is to develop an elasto-viscoplastic constitutive model in order to characterize the deformation behaviour of polymeric materials. More importantly, it is intended to incorporate several numerical techniques to solve the equilibrium equation obtained from kinematic and material constitutive relations. For this purpose, the thesis starts with a global overview of the state-of-art on the present topic in Chapter 1. In Chapter 2, the main concepts of the Continuum Mechanics Theory and the Finite Element Method are briefly reviewed. This chapter is followed by development of the BPA based constitutive model in Chapter 3. In Chapter 4, the most important contribution of this work is introduced, which is derivation of the integration algorithm and finite element implementation of the model. In order to assess the capability of the model to predict the typical deformation behaviour of polymeric materials and the efficiency of the integration algorithm and the code, some numerical examples are provide in Chapter 5. The document finishes with the presentation of the main conclusions of this work in Chapter 6and, in addition, suggestions for future research are made. Contents Nomenclature vii 1 Introduction 1 1.1 Generalintroduction ............................ 1 1.2 Modeling of the mechanical behavior of polymers . . . . . . . . . . . . 2 1.3 Scopeandoutline ............................. 2 2 Continuum mechanics and finite element method 4 2.1 ContinuumMechanics ........................... 4 2.1.1 Kinematics of deformation . . . . . . . . . . . . . . . . . . . . . 4 2.1.1.1 Motion .......................... 4 2.1.1.2 The deformation gradient . . . . . . . . . . . . . . . . 5 2.1.1.3 Strain measures . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1.4 The velocity gradient . . . . . . . . . . . . . . . . . . . 7 2.1.1.5 Stress measures . . . . . . . . . . . . . . . . . . . . . . 7 2.1.2 Fundamentallaws.......................... 8 2.1.2.1 Conservation of mass . . . . . . . . . . . . . . . . . . . 9 2.1.2.2 Momentum balance . . . . . . . . . . . . . . . . . . . 9 2.1.2.3 The first and second principles of thermodynamics . . 9 2.1.2.4 The Clausius-Duhem inequality . . . . . . . . . . . . . 10 2.1.3 The quasi-static initial boundary value problem . . . . . . . . . 11 2.1.3.1 The material version . . . . . . . . . . . . . . . . . . . 12 2.2 Displacement-based finite elements . . . . . . . . . . . . . . . . . . . . 12 2.2.1 Spatial discretisation . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.2 Temporal discretisation. The non-linear incremental finite elementprocedure........................... 14 vCONTENTS 3 BPA based model 19 3.1 Introduction................................. 19 3.2 Formulation................................. 19 3.2.1 Multiplicative kinematics . . . . . . . . . . . . . . . . . . . . . . 20 3.2.2 Strainmeasure ........................... 22 3.2.3 Flowpotential............................ 23 3.2.4 Stress decomposition . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3 Conclusions ................................. 24 4 The integration algorithm of the model 25 4.1 Introduction................................. 25 4.2 Stateupdate................................. 26 4.2.1 Returnmapping........................... 26 4.3 Consistent tangent operator . . . . . . . . . . . . . . . . . . . . . . . . 35 5 Numerical examples 42 5.1 Compression of a cube (plane strain compression) . . . . . . . . . . . . 42 5.2 Cylinder upsetting (axisymmetric compression) . . . . . . . . . . . . . 44 5.3 Notched bar compression . . . . . . . . . . . . . . . . . . . . . . . . . . 45 6 Conclusions and suggestions for future research 48 References 52 List of Figures 3.1 Elastic-Plastic decomposition of the total deformation gradient-Reproduced from De Souza Neto et al. (2008) ..................... 20 3.2 Polar decomposition of the deformation gradient-Reproduced from De Souza Neto et al. (2008) .......................... 21 5.1 The geometry and mesh of the plane strain compression test. . . . . . . 43 5.2 The stress-strain curve for the compression of the cube. . . . . . . . . . 43 5.3 The stress-strain curve for the compression of the cylinder. . . . . . . . 44 5.4 The geometry and mesh of the compression test on a notched bar . . . 45 5.5 Contour plot of effective stress at ten percent of the deformation. . . . 46 5.6 Contour plot of effective stress at twenty percent of the deformation. . 46 5.7 Contour plot of effective stress at thirty percent of the deformation. . . 46 5.8 Contour plot of effective stress at half of the deformation. . . . . . . . . 47 5.9 Contour plot of effective stress at the end of the deformation. . . . . . . 47 List of Tables 5.1 Material properties for PMMA required for the model . . . . . . . . . . 42 5.2 Local convergence table for cube compression simulation . . . . . . . . 44 5.3 Local convergence table for cylinder compression simulation . . . . . . 45 72.1 Continuum Mechanics tensors, called Eulerian strain tensors, are defined by using the left stretch tensor, V. ε(m)=(1 m(Vm−I), m 6= 0 ln[V], m = 0 (2.14) In this work, Eulerian strain tensor is used and the logarithmic strain tensor (m= 0) is considered. It is also worth noting that when the deformation is just a rotation (i.e., F=R), Eulerian strain tensor and Lagrangian strain tensor are null. U=I=⇒E(m)= 0,(2.15) V=I=⇒ε(m)= 0.(2.16) 2.1.1.4 The velocity gradient Velocity gradient is defined by: L=Oxυ,(2.17) or, alternatively by L=˙ FF−1.(2.18) The velocity gradient is additively composed of symmetric and skew part: L=D+W,(2.19) where, D=sym(L) = 1 2L+LT,W=skew(L) = 1 2L−LT.(2.20) The symmetric part of the deformation gradient, D, is called rate of deformation (stretching) tensor and the skew part, W, is called spin tensor. Stretching tensor is associated with straining whereas spin tensor is associated with rigid velocities. 2.1.1.5 Stress measures The Cauchy stress, σ, the Kirchhoff stress, τ, and the first Piola-Kirchhoff stress, P, are three important stress measures. Which, Cauchy stress tensor or true stress is defined as: t=σn,(2.21) where tis the surface traction and nis its associated normal vector. The Cauchy stress tensor is additively composed of hydrostatic term and deviatoric part σ=s+pI,(2.22) 82.1 Continuum Mechanics where sis the deviatoric stress and pis hydrostatic pressure: s=dev (σ) = Id:σ,(2.23) where Idis the deviatoric fourth order identity tensor: Id=Is−1 3(I⊗I),(2.24) where Isis the symmetric fourth order identity and Iis the second order identity. Is and Iin component form can be shown as: Iijkl =1 2δikδjl +δilδjk,(2.25) Iij =δij,(2.26) where, δis the Kronecker delta. The hydrostatic pressure is given by: p=1 3tr (σ).(2.27) The second stress measure is the Kirchhoff stress tensor denoted by τdefined as τ=Jσ.(2.28) Similarly to the Cauchy stress, the Kirchhoff stress tensor can also be split into two parts, i.e., τ=τd+τhI,(2.29) where τd=dev[τ] and τh=1 3tr[τ] are, respectively, the deviatoric and hydrostatic (or spherical) parts. The last stress measure mentioned here is first Piola-Kirchhoff stress tensor, denoted by P, also known as nominal stress. The Piola-Kichhoff stress is defined by: P=JσF−T=τF−T.(2.30) 2.1.2 Fundamental laws The quantities and definitions mentioned above are essential mathematical relations of deformation. The fundamental laws of continuum mechanics governing the physical phenomena of deformation provide relations between defined and expressed quantities. 92.1 Continuum Mechanics 2.1.2.1 Conservation of mass The conservation of mass postulate implies: ˙ρ+divx(ρ˙ u)=0,(2.31) where ρis the density at the deformed configuration. 2.1.2.2 Momentum balance The momentum balance, also referred to as strong form of the equilibrium equation, of any given body can be expressed as divxσ+b=ρ¨ u,(2.32) where bdenotes the body force vector in the deformed configuration. The equilibrium equation (2.32) needs to fulfill the following boundary condition: t=σn,(2.33) where tis a tension vector applied on the boundary of the body. 2.1.2.3 The first and second principles of thermodynamics The first principle of thermodynamics postulates that the energy must be conserved. This can be mathematically expressed as ρ˙e=σ:D+ρr −divxq,(2.34) where e,rand qare, respectively, the specific internal energy, the density of heat production and the heat flux. Throughout this work, only processes with constant temperature will be considered. In this case, the first principle reduces to ρ˙e=σ:D.(2.35) The equation above states that the rate of internal energy per unit deformed volume must equal the stress power, σ:D, per unit deformed volume. Making use of the relation below, ¯ρ=Jρ, (2.36) 10 2.1 Continuum Mechanics where ¯ρdenotes the reference density, it is possible to rewrite equation (2.35) as ¯ρ˙e=τ:D.(2.37) The second principle of thermodynamics is associated with the so-called irreversibility of entropy production, expressed by the following inequality: ρT ˙s+divxq−ρr ≥0,(2.38) where sdenotes the entropy and Tis the temperature. Similar to the case of the first principle of thermodynamics, if only isothermal processes are considered, Equation (2.38) is then given by ρT ˙s≥0.(2.39) Pre-multiplying the equation (2.39) by J, we have ¯ρT ˙s≥0.(2.40) 2.1.2.4 The Clausius-Duhem inequality Firstly, we introduce the Helmholtz free energy, ψ, defined by ψ=e−Ts (2.41) Re-arranging the equation (2.41) and differentiating with respect to time, we have T˙s= ˙e−˙ ψ. (2.42) We remark that the temperature has been assumed constant, thus, ˙ T= 0 Using the equation (2.35), we conclude that ¯ρT ˙s=τ:D−¯ρ˙ ψ, (2.43) with τ:D−¯ρ˙ ψ≥0.(2.44) The fundamental inequality above is called Clausius-Duhem inequality. 11 2.1 Continuum Mechanics 2.1.3 The quasi-static initial boundary value problem The fundamental laws presented in the last section allow us to define an initial boundary value problem (IBVP) associated with the description of a given deformation process. The solution of the IBVP delivers the prediction of how a given solid will mechanically behave when subjected to certain boundary conditions. Within the scope of this work, only quasi-static problems will be addressed, hence, all inertia effects will be neglected. Thus, the equilibrium equation (2.32), in the strong form, is re-written to be given by divxσ= 0.(2.45) Multiplying the equation (2.45) by a virtual displacement, η, and integrating over the volume, we have Zϕ(Ω) (divxσ)TηdV = 0.(2.46) After some straightforward operations, (2.46) becomes Zϕ(Ω) [σ:Oxη−(divxσ·η)]dV = 0.(2.47) Making use of the divergence theorem, we have Zϕ(Ω) σ:OxηdV −Zϕ(∂Ω) (σ·n)TηdA = 0.(2.48) Finally, substituting Equation (2.33) into Equation (2.48) leads to Zϕ(Ω) σ:OxηdV −Zϕ(∂Ω) t·ηdA = 0.(2.49) Equation (2.49) is called weak form of the equilibrium equation. The use of the weak form can significantly facilitate the use of efficient numerical methods for the solution of the structural IBVP. With the definition of weak equilibrium at hand, we can define the quasi-static IBVP, in the spatial description, as follows. Problem 2.1 Given a prescribed deformation gradient history F(t) = I+Opu(p, t) (2.50) and the Cauchy stress, at each point of the body expressed as σ(t) = σ(F(t, α(t)) (2.51) 12 2.2 Displacement-based finite elements obtained from the solution of the constitutive initial boundary value problem where αis the set of internal variables associated with the material, find a kinematically admissible displacement function, u∈Ksuch that the equation Zϕ(Ω) σ(t) : OxηdV −Zϕ(∂Ω) t(t)·ηdA = 0,(2.52) is satisfied for all t∈[t0,tn] and for all η∈νt. The set of kinematically admissible displacements, K, and the space of virtual displacements at time t, νt, are respectively given by K={u: Ω →U|u(p,t) = ¯ u(p,t),t∈[t0,tn],p∈∂ϕu},(2.53) νt={η:ϕ→U|η= 0 ∈∂ϕu(t)}.(2.54) Unfortunately, analytical solutions for the problem defined above exist only for a restricted set of special cases. For accurate predictions of the mechanical behavior of solids in the general case, the use of numerical methods is therefore indispensable. 2.1.3.1 The material version For reference, the material version of the weak form of the equilibrium equation is also herein provided, which reads Zϕ P:∇pηdV −Z∂ϕ ¯ t·ηdA= 0,(2.55) where ¯ tis the surface traction per unit reference area. 2.2 Displacement-based finite elements In this section, the general concepts of the finite element method formulated with a displacement-based approach will be briefly addressed. The finite element method has been chosen in this work as the base numerical tool mainly due to its versatility and its high effectiveness when adopted for the simulation of deformation processes. 2.2.1 Spatial discretisation As stressed out above, the solution of Problem 2.1 often requires the use of some sort of numerical strategy. Within a typical finite element framework, the field variables 13 2.2 Displacement-based finite elements are discretised through the so-called interpolation or shape functions. In the case of displacement-based finite elements, the interpolated field variable are the displacements. Within a given element e, the interpolation is assumed to be u(x)≡ nnode X i=1 N(e) i(x)ui,(2.56) where N(e) i(x) is the shape function associated with node i(evaluated at x) and nnode is the number of nodes of the element. In similar manner, a global interpolation function can also be set, that is, u(x)≡ npoin X i=1 Ng i(x)ui,(2.57) where npoin is the total number of nodes of the finite element mesh and Ng i(x) is the global interpolation matrix which can be represented by Ng(x) = [diag[Ng 1(x)] diag[Ng 2(x)] ··· diag[Ng npoin (x)]] (2.58) where diag[Ng i] denotes a ndim ×ndim diagonal matrix defined diag[Ng i] =      Ng i0··· 0 0Ng i··· 0 . . .. . ..... . . 0 0 ··· Ng i     .(2.59) At this point, it is also convenient to define the global vector of nodal displacements, given by u= [u1 1,···undim 1,······ u1 npoin ,···unpoin ndim ]T.(2.60) With the above matrix notation at hand, equation (2.55) can be rephrased to be given by u(x) = Ng(x)u,(2.61) where the equation above represents the interpolation of the displacement field by means of discrete functions. Analogously, we can write the field of virtual displacements to be given by η(x) = Ng(x)η.(2.62) We also define the global discrete symmetric gradient matrix, Bg, which in the case of plane stress and plane strain problems assumes the form Bg=  Ng 1,10Ng 2,1··· Ng npoin,10 0Ng 1,20··· 0Ng npoin,2 Ng 1,2Ng 1,1Ng 2,2··· Ng npoin,2Ng npoin,1  ,(2.63) 14 2.2 Displacement-based finite elements where use of the following notation has been made: (·)i,j =∂(·)i ∂xj .(2.64) For completeness, the global discrete full gradient operator,Gg, is also provided herein, whose format in plane stress and plane strain analyses is given by Gg=    Ng 1,10Ng 2,10··· Ng npoin,10 0Ng 1,10Ng 2,1··· 0Ng npoin,1 Ng 1,20Ng 2,20··· Ng npoin,20 0Ng 1,20Ng 2,2··· 0Ng npoin,2    .(2.65) 2.2.2 Temporal discretisation. The non-linear incremental finite element procedure In practical engineering applications, it is often required the modeling of materials that are dependent of the deformation history. Such materials are called path dependent and, regardless whether they take strain rate effects into account or not, a suitable temporal discretisation needs to be performed. Within the context of most of this work, a pseudo-time discretisation between the time increments [tn,tn+1] will be considered for which a fully implicit scheme is adopted. For a generic path-dependent material model, an incremental constitutive function, ˆ σ, is assumed to exist, i.e., σn+1 =ˆ σ(Fn+1,αn).(2.66) In practice, the function ˆ σis associated with an integration algorithm that computes the material behavior for a given deformation gradient Fn+1 and set of internal variables, αn, which remain constant within [tn,tn+1]. Accordingly, a similar incremental function for the set of internal variables is also defined: αn+1 =ˆ α(Fn+1,αn).(2.67) Up to this point, it suffices to leave the incremental constitutive functions ˆ σand ˆ α unspecified for the sake of generality. The non-linear incremental finite element equation The definition of the incremental constitutive function of Section 2.2.2, combined with the spatial discretisation presented in Section 2.2.1, allows the definition of the 15 2.2 Displacement-based finite elements incremental finite element equilibrium equation, obtained after some straight forward substitutions and rearrangements from Equation (2.49): r(un+1)≡fint(un+1)−fext = 0,(2.68) where fint and fext are, respectively, the internal and external force vectors, defined for a given element eas fint =Zϕn+1(ϕ(e)) BTˆ σ(Fn+1,α)dV, (2.69) fext =Zϕn+1(∂ϕ(e)) NTtn+1dA, (2.70) Within the scope of this work, Equation (2.64) will be generally non-linear and therefore requires an adequate (numerical) method for its solution. Considering an incremental scheme, in which a given fraction of the external prescribed load is applied at each increment, Equation (2.64) is solved as summarized in Box 1. Remark 2.1. In practice, the external force is computed by the expression fext n+1 =λn+1¯ fext,(2.71) where λn+1 is the prescribed load factor at time tn+1 and ¯ fext is computed only once at the first iteration of the incremental procedure through the expression ¯ fext =Zϕn+1(∂Ω(e)) NT¯ tdA (2.72) where ¯ tis a prescribed field, which remains constant through the incremental procedure. Numerical integration of fint and fext One important aspect of the finite element implementation is the substitution of the exact integrals of Equations (2.69) and (2.70) by some numerical approximations. In this work, both the internal and external force vectors will be integrated using standard Gaussian quadratures. For instance, in the case of the internal force vector, fint is approximated by the following expression: fint (e)=Zϕn+1(Ω(e)) BTˆ σdV ≈ ngp X i=1 wiJiBT iˆ σi(2.73) 16 2.2 Displacement-based finite elements Box 1: The incremental non-linear finite element scheme - implicit solution. 1. Assemble the global external force vector,¯ fext 2. Initialize increment counter i= 1 3. Set load factor ,λi 4. Solve the non-linear equilibrium equation r(un+1) = fint(un+1)−λi¯ fext = 0 5. Update increment counter i=i+ 1 6. Check if prescribed number of increments has been achieved where wiand Jiare, respectively, the Gaussian weight and the Jacobian at the ith integration point. A similar procedure is carried out for the external force vector. Newton-Raphson method As stressed out before, the equilibrium equation (2.68) is generally non-linear and demands an appropriate solution method. We will adopt herein the Newton-Raphson method, which is particularly attractive due to its quadratic rates of convergence. Following standard procedures of the method and particularising for the case of the finite element framework presented in this chapter, the displacements are updated according to: uk+1 n+1 =uk n+1 −∂r ∂un+1 uk n+1 −1 ruk n+1(2.74) Equation (2.74) can be more conveniently written as KTδu=−r(uk n+1),(2.75) where u=uk+1 n+1 −uk n+1,(2.76) and KTis called global tangent stiffness matrix, given by KT=∂r ∂un+1 uk n+1 .(2.77) The correct derivation of the tangent stiffness is crucial to guarantee the quadratic rates of convergence of the Newton-Raphson method. In the case of finite strains 23 3.2 Formulation 3.2.3 Flow potential The dissipation potential and flow rule of the model are defined using the Argon theory. The plastic stretching tensor at the current configuration is given by: dp= ˙γp∂Ψ ∂σ= ˙γpN,(3.13) where, Ψ is a dissipation potential, Nis a flow vector and dpis spatial plastic stretching tensor: dp≡ReDpReT.(3.14) Different terms in relation (3.14) could be explained and interpreted as follows: •Dpis plastic stretching tensor in the relaxed configuration. •dpis plastic stretching tensor in the current configuration (rotated from intermediate configuration using elastic rotation). As stated above, the flow vector Nis the derivative of dissipation potential in order to stress. Hence, the dissipation potential should be defined in order to obtain the flow vector. The dissipation potential is assumed to be: Ψ = r1 2s:s.(3.15) As a result, the flow vector is obtained: N=∂Ψ ∂σ=r1 2 s ||s||,(3.16) where, ||s|| is the norm of sdefined by: ||s|| =√s:s.(3.17) Combining relations (3.13) and (3.16), the multi-dimensional plastic flow rule of the model is obtained as: dp= ˙γpr1 2 s ||s||.(3.18) Plastic shear strain rate, ˙γp, is given by: ˙γp= ˙γ0exp −∆G∗ kT ,(3.19) where, ˙γ0is pre-exponential shear shear strain rate factor (one of the material properties); kis Boltzmanns constant; Tis the absolute temperature and ∆G∗is the free energy change. The free energy change is given by: ∆G∗=3πµω2α3 16(1 −ν) 1− τ 0.077µ 1−ν!(5/6) ,(3.20) 24 3.3 Conclusions where, µis shear modulus; νis poissons ratio; ωis the net angle of rotation of the molecular segment between the initial configuration and the activated configuration and αis the mean molecular radius. Relation (3.19) can be rewritten as follows: ˙γp= ˙γ0exp −Aes T1−τ es(5/6),(3.21) where, τis the shear stress. The parameter esis given by: es=s+αp, (3.22) where, pis hydrostatic pressure. The evolution of sis characterized by: ˙s=h1−s sss ˙γp,(3.23) s=s0+αp (3.24) Aand s0are defined by relations (3.25) and (3.26), respectively. A≡39πω2α3 16K,(3.25) s0≡0.077µ 1−ν.(3.26) 3.2.4 Stress decomposition The total stress, σtotal, of the model, is additively composed of driving stress and hardening stress: σtotal =σdriving +σhardening.(3.27) In this study, the hardening stress is not taken into account. Since the main objective of this thesis is to use numerical methods, Finite Element Method and NewtonRaphson method, to solve the nonlinear system of equations of equilibrium, the total stress is assumed to be the driving component of the stress. 3.3 Conclusions In this chapter a finite strain elasto-viscoplastic constitutive model was introduced. The kinematic and material contribution of the constitutive model were introduced. The description of the flow potential, which is in fact the core of the model, was introduced by Boyce et al. (1988) and the model is known as BPA model. Chapter 4 The integration algorithm of the model 4.1 Introduction This chapter deals with the integration algorithm of the model introduced in chapter (3). Stress integration algorithm and consistent tangent operator are explained. The numerical implementation of of stress integration algorithm and tangent operator is demonstrated. (Ortiz et al. (1983); Zienkiewicz and Taylor (1991); Owen and Hinton (1980); De Souza Neto et al. (2008); (Simo and Hughes ,1998); Simo et al. (1998)) studied numerical solution of constitutive models within finite element approach which basically is solution of a set of evolutionary equations in an iterative fashion. The solution pursued in this study has a strain driven structure and is performed at each Gauss point of the finite element mesh. Stress and internal variables are updated according to the the level of current strain and the values of the internal variables in the previous increment. In order to have the solution of the global boundary value problem, using NewtonRaphson method, efficiently converged, the consistent linearization of the time discritized constitutive equations is very important. Using operator split algorithms are nowadays standard procedure for numerical integration of elasto-plastic and elasto-viscoplastic constitutive equations (Criesfield (1997); Simo and Hughes (1998)). In order to implement the constitutive model into a finite element code, it is required to derive the state update equations and also consistent tangent operator. The state update procedure is derived using operator split algorithm and as a result an elastic predictor/return mapping algorithm is obtained for the state update procedure. 26 4.2 State update It is worth emphasizing that state update and consistent tangent are derived on the small strain format of the constitutive equations. Then they are extended to finite strain as explained (De Souza Neto et al.,2008). 4.2 State update In this section, the constitutive equations of the model introduced in chapter (3) are time discretized and operator split algorithm is applied to the equations in order to derive the elastic predictor/return mapping algorithm of the state update procedure. Initially, the total strain at the beginning of state update is assumed elastic (this stage is called elastic predictor) and then using the return mapping equations (derived form the flow potential and other constitutive equations), the inelastic contribution of the total strain would be determined through an iterative procedure (this stage is called plastic corrector or return mapping). 4.2.1 Return mapping In order to update the deformation gradient and obtain the deformation gradient at time step tn+1, we use: Fn+1 =F∆Fn.(4.1) where, Fn+1 is the deformation gradient at tn+1,Fnis the deformation gradient at tn and F∆is the incremental deformation gradient. Considering incremental displacement known, the incremental deformation gradient is obtained by: F∆=I+5n(∆u),(4.2) where, 5n(∆u) is the gradient of the incremental displacement. We need to obtain the elastic trial state so that we will then apply return mapping procedure in order to determine which part of the strain has been elastic and which part has been plastic strain. As we know the total elastic strain at tnwhich is εe n, the elastic left Cauchy-Green deformation tensor is obtained by: Be n= exp[2 εe n].(4.3) The elastic trial left Cauchy-Green deformation tensor at tn+1 is given by: Be trial n+1 =F∆Be n(F∆)T.(4.4) 27 4.2 State update Having elastic trial left Cauchy-Green deformation tensor computed, the driving parameter of the study which is the elastic trial strain at tn+1 i.e. εe trial n+1 is given by: εe trial n+1 = ln [Ve trial n+1 ] = 1 2ln [Be trial n+1 ].(4.5) It is worth noting that so far, everything is independent of the material models i.e. everything is done at the kinematic level and for whatever constitutive model could be used. At this stage, when we have the elastic trial strain determined, the stresses and state variables should be updated using the constitutive relations. In order to derive the integration algorithm of the model, the small strain counterpart of relation (3.18) is used: ˙ εp= ˙γpr1 2 s ||s||.(4.6) where, sis the deviatoric part of the Kirchhof stress. s=Id:τ,(4.7) and ||s|| is the norm of the deviatoric stress. ||s|| =√s:s.(4.8) Integration of relation (4.6) over the time step [tn, tn+1] gives: εp n+1 −εp n= ˙γp n+1r1 2 ∆t ||sn+1||sn+1 (4.9) The time discretized form of relation (3.21) is given by: ˙γp n+1 = ˙γ0exp "−Aesn+1 T 1−τ esn+1 (5/6)!#,(4.10) where, esn+1 =sn+1 +αpn+1 (4.11) Integrating relation (3.23) over time step [tn, tn+1] gives relation (4.12): sn+1 −sn−h1−sn+1 sss ˙γp n+1∆t= 0 (4.12) 28 4.2 State update From relation (4.9), the following system of four non-linear equations are obtained:                            ∆εP n+1(1) = ˙γp n+1q1 2 ∆t ||sn+1|| sn+1(1) ∆εP n+1(2) = ˙γp n+1q1 2 ∆t ||sn+1|| sn+1(2) ∆εP n+1(3) = ˙γp n+1q1 2 ∆t ||sn+1|| sn+1(3) ∆εP n+1(4) = ˙γp n+1q1 2 ∆t ||sn+1|| sn+1(4) (4.13) Considering relations (4.10)-(4.12), and the system of equations (4.13), the following system of six algebraic equations are obtained, in the two dimensional space, for the finite element implementation of the model.                                                  R1:= ∆εP n+1(1) −˙γp n+1q1 2 ∆t ||sn+1|| sn+1(1) = 0 R2:= ∆εP n+1(2) −˙γp n+1q1 2 ∆t ||sn+1|| sn+1(2) = 0 R3:= ∆εP n+1(3) −˙γp n+1q1 2 ∆t ||sn+1|| sn+1(3) = 0 R4:= ∆εP n+1(4) −˙γp n+1q1 2 ∆t ||sn+1|| sn+1(4) = 0 R5:= ˙γp n+1 −˙γ0exp −Aesn+1 T1−τ esn+1 (5/6)= 0 R6:= sn+1 −sn−h1−sn+1 sss ˙γp n+1∆t= 0 (4.14) The unknowns of the system of equations (4.14) are ∆εP n+1(1), ∆εP n+1(2), ∆εP n+1(3), ∆εP n+1(4), ˙γp n+1 and sn+1. The system of equation will be solved suing the well-known iterative Newton-Raphson method. The unknowns are called u1,u2,u3,u4,u5and u6, respectively. In order to solve the system of equations using the Newton-Raphson method in an iterative fashion, the derivatives of all equations in order to all unknowns should be calculated:      ∂R1 ∂u1 ∂R1 ∂u2 ∂R1 ∂u3··· ∂R1 ∂u6 ∂R2 ∂u1 ∂R2 ∂u2 ∂R2 ∂u3··· ∂R2 ∂u6 . . .. . ..... . . ∂R6 ∂u1 ∂R6 ∂u2 ∂R6 ∂u3··· ∂R6 ∂u6              δu1 δu2 δu3 δu4 δu5 δu6         k =−         R1 R2 R3 R4 R5 R6         k−1 .(4.15) 29 4.2 State update Where, superscripts k−1 and kstand for two consecutive Newton-Raphson iterations. The required derivatives in relation (4.15) are provided below. Derivatives of the first residual equation Derivative of the first residual equation in order to first unknown: ∂R1 ∂u1 = 1 −˙γp n+1r1 2∆t C1,(4.16) where, C1=||sn+1||C2−sn+1(1)C3 (||sn+1||)2,(4.17) where, C2and C3are given by: C2=−2G , C3=−2Gsn+1(1) ||sn+1|| .(4.18) Derivative of the first residual equation in order to second unknown: ∂R1 ∂u2 =−˙γp n+1r1 2∆t C4,(4.19) where, we have: C4=−sn+1(1)C5 (||sn+1||)2, C5=−2Gsn+1(2) ||sn+1||.(4.20) Derivative of the first residual equation in order to third unknown: ∂R1 ∂u3 =−˙γp n+1r1 2∆t C6,(4.21) where, C6=−sn+1(1)C7 (||sn+1||)2, C7=−4Gsn+1(3) ||sn+1||.(4.22) Derivative of the first residual equation in order to fourth unknown: ∂R1 ∂u4 =−˙γp n+1r1 2∆t C8,(4.23) where, C9=−sn+1(1)C9 (||sn+1||)2, C9=−2Gsn+1(4) ||sn+1||.(4.24) Derivative of the first residual equation in order to fifth unknown: ∂R1 ∂u5 =−r1 2∆tsn+1(1) ||sn+1||.(4.25) 30 4.2 State update Derivative of the first residual equation in order to sixth unknown: ∂R1 ∂u6 = 0.(4.26) Derivatives of the second residual equation Derivative of the second residual equation in order to first unknown: ∂R2 ∂u1 =−˙γp n+1r1 2∆t C10,(4.27) where, C10 =−sn+1(2)C3 (||sn+1||)2.(4.28) Derivative of the second residual equation in order to second unknown: ∂R2 ∂u2 = 1 −˙γp n+1r1 2∆t C11,(4.29) where, C11 =||sn+1||C12 −sn+1(2)C5 (||sn+1||)2,(4.30) where, C12 and is given by: C12 =−2G. (4.31) Derivative of the second residual equation in order to third unknown: ∂R2 ∂u3 =−˙γp n+1r1 2∆t C13,(4.32) where, C13 =−sn+1(2)C7 (||sn+1||)2.(4.33) Derivative of the second residual equation in order to fourth unknown: ∂R2 ∂u4 =−˙γp n+1r1 2∆t C14,(4.34) where, C14 =−sn+1(2)C9 (||sn+1||)2.(4.35) Derivative of the second residual equation in order to fifth unknown: ∂R2 ∂u5 =−r1 2∆tsn+1(2) ||sn+1||.(4.36) 31 4.2 State update Derivative of the second residual equation in order to sixth unknown: ∂R2 ∂u6 = 0.(4.37) Derivatives of the third residual equation Derivative of the third residual equation in order to first unknown: ∂R3 ∂u1 =−˙γp n+1r1 2∆t C15,(4.38) where, C15 =−sn+1(3)C3 (||sn+1||)2.(4.39) Derivative of the third residual equation in order to second unknown: ∂R3 ∂u2 =−˙γp n+1r1 2∆t C16,(4.40) where, C16 =−sn+1(3)C5 (||sn+1||)2.(4.41) Derivative of the third residual equation in order to third unknown: ∂R3 ∂u3 = 1 −˙γp n+1r1 2∆t C17,(4.42) where, C17 =||sn+1||C18 −sn+1(3)C7 (||sn+1||)2,(4.43) where, C18 and is given by: C18 =−2G. (4.44) Derivative of the third residual equation in order to fourth unknown: ∂R3 ∂u4 =−˙γp n+1r1 2∆t−sn+1(3)C9 (||sn+1||)2.(4.45) Derivative of the third residual equation in order to fifth unknown: ∂R3 ∂u5 =−r1 2∆tsn+1(3) ||sn+1||.(4.46) Derivative of the third residual equation in order to sixth unknown: ∂R3 ∂u6 = 0.(4.47) 32 4.2 State update Derivatives of the fourth residual equation Derivative of the fourth residual equation in order to first unknown: ∂R4 ∂u1 =−˙γp n+1r1 2∆t C19,(4.48) where, C19 =−sn+1(4)C3 (||sn+1||)2.(4.49) Derivative of the fourth residual equation in order to second unknown: ∂R4 ∂u2 =−˙γp n+1r1 2∆t C20,(4.50) where, C20 =−sn+1(4)C5 (||sn+1||)2.(4.51) Derivative of the fourth residual equation in order to third unknown: ∂R4 ∂u3 =−˙γp n+1r1 2∆t C21,(4.52) where, C21 =−sn+1(4)C7 (||sn+1||)2.(4.53) Derivative of the fourth residual equation in order to fourth unknown: ∂R4 ∂u4 = 1 −˙γp n+1r1 2∆t C22,(4.54) where, C22 =||sn+1||C23 −sn+1(4)C9 (||sn+1||)2.(4.55) where, C18 and is given by: C23 =−2G(4.56) Derivative of the third residual equation in order to fifth unknown: ∂R4 ∂u5 =−r1 2∆tsn+1(4) ||sn+1||.(4.57) Derivative of the fourth residual equation in order to sixth unknown: ∂R4 ∂u6 = 0.(4.58) 39 4.3 Consistent tangent operator In order to compute the required unknowns for the tangent operator, D, we need to have the derivatives of the residual equations, R1, R2, R3, R4, R5, R6, R7, in order to unknowns,∆εP n+1(1), ∆εP n+1(2), ∆εP n+1(3), ∆εP n+1(4), ˙γp n+1, sn+1 and pn+1 and also in order to strain εn+1. For the state update stage of the solution the derivatives of the first six residual equations, R1, R2, R3, R4, R5, R6, in order to first six unknowns, ∆εP n+1(1), ∆εP n+1(2), ∆εP n+1(3), ∆εP n+1(4), ˙γp n+1, sn+1, are computed and provided in relations (4.16)-(4.81). We need to Compute: •The derivatives of the the first six residual equations in order to hydrostatic pressure, pn+1: •The derivatives of seventh equation, R7, in order to all unknowns: •The derivatives of all equations in order to strain to strain. The first set of required calculations are the derivatives in order to pressure: ∂R1 ∂pn+1 =∂R2 ∂pn+1 =∂R3 ∂pn+1 =∂R4 ∂pn+1 = 0,(4.109) ∂R5 ∂pn+1 =−˙γ0exp (C24)K1,(4.110) where, K1=∂C24 ∂pn+1 ,(4.111) which is given by: K1=−Aα T 1−τn+1 esn+1 (5/6)!−Aesn+1 T 5 6τn+1 esn+1 (−1/6) K2!,(4.112) where, K2=−ατn+1 esn+1 .(4.113) The other derivatives in order to pressure is given below: ∂R6 ∂pn+1 = 0 ,∂R7 ∂pn+1 = 1.(4.114) The last equation, R7, is only dependent on pressure, pn+1, and consequently the derivatives of the last residual equation in order to all unknowns but pressure, pn+1, are zero. ∂R7 ∂u1 =∂R7 ∂u2 =∂R7 ∂u3 =∂R7 ∂u4 =∂R7 ∂u5 =∂R7 ∂u6 = 0.(4.115) 40 4.3 Consistent tangent operator The last set of required relations to be determined are the derivatives of the all residual equations in order to strain. It should be emphasized that since strain, εn+1, is a second order tensor and the equations are scalar quantities, the derivatives of each residual equation in order to strain result in a second order tensor. Now, we proceed with the derivatives R1, R2, R3, R4in order to strain: ∂R1 ∂εn+1 =? ∂R1 ∂εn+1 = ˙γP n+1r1 2∆t∂ ∂εn+1 (sn+1(1) ||sn+1||),(4.116) where, ∂sn+1(1) ∂εetrial n+1 = 2G       2 3 −1 3 0 −1 3        ,(4.117) and the next one ∂R2 ∂εn+1 =?, ∂R2 ∂εn+1 = ˙γP n+1r1 2∆t∂ ∂εn+1 (sn+1(2) ||sn+1||),(4.118) where, we have: ∂sn+1(2) ∂εetrial n+1 = 2G       −1 3 2 3 0 −1 3        .(4.119) Now, it is the turn of ∂R3 ∂εn+1 =? ∂R3 ∂εn+1 = ˙γP n+1r1 2∆t∂ ∂εn+1 (sn+1(3) ||sn+1||),(4.120) where, ∂sn+1(3) ∂εetrial n+1 = 2G       0 0 1 2 0        .(4.121) The derivative ∂R4 ∂εn+1 is given in relation (4.122): ∂R4 ∂εn+1 = ˙γP n+1r1 2∆t∂ ∂εn+1 (sn+1(4) ||sn+1||),(4.122) where, ∂sn+1(4) ∂εetrial n+1 = 2G       −1 3 −1 3 0 2 3        .(4.123) 41 4.3 Consistent tangent operator The norm of stress deviator is given by: ||sn+1|| =psn+1(1)sn+1(1) + sn+1(2)sn+1(2) + 2sn+1(3)sn+1(3) + sn+1(4)sn+1(4). (4.124) The derivative of the fifth equation in order to strain is given in: ∂R5 ∂εn+1 =−˙γ0exp −Aesn+1 θ(1 −(τn+1 esn+1 ) 5 6)×, ∂ ∂εn+1 [−Aesn+1 θ(1 −(τn+1 esn+1 ) 5 6)].(4.125) Relation (4.125) could be rewritten as: ∂R5 ∂εn+1 =−˙γ0exp −Aesn+1 θ(1 −(τn+1 esn+1 ) 5 6)×, −Aesn+1 θ×∂ pεn+1 [−(τn+1 esn+1 ) 5 6],(4.126) or in a more simplified form: ∂R5 ∂εn+1 =−˙γ0exp −Aesn+1 θ(1 −(τn+1 esn+1 ) 5 6)×, +Aesn+1 θ×∂ ∂εn+1 [(τn+1 esn+1 ) 5 6],(4.127) where τn+1 is an equivalent stress defined by: τn+1 =r1 2sn+1 :sn+1.(4.128) Using relations (4.127) and (4.128), we can have: ∂R5 ∂εn+1 =−5√2 6 ˙γ0AG θ||sn+1||(τn+1 esn+1 )−1 6sn+1 exp{−Aesn+1 θ(1 −(τn+1 esn+1 ) 5 6)}.(4.129) The other derivatives are provided below: ∂R6 ∂εn+1 = 0,(4.130) ∂R7 ∂εn+1 =KI.(4.131) Chapter 5 Numerical examples This chapter presents some numerical examples through which the capability of the model to characterize the typical deformation behaviour of polymers is assessed and also the efficiency of the derived and implemented algorithm is shown. The material under study is polymethylmethacrylate (PMMA). The material properties required for performing the simulations are taken from (Boyce et al. ,1988). The material properties are tabulated in Table (5.1). Different examples including Table 5.1: Material properties for PMMA required for the model E ν s0α˙γ0A h sss R PMMA 2300 0.37 88E+6 0.2 1.13E+11 167E-06 900E+06 77E+06 8.3143 compression on a cube, cylinder compression, notched bar compression and necking of a round bar are selected to be analysed. 5.1 Compression of a cube (plane strain compression) In this section, a compression test on a cube is performed which can be approximated in 2D as plane strain compression. The geometry and mesh of the problem is shown in Figure (5.1). The specimen is spatially discretized with 128 quadrilateral eight noded elements with reduced four integration Gauss points. The stress strain curve for the plane strain compression test is given in Figure (5.2). In Figure (5.2), the typical deformation behaviour of glassy polymers which includes initial linear elastic, yield, post-yield softening could be observed. It should be emphasized that the reason for not-seeing the final hardening regime in the deformation is that the hardening stress 43 5.1 Compression of a cube (plane strain compression) 8 mm 8 mm Figure 5.1: The geometry and mesh of the plane strain compression test. 0,00E+00 1,00E+07 2,00E+07 3,00E+07 4,00E+07 5,00E+07 6,00E+07 7,00E+07 0 0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 True compressive stress True compressive strain Figure 5.2: The stress-strain curve for the compression of the cube. is not considered in this study. In order to assess the performance and numerical efficiency of the derived and implemented algorithm, the convergences are going to to be checked. In Table (5.2), the local convergence of the problem at two different load increments is shown. From Table (5.2), it can be concluded that the state update algorithm is derived 44 5.2 Cylinder upsetting (axisymmetric compression) Table 5.2: Local convergence table for cube compression simulation The value of the residual at local level Iteration number Increment (6) Increment (81) 1 0.383915E-05 1.885185 2 0.295615E-10 0.194465 3 2.044682E-03 4 4.130582E-06 5 8.780845E-09 properly. 5.2 Cylinder upsetting (axisymmetric compression) This section presents results of a compression test on a cylinder. The Geometry is the same as the one shown in Figure (5.1) and the same material properties but under axisymmetric condition is used for the simulation. Figure (5.3), shows stressstrain curve for the cylinder upsetting simulation. In order the check the numerical 0,00E+00 1,00E+07 2,00E+07 3,00E+07 4,00E+07 5,00E+07 6,00E+07 0 0,05 0,1 0,15 0,2 0,25 0,3 0,35 0,4 True compressive stress True compressive strain Figure 5.3: The stress-strain curve for the compression of the cylinder. efficiency, the convergence of the cylinder upsetting simulation are investigated. Table (5.3) presents local convergence (state update) for two different load increments. 45 5.3 Notched bar compression Table 5.3: Local convergence table for cylinder compression simulation The value of the residual at local level Iteration number Increment (6) Increment (35) 1 0.779533 0.532699 2 0.010266E-04 0.880138E-02 3 0.312095E-08 0.059476E-04 4 0.740346E-05 5 0.737990E-07 5.3 Notched bar compression The next example is a notched bar under uniaxial compression. The geometry and mesh of the problem is shown in Figure (5.4). The specimen is spatially discretized Figure 5.4: The geometry and mesh of the compression test on a notched bar with 350 eight noded quadrilateral elements with reduced four integration Gauss points. Figure (5.5) depicts contour plot of the effective stress in whole specimen when ten percent of the deformation is applied. Figure (5.6) shows contour plot of the effective stress in whole specimen at twenty percent of the deformation. Effective stress in the specimen when 30 percent of the deformation is applied is shown in Figure (5.7). When half of the deformation has proceed, the effective stress is given in Figure (5.8). The final contour plot ,depicted in Figure (5.9), shows the effective stress at the end of the deformation. 46 5.3 Notched bar compression 3.87e+004 3.45e+007 X Y Equivalent Stress (0.1) Z 6.89e+007 Figure 5.5: Contour plot of effective stress at ten percent of the deformation. 3.87e+004 3.45e+007 X Y Equivalent Stress (0.2) Z 6.89e+007 Figure 5.6: Contour plot of effective stress at twenty percent of the deformation. 3.87e+004 3.45e+007 X Y Equivalent Stress (0.3) Z 6.89e+007 Figure 5.7: Contour plot of effective stress at thirty percent of the deformation. 47 5.3 Notched bar compression 3.87e+004 3.45e+007 X Y Equivalent Stress (0.5) Z 6.89e+007 Figure 5.8: Contour plot of effective stress at half of the deformation. 3.87e+004 3.45e+007 X Y Equivalent Stress (1) Z 6.89e+007 Figure 5.9: Contour plot of effective stress at the end of the deformation. Chapter 6 Conclusions and suggestions for future research In this thesis, a constitutive model was developed based on BPA model, which is one of the most well-known constitutive models proposed for polymers. Finite Element Method was used to solve the equilibrium equation. The integration algorithm of the model was presented for finite element implementation. The well-known NewtonRaphson method is used at state update (local) level and also equilibrium (global) level to solve the non-linear system of equations. Hence, in this thesis we have used two important numerical methods, namely finite element method and Newton-Raphson method, to deal with equilibrium equations of solid polymers. It is worth emphasizing that, there might be a possibility to improve the integration algorithm of the model by reducing the system of equations to another system of equations with reduced number of equations or even to a single scalar equation. According to the results presented in chapter five of this work, it can be concluded that the BPA based model can predict the typical polymeric materials behavior. In order to have an idea of the capabilities and drawbacks of the model predictions, it is required to have a wide range of experimental results under different deformation modes. Comparing simulations and experimental results will help to realize the accuracy of the model predictions and also the necessity of required changes and modifications to the model. Considering all the remarks provided above, the suggestions for the future research can be summarized as follows: •Checking the possibility of modification of the integration algorithm •Making some comparisons with existing experimental results •Proposing some applicable modifications to the constitutive relations in order to improve model predictions