Numerical modelling of complex geomechanical problems

NUMERICAL MODELLING OF COMPLEX GEOMECHANICAL PROBLEMS by AGUSTÍ P ÉREZ FOGUET UNIVERSITAT POLITÈCNICA DE CATALUNYA ESCOLA TÈCNICA SUPERIOR D’ENGINYERS DE CAMINS, CANALS I PORTS DE BARCELONA DEPARTAMENT DE MATEM ÀTICA APLICADA III Advisors: Antonio Huerta Antonio Rodr íguez-Ferran Barcelona, October 2000 Doctoral Thesis To Paloma and my parents Contents Acknowledgements xv 1 Introduction 1 1.1 Kinematicformulation.............................. 2 1.2 Nonlinearsolvers................................. 3 1.3 Constitutive modelling .............................. 4 2 Arbitrary Lagrangian–Eulerian formulations 7 2.1 Analysisofthevanetestconsideringsizeandtimeeffects .......... 8 2.1.1 Introduction ............................... 8 2.1.2 Mainfeaturesoffieldvanetest ..................... 10 2.1.3 Constitutive laws . . . .......................... 16 2.1.4 BasicequationsandALEformulation ................. 20 2.1.5 Analyses using theoretical constitutive laws .............. 24 2.1.6 Applicationtorealmaterials ...................... 27 2.1.7 Sizeandtimeeffects........................... 31 2.1.8 Concludingremarks ........................... 34 2.2 Arbitrary Lagrangian–Eulerian formulation for hyperelastoplasticity . . . . 37 2.2.1 Introduction ............................... 37 2.2.2 Multiplicative finite–strain plasticity in a Lagrangian setting . . . . 39 2.2.3 Multiplicativefinite–strainplasticityinanALEsetting ....... 41 2.2.4 Numericalexamples ........................... 46 2.2.5 Concludingremarks ........................... 59 3 Non-standard consistent tangent operators 61 3.1 Numerical differentiation for local and global tangent operators in computationalplasticity................................. 62 3.1.1 Introduction ............................... 62 3.1.2 Problemstatement............................ 63 3.1.3 Numerical differentiation . . ...................... 66 i ii Contents 3.1.4 Examples ................................. 70 3.1.5 Concludingremarks ........................... 84 3.2 Consistenttangentmatricesforsubsteppingschemes............. 86 3.2.1 Introduction ............................... 86 3.2.2 Problemstatement............................ 88 3.2.3 The consistent tangent matrix for the substepping technique . . . . 90 3.2.4 Examples ................................. 93 3.2.5 Concludingremarks ...........................104 3.2.6 Appendix: Consistent tangent moduli for substepping with the generalizedmidpointrule...........................104 3.2.7 Appendix: Computationally efficient expression of the consistent tangent moduli for substepping with the backward Euler method . . 107 4 Elastoplastic models for granular materials 111 4.1 PlasticflowpotentialfortheconeregionoftheMRS–Lademodel......112 4.1.1 Introduction ...............................112 4.1.2 MRS–Lademodel ............................112 4.1.3 Modifiedplasticflowpotential .....................115 4.1.4 Corrected plastic flow potential . . . ..................115 4.1.5 Concludingremarks ...........................116 4.2 Numerical differentiation for non-trivial consistent tangent matrices: an applicationtotheMRS-Lademodel ......................118 4.2.1 Introduction ...............................118 4.2.2 Problemstatement............................120 4.2.3 Proposedapproach............................123 4.2.4 Numerical differentiation . . ......................124 4.2.5 Examples .................................125 4.2.6 Concludingremarks ...........................142 4.2.7 Appendix:MRS–Lademodeldefinition ................142 4.3 TheMRS-Lademodelforcohesivematerials..................144 4.3.1 Introduction ...............................144 4.3.2 TheoriginalMRS–Lademodelforcohesionlessmaterials.......144 4.3.3 TheproposedMRS–Lademodelforcohesivematerials........146 4.3.4 Recoveringthemeaningofcapparameters...............149 4.3.5 Yield function ..............................149 4.3.6 Concludingremarks ...........................150 4.4 Consistent tangent matrices for density–dependent finite plasticity models . 151 4.4.1 Introduction ...............................151 4.4.2 Problemstatement............................153 Contents iii 4.4.3 Consistent tangent moduli for density–dependent finite plasticity models...................................157 4.4.4 Examples .................................160 4.4.5 Concludingremarks ...........................176 4.4.6 Appendix:Density–dependentfiniteplasticitymodels ........177 5 An application to powder compaction processes 181 5.1 Introduction....................................182 5.2 Problemstatement................................183 5.2.1 Kinematics ................................183 5.2.2 Constitutive model . . ..........................184 5.2.3 Numericaltime–integration .......................186 5.3 Numericalsimulations ..............................188 5.3.1 Homogeneoustests............................189 5.3.2 Aplainbushcomponent.........................190 5.3.3 Arotationalflangedcomponent.....................192 5.3.4 Amulti–levelcomponent ........................200 5.4 Concludingremarks ...............................203 6 Summary and future developments 205 iv Contents List of Tables 2.1 Some generalized Newtonian fluid models defined in terms of viscosity. A simplified 1D representation of the models is included. (After Huerta and Liu 1988). . . . .................................. 20 2.2 Numerical values of N1and N2for different materials, vane sizes and angular velocities. ..................................... 21 2.3 Numerical values of the analyses of the vane test using theoretical constitutive laws. . . .................................. 26 2.4 Numerical values of the analyses of the vane test applied to Red Mud. . . . 27 2.5 Numerical values of the analyses of the vane test applied to soft clay. . . . . 32 2.6 The proposed ALE approach for hyperelastoplasticity . . .......... 45 2.7 Materialparametersfortheneckingtest..................... 47 3.1 Numericalapproximationstothederivativesoftheflowvector........ 67 3.2 Relative stepsizes that give the same convergence results as analytical derivatives,forthevonMisesperfectplasticityglobalproblem. ......... 80 3.3 Relative stepsizes that give the same convergence results as analytical derivatives,forthevonMisesexponentialhardeningglobalproblem. ...... 80 3.4 The convergence results for sixth load step of the pile problem: (a) 1NDO(h), (b) 2ND-O(h2). .............................. 81 3.5 Relative stepsizes that give the same convergence results as analytical derivatives,forthepileproblem. .......................... 82 3.6 Relative stepsizes that give the same convergence results as analytical derivatives,fortherigidfootingproblem....................... 83 3.7 Time overheads of the numerical approximations, for the rigid footing problem. ........................................ 84 3.8 Rigid footing problem. Relationship between number of global load incrementsandrelativeCPUtime........................... 96 3.9 Triaxial test problem. Relationship between number of global load incrementsandrelativeCPUtime...........................100 4.1 Numericalapproximationstofirstderivatives..................125 v xii List of Figures 4.31 Powder–A triaxial test (150 MPa initial isostatic compaction). Convergence results for different load levels and with (a,c) unsymmetric and (b,d) symmetriclinearsolvers,leftandrightrespectively...............168 4.32 Results of isostatic and uniaxial tests with powder–B material, see table 4.12. Experimental data from Ernst and Barnekov (1994). ..........169 4.33 Powder–B isostatic test. Convergence results for different load levels. . . . . 169 4.34 Frictionless compaction of the flanged component with the elliptic model. Final relative density, η, distribution after a top displacement of 6.06 mm. 170 4.35 Frictionless compaction of the flanged component with the elliptic model. Final distribution of Cauchy stresses at Gauss–point level on the meridian plane q¯ σ—p¯ σ. Blue marks indicate plastic steps and green marks indicate elasticsteps.....................................170 4.36 Frictionless compaction of the flanged component with the elliptic model. Convergence results for different load levels. ..................171 4.37 Traces of the cone–cap yield function on the meridian plane q¯ τ–p¯ τfor different relative densities, η. Powder–A material parameters, see table 4.12. . 172 4.38 Frictionless compaction of the flanged component with the cone–cap model. Distribution of Cauchy stresses at Gauss–point level on the meridian plane q¯ σ–p¯ σfor different load levels. Blue marks indicate plastic steps on the cap region, red marks plastic steps on the cone region and green marks elastic steps. .......................................173 4.39 Frictionless compaction of the flanged component with the cone–cap model. Final relative density, η, distribution after a top displacement of 6.06 mm. . 174 4.40 Frictionless compaction of the flanged component with the cone–cap model. Convergence results for different load levels. ..................175 5.1 Trace of the elliptic yield function on the meridian plane qτ–pτfor different relative densities, η. Parameters of Powder–C material, see table 5.2. . . . . 185 5.2 Dependence of parameters a1(η)anda2(η) on the relative density, η.ParametersofPowder–Cmaterial,seetable5.2..................185 5.3 Plain bush component. Problem definition (after Lewis and Khoei 1998) and computational mesh. . . ..........................190 5.4 Plain bush component. Relative density profile at radius equal to 10.5 mm. 191 5.5 Plain bush component. Relationship between top punch vertical reaction anditsverticaldisplacement. ..........................191 5.6 Plain bush component. Relative density distribution for different top punch movements. Note that two different scales are used. . . . ..........193 List of Figures xiii 5.7 Flanged component. Relationships between the vertical reactions of the punches and their vertical movements: top reaction for top punch compaction, bottom reaction for bottom punch compaction, and top and bottom reactions for double–punch compaction. ..................194 5.8 Relative mass variation during the three compaction processes of the flanged component. The load levels are referred to the punch displacements imposed attheendofeachtest...............................195 5.9 Top punch compaction. (a–d) Relative density distribution for different top punch movements and (e) relative density profile at 1.88 mm from line GD (section 1–1’). ..................................196 5.10 Bottom punch compaction. (a–d) Relative density distribution for different bottom punch movements and (e) relative density profile at 3.9 mm from line GD (section 1–1’). . . . ..........................197 5.11 Double–punch compaction. Relative density distribution for different movements of the top and bottom punches. (a) common density scale and (b–e) different scales. ..................................198 5.12 Double–punch compaction. Relative density profiles at 3.47 mm (section 1–1’) and 1.3 mm (section 2–2’) from the line GD, and at radii equal to 9.37 mm (section 3–3’) and 8.77 mm (section 4–4’). ..............199 5.13 Multi–level component. Problem definition (after Khoei and Lewis 1999) and computational mesh. . . ..........................201 5.14 Multi–level component. Relative density distribution for three different compactionprocesses. ..............................202 5.15 Multi–level component. Relative density profile at radius equal to 30 mm for the three different compaction processes. . . . ..............202 xiv List of Figures Acknowledgements I would like to thank Antonio Rodr´ıguez and Antonio Huerta, my supervisors, for their suggestions and support during this research. I am also thankful to my colleagues of the Applied Mathematics Department III and the Faculty of Civil Engineering of the Technical University of Catalonia, in particular those of my research area. I would not been able to arrange teaching and research during this last four years without them. I am in debt with Sergio Oller, who introduced me in the powder compaction field, and Francisco Armero, who received me at Berkeley during my firsts steps into finite strain plasticity. Of course, I am also grateful to my family for their patience, support, and love, with a particular mention to Paloma, my wife. She is the only one that knows the cost of this thesis. Finally, I wish to thank all the people of ISF for helping me to not forget how the real world is, and all the friends of the Raval adventure for all the good times we had together. Partial financial support to this research has been provided by the Ministerio de Educaci´on y Cultura (grant number: TAP98–0421) and the Comisi´on Interministerial de Ciencia y Tecnolog´ıa (grant number: 2FD97–1206). Their support is gratefully acknowledged. Barcelona Agust´ıP´erez–Foguet October, 2000 xv xvi Chapter 1 Introduction Numerical modelling of problems involving geomaterials (i.e. soils, rocks, concrete and ceramics) has been an area of active research over the past few decades. This fact is probably due to three main causes: the increasing interest of predicting the material behaviour in practical engineering situations, the great change of computer capabilities and resources, and the growing interaction between computational mechanics, applied mathematics and different engineering fields (concrete, soil mechanics ...). This thesis fits within this last multidisciplinary approach. Based on constitutive modelling and applied mathematics and using both languages the numerical simulation of some complex geomechanical problems has been studied. The state of the art regarding experiments, constitutive modelling, and numerical simulations involving geomaterials is very extensive. The thesis focuses in three of the most important and actual ongoing research topics within this framework: 1. The treatment of large boundary displacements by means of Arbitrary Lagrangian– Eulerian (ALE) formulations. 2. The numerical solution of highly nonlinear systems of equations in solid mechanics. 3. The constitutive modelling of the nonlinear mechanical behaviour of granular materials. The three topics have been analyzed and different contributions for each one of them have been developed. The new developments are presented in chapters 2 to 4, which are related, respectively, to the three topics outlined above. Several applications have been included in these chapters in order to show the main features of each contribution. After that, in chapter 5, some of the new developments have been applied to the numerical modelling of cold compaction processes of powders. These processes consist in the vertical compaction through the movement of a set of punches of a fine powder material at room temperature. The process transforms the loose powder into a compacted sample through a large volume reduction. This problem has been chosen as a reference application 1 2 Introduction of the thesis because it involves large boundary displacements, finite deformations and highly nonlinear material behaviour. Therefore, it is a challenging geomechanical problem from a numerical modelling point of view. Finally, in chapter 6, a brief summary of the thesis is presented, together with an overview of the possible future developments. In the following, a general description of the three main research topics is presented. 1.1 Kinematic formulation In many geomechanical problems, the domain of interest is subjected to large boundary displacements. Moreover, these problems usually lead to non-uniform large–strain solutions. An adaptive strategy for space discretization can help to handle these problems. Several alternatives are available (Huerta, Rodr´ıguez-Ferran, D´ıez and Sarrate 1999). Two common approaches in finite element simulations are h-remeshing techniques and ALE schemes. The h-remeshing techniques imply, usually, a high computational effort and the loss of accuracy due to interpolations from the old mesh to the new mesh. For these reasons the thesis focuses in the formulation of different ALE schemes. ALE formulations reduce the drawbacks of purely Lagrangian or Eulerian formulations. In a Lagrangian formulation, mesh points coincide with material particles. Each element contains always the same amount of material and no convective effects are generated. In this case the resulting governing equations are simple, but it is difficult to deal with large deformations. On the other hand, an Eulerian formulation takes a fixed mesh, and the particles move through it. Now important convective effects appear due to the relative motion between the grid and the particles, but it is possible to simulate large strains. In ALE formulations, the mesh and material deformations are uncoupled, and a computationally efficient compromise between Lagrangian and Eulerian formulations is achieved. The ALE approach was first proposed for fluid problems with moving boundaries (Don´ea, Fasoli-Stella and Giuliani 1977, Donea 1983, Hughes, Liu and Zimmermann 1981, Huerta and Liu 1988). Nowadays, in some fields of solid mechanics, such as the modelling of forming processes, ALE fluid–based formulations are widely used. In the thesis, an ALE fluid–based formulation has been applied to quasistatic and dynamic simulations of the vane test for soft materials, see section 2.1 (P´erez-Foguet, Ledesma and Huerta 1999). This particular application is an example of how an ALE formulation allows to manage the movement of the boundary in a straightforward manner. Several detailed analyses related with the test are included in section 2.1. On the other hand, the ALE approach has been extended and successfully employed in nonlinear solid mechanics (Liu, Belytschko and Chang 1986, Benson 1986, Hu´etink, Vreede Introduction 3 and van der Lugt 1990, Ghosh and Kikuchi 1991, Rodr´ıguez-Ferran, Casadei and Huerta 1998) and structural mechanics (Askes, Rodr´ıguez-Ferran and Huerta 1998, Huerta et al. 1999). However, this is still an important research area. Recently, a specific ALE scheme for hyperelastoplasticity has been presented (Armero and Love 2000). This approach, as a previous formulation for hyperelasticity (Yamada and Kikuchi 1993), is based on a total Lagrangian formulation of the problem. In both cases, the distortion on both the material mesh and the spatial mesh have must be kept under control. Here, in section 2.2, a new ALE scheme for hyperelastoplasticity based on an updated Lagrangian approach is presented (Rodr´ıguez-Ferran, P´erez-Foguet and Huerta 2000). The deformed configuration at the beginning of the time–step is chosen as the reference configuration; consequently, only the quality of the spatial mesh must be ensured by the ALE remeshing strategy. Several applications of the proposed approach are presented in sections 2.2, 4.4 and in chapter 5. Most of them correspond to powder compaction problems, a field where the ALE approach shows to be crucial for an accurate numerical simulation. 1.2 Nonlinear solvers Numerical simulation in solid mechanics usually leads to highly (geometrically and materially) nonlinear problems. The thesis focuses in quasistatic problems. Two kinds of methods can be applied to solve them: implicit and explicit. Implicit methods are preferred, because the are unconditionally stable. Nevertheless, for large scale problems explicit methods are widely used, because the computational effort is lower. The adequacy of explicit solutions for simulations in solid mechanics is still a subject of debate (Owen, Peri´c, de Souza Neto, Yu and Dutko 1995). It is expected that implicit methods will become the standard approach for large scale problems (like during the later eighties they became for 2-D elastoplastic problems). The thesis focuses on implicit methods. Although various nonlinear solvers may be used, see for instance the recent works of Alfano, Rosati and Valoroso (1999) and Sloan, Sheng and Abbo (2000), a common choice is the full Newton-Raphson method. In the context of solid mechanics, the consistent tangent operators ensure quadratic convergence of the Newton-Raphson method. The concept of consistent tangent operator was introduced for simple elastoplastic models (Simo and Taylor 1985, Runesson, Samuelsson and Bernspang 1986). Since then, it has been systematically applied to a broad class of constitutive models of inelastic behaviour (Simo, Kennedy and Govindjee 1988, Ramm and Matzenmiller 1988, Simo 1992, Hofstetter, Simo and Taylor 1993, de Souza Neto, Peri´c and Owen 1994, Li 1995, Crisfield 1997, Simo 1998, Armero 1999, Belytschko, Liu and Moran 2000)). In section 4.4, the expression for finite strain density–dependent plastic models is developed (P´erez-Foguet, Rodr´ıguez-Ferran and Huerta 2000b). Several applications of these models to powder compaction processes are shown in sections 2.2, 4.4 and 4 Introduction in chapter 5. Consistent tangent operators depend on the numerical scheme used for the time– integration of the constitutive equations. The expressions of the consistent tangent operators for the most usual time–integration schemes (such as the backward Euler method or the midpoint rule) are available elsewhere. Other time–integration schemes are those based on substepping techniques (Sloan 1987, Sloan and Booker 1992, Potts and Ganendra 1994). They are less common because, among other things, they lack of the corresponding consistent tangent operators. For this reason, in section 3.2, the consistent operators for different substepping techniques are presented. Moreover, one of them is applied together with an adaptive time–integration strategy. This approach can result in a large reduction of the computational cost for some complex geomechanical problems. This fact is illustrated with the simulation of the rigid footing problem on a frictional material in section 3.2. Independently of the time–integration scheme, there are some difficulties in the analytical definition and the computation of consistent operators for non-trivial constitutive laws. Because of this, in some cases the full Newton-Raphson method is abandoned and specific nonlinear solvers are devised to integrate the constitutive model (Pramono and Willam 1989b, Etse and Willam 1994, Jeremi´c and Sture 1997, Macari, Weihe and Arduino 1997). However, quadratic convergence is not achieved with these methods, because they are not based on a consistent linearization of all the equations with respect to all the unknowns. Here, numerical differentiation has been applied to compute consistent tangent operators with the goal of precluding these difficulties. The proposed approach (P´erez-Foguet, Rodr´ıguez-Ferran and Huerta 2000d) is presented in section 3.1. Several applications are shown in different parts of the thesis: in section 4.2 it is applied to a work hardening cone–cap model for sands, in sections 2.2, 4.4 and in chapter 5 to density– dependent plastic models for powder compaction simulations, and in section 3.2 to several problems solved with the substepping scheme. 1.3 Constitutive modelling Two different approaches are usually followed to model the mechanical behaviour of granular materials: micromechanical and macroscopic. The micromechanical approach consists in modelling each particle separately and computing the behaviour at macroscopic level through the relative interaction between many particles (Borja and Wren 1995, Wren and Borja 1997). The macroscopic approach is based on model the material as a continuous medium. Several continuous constitutive laws have been used for granular materials (typically elastoplastic and viscoplastic models). Elastoplastic models are the most used for isothermal modelling. In fact, nowadays, elastoplastic models are widely used to model many different geomaterials in small strain problems, ranging from virgin noncohesive Introduction 5 sand (Sture, Runesson and Macari-Pasqualino 1989, di Prisco, Nova and Lanier 1993) to hard rocks and concrete (Pramono and Willam 1989a). Recently, several finite strain hyperelastic–plastic models specific for geomaterials have also been developed (Simo and Meschke 1993, Meschke, Liu and Mang 1996, Oliver, Oller and Cante 1996, Borja and Tamagnini 1998, Callari, Auricchio and Sacco 1998). Moreover, an important issue related with the numerical modelling of granular materials is the mathematical notion of well–posedness. This requires the use of regularized constitutive laws to model properly failure situations, see for instance Askes (2000) and references therein. This is an actual important area of research. The thesis focuses in classical elastoplastic modelling. Two different approaches are followed: a work hardening–softening cap model for small strain problems and density– dependent hyperelastic–plastic models for finite strain simulations. Two issues are common to both approaches: the treatment of non-smooth plastic equations and the proper computation of tangent operators. In many elastoplastic models (such as Tresca, Mohr Coulomb and cone–cap models) the yield function is defined by parts, leading to a non-smooth transition of the flow equations at the intersection. The problem is solved, typically, by means of corner return– mapping algorithms based on the Koiter’s rule, see for instance Simo et al. (1988) and references therein. However, in order to reduce the computational effort of the time– integration algorithm, smoothing approaches are preferred in some cases. Two alternatives can be devised: smoothing the yield function (and as a result the flow equations) or smoothing only the flow equations. In both cases corner return–mapping algorithms are not needed. The von Mises–Tresca model (Miehe 1996, P´erez-Foguet and Armero 2000) and the Rounded Hyperbolic Mohr Coulomb model (Sloan and Booker 1986, Abbo and Sloan 1995) are examples of the first alternative. Here, both models are used in several applications in section 2.2 and in chapter 3. The second alternative, to smooth only the flow equations, is more recent (Macari, Runesson and Sture 1994, Macari et al. 1997). In sections 4.1 (P´erez-Foguet and Huerta 1999) and 4.3 (P´erez-Foguet and Rodr´ıguez-Ferran 1999) the second alternative has been applied to work hardening–softening cap models and in section 4.4 to a density–dependent hyperelastic–plastic model. On the other hand, especial attention has been dedicated to the computation of tangent operators for both types of models in sections 4.2 and 4.4. Both cases are examples of non-trivial elastoplastic models with coupled dependence between the plastic strain flow and the hardening–softening laws (which are based on the plastic work and the relative density respectively). As a result of the use of full tangent operators together with the new developments presented in chapters 2 and 3 (ALE formulations, numerical differentiation and consistent operators for substepping techniques) it has been possible to simulate numerically and analyze three complex geomechanical problems involving granular materials: 12 ALE formulations Figure 2.2: Measured stress distributions at vane blades (after Menzies and Merrifield 1980). or for a vane of height equal to twice its diameter su=T 3.66D3and Th Tv =1 6.(2.1.3) The stress distributions obtained from experiments or from numerical analyses are partly different from the assumptions above considered. For instance, figure 2.2 shows the stress distributions on the failure surfaces measured in the blades of an instrumented vane (Menzies and Merrifield 1980), and it may be seen that the distribution of shear stresses on the top is very different from the uniform assumption. Numerical results using an elastic constitutive law (Donald et al. 1977) already suggested a nonlinear distribution of stresses in the top of the vane (figure 2.3a). From these results Worth (1984) proposed a polynomial function to represent the shear stress distribution τ=τ(r) at the top and bottom surfaces, and therefore τ=sur D/2n and Th=πD3su 2(n+3) (2.1.4) where ris indicated in figure 2.3a. Worth suggested a value of n= 5 for London clay, based on the results from Menzies and Merrifield (1980). For this value, the torque ratio becomes smaller: Th/Tv=1/16. Hence, the contribution of the horizontal failure surfaces to the total torque seems to be less significant in practice. That is, almost 94% of the resistance to torque is provided by the vertical failure surface. As a consequence of that, classical expressions to obtain suwould underestimate the actual value of shear strength and that has been reported by some authors (Worth 1984, Eden and Law 1980). Equation (2.1.4) is quite general as according to the value of ndifferent stress distributions for the top and bottom of the vane can be considered. From that results it seems that a value about n= 5 could be appropiate. However, some authors have confirmed recently values close to n= 0 for different soils (Silvestri et al. 1993), which corresponds again to a uniform stress distribution. A finite element analysis presented by Griffiths ALE formulations 13 Figure 2.3: Shear stress distributions on sides and top of vane obtained from numerical simulations. a) Elastic model (after Donald et al. 1977) b) Using an elastoplastic model and a strain softening model including anisotropy (after Griffiths and Lane 1990). and Lane (1990) confirms that for elastoplastic materials the shear stress can be close to a constant value on the top of the vane (figure 2.3b). They also showed an elastic analysis which is consistent with that presented by Donald et al. (1977). Therefore, the value of nwill depend on the stress state reached on the top and bottom of the vane, and that is difficult to predict in advance. This conclusion assumes that soil is isotropic, which is not always the case. When the soil is anisotropic, the interpretation of the test becomes more difficult, as, for instance, maximum shear stress can be reached in the vertical surface whereas the situation on the top is still elastic. As the result used from the test is the peak of the curve torque–rotated angle, which is in fact an integral of all these stresses, it is difficult to distinguish all these effects from just one measured value. As the vane includes vertical and horizontal failure surfaces, some attempts have been made to identify anisotropy by means of vanes with different dimensions and shapes in order to estimate Th and Tvseparately (Aas 1965, Wiesel 1973, Donald et al. 1977). Bjerrum (1973) proposed a correction factor to account for the anisotropy that has been critiziced in some cases (Garga and Khan 1994, Tanaka 1994). When the soil is isotropic and is not strain softening, as the maximum shear strains are produced at r=Rv, where Rvis the vane radius, it is expected to reach the maximum shear stress at the vertical surface failure. If this value is kept constant, then plastification of the top and bottom vane will occur and the peak measured torque will correspond to a uniform distribution of shear stresses in all surfaces. However, if the soil has a strain 14 ALE formulations softening constitutive law, the shear stress on the vertical failure will decrease and the peak torque will correspond to an intermediate situation and n>0. Moreover, when strain softening occurs, the shear stress is not constant, which makes the result of the vane test insufficient to estimate su. These arguments are consistent with the conclusions obtained by de Alencar et al. (1988) in a 2D numerical analysis. They simulated the vane using different strain softening constitutive laws (but all of them with the same peak shear strength). The torque–rotation curve was totally dependent on the constitutive law employed. Numerical simulations presented by Griffiths and Lane (1990) (figure 2.3b) present the same dependence. All those results showed the influence of progressive failure on the final interpretation of the test. A consequence of all the works involved in the study of the interpretation of the vane is that the complete stress–strain curve of the material and its anisotropy must be known in advance in order to explain correctly the results of the test. However, for isotropic soft materials it seems to be an appropiate test, and the vertical failure surface would be predominant in that case. Time effects The influence of time on the results of the test has two different aspects: the delay between insertion and rotation of vane, and the rate of vane rotation. The disturbance originated by the vane insertion and the consolidation following that insertion are difficult to predict in general. There are a few experimental studies about these effects. They suggest that in order to reduce the vane insertion effects, blade thickness related to vane size must be as small as possible (Rochelle, Roy and Tavenas 1973, Torstensson 1977). On the other hand, the delay on carrying out the test after vane insertion increases the measured shear strength, due to the dissipation of pore water pressures originated by the insertion and also due to thixotropic effects. This effect is usually not considered in vane analyses, but there is experimental evidence on the high pore pressure developed by vane insertion (Kimura and Saitoh 1983) and on the microestructural changes due to thixotropic phenomena (Osipov, Nikolaeva and Sokolov 1984). Results from Torstensson (1977) show that within five minutes after insertion measured shear strength does not change. It must be pointed out that both effects depend on the type of clay involved in the test (sensitivity, consolidation coefficient cv, etc.). Fabric disturbance due to insertion reduces the true undrained strength in about 10%, but if consolidation after insertion is permitted a 20% increase on strength is produced (Chandler 1988). The standard vane test is usually performed 1 minute after the insertion of the blades, which is the maximum delay value suggested by Roy and Leblanc (1988). In that case, no consolidation is allowed. The effect of the rate of vane rotation on the interpretation of the test is also important. The standard rate is 6 — 12◦/minute. That produces failure in about 30 — 60 seconds, a ALE formulations 15 Figure 2.4: Shear stress–angular rotation obtained using different testing rates on B¨ackebol clay, Sweden (after Torstensson 1977). shorter time than in classical triaxial tests or shear tests. Due to this difference, undrained strengths from vane tests are overestimated if compared with that obtained from classical laboratory tests. This effect can be compensated with the underestimation of suprovided by other effects (fabric disturbance, stress uniformity,...), but their magnitude is difficult to estimate. As the vane is an undrained test, some recommendations regarding a minimum angular velocity are defined in the codes of practice. Assuming a time to failure of one minute, undrained conditions can be assured if the consolidation coefficient of the material is : cv≤3.5·10−2cm2/s (Chandler 1988), which is usually the case when soft clays are tested. However, as in other undrained tests, measured strength increases with the velocity of the load application and this effect leads to difficulties in the interpretation of results. To consider that, Bjerrum (1973) proposed a reduction of the shear strength measured with the vane according to the plasticity index of the clay, as the usual values of shear strength obtained in the laboratory correspond to slower experiments than field vane. Measurements of vane shear strength for different velocities of rotation have been published by Wiesel (1973) and Torstensson (1977). Both presented a potential expression between shear strength, su, and angular velocity, ω, (or time to failure) from interpolation of their results (su)vane =k1ωk2(2.1.5) where k1and k2are constants. The value of k2ranged from 0.02 to 0.07. Figure 2.4 presents shear strength versus rotation angle for different durations of the test (Torstensson 1977). The comparison of results from vanes of different shapes and different strain rates has been difficult as these effects are only related to the shear strength on the basis of empirical relations which may depend on the soil considered. A reference shear velocity: 16 ALE formulations v=ωRvas a variable to compare results from different vanes was proposed by Perlow and Richards (1977). They obtained an almost linear relationship between vane shear strength and shear velocity vfor two marine sediments, but they did not have enough experimental data to propose a definite relationship. In fact some results reported by other authors (Tanaka 1994), show no influence of the vane radius on the measured shear strength. These differences may be due to side effects as sampling disturbance or stress relief for soils used in laboratory vane tests, which would reduce the apparent shear strength (Kirkpatrick and Khan 1984). However, in general, that is taken into account when estimating su. This is still a controversial issue, and it will be considered later using the formulation presented in following sections. The shortcomings presented have been extensively studied by many authors, but still there are contradictions on results and on interpretations of vane measurements. This is due to the fact that vane test is a model test rather than an element test (Morris and Williams 1993). Apart from that, other effects on vane strength have been rarely studied: for instance, influence of the stress state and K0(coefficient of lateral earth pressure at rest) on the result (Law 1979, Garga and Khan 1994). As a consequence of the drawbacks above mentioned, the vane test seems indicated for very soft isotropic materials. Also the dominant failure surface is the vertical one. Therefore, it is reasonable to perform numerical simulations by means of two–dimensional analyses, in order to study rate effects and stress distributions. Hence the use of plastic fluid constitutive laws and fluid mechanics equations may be appropiate to analyze time effects and to deal with this soft isotropic materials. The disturbance due to vane insertion and the 3D effects of the test are not considered within this approach. 2.1.3 Constitutive laws Soft materials have been successfully modelled by means of fluid constitutive laws to simulate classical geomechanics problems like landslides or debris flows (Vulliet and Hutter 1988, Rickenmann 1991, Gili, Huerta and Corominas 1993). The model is defined in terms of a shear stress–shear strain rate relationship, instead of a stress–strain one. The study of strain rate effects on soil behaviour has been considered in many works (Berre and Bjerrum 1973, Komamura and Huang 1974, Cheng 1981, Mesri, Febres-Cordero, Shields and Castro 1981, Leroueil, Kabbaj, Tavenas and Bouchard 1985). Some of them have proposed constitutive laws relating stresses with strains and strain rates to account for time influence, based on a visco–elasto–plastic theory. Bjerrum (1973) indicates that time effects in soft clays are associated with the cohesive component of the shear strength which is of a viscous nature; the frictional component of the shear strength would be further mobilized. Therefore it is not absurd to study the vane test assuming that the material involved is a plastic fluid, particularly at the beginning of the test, when viscous ALE formulations 17 Figure 2.5: Rheological state of soil in accordance with water content for some Japanese clays (after Komamura and Huang 1974). effects are definite. These effects depend on the type of clay considered, and plasticity index has been used historically to distinguish between different behaviours of normally consolidated clays. However, Tavenas and Leroueil (1980) propose to use the limit liquid instead, because it requires only a single test and in fact the resulting correlations are essentially the same. Let us assume that the clay involved in the test is saturated and normally consolidated. When viscous effects are studied in detail, it is found that for a particular soil (that is for a particular liquid limit wl), the soil behaviour depends on the water content w,as shown in figure 2.5, from Komamura and Huang (1974). According to them, when w>w l the behaviour is viscous, that is, close to a newtonian fluid; whereas when w<w vp the behaviour is visco–elasto–plastic. The value wvp was defined as viscoplastic limit, between plastic and liquid limits. Results obtained in a torsional hollow cylinder are reproduced in figure 2.6, from Cheng (1981). The material was a Mississippi Buckshot clay with a water content close to its plastic limit. It may be seen that its behaviour is consistent with the trends established by Komamura and Huang (1974). However, further works have shown that some clays may have a visco–elasto–plastic behaviour with water content above their liquid limit 18 ALE formulations Figure 2.6: Effect of strain rate on undrained shear stress obtained using torsional hollow cylinder (after Cheng 1981). (Bentley 1979, Torrance 1987, Locat and Demers 1988). Curves obtained by Locat and Demers (1988) using a viscosimeter device are presented in figure 2.7. Note, nevertheless, that the shear rate range is different for the results reproduced in figure 2.6 and for those indicated in figure 2.7. Some consequences of this difference will be treated later. The works indicated show that fluid constitutive laws for modelling the soil behaviour have been successfully employed to account for the time effects which are supposed to be important when water content is high, but it is difficult to establish a particular behaviour for each soil state in advance. The general constitutive equation to be used in this work is a general relationship between stresses and strain rates: σij =f(dij) (2.1.6) where σij is Cauchy’s stress tensor and the strain rate tensor is defined as dij =1 2∂vi ∂xj +∂vj ∂xi(2.1.7) where xiand viare position and velocity vectors respectively. A common expression for equation (2.1.6) is: σij =−pδij +2µdij (2.1.8) where pis the hydrostatic pressure (tension positive) and µis the dynamic viscosity. Equation (2.1.8) can be rewritten as σd ij =2µdij (2.1.9) ALE formulations 19 Figure 2.7: Shear stress–shear strain rate obtained from viscosimeter experiments with St. Alban–1 marine clay with a salt content of 0.2 g/l; τyis yield stress, ILis liquidity index and µis viscosity (after Locat and Demers 1988). where σd ij is the deviatoric stress tensor. When viscosity is assumed constant, the fluid is called newtonian. A generalized newtonian fluid is defined by a viscosity which depends on the strain rate tensor. Also, some plastic fluids show a “yield stress”, that is, below that value no flow is observed, which is equivalent to an infinite viscosity in terms of equation (2.1.9). However, some authors (Barnes and Walters 1985) indicate that the yield stress does not exist, provided that accurate measurements are performed. That is, “yield stress” is just an idealization of the actual behaviour. Table 2.1 shows a list of constitutive laws available for fluids (Huerta and Liu 1988). From that list, the models by Bingham, Casson and Herschel and Bulkley seem to be most appropiate for soft materials. Those models simulate a “yield stress” value up to which no velocities or displacements occur. According to the yield stress magnitude, it is possible to reproduce effects observed in the actual behaviour of soft materials. For instance, the amplitude of the strain localized zone, the influence of the progressive failure or the brittleness of the material are supposed to be determined by the constitutive law used and specially by the existence of the “yield stress”. The fluid models indicated above are consistent with experimental results like the ones depicted in figures 2.6 and 2.7, and they will be used as constitutive laws in this simulation. However, any relation between stresses and strain rates could be implemented in the formulation. 20 ALE formulations Newtonian µ=µ0 γγ ττ . Carreau µ=(µ0−µ∞)1+(λ˙γ)2(n−1)/2+µ∞ γγ ττ . Bingham µ=∞if τ≤τ0 µ=µp+τ0/˙γif τ>τ 0γγ ττ . Herschel–Bulkley µ=∞if τ≤τ0 µ=µp˙γn−1+τ0/˙γif τ>τ 0γγ ττ . Casson µ=∞if τ≤τ0 √µ=√µp+τ0/˙γif τ>τ 0γγ ττ . Table 2.1: Some generalized Newtonian fluid models defined in terms of viscosity. τ= σd ijσd ij/2and ˙γ=2dijdij. A simplified 1D representation of the models is included. (After Huerta and Liu 1988). 2.1.4 Basic equations and ALE formulation Basic equations Fluid movement is described by two basic equations: mass conservation equation and equilibrium equation. These equations are, respectively: ∂ρ ∂t +ρ∂vi ∂xi =0 in Ω,(2.1.10a) ρ∂vi ∂t +ρcj ∂vi ∂xj =ρbi+∂σji ∂xj in Ω ,(2.1.10b) where ρis the density of the material, tis time, cj=ˆvj−vjwhere ˆvjis the velocity of the reference system, biis the mass forces vector and Ω is the domain of study. Also, repeated index means summation. As incompressible flow is assumed, density is constant and expression (2.1.10a) leads to divv= 0 which is equivalent to the undrained condition assumed in the standard field vane test. Introducing this result and the constitutive law presented above in (2.1.10b) gives ρ∂vi ∂t +ρcj ∂vi ∂xj =ρbi−∂p ∂xi +∂ ∂xjµ[˙γ]∂vi ∂xj +∂vj ∂xiin Ω (2.1.11) ALE formulations 21 ρ R∗τ∗ω∗N1min(N2) Kg/m3mPa s−1 Red Mud 1200 0.013 126 0.01 1.6E-7 2.E-2 1200 0.013 126 0.21 7.1E-5 2.E-2 Soft Clay 1400 0.0325 2000 0.0017 2.1E-9 4.E-2 1400 0.0325 2000 0.0035 9.1E-9 4.E-2 1400 0.0325 2000 0.0070 3.6E-8 4.E-2 Table 2.2: Numerical values of N1and N2for different materials, vane sizes and angular velocities. where µ[˙γ] is the dynamic viscosity, function of the shear strain rate, ˙γ=2dijdij. The boundary conditions applied are vx=vy= 0 at the outer boundary and vx= ωr sin(θ◦+ωt)andvy=−ωr cos(θ◦+ωt) at the blade contours, with ωthe angular velocity of the vane and (r, θ◦) the polar coordinates of the blade nodes at t= 0. The initial stresses and velocities are considered zero in all the points and in the boundaries. In order to find out the parameters that govern the problem, a set of dimensionless variables ¯x=x/R∗,¯p=p/τ∗,¯v=v/R∗ω∗and ¯ t=tω∗(2.1.12) is substituted in equation (2.1.11), where R∗,τ∗,ω∗are characteristic length, stress and angular velocity of the test respectively. Usually, the vane radius, Rv, is adopted for R∗, the angular velocity of the test for ω∗, and the material yield stress for τ∗. Then equation (2.1.11) is transformed into a dimensionless expression: 1 Ne∂¯vi ∂¯ t+¯cj ∂¯vi ∂¯xj=−∂¯p ∂¯xi +∂ ∂¯xj1 NeRe[˙γ]∂¯vi ∂¯xj +∂¯vj ∂¯xiin Ω (2.1.13) with Re and Ne equal to Reynolds and Newton number. The Reynolds number is related to viscous forces and the Newton one to inertial forces. They are defined as Ne =τ∗ ρ(R∗ω∗)2and Re =ρω∗(R∗)2 µ.(2.1.14) The influence of Re and Ne in equation (2.1.13) is in the form of N1=1 Ne =ρ(R∗ω∗)2 τ∗and N2[˙γ]= 1 NeRe[˙γ]=µ[˙γ]ω∗ τ∗,(2.1.15) and their characteristic values for some typical vane tests are shown in table 2.2. As N1 is much less than N2and accelerations usually are not large enough to compensate this difference, the inertial terms can be neglected and the problem becomes quasi-static. In these cases the problem depends just on N2, and therefore the test will be independent of the vane radius or the fluid density. The use of those dimensionless numbers may be useful when comparing different vanes, and they will be considered later to account for size and time effects. 28 ALE formulations 0 100 200 300 400 500 600 0 50 100 150 200 250 300 350 400 γγ [s-1] ττ [Pa] Experimental data Casson Herschel - Bulkley Bin g ham Zoom 75 100 125 150 175 0 0,1 0,2 0,3 0,4 . Figure 2.13: Shear stress [Pa] versus shear strain rate [s−1] for the Red Mud constitutive laws. ported by Nguyen and Boger (1985). Vane tests performed with vane radius equal to 1.3 cm and angular velocity equal to 0.1 cycles/min gave (su)vane of 126 Pa. Viscosimeter experimental data and least square aproximation by means of three different constitutive models are presented in figure 2.13. The models used are Bingham, Casson and HerschelBulkley. Table 2.4 shows some numerical results of the analyses in order to compare the models employed (with R∗=1.3cm,ω∗=0.1 cycles/min and τ∗= 133.5 Pa). Note that again r=1.01 is the radius at which the maximum shear strain rate and shear strength is produced irrespective of the model used. The applied torque values are very different. The reason for that is the range of shear strain rate mobilized during the test: 0 to 40 (dimensionless values) which correspond to ˙γ=0to0.42 s−1. In the zoom of figure 2.13, the different stress level for this range of ˙γis clearly highlighted. That zoom shows that a correct interpolation must be used to perform a correct analysis. The value of (su)vane calculated with Bingham model simulations, 136 Pa, is quite similar to (su)vane measured by Nguyen and Boger (1985), 126 Pa. Note that Bingham model has been approximated from experimental viscosimeter data at low shear strain rate, but still 10 times larger than the shear strain rate mobilized in the vane test. In general, extrapolation of experimental values from viscosimeters must be used carefully because the range applied in a vane test is very small when compared to that of a viscosimeter. ALE formulations 29 Figure 2.14: Summary of undrained rate effects in isotropically consolidated soils of different composition (after Lacasse 1979). Continuous line corresponds to the case analized in numerical simulations. Soft Clay Soft clays have been extensively studied also by many researchers, but usually a soil mechanics point of view has been employed to define their behaviour. As a consequence of that, constitutive laws suggested for soft clays are presented usually in terms of a stress–strain relationship instead of a stress–strain rate one, even for clays with a high liquidity index. Most of the studies have been carried out using the conventional triaxial test, which can take from a few minutes to one or two hours. Hence rate effects are expected to be less important than in the vane test, where failure is reached in about one minute. Nevertheless, experimental results obtained in viscosimeters have also been published (figures 2.5 and 2.7), although its shear strain range is different from that used in the vane. In order to obtain a general form of a shear stress–strain rate relationship, information about duration of triaxial tests and rate effects has been used (Torstensson 1977, Lacasse 1979, Hight, Jardine and Gens 1987). A logarithmic relation between shear stress and time to failure has been proposed by many researchers. Typical curves from Lacasse (1979) are presented in figure 2.14. Using strain rate as main variable, this relationship may be expressed as τ/τr=alog( ˙γ/˙γr)+b, (2.1.19) where τris the failure shear strength at reference time, ˙γris the failure strain rate at same reference time and a,bare constants. If the test is performed at a constant strain rate, ˙γr can be computed from the ratio shear deformation at failure vs. time to failure. Equation (2.1.19) shows the effect of increase in shear strength when the test is faster, which is a known behaviour for soft clays. 30 ALE formulations 1,3 1,4 1,5 1,6 1,7 0 10203040 γγ / ωω∗ ττ / ττ∗ Bin g ham Carreau Lo g arithmic . Figure 2.15: Dimensionless shear stress versus dimensionless shear strain rate for the soft clay constitutive laws. From the curves of figure 2.14, a clay with the following constitutive law has been considered (assuming the reference time equal to 140 min and a 3% of the failure shear deformation at that time): Logarithmic model : τ/τ∗=0.13 log(˙γ/ω∗)+1.39 ,(2.1.20a) where τ∗is a shear strength reference value and ω∗the rotation velocity that has been fixed to 12◦/min. Also a Bingham and a Carreau models have been considered interpolating the Logarithmic one in the interval ˙γ/ω∗∈[0.1,40]: Bingham model : τ/τ∗=1.46 + 0.0042 ˙γ/ω∗,(2.1.20b) Carreau model : µω∗/τ∗= 100 1 + 7200(˙γ/ω∗)2−0.48 .(2.1.20c) As in previous sections, Bingham and Logarithmic models have been approximated with a initial dimensionless viscosity of 100, to avoid infinite viscosity values. Figure 2.15 shows the shear strain–strain rate curves corresponding to these models. Note that the Carreau model interpolates so well the Logarithmic model that no difference can be observed over the interesting interval. The results obtained with these models are represented in figures 2.16, 2.17 and 2.18, and they are compared in table 2.5. As could be expected, Logarithmic and Carreau results are exactly the same, and the shear strain localizated zone originated using the Logarithmic and Carreau models is slightly wider than the one obtained with the Bingham model. However, the value of the radius at which maximum shear strain rate occurs in an intermediate plane is again 1.01 for all the models. Calculated torques are the same for the three models, so the (su)vane associated to the simulations will be also the ALE formulations 31 0.143E+01 0.429E+01 0.714E+01 0.100E+02 0.129E+02 0.157E+02 0.186E+02 0.214E+02 0.243E+02 0.271E+02 0.300E+02 0.329E+02 0.357E+02 0.386E+02 0.571E-01 0.171E+00 0.286E+00 0.400E+00 0.514E+00 0.629E+00 0.743E+00 0.857E+00 0.971E+00 0.109E+01 0.120E+01 0.131E+01 0.143E+01 0.154E+01 (b) (a) Figure 2.16: Shear strain rate and shear stress distributions using the soft clay constitutive laws: Bingham, a), and Logarithmic, b). same. Nevertheless if these laws are extrapolated to other strain rate ranges, the predicted behaviour could be very different, as shown in figure 2.17. After comparing experimental results (figures 2.6 and 2.7) with figure 2.19, it must be pointed out that Bingham models must be used with caution. It is very important to choose the correct range of applicability when one desires to approximate the actual behaviour of soft clays by Bingham models. However Carreau and Logarithmic models seem to capture better the changing scale between the results of vane test and viscosimeter ones. Therefore, Carreau and Logarithmic models seem to give continuity from the results of the triaxial undrained tests to that obtained from viscometers. 2.1.7 Size and time effects The dimensionless numbers defined in (2.1.15) are very useful to analyze the factors that influence the results of the test. For instance, size effects that have traditionally been considered from an empirical point of view, can be studied in a more objective manner with this approach. A relationship between vane shear strength and shear velocity computed as v=ωRvwas depicted from experimental results by Perlow and Richards (1977). However, 32 ALE formulations 0,0 0,2 0,4 0,6 0,8 1,0 v / R∗ωω∗ Lo g arithmic Bin g ham 0 2 4 6 8 10 0,0 0,5 1,0 1,5 2,0 2,5 3,0 r / R∗ γγ / ω ω∗ . Figure 2.17: Dimensionless velocity and shear strain rate between blades using soft clay constitutive laws. Logarithmic Carreau Bingham maxθ=0◦(v/R∗ω∗)1.00 1.00 1.00 maxθ=45◦(v/R∗ω∗)0.88 0.88 0.94 r/R∗|maxθ=45◦(v/R∗ω∗)0.89 0.89 0.94 maxθ=0◦(˙γ/ω∗)38.39 38.34 39.84 maxθ=45◦(˙γ/ω∗)6.70 6.73 10.00 r/R∗|maxθ=45◦(˙γ/ω∗)1.01 1.01 1.01 maxθ=0◦(τ/τ∗)1.60 1.60 1.63 maxθ=45◦(τ/τ∗)1.50 1.50 1.50 r/R∗|maxθ=45◦(τ/τ∗)1.01 1.01 1.01 T/τ∗(R∗)29.81 9.82 9.85 Table 2.5: Numerical values of the analyses of the vane test applied to soft clay. that result was based on just few data, and it was not definite. In fact, other authors (Tanaka 1994) have reported completely different results. Based on many laboratory small vane tests and field vane tests on Japanese clays, they did not find any sustantial difference between vane sizes, when the same angular velocity, ω, was used. In fact, according to (2.1.15) the mathematical problem is controlled by N2=µω∗/τ∗, as N1is usually very small. Thus the problem does not depend on the vane size, but on the rotation velocity, the yield stress and the viscosity. For a particular soil, ωis the fundamental parameter. That applies for soft materials for which the constitutive laws used are reasonable. And it is confirmed by experimental evidence, as the results presented by Tanaka (1994). An attempt has also been made to use the formulation presented above to reproduce the effect of vane rotation (and therefore the time to failure) on the shear strength provided by the vane. ALE formulations 33 5 6 7 8 9 10 11 0,000 0,005 0,010 0,015 0,020 t ωω∗ T / ττ∗ (R∗)2 Carreau & Lo g arithmic Bin g ham Figure 2.18: Dimensionless torque versus dimensionless time for the soft clay constitutive laws. 0 1 2 3 4 03691215 γγ [s-1] ττ / ττ∗ Bin g ham Lo g arithmic Carreau . Figure 2.19: Dimensionless shear stress versus shear strain rate (1/s) for the soft clay constitutive laws. If different rotation velocities are used in the simulation, a relationship similar to equation (2.1.5) is expected to be found. Figure 2.20 presents the results obtained using two models employed previously (Bingham and Logarithmic). If the computed torque is expressed as (su)vane/τ∗=k1ωk2,(2.1.21) the values obtained for Bingham model are (k1)Bin =1.358 and (k2)Bin =0.052 ,with r2=0.95 ,(2.1.22a) where r2is the correlation coefficient, and for Logarithmic model are (k1)Log =1.395 and (k2)Log =0.037 ,with r2=0.99 .(2.1.22b) These values are consistent with the experimental ones provided Wiesel (1973) and Torstensson (1977) indicated in equation (2.1.5). In particular, Logarithmic model seems to 34 ALE formulations 9,0 9,5 10,0 10,5 11,0 11,5 12,0 0 20406080100 ωω [o/min] T / τ τ∗(R∗)2 Bin g ham Lo g arithmic Potential ( Bin g ham ) Potential ( Lo g arithmic ) Figure 2.20: Simulation results and potential interpolation of relationship between dimensionless torque and angular velocity (◦/min) for different soft clay constitutive laws. be specially designed to approximate experimental relationships expressed as equation (2.1.21). Therefore, time effects seem to be simulated correctly by means of this model based on fluid mechanics principles. This was expected as those effects were measured on soft clays where viscous phenomena are supposed to be important. In conventional vane tests, viscous forces dominate inertial ones, and the problem depends on N2(equation 2.1.15). When a fluid constitutive law is used, the torque increases with time up to a limit value, as in figure 2.18. It is not possible to reproduce a peak in the torque–time curve in this way, unless inertial forces become important. For usual vane velocities, this is not the case, but some measurements at high velocities have also been reported in the literature (Torstensson 1977). Figure 2.21 shows the effect of inertial forces on the shape of the torque–time curve, in terms of N1dimensionless value, using the Bingham–1 model from figure 2.9. Thus, even for that model, at high rotation velocities, inertial effects produce a peak on that curve. This is consistent with the measurements presented in figure 2.4 (Torstensson 1977), where the peak of the curve torque–rotation (or time) is more pronounced when the test is faster, as inertial forces become more important. For the normal velocity range however, a Bingham model can not produce such a peak, and an explanation in terms of softening of the material could be appropiate. 2.1.8 Concluding remarks A simulation of the vane test using an Arbitrary Lagrangian–Eulerian formulation and appropiate for soft clays has been presented. As the dominant failure surface is the vertical one, a 2D analysis has been useful enough to study stress distributions. Also, the use of fluid mechanics principles and fluid mechanics constitutive laws have permitted the characterization of time effects in a natural manner, as velocities instead of displacements are the main variables. ALE formulations 35 0,0 0,5 1,0 1,5 2,0 2,5 3,0 3,5 4,0 4,5 0,000 0,002 0,004 0,006 0,008 0,010 t ωω∗ T / Tquasi-static N1 = 10-1 N1 = 10-2 N1 = 103 N1 = 10-4 ω = ct. Figure 2.21: Dimensionless torque versus dimensionless time for different inertial forces, using Bingham–1 model. The mathematical problem is governed by two dimensionless numbers. They are related to inertial forces (Newton number) and to viscous forces (Reynolds number). Two tests performed with vanes of different sizes and different rotation velocities, can only be compared by means of these numbers. In most cases, at a typical angular velocity and with usual vane dimensions, the problem becomes quasi-static and independent from inertial forces. Thus in this case, the problem is independent from vane radius and density, and it is controlled by the value of N2=µω∗/τ∗(µ:viscosity,ω∗: rotation velocity, τ∗: characteristic yield stress). Therefore, rotation velocity should be used as main variable to compare different vanes tested on the same soil, provided that the assumptions considered apply (i.e. when 2D conditions are predominant and soft materials are tested). For soft clays and usual vane conditions, the simulated torque is always increasing with time. Thus a peak in the curve torque–time should be related to other effects as strain softening of the material tested. However, if rotation velocity is increased, inertial effects become more important, and a peak in that curve is always obtained. Stress and strain rate distributions on the failure surface depend on the constitutive laws adopted for the material, as stated in previous works. However, the position of the failure surface (where maximum shear stresses are developed) has been found to be always at 1. to 1.01 times the vane radius. Differences around 10% have been found in the shear stress distribution along the failure surface, depending on the material model. Also, amplitude of the shear band is related to that: Bingham models tend to produce more definite shear bands. On the other hand, a yield stress is clearly obtained when results of the vane test are independent of 36 ALE formulations its rotation velocity. That is, when shear stress is constant irrespective of the shear strain rate reached in the test. To compare triaxial, vane and viscosimeters results, it is necessary to take into account the different shear strain rate mobilized in each test. The same model will give different strengths in each case. Carreau and Logarithmic models seem to reproduce well that change of scale. Also, the effect of rotation velocity in shear strength has been simulated using this approach. In fact, shear strength increase associated to rotation velocity increase is directly related to the increment of shear strength due to shear strain increments. The experimental relation between these variables that has been reported in the literature, has been reproduced by means of this approach. Thus time effects, defined in terms of rotation velocity (or time to failure), have been studied in this manner. Finally, it can be concluded that the use of a fluid mechanics approach has proved to be appropiate for the interpretation of this test when soft materials are involved. ALE formulations 37 2.2 Arbitrary Lagrangian–Eulerian formulation for hyperelastoplasticity The Arbitrary Lagrangian–Eulerian (ALE) description in nonlinear solid mechanics is nowadays standard for hypoelastic–plastic models. An extension to hyperelastic–plastic models is presented here. A fractional–step method —a common choice in ALE analysis— is employed for time–marching: every time–step is split into a Lagrangian phase, which accounts for material effects, and a convection phase, where the relative motion between the material and the finite element mesh is considered. In contrast to previous ALE formulations of hyperelasticity or hyperelastoplasticity, the deformed configuration at the beginning of the time–step, not the initial undeformed configuration, is chosen as the reference configuration. As a consequence, convecting variables is required in the description of the elastic response. This is not the case in previous formulations, were only the plastic response contains convection terms. In exchange for the extra convective terms, however, the proposed ALE approach has a major advantage: only the quality of the mesh in the spatial domain must be ensured by the ALE remeshing strategy; in previous formulations, it is also necessary to keep the distortion of the mesh in the material domain under control. Thus the full potential of the ALE description as an adaptive technique can be exploited here. These aspects are illustrated in detail by means of three numerical examples: a necking test, a coining test and a powder compaction test. 2.2.1 Introduction The Arbitrary Lagrangian–Eulerian (ALE) formulation is a standard approach in large strain solid mechanics to keep mesh distortion and element entanglement under control (Liu et al. 1986, Benson 1986, Hu´etink et al. 1990, Ghosh and Kikuchi 1991, Huerta and Casadei 1994, Rodr´ıguez-Ferran et al. 1998). The basic idea of the ALE formulation is the use of a referential domain for the description of motion, different from the material domain (Lagrangian description) and the spatial domain (Eulerian description). When compared to fluid dynamics, where the ALE formulation originated (see Donea 1983 and references therein), the main difficulty of ALE solid mechanics is the path–dependent behaviour of plasticity models. The relative motion between the mesh and the material must be accounted for in the treatment of the constitutive equation. Two approaches may be used to describe large (elastic and inelastic) strains. In hypoelastic–plastic models (Bonet and Wood 1997, Simo and Hughes 1998, Belytschko et al. 2000), the evolution of stresses is expressed in rate format, relating an objective stress rate with a rate of deformation. In hyperelastic–plastic models (Bonet and Wood 1997, Simo and Hughes 1998, Simo 1998, Belytschko et al. 2000), on the contrary, there is no rate equation for stresses: they can be computed from the deformation gradient 44 ALE formulations mesh (i.e. the mapping of the finite element mesh in the referential domain Rχinto the spatial domain Rx) must be ensured by the ALE remeshing strategy. This is standard in ALE analysis of fluid dynamics and hypoelastoplasticity, and is in sharp contrast with the situation in Yamada and Kikuchi (1993), where both the spatial mesh and the material mesh (i.e. the mapping of the FE mesh in Rχinto RX) must be kept undistorted, thus seriously limiting the potential of the ALE description as an r–adaptive technique. Moreover, the need of convection cannot be regarded as a significant drawback of the proposed approach, since convecting variables is needed anyway in the ALE description of the plastic response. Hyperelastoplasticity Indeed, if plastic strains are considered, a convective term appears in the equation that describes the evolution of plastic variables. For instance, Armero and Love (2000) work with the strain measure Gp=(FpTFp)−1, which is related to the elastic left Cauchy-Green tensor bethrough be=FG pFT, and rewrite equation (2.2.10)1 into ˙ Gp=−2˙γ(F−1mτF)Gp(2.2.27) in a Lagrangian setting, or ∂Gp ∂t |χ−∇χGp·FΨ −1∂X ∂t |χ=−2˙γ(F−1mτF)Gp(2.2.28) in an ALE setting. In equation (2.2.28), |χdenotes that the referential time derivative is computed holding the grid point χfixed, and ∇χis the gradient operator with respect to referential coordinates. As expected, the material time derivative in the left-hand-side of equation (2.2.27) is transformed, in equation (2.2.28), into a referential time derivative and a convective term which accounts for the relative motion between material particles Xand grid points χ. A similar result is obtained for equation (2.2.10)2. Alternatively, equation (2.2.10) can be reformulated into an ALE setting as ∂be ∂t |χ+c∇xbe=lbe+belT−2˙γmτ(τ,q)be(2.2.29) ∂p ∂t |χ+c∇xp=˙γmq(τ,q) (2.2.30) In equations (2.2.29) and (2.2.30), the material time derivative, see equations (2.2.10) and (2.2.11), has again been replaced by a referential time derivative and a convective term. The main difference between this approach and the one represented by equation (2.2.28) resides in the convective term: the relative motion is now expressed by the so– called convective velocity c(i.e. the difference between the particle velocity and the mesh velocity), and the gradient operator is now with respect to spatial, not referential, coordinates. In fact, equations (2.2.29) and (2.2.30) are in the “quasi–Eulerian” format ALE formulations 45 FOR EVERY TIME–STEP [nt, n+1t]: Material phase •Neglect convective terms •Advance the solution in an updated Lagrangian fashion: compute the increment of particle displacements n+1∆uand quantities Lbe,Lpand det(LF) (superscript L denotes Lagrangian) Remeshing •Compute the increment of mesh displacements n+1∆uΦand the increment of convective displacements n+1∆uconv by means of a remeshing algorithm that reduces element distortion •Compute the convective velocity c=n+1∆uconv/∆t Convection phase •Account for convective terms •Use the Godunov-type technique to convect quantities Lbe, Lpand det(LF)inton+1be,n+1pand det(n+1F) •Compute stresses n+1τand n+1σ Table 2.6: The proposed ALE approach for hyperelastoplasticity commonly encountered in ALE fluid dynamics and ALE hypoelastoplasticity (Huerta and Liu 1988). Similarly, if the general case of non–isochoric response is considered det(F)mustalso be convected, because it cannot be computed solely from be. For doing so, it is convenient to rewrite equation (2.2.8) into ∂|F| ∂t |χ+c∇x|F|=|F|∇x·v.(2.2.31) Note that equation (2.2.31) has the same structure that equations (2.2.29) and (2.2.30): in the left–hand–side, a referential time derivative and a convective term; in the right– hand–side, the material terms. In summary, the quantities to convect in the proposed ALE approach are be,pand det(F). This is done by means of a fractional–step method, a very common strategy to treat ALE convective terms (Hu´etink et al. 1990, Baaijens 1993, Huerta and Casadei 1994, Rodr´ıguez-Ferran et al. 1998, Askes et al. 1998, Armero and Love 2000). Every time–step is divided into two phases: the Lagrangian phase and the convection phase. During the Lagrangian phase, convection is neglected and the increment of particle 46 ALE formulations displacements n+1∆uis computed in the usual Lagrangian fashion (i.e. elastic predictor and plastic corrector). After that, an ALE remeshing algorithm is employed to compute the increment of mesh displacements n+1∆uΦ. During the convection phase, the convective term is taken into account. A Godunov–like technique (Huerta, Casadei and Donea 1995, Rodr´ıguez-Ferran et al. 1998) is used for that purpose. The proposed ALE approach for hyperelastoplasticity is summarized in table 2.6. Remark 2.2.4.In the proposed approach, the time–integration of the elastic response is not exact. During the convection phase, truncation errors are introduced. It must be noted, however, that these numerical errors can be controlled by the time–step ∆tand the mesh size. Moreover, these errors are completely unrelated with the drawbacks of hypoelastic–plastic models (i.e. elastic dissipation in a closed path). From the viewpoint of modelling, the proposed approach is fully hyperelastic–plastic. ✷ 2.2.4 Numerical examples The proposed ALE approach is illustrated and validated here by means of three representative numerical examples: a necking test, a coining test and a powder compaction test. Eight–noded quadrilateral elements with 2 ×2 Gauss points are employed for all the computations. Necking test The necking problem is a well–known benchmark test in large–strain solid mechanics (Simo 1988, Rodr´ıguez-Ferran et al. 1998, Miehe 1998, Peri´c and de Souza Neto 1999). A cylindrical bar with circular cross–section, with a radius of 6.413 mm and 53.334 mm length, is subjected to uniaxial extension. A slight geometric imperfection (1% reduction in radius), see figure 2.23, induces necking in the central part of the bar. An axisymmetric analysis is carried out with the mesh of 5 ×10 finite elements shown in figure 2.23. Two hyperelastic–plastic models have been used: the classical von Mises model (Simo and Hughes 1998, Simo 1998) and a Tresca–type model (Miehe 1998). The two models are isochoric and exhibit hardening. In consequence, the quantities to transport in the convection phase are beand p, which contains one internal plastic variable. The material parameters for both models are summarized in table 2.7, see Simo (1988) and P´erez-Foguet and Armero (2000) for further details. For comparative purposes, both Lagrangian and ALE analyses have been performed. A very simple ALE remeshing strategy is used (Rodr´ıguez-Ferran et al. 1998): the outer region BDEG of the mesh is Lagrangian, and equal height of elements is prescribed in the central region ABGH. The results with the von Mises model are discussed first. Figure 2.24 depicts the deformed shapes up to an elongation dof 8 mm for half the bar. As expected, the elements ALE formulations 47 A B C D E F H G AD = EH = 26.667 mm AB = GH = 0.167 AD CD = EF = 0.5 AD DE = FC = 6.413 mm HA = 0.99 DE . . . Figure 2.23: Necking of a cylindrical bar. Problem definition and computational mesh. in the neck zone become very distorted with the Lagrangian description, see figure 2.24(a– c). The distortion is highly reduced with the ALE description, see figure 2.24(d–f). In the proposed ALE approach, the quality of these spatial meshes is the only concern of the ALE remeshing strategy. There is no need to ensure the quality of the material meshes (Yamada and Kikuchi 1993, Armero and Love 2000). A more quantitative comparison is offered in figure 2.25, which shows the evolution of the vertical reaction and dimensionless radius (ratio of current radius to initial radius) with elongation. Very similar results are obtained up to an elongation of 6.5–7 mm. If pulling proceeds, however, discrepancies arise between the Lagrangian and ALE solutions. With only one row of (very distorted) elements in the necking zone, the Lagrangian simulation on the coarse mesh does not fully capture the plastification process, and this results in less necking. The ALE response is in much better agreement with a reference Lagrangian solution with a very fine mesh (not shown in the figure). Bulk modulus 164.206 GPa Shear modulus 80.1938 GPa Initial flow stress 0.45 GPa Residual flow stress 0.715 GPa Linear hardening coefficient 0.12924 GPa Saturation exponent 16.93 Yield shape exponent (Tresca–type model) 20 Table 2.7: Material parameters for the necking test. 48 ALE formulations Lagrangian formulation ALE formulation d = 7 mm d = 8 mm (a) (b) (c) d = 6 mm d = 7 mm d = 8 mm (e)(d) (f) d = 6 mm Figure 2.24: Necking test with the von Mises model. Mesh configurations for different top displacements d: (a-c) Lagrangian formulation and (d-f) ALE formulation. (a) (b) 0 10 20 30 40 50 60 70 80 012345678 Edge displacement [mm] Edge reaction [kN] 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 012345678 Edge displacement [mm] Dimesionless radius ALE Lag. ALE Lag. Figure 2.25: Necking test with the von Mises model. Lagrangian and ALE formulations. Global response: (a) vertical edge reaction and (b) dimensionless radius in the necking zone versus vertical edge displacement. ALE formulations 49 0.04 0.11 0.18 0.25 0.32 0.39 0.46 0.53 0.60 0.67 0.74 0.81 0.88 0.95 Hyperelastic model Hypoelastic model Lagrangian formulation ALE formulation Lagrangian formulation ALE formulation (c) (a) (b) (d) [GPa] q Figure 2.26: Necking test with the von Mises model. Distribution of the von Mises stress in the necking zone. Hyperelastic–plastic model: (a) Lagrangian formulation, (b) ALE formulation. Hypoelastic–plastic model (after Rodriguez–Ferran et al. 1998): (c) Lagrangian formulation, (d) ALE formulation. 50 ALE formulations (a) (b) 0 10 20 30 40 50 60 70 80 012345678 Edge displacement [mm] Edge reaction [kN] 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 012345678 Edge displacement [mm] Dimesionless radius ALE Lag. ALE Lag. Figure 2.27: Necking test with the Tresca–type model. Lagrangian and ALE formulations. Global response: (a) vertical edge reaction and (b) dimensionless radius in the necking zone versus vertical edge displacement. 0.04 0.11 0.18 0.25 0.32 0.39 0.46 0.53 0.60 0.67 0.74 0.81 0.88 0.95 Lagrangian formulation ALE formulation (a) (b) [GPa] q Figure 2.28: Necking test with the Tresca–type model. Distribution of the von Mises stress in the necking zone: (a) Lagrangian formulation, (b) ALE formulation. ALE formulations 51 Lagrangian formulation ALE formulation d = 7 mm d = 8 mm (a) (b) (c) d = 6 mm d = 7 mm d = 8 mm (e)(d) (f) d = 6 mm Figure 2.29: Necking test with the Tresca–type model. Mesh configurations for different top displacements, d: (a-c) Lagrangian formulation and (d-f) ALE formulation. A final, qualitative assessment of the two simulations can be made by looking at the distribution of the von Mises stress. According to some empirical and semianalytical studies (Goicolea 1985) this field is constant in the neck zone, along the z= 0 plane of symmetry. Figure 2.26 shows the distribution of the von Mises stress after an elongation of 7 mm. It can be seen that the Lagrangian analysis, figure 2.26(a), does not provide a constant value along z= 0, while the ALE analysis does, figure 2.26(b). The results reported in Rodr´ıguez-Ferran et al. (1998) for the hypoelastic–plastic von Mises model are shown in figures 2.26(c) and 2.26(d). Since elastic strains are small compared to elastic strains for this test, the hyperelastic and hypoelastic approaches yield very similar results. The same test is performed with the Tresca–type model. Figure 2.29 shows the deformed shapes as elongation proceeds. The highly distorted Lagrangian meshes show the different failure pattern with respect to the von Mises model, compare figures 2.24(a–c) and 2.29(a–c). The ALE meshes, however, are very similar, compare figures 2.24(d–f) and 2.29(d–f), because the same ALE remeshing strategy is used. Thanks to the ALE description, the quality of the spatial mesh can be ensured independently of the material deformation. Regarding the structural response and the von Mises stress distribution, results are qualitatively very similar to those with the von Mises model. Figure 2.27 shows that the global behaviour of the piece is captured correctly with the Lagrangian formulation only up to an elongation of 6.5–7 mm; figure 2.28 shows that only the ALE formulation correctly captures the constant stress distribution in the neck zone. 52 ALE formulations ABC D E AC = DF = 30 mm CD = FA = 10 mm AB = EF = 12 mm F Figure 2.30: Coining test. Problem definition and computational mesh. Coining test As a second example, a coining process is simulated (Rodr´ıguez-Ferran et al. 1998). A metallic disk, with a radius of 30 mm and a height of 10 mm, is deformed by a punch 12 mm in radius, see figure 2.30. The von Mises hyperelastic–plastic model is employed, with elastic modulus E= 200 GPa, Poisson’s coefficient ν=0.3, yield stress σy= 250 MPa and plastic modulus Ep= 1 GPa. Both the punch and the die are rigid. Perfect friction (stick) conditions are assumed in the punch–disk and disk–die interfaces. An axisymmetric analysis is performed to model a 60% height reduction with a mesh of 20×8 finite elements. Again, both a Lagrangian and an ALE analysis have been performed. The evolution of the deformed shape and the von Mises stress field is depicted in figure 2.31. Due to the stick conditions in the two interfaces, the material in the central part of the disk flows outward. This leads to a highly distorted Lagrangian mesh, especially under the punch corner and in the contact with the die. Remarkably, the convergence of the Lagrangian analysis is not disturbed by the mesh distortion, thanks to the use of consistent tangent matrices and a 2 ×2 Gauss–point quadrature, with negative Jacobians less likely than with a 3 ×3 quadrature. The mesh distortion is greatly reduced in the ALE analysis. The following remeshing strategy is used: (1) line AB is Eulerian, lines CD and EF are Lagrangian and equal length of elements is prescribed in lines BC, CF and DA; (2) parabolic profiles of horizontal mesh displacements are prescribed in region ABEF; (3) mesh displacements in region BCDE are interpolated from the contour values. As expected, the mesh distortion clearly affects the numerical solution. In figure 2.32, the Lagrangian and ALE solutions are compared in terms of the displacement of points C and D and the punch reaction (radial and vertical components shown separately). For all the outputs, the two descriptions provide very similar results up to a height reduction of 20%–30%, which induces a relatively small distortion in the Lagrangian mesh. Significant differences, however, are encountered for increasing height reductions and a more distorted Lagrangian mesh. The same type of behaviour is obtained with the Tresca–type model (with the same material parameters used for the von Mises analysis plus the yield shape exponent equal to 20). ALE formulations 53 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95 1.05 1.15 1.25 1.35 1.45 1.55 1.65 1.75 Lagrangian formulation ALE formulation 12% 36% 24% 48% 60% [GPa] q Figure 2.31: Coining test with the von Mises model. Mesh configurations and distribution of the von Mises stress for different height reductions. Lagrangian and ALE formulations. 60 ALE formulations Chapter 3 Non-standard consistent tangent operators Two topics related with the computation of consistent tangent operators in elastoplasticity are developed in this chapter. Both are of special interest for the solution of problems involving highly nonlinear constitutive models. First, in section 3.1, numerical differentiation is applied to integrate the plastic constitutive laws and to compute the corresponding consistent tangent operators. In many cases, the plastic constitutive laws used in geomaterial modelling have complicated expressions of the flow vector and the internal variable flow direction. The derivatives of both are needed to achieve quadratic convergence in the numerical time–integration of the constitutive equations at Gauss–point level and in the solution of boundary value problems with an incremental–iterative strategy. In this context, it is shown that the numerical approximation of first derivatives is a simple, robust and competitive alternative to analytical derivatives. The behaviour of several numerical differentiation techniques is analyzed. On the other hand, in section 3.2, a very simple expression of consistent tangent operators for substepping time–integration schemes is presented. Substepping schemes based on the forward Euler, the backward Euler and the midpoint rules are included. An adaptive strategy which activates a substepping scheme at Gauss points where the standard return mapping algorithms does not converge is used as an example. With this strategy, the most restrictive local problems no longer control the load increment of the global problem. This results in a large reduction of the number of load increments of the global problem, which typically implies a large reduction of the computational cost. Examples involving the adaptive substepping scheme jointly with the proposed numerical differentiation approximations are included. It is shown that the use of both strategies lead to quadratic convergence in complex nonlinear elastoplastic problems. 61 62 Consistent tangent operators 3.1 Numerical differentiation for local and global tangent operators in computational plasticity In this section, numerical differentiation is applied to integrate plastic constitutive laws and to compute the corresponding consistent tangent operators. The derivatives of the constitutive equations are approximated by means of difference schemes. These derivatives are needed to achieve quadratic convergence in the integration at Gauss–point level and in the solution of the boundary value problem. Numerical differentiation is shown to be a simple, robust and competitive alternative to analytical derivatives. Quadratic convergence is maintained, provided that adequate schemes and stepsizes are chosen. This point is illustrated by means of some numerical examples. 3.1.1 Introduction The main goal of this section is to show that numerical differentiation is a useful tool for achieving quadratic convergence in computational plasticity. Two different problems must be solved: 1) the integration of the constitutive law at each Gauss point (the local problem) and 2) the boundary value problem (the global problem); see, for instance, Crisfield (1991, 1997). If the constitutive law is integrated with an implicit method, the local problem is nonlinear. To solve it with the Newton-Raphson method, thus attaining quadratic convergence (Dennis and Schnabel 1983), it is necessary to compute the Jacobian of the residual at the Gauss–point level. The global problem, on the other hand, is typically solved via an incremental/iterative approach (Crisfield 1991). At each load increment, a nonlinear system of equations must be solved. To do it with the Newton-Raphson method and achieve quadratic convergence, the consistent tangent matrix (computed with the consistent elastoplastic moduli) must be used (Simo and Taylor 1985, Runesson et al. 1986). In both problems (local and global), the derivatives of the constitutive equation are needed. These derivatives are a key ingredient of both the Jacobian of the residual —local problem— and the consistent elastoplastic moduli —global problem. Various approaches for computing these derivatives can be found in the literature. For simple plasticity models, analytical derivatives are readily available, and this leads to closed–form return mapping algorithms for the local problem and compact, explicit expressions of the consistent elastoplastic moduli for the global problem (Simo and Taylor 1985, Simo and Hughes 1998). In more complicated models, analytical differentiation is rather more cumbersome. Algebraic manipulators such as Maple or Mathematica can be a very effective tool for obtaining analytical derivatives. Here a different approach is proposed: derivatives are approximated by means of classical difference schemes (Isaacson and Keller 1966, Dennis and Schnabel 1983, Hoffman Consistent tangent operators 63 1982). The approximated derivatives are used both for the integration of the constitutive equations (local problem) and for the computation of the consistent tangent moduli (global problem). The resulting algorithm is both robust and computationally efficient. It maintains the characteristic quadratic convergence of the Newton-Raphson method, provided that adequate difference schemes and stepsizes (i.e. the perturbation in the difference scheme) are chosen. Some applications of numerical differentiation to other problems of computational plasticity can be found in the literature. For instance, a first–order forward difference scheme is used by Miehe (1996) to compute the consistent tangent moduli needed in large–strain inelasticity. There, simple material models are used and closed expressions for the return mapping are available, so numerical derivatives are only applied to the global problem. The proposed numerical approach computes the derivatives of each stress component with respect to each strain component to get directly the consistent moduli. This implies a computational overhead (with respect to analytical derivatives) that ranges from 40% to 80% and, perhaps more importantly, its robustness is limited by the choice of the stepsize. More complicated models are considered by Jeremi´c and Sture (1994, 1997), but it is also concluded that analytical derivatives are clearly superior to numerical derivatives regarding computational cost. In constrast to these approaches, the strategy proposed here, concerned with small strains, combines the three following features: 1) it can be used for both local and global problems, with complicated material models, see section 4.2 (P´erez-Foguet, Rodr´ıguez-Ferran and Huerta 2000e), 2) it has a marginal computational overhead (around 1 — 2 %) and 3) it is very robust (i.e. insensitive to the choice of the stepsize). An outline of this section follows. The local and the global problems of computational plasticity are briefly reviewed in subsection 3.1.2, with emphasis on the crucial role of the derivatives of the constitutive equation. In subsection 3.1.3 the various numerical differentiation techniques are summarized, and a simple rule for selecting the stepsize is presented. These numerical approximations are applied to several problems in subsection 3.1.4. Finally, some concluding remarks are made on subsection 3.1.5. For the sake of simplicity, the simple case of a single yield surface and the backward Euler integration rule is considered. However, the same derivatives are needed in multisurface plasticity (Simo et al. 1988, Pramono and Willam 1989b) or with other implicit integration rules (Ortiz and Popov 1985, Runesson, Sture and Willam 1988, Chaboche and Cailletaud 1996), so numerical differentiation can also be applied in these general cases. 3.1.2 Problem statement In small strain elastoplasticity, the stress tensor σand the stress–like internal variables qcan be related with the small strain tensor εand the strain–like internal variables p 64 Consistent tangent operators (Lubliner 1990, Simo and Hughes 1998) through σ=∇We(εe)=∇We(ε−εp)andq=−∇Wp(p),(3.1.1) where εeand εpare the elastic and the plastic strain tensors, Weand Wpare the elastic and the plastic part of the free–energy function per unit of volume, and ∇Ψ(χ) denotes the gradient of Ψ with respect to χ. The yield function Fdefines the admissible stress states, F(σ,q)≤0. Its derivatives are denoted by nσ=∂F/∂σand nq=∂F/∂q. The equations of evolution for εpand pare ˙ εp=˙ λmσand ˙ p=˙ λmq,(3.1.2) where mσand mqare the corresponding flow directions, and ˙ λis the plastic consistency parameter. For convenience, the terms generalized stress and generalized flow vector are used to refer to (σ,q)and(mσ,mq) respectively (Ortiz and Martin 1989, Simo and Hughes 1998). For some constitutive models, the generalized flow vector is expressed as the derivative of a generalized flow potential G(σ,q) (Lubliner 1990, Runesson and Larsson 1993): mσ=∂G/∂σand mq=∂G/∂q. The flow potential Gmay either coincide with (associate plasticity) or differ from (non-associate plasticity) the yield function F. Local problem: integration of the constitutive law The integration of the constitutive law with an implicit rule results in a nonlinear problem. If the backward Euler rule is used, the equations are n+1σ=∇We(n+1ε−nεp−n+1λn+1mσ) (3.1.3) n+1q=−∇Wp(np+n+1λn+1mq) (3.1.4) F(n+1σ,n+1q) = 0 (3.1.5) where the superscripts nand n+ 1 refer to instants ntand n+1t=nt+∆trespectively. There are several ways to solve this nonlinear problem. For some very simple models, it can be solved analytically (Crisfield 1991, Simo and Hughes 1998). However, for general models an analytical solution is not possible, so a numerical method must be used. A typical choice is the Newton-Raphson method, because it converges quadratically (Dennis and Schnabel 1983). To use this method, the derivatives of equations (3.1.3–3.1.5) with respect to the unknowns n+1σ,n+1q,andn+1λare needed. With the standard vector notation of computational mechanics (Zienkiewicz and Taylor 1988, Crisfield 1991), this Jacobian is n+1J=    Inσ+λE∂mσ ∂σλE∂mσ ∂qEmσ λH∂mq ∂σInq+λH∂mq ∂qHmq nT σnT q0     t=n+1t ,(3.1.6) Consistent tangent operators 65 where E=∇2We(εe) is the tensor of elastic moduli, H=∇2Wp(p) is the tensor of plastic moduli, Tdenotes transpose, nσis the dimension of the stress vector, nqis the number of internal variables and I∗is the identity matrix of order ∗. The expression of the Jacobian given in equation (3.1.6) can be compacted by multiplying the first row by E−1and the second row by H−1(Simo and Hughes 1998). This transformation (which must, of course, also be performed on the RHS residual vector of the Newton-Raphson iteration) is computationally appealing if matrices Eand Hare constant. However, the ideas presented in this section do not rely at all in this transformation. For this reason, the original expression of the Jacobian is retained here. Of the terms in equation (3.1.6), the derivatives of the generalized flow vector with respect to the generalized stresses are usually the most difficult to compute. For complex material models, these derivatives are either not available or computationally too expensive (Pramono and Willam 1989b, Sture et al. 1989, Etse and Willam 1994, Jeremi´c and Sture 1997). The standard approach in these cases is to use nonlinear solvers (different from the Newton-Raphson method, i.e. a fully tangent approach) that do not need to compute all the derivatives of the generalized flow vector. Some of these alternatives are: a tangent approach for the plastic multiplier, equation (3.1.5), with an explicit expression for internal variables, equation (3.1.4), and a secant approach for stresses, equation (3.1.3) (Pramono and Willam 1989b); a tangent approach for the stresses, equation (3.1.3), and a direct substitution of the internal variable equations, equation (3.1.4) (Jeremi´c and Sture 1994, 1997); a two–level technique with a tangent approach for the stress invariants and a Picard iteration with an adaptive order inverse interpolation for the internal variables (Macari et al. 1997). However, quadratic convergence cannot be achieved with these methods, because they are not based on a consistent linearization of all the equations and unknowns. The goal of this section is to show that numerical differentiation is a valid alternative to these approaches. By approximating numerically the derivatives of the generalized flow vector, the standard full Newton-Raphson method can be applied to equations (3.1.3– 3.1.5). In this manner, quadratic convergence is obtained. Global problem: the consistent tangent matrix Various nonlinear solvers may be used for the global problem (Dennis and Schnabel 1983, Crisfield 1991). Again, one of the best choices is the full Newton-Raphson method. To achieve quadratic convergence, the consistent tangent matrix must be used (Simo and Taylor 1985, Runesson et al. 1986). To compute this matrix, the consistent tangent moduli dn+1σ/dn+1εare needed. They can be computed by linearizing the discrete constitutive 66 Consistent tangent operators equations (3.1.3–3.1.5). This linearization can be represented as n+1J     dn+1σ dn+1q dn+1λ      =     Edn+1ε 0 0      (3.1.7) where n+1Jis the Jacobian of the local problem defined in equation (3.1.6). The upper– left block–matrix of the inverse of n+1Jcontains the consistent tangent moduli. In the literature, compact expressions (after inverting n+1Jand taking the upper–left block– matrix) of these moduli can be found for particular models (Ortiz and Martin 1989, Crisfield 1991, Simo and Hughes 1998). However, the more general expression given in equation (3.1.7) is preferred here because it highlights an important fact in the context of this work: the derivatives of the generalized flow vector needed to solve the local problem are also required in the computation of the consistent tangent matrix. Numerical differentiation is used by Miehe (1996) to approximate directly the consistent tangent moduli dn+1σ/dn+1ε(or, more precisely, an equivalent expression for the large– strain problems treated there). The resulting algorithm has a considerable CPU overhead in comparison to analytical derivation (40% to 80%) because the derivatives of all the stress components with respect to all the strain components are approximated numerically. 3.1.3 Numerical differentiation The derivatives of the generalized flow vector with respect to the generalized stresses are approximated by means of classical difference schemes. The approximation will be used for both the local and the global problems defined in the previous section. Thus, it must be accurate enough to maintain the characteristic quadratic convergence of the Newton-Raphson method. Some authors (Jeremi´c and Sture 1994, 1997) suggest to approximate numerically the second derivatives of the flow potential G(recall that the flow vector is the derivative of the flow potential). The standard approach to obtain second order of accuracy is the typical centered difference scheme (Isaacson and Keller 1966, Dennis and Schnabel 1983, Hoffman 1982) applied to a general n–dimensional function, f(x): ∂2f ∂x2 i =f(x+hiei)−2f(x)+f(x−hiei) h2 i +O(h2 i).(3.1.8) In equation (3.1.8), xiis the ith component of x,eithe ith unit vector, hithe stepsize in the ith direction and the Odenotes the order of convergence. The scheme represented by equation (3.1.8) will be denoted by 2ND-O(h2), see table 3.1. The approach used in this section consists on approximating numerically the first derivatives of (the analytical expression of) the flow vector. That is, the flow vector can be obtained via analytical differentiation of the flow potential (this step is relatively Consistent tangent operators 67 Notation Description 1ND-O(h)Forward difference scheme for first derivatives of the flow vector 1ND-O(h2) Centered difference scheme for first derivatives of the flow vector 1CND-O(h2) Approximation to first derivatives of the flow vector using complex variable 2ND-O(h2) Centered difference scheme for second derivatives of the flow potential Table 3.1: Numerical approximations to the derivatives of the flow vector. simple, even for complex constitutive laws) or it can be an input of the model. Then, numerical differentiation is applied to approximate the derivatives of the flow vector (which is the computationally involved step for complex models). Standard forward or centered difference schemes are used: ∂f ∂xi =f(x+hiei)−f(x) hi +O(hi) (3.1.9) ∂f ∂xi =f(x+hiei)−f(x−hiei) 2hi +O(h2 i).(3.1.10) The schemes represented by equations (3.1.9) and (3.1.10) will be denoted by 1ND-O(h) and 1ND-O(h2) respectively, see table 3.1. It must be noted that the generic function fplays the role of a component of the flow vector in equations (3.1.9) and (3.1.10), whereas it denotes the flow potential in equation (3.1.8). Error analysis The key issue in numerical differentiation is the choice of the stepsize hi. Approximations based on difference schemes are affected by truncation and rounding errors (Stepleman and Winarsky 1979, Dennis and Schnabel 1983). The truncation errors (represented in equations (3.1.8–3.1.10) with the Osymbol) decrease as the stepsize tends to zero. The rounding errors, on the contrary, increase as the stepsize tends to zero. Therefore, there is an optimal stepsize hopt that minimizes the summation of both errors. Dennis and Schnabel (1983) present an expression of this optimal stepsize for first and second–order approximations to first derivatives (1ND-O(h)and1ND-O(h2)) and for first–order approximation to second derivatives (not used in this section). Their work has been extended here to the case of second–order approximation to second derivatives (2ND-O(h2)). The optimal stepsize hopt can be written as hopt =hopt rmax{|x|,typx},(3.1.11) 68 Consistent tangent operators where hopt ris the optimal relative stepsize and typxis a typical value of xused to avoid choosing a null (or extremely small) hopt for null (or extremely small) x. Numerical experimentation shows that typxcan be chosen in a rather arbitrary manner, because it has a very small influence on the results (in all the numerical examples of section 4, typx= 1). The main idea behind equation (3.1.11) is that hopt ris independent of x. This means that a constant value of hopt rcan be used all over the domain, for every load step, and all the stress components. In general, it is not possible to compute the exact value of hopt r. This would require the rigorous minimization of the sum of truncation and rounding errors. However, the following expressions can be found after some simplifying assumptions: 1ND −O(h): hopt r=√rf(3.1.12) 1ND −O(h2): hopt r=3 √rf(3.1.13) 2ND −O(h2): hopt r=4 √rf.(3.1.14) In equations (3.1.12–3.1.14), rfis the accuracy in the evaluation of f. The expressions (3.1.12) and (3.1.13) are given by Dennis and Schnabel (1983). Expression (3.1.14) has been derived following the same arguments, as described next. Consider equation (3.1.8). A standard error propagation analysis renders the following bound on the rounding error ERo in the approximation to the second derivative, i.e. errors induced by finite precision computations of the first term in the RHS of equation (3.1.8): |ERo|≤(4rf+12r)ˆ f h2.(3.1.15) In equation (3.1.15), ris the machine precision (r≈10−16 in IEEE double precision), ˆ fis an upper bound of |f|in the neighborhood of xand the subscript iis ommited from the stepsize hto ease the notation. Regarding the truncation error ETr , i.e. the second term in the RHS of equation (3.1.8), a bound can be easily derived (Isaacson and Keller 1966, Hoffman 1982) as |ETr|≤γh2 12 ,(3.1.16) where γis a bound on the fourth derivative of f. Putting equations (3.1.15) and (3.1.16) together, the following bound on the total error ETot is obtained: |ETot |≤|ETr |+|ERo|≤γh2 12 +(4rf+12r)ˆ f h2.(3.1.17) The effect of the truncation error and the rounding error on the total error is clearly highlighted in a double logarithmic graph of ETot versus h: the rounding error contributes with a line of slope 2 (because second derivatives are approximated) and the truncation error with a line of slope −2 (because a second–order scheme is used). Similar results are Consistent tangent operators 69 obtained for the other schemes, see Dennis and Schnabel (1983) and figures 3.1 and 3.3 of subsection 3.1.4. The optimal value of the stepsize, which minimizes the bound of |ETot |given by equation (3.1.17) is hopt =4 12(4rf+12r)ˆ f γ.(3.1.18) The value of hopt cannot be computed from equation (3.1.18): the values of ˆ fand γ depend on x, and, in general, they are very difficult to approximate. Following Dennis and Schnabel (1983), it is assumed that ˆ f γ∝(max{|x|,typx})4,(3.1.19) where the fourth power is obtained as the sum of the degree of derivation and the order of the approximation (2 and 2), and ∝means “is proportional to”. This assumption is rather reasonable; it is verified, for instance, for monomials f(x)=xα. Combining equations (3.1.18) and (3.1.19) yields hopt ∝4 max{rf,r}max{|x|,typx}.(3.1.20) Two further assumptions are needed: 1) rfis larger or equal to r(if function fis simple, then rfis similar to r, but if many operations are involved, rfcan be quite larger than r) and 2) the proportionality constant can be taken as 1 (again a reasonable assumption, as illustrated by monomials f(x)=xα). With these assumptions, one finally gets the expression of the relative optimal stepsize given in equation (3.1.14). The concept of relative stepsize is essential in computational plasticity. In a global problem the range of values of the generalized stresses is usually very large, and the values of the different components at a certain Gauss point are very different. For this reason, the optimal stepsize hopt can show huge variations. However, the previous analysis shows that the optimal relative stepsize hopt rcan be assumed to be constant. Since a simple technique is wanted, the same value of hopt r(obtained from equations 3.1.12, 3.1.13 or 3.1.14, depending on the scheme) will be used for all the computation. In fact, equations (3.1.12–3.1.14) are only used to select the order of magnitude of hopt r, because there is a wide range of relative stepsizes for which the difference schemes (3.1.8–3.1.10) are accurate enough to attain quadratic convergence. In this context, the assumptions needed to deduce the simplified expressions of hopt r, equations (3.1.12–3.1.14), do not restrict at all the applicability of the proposed approach. This point is illustrated in the next section by means of several numerical examples. Remark 3.1.1.An unconventional approximation of first derivatives has also been used in this work. It is based on the theory of functions of complex variable (Lyness and Moler 1967, Squire and Trapp 1998). If f:C→Cis analytic in a neighborhood that 76 Consistent tangent operators Figure 3.5: Perforated strip under traction (after Simo 1985). Due to symmetry, only one quarter is considered. Application to the global problem Numerical differentiation is applied to solve several boundary value problems (i.e., global problems). That is, the numerical approximations of table 3.1 are employed to compute consistent tangent matrices. Moreover, in the examples with the RHMC model, they are also used to solve the local problem. For all the stress components and over the whole domain a constant relative stepsize has been used. Four examples are presented: two with von Mises plasticity and two more with the RHMC model. These examples will illustrate that any numerical approximation to first derivatives of the flow vector (1ND-O(h), 1ND-O(h2) and 1CND-O(h2)) is useful to solve the global problem with quadratic convergence. On the contrary, approximating numerically the second derivatives of the flow potential (2ND-O(h2)) is not robust enough to achieve quadratic convergence. Von Mises model with perfect plasticity First, a perforated strip under uniaxial traction is analyzed with von Mises perfect plasticity. This test is presented by Simo and Taylor (1985), and it is depicted in figure 3.5. A total displacement of 0.2misimposed in ten steps. The convergence results for the fourth and the eighth load steps obtained with the analytical consistent tangent matrix and with numerical differentiation are compared in figure 3.6. The ranges of relative stepsizes that give the same convergence results that analytical consistent tangent matrix (up to a relative error in energy of 10−12) are presented in table 3.2 (convergence results are considered to be the same if the energy error has the same order of magnitude during all the iterations of all the load steps). Also in table 3.2, Consistent tangent operators 77 between parenthesis, the values of hrthat give almost the same results as the analytical derivatives are indicated (convergence results are considered to be almost the same if the energy error has the same order of magnitude during all the iterations except the last one of each load step). It can be seen that the difference schemes for first derivatives, 1ND-O(h) and 1ND-O(h2), give the same result that analytical differentiation for a wide range of relative stepsizes, even though a very strict tolerance has been used. As expected, second– order of accuracy presents a wider range of adequate relative stepsizes than first–order. The 1CND-O(h2) approximation presents the same convergence results as the 1ND-O(h2) one for relative stepsizes higher than the optimal one. Moreover, in agreement with the previous results (see figure 3.1), it maintains the quadratic convergence for arbitrarily small relative stepsizes. On the other hand, the numerical difference scheme for second derivatives of the flow potential, 2ND-O(h2), does not give very good results even in this simple boundary value problem (the typical quadratic convergence is lost). Von Mises model with exponential isotropic hardening The perforated strip under traction has also been simulated using the von Mises model with exponential isotropic hardening. The problem definition is the same of Simo and Taylor (1985), except for the plastic parameters: the initial yield stress is equal to 0.243 MPa, the yield stress at infinite equivalent plastic strain is 0.729 MPa and the exponential parameter is 0.1 MPa. In figure 3.7, the convergence results for the fourth and the eighth load steps are depicted, and in table 3.3 the ranges of relative stepsizes that give the same convergence results that analytical consistent tangent matrices are presented. The results are more strict than those obtained with perfect plasticity: numerical first derivatives present a narrower range of relative stepsizes that give the same results as analytical derivatives, and the numerical second derivatives do not attain quadratic convergence. Moreover, note that the optimal relative stepsize is higher than the one of figure 3.1. This is in agreement with the previous comments about the value of hopt rwhen the components of the stress tensor have very different values. On the other hand, the 1CND-O(h2) approximation presents the same behavior as in perfect plasticity: the results are the same as the 1ND-O(h2) approximation for hrhigher than the optimal one, and quadratic convergence is achieved for arbitrarily small hr. The main conclusion of the two examples with von Mises plasticity is that numerical first derivatives of the flow vector do not affect the properties of convergence of the NewtonRaphson method (i.e., the global consistent tangent matrix is accurately approximated). The typical difference schemes present an adequate behavior for a wide range of relative stepsizes. And the unconventional approximation based on complex variable theory allows to use stepsizes as small as wanted. On the other hand, it has been shown that the numerical second derivatives of flow potential are not robust enough: they do not work properly in demanding boundary value problems. 78 Consistent tangent operators 1ND-O(h) 4th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 12345678 Iteration Relative error Analytical hr = 1.E-2 hr = 1.E-7 hr = 1.E-11 1ND-O(h2) 4th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 12345678 Iteration Relative error Analytical hr = 1.E-2 hr = 1.E-7 hr = 1.E-11 2ND-O(h2) 4th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 12345678 Iteration Relative error Analytical hr = 1.E-1 hr = 1.E-2 hr = 1.E-4 1ND-O(h) 8th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1234567 Iteration Relative error Analytical hr = 1.E-2 hr = 1.E-7 hr = 1.E-11 1ND-O(h2) 8th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1234567 Iteration Relative error Analytical hr = 1.E-2 hr = 1.E-7 hr = 1.E-11 2ND-O(h2) 8th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1234567 Iteration Relative error Analytical hr = 1.E-1 hr = 1.E-2 hr = 1.E-4 Figure 3.6: Convergence results for the fourth and the eighth load steps. Von Mises perfect plasticity. Consistent tangent operators 79 1ND-O(h) 4th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 123456789 Iteration Relative error Analytical hr = 1.E-2 hr = 1.E-7 hr = 1.E-11 1ND-O(h2) 4th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 123456789 Iteration Relative error Analytical hr = 1.E-2 hr = 1.E-7 hr = 1.E-11 2ND-O(h2) 4th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 123456789 Iteration Relative error Analytical hr = 1.E-1 hr = 1.E-2 hr = 1.E-4 1ND-O(h) 8th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 123456 Iteration Relative error Analytical hr = 1.E-2 hr = 1.E-7 hr = 1.E-11 1ND-O(h2) 8th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 123456 Iteration Relative error Analytical hr = 1.E-2 hr = 1.E-7 hr = 1.E-11 2ND-O(h2) 8th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 123456 Iteration Relative error Analytical hr = 1.E-1 hr = 1.E-2 Figure 3.7: Convergence results for the fourth and the eighth load steps. Von Mises exponential hardening plasticity. 80 Consistent tangent operators Num. approx. Range of hr 1ND-O(h)(10−4)—10 −5—10 −7—(10 −8) 1ND-O(h2)(10−2)—10 −3—10 −8—(10 −9) 1CND-O(h2)(10−2)—10 −3—... 2ND-O(h2)(10−2) Table 3.2: Relative stepsizes that give the same convergence results as analytical derivatives, for the von Mises perfect plasticity global problem. Num. approx. Range of hr 1ND-O(h)(10−4)—10 −5—(10 −6) 1ND-O(h2) 10−3—10 −5—(10 −6)—(10 −7) 1CND-O(h2) 10−3—... 2ND-O(h2) — Table 3.3: Relative stepsizes that give the same convergence results as analytical derivatives, for the von Mises exponential hardening global problem. Figure 3.8: Pile problem (after Potts and Gens 1985). In the following, the applicability of the approximations 1ND-O(h), 1ND-O(h2)and 2ND-O(h2) to the Rounded Hyperbolic Mohr-Coulomb (RHMC) model is assessed by means of two boundary value problems. In both problems numerical differentiation is applied to the local and the global problem simultaneously. RHMC model: vertical displacement of a pile The first problem is the vertical displacement of a pile. The definition of the problem is presented by Potts and Gens (1985), and it is only summarized here. Figure 3.8 shows the finite element mesh. It corresponds to a horizontal disc of soil. The thickness of the disc is 5 units of length (u.l.) and the pile radius is 7.5 u.l. To model the loading of the pile, a vertical displacement of 2 u.l. is imposed over the boundary AF in 20 load steps. To model the infinite extension of the disc, zero vertical displacement is imposed over the boundary CD. Due to the essentially one–dimensional nature of the problem, vertical lines (such as EB) are prescribed to remain vertical during loading. Table 3.4 shows the convergence results for the tenth load step, and in table 3.5 the ranges of relative stepsizes that give the same convergence results that analytical consistent tangent matrices are presented. It is clear that the numerical approximations to the second derivatives of the flow potential, 2ND-O(h2), do not yield quadratic convergence. Consistent tangent operators 81 (a) Analytic 10−310−510−710−9 12.2 E+00 2.2 E+00 2.2 E+00 2.2 E+00 2.2 E+00 21.5 E-03 1.5 E-03 1.5 E-03 1.5 E-03 1.5 E-03 33.1 E-05 3.1 E-05 3.1 E-05 3.1 E-05 3.1 E-05 48.6 E-09 1.3 E-07 8.6 E-09 8.6 E-09 8.6 E-09 54.6 E-15 1.7 E-09 4.7 E-15 1.4 E-13 8.0 E-12 62.6 E-11 8.7 E-12 73.4 E-13 1.1 E-11 8... (b) Analytic 10−210−3 12.2 E+00 2.2 E+00 2.2 E+00 21.5 E-03 1.5 E-03 1.5 E-03 33.1 E-05 4.0 E-05 3.0 E-05 48.6 E-09 7.2 E-06 1.8 E-07 54.6 E-15 8.6 E-07 1.7 E-08 61.3 E-07 1.5 E-09 71.7 E-08 1.4 E-10 82.4 E-09 1.3 E-11 93.2 E-10 1.2 E-12 10 4.5 E-11 2.1 E-13 11 ... Table 3.4: The convergence results for sixth load step of the pile problem: (a) 1ND-O(h), (b) 2ND-O(h2). To reach quadratic convergence it would be necessary to do a supplementary effort to choose a correct stepsize at each Gauss point. The results of the approximation 1NDO(h) are quite similar to the previous example, see table 3.3. The hrthat gives quadratic convergence in the global problem is higher than the estimated value from an analysis of the approximations to the flow vector derivatives (see figure 3.3). This is in agreement with the previous analysis of the influence of the stresses on the value hopt r. On the other hand, although one may conclude from table 3.5 that the range of adequate hris quite small, it must be pointed out that the comparison has been done up to a very strict tolerance of 10−12. If the results are compared to a tolerance of 10−8(that for a relative energy error is still a strict tolerance), the range becomes 10−4to 10−9. Thus, the 1ND-O(h) approximation is accurate enough. The results of 1ND-O(h2) are also similar to previous ones: the range of adequate hris larger than for 1ND-O(h), see table 3.5. Nevertheless, if realistic tolerances are considered (i.e., 10−8or higher), the advantage of second order approximation is lost: the two alternatives give good results for a very wide range of relative stepsizes. Moreover, since the computational cost (the number of evaluations of the flow vector) of the first– order approximation is half of the second–order one, it can be concluded that the first– order scheme is more suitable. 82 Consistent tangent operators Num. approx. Range of hr 1ND-O(h)(10−4)—10 −5—(10 −6) 1ND-O(h2) 10−4—10 −6—(10 −7) 2ND-O(h2) — Table 3.5: Relative stepsizes that give the same convergence results as analytical derivatives, for the pile problem. (a) (b) Smooth Smooth Smooth 605 four-noded elements B/2 5B Smooth Figure 3.9: Rigid footing: (a) problem definition, (b) mesh and final plastic strains. Due to symmetry, only one half is considered. RHMC model: vertical displacement of a rigid footing The second global problem solved using the RHMC model is the vertical displacement of a rigid footing (Abbo and Sloan 1996). The scheme of the problem is depicted in figure 3.9(a), and the mesh and the final distribution of plastic strains are shown in figure 3.9(b). The dimension of the rigid footing, B, is 20 units of length (u.l.) and a vertical displacement of 0.02 u.l. is imposed in 20 increments. The convergence results for the load steps fifteen and twenty are shown in figure 3.10, and the ranges of relative stepsizes that give the same convergence results that analytical consistent tangent matrices are presented in table 3.6. The results are similar to previous ones: the 1ND-O(h) approximation gives quadratic convergence with a small range of relative stepsizes if a very strict tolerance is used (10−12), and with second–order of accuracy, 1ND-O(h2), the range of relative stepsizes is wider. However, as in the pile problem example, with a tolerance of 10−8almost any relative stepsize gives good results for both approximations. Therefore, for practical applications the first–order difference Consistent tangent operators 83 1ND-O(h) 15th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 12345678 Iteration Relative error Analitycal hr = 1.E-3 hr = 1.E-5 hr = 1.E-7 hr = 1.E-9 1ND-O(h) 20th Load step 1,E-17 1,E-13 1,E-09 1,E-05 1,E-01 1234567 Relative error Analitycal hr = 1.E-3 hr = 1.E-5 hr = 1.E-7 hr = 1.E-9 1ND-O(h2) 15th Load step 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 12345678 Iteration Relative error Analitycal hr = 1.E-2 hr = 1.E-4 hr = 1.E-6 hr = 1.E-8 1ND-O(h2) 20th Load step 1,E-17 1,E-13 1,E-09 1,E-05 1,E-01 1234567 Iteration Relative error Analitycal hr = 1.E-2 hr = 1.E-4 hr = 1.E-6 hr = 1.E-8 Iteration Figure 3.10: Convergence results for the load steps fifteen and twenty. Rigid footing with RHMC model. Num. approx. Range of hr 1ND-O(h)(10−4)—10 −5—(10 −6)—(10 −7) 1ND-O(h2)(10−2)—10 −3—10 −5—(10 −6)—(10 −8) Table 3.6: Relative stepsizes that give the same convergence results as analytical derivatives, for the rigid footing problem. scheme approximation, 1ND-O(h), is accurate enough. Moreover, the choice of the relative stepsize is not a problem (as in the previous example, the values of hrbetween 10−4and 10−9give quadratic convergence up to a relative error less than 10−8). In this global problem, the computing time of the alternatives 1ND-O(h) and 1NDO(h2) has been analyzed. In table 3.7 there are the time overheads with respect to the problems solved with analytical consistent tangent matrices. All the approximations, except 1ND-O(h)withhr=10 −3, have needed the same number of iterations. For this particular case, the 1ND-O(h) approximation is between 1 — 1.5% more expensive than the one computed with the analytical derivatives, and the 1ND-O(h2) is between 2.6 — 2.8%. Therefore, the computational cost of the proposed approach is marginal, even in the context of material models (von Mises and RHMC) where analytical derivatives are 84 Consistent tangent operators hr1ND-O(h)1ND-O(h2) 10−34.1% 2.8% 10−42.5% 3.4% 10−51.5% 2.6% 10−61.0% 2.8% 10−71.4% 2.6% 10−81.2% 3.8% Table 3.7: Time overheads of the numerical approximations, for the rigid footing problem. available and relatively simple to compute. The main conclusion of these last two examples is that typical difference schemes applied to first derivatives of the flow vector are adequate to approximate the consistent tangent matrix and the Jacobian of the residual of the integration rule for highly nonlinear flow potentials. The convergence properties are the same that with analytical derivatives for a wide range of relative stepsizes. The first–order approximation is enough to maintain quadratic convergence up to a tolerance of less than 10−8. Nevertheless, for tolerances less or equal to 10−12, second–order accuracy is needed in order to have a wide range of relative stepsizes. Finally, for complex models the numerical second derivatives of the flow potential are not robust enough to maintain quadratic convergence, even in simple boundary value problems. 3.1.5 Concluding remarks It has been shown that numerical differentiation is a useful tool in computational plasticity, both in the integration of the constitutive law (local problem) and in the computation of the consistent tangent matrix (global problem). The analytical derivatives of the generalized flow vector with respect to the generalized stresses are the components most difficult to compute in the consistent tangent matrix and the Jacobian of the local residual. In some cases, they are not even available. The main conclusion of this section is that numerical differentiation is a valid alternative to analytical derivatives. Two approaches are possible: 1) approximating the second derivatives of the flow potential (recall that the flow vector is the first derivative of the flow potential) or 2) approximating the first derivatives of the (analytical) flow vector. The first approach, suggested in the literature, is not robust enough. It can be used to solve the local problem, but not the global problem. The resulting consistent tangent matrices are not accurate enough and the Newton-Raphson method looses its characteristic quadratic convergence. The second approach, on the contrary, is a simple and robust alternative to analytical differentiation. Quadratic convergence is achieved, both for simple (von Mises) and more complicated (Rounded Hyperbolic Mohr-Coulomb) material models. Consistent tangent operators 85 Various schemes may be used to approximate the first derivatives of the flow vector: the classical first–order and second-order difference schemes, and an unconventional secondorder approximation based on complex variable theory. The first–order approximation maintains quadratic convergence up to a tolerance of less than 10−8, and the second–order approximations up to less than 10−12. Thus, all of them are accurate enough for any practical application. The choice of an adequate stepsize (a typical problem of difference schemes) does not present any difficulty. The concept of the relative stepsize has been presented, and it has been verified that the three approximations to first derivatives of the flow vector give good results with a wide range of relative stepsizes. That is, the proposed strategy is very robust in the sense that the choice of the stepsize has a very small influence on its performance. Moreover, the complex variable approximation presents a very interesting property: the rounding errors are very small and constant (i.e., they do not increase as the stepsize tends to zero), so arbitrarily small stepsizes may be used. The computational overhead of the proposed strategy (with repect to analytical derivatives) is marginal, even for material models where the analytical derivatives have relatively simple expressions. This result is in sharp contrast with previous applications of numerical differentiation to computational plasticity (Jeremi´c and Sture 1994, 1997, Miehe 1996). In section 4.2 (P´erez-Foguet et al. 2000e), numerical differentiation is applied to more complex constitutive models, where analytical derivatives are not available. 92 Consistent tangent operators which can be rewritten as d( n+k mσ,n+k mκ,λ k) d∆ε=n+k mAαkPE +d( n+k−1 mσ,n+k−1 mκ,λ k−1) d∆ε(3.2.14) using PT=(Inσ,0nσ,nκ+1)and n+k mA=n+k mJ−1  Inσ+nκ0nσ+nκ,1 01,nσ+nκ0 .(3.2.15) Equation (3.2.14) is valid for k=m,...,2. For k= 1, the first substep, equation (3.2.8) is linearized as follows: d( n+1 mσ,n+1 mκ,λ 1) d∆ε=α1n+1 mJ−1     E 0nκ,nσ 01,nσ      =α1n+1 mAPE .(3.2.16) The final expression is obtained after substitution of equation (3.2.16) into equation (3.2.14) particularized at k= 2, and a recursive use of (3.2.14) from k= 3 up to k=m. Finally, d( n+1σ,n+1κ,λ m) d∆ε=n+1AαmPE +n+m−1 mAαm−1PE +··· ···α2PE +α1n+1 mAPE ···   m−1 .(3.2.17) The consistent tangent moduli are in the leading principal minor of order nσof the LHS of equation (3.2.17). They are obtained by means of the projection matrix P, and have the compact expression dn+1σ d∆ε=PT  m  i=1  αi i  j=m n+j mA  PE .(3.2.18) Equation (3.2.18) has the same structure of equation (3.2.6). The summation in brackets of equation (3.2.18) plays the same role that the inverse of the Jacobian of equation (3.2.6). In fact, equation (3.2.18) is identical to equation (3.2.6) if only one substep is considered, m=1. Note that the computation of the consistent tangent moduli involves matrix inversions, even in the case of single–step integration rules with no substepping, equation (3.2.6). It must be remarked, however, that small matrices are considered, so numerical inversion poses no difficulties. The order of the Jacobian matrices to be inverted is only nσ+nκ+1, see equations (3.2.4) and (3.2.11), and it can be further reduced. A computationally efficient expression of the consistent tangent moduli, with smaller matrices and a recursive structure, is presented in subsection 3.2.7. Consistent tangent operators 93 < 3 3 5 6 4 5 < 3 36 >1 0 536 >1 0>1 0 5 I1/3=−2ccotφI1/3=0 I 1/3=2ccotφ θ −30◦30◦ θ −30◦30◦ θ −30◦30◦ J2 0 10 Figure 3.11: Rounded Hyperbolic Mohr-Coulomb model. Number of iterations of the local problem for trial stresses in the deviatoric plane (θ–J2) at various levels of confinement (I1/3). 3.2.4 Examples In this subsection, two global problems are solved using the consistent tangent matrix for substepping presented in subsection 3.2.3. Quadratic convergence is attained. The first example is the simulation of a rigid footing with the Rounded Hyperbolic Mohr-Coulomb model. When solving this problem with single–step integration rules, the Gauss points located under the corner of the footing restrict the load increment to a forbidding small value. When substepping is used, on the contrary, the load increment no longer depends on local demands. As a consequence, larger load increments can be used. Until now, the use of substepping was incompatible with quadratic convergence at the global level. However, thanks to the consistent tangent matrix presented in subsection 3.2.3, quadratic convergence results are presented in the following. The second example is the simulation of a triaxial test with the MRS–Lade model. In this example substepping is combined with numerical differentiation. For the MRS–Lade model, not all the derivatives needed to compute the Jacobian of the residual, equation (3.2.4), have a readily available analytical expression. Because of this, and following sections 3.1 and 4.2, numerical differentiation is used to approximate the Jacobian. In these references, numerical differentiation is combined with the backward Euler integration rule and quadratic convergence is obtained at local and global level. In this subsection, numerical differentiation is combined with substepping and quadratic results are also obtained. Simulation of a rigid footing with the Rounded Hyperbolic Mohr-Coulomb model In the following a rigid footing is analysed with the Rounded Hyperbolic Mohr-Coulomb model and substepping (Abbo and Sloan 1995, 1996). The consistent tangent matrix developed in subsection 3.2.3 is employed to achieve quadratic convergence at the global 94 Consistent tangent operators (a) (b) Smooth Smooth Smooth 1784 eight-noded elements B/2 10 B Smooth Figure 3.12: Rigid footing: (a) problem definition, (b) mesh. level. The model is presented by Abbo and Sloan (1995). The flow potential–yield surface is divided into two regions: one corresponding to the Hyperbolic Mohr-Coulomb zone that smoothes the apex of the classical Mohr-Coulomb model on the hydrostatic axis, and the other corresponding to the rounding zones, that smoothes the corners present on the deviatoric plane. The traces of the flow potential–yield surface on the meridian and the deviatoric planes are presented in figure 3.2. The dimensionless material parameters are a cohesion of 1, a friction angle of 30◦, a Young modulus of 3000 and a Poisson coefficient of 0.3. Associated plasticity is considered. First, the convergence of the local problem is analysed. Figure 3.11 depicts the number of iterations for convergence (up to a tolerance of 10−12) with the standard initial approximation, equation (3.2.5). Note that the Newton-Raphson method does not converge in 10 iterations in some regions of the stress space. This is due to the high curvature of the yield function close to the apex (I1/3≥0) combined with the non-differentiability of the flow vector at |θ|=25 ◦. In order to avoid these regions of non-convergence of the local problem, the adaptive substepping scheme is activated at the Gauss points that require more than 6 iterations at the local level for convergence (to a tolerance of 10−12). Due to symmetry, only one half of domain is considered for the global problem, see figure 3.12. The soil mass is modelled as a square of 20 units of length (u.l.), twenty times the footing half–width, B = 2 u.l. It has been checked that this domain is large enough to preclude any undesired influence of the boundary on the results. An unstructured mesh of 1784 quadrilateral eight–noded elements is used. A vertical displacement of the footing Consistent tangent operators 95 0 10 20 30 40 0,000 0,025 0,050 0,075 0,100 Vertical displacement [u.l.] Dimensionless force Figure 3.13: Rigid footing problem. Dimensionless force versus vertical displacement. of 0.1 u.l. is prescribed in 100, 200, 400 and 800 uniform increments. The relationship between force and vertical displacement is depicted in figure 3.13. The computed limit dimensionless force is 30.71, 2% above the exact Prandtl collapse dimensionless force of 30.14. Figure 3.14 shows the distribution of equivalent plastic strain for different values of the load level (i.e. fraction of the total prescribed displacement). Note that the failure mechanism is well captured. Very similar results are obtained for the four problems, with 100, 200, 400 and 800 load increments (l.i.). The substepping has been activated just under the right corner of the footing. This agrees with the fact that the non-convergence regions of the local problem are close to the rounded apex. Figure 3.15 shows the evolution of the number of Gauss points with substepping activated, the sum of substeps over all the domain and the maximum number of substeps for all the Gauss points, for the problem solved with 100, 200, 400 and 800 l.i. The number of Gauss points with substepping activated is not very different in the four problems (it is reduced by a factor of less than 2 when the number of load increments increases by a factor of 8). On the other hand, the total number of substeps and the maximum number of substeps are in inverse proportion to the number of steps; if the number of load increments is doubled, the maximum number of substeps is divided by two. This indicates that a very large number of steps would be needed to solve the problem without substepping (extrapolating the results of figure 3.15, the number of uniform l.i. would be greater than 50 000). The convergence results for several load levels and for the four problems (with 100, 200, 400 and 800 l.i.) are shown in figure 3.16. All the results are quadratic. As expected, the number of iterations per load increment decreases as the number of load increments increases. In fact, the problem with 100 l.i. requires up to 11 Newton-Raphson iterations at the increments previous to the plateau in the load–displacement curve. This indicates that larger increments should not be used in this part of the problem. It has been checked that the influence of the substepping criterion is marginal. If the threshold for activating the substepping is set at 12 iterations (instead of 6), the same results are found (except 96 Consistent tangent operators 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 Load level 0.25 Load level 0.50 Load level 0.75 Load level 1.00 Figure 3.14: Rigid footing problem. Equivalent plastic strain for different load levels. for the sum of substeps over the domain, which is a little lower). The computational cost of the four load discretizations is compared in table 3.8 (relative CPU time) and in figure 3.17 (accumulated iterations). The computational cost increases with the number of load increments. Within the range presented in table 3.8, twice the number of load increments implies a computational cost 1.6 times greater. The case with no substepping (i.e. 50 000 uniform increments) is also shown in figure 3.17. Note that the computational cost is much higher: the number of accumulated iterations exceeds 10 000 after only one–eighth of the analysis. This clearly illustrates the computational efficiency of the substepping scheme with consistent tangent matrix. Load increments 100 200 400 800 Relative CPU time 100% 160% 263% 432% Table 3.8: Rigid footing problem. Relationship between number of global load increments and relative CPU time. Consistent tangent operators 97 0 10 20 30 40 0,00 0,25 0,50 0,75 1,00 Load level Num ber of Gauss points 100 load incr. 200 load incr. 400 load incr. 800 load incr. 0 50 0 1000 1500 2000 2500 0,00 0,25 0,50 0,75 1,00 Load level Sum of substeps over the dom ain 100 load incr. 200 load incr. 400 load incr. 800 load incr. 0 100 200 30 0 400 50 0 600 0,00 0,25 0,50 0,75 1,00 Load level Max. Substeps at Gauss point 100 load incr. 200 load incr. 400 load incr. 800 load incr. Figure 3.15: Rigid footing problem. Evolution of the number of Gauss points with substepping activated (top), the sum of substeps over the domain (center) and the maximum number of substeps for all the Gauss points (bottom). 98 Consistent tangent operators ORDGLQFUHPHQWV ORDGLQFUHPHQWV    ORDG  LQFUHPHQWV    ORDG  LQFUHPHQWV 1,E-17 1,E-14 1,E-11 1,E-08 1,E-05 1,E-02 1234567891011 ,WHUDWLRQ 5 H O D W L Y H  H U U R U Load level 0.2 Load level 0.5 Load level 0.6 Load level 0.7 Load level 0.8 1,E-17 1,E-14 1,E-11 1,E-08 1,E-05 1,E-02 1234567891011 ,WHUDWLRQ 5 H O D W L Y H  H U U R U Load level 0.2 Load level 0.5 Load level 0.6 Load level 0.7 Load level 0.8 1,E-17 1,E-14 1,E-11 1,E-08 1,E-05 1,E-02 1234567891011 ,WHUDWLRQ 5 H O D W L Y H  H U U R U Load level 0.2 Load level 0.5 Load level 0.6 Load level 0.7 Load level 0.8 1,E-17 1,E-14 1,E-11 1,E-08 1,E-05 1,E-02 1234567891011 ,WHUDWLRQ 5 H O D W L Y H  H U U R U Load level 0.2 Load level 0.5 Load level 0.6 Load level 0.7 Load level 0.8 Figure 3.16: Rigid footing problem. Convergence results for different load levels. Triaxial test with the MRS–Lade model In the following, a triaxial test with end–platen friction is analysed with the MRS–Lade model (Sture et al. 1989, P´erez-Foguet and Huerta 1999, P´erez-Foguet et al. 2000e)and substepping. The MRS–Lade model is used to simulate the behaviour of granular materials under both low and high confinement stresses (Macari et al. 1994, Macari et al. 1997). It features 1) a two–surface yield function, comprising a smooth cone surface and a smooth cap surface, 2) hardening and softening variables that depend on dissipated plastic work, and 3) a non-associated flow rule in the meridian plane of the cone region. Several slight modifications to the original formulation of the model have been devised (Jeremi´cand Sture 1994, Macari et al. 1997, P´erez-Foguet and Huerta 1999). In this section the modification presented in section 4.1 (P´erez-Foguet and Huerta 1999) is used. It consists on a new definition of the flow vector for the cone region that avoids the corner problem or the flip–over of previous formulations. The traces of the yield surface on the meridian plane and on the deviatoric plane are Consistent tangent operators 99 0 2000 4000 6000 8000 10000 0,00 0,25 0,50 0,75 1,00 Load level Accum ulated iterations 100 load incr. 200 load incr. 400 load incr. 800 load incr. W ithout Subs. Figure 3.17: Rigid footing problem. Relationship between accumulated iterations and load level. depicted in figure 4.4, and the hardening–softening function ηcon(κcon) of the cone in figure 4.5. The value of ηcon is directly related with the friction angle (slope of the cone in the meridian plane) and κcon is the cone internal variable, which depends on the plastic work. The softening at Gauss point level starts for κcon =1. The MRS–Lade model exhibits a high coupling between the flow vector and the plastic moduli. This coupling makes the analytical computation of the derivatives a very cumbersome task. However, using any of the numerical differentiation techniques presented in sections 3.1 and 4.2, all the derivatives are computed in a simple and efficient way. Going beyond sections 3.1 and 4.2, substepping and numerical differentiation are combined here. Like in the first example, the adaptive substepping technique is activated at the Gauss points where the local problem requires more than 6 iterations for convergence (up to a tolerance of 10−12). As shown in section 4.2, quadratic convergence at the local level is obtained in all the stress–internal variable space, even with large total strain increments and no substepping. Substepping is used for this problem to ensure proper time–integration of the constitutive law and a reduction of the computational cost, not to avoid non-convergence regions as with the Rounded Hyperbolic Mohr-Coulomb model. A structured mesh of 1350 (30×45) quadrilateral eight–noded elements has been used. Due to double symmetry only the upper right quarter of the sample is modelled. The end– platen friction is imposed by restraining the radial displacement of the sample top during loading. The same material parameters used by Macari et al. (1997) to simulate the triaxial test at local level of a Sacramento River sand are used in this example. Since the MRS– Lade model is not regularized and includes non-associated plasticity and softening, the problem can localize. However, the low degree of softening of the material parameters used and the axisymmetric nature of the test prevent localization (Rudnicki and Rice 1975). 100 Consistent tangent operators 0 2500 5000 7500 10000 0% 2% 4% 6% 8% 10% Vertical displacement / initial height Force [kN] 10 load incr. 100 load incr. Figure 3.18: Triaxial test problem. Force versus relative vertical displacement. The results do not depend significantly on the space or time discretization. The problem has been solved with 10, 20, 50 and 100 load increments. The curves of force versus relative vertical displacement (i.e. vertical displacement over initial height) for 10 and 100 l.i. are depicted in figure 3.18. The results are almost identical for all load discretizations (the relative error of the force at the end of the simulation computed with 10 l.i. is less than 0.6%). Figure 3.19 shows the evolution of the distribution of the cone internal variable. As expected, the material response is clearly non-homogeneous. Note that a wide region at the top of the sample does not enter in the softening regime (κcon <1), even for large vertical displacements, while in the center of the sample (lower left corner of the computational domain) softening starts before the global limit force is reached. Figure 3.20 shows the distribution of the number of substeps required at various load levels. Figure 3.21 shows the evolution of the number of Gauss points with substepping activated, the sum of substeps over the domain and the maximum number of substeps for the four problems (with 10, 20, 50 and 100 l.i.). In the problems with 10 and 20 l.i., substepping is activated from the beginning of the analysis and almost everywhere in the domain. At the end of the test, part of the domain changes from plastic loading to elastic unloading. In that region the substepping is deactivated. This results in a decrease in the number of Gauss points with substepping activated, see figures 3.21(top) and 3.20(a). However, the number of substeps at each Gauss point is very low, compared with the previous example. On the other hand, with 50 l.i., substepping is activated only in the region of the domain where the local problems are more demanding. Finally, note that with 100 l.i. there are only a few Gauss points with substepping activated at the end of the test, and they require only two substeps. Load increments 10 20 50 100 Relative CPU time 119% 100% 120% 215% Table 3.9: Triaxial test problem. Relationship between number of global load increments and relative CPU time. Consistent tangent operators 101 0. 1. 2. 3. 4. 2.5% 5% 7.5% 10% Figure 3.19: Triaxial test problem. Distribution of the cone internal variable for different values of the relative vertical displacement. 1 2 3 4 6 7 8 9 10 5 2.5% 5% 7.5% 10% (a) 1 2 3 4 6 7 8 9 10 5 2.5% 5% 7.5% 10% (b) Figure 3.20: Triaxial test problem. Distribution of the number of substeps required for different values of the relative vertical displacement. Problem solved with (a) 20 l.i., and (b) with 50 l.i. 108 Consistent tangent operators First, the equivalence      Inσ0nσ,nκ0nσ,1 0nκ,nσInκ0nκ,1 01,nσ01,nκ0      =     Inσ0nσ,nκ 0nκ,nσInκ 01,nσ01,nκ       Inσ0nσ,nκ0nσ,1 0nκ,nσInκ0nκ,1 (3.2.34) is employed in order to rewrite equation (3.2.13) into d( n+k mσ,n+k mκ,λ k) d∆ε=n+k mJ−1     Inσ0nσ,nκ 0nκ,nσInκ 01,nσ01,nκ         Inσ0nσ,nκ0nσ,1 0nκ,nσInκ0nκ,1 d( n+k−1 mσ,n+k−1 mκ,λ k−1) d∆ε+αk  E 0nκ,nσ  .(3.2.35) Then, equation (3.2.35) is pre–multiplied by   Inσ0nσ,nκ0nσ,1 0nκ,nσInκ0nκ,1 (3.2.36) in order to get d( n+k mσ,n+k mκ) d∆ε=n+k mAcd( n+k−1 mσ,n+k−1 mκ) d∆ε+αkPcE(3.2.37) where Pc=  Inσ 0nκ,nσ , n+k mAc=  Inσ0nσ,nκ0nσ,1 0nκ,nσInκ0nκ,1 n+k mJ−1     Inσ0nσ,nκ 0nκ,nσInκ 01,nσ01,nκ      , (3.2.38) and use is made of the relation   Inσ0nσ,nκ0nσ,1 0nκ,nσInκ0nκ,1 d( n+k mσ,n+k mκ,λ k) d∆ε=d( n+k mσ,n+k mκ) d∆ε.(3.2.39) Finally, following the same process of section 3.2.3, the consistent tangent moduli are obtained. The final expression is dn+1σ d∆ε=PcT  m  i=1  αi i  j=m n+j mAc  PcE.(3.2.40) Consistent tangent operators 109 Equation (3.2.40) has the same structure of equation (3.2.18). However, the matrices involved here are smaller. In fact, the consistent tangent moduli are computed with the following expression: dn+1σ d∆ε=PcTn+1AcαmPc+n+m−1 mAcαm−1Pc+··· ···α2Pc+α1n+1 mAcPc···   m−1 E,(3.2.41) which is equivalent to equation (3.2.40). As suggested by equation (3.2.41), the consistent tangent moduli are computed recursively during time–integration: when the substep jis integrated, the matrix n+j mAcis computed and the corresponding part of the consistent tangent moduli is calculated. The process is always the same, except for the first and the last substeps. Moreover, because of the special structure of the Jacobian, the matrices n+j mAcare computed inverting only the leading principal minor of order nσ+nκof the Jacobian. This result can be obtained using Sherman and Morrison’s lemma (Dennis and Mor´e 1977). 110 Consistent tangent operators Chapter 4 Elastoplastic models for granular materials In this chapter, the highly nonlinear elastoplastic behaviour of granular materials is modelled following two different approaches. In both cases the attention is focused in the efficient solution of the nonlinear constitutive equations and the boundary value problems by means of the corresponding tangent operators. The first part, sections 4.1 to 4.3, is devoted to the analysis of a work hardening cone–cap model for small strain problems. In section 4.1 a new formulation of the plastic potential of the model is presented. This formulation avoids the grey zone and the flip– over found at cone–cap intersection of previous ones. A detailed analysis of the model features and its application to various boundary value problems is presented in section 4.2. The consistent linearization of all equations with respect to all unknowns and the application of the numerical differentiation techniques presented in section 3.1 lead to quadratic convergence results in both the time–integration of constitutive laws and the solution of boundary value problems. Previous applications found in the literature do not exhibit quadratic convergence. Finally, in section 4.3, an extension of the model which includes the modelling of cohesive behaviour is discussed. In the second part, section 4.4, several examples involving different density–dependent models within the framework of isotropic finite strain multiplicative plasticity are presented. The flow directions and yield functions of these models are expressed in terms of the Kirchhoff stresses and the relative density. This type of models include the finite strain multiplicative plasticity based on the Cauchy stresses as a particular case. The consistent linearization of the constitutive equations, including the density influence, is presented. Moreover, it is shown that the standard numerical time–integration based on the exponential return mapping also apply, without any modification, to this type of models. 111 112 Elastoplastic models 4.1 Plastic flow potential for the cone region of the MRS–Lade model The original formulation of the MRS–Lade model, with non-associated flow rule on the meridian plane in the cone region, has a corner. In order to reduce the computational effort of corner solution algorithms, a modified plastic flow potential for the cone part is found in the literature. This modification may have a non-admissible flip over of the flow vector in the cone–cap intersection if the plastic flow potential is not correctly defined. Here a corrected plastic flow potential for the cone region is defined to obtain a continuous transition of the flow vector. 4.1.1 Introduction In computational plasticity, the definition of the plastic flow vector is more useful than the definition of the plastic flow potential. The flow vector is needed for the integration of the constitutive law and for the resolution of the global finite element problem, see Ortiz and Popov (1985), Simo and Taylor (1985), Runesson et al. (1988) and Crisfield (1991, 1997) among others. In fact, the flow potential is hardly ever employed, it is defined mainly for convenience (Lubliner 1990). Nevertheless the flow potential is useful in theoretical analysis (Kim and Lade 1988, Lubliner 1990, Lade 1994) and in the formal description of the model, for instance, Pramono and Willam (1989a), Etse and Willam (1994) and Khan and Huang (1995). In order to implement a non-associated flow rule three approaches are possible. The first one defines the flow rule (usually the flow potential) directly from experimental analysis, independently of other characteristics of model, see Lade and Duncam (1975), Nova and Wood (1979) and Lade and Kim (1988b) among others. The second one prescribes the flow vector modifying the normal to the yield function (the corresponding flow potential is obtained by integration), see for instance, Runesson (1987), Alawaji, Runesson and Sture (1992) and Larsson and Runesson (1996). The third one defines the flow potential as a direct modification of the yield function, Macari et al. (1994), Jeremi´c and Sture (1994), Macari et al. (1997); this must be done carefully in order to get the desired properties in the flow vector. 4.1.2 MRS–Lade model This work focuses in the non-associated flow rule for the cone region of the MRS–Lade model. This model has been developed at the University of Colorado by Macari-Pasqualino, Runesson and Sture (Sture et al. 1989, Jeremi´c and Sture 1994, Macari et al. 1997) and it is a further development of Lade’s three-invariant model for cohesionless soils (Lade and Duncam 1975, Lade 1989). Elastoplastic models 113 The model has been used to simulate the behavior of granular materials, such as sand, under both low and high confinement stresses (Macari et al. 1994, Macari et al. 1997). Quoting Jeremi´c and Sture (1997), the MRS–Lade model features: •a two-surface formulation, comprising a smooth cone surface and a smooth cap surface intersecting in plane curve (ellipse segment) in the deviatoric plane, •hardening and softening variables for both surfaces are based on dissipated plastic work, •a non-associated flow rule in the meridian plane and an associated flow rule in the deviatoric plane of the cone region, and an associated flow rule in the cap region, •ability to model cohesive strength and a curved meridian in the cone region. In order to center in the essential issue of this work, a simplified version of the original model is described here. Detailed discussion is presented by Sture et al. (1989) and Jeremi´c and Sture (1994). The following expressions define the yield function Fcon =qg(θ)−ηconep=0 Fcap =p−αpcap (1 −α)pcap 2 +qg(θ) ηconeαpcap 2 −1=0 , (4.1.1) where p,qand θare functions of the stress invariants, g(θ) is a function that defines the shape in the deviatoric plane, ηcone and pcap are the hardening–softening functions, and αis a parameter of the model. A scheme of the trace of the yield function in the p—q plane is depicted in figure 4.1. The flow rule is associated on the entire cap region and in the deviatoric plane of the cone region. Thus, in the cap region, the flow potential is equal to Fcap, and, in the cone region, it has the same dependence of qand θthat Fcon does. Following the original formulation of the model, Sture et al. (1989), the non-associated flow defined in correspondence with the expansive behavior in the cone region is represented by a plastic potential function of the form Gcone =qg(θ)−nηconep, (4.1.2) where nis a non-negative constant. Typical values for nare close to 0.1. This potential reduces the dilatancy of the associated flow rule in the cone region. Note that for nequal to zero incompressibility is enforced for values of pbetween zero and αpcap. The components of the flow vector on the meridian plane, m=(mp,m q), are mp=−nηcone mq=g(θ).(4.1.3) 114 Elastoplastic models (a) p q G rey region αpcap pcap (b ) p q F lip over αpcap pcap (c) p qSm ooth tra n s itio n αpcap pcap Figure 4.1: Trace of the simplified MRS–Lade yield function and characteristic flow vectors on the p−qplane. Original formulation (a), modified flow potential (b) and corrected flow potential (c). Elastoplastic models 115 4.1.3 Modified plastic flow potential For the usual case, ndifferent than zero, this flow rule has a grey region at the intersection of cone and cap surfaces, see figure 4.1(a). In this region there is not continuity of the flow vector. This implies that the Koiter’s rule must be applied: the direction of the plastic strain rate is defined as a linear combination of the cone and cap flow vectors. Therefore, a corner solution algorithm is needed (Simo et al. 1988, Pramono and Willam 1989b, Hofstetter et al. 1993, Jeremi´c and Sture 1994). Such algorithms are usually expensive from a computational point of view. Thus, in order to reduce the computational effort, an alternative definition of the flow potential for the MRS–Lade model is presented in Macari et al. (1994) and Macari et al. (1997). To avoid the grey zone, the flow vectors corresponding to cone and cap regions at the corner, i.e. p=αpcap, must have the same direction, see figure 4.1(c). Thus the p-component of the cone flow vector must be zero at the corner. The previously cited references propose the use of a pressure dependent n, n=−γp−αpcap p+αpcap ,(4.1.4) where γis a non-negative constant. If equation (4.1.4) is used in (4.1.3) and the expression of the flow potential is not necessary at all, the grey zone disappears and a continuous variation of the flow vector is obtained. However, if equation (4.1.4) is directly substituted in (4.1.2), as Macari et al. (1994) and Macari et al. (1997) seem to indicate, a modified flow potential is defined which induces the following components of the flow vector on the meridian plane, mp=γηcone2p−αpcap p+αpcap −(p−αpcap)p (p+αpcap)2mq=g(θ).(4.1.5) Then, mpvaries from a negative value, −γηcone,atp= 0 to a positive one, γηcone/2 at the corner, p=αpcap. This variation is illustrated in figure 4.2, where characteristic values of γ=0.125, ηcone =0.2andαpcap = 640 are used, Macari et al. (1997). Thus the proposed objective is not attained. Moreover, since mpat αpcap is strictly positive, there is a flip over of the vector flow at the corner, see figure 4.1(b). This situation is not desired and, in general, not admissible. 4.1.4 Corrected plastic flow potential In order to obtain the desired flow vector a different flow potential must be defined. With the desired expression of the flow vector a flow potential is obtained by integration. After substitution of equation (4.1.4) into (4.1.3) the components of the flow vector, m#= (m# p,m # q), become m# p=γηcone p−αpcap p+αpcap m# q=mq=g(θ).(4.1.6) 116 Elastoplastic models -0,025 -0,020 -0,015 -0,010 -0,005 0,000 0,005 0,010 0,015 0 200 400 600 p mp Figure 4.2: p-component of modified flow vector as function of p. As stated previously this flow vector induces the desired behavior. In figure 4.3, the variation of m# pwith respect to pis presented. In this case, the component m# premains negative for every pless then αpcap and reaches zero at this limit. Therefore, the grey region disappears, without any flip over of the flow vectors, and a specially designed corner algorithm is precluded. Note that other expressions for the evolution of m# preaching zero from below at p=αpcap can be used, if they do conform with experimental results. The plastic flow potential corresponding to the flow vector defined in (4.1.6) is obtained by integration: G# cone =qg(θ)+γηcone(p−2αpcap ln(p+αpcap)) ,(4.1.7) where the integration constants are taken equal to zero since the purpose of (4.1.7) is simply the definition of a potential for the plastic strain rate. Equation (4.1.7) represents a new plastic flow potential for the cone region of the MRS–Lade model. The corresponding flow vector components are defined in equation (4.1.6). This flow vector is justified from a physical view point in Macari et al. (1994). But, the flow potential presented here will not induce the undesired features of the one presented in Macari et al. (1994). 4.1.5 Concluding remarks Flow potentials are seldom employed in computational plasticity if the flow vector is known a priori. Nevertheless, if they are needed for theoretical or verification purposes they must agree with the desired behavior of the flow vectors. If the flow rule is modified by acting on the flow vector, the corresponding flow potential should be obtained by simple Elastoplastic models 117 -0,025 -0,020 -0,015 -0,010 -0,005 0,000 0,005 0,010 0,015 0 200 400 600 p mp   Figure 4.3: p-component of corrected flow vector as function of p. integration. Here, a new plastic flow potential for the cone region of the MRS–Lade model is presented. This potential induces a flow vector with continuous transition between cone and cap regions. Thus, the corner problem (grey region) inherent to the original formulation and the non-admissible flip over of previously published modifications, are avoided. 124 Elastoplastic models where nand ξare the derivatives of F(σ,κ) with respect to σand κrespectively, and nκ is the number internal variables (i.e. nκ= 2 for the MRS–Lade model). On the other hand, to solve the global problem with quadratic convergence it is necessary to use the consistent tangent matrix (Simo and Taylor 1985, Runesson et al. 1986). To compute this matrix, the consistent tangent moduli dn+1σ/dn+1εat each Gauss point are needed. They are obtained by linearizing equation (4.2.10). This linearization can be represented in a compact form as (Ortiz and Martin 1989) PTn+1J−1PE ,(4.2.12) where PT=(Inσ,0nκ+1) is the projection matrix on stress space. Therefore, the Jacobian matrix, equation (4.2.11), is needed for both the local and the global problems. The most difficult components to compute of the Jacobian are typically the derivatives of mand hwith respect to σand κ. This is the case for the MRS–Lade model, which exhibits a high coupling of all the components. Note, for instance, that the hardening moduli hare defined in terms of the flow vector m. This is caused by the fact that plastic work drives the hardening, see equation (4.2.7). Thus hdepends on σand κ both explicitly and through m, see equation (4.2.9). This coupling makes the analytical computation of the derivatives a very cumbersome task. Jeremi´c and Sture (1994), for instance, only present the analytical expression of some of the required derivatives for the original formulation of the model (that is, for Aconstant). They use these derivatives to solve the local problem and to compute an approximation to the consistent tangent matrix for the global problem. However, quadratic convergence is not achieved, because not all the required derivatives are used. 4.2.4 Numerical differentiation Indeed, quadratic convergence can only be attained by means of a full Newton-Raphson method. That is, all the derivatives of mand hwith respect to σand κare needed. One possibility would be to obtain the analytical expression of the missing derivatives. However, this is rather involved, even with the help of an algebraic manipulator. For this reason, a different course is followed here: all required derivatives are approximated numerically. Three of the techniques discussed in section 3.1 will be employed: the forward difference scheme, 1ND-O(h), the centered difference scheme, 1ND-O(h2), and the scheme based on complex variables, 1CND-O(h2), see table 4.1. The forward difference scheme is first–order accurate, and the other two schemes are second–order accurate. With these schemes, the derivative of miwith respect to κj(recall that vector notation is used), for Elastoplastic models 125 instance, is approximated either by 1ND −O(h)∂mi ∂κj (σ,κ)=mi(σ,κ+hej)−mi(σ,κ) h, 1ND −O(h2)∂mi ∂κj (σ,κ)=mi(σ,κ+hej)−mi(σ,κ−hej) 2h, 1CND −O(h2)∂mi ∂κj (σ,κ)=Immi(σ,κ+√−1hej) h, (4.2.13) where his the stepsize and ejis the jth unit vector. Similar expressions are used for ∂m/∂σ,∂h/∂σand ∂h/∂κ. The approximated derivatives are then used to solve the local and the global problems. Notation Description 1ND-O(h) Forward difference scheme (1st order accurate) 1ND-O(h2) Centered difference scheme (1st order accurate) 1CND-O(h2) Approximation based on complex variables (2nd order accurate) Table 4.1: Numerical approximations to first derivatives. A crucial issue in numerical differentiation is the choice of the stepsize h. In this work it is selected as shown in section 3.1, by using the concept of relative stepsize, hr. The optimal value of the relative stepsize can be approximated by √macheps for the first–order scheme, 1ND-O(h), and by 3 √macheps for the second–order accurate schemes, 1ND-O(h2) and 1CND-O(h2), with macheps the machine precision. Numerical experiments reveal a good behaviour of numerical differentiation (that is, quadratic convergence for both the local and the global problems) for a wide range of relative stepsizes, hr. In order to reduce the effect of rounding errors, hris taken as a negative power of 2 (hr=2 −k), not of 10 (hr=10 −k). This choice is relevant in some critical zones, as illustrated in next subsection. 4.2.5 Examples In this subsection, several local and global problems are solved quadratically with numerical differentiation. The three techniques of table 4.1 are compared and the main features of each one are remarked. Two sets of parameters have been used, see table 4.2. Soil S1 is a dense Sacramento River sand (Macari et al. 1997). Soil S2 is a small modification of soil S1. The modifications are 1) a smaller value of ¯ηcon, which reduces the size of the elastic domain (¯ηcon is the maximum value of ηcon, see equations (4.2.17) and (4.2.18) in subsection 4.2.7, and 2) different values of ccon and ;, which result in a more nonlinear evolution of the cone internal variable, κcon, see equations (4.2.9), (4.2.17) and (4.2.18). With these two modifications, soil S2 is quite more demanding from a numerical point of view than soil S1. 126 Elastoplastic models S1 S2 S1 S2 E[MPa] 1.46 E5 1.46 E5 ¯ηcon 2.8499 1.2 ν0.2 0.2 v1.15 1.15 pa[kPa] 1. 1. k10.2 0.2 qa[kPa] 1. 1. k20.7256 0.7256 pcap,0 [kPa] 5. E3 5. E3 ccon 4.3067 E-2 4. E-3 e0.7 0.7 ccap 1.59 E-4 1.59 E-4 m7.423 E-2 7.423 E-2 l1.0867654 1.0867654 γ0.5 0.5 r1.592 1.592 α0.8 0.8 ;7.5 E-5 7.5 E-1 Table 4.2: Sets of material parameters. S1 is a Sacramento River sand (Macari et al. 1997). S2 is a modification of S1. First part of subsection 4.2.5 deals with local problems. The relative error of the vector of unknowns x=(σ,κ,λ)measuredinthemaximumnormisusedtocontrolthe convergence. Global problems are treated in second part of subsection 4.2.5. Convergence is checked with the relative error in energy norm. All the computations (except where the opposite is explicitly stated) have been performed by using a negative power of 2, not of 10, as the relative stepsize (hr=2 −k), to reduce the effect of rounding errors. However, relative stepsizes are expressed as powers of 10 to indicate clearly the order of magnitude. For instance, hr=10 −6in the text or in a table means that the actual computation is performed with hr=2 −19. Strict tolerances have been used: 10−14 for local problems and 10−8to 10−10 for global problems. This allows for a comparative assessment of the three differentiation techniques. Quite larger values may be chosen in practice. Local problems In the local problem, numerical differentiation is applied to compute the Jacobian shown in equation (4.2.11) at each Gauss point. In order to show that quadratic convergence is obtained in all stress–internal variable space, three different deformation paths are considered. The paths are characterized by an initial stress–internal variable state, σini and κini, and a total strain increment, ∆ε (applied in 50 steps) see table 4.3. The material parameters of soil S2 have been used. Path A Path B Path C σini (1000,1000,1000,0) (4800,4800,4800,0) (4800,4800,4800,0) κini (0,0) (0,0) (0,0) ∆ε(0,0,0,0.2) (0,0,0,0.2) (−0.1,0,0,0) n50 50 50 Table 4.3: Definition of the three stress paths for the local problems. Elastoplastic models 127 0 1000 2000 3000 0 2500 5000 7500 10000 S T 0 1000 2000 3000 4000 5000 0 2500 5000 7500 10000 12500 15000 S T 0 1000 2000 3000 4000 5000 0 2500 5000 7500 10000 12500 15000 S T $ % & Figure 4.6: Trace of the three paths and of their initial and final yield criteria on the meridian plane. In figure 4.6 the trace on the meridian plane of the three stress paths and the trace of the initial and final yield criteria are depicted. Paths A and B correspond to pure shear deformation, see ∆εin table 4.3. Path A develops in the cone region, and path B starts in the cap region and then changes to the cone region. Path C corresponds to uniaxial compression, see table 4.3, and it develops in the cap region. In table 4.4 the evolution of the Lode angle, θ, is shown. Note that, in general, the paths are three–dimensional curves in the three–invariant space (p,q,θ). Indeed, θchanges during loading in paths A and B. For path C, on the other hand, θremains constant and equal to π/3. Finally, note that the three paths start in hardening regime and finish during softening. Therefore, a wide range of different local problems is covered. 128 Elastoplastic models C 1,E-17 1,E-13 1,E-09 1,E-05 1,E-01 1 2 34 56 7 Iteration Relative error Step 11 Step 21 Step 31 Step 41 A 1,E-17 1,E-13 1,E-09 1,E-05 1,E-01 1234567 Iteration Relative error Step 1 Step 11 Step 31 Step 41 B 1,E-17 1,E-13 1,E-09 1,E-05 1,E-01 1 2 34 56 7 Iteration Relative error Step 11 Step 21 Step 31 Step 41 Figure 4.7: Convergence results for various steps of paths A, B and C In the following, quadratic convergence results are presented and analyzed for moderate strain increments. After that, the behaviour of the three numerical differentiation schemes, see table 4.1, is compared. Finally, quadratic convergence results for large excursions outside the elastic domain are also shown. Convergence illustration The convergence results for different steps of the three stress paths are depicted in figure 4.7. These results correspond to local problems in the cone and the cap regions and to hardening and softening regimes. All the convergence results are quadratic up to a (very strict) tolerance of 10−14. These results have been obtained with the approximation 1ND-O(h) and with hr=10 −6, see equation (4.2.13) and table Elastoplastic models 129 Step 1 5 10 25 50 Path A 40.3◦51.5◦53.1◦54.1◦53.9◦ Path B 30.0◦35.8◦42.2◦54.7◦54.1◦ Path C 60.0◦60.0◦60.0◦60.0◦60.0◦ Table 4.4: Evolution of the Lode angle θ(in degrees) during the three stress paths defined in table 4.3. 1,E-20 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 1,E+04 1234567 Iteration Increm ent Sigm a_xx Sigm a_yy Sigm a_zz Sigm a_xy k_cone k_cap lam bda Figure 4.8: Convergence of the different unknowns of the local problem for step 31 of path A. 4.1. The convergence is also quadratic if checked independently for each unknown of the local problem (σ,κand λ), see figure 4.8. In figure 4.9 the stress invariants and the yield criterion during the iterations are depicted. In three iterations the approximations are very close to the final result. The remaining iterations are just to improve the accuracy. Comparison of numerical differentiation schemes The three numerical differentiation schemes of table 4.1 have been compared through the integration of paths A, B and C defined in table 4.3. The convergence results of step 31 of path A with different relative stepsizes hrare depicted in figure 4.10. They are quadratic up to a tolerance of 10−14 for a wide range of hrwith the three schemes. The same results are obtained with the other steps of paths A and B. In table 4.5, the ranges of hrthat give quadratic convergence during all the steps of paths A and B are summarized. The main difference between Num. approx. Path A Path B 1ND-O(h) 10−6—10 −910−6—10 −10 1ND-O(h2) 10−4—10 −10 10−4—10 −10 1CND-O(h2) 10−4—10 −11 10−4—10 −11 Table 4.5: Range of relative stepsizes hrthat give quadratic convergence in the local problem, stress paths A and B. 130 Elastoplastic models 800 900 1000 1100 1200 1300 1000 1100 1200 1300 1400 1500 p q 850 950 1050 1150 1250 4500 4600 4700 4800 4900 p q A step 11 B step 3 1 2 2,3 1 1 1 2,3 2,3 3 trial = 45.46o 1 = 56.41o 2 = 53.25o 3 = 53.32o trial = 30.92o 1 = 32.34o 2 = 32.95o 3 = 33.04o nFcon nFcap trial trial Figure 4.9: Evolution of the stress invariants and the yield criterion during the iterations of step 11 of path A and step 3 of path B. the three techniques is that second order of accuracy provides quadratic convergence with larger hr. This is in agreement with section 3.1 and it is due to the fact that the truncation error of the second–order schemes is lower than for the first–order scheme. The ranges of hrthat give quadratic convergence during all the steps of path C are summarized in table 4.6. Both the ranges obtained using hr=10 −kand hr=2 −kare indicated. Two aspects are important: first, the ranges are quite narrower than for paths A and B; and second, with the approximation 1ND-O(h) the improvement of using hr=2 −k is notorious. This is because path C develops at Lode angle equal to π/3. In this zone, the influence of the rounding errors is quite more important than in the other regions of the stress space. Nevertheless, the range in which one can choose hris still wide enough and includes the approximation indicated before. Elastoplastic models 131 Num. approx. hr=10 −khr=2 −k 1ND-O(h) 10−5—10 −610−5—10 −8 1ND-O(h2) 10−3—10 −710−3—10 −7 1CND-O(h2) 10−4—10 −910−4—10 −9 Table 4.6: Range of relative stepsizes hrthat give quadratic convergence in the local problem, stress path C. 1CND-O(h2) 1,E-17 1,E-13 1,E-09 1,E-05 1,E-01 12345678 Iteration Relative error hr = 1.E-3 hr = 1.E-6 hr = 1.E-8 hr = 1.E-12 1ND-O(h) 1,E-17 1,E-13 1,E-09 1,E-05 1,E-01 12345678 Iteration Relative error hr = 1.E-3 hr = 1.E-6 hr = 1.E-8 hr = 1.E-12 1ND-O(h2) 1,E-17 1,E-13 1,E-09 1,E-05 1,E-01 12345678 Iteration Relative error hr = 1.E-3 hr = 1.E-6 hr = 1.E-8 hr = 1.E-12 Figure 4.10: Convergence results for step 31 of path A using the approximations defined in table 4.1 with several relative stepsizes hr 132 Elastoplastic models 1,E-16 1,E-12 1,E-08 1,E-04 1,E+00 123456789101112 Iteration Relative error Step 1 Step 3 Step 5 Step 7 Step 9 Figure 4.11: Convergence results for path A with only 10 steps Large excursions outside the elastic domain In order to show that quadratic convergence is also attained for large excursions outside the elastic domain, path A defined in table 4.3 is solved with only 10 steps. The convergence results for various steps, depicted in figure 4.11, are again quadratic. On the other hand, note that, as expected, more iterations than in the original integration of path A with 50 steps are needed (compare Figures 4.11 and 4.7). This is clearly due to the use of the solution of one step as the initial approximation for the next step. Larger steps mean worse initial approximations and, thus, more iterations. To summarize these examples on local problems: quadratic convergence can be attained in a simple manner with any of the three techniques of numerical differentiation. This is valid for any stress path (cone and cap regions, hardening and softening regimes), and for both moderate and large steps. The choice of the stepsize presents no difficulties, because quadratic convergence is obtained for a wide range of relative stepsizes. Global problems In the following, numerical differentiation is applied to solve several boundary value problems (i.e., global problems). That is, the numerical approximations of table 4.1 are employed to compute consistent tangent matrices, see equation (4.2.12). Moreover, they are also used to solve the corresponding local problems. Three examples are presented: the vertical displacement of a pile, a triaxial test with an homogeneous sample and a triaxial test with a non-homogeneous sample. These examples illustrate that the three numerical approximations to first derivatives of the flow vector and the hardening moduli, see table 4.1, are useful to solve the global problem with quadratic Elastoplastic models 133 0 200 400 600 800 1000 1200 1400 1600 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 Displacem ent [cm ] Load [kN] Depicted convergence results Figure 4.12: Load versus displacement curve for the pile problem. convergence. Moreover, their main features (range of adequate relative stepsizes and computational cost) are compared. Vertical displacement of a pile The first example is the vertical displacement of a pile. The definition of the problem is presented by Potts and Gens (1985), and it is only summarized here. Figure 3.8 shows the finite element mesh. It corresponds to a horizontal disc of soil. The thickness of the disc is 5 cm and the pile radius is 7.5 cm. A hydrostatic initial stress state of −250 kPa is imposed. To model the loading of the pile, a vertical displacement of 0.625 cm is prescribed over the boundary AF in 25 load steps. To model the infinite extension of the disc, zero vertical displacements are prescribed over the boundary CD. Due to the essentially one–dimensional nature of the problem, vertical lines (such as EB) are prescribed to remain vertical during loading. The material parameters correspond to the dense Sacramento River sand, see table 4.2. The load versus displacement curve is depicted in figure 4.12. The most stressed points are those next to the boundary AF. Because of the one–dimensional behaviour of the problem, the limit state is reached when the integration points close to AF start the softening regime. After that, stresses are no longer transferred to the rest of the disc. Thus, this simple example only tests the behaviour during hardening. The convergence results for several load steps (indicated in figure 4.12) are shown in figure 4.13. Convergence is quadratic up to a strict tolerance of 10−10. These particular results have been obtained with the approximation 1CND-O(h2) and with hr=10 −5. However, similar results are obtained with the other techniques and other relative stepsizes. In table 4.7, the ranges of hrthat give quadratic convergence during all the test are summarized.