- Email: [email protected]
Cyclic plasticity using Prager’s translation rule and both nonlinear kinematic and isotropic hardening: Theory, validation and algorithmic implementation
Accepted Manuscript Cyclic plasticity using Prager’s translation rule and both nonlinear kinematic and isotropic hardening: Theory, validation and algorithmic implementation Meijuan Zhang, Jos´e Mar´ıa Ben´ıtez, Francisco J. Mont´ans
PII: DOI: Reference:
S0045-7825(17)30655-2 https://doi.org/10.1016/j.cma.2017.09.028 CMA 11622
To appear in:
Comput. Methods Appl. Mech. Engrg.
Received date : 24 May 2017 Revised date : 5 August 2017 Accepted date : 18 September 2017 Please cite this article as: M. Zhang, J.M. Ben´ıtez, F.J. Mont´ans, Cyclic plasticity using Prager’s translation rule and both nonlinear kinematic and isotropic hardening: Theory, validation and algorithmic implementation, Comput. Methods Appl. Mech. Engrg. (2017), https://doi.org/10.1016/j.cma.2017.09.028 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
*Manuscript Click here to download Manuscript: Zhang_Benitez_Montans_CMAME1_R1_vf.pdf
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Click here to view linked References
Cyclic plasticity using Prager’s translation rule and both nonlinear kinematic and isotropic hardening: theory, validation and algorithmic implementation Meijuan Zhanga , Jos´e Mar´ıa Ben´ıteza , Francisco J. Mont´ansa,∗ a
Escuela T´ecnica Superior de Ingenier´ıa Aeron´ autica y del Espacio Universidad Polit´ecnica de Madrid Plaza Cardenal Cisneros, 3, 28040-Madrid, Spain
Abstract Finite element analysis of structures under elasto-plastic nonproportional cyclic loadings is useful in seismic engineering, fatigue analysis and ductile fracture. Usual models with nonlinear stress-strain curves in cyclic behavior are based on Mr´oz multisurface plasticity, bounding surface models or models derived from the ArmstrongFrederick rule. These models depart from the associative Prager’s rule with the main purpose of modeling aspects of cyclic nonlinear hardening. In this paper we develop a model for cyclic plasticity within the framework of the associative classical plasticity theory using Prager’s rule accounting for anisotropic nonlinear kinematic hardening coupled with nonlinear isotropic hardening. We include the validation of the theory against several uniaxial and multiaxial cyclic experiments and an efficient fully implicit radial return algorithm. The parameters of the model are obtained directly by a discretization of the uniaxial stress-strain behavior. Remarkably, both the presented theory and the computational algorithm automatically recover classical bilinear plasticity and the Krieg and Key algorithm if the user-prescribed stress-strain curve is bilinear. Keywords: Cyclic plasticity, Associative plasticity, Radial return algorithm, Multisurface model, Nonlinear kinematic hardening
∗
Corresponding author Email addresses: [email protected] (Meijuan Zhang), [email protected] (Jos´e Mar´ıa Ben´ıtez), [email protected] (Francisco J. Mont´ ans )
Preprint submitted to Computer Methods in Applied Mechanics and Engineering
August 5, 2017
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
1. Introduction The finite element analysis of structures under nonproportional and cyclic nonlinear behavior is very important in many applications, among them seismic engineering [1], [2], [3], fatigue and fracture analyses [4], [5], [6], [7] and plastic forming [8], [9], [10], [11]. The correct description of multiaxial hardening effects have proved to be critical for an accurate prediction of the effective displacements, accelerations and safety of the structures [12]. In fatigue analysis of notched specimens, the plastic loading at the notch edge induces nonproportional multiaxial loading [13], frequently studied through cyclic plasticity models which, thereafter, are related to nominal strains and stresses through incremental Neuber rules or related strain or energy methods [5], [14], [15], [16], [17], [18]. In cyclic behavior, the accurate representation of the Bauschinger effect is crucial, a behavior which is approximately represented by Masing’s rules. The classical linear kinematic hardening reproduces accurately such rules, has a very efficient integration algorithm and is available in most finite element programs. This model uses Prager’s translation rule for the backstress, which is consistent with the principle of maximum dissipation [19]. However, a priori, the direct inclusion of a nonlinear anisotropic kinematic hardening function produces unphysical loops [20]. Therefore, finite element programs use Prager’s rule exclusively in the case of linear kinematic hardening, and for nonlinear kinematic hardening, which is needed to better describe the cyclic loops, they resort to other types of formulations (if available to the user), typically based on the Armstrong and Frederick rule [21], [22]. One of the procedures to account for anisotropic, history-dependent nonlinear kinematic hardening is based on the history of observable quantities by hereditary integrals, as in the endochronic theory, similar to those used in viscoelasticity [23]. However, since this approach is more cumbersome for computational application and for finite element analysis, the approach based on internal variables [20], where only information of the previous step is needed, is usually preferred. Remarkably, the theories to model nonlinear kinematic hardening based on internal variables have departed from the classical linearly hardened framework, introducing different non-associative translation rules for the backstress, as Mr´oz’s rule [24], Garud’s rule [25] or the commented Armstrong and Frederick rule [26]. These rules have resulted in different models as the Mr´oz model [24], Chaboche’s model [20], [27], bounding surface models [28], [29], [1], and nonlinear kinematic hardening models with the addition of multiple backstresses [20], [30]. The algorithmic implementation of some of these models is usually more elaborate than that of classical plasticity [31], [32], [33], [34], [35], and even though some recent efficient algorithms are available for some cases [21], [22], they do not constitute a natural extension of 2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
classical Prager’s plasticity. To improve the cyclic multiaxial plastic behavior, sometimes non-proportionality parameters, obtained for certain materials under certain loading paths by fitting experiments, are employed, like in Refs. [36], [37], [38], [39], [40], [41], [42], [35], [43]. However, despite the common belief [20, p. 206], it is possible to develop nonlinear kinematic hardening models using Prager’s translation rule and preserving Masing’s rules [44], [45]. These models are similar in conception to Mroz’s model, but with some crucial diferences, as the preservation of Prager’s rule, the use of a single yield surface (outer surfaces are hardening surfaces) and the simplicity of the integration algorithms. Furthermore, they lack the inconsistencies present in Mr´oz’s model under multiaxial loading [46], [47]. As we show in this paper, the predictions for multiaxial plastic behavior using Prager’s rule are similar to the behavior observed in the experiments, at least for the different materials addressed below. Despite of the employed kinematic hardening rule, mixed kinematic-isotropic hardening models have been increasingly used in cyclic loading analysis [48], [49], [50], [51], [52], [53], [54], [55], [10]. For a better description of cyclic hardening/softening in certain materials, the effect of a memory for the strain amplitude can be included in the isotropic hardening rule [27], [56], [57], [42], and the inclusion of isotropic cyclic softening may be used to model the effects of damage by fatigue. However, as seen below, the presence of isotropic hardening also modifies the anisotropic kinematic hardening moduli, an observation that must be considered in the formulation and computational algorithm. Therefore, it is valuable, and the purpose of this paper, to extend the classical von Mises associative theory of plasticity using Prager’s translation rule to account for history-dependent nonlinear anisotropic kinematic hardening combined with cyclic isotropic hardening/softening. In the next sections we first introduce the theory, showing that we just propose a simple extension of the classical framework by allowing a dependence of the effective hardening moduli on some internal variables. Then we introduce the method to compute such dependence, i.e. the effective hardening modulus at each instant. Thereafter we explain the efficient radial return algorithm, which for the bilinear case naturally recovers the solution from Krieg and Key [58]. With some examples we show that the theory predicts rather well the experimentally observed multiaxial behavior in several materials and that the finite element implementation preserves the asymptotic quadratic convergence of Newton schemes.
3
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
2. Classical J2 –plasticity using Prager’s rule with nonlinear anisotropic mixed hardening 2.1. Classical associative plasticity with anisotropic, cyclically modified hardening The stored “observable” (elastic) energy depends on the elastic strains εe , which are a function of the total strains and some internal strains εp (the plastic strains), i.e. εe (ε, εp ). The chain rule derived from the explicit dependencies εe (ε, εp ), gives immediately the rate as ∂εe ∂εe : ε˙ + : ε˙p = ε˙e |ε˙p =0 + ε˙e |ε=0 (1) ε˙e = ˙ ∂ε ε˙p =0 ∂εp ε=0 ˙ ≡
tr
ε˙e +
ct
ε˙e ≡ ε˙ − ε˙p
(2)
where the last identity accomodates the rate form of the usual additive decomposition ansatz εe = ε − εp . We defined the trial and corrector partial derivative contributions tr ε˙e ≡ ε˙e |ε˙p =0 and ct ε˙e ≡ ε˙e |ε=0 ˙ . In the new formulation of large strain elasto-plasticity presented in [59] and [60], this framework based on the chain rule remains additive and unaltered at large strains even when using the multiplicative decomposition. The stored energy is —see Figure 1 Ψ (εe , ξ, γ) = W (εe ) + H (ξ, γ) = 12 εe : C : εe + H (ξ, γ)
(3)
where W (εe ) is the elastic energy, immediately recovered upon unloading, and H (ξ, γ) is the internal, microstructurally-blocked “elastic” energy of which part can be recovered following certain plastic deformation paths. This energy may be function of a scalar variable γ (accounting for the cummulative plastic strain) and a internal variable ξ of tensor nature; and in turn ξ itself may be a function of some internal variables ξi of microscopic origin, mapped to the continuum level. For elastically isotropic materials, and using the usual observation that plastic flow is isochoric, we write (4) Ψ (εe , ξ, γ) = 21 κ (εv )2 + 21 2µεde : εde + H (ξ, γ)
where εv are the volumetric strains (elastic), εde are the elastic deviatoric strains, µ is the shear modulus, κ is the bulk modulus. We define the backstress α and the overstress rˆ as α (ξ, γ) :=
∂H (ξ, γ) ∂H (ξ, γ) , and rˆ (ξ, γ) := ∂ξ ∂γ
With these definitions, the rate of the stored energy is 4
(5)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 1: Decomposition of the input energy
˙ = σ : ε˙e + α : ξ˙ + rˆγ˙ = κεv ε˙v + 2µεd : ε˙d + α : ξ˙ + rˆγ˙ Ψ e e
(6)
The rate ξ˙ = − ct ε˙e = ε˙p , accounting for the evolution of internal variables, only takes place during plastic flow. The definition ξ˙ = − ct ε˙e is just a traditional, convenient definition so α results in the usual macroscopic, continuum backstress, that can be directly subtracted from the stress in the yield function as seen below. In ˙ a similar we take γ˙ as the norm of ξ,
definitions,
fashion, to accommodate the usual
˙ ˆ and we define n ˆ = ξ/ γ˙ = ξ˙ , so ξ˙ = γ˙ n
ξ˙ , and the last term in Eq. (6) adopts the familiar form for isotropic hardening. Note also that sometimes the sign of both ξ and α are opposite than in our definition here to understand ξ as internal elastic strains mapped to the continnum (because of the “elastic” nature of H (ξ, γ)), but in practice the numerical result remains to be the same. Denoting the stress tensor by σ, the stress power is P = σ : ε. ˙ Using ε˙ ≡ tr ε˙e , the dissipation is ˙ =σ: D˙ = P − Ψ
ε˙e − κ (εv ) ε˙v − 2µεde : ε˙de − α : ξ˙ − rˆγ˙ ≥ 0 (7) ˙ = tr Ψ− ˙ ˙ ≡ − ct Ψ which can be interpreted as D˙ ≡ − ∂Ψ (ε, εp ) /∂εp |ε=0 : ε ˙ = − Ψ p ˙ ε=0 ˙ ˙ Following the usual Coleman procedure, for purely elastic incremental deformaΨ. tions, when D˙ = 0, we also have ε˙p = 0 and by definition ct ε˙e = −ξ˙ = 0. Thus by chain differentiation we immediately identify the necessary condition σ = σv + σd = tr
5
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Kεv I + 2µεde , where σv and σd are the volumetric and deviatoric parts, respectively. Then, since γ˙ = (dγ/dξ) : ξ˙ = n ˆ : ξ˙ = − n ˆ : ct ε˙e and tr ε˙e = ε˙e − ct ε˙e , for the plastic case, the inequality D˙ > 0 can be written as D˙ = − σd − α − rˆn ˆ : ct ε˙e > 0 (8)
Hence, the ddissipation is maximum, and the inequality is always guaranteed when ct ε˙e and σ − α − rˆn ˆ have the opposite direction, i.e. in this case n ˆ may also be defined as σd − α (9) n ˆ= kσ d − αk
so the flow rule − ct ε˙e ≡ ε˙p = γ˙ n ˆ and the hardening rule ξ˙ = − ct ε˙e = γ˙ n ˆ are defined in terms of this tensor n ˆ . The multiplier γ˙ ≥ 0 is now interpreted as the norm of the plastic strain rate (the equivalent plastic strain rate, up to a constant). This variable is non-negative because D˙ ≥ 0. Then we can write
D˙ = σd − α − rˆn ˆ γ˙ ≡ 0 r γ˙ ≡ σY p˙ ≥ 0 (10)
where we have defined 0 r := σd − α − rˆ ≥ 0, and k·k denotes the 2-norm. The p term σY p˙p is the uniaxial equivalence, where σY = 3/2 0 r is the uniaxial yield stress and p˙ = 2/3γ˙ is the uniaxial plastic strain rate, obtained from r Z t r 2 2 p= kε˙p k dt = γ (11) 3 0 3 Then, using σ = α + (0 r + rˆ) n ˆ and ε˙p = γ˙ n ˆ , the stress power P may be written as the sum of the following terms σ : ε≡σ ˙ : tr ε˙e = σ : ε˙e − ct ε˙e = σ : ε˙e + σ : ε˙p (12) ˙ Ψ D˙ }| { z z}|{ (13) ≡ P = σ : ε˙e + α : ξ˙ + rˆγ˙ + 0 r γ˙ | {z } | {z } |{z} ˙ W H˙ kin H˙ iso | {z } ˙ H
6
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
so the integral of these terms are the energies shown in Fig. 1: Z Z Z Z ˙ W = σ : ε˙e dt, Hkin = α : ξdt, Hiso = rˆγdt, ˙ D=
0
r γdt ˙
(14)
R and the area below the stress-strain curve, Pdt, is the external (“stress”) work. The energy W is immediately recovered upon unloading. The energy Hkin may be recovered when plastic flow happens in the reverse direction; note that the product α : ξ˙ may have positive or negative sign. The energy Hiso is usually not recovered during cyclic loading (isotropic hardening) because rˆγ˙ is positive for hardening (as seen below); however, it may be recovered through thermal treatments. Finally, the dissipation energy D is never recovered; it is dissipated at each instant, typically, into heat. Sometimes Hiso is considered as part of the dissipation (phenomenologically, if thermal effects are neglected, there is no difference), so in such case we could define D˙ = (0 r + rˆ) γ. ˙ Equation (10) may be written in the following convenient format for identification of the yield surface
D˙ = σ d − α − rˆ − 0 r γ˙ + 0 r γ˙ ≥ 0 (15) which implies that (case of plastic loading)
f := σd − α − r = 0 if γ˙ > 0
(16)
where we defined the current radius r := 0 r + rˆ, but we can have (case of elastic loading or neutral loading) f ≤ 0 if γ˙ = 0 (17) During plastic deformations, the consistency parameter γ˙ > 0 can be obtained from the consistency condition which must be f = 0 during plastic flow (γ˙ > 0), so f˙ = 0 (or equivalently γ˙ f˙ = 0 in any case, elastic or plastic) tr d (18) +n ˆ : ct σ˙ d − α˙ − r˙ = 0 f˙ (ε, εp ) = n ˙} |ˆ : {z σ {z } | tr ˙ ct ˙ f ≡ f˙ f ≡ f˙ ε˙p =0 ε=0 ˙
where we define the trial deviatoric stress rate as tr
σ˙ d := 2µtr ε˙d ≡ 2µε˙d
7
(19)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
and the corrector one as ct
σ˙ d = 2µ ct ε˙e ≡ −2µε˙p = −2µγ˙ n ˆ
(20)
The latter is the rate of the stress tensor during plastic flow which produce the maximum dissipation. Note that this novel continuum framework presented in [59], [60] mimics the usual algorithmic implementation, even at large strains, but it is just an immediate application of the chain rule when considering the dependencies σ(ε, εp ), i.e. ∂σ ∂σ : ε˙ + : ε˙p = tr σ˙ + ct σ˙ (21) σ˙ (ε, εp ) = ∂ε ε˙p =0 ∂εp ε=0 ˙
This type of decomposition is also recovered at large strains and its backward-Euler discretization results immediately in the typical Closest Point Projection algorithms. The component of the rate of α perpendicular to n ˆ does not affect the yield condition and does not change the dissipation, so as it is well known, the direction of minimum accumulation of internal elastic energy, and maximum dissipation for α, ˙ is the one given by Prager’s rule; see e.g. [19]. Accordingly, we assume Prager’s translation rule, so α ˙ is proportional to n ˆ ∂ 2 H (ξ, γ) ˙ ∂ 2 H (ξ, γ) α˙ = :ξ+ γ˙ ∂ξ∂ξ ∂ξ∂γ ∂α ˙ ∂α ˙ n =: 2 H ¯ (ξ, γ) γ˙ n :ξ+ γ˙ =: λˆ ˆ = 3 ∂ξ ∂γ
(22)
¯ γ, ˙ where ∂α/∂γ = ∂r/∂ξ. The backstress rate modulus is given as λ˙ = kαk ˙ =: 32 H or equivalently, multiplying Eq. (22) by n ˆ γ˙ =
hα˙ : n ˆi (ξ, γ)
2 ¯ H 3
(23)
and we defined 2 ¯ H 3
(ξ, γ) := n ˆ:
∂ 2 H (ξ, γ) ∂ 2 H (ξ, γ) :n ˆ+ :n ˆ ∂ξ∂ξ ∂γ∂ξ
(24)
as the effective kinematic hardening. The Macaulay bracket in hα˙ : n ˆ i implies that γ˙ ≥ 0. Then, from the consistency condition, Eq. (18), we restrict the formulation to a dependence r (γ) —a dependence of r on ξ would include a path-dependent 8
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
isotropic hardening, found in some materials [36], but which for simplicity we do not consider in this work— so r˙ :=
∂ 2 H (ξ, γ) ∂ 2 H (ξ, γ) dr γ˙ ≡ r ′ (γ) γ˙ := : γ˙ n ˆ+ γ˙ dγ ∂γ∂ξ ∂γ∂γ
(25)
where, note, both types of hardening are implicitly coupled through the term ∂ 2 H/∂γ∂ξ as it will be seen below. Then, we can factor out the consistency parameter as usual γ˙ =
2µ d n ˆ : ε ˙ ¯ (ξ, γ) + r ′ (γ) 2µ + 23 H
(26)
¯ (ξ, γ) and r ′ (γ) are the functions to be obtained from experiments. Then, where H it is immediate to obtain the continnum tangent Cep
(2µ)2 = κI ⊗ I + 2µP − ˆ ⊗n ˆ ¯ (ξ, γ) + r ′ (γ) n 2µ + 32 H
(27)
Remarkably, these are the classical J2 plasticity expresions given, for example, in Sec. 2.3.2 of [19] or in Sec. 4.3 of [63], where the previous equations are obtained ¯ Then, the single important from an equivalent procedure but assuming a constant H. difference in Eqs. (26) and (27) respect to the classical ones is that we do not assume ¯ (ξ, γ) is constant, but instead that it changes according to the flow direction that H and the history of deformation by way of the implicit dependences on ξ and γ, as ¯ results constant, both the shown below. However, in the particular case that H theory and computational algorithms in [19] (see also [64]) are naturally recovered regardless of the number or surfaces employed, and without user intervention; i.e. the solution is reached in just one iteration if r ′ is also constant. Therefore, our model constitutes a natural extension of the classical computational algorithms in finite element software to account for history-dependent non-linear kinematic hardening and anisotropic hardening effects. 2.2. Description of anisotropic hardening The anisotropic hardening is described in our model with the aid of different, initially concentric, hardening surfaces (i.e. not macroscopic yield surfaces). We remark that we follow the idea of Mr´oz and Iwan in using concentric surfaces for describing the hardening, but our theory is different from Mr´oz’s theory because our hardening surfaces are not yield surfaces at the macroscopic level, but just a tool ¯ which take to compute an effective macroscopic anisotropic hardening modulus H into account microstructural changes due to the deformation history in the material. 9
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
a.
b.
Yield surface
Hardening surfaces
Mróz-like view of the theory
Hereditary view of the theory
c. Trial stress
Yield surface (final stress)
Backward-Euler / radial-return stress integration
Figure 2: Geometric description of the model and the integration algorithm in the deviatoric plane. (a) Multisurface (Mr´ oz-like) view — Yield and hardening surfaces fi , flow direction n ˆ , backstresses αi , contact points σ ˆ i and hardening directions m ˆ i . (b) Hereditary (endochronic-like) view — Yield surfaces f¯i for internal variables αi . (c) Backward-Euler, radial-return integration algorithm; for simplicity only kinematic hardening is considered in the figure.
10
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Furthermore, we keep Prager’s hardening rule instead of changing to Mr´oz’s or other alternative rule. These hardening surfaces may be written as, see Figure 2a fi = kˆ σ i − αi k − ri = 0, i = 1, ..., N
(28)
where N is the number of surfaces and σ ˆ i is the contact point in the deviatoric stress space of surface i with surface (i − 1), and we define the first, innermost surface as the (macroscopic) yield surface, i.e. f1 ≡ f , σ ˆ 1 ≡ σd , α1 ≡ α, r1 ≡ r
(29)
As another phenomenological difference with Mr´oz’s proposal, the contact point between surfaces is not, in general, the same point, i.e. σ ˆ i 6= σ ˆ j for i 6= j. Equations (28) are convenient for visualization. However, they are fully equivalent to the following equations, more convenient for mathematical derivations, see Figure 2b
f¯1 ≡ f1 = σd − α1 − r1 = 0 (30) f¯i = kαi − αi−1 k − ρi = 0, i = 2, ..., N
where we have defined ρi = ri − ri−1 > 0 (and ρ1 = r1 ). Hence, the hardening surfaces can be alternatively interpreted as yield surfaces for internal microstructural stresses αi . These internal stresses are a result of internally stored (blocked) energy, which is returned upon stress reversal. Interestingly, the view in Figure 2b shows an interpretation of the theory close to that of hereditary models, as for example the endochronic theory [23], [67]. 2.3. Consistency condition in hardening surfaces For those surfaces such that fi < 0 (or equivalently f¯i < 0), the hardening surfaces remain fixed because the inner surfaces are fully inside, not touching surface i. However, if fi = 0 for any i > 1, this implies that surface (i − 1) is touching surface i. In such a case, we must guarantee that f˙i ≤ 0 (or equivalently df¯i /dt ≤ 0) in order to avoid overlapping surfaces which would yield an inconsistent definition of the hardening function in the stress domain. From Eq. (302 ) αi − αi−1 df¯i = : (α ˙i −α ˙ i−1 ) − ρ˙ i = 0 dt kαi − αi−1 k
11
(31)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Define —this tensor can be interpreted as the flow direction for the internal stresses αi of microstructural origin, see Fig. 2b m ˆ i :=
αi−1 − αi mi := kmi k kαi−1 − αi k
and m ˆ1 ≡n ˆ
(32)
so
df¯i =m ˆ i : (α ˙i −α ˙ i−1 ) + ρ˙ i = 0 (33) dt The first surface is the actual yield surface, so as seen above, it translates according to Prager’s rule, α ˙ 1 = kα˙ 1 k n ˆ . Then the consistency on the second surface is −
m ˆ 2 : α˙ 2 − kα ˙ 1 k (ˆ n:m ˆ 2 ) + ρ˙ 2 = 0
(34)
Note that as in the case of the yield surface, the component of α˙ 2 perpendicular to m ˆ 2 is irrelevant here, but the minimum translation for a given kα˙ 1 k takes place for translation in direction m ˆ 2 , i.e. α˙ 2 = kα ˙ 2k m ˆ 2 —this can be seen as Prager’s rule ¯ applied to fi , interpreted as yield conditions at different microstructural levels kα ˙ 2 k = kα˙ 1 k (ˆ n:m ˆ 2 ) − ρ˙ 2
(35)
kα˙ i k = kα ˙ i−1 k (m ˆ i−1 : m ˆ i ) − ρ˙ i ≥ 0
(36)
In general, It is important to note that if the position of the inner surface is known and contact is taking place, from f¯i = 0, Eq. (302 ), and using α˙ i = kα ˙ ik m ˆ i , we can obtain the position of the outer surface in an explicit manner, see Figure 2 αi = αi−1 − ρi m ˆi
(37)
which is a condition easier to use in the algorithmic (discrete) implementation. Note that defining ρi = ρi m ˆ i , we can also have the interpretation of the backstress as the addition of several overstresses (which could in turn be interpreted as an internal equilibrium equation) α = αa+1 +
a X
ρi and σ = αa+1 +
a X i=1
i=2
where a is the outermost active surface.
12
ρi
(38)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
2.4. Consistency parameter from multiple hardening surfaces Equation (23) is directly applicable for one surface. In the presence of N hardening surfaces, we assume the existence of some strain-like internal variables ξi , i = 1, ..., N (of microstructural origin) describing an anisotropic hardening field, such that we have the dependency ξ (ξi ), and in particular we assume that they follow the compatibility equation ξ˙ =
N X
ξ˙i
(39)
i=1
In this work we do not include ratcheting effects. However, these effects may be included in the theory by the addition of ratchetting internal strains in the previous equation. From the compatibility equation of internal strains, we have at the continuum observable level that the consistency parameter rate is γ˙ = ξ˙ : n ˆ=
a X
ξ˙i : n ˆ=
a X
γ˙ i
(40)
i=1
i=1
where we defined γ˙ i := ξ˙i : n ˆ and a ≤ N is the outermost active surface, i.e. a = arg max (fi = 0) i
(41)
If we associate a hardening modulus Hi (γ) to each one of the hardening surfaces such that α ˙i (42) α˙ i = 32 Hi (γ) ξ˙i so ξ˙i = 2 H (γ) 3 i
then by Eq. (40)
a
a
X X hα˙ i : n hα ˙ :n ˆi ˆi γ˙ = 2 ¯ = γ˙ i = 2 H (ξ, γ) Hi (γ) 3
i=1
(43)
i=1 3
In this work we restrict the dependence of the moduli Hi (γ) to deppend only on the ¯ (ξ, γ) is implicitly defined by cummulative parameter γ. The other dependency in H
13
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
the history of the backstresses αi and their flows ξ˙i . Using Eq. (36) recursively γ˙ i =
hα ˙i : n ˆi kα ˙ ik = 2 hm ˆi :n ˆi 2 H H 3 i 3 i i
Qi hm ˆi :n ˆ i X Qi−1 hm ˆi :n ˆi kαk ˙ ˆj :m ˆ j−1 i − 2 ρ˙ j k=j (m ˆk :m ˆ k+1 ) = 2 j=2 hm H H 3 i 3 i j=2
(44)
The first term in Eq. (44) is the most important one, whereas the other terms represent the coupling between isotropic and kinematic hardening, and vanish even with isotropic hardening in the case that r˙j = r˙j−1 . Since i X j=2
ρ˙ j = r˙i − r˙1
(45)
for the case of uniaxial loading, we have m ˆi =n ˆ for all i, kαk ˙ r˙i − r˙1 − 2 2 H H 3 i 3 i
(46)
kαk ˙ (r˙i − r˙1 ) / kαk ˙ kαk ˙ − kαk ˙ = 2 [1 − ζia (p)] 2 2 H H H 3 i 3 i 3 i
(47)
γ˙ i =
γ˙ i =
where we have defined ζia (p) :=
r˙i − r˙1 σˆ ′ (p) − σ ˆ1′ (p) = i ¯a kαk ˙ H
(48)
in which σ ˆi′ (p) is the derivative of the uniaxial-equivalent size of surface i as a function ¯ a (with superindex a) is the uniaxial equivalent of the uniaxial plastic strain p, and H hardening modulus when the outermost active surface is surface a. Then kα ˙ ik z }| { 2 kαk ˙ = Hi γ˙ i + (r˙i − r˙1 ) 3
(49)
and, from Eq. (43), in the particular case of proportional loading a
a
X 1 X 3 1 3 γ˙ = kαk ˙ ¯a = ˙ (1 − ζia) γ˙ i = kαk 2 2 H H i i=1 i=1 14
(50)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
This equation results in
a
X 1 − ζa 1 i = a ¯ Hi H i=1
(51)
which is the expression we may use to determine the hardening moduli associated to the surfaces from a uniaxial stress-strain curve. We note that during a monotonic loading it is impossible to distinguish which part of the hardening corresponds to isotropic and which part to kinematic hardening. Only upon stress reversal, identifying the position of the backstress (at half the elastic region), both components can be identified. For a given surface a−1
X 1 − ζa 1 1 1 i = for a = 1, ..., N − aH a a ¯ Ha 1 − ζ (1 − ζa ) H i a i=1
(52)
This expression simplifies in the usual case where isotropic hardening is small respect to kinematic hardening or all the surfaces harden equally, i.e. ζia ≪ 1, to a−1
X 1 1 1 = ¯a − for a = 1, ..., N Ha H H i i=1
(53)
In this case the modulus is given by the uniaxial equivalence ¯a σˆa+1 − σˆa EH (a) ≃ Eep = a ¯ εˆa+1 − εˆa E +H
(54)
(a)
in which E is the elastic Young modulus, Eep is the uniaxial elastoplastic modulus while surface a is the outermost active surface, and σ ˆa and εˆa are the stress-strain data points in which the uniaxial stress-strain curve has been discretized, see Figures 3 and 4. Then (a) E − Eep a ¯ H = (55) (a) EEep (N )
For the case a = N, the value of Eep , and hence HN , is given by the desired limiting hardening modulus. Remarkably, if we insert a hardening surface in a straight section of the stress(i+1) (i) ¯ i+1 = H ¯ i and 1/Hi+1 = 0. Then note that strain curve, we have Eep = Eep , so H γ˙ i+1 = 0 and the surface becomes irrelevant in terms of plastic work. Furthermore, if the uniaxial stress-strain curve happens to be bi-linear, naturally all 1/Hi = 0 except 1/H1 , and γ˙ i = 0 except for γ˙ 1 = γ. ˙ Therefore, the classical linearly-hardened 15
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 3: Example of a discretization of a stress-strain curve and the generation of the hardening √ surfaces; depicted in the 3τ − σ plane.
J2 -plasticity is recovered naturally without explicitly accounting for it, just as a consequence of the prescribed stress-strain curve. This is extremely important in order to obtain consistent predictions for different, arbitrary discretizations, a property not enjoyed, for example, by Mr´oz’s model [68]. 2.5. Cyclic isotropic hardening/softening There are many possible implementations for isotropic hardening/softening in the multisurface model. Since we are dealing with multiple surfaces to describe hardening, these surfaces may be used even to change the shape of the stress-strain curve during cyclic loading, for example from the monotonic backbone curve to a very different skeleton one. However, the purpose of the present paper is to show the capabilities of the model to describe multiaxial loading with cyclic hardening/softening. Therefore, since in the examples we will not account for a relevant change of shape in the skeleton curve, to keep the presentation simple, we consider an isotropic hardening function such that all surfaces remain proportional in size ri (p) = R (p) 0 ri
16
(56)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 4: Procedure for modification of the skeleton curve due to cyclic hardening.
where 0 ri = ri (0), R (p) is the isotropic hardening function and we made, as usual, an abuse of notation in ri (p) = ri (γ) to avoid proliferation of symbols. Then r dri 2 ′ 0 0 ′ r˙i (p) = p˙ = ri R (p) p˙ = R (p) ri γ˙ (57) dp 3 and obviously ρ˙ i (p) = r˙i (p) − r˙i−1 (p) =
0
ρi R′ (p) p˙
(58)
An scheme of the surface growth may be seen in Figure 4. The procedure that we use to determine the function R (p) is the follow-up of the peak stresses during subsequent cycles, which is a procedure useful in fatigue analysis. With the data from the change of peak stresses and the accumulated plastic strain during a cycle, we obtain the relation between R and p. With the knowledge of R (p), and hence of the change in the stresses σ ˆa , the derivative of the radius of any hardening surface dri /dp can be obtained. Hence, the relevant key function is R (p). A popular way of expressing the function R (p) is in relative terms, i.e. R/Rs ([27], [57]), where Rs is the saturated value. During cyclic loading, we can approximate R (p) measuring the peak stresses in succesive strain-controlled cycles. Assume that
17
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
when the strain range of a loop is known, the accumulative plastic strain of that first loop is approximately p0 . Then, during cyclic loading, the accumulated plastic strain is approximately pn = 0.25p0 + np0 , where n is the number of loops. The change of R (p) at different loop counting n, gives pairs {pn , Rn } , which can be immediately used to establish the relation R (p) in any desired format (piecewise linear, spline or adjusting a model function). Once the relation R (p), or more generally, the relation ri (p) is known, the expression for the derivative of the kinematic hardening moduli Hi can be obtained, see Figure 4. The inverse of the uniaxial elastoplastic tangent when surface a is active is, for a given p εa+1 (p) − εa (p) 1 1 1 = = + (59) (a) ¯ a (p) E H σ ˆa+1 (p) − σ ˆa (p) Eep (p) Since the elastic behavior is constant, the variations during a cycle are related by— see Figure 4 1 ∆ [εa+1 (p) − εa (p)] = (60) ∆ [ˆ σa+1 (p) − σ ˆa (p)] E Then, using the chain rule, taking the derivative of Eq. (59) respect to p—note that as it should be expected, this derivative vanishes if all surfaces harden equally in an isotropic manner, so dˆ σa+1 /dp = dˆ σa /dp— ¯ a (p) d 1/H dˆ σa+1 (p) dˆ 1 εa+1 (p) − εa (p) σa (p) (61) = − − dp E [ˆ dp dp σa+1 (p) − σ ˆa (p)]2
where dˆ σa = dp
r
2 dra 3 dp
(62)
Thus, from Eq. (51) we can compute in a recursive manner a−1 ¯ a (p) X d (1/Hi (p)) d 1/H d (1/Ha (p)) = − dp dp dp i=1
(63)
d (1/Hi (p)) dHi = −Hi2 (p) dp dp
(64)
and, obviously
In summary, from R (p) we can immediately compute ri (p), dri /dp and dHi /dp for all surfaces (the hardening of the last surface may be set constant). If a change of the shape of the stress-strain curve from monotonic to cyclic behavior is desired, 18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
different functions for ri (p) can be selected to model such effect. 3. Fully implicit algorithmic implementation In this section we introduce a fully implicit integration algorithm based on the two-parameter Closest Point Projection (CPP) Algorithm [69]. 3.1. Incremental formulation The stress integration algorithm computes the stress t+∆t σ at “time” (or step) t + ∆t with the knowledge of: the elastic strains t εe at time t, the increment of total deformations ∆ε during the step, the history variables t αi , and either t p or t γ. The time-integration scheme, based on the classical radial return algorithm is depicted in Figure 2c. Using the trial stress rate Eq. (19), and defining the algorithmic trial stress tr σ = t σ + ∆tr σ, tr σ ≡ σ t εe + ∆ε = t σ + κ∆εv I + 2µ∆εd (65) Since during the absence of plastic flow the backstresses do not change (i.e. we can define the algorithmic backstresses tr
tr
α˙ i = 0),
αi := t αi , i = 1, ..., N
(66)
The corrector stress increment, using a Backward-Euler (radial-return) integration algorithm, is ∆ct σ = 2µ∆ctεe = −2µ∆γ t+∆t n ˆ (67) where
t+∆t
n ˆ is the normal to the yield surface at the end of the step. Then t+∆t
σ = tr σ − 2µ∆γ
t+∆t
n ˆ
(68)
Using the same integration rule, from Eq. (22) for surface 1 t+∆t
α = tα +
t+∆t
˜ 1 ∆γ1 H
t+∆t
n ˆ
(69)
where for convenience, we have absorbed the uniaxial equivalence factor 2/3 in t+∆t
˜ i := 2 H 3
19
t+∆t
Hi
(70)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Then, note that as in classical plasticity, we can obtain the flow direction directly from the deviatoric trial state—recall that the volumetric part is elastic t+∆t
tr
n ˆ=
t+∆t σd − t α n =: tr d t t+∆t k σ − αk k nk
(71)
However, to obtain t+∆t α in Eq. (69), which accounts for the correction to the dissipation (the energy not dissipated but microstructurally stored), we need to compute the consistency parameters which guarantee material compatibility. In order to do so, we need to develop an integration algorithm. We will be able to establish below a very efficient iterative algorithm based just on two scalars, namely ∆γ1 and ∆γ (i.e. a two-parameter CPP iterative algorithm), regardless of the number of surfaces being employed. The two functions to be minimized are the yield function f and the compatibility requirement g (∆γ, ∆γ1 ) = ∆γ − ∆γ a → 0
(72)
where we defined —note the different meaning of sub- and superindices ∆γ a (∆γ1 ) :=
a X
∆γi (∆γ1 ) = ∆γa (∆γ1 ) + ∆γ a−1 (∆γ1 )
(73)
i=1
The last recursive identity is given for future reference. In order to minimize the function, we need the derivatives respect to these two basic iterative variables, namely ∆γ and ∆γ1 . From Equation (37), using also Eq. (32), after some straightforward algebra, we can obtain dt+∆t αi−1 d∆γ a dt+∆t αi = t+∆t D1i : − t+∆t D2i (74) d∆γ1 d∆γ1 d∆γ1 where the algorithmic projectors and algorithmic radii vector changes are t+∆t t+∆t ρi ρi t+∆t t+∆t D1i = 1 − t+∆t P + t+∆t m ˆ i ⊗ t+∆t m ˆi k mi k k mi k r t+∆t 2 dt+∆t ρi t+∆t dt+∆t ρi R t+∆t 0 d D2i = m ˆ with = ρ i i t+∆t t+∆t t+∆t 3d p d p d p
20
(75) (76)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
in which we defined—note that t+∆t
mi :=
t+∆t
t+∆t
m ˆ i is also recursively defined
αi−1 − t αi and
t+∆t
t+∆t
m ˆ i :=
k
mi t+∆t m k i
(77)
and P is the fourth order deviatoric projector tensor; i.e. P = I− 31 I ⊗ I, where I is the fully symmetric 4th-order identity tensor. Then, taking the derivative of Eq. (43) respect to ∆γ1 , we can compute d∆γ a /d∆γ1 in a recursive manner using Eq. (73)2 . Consider ! b−1 X d∆γb d∆γi d∆γ b−1 d∆γb d∆γ b + = = + (78) d∆γ1 d∆γ1 d∆γ1 d∆γ1 d∆γ1 i=1 After some straightforward algebra we obtain t+∆t
dt+∆t αb−1 ∂∆γ b−1 + d∆γ b d∆γ1 ∂∆γ1 = q t+∆t t+∆t ˜b d∆γ1 n ˆ : t+∆t m ˆ b d t+∆t H Zb t+∆t 1+ n ˆ : t+∆t D2b + t+∆t Zb ∆γb 23 t+∆t H t+∆t H ˜b ˜b d t+∆t p (79) which is applied repeatedly for b = 1, ..., a. In this equation, to account in Eq. (43) for γ˙ i ≥ 0, we defined for convenience the Heaviside function t+∆t t+∆t 1 + sign m ˆ : n ˆ b t+∆t Zb = ∈ {0, 1} (80) 2 Zb t+∆t H ˜b
t+∆t
n ˆ:
t+∆t
D1b :
(b) (b−1) ¯b = H ¯ b−1 and H ˜ b → ∞; then ∆γ b = ∆γ b−1 , Remarkably, if Eep = Eep we have H d∆γ b /d∆γ1 = d∆γ b−1 /d∆γ1 and surface b has no effect in the formulation or the iterations.
3.2. Local iterative algorithm As mentioned, the implicit algorithm is based on the search for the values of ∆γ1 and ∆γ such that the following conditions are met t+∆t ˜ 1 ∆γ1 − t+∆t r1 → 0 f (∆γ, ∆γ1 ) := ktr σ − t αk − 2µ∆γ − t+∆t H (81) t+∆t g (∆γ, ∆γ1 ) := ∆γ a (∆γ1 ) − ∆γ → 0
21
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
The solution for [∆γ, ∆γ1 ]T is obtained by a Newton-Raphson algorithm until convergence is achieved in both f and g,
∆γ ∆γ1
(j+1)
=
∆γ ∆γ1
(j)
−
∂∆γ/∂f ∂∆γ/∂g ∂∆γ1 /∂f ∂∆γ1 /∂g
(j)
f g
(j)
(82)
where [·](j) implies values at iteration j. The derivatives in the tangent matrix in Eq. (82) are obtained as usual from its inverse matrix, whose components are ∂f ˜1 = − t+∆t H ∂∆γ 1 q ˜ 1 q d t+∆t r1 d t+∆t H ∂f = −2µ − 23 ∆γ1 t+∆t − 23 t+∆t ∂∆γ d p d p (83) a d∆γ ∂g = ∂∆γ1 d∆γ1 ∂g = −1 ∂∆γ which right-hand-side terms have already been computed above, see Eqs. (79) and (57). Remarkably, if linear hardening is considered, in which case only one surface ˜ i /dp = 0, dri /dp is constant and ∆γ a = ∆γ1 . Then, the function g is is relevant, dH identically satisfied and the algorithm gives the exact solution in just one iteration because all coefficients in the tangent are constants. As expected, this result is in line with the well-known classical radial return algorithm of Krieg and Key [58], so both the continuum and algorithmic classical plasticity are recovered. In the case when only isotropic hardening is nonlinear, g is also identically satisfied and the classical scalar-based algorithm is also recovered, see Ref. [19]. For the case when multiple surfaces are present, a guess (·)(0) is needed to start (0) the iterations. We always take ∆γ1 = 0. A simplest guess for the active surface is to assume a(0) = 1. If after the computation of the solution we obtain f¯2 > 0, then the index is increased a ← a+1, and the solution is recomputed; and so on. However, if a large number of surfaces is employed, a better guess to save computational time is to set a for the first iterations as the outermost surface such that f¯a = 0 and
t+∆t
n ˆ : tm ˆa >0
(84)
If after reaching a solution f¯a+1 > 0, then the active surface index is increased. On the contrary, if a solution with ∆γa > 0 is not obtained, the active surface index is 22
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
decreased. A layout of the computational algorithm is shown in Table 1 3.3. Consistent tangent for the global iterative algorithm To preserve the efficient asymptotic quadratic convergence rate of implicit Newton algorithms in global equilibrium iterations in finite element simulations, the algorithmic (consistent) tangent moduli must be used. These moduli relate the increment of stresses to the increment of strains according to the algorithmic implementation, i.e. t+∆t
dt+∆t σd dεd dt+∆t σ C = t+∆t = κI ⊗ I + t+∆t d : d ε d ε dε
(85)
where κI ⊗ I is the volumetric part (purely elastic and constant), dεd /dε ≡ P is the also constant deviatoric projector tensor and dt+∆t σd /dt+∆t εd is the relevant nontrivial algorithmic part. This tensor may be obtained taking the derivative of—see Figure 2c t+∆t d ˜ 1 ∆γ1 + t+∆t r1 ) t+∆t n σ = t α + (t+∆t H ˆ (86) as dt+∆t σ d dt+∆t n ˆ d∆γ1 t+∆t ˜ t+∆t ˜ 1 t+∆t n = ( H ∆γ + r ) + t+∆t H ˆ ⊗ t+∆t d 1 1 1 t+∆t d t+∆t d d ε d ε d ε t+∆t ˜ t+∆t d H1 d r1 + ∆γ1 t+∆t n ˆ ⊗ t+∆t d + t+∆t n ˆ ⊗ t+∆t d d ε d ε p Taking into account that ∆εp = ∆γ t+∆t n ˆ = 3/2∆p t+∆t n ˆ , we obtain dt+∆t n ˆ d∆γ dt+∆t εp = ∆γ + t+∆t d t+∆t d d ε d ε d∆γ1
and
dt+∆t p = dt+∆t εp
r
2 3
t+∆t
n ˆ , so
t+∆t
n ˆ⊗
d∆γ1 dt+∆t εd
dt+∆t p dt+∆t p dt+∆t εp : = dt+∆t εd dt+∆t εp dt+∆t εd
(87)
(88)
(89)
Using the chain rule ˜ 1 dt+∆t p ˜1 dt+∆t H dt+∆t r1 dt+∆t r1 dt+∆t p dt+∆t H = and = dt+∆t εd dt+∆t p dt+∆t εd dt+∆t εd dt+∆t p dt+∆t εd
23
(90)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Local iterative algorithm 1. Compute trial stress tr σ, Eq. (65) and both t+∆t n and t+∆t n ˆ , Eq. (71)
2. If tr f := t+∆t n − t r1 ≤ 0, step is elastic; accept trial solution and exit. Else, step is plastic: continue
3. Set ∆γ = ∆γi = 0 and guess a; either a = 1 or by Eq. (84) 4. For each iteration (j) Compute Set
t+∆t
t+∆t
˜i H
t+∆t
α1 = t α1 +
p and
t+∆t
t+∆t
˜ 1 ∆γ1 H
ri
t+∆t
t+∆t
p
with
t+∆t
p=
p
2/3 (t γ + ∆γ)
n ˆ and initialize ∆γ a = ∆γ1
For each hardening surface i = 2, ..., a Compute
t+∆t
mi and
t+∆t
m ˆ i via Eq. (77)
By Eq. (37) t+∆t αi = t+∆t αi−1 − t+∆t ρi t+∆t m ˆi ˜ i t+∆t αi − t αi : t+∆t n ∆γi = 1/t+∆t H ˆ
∆γ a ←− ∆γ a + ∆γi
Compute t+∆t f = t+∆t n − 2µ∆γ −
t+∆t
˜ 1 ∆γ1 − H
t+∆t
r1
Compute t+∆t g = ∆γ a − ∆γ If t+∆t f < tolf and t+∆t g < tolg accept the solution, go to step 5. Perform loop i = 2, ..., a to compute d∆γ a /d∆γ1 via Eq. (79) Compute the local tangent matrix, Eqs. (83) and (79) Update solution solving 2 × 2 linear system of Equations, e.g. Eq. (82) (j) 5. Solution: t+∆t εe = tr εe − ∆γ t+∆t n ˆ , t+∆t σ = σ t+∆t εe and t+∆t αi = t+∆t αi 6. If f¯a+1 > 0, a ←− a + 1 and go to step 4.
If solution needs ∆γ a < 0, a ←− a − 1 and go to step 4
7. Compute global tangent; see Section 3.3 Table 1: Local stress integration algorithm
24
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Taking into account that by Eq. (71)
d t+∆t n dt+∆t εd
= 2µ
2µ dt+∆t n ˆ I− = tr t+∆t d d ε k σ − t αk
t+∆t
n ˆ
t+∆t
(91)
n ˆ⊗
t+∆t
n ˆ
and that by Eq. (81) 0≡
dt+∆t f d∆γ = 2µ t+∆t n ˆ − 2µ t+∆t d − t+∆t d d ε d ε t+∆t ˜ t+∆t d H1 d r1 − t+∆t d ∆γ1 − t+∆t d d ε d ε
t+∆t
(92)
˜ 1 d∆γ1 H dt+∆t εd (93)
after some lengthy but straightforward algebra, we can factor-out the key derivative d∆γ1 = dt+∆t εd
2µ d∆γ + (2µ + t+∆t θ) d∆γ1
where t+∆t
θ=
r
2 3
t+∆t t+∆t H ˜
n ˆ
(94)
1
˜ 1 dt+∆t r1 dt+∆t H ∆γ1 t+∆t + t+∆t d p d p
!
(95)
Because at local convergence ∆γ a ≡ ∆γ, then d∆γ/d∆γ1 is given by Eq. (79). ˜ 1 = 2/3H ¯ and we Remarkably if only one surface is used we have d∆γ/d∆γ1 = 1, H naturally recover the usual expression for classical von Mises plasticity with mixed ¯ hardening if kinematic hardening is linear (dH/dp = 0). With the key Equation (94), all the quantities in the tangent Eq. (87) are known. The flow of the computational procedure for the tangent is Eq. (79) d∆γ Eq. (94) d∆γ1 −→ −→ t+∆t εp Eq. (89) Eq. (88) d d∆γ 1 dt+∆t εd −→ t+∆t d −→ t+∆t d n ˆ Eq. (92) d ε −→ t+∆t d d ε t+∆t ˜ 1 dt+∆t r1 Eq. (87) dt+∆t σd Eq. (85) p Eq. (90) dt+∆t H Eq. (89) d −→ t+∆t d −→ −→ t+∆t d −→ t+∆t C (96) , d ε dt+∆t εd dt+∆t εd d ε Obviously this diagram can also be employed to obtain the tangent in closed form. As it should be expected, it is straightforward to verify that if only one surface is 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
present, the classical tangent in Ref. [19] is recovered. 4. Large strains formulation As mentioned, the current formulation is a special case of classical von Mises plasticity with mixed hardening in which the effective kinematic hardening modulus is computed with the aid of internal variables. Therefore, as a remarkable advantage, the large strains implementation is immediate. To this end, we can use the algorithm of Eterovic and Bathe [65] based on logarithmic strains. This algorithm formulated in the full tensorial space is valid for kinematic hardening and uses a purely geometric pre- and post-processor, keeping the infinitesimal plastic stress integration unmodified. Using the multiplicative decomposition of the deformation gradient X = X e X p into elastic X e plastic X p parts, the elastic logarithmic strain in the intermediate configuration is obtained from the stretch tensor U e of the right polar decomposition X e = Re U e as E e = ln (U e )
(97)
The resulting incrementally additive framework using logarithmic strains mimics the infinitesimal theory, in which small elastic strains εe are interpreted as logarithmic elastic strains E e in the intermediate configuration, and Cauchy stresses σ as generalized Kirchhoff stresses T (or rotated Kirchhoff stresses in the elastic-isotropic case at hand). The hyperelastic stored energy is written in terms of these stress-strain measures as T = JU ′ (J) I + 2µE de (98) where JU ′ (J) and E de are respectively the volumetric and distorsional parts, the former obtained from the volumetric stored energy U ′ (J) as a function of the Jacobian determinant J. Algorithmically, the trial state is given by tr
Xe =
t+∆t t tX 0X e
(99)
and tr T d = 2µ tr E de with tr E de being the distorsional component of the trial elastic stretch tensor. Then, interpreting the backstress as a Kirchhoff-alike backstresses B, the small strains algorithm may be used without modification. Upon convergence, the final geometric update is performed from the converged elastic strains t+∆t t+∆t εe , as 0Ee ↔ tr t+∆t Re exp t+∆t0 E e (100) 0 Xe = 26
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
where we used the usual assumption that tr Re ≈ t+∆t0 Re [65], [63], [66], [59]. By systematic use of the chain rule, the tangent of the small strains algorithm may be converted to the typical tangent employed in finite element codes which relates second Piola-Kirchhoff stresses with Green-Lagrange strains in the reference configuration. Since the large strains pre- and post-processors are unchanged from the well-established classical plasticity and the experimental data used below are in the range of small strains, we refer for example to [65], [63] (Table 7.3.1), [66], [60], [59] for further details on the large strain numerical implementation using geometrical mappings. 5. Numerical examples In this section we include several numerical examples to demonstrate the predictive capabilities of the model and to demonstrate the performance of the fully implicit computational algorithm in actual finite element simulations. 5.1. Uniaxial example with cyclic isotropic hardening The purpose of this example is to show the capabilities to model cyclic hardening/softening during uniaxial cycles. The material for the simulation is steel BLY160. We used the material data from [57] and we predict the uniaxial cyclic experiments performed in Ref. [57]. In [57], the isotropic hardening is described with the relation between the normalized isotropic hardening parameter R/Rs and the accumulated plastic strain p, where Rs is the stabilized value. This function R/Rs is obtained in [57] by curve fitting from cyclic experiments shown in Fig. 18a of Ref. [57]. Since the saturated value Rs is also needed for our simulations, we take that data from Fig. 4 of Ref. [57], which contains experimental data of saturated stress values. Hence, from both figures in Ref. [57] we can obtain the needed relation R (p). However, we note that the shape of the skeleton curve is a saturated one with strain range 3% so the stress values have been adjusted by a ratio that is obtained by comparing the estimated yield stress of the saturated curve with the estimated yield stress of the first half loop curve. The resultant curve is the basic backbone which we used to obtain the initial kinematic hardening moduli field. With this data from Ref. [57], we performed simulations to predict the experimental observations shown in Fig. 3 of the same paper, which we reproduce in Figure 5a for the reader’s convenience. These experiments on three different specimens show the cyclic hardening behavior of the material at different deformation levels. In Figure 5b we show the predictions from our model. It can be seen that the predictions correlate quite well with the experimental observations. 27
a.
400
strain %
b.
300 200
stress (MPa)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
100 0 −100 −200 −300 −400 −0.04
−0.03
−0.02
−0.01
0
0.01
0.02
0.03
0.04
strain
Figure 5: (a) Experimental stress-strain cycles for prescribed strains of ∆ε = 1%, 2% and 3% for steel BLY160; adapted from [57], with permision form Elsevier. (b) Predictions of our model for the same experiments.
28
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
5.2. Examples under multiaxial, nonproportional loading. 5.2.1. Self-consistency of the formulation and recovery of classical bi-linear plasticity In order to show that classical von Mises plasticity with linear kinematic hardening is recovered in multiaxial cases regardless of the number of surfaces employed, we perform multiaxial simulations using as uniaxial stress-strain curve the bilinear one shown in Fig. 6a. In this figure we also show the locations of the different surfaces employed in the discretization and the slight modification to better appreciate the three cases in the simulations. We have prescribed the deformation path shown in Figure 6b. The result of the simulations using the three cases are shown in Figure 6c. It is seen that there is no relevant difference in the results. In all cases, the solution ˜ i /H ˜ i−1 is high enough to overcome the required is obtained in just one iteration if H tolerances in the local system of equations, except when a change of active surface is encountered. In such latter case, an additional iteration is needed just to update the position of the surface. 5.2.2. Experiments of Lamba and Sidebottom [61] We now show the predictions for the well-known nonproportional experiments of Lamba and Sidebottom [61] on oxygen-free high-conductivity (OFHC) copper. The nonproportional loading sequence (0 → 1 → · · · → 8) is shown in Figure 7. In this path, the positive and negative strain peaks are equal, and the angles between the different segments take very different values in order to test the response under change of direction, including a variety of angles [61]. To model these experiments we did not consider an isotropic hardening contribution. According to the loading sequence of Lamba and Sidebottom, the first segment 0 − 1 of the path is uniaxial shear loading, so the experimental stress-strain curve of this segment was used for generating the kinematic-hardening parameters of our model, see Eqs. (53) and (54), and the Poisson ratio to ν = 0.3. Figure 8 shows the comparison between our predictions (a) and the experimental results (b), the latter redrawn from Lamba and Sidebottom [61]. From the comparison of predicted and experimental results we can conclude that the model is capable of capturing most of the details of the loading history. Recall that the only data employed in the model is the proportional stress-strain curve for segment 0 − 1, which is shown in Figure 8 (1.a). A quantitative comparison of the von Mises stress in the characteristic points of the stress path is given in Table 2, where experimental values have been estimated from the plots in Ref. [61]. It is seen that the maximum differences are about 10%.
29
Uniaxial stress σ (GPa)
1
a
0.8 0.6
2 surfaces
0.4
4 surfaces
0.2
6 surfaces
0 0
0.2
0.4
0.6
0.8
Uniaxial strain ε (%) 1
1.5
1.2
1.4
1.6
c
b 1
Shear stress τ (GPa)
Shear strain γ/2 (%)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
0.5 0 −0.5 −1
0.5
0
−0.5 2 surfaces
−1
4 surfaces 6 surfaces
−1.5 −1
−0.5
0
0.5
−1.5
1
−1
−0.5
0
0.5
1
Axial stress σ (GPa)
Axial strain ε (%)
Figure 6: (a) Uniaxial stress-strain curve with different discretizations. Additional data: Poisson ˜ i = 103 H ˜ i−1 . (b) Prescribed strain path. (c) Predicted stress paths ratio ν = 0.3, H
30
Engineering shear strain γ %
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
1.5 1
6
1.0
5
0.5 0
0 −0.5
8 4 2
−1.0 −1.5 −0.4
7 3 −0.2
0
0.2
Axial strain ε %
0.4
a
Figure 7: Prescribed strain path in the tests from Lamba and Sidebottom [61]
Point 1 2 3 4 5 6 7 8
VM stress [MPa] (exp., approx.) 180.3 177.1 177.4 201.0 182.9 203.4 183.8 162.3
VM stress [MPa] (model) 185.5 185.1 166.1 178.5 164.5 179.3 198.8 164.8
% difference 2.88% 4.51% 6.37% 11.2% 10.1% 11.8% 8.16% 1.54%
Table 2: Approximate differences between experimental and computed von Mises stress at the turning points labelled in Fig. 7.
31
Shear Stress τ (MPa)
150 100
150 1(a)
5
50
4
0 −50
Axial Stress σa (MPa)
50
0
0
2
−150
100
100
−50
3 7
−100
1(b)
1 6
8
−100
Engrg.Shear Strain −0.01
0
γ
γ
Engrg.Shear Strain
0.01
−0.01
5 2
2(a)
0
−150
0.01
2(b)
100 0
0,1
0
7,8 3 6
−100
4
−200
100
3(a)
ε
0
0.004 −0.004 −0.002
0.002
Axial Strain
a
0
ε
−200
100
5
50
4
0
0
−50
−50
3
2
−100
7
−150
Axial Stress
−200
a
0.002 0.004
3(b)
1
8 6
50 0
−100
Axial Strain
−0.004 −0.002
Shear stress τ (MPa)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
−100
σ
a
0
−100 σ (MPa)
Axial Stress
(MPa) −200
100
−100
a
0
−150
100
Figure 8: Left column (a): Simulations of the Lamba-Sidebottom experiments using our model. The prescribed strain path is given in Figure 7. Right column (b): Experimental data from Lamba and Sidebottom, redrawn from [61]. The input data for the model is the monotonic curve of section 0-1
32
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
5.2.3. Experiments of Hamala et al [35] For further demonstration of the ability of the model for the description of multiaxial loading in other materials, we predicted the behavior of different specimens under four different strain-controlled tension-torsion loading paths given in Hamala et al [35], see also therein references. Then, the results of our simulations are compared with the experimental results in [35]. The uniaxial cyclic curve of the material (2124-T851 aluminum) given in Figure 8 of the same paper is used to determine our model, again using Eqs. (53) and (54). The four different multiaxial strain-loading paths (rhombus, circular, square and two-block) are simulated with a strain-driven procedure. The maximum amplitude in the prescribed √strain cycles shown in the small details in Fig. 9 is 1% (either axial ε or shear γ/ 3) in all cases. The comparison of the results from simulations and experiments along with prescribed strain paths are shown in Fig. 9. Since we did not consider isotropic hardening in this example, the first full cycle corresponds to a stabilized cycle in our model, and it is compared with the stabilized experimental results given in [35]. As shown in Fig. 9, the model predicts well the observed experimental stresspaths in all cases, in both the amplitudes and the different experimental shapes of the axial-shear stress√paths observed in the experiments. Since the plots given in the figure are in σ − 3τ axes, the distance to the origin gives a measure of the von Mises stress. The maximum differences in von Mises stress are found in the square and rhombus paths. In the corners of the square path, the maximum differences are about 15% at point labelled 2 in Fig. 9a, and about 10% at corner labelled 1. For the rhombus path, the maximum differences are found at corners labelled 1 and 3, with differences of about 10% and 18%, respectively. The maximum differences in the other two paths are about 10%. 5.2.4. Experiments of Tanaka et al [39] A series of experiments with nonproportional loadings were carried out also by Tanaka et al [36], [39]. Material parameters obtained from the experiments of Tanaka et al with loading paths of certain shapes have been used in different constitutive models for cyclic plasticity to derive a better description of the multiaxial behaviors of materials under those loading paths, see for example [39], [40], [41]. The prescribed loading path is given in Fig. 10a and the experimental result of the initial loading loop from [36] is redrawn in Fig 10b for comparison purposes. We note that the experimentally prescribed path given in [36], idealized in Fig. 10a, is given in terms of plastic strains, from machine-prescribed total strains with an iterative estimation therein done for the elastic part. We have performed our simulations using total strains, so an estimation of the maximum elastic contribution has been 33
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 9: Comparison of the predictions from our model against experimental cyclic data published in [35] for several prescribed multiaxial tension-torsion strain paths. (a) Square path. (b) Two-block path. (c) Rhombic path. (d) Circular path.
34
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
performed from stress values and the reported mean Young and shear modulus of the specimens, given in [36] as E = 203 MPa and G = 78 MPa. From Fig. 10b, it can be seen that the experimentally measured stresses in the initial strain-proportional path are not proportional, an effect also observed in other experiments in [36]. This can be attributed to experimental errors or to lack of compliance of the material state to the initial isotropic assumption. We show in Figure 10c the stress-strain curve employed in the simulations taken from Fig. 2 of Ref. [36], and in Fig. 10d the result of the simulations for the first cycle. Whereas there are several commented uncertainties that may affect the accuracy of the simulations, it is remarkable that the model captures well the shape of the stress path, including the special corner√effects in the path. The von Mises stress is the distance to the origin in the σ − 3τ plot. In Fig. 10d it is observed that the maximum difference is about 18% at the corner labelled 3, where isotropic hardening is already relevant. To include cyclic isotropic hardening in the simulations, from the peak stresses in the cycles in Fig 2 of Ref. [36] we have adjusted a function R(p), shown in Fig 10e. This curve, along with curve 10c was used for the simulations in Fig. 10d, which can be compared to those in Fig.3d of Ref. [36]. When including isotropic hardening, the error is about 10% at the mentioned corner. We note that we have not included path-dependent (directional) isotropic hardening. However, as commented, such effect could be included in the model allowing for a dependency r(γ, ξ). 5.2.5. Ramp-up loading experiments of Madrigal et al [71] In [71], a 90◦ tension-torsion out-of-phase experiment with continuously increasing strain amplitude is performed in AISI 4140 quenched and tempered steel, which is a material frequently used in the aerospace and automotive industries in elements that are subject to fatigue [71]. These experiments, along those from Lamba and Sidebottom, were used to test a previously proposed model [71], [72], [73], [74], similar to bounding surface models, but including a metric tensor which allows for the inclusion of different yield criteria, e.g. Tresca’s criterion [71]. This metric tensor also allows them to loosen the correlation between proportional axial and shear behaviors through additional material parameters in order to obtain better predictions for multiaxial loading. The 90o out-of-phase loading experiments with continuously increasing amplitude are challenging because they produce a large amount of additional multiaxial hardening and predictions are very sensitive to the uniaxial stress-strain curve and the size of the yield surface [71]. The experimental prescribed strain path of [71] is shown in Figure 11a. The experimental measured stress path obtained by Madrigal et al [71] is shown in Figure 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 10: Multiaxial cyclic experiments from [36]. (a) Prescribed strain path, 0 → 1 → 2 → 3 → 4. (b) Measured axial-shear stress path for the first cycle, redrawn from [36]. (c) Uniaxial monotonic stress-strain curve, obtained from Fig. 2 of [36] using a Young modulus of 203GP a and a shear modulus of 78GP a, which are average values given in Eq. (4) of [36]. (d) Predicted stress path for the first cycle. (e) Cyclic hardening curve R(p) − 1, obtained fitting approximate values of peak stress for subsequent cycles in Fig. 2 of [36]. (f) Predictions for several cycles including isotropic hardening (compare to Fig. 3d of [36]
36
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
11b. In order to perform the simulations, we have discretized the monotonic uniaxial stress-strain curve, given in Figure 5 of Ref. [71] and reproduced in Fig. 11c. We note that this material has a slightly different behavior in axial than in shear, which was the reason because of which Madrigal et al [71] developed a model including tensorial metrics. However, we performed the simulations using a von Mises criterion and taking the monotonic stress-strain curve for the model from axial loading; e.g. from Fig. 5 of Ref. [71]. In Figure 11d we show the results from our simulations for the ramp-up experiments with increasing amplitude. It is seen that the model reproduces the overall behavior, and quite accurately the axial stress values, having a difference precisely in the shear stress values, specially in the first loop. We attribute this difference in shear values to the commented slightly different behavior of the material in shear, which was the reason for the constitutive model proposal in Ref. [71]. 5.3. Finite element example One of the main purposes of multiaxial cyclic models is to compute stress redistribution during loading in a component subjected to multiaxial plastic cyclic loading. Then, an efficient finite element implementation of any proposed model, that allows for relatively large steps, still giving accurate predictions, is needed [5]. Accordingly, we have implemented the described implicit computational algorithm in our in-house finite element code Dulcinea, already used in simulations with other constitutive models [75], [76]. In order to verify the performance of the algorithm and, at the same time, to show a comparison with finite element simulations using the model of Ohno, we simulate a finite element example taken from the paper of Ohno et al [34], who used data from annealed OFHC copper. Ohno et al used this example to show the performance of the stress integration algorithm presented therein for the well-known model of Ohno. The simulation is performed on the finite element model of half of a singlehole plate fixed at one end and subjected to cyclic bending and a constant tension at the opposite end. The cyclic bending value uend y , which is varied from 0 mm to end 5 mm, is added after the axial tension σxx is increased from 0 MPa to 5 MPa. All the displacements are imposed in our model by means of a penalty formulation. The axial stresses are imposed by equivalent uniformly distributed loads at the nodes. Details of the example are given in Figure 12. In our simulations we obtain the hardening surfaces from the stress-strain curve shown in Figure 2 of Ref. [34] following the same procedure as in the previous examples. As in the simulations performed by Ohno, we include in this example 37
0.02 1000 Shear stress τ√3 (MPa)
Shear strain γ/√3
a 0.01
0
−0.01
b
500 0 −500
−1000 −0.02 −0.02
−0.01
0
0.01
−1000 −500
0.02
Axial strain ε
1000 Shear stress τ√3 (MPa)
1000
750
500
250
0
500
1000
Axial stress σ (MPa)
c Axial stress σ (MPa)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
d
500 0 −500
−1000 0 0
0.005 0.01 Axial strain ε
0.015
−1000 −500
0
Axial stress σ
500
1000
(MPa)
Figure 11: Experimental and simulated ramp-up 90◦ out-of-phase strain path with continuously increasing amplitude. (a) Prescribed experimental axial-shear strain path [71]. (b) Experimental axial-shear stress path [71]. (c) Axial monotonic stress-strain curve employed in the simulation, obtained from Fig. 5 of [71]. (d) Predicted stress path
38
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
a 100 Ø25
75
25
end xx
s
z
x
3.5
end
uy
y Time (step)
b end
uy
0 0.7
5
15
25
35
45
55
2.1 3.5
Figure 12: Plate with a hole, clamped at one end and subjected to an axial constant tension and a cyclic bending displacement at the other end; example taken from Ref. [34]. (a) Details of the geometry with dimensions in mm. (b) Prescribed cyclic bending displacement applied to the free end.
39
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
cyclic isotropic hardening obtained from Table 4 of Ref. [34], so all the presented capabilities are active and the efficiency of the consistent tangents may be verified in the general case. In Ref. [34] the authors made the implementation of their model in the commercial code Abaqus via user subroutine. The elements they used are the typical 20-node serendipity elements C3D20R with reduced integration to avoid volumetric locking. However, in our case, in order to include more integration points per element (important in stress concentration cases as the one at hand) and to use a mixed u/p formulation, we selected fully integrated 27/4 tri-quadratic elements with linear pressure interpolation. These are reliable elements which are known to pass the inf-sup condition and give optimal performance [77]. The finite element mesh is shown in Figure 13. We included a finer layer of elements as to compare our integration-point stress-strain results with the nodal ones in Ref. [34]. Since both the specimen and the loading pattern are symmetric, we used symmetry conditions in the finite element model, so only one half of the specimen has been modelled. In our simulations we have used the same step increment for cyclic bending as the one used in Ref. [34]. Even though this time step is too large for an accurate description of the behavior of the structure, this step size shows the robustness of the numerical implementation and facilitates the comparison with [34]. However, for the initial tension loading, to better describe the incursion from the elastic into the plastic regime and to show the effect in the convergence rate, we used 5 steps instead of the single step used by Ohno et al [34]. In Figure 13 we show the von Mises stress distribution in the upper face or the specimen (in MPa) at the end of the loading cycle. Note that the band plot results shown in the figure are unaveraged, nonsmoothed stresses extrapolated at nodes in order to show the quality of the mesh. The stress jumps at both ends of the specimen are due to the approximate approach used in prescribing the boundary conditions, and for our purposes those jumps are irrelevant in the simulation according to the Saint Venant principle. The stress-strain cyclic results shown in Figure 13c are obtained at the integration point closest to the middle point of the half circle in the direction of width and closest to the upper surface in thickness. As shown in Figure 13c, our predictions of the stress-strain hysteresis loops show the feature of cyclic hardening. These predictions are very similar to those given by Ohno et al in Figure 8b of Ref. [34] using their model. Apart from being different models (although both use a multilinear approximation obtained from stress-strain curves), the small discrepancies may also be attributed to the ratchetting effect included in Ref. [34], the different element type used, the finer mesh we have used and the actual location of the corresponding integration points, which is specially important in that zone due to the stress concentration promoted by the hole. 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 13: Finite element mesh used for our simulations of the example given in [34] (a); elements are 27/4 u/p mixed elements (tri-quadratic interpolation in displacements and linear in pressure) with full 3 × 3 × 3 integration. Band plots correspond to nodal-unaveraged von-Mises stresses at the top layer (band-plot jumps are an estimate of the local error in the solution) (b). It is also indicated the location of node nearest to the integration point which corresponds to the stress-strain history plot (c)
41
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Iteration (1) (2) (3) (4)
Step 14 4.366E + 06 4.706E + 02 1.711E + 01 1.973E − 01
Iteration (1) (2) (3) (4)
Step 14 3.063E + 07 6.833E − 01 1.713E − 03 2.282E − 07
Residual norm Step 34 4.366E + 06 4.270E + 02 1.539E + 01 1.021E − 01 Energy norm Step 34 3.063E + 07 6.133E − 01 1.302E − 03 8.144E − 08
Step 54 4.366E + 06 4.839E + 02 1.663E + 01 2.542E − 01 Step 54 3.063E + 07 6.889E − 01 1.125E − 03 1.213E − 07
Table 3: Global equilibrium convergence values at steps 14, 34, 54 of the finite element simulation
The number of global equilibrium iterations employed for each step is shown in Figure 14. As usual, the first iteration corresponds to the imposition of the incremental loads of the step, giving the reference non-equilibrium values, whereas the remaining iterations are the actual iterative search for the equilibrium configuration. Since the first two steps are elastic, just two iterations (one for imposing load and one to regain equilibrium) give the solution to machine precision. The rest of the steps involve plastic strains, so more iterations are needed to reach the prescribed relative tolerances in forces of 10−7 and in energy of 10−14 . However note that usually only 4 iterations (3 equilibrium iterations) are needed to reach such tolerances, which is a performance similar to that of classical von Mises plasticity. Furthermore, in Table 3 we show force and energy convergence values for three typical steps. It can be observed that asymptotic quadratic convergence is obtained. 6. Conclusion In this paper we have presented a new model for multiaxial cyclic plasticity which follows the classical Prager’s formulation, including nonlinear kinematic hardening and cyclic hardening/softening. We also developed a fully implicit finite element implementation of the model. The model uses the idea from Mr´oz of discretizing the uniaxial stress-strain curve in several segments, resulting in a field of hardening surfaces. This approach is very simple in obtaining material parameters of the model. However, in contrast to Mr´oz’s proposal, our model is based on the traditional associative Prager rule and the yield surface is always the innermost surface. The 42
Number of global iterations per step
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
7 6 5 4 3 2 1 0
10
20
30
40
50
step number
Figure 14: Number of global iterations used for each loading step
outer surfaces are just hardening surfaces, a tool to describe directional and historydependent anisotropic hardening moduli. Our formulation and numerical algorithm naturally recover classical von Mises plasticity with mixed hardening in the case that kinematic hardening is linear in practice, regardless of the number of surfaces being employed. Hence, the predictions are self-consistent. Furthermore, since the model is in essence classical plasticity, the large strain implementation is straightforward and unmodified from the classical plasticity using the typical geometric pre- and post-processors in terms of logarithmic strains. We have validated the model with predictions of experimental results from both uniaxial and multiaxial tests taken from the literature. It is seen that the model is capable of giving good predictions in all cases, which include different materials and multiaxial loading paths. Finally, we have demonstrated the numerical efficiency of the fully implicit, Backward-Euler algorithmic formulation by the implementation in a finite element program and the simulation of a cyclic nonproportional loading example from the literature. Two effects are still not included in the formulation herein presented: the ratchetting effect and the path-dependent isotropic hardening. These effects can be accommodated in the present theory preserving Prager’s rule, respectively, by the inclusion of a ratchetting internal strain in the internal compatibility equation, and by allowing the full dependency of the isotropic hardening on ξ.
43
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Acknowledgments Partial funding from grant DPI2015-69801-R from the Direcci´on General de Proyectos de Investigaci´on of the Ministerio de Econom´ıa y Competitividad of Spain is acknowledged. We thank the authors of Ref. [71] for the access to the digital data of the experiment in their paper, shown in Fig. 11. FJM also acknowledges the support of the Department of Mechanical and Aerospace Engineering of University of Florida during the sabbatical period in which part of this work was done, and Ministerio de Educaci´on Cultura y Deporte of Spain for the financial support for that stay under grant PRX15/00065. References [1] Borja RI, Amies AP. Multiaxial cyclic plasticity model for clays. J Geotech Eng ASCE 1994; 120(6):1051-1070. 1 [2] Darwin D, Nmai CK. Energy dissipation in RC beams under cyclic load. J Struct Eng 1986; 112(8):1829-1846. 1 [3] Benavent-Climent, A. Seismic behavior of RC wide beam-column connections under dynamic loading. J Earthquake Eng 2007; 11(4):493-511. 1 [4] Gates NR, Fatemi A. A simplified cyclic plasticity model for calculating stressstrain response under multiaxial non-proportional loadings. Eur J Mech-A/Solid 2016; 59:344-355. 1 [5] Madrigal C, Navarro A, Chaves V. Numerical implementation of a multiaxial cyclic plasticity model for the local strain method in low cycle fatigue. Theor Appl Fract Mech 2015; 80:111-119. 1, 5.3 [6] Navarro A, Vallellano C, Chaves V, Madrigal C. A microstructural model for biaxial fatigue conditions. Int J Fatigue 2011;33(8):1048-1054. 1 [7] Teixeira P, Santos AD, Andrade Pires FM, C´esar de S´a JMA. Finite element prediction of ductile fracture in sheet metal forming processes. J Mater Proc Technol 2006;177:278-281. 1 [8] Khoei A. Computational plasticity in powder forming processes. Elsevier; 2010. 1
44
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[9] Cvitani´c V, Vlak F, Lozina Z. A finite element formulation based on nonassociated plasticity for sheet metal forming. Int J Plasticity 2008; 24(4):646687. 1 [10] Vladimirov IN, Pietryga MP, Reese S. Prediction of springback in sheet forming by a new finite strain model with nonlinear kinematic and isotropic hardening. J Mater Proc Technol 2009; 209:4062-4075 1 [11] Levkovitch V, Svendsen B. Accurate hardening modeling as basis for the realistic simulation of sheet forming processes with complex strain-path changes. AIP Conf Proc 2007; 908:1331-1336 1 [12] Milne I, Ritchie RO, Karihaloo BL (Eds). Comprehensive structural integrity: Cyclic loading and fatigue (vol 4). Elsevier; 2003. 1 [13] Li B, Reis L, De Freitas M. Simulation of cyclic stress/strain evolutions for multiaxial fatigue life prediction. Int J Fatigue, 2006;28(5):451-458. 1 [14] Ince A, Glinka G, Buczynski A. Computational modeling of multiaxial elastoplastic stress–strain response for notched components under non-proportional loading. Int J Fatigue 2014; 62:42-52. 1 [15] Rezaiee-Pajand M, Sharifian M. Accurate and approximate integrations of Drucker-Prager plasticity with linear isotropic and kinematic hardening. Eur J Mech A/Solid 2011; 30:345-361 1 [16] Moftakhar A, Buczynski A, Glinka G. Calculation of elasto-plastic strains and stresses in notches under multiaxial loading. Int J Fract 1994; 70(4):357-373. 1 [17] Zeng Z, Fatemi A. Elasto-plastic stress and strain behavior at notch roots under monotonic and cyclic loadings. J Str Anal Eng Des 2001; 36(3):287-300. 1 [18] Singh MNK, Glinka G, Dubey RN. Elastic-plastic stress-strain calculation in notched bodies subjected to non-proportional loading. Int J Fract 1986; 76(1):39-60. 1 [19] Simo JC, Hughes TJR. Computational Inelasticity. Springer; 1998. 1, 2.1, 2.1, 3.2, 3.3 [20] Lemaitre J, Chaboche JL. Mechanics of solid materials. Cambridge University Press; 1994. 1
45
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[21] De Angelis F, Taylor RL. An efficient return mapping algorithm for elastoplasticity with exact closed form solution of the local constitutive problem. Eng Comput 2015; 32(8): 2259-2291 1 [22] De Angelis F, Taylor RL. A nonlinear finite element plasticity formulation without matrix inversions. Fin Elem Anal Des 2016; 112:11-25. 1 [23] Valanis KC. Fundamental consequences of a new intrinsic time measure. Plasticity as a limit of the endochronic theory. Arch Mech 1978; 32:171. 1, 2.2 [24] Mr´oz Z. On the description of anisotropic workhardening. J Mech Phys Solid 1967; 15(3): 163-175. 1 [25] Garud YS. A new approach to the evaluation of fatigue under multiaxial loadings. J Eng Mater Technol ASME 1981; 103(2):118-125. 1 [26] Armstrong PJ, Frederick CO. A mathematical representation of the multiaxial Bauschinger effect. GEGB Report RD/B/N/731, Berkeley Nuclear Laboratories, UK; 1966. 1 [27] Chaboche JL. Time-independent constitutive theories for cyclic plasticity. Int J Plasticity 1986; 2(2):149-188 1, 2.5 [28] Dafalias YF, Popov EP. Plastic internal variables formalism of cyclic plasticity. J Appl Mech ASME 1976; 43: 645. 1 [29] Dafalias YF. The Concept and Application of the Bounding Surface in Plasticity Theory. In: Hult J, Lemaitre J (eds) Physical Non-Linearities in Structural Analysis. International Union of Theoretical and Applied Mechanics. Springer; 1981. 1 [30] Ohno N, Wang JD. Kinematic hardening rules with critical state of dynamic recovery. Part I: formulation and basic features for ratchetting behavior. Int J Plasticity 1993; 9:375–390. 1 [31] Mont´ans FJ. Implicit algorithms for multilayer J 2-plasticity. Comput Meth Appl Mech Eng 2000; 189(2):673-700. 1 [32] Caminero MA, Mont´ans FJ. An enhanced algorithm for nested surfaces plasticity using the implicit Mr´oz translation rule. Comput Struct 2006; 84(26):16841695. 1
46
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[33] Khoei AR, Jamali N. On the implementation of a multi-surface kinematic hardening plasticity and its applications. Int J Plasticity 2005; 21(9):1741-1770. 1 [34] Ohno N, Tsuda M, Kamei T. Elatoplastic implicit integration algorithm applicable to both plane stress and three-dimensional stress states. Finite Elem Anal Des 2013; 66:1-11 1, 5.3, 5.3, 12, 5.3, 13 [35] Halama R, Markopoulos A, Janco R, Bartecky M. Implementation of MAKOC cyclic plsticity model with memory. Adv Eng Soft 2017; in press, doi: 10.1016/j.advengsoft.2016.10.009 1, 5.2.3, 9 [36] Tanaka E, Murakami S, Ooka M. Effects of strain path shapes on nonproportional cyclic plasticity. J Mech Phys Solid 1985; 33(6):559-575 1, 2.1, 5.2.4, 5.2.4, 10 [37] Benallal A, LeGallo P, Marquis D. An experimental investigation of cyclic hardening of 316 stainless steel and of 2024 aluminum alloy under multiaxial loadings. Nucl Eng Des 1989; 114:345-353 1 [38] Calloch S , Marquis D. Triaxial tension–compression tests for multiaxial cyclic plasticity. Int J Plasticity. 1999; 5(5):521–49. 1 [39] Tanaka E. A nonproportionality parameter and a cyclic viscoplastic constitutive model taking into account amplitude dependences and memory effects of isotropic hardening. Eur J Mech/A Solid 1994; 13(2):155-173 1, 5.2.4 [40] Gates NR, Fatemi A. A simplified cyclic plasticity model for calculating stressstrain response under multiaxial non-proportional loadings. Eur J Mech A/Solid 2016; 59:344-355 1, 5.2.4 [41] Meggiolaro MA, Wu H, Pinho de Castro, JT. Non-proportional hardening models for predicting mean and peak stress evolution in multiaxial fatigue using Tanaka’s incremental plasticity concepts. Int J Fatigue 2016; 82:146-157 1, 5.2.4 [42] Taleb L, Cailletaud G. An updated version of the multimechanism model for cyclic plasticity. Int J Plasticity 2010; 26:859-874 1 [43] Cardoso RPR, Yoon JW. Stress integration method for a nonlinear kinematic/isotropic hardening model and its characterization based on polycrystal plasticity. Int J Plasticity 2009; 25:1684-1710 1 47
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[44] Mont´ans FJ. Implicit multilayer J2-plasticity using Prager’s translation rule. Int J Num Meth Eng 2001; 50:347-375 1 [45] Mont´ans FJ. Implicit plane stress algorithm for multilayer J2-plasticity using the Prager–Ziegler translation rule. Int J Num Meth Eng 2004; 59(3):409-418. 1 [46] Jiang Y, Sehitoglu H. Comments on the Mr´oz multiple surface type plasticity models. Int J Solid Struct 1996; 33(7):1053-1068. 1 [47] Mont´ans FJ, Caminero MA. On the consistency of nested surfaces models and their kinematic hardening rules. Int J Solid Struct 2007; 44(14):5027-5042. 1 [48] Chaboche JL. On some modifications of kinematic hardening to improve the description of ratchetting effects. Int J Plasticity 1991; 7:661-678. 1 [49] Abdel-Karim M. An extension for the Ohno-Wang kinematic hardening rules to incorporate isotropic hardening. Int J Pres Vessel Piping 2000; 87(4):170-176. 1 [50] Kossa A, Szabo L. Exact integration of the Von Mises elastoplasticity model with combined linear isotropic-kinematic hardening. Int J Plasticity 2009; 25:10831106 1 [51] Berisha B, Hora P, Wahlen A, Tong L. A combined isotropic-kinematic hardening model for the simulation of warm forming and subsequent loading at room temperature. Int J Plasticity 2009; 26:126-140 1 [52] Cao J, Lee W, Cheng HS, Seniw M, Wang HP, Chung K. Experimental and numerical investigation of combined isotropic-kinematic hardening behavior of sheet metals. Int J Plasticity 2009; 25:942-972 1 [53] Firat M. Cyclic plasticity modeling and finite element analyzes of a circumferentially notched round bar under combined axial and torsion loadings. Mater Design 2012; 34:842-852 1 [54] Hong PY, Pereira JM, Cui YJ, Tang AM, Collin F, Li XL. An elastoplastic model with combined isotropic-kinematic hardening to predict the cyclic behavior of stiff clays. Comput Geotech 2014; 62:193-202 1 [55] Aahnken R. Improved implementation of an algorithm for non-linear isotropic/kinematic hardening in elstoplasticity. Commun Num Meth Eng 1999; 15:745-754 1 48
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[56] Ohno N. A Constitutive model of cyclic plasticity with a non-hardening strain region. J Appl Mech. 1982; 49:721-727 1 [57] Xu L, Nie X, Fan J, Tao M, Ding R. Cyclic hardening and softening behavior of the low yield point steel BLY160: Experimental response and constitutive modeling. Int J Plasticity 2016; 78:44-63 1, 2.5, 5.1, 5 [58] Krieg RD, Key SW. Implementation of a time dependent plasticity theory into structural computer programs. In Constitutive Equations in Viscoplasticity: Computational and Engineering Aspects, JA Stricklin and KJ Saczalski eds, Appl Mech Div AMD-20, ASME, New York 1976. 1, 3.2 [59] Sanz MA, Mont´ans FJ, Latorre M. Computational anisotropic hardening multiplicative elstoplasticity based on the corrector elastic logarithmic strain rate. Comput Meth Appl Mech Eng 2017; 320:82-121 2.1, 2.1, 4 [60] Latorre M, Mont´ans FJ. A new class of plastic flow evolution equations for anisotropic multiplicative elastoplasticity based on the notion of a corrector elastic strain rate. Under review. arXiv:1701.00095 [cond-mat.soft] 2016 2.1, 2.1, 4 [61] Lamba HS, Sidebottom OM. Cyclic plasticity for nonproportional paths: Part 1 — Cyclic hardening, erasure of memory, and subsequent strain hardening experiments. Part 2 — Comparison with predictions of three incremental plasticity models. J Eng Mater Technol 1978; 100:93-111 5.2.2, 7, 8 [62] Zhang M, Ben´ıtez JM, Mont´ans. Capturing yield surface evolution with a multilinear anisotropic kinematic hardening model. Int J Solid Struct 2016; 81:329– 336 [63] Koji´c M, Bathe KJ. Inelastic Analysis of Solids and Structures. Springer; 2005. 2.1, 4 [64] Borja RI. Plasticity. Springer; 2013. 2.1 [65] Eterovic AL, Bathe KJ. A hyperelastic-based large strain elasto-plastic constitutive formulation with combined isotropic-kinematic hardening using the logarithmic stress and strain measures. Int J Num Meth Engrg 1990; 30(6):10991114. 4, 4
49
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[66] Caminero MA, Mont´ans FJ, Bathe KJ. Modeling large strain anisotropic elastoplasticity with logarithmic strain and stress measures. Comput Struct 2011; 89(11): 826-843. 4 [67] Valanis KC, Wu HC. Endochronic representation of cyclic creep and relaxation of metals. J Appl Mech ASME 1975; 42:67-73. 2.2 [68] Mont´ans FJ, Caminero MA. On the consistency of nested surfaces models and their kinematic hardening rules. Int J Solid Struct 2007; 44(14): 5027-5042. 2.4 [69] Mi˜ nano M, Caminero MA, Mont´ans FJ. On the numerical implementation of the Closest Point Projection algorithm in anisotropic elasto-plasticity with nonlinear mixed hardening. Fin Elem Anal Des 2016; 121:1-17. 3 [70] Mont´ans FJ, Borja RI. Implicit J2-bounding surface plasticity using Prager’s translation rule. Int. J Num Meth Eng 2002; 55:1129:1166. [71] Madrigal C, Navarro A, Chaves V. Biaxial cyclic plasticity experiments and application of a constitutive model for cyclically stable material behaviour. Int J Fatigue 2016; 83(2):240-252 5.2.5, 11, 6 [72] Navarro A, Brown MW. A constitutive model for elastic-plastic deofrmation under cyclic multiaxial straining. Fatigue Fract Eng Mater Struct 1997; 20(5):747758 5.2.5 [73] Navarro A, Gir´aldez JM, Vallellano C. A constitutive model for elastoplastic deformation under variable amplitude multiaxial cyclic loading. Int J Fatigue 2005; 27:838-846 5.2.5 [74] Navarro A, Brown MW, Miller KJ. A multiaxial stress-strain analysis for proportional cyclic loading. J Strain Anal 1993; 28(2):125-133 5.2.5 [75] Latorre M, Mont´ans FJ. Fully anisotropic finite strain viscoelasticity based on a reverse multiplicative decomposition and logarithmic strains. Comput Struct 2016; 163:56-70 5.3 [76] Crespo J, Latorre M, Mont´ans FJ. WYPIWYG hiperelasticity for isotropic, compressible materials. Comput Mech 2017; 59:73-92. 5.3 [77] Bathe KJ. Finite Element Procedures, 2nd Ed. Klaus-J¨ urgen Bathe; 2014. 5.3
50
PII: DOI: Reference:
S0045-7825(17)30655-2 https://doi.org/10.1016/j.cma.2017.09.028 CMA 11622
To appear in:
Comput. Methods Appl. Mech. Engrg.
Received date : 24 May 2017 Revised date : 5 August 2017 Accepted date : 18 September 2017 Please cite this article as: M. Zhang, J.M. Ben´ıtez, F.J. Mont´ans, Cyclic plasticity using Prager’s translation rule and both nonlinear kinematic and isotropic hardening: Theory, validation and algorithmic implementation, Comput. Methods Appl. Mech. Engrg. (2017), https://doi.org/10.1016/j.cma.2017.09.028 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
*Manuscript Click here to download Manuscript: Zhang_Benitez_Montans_CMAME1_R1_vf.pdf
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Click here to view linked References
Cyclic plasticity using Prager’s translation rule and both nonlinear kinematic and isotropic hardening: theory, validation and algorithmic implementation Meijuan Zhanga , Jos´e Mar´ıa Ben´ıteza , Francisco J. Mont´ansa,∗ a
Escuela T´ecnica Superior de Ingenier´ıa Aeron´ autica y del Espacio Universidad Polit´ecnica de Madrid Plaza Cardenal Cisneros, 3, 28040-Madrid, Spain
Abstract Finite element analysis of structures under elasto-plastic nonproportional cyclic loadings is useful in seismic engineering, fatigue analysis and ductile fracture. Usual models with nonlinear stress-strain curves in cyclic behavior are based on Mr´oz multisurface plasticity, bounding surface models or models derived from the ArmstrongFrederick rule. These models depart from the associative Prager’s rule with the main purpose of modeling aspects of cyclic nonlinear hardening. In this paper we develop a model for cyclic plasticity within the framework of the associative classical plasticity theory using Prager’s rule accounting for anisotropic nonlinear kinematic hardening coupled with nonlinear isotropic hardening. We include the validation of the theory against several uniaxial and multiaxial cyclic experiments and an efficient fully implicit radial return algorithm. The parameters of the model are obtained directly by a discretization of the uniaxial stress-strain behavior. Remarkably, both the presented theory and the computational algorithm automatically recover classical bilinear plasticity and the Krieg and Key algorithm if the user-prescribed stress-strain curve is bilinear. Keywords: Cyclic plasticity, Associative plasticity, Radial return algorithm, Multisurface model, Nonlinear kinematic hardening
∗
Corresponding author Email addresses: [email protected] (Meijuan Zhang), [email protected] (Jos´e Mar´ıa Ben´ıtez), [email protected] (Francisco J. Mont´ ans )
Preprint submitted to Computer Methods in Applied Mechanics and Engineering
August 5, 2017
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
1. Introduction The finite element analysis of structures under nonproportional and cyclic nonlinear behavior is very important in many applications, among them seismic engineering [1], [2], [3], fatigue and fracture analyses [4], [5], [6], [7] and plastic forming [8], [9], [10], [11]. The correct description of multiaxial hardening effects have proved to be critical for an accurate prediction of the effective displacements, accelerations and safety of the structures [12]. In fatigue analysis of notched specimens, the plastic loading at the notch edge induces nonproportional multiaxial loading [13], frequently studied through cyclic plasticity models which, thereafter, are related to nominal strains and stresses through incremental Neuber rules or related strain or energy methods [5], [14], [15], [16], [17], [18]. In cyclic behavior, the accurate representation of the Bauschinger effect is crucial, a behavior which is approximately represented by Masing’s rules. The classical linear kinematic hardening reproduces accurately such rules, has a very efficient integration algorithm and is available in most finite element programs. This model uses Prager’s translation rule for the backstress, which is consistent with the principle of maximum dissipation [19]. However, a priori, the direct inclusion of a nonlinear anisotropic kinematic hardening function produces unphysical loops [20]. Therefore, finite element programs use Prager’s rule exclusively in the case of linear kinematic hardening, and for nonlinear kinematic hardening, which is needed to better describe the cyclic loops, they resort to other types of formulations (if available to the user), typically based on the Armstrong and Frederick rule [21], [22]. One of the procedures to account for anisotropic, history-dependent nonlinear kinematic hardening is based on the history of observable quantities by hereditary integrals, as in the endochronic theory, similar to those used in viscoelasticity [23]. However, since this approach is more cumbersome for computational application and for finite element analysis, the approach based on internal variables [20], where only information of the previous step is needed, is usually preferred. Remarkably, the theories to model nonlinear kinematic hardening based on internal variables have departed from the classical linearly hardened framework, introducing different non-associative translation rules for the backstress, as Mr´oz’s rule [24], Garud’s rule [25] or the commented Armstrong and Frederick rule [26]. These rules have resulted in different models as the Mr´oz model [24], Chaboche’s model [20], [27], bounding surface models [28], [29], [1], and nonlinear kinematic hardening models with the addition of multiple backstresses [20], [30]. The algorithmic implementation of some of these models is usually more elaborate than that of classical plasticity [31], [32], [33], [34], [35], and even though some recent efficient algorithms are available for some cases [21], [22], they do not constitute a natural extension of 2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
classical Prager’s plasticity. To improve the cyclic multiaxial plastic behavior, sometimes non-proportionality parameters, obtained for certain materials under certain loading paths by fitting experiments, are employed, like in Refs. [36], [37], [38], [39], [40], [41], [42], [35], [43]. However, despite the common belief [20, p. 206], it is possible to develop nonlinear kinematic hardening models using Prager’s translation rule and preserving Masing’s rules [44], [45]. These models are similar in conception to Mroz’s model, but with some crucial diferences, as the preservation of Prager’s rule, the use of a single yield surface (outer surfaces are hardening surfaces) and the simplicity of the integration algorithms. Furthermore, they lack the inconsistencies present in Mr´oz’s model under multiaxial loading [46], [47]. As we show in this paper, the predictions for multiaxial plastic behavior using Prager’s rule are similar to the behavior observed in the experiments, at least for the different materials addressed below. Despite of the employed kinematic hardening rule, mixed kinematic-isotropic hardening models have been increasingly used in cyclic loading analysis [48], [49], [50], [51], [52], [53], [54], [55], [10]. For a better description of cyclic hardening/softening in certain materials, the effect of a memory for the strain amplitude can be included in the isotropic hardening rule [27], [56], [57], [42], and the inclusion of isotropic cyclic softening may be used to model the effects of damage by fatigue. However, as seen below, the presence of isotropic hardening also modifies the anisotropic kinematic hardening moduli, an observation that must be considered in the formulation and computational algorithm. Therefore, it is valuable, and the purpose of this paper, to extend the classical von Mises associative theory of plasticity using Prager’s translation rule to account for history-dependent nonlinear anisotropic kinematic hardening combined with cyclic isotropic hardening/softening. In the next sections we first introduce the theory, showing that we just propose a simple extension of the classical framework by allowing a dependence of the effective hardening moduli on some internal variables. Then we introduce the method to compute such dependence, i.e. the effective hardening modulus at each instant. Thereafter we explain the efficient radial return algorithm, which for the bilinear case naturally recovers the solution from Krieg and Key [58]. With some examples we show that the theory predicts rather well the experimentally observed multiaxial behavior in several materials and that the finite element implementation preserves the asymptotic quadratic convergence of Newton schemes.
3
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
2. Classical J2 –plasticity using Prager’s rule with nonlinear anisotropic mixed hardening 2.1. Classical associative plasticity with anisotropic, cyclically modified hardening The stored “observable” (elastic) energy depends on the elastic strains εe , which are a function of the total strains and some internal strains εp (the plastic strains), i.e. εe (ε, εp ). The chain rule derived from the explicit dependencies εe (ε, εp ), gives immediately the rate as ∂εe ∂εe : ε˙ + : ε˙p = ε˙e |ε˙p =0 + ε˙e |ε=0 (1) ε˙e = ˙ ∂ε ε˙p =0 ∂εp ε=0 ˙ ≡
tr
ε˙e +
ct
ε˙e ≡ ε˙ − ε˙p
(2)
where the last identity accomodates the rate form of the usual additive decomposition ansatz εe = ε − εp . We defined the trial and corrector partial derivative contributions tr ε˙e ≡ ε˙e |ε˙p =0 and ct ε˙e ≡ ε˙e |ε=0 ˙ . In the new formulation of large strain elasto-plasticity presented in [59] and [60], this framework based on the chain rule remains additive and unaltered at large strains even when using the multiplicative decomposition. The stored energy is —see Figure 1 Ψ (εe , ξ, γ) = W (εe ) + H (ξ, γ) = 12 εe : C : εe + H (ξ, γ)
(3)
where W (εe ) is the elastic energy, immediately recovered upon unloading, and H (ξ, γ) is the internal, microstructurally-blocked “elastic” energy of which part can be recovered following certain plastic deformation paths. This energy may be function of a scalar variable γ (accounting for the cummulative plastic strain) and a internal variable ξ of tensor nature; and in turn ξ itself may be a function of some internal variables ξi of microscopic origin, mapped to the continuum level. For elastically isotropic materials, and using the usual observation that plastic flow is isochoric, we write (4) Ψ (εe , ξ, γ) = 21 κ (εv )2 + 21 2µεde : εde + H (ξ, γ)
where εv are the volumetric strains (elastic), εde are the elastic deviatoric strains, µ is the shear modulus, κ is the bulk modulus. We define the backstress α and the overstress rˆ as α (ξ, γ) :=
∂H (ξ, γ) ∂H (ξ, γ) , and rˆ (ξ, γ) := ∂ξ ∂γ
With these definitions, the rate of the stored energy is 4
(5)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 1: Decomposition of the input energy
˙ = σ : ε˙e + α : ξ˙ + rˆγ˙ = κεv ε˙v + 2µεd : ε˙d + α : ξ˙ + rˆγ˙ Ψ e e
(6)
The rate ξ˙ = − ct ε˙e = ε˙p , accounting for the evolution of internal variables, only takes place during plastic flow. The definition ξ˙ = − ct ε˙e is just a traditional, convenient definition so α results in the usual macroscopic, continuum backstress, that can be directly subtracted from the stress in the yield function as seen below. In ˙ a similar we take γ˙ as the norm of ξ,
definitions,
fashion, to accommodate the usual
˙ ˆ and we define n ˆ = ξ/ γ˙ = ξ˙ , so ξ˙ = γ˙ n
ξ˙ , and the last term in Eq. (6) adopts the familiar form for isotropic hardening. Note also that sometimes the sign of both ξ and α are opposite than in our definition here to understand ξ as internal elastic strains mapped to the continnum (because of the “elastic” nature of H (ξ, γ)), but in practice the numerical result remains to be the same. Denoting the stress tensor by σ, the stress power is P = σ : ε. ˙ Using ε˙ ≡ tr ε˙e , the dissipation is ˙ =σ: D˙ = P − Ψ
ε˙e − κ (εv ) ε˙v − 2µεde : ε˙de − α : ξ˙ − rˆγ˙ ≥ 0 (7) ˙ = tr Ψ− ˙ ˙ ≡ − ct Ψ which can be interpreted as D˙ ≡ − ∂Ψ (ε, εp ) /∂εp |ε=0 : ε ˙ = − Ψ p ˙ ε=0 ˙ ˙ Following the usual Coleman procedure, for purely elastic incremental deformaΨ. tions, when D˙ = 0, we also have ε˙p = 0 and by definition ct ε˙e = −ξ˙ = 0. Thus by chain differentiation we immediately identify the necessary condition σ = σv + σd = tr
5
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Kεv I + 2µεde , where σv and σd are the volumetric and deviatoric parts, respectively. Then, since γ˙ = (dγ/dξ) : ξ˙ = n ˆ : ξ˙ = − n ˆ : ct ε˙e and tr ε˙e = ε˙e − ct ε˙e , for the plastic case, the inequality D˙ > 0 can be written as D˙ = − σd − α − rˆn ˆ : ct ε˙e > 0 (8)
Hence, the ddissipation is maximum, and the inequality is always guaranteed when ct ε˙e and σ − α − rˆn ˆ have the opposite direction, i.e. in this case n ˆ may also be defined as σd − α (9) n ˆ= kσ d − αk
so the flow rule − ct ε˙e ≡ ε˙p = γ˙ n ˆ and the hardening rule ξ˙ = − ct ε˙e = γ˙ n ˆ are defined in terms of this tensor n ˆ . The multiplier γ˙ ≥ 0 is now interpreted as the norm of the plastic strain rate (the equivalent plastic strain rate, up to a constant). This variable is non-negative because D˙ ≥ 0. Then we can write
D˙ = σd − α − rˆn ˆ γ˙ ≡ 0 r γ˙ ≡ σY p˙ ≥ 0 (10)
where we have defined 0 r := σd − α − rˆ ≥ 0, and k·k denotes the 2-norm. The p term σY p˙p is the uniaxial equivalence, where σY = 3/2 0 r is the uniaxial yield stress and p˙ = 2/3γ˙ is the uniaxial plastic strain rate, obtained from r Z t r 2 2 p= kε˙p k dt = γ (11) 3 0 3 Then, using σ = α + (0 r + rˆ) n ˆ and ε˙p = γ˙ n ˆ , the stress power P may be written as the sum of the following terms σ : ε≡σ ˙ : tr ε˙e = σ : ε˙e − ct ε˙e = σ : ε˙e + σ : ε˙p (12) ˙ Ψ D˙ }| { z z}|{ (13) ≡ P = σ : ε˙e + α : ξ˙ + rˆγ˙ + 0 r γ˙ | {z } | {z } |{z} ˙ W H˙ kin H˙ iso | {z } ˙ H
6
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
so the integral of these terms are the energies shown in Fig. 1: Z Z Z Z ˙ W = σ : ε˙e dt, Hkin = α : ξdt, Hiso = rˆγdt, ˙ D=
0
r γdt ˙
(14)
R and the area below the stress-strain curve, Pdt, is the external (“stress”) work. The energy W is immediately recovered upon unloading. The energy Hkin may be recovered when plastic flow happens in the reverse direction; note that the product α : ξ˙ may have positive or negative sign. The energy Hiso is usually not recovered during cyclic loading (isotropic hardening) because rˆγ˙ is positive for hardening (as seen below); however, it may be recovered through thermal treatments. Finally, the dissipation energy D is never recovered; it is dissipated at each instant, typically, into heat. Sometimes Hiso is considered as part of the dissipation (phenomenologically, if thermal effects are neglected, there is no difference), so in such case we could define D˙ = (0 r + rˆ) γ. ˙ Equation (10) may be written in the following convenient format for identification of the yield surface
D˙ = σ d − α − rˆ − 0 r γ˙ + 0 r γ˙ ≥ 0 (15) which implies that (case of plastic loading)
f := σd − α − r = 0 if γ˙ > 0
(16)
where we defined the current radius r := 0 r + rˆ, but we can have (case of elastic loading or neutral loading) f ≤ 0 if γ˙ = 0 (17) During plastic deformations, the consistency parameter γ˙ > 0 can be obtained from the consistency condition which must be f = 0 during plastic flow (γ˙ > 0), so f˙ = 0 (or equivalently γ˙ f˙ = 0 in any case, elastic or plastic) tr d (18) +n ˆ : ct σ˙ d − α˙ − r˙ = 0 f˙ (ε, εp ) = n ˙} |ˆ : {z σ {z } | tr ˙ ct ˙ f ≡ f˙ f ≡ f˙ ε˙p =0 ε=0 ˙
where we define the trial deviatoric stress rate as tr
σ˙ d := 2µtr ε˙d ≡ 2µε˙d
7
(19)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
and the corrector one as ct
σ˙ d = 2µ ct ε˙e ≡ −2µε˙p = −2µγ˙ n ˆ
(20)
The latter is the rate of the stress tensor during plastic flow which produce the maximum dissipation. Note that this novel continuum framework presented in [59], [60] mimics the usual algorithmic implementation, even at large strains, but it is just an immediate application of the chain rule when considering the dependencies σ(ε, εp ), i.e. ∂σ ∂σ : ε˙ + : ε˙p = tr σ˙ + ct σ˙ (21) σ˙ (ε, εp ) = ∂ε ε˙p =0 ∂εp ε=0 ˙
This type of decomposition is also recovered at large strains and its backward-Euler discretization results immediately in the typical Closest Point Projection algorithms. The component of the rate of α perpendicular to n ˆ does not affect the yield condition and does not change the dissipation, so as it is well known, the direction of minimum accumulation of internal elastic energy, and maximum dissipation for α, ˙ is the one given by Prager’s rule; see e.g. [19]. Accordingly, we assume Prager’s translation rule, so α ˙ is proportional to n ˆ ∂ 2 H (ξ, γ) ˙ ∂ 2 H (ξ, γ) α˙ = :ξ+ γ˙ ∂ξ∂ξ ∂ξ∂γ ∂α ˙ ∂α ˙ n =: 2 H ¯ (ξ, γ) γ˙ n :ξ+ γ˙ =: λˆ ˆ = 3 ∂ξ ∂γ
(22)
¯ γ, ˙ where ∂α/∂γ = ∂r/∂ξ. The backstress rate modulus is given as λ˙ = kαk ˙ =: 32 H or equivalently, multiplying Eq. (22) by n ˆ γ˙ =
hα˙ : n ˆi (ξ, γ)
2 ¯ H 3
(23)
and we defined 2 ¯ H 3
(ξ, γ) := n ˆ:
∂ 2 H (ξ, γ) ∂ 2 H (ξ, γ) :n ˆ+ :n ˆ ∂ξ∂ξ ∂γ∂ξ
(24)
as the effective kinematic hardening. The Macaulay bracket in hα˙ : n ˆ i implies that γ˙ ≥ 0. Then, from the consistency condition, Eq. (18), we restrict the formulation to a dependence r (γ) —a dependence of r on ξ would include a path-dependent 8
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
isotropic hardening, found in some materials [36], but which for simplicity we do not consider in this work— so r˙ :=
∂ 2 H (ξ, γ) ∂ 2 H (ξ, γ) dr γ˙ ≡ r ′ (γ) γ˙ := : γ˙ n ˆ+ γ˙ dγ ∂γ∂ξ ∂γ∂γ
(25)
where, note, both types of hardening are implicitly coupled through the term ∂ 2 H/∂γ∂ξ as it will be seen below. Then, we can factor out the consistency parameter as usual γ˙ =
2µ d n ˆ : ε ˙ ¯ (ξ, γ) + r ′ (γ) 2µ + 23 H
(26)
¯ (ξ, γ) and r ′ (γ) are the functions to be obtained from experiments. Then, where H it is immediate to obtain the continnum tangent Cep
(2µ)2 = κI ⊗ I + 2µP − ˆ ⊗n ˆ ¯ (ξ, γ) + r ′ (γ) n 2µ + 32 H
(27)
Remarkably, these are the classical J2 plasticity expresions given, for example, in Sec. 2.3.2 of [19] or in Sec. 4.3 of [63], where the previous equations are obtained ¯ Then, the single important from an equivalent procedure but assuming a constant H. difference in Eqs. (26) and (27) respect to the classical ones is that we do not assume ¯ (ξ, γ) is constant, but instead that it changes according to the flow direction that H and the history of deformation by way of the implicit dependences on ξ and γ, as ¯ results constant, both the shown below. However, in the particular case that H theory and computational algorithms in [19] (see also [64]) are naturally recovered regardless of the number or surfaces employed, and without user intervention; i.e. the solution is reached in just one iteration if r ′ is also constant. Therefore, our model constitutes a natural extension of the classical computational algorithms in finite element software to account for history-dependent non-linear kinematic hardening and anisotropic hardening effects. 2.2. Description of anisotropic hardening The anisotropic hardening is described in our model with the aid of different, initially concentric, hardening surfaces (i.e. not macroscopic yield surfaces). We remark that we follow the idea of Mr´oz and Iwan in using concentric surfaces for describing the hardening, but our theory is different from Mr´oz’s theory because our hardening surfaces are not yield surfaces at the macroscopic level, but just a tool ¯ which take to compute an effective macroscopic anisotropic hardening modulus H into account microstructural changes due to the deformation history in the material. 9
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
a.
b.
Yield surface
Hardening surfaces
Mróz-like view of the theory
Hereditary view of the theory
c. Trial stress
Yield surface (final stress)
Backward-Euler / radial-return stress integration
Figure 2: Geometric description of the model and the integration algorithm in the deviatoric plane. (a) Multisurface (Mr´ oz-like) view — Yield and hardening surfaces fi , flow direction n ˆ , backstresses αi , contact points σ ˆ i and hardening directions m ˆ i . (b) Hereditary (endochronic-like) view — Yield surfaces f¯i for internal variables αi . (c) Backward-Euler, radial-return integration algorithm; for simplicity only kinematic hardening is considered in the figure.
10
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Furthermore, we keep Prager’s hardening rule instead of changing to Mr´oz’s or other alternative rule. These hardening surfaces may be written as, see Figure 2a fi = kˆ σ i − αi k − ri = 0, i = 1, ..., N
(28)
where N is the number of surfaces and σ ˆ i is the contact point in the deviatoric stress space of surface i with surface (i − 1), and we define the first, innermost surface as the (macroscopic) yield surface, i.e. f1 ≡ f , σ ˆ 1 ≡ σd , α1 ≡ α, r1 ≡ r
(29)
As another phenomenological difference with Mr´oz’s proposal, the contact point between surfaces is not, in general, the same point, i.e. σ ˆ i 6= σ ˆ j for i 6= j. Equations (28) are convenient for visualization. However, they are fully equivalent to the following equations, more convenient for mathematical derivations, see Figure 2b
f¯1 ≡ f1 = σd − α1 − r1 = 0 (30) f¯i = kαi − αi−1 k − ρi = 0, i = 2, ..., N
where we have defined ρi = ri − ri−1 > 0 (and ρ1 = r1 ). Hence, the hardening surfaces can be alternatively interpreted as yield surfaces for internal microstructural stresses αi . These internal stresses are a result of internally stored (blocked) energy, which is returned upon stress reversal. Interestingly, the view in Figure 2b shows an interpretation of the theory close to that of hereditary models, as for example the endochronic theory [23], [67]. 2.3. Consistency condition in hardening surfaces For those surfaces such that fi < 0 (or equivalently f¯i < 0), the hardening surfaces remain fixed because the inner surfaces are fully inside, not touching surface i. However, if fi = 0 for any i > 1, this implies that surface (i − 1) is touching surface i. In such a case, we must guarantee that f˙i ≤ 0 (or equivalently df¯i /dt ≤ 0) in order to avoid overlapping surfaces which would yield an inconsistent definition of the hardening function in the stress domain. From Eq. (302 ) αi − αi−1 df¯i = : (α ˙i −α ˙ i−1 ) − ρ˙ i = 0 dt kαi − αi−1 k
11
(31)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Define —this tensor can be interpreted as the flow direction for the internal stresses αi of microstructural origin, see Fig. 2b m ˆ i :=
αi−1 − αi mi := kmi k kαi−1 − αi k
and m ˆ1 ≡n ˆ
(32)
so
df¯i =m ˆ i : (α ˙i −α ˙ i−1 ) + ρ˙ i = 0 (33) dt The first surface is the actual yield surface, so as seen above, it translates according to Prager’s rule, α ˙ 1 = kα˙ 1 k n ˆ . Then the consistency on the second surface is −
m ˆ 2 : α˙ 2 − kα ˙ 1 k (ˆ n:m ˆ 2 ) + ρ˙ 2 = 0
(34)
Note that as in the case of the yield surface, the component of α˙ 2 perpendicular to m ˆ 2 is irrelevant here, but the minimum translation for a given kα˙ 1 k takes place for translation in direction m ˆ 2 , i.e. α˙ 2 = kα ˙ 2k m ˆ 2 —this can be seen as Prager’s rule ¯ applied to fi , interpreted as yield conditions at different microstructural levels kα ˙ 2 k = kα˙ 1 k (ˆ n:m ˆ 2 ) − ρ˙ 2
(35)
kα˙ i k = kα ˙ i−1 k (m ˆ i−1 : m ˆ i ) − ρ˙ i ≥ 0
(36)
In general, It is important to note that if the position of the inner surface is known and contact is taking place, from f¯i = 0, Eq. (302 ), and using α˙ i = kα ˙ ik m ˆ i , we can obtain the position of the outer surface in an explicit manner, see Figure 2 αi = αi−1 − ρi m ˆi
(37)
which is a condition easier to use in the algorithmic (discrete) implementation. Note that defining ρi = ρi m ˆ i , we can also have the interpretation of the backstress as the addition of several overstresses (which could in turn be interpreted as an internal equilibrium equation) α = αa+1 +
a X
ρi and σ = αa+1 +
a X i=1
i=2
where a is the outermost active surface.
12
ρi
(38)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
2.4. Consistency parameter from multiple hardening surfaces Equation (23) is directly applicable for one surface. In the presence of N hardening surfaces, we assume the existence of some strain-like internal variables ξi , i = 1, ..., N (of microstructural origin) describing an anisotropic hardening field, such that we have the dependency ξ (ξi ), and in particular we assume that they follow the compatibility equation ξ˙ =
N X
ξ˙i
(39)
i=1
In this work we do not include ratcheting effects. However, these effects may be included in the theory by the addition of ratchetting internal strains in the previous equation. From the compatibility equation of internal strains, we have at the continuum observable level that the consistency parameter rate is γ˙ = ξ˙ : n ˆ=
a X
ξ˙i : n ˆ=
a X
γ˙ i
(40)
i=1
i=1
where we defined γ˙ i := ξ˙i : n ˆ and a ≤ N is the outermost active surface, i.e. a = arg max (fi = 0) i
(41)
If we associate a hardening modulus Hi (γ) to each one of the hardening surfaces such that α ˙i (42) α˙ i = 32 Hi (γ) ξ˙i so ξ˙i = 2 H (γ) 3 i
then by Eq. (40)
a
a
X X hα˙ i : n hα ˙ :n ˆi ˆi γ˙ = 2 ¯ = γ˙ i = 2 H (ξ, γ) Hi (γ) 3
i=1
(43)
i=1 3
In this work we restrict the dependence of the moduli Hi (γ) to deppend only on the ¯ (ξ, γ) is implicitly defined by cummulative parameter γ. The other dependency in H
13
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
the history of the backstresses αi and their flows ξ˙i . Using Eq. (36) recursively γ˙ i =
hα ˙i : n ˆi kα ˙ ik = 2 hm ˆi :n ˆi 2 H H 3 i 3 i i
Qi hm ˆi :n ˆ i X Qi−1 hm ˆi :n ˆi kαk ˙ ˆj :m ˆ j−1 i − 2 ρ˙ j k=j (m ˆk :m ˆ k+1 ) = 2 j=2 hm H H 3 i 3 i j=2
(44)
The first term in Eq. (44) is the most important one, whereas the other terms represent the coupling between isotropic and kinematic hardening, and vanish even with isotropic hardening in the case that r˙j = r˙j−1 . Since i X j=2
ρ˙ j = r˙i − r˙1
(45)
for the case of uniaxial loading, we have m ˆi =n ˆ for all i, kαk ˙ r˙i − r˙1 − 2 2 H H 3 i 3 i
(46)
kαk ˙ (r˙i − r˙1 ) / kαk ˙ kαk ˙ − kαk ˙ = 2 [1 − ζia (p)] 2 2 H H H 3 i 3 i 3 i
(47)
γ˙ i =
γ˙ i =
where we have defined ζia (p) :=
r˙i − r˙1 σˆ ′ (p) − σ ˆ1′ (p) = i ¯a kαk ˙ H
(48)
in which σ ˆi′ (p) is the derivative of the uniaxial-equivalent size of surface i as a function ¯ a (with superindex a) is the uniaxial equivalent of the uniaxial plastic strain p, and H hardening modulus when the outermost active surface is surface a. Then kα ˙ ik z }| { 2 kαk ˙ = Hi γ˙ i + (r˙i − r˙1 ) 3
(49)
and, from Eq. (43), in the particular case of proportional loading a
a
X 1 X 3 1 3 γ˙ = kαk ˙ ¯a = ˙ (1 − ζia) γ˙ i = kαk 2 2 H H i i=1 i=1 14
(50)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
This equation results in
a
X 1 − ζa 1 i = a ¯ Hi H i=1
(51)
which is the expression we may use to determine the hardening moduli associated to the surfaces from a uniaxial stress-strain curve. We note that during a monotonic loading it is impossible to distinguish which part of the hardening corresponds to isotropic and which part to kinematic hardening. Only upon stress reversal, identifying the position of the backstress (at half the elastic region), both components can be identified. For a given surface a−1
X 1 − ζa 1 1 1 i = for a = 1, ..., N − aH a a ¯ Ha 1 − ζ (1 − ζa ) H i a i=1
(52)
This expression simplifies in the usual case where isotropic hardening is small respect to kinematic hardening or all the surfaces harden equally, i.e. ζia ≪ 1, to a−1
X 1 1 1 = ¯a − for a = 1, ..., N Ha H H i i=1
(53)
In this case the modulus is given by the uniaxial equivalence ¯a σˆa+1 − σˆa EH (a) ≃ Eep = a ¯ εˆa+1 − εˆa E +H
(54)
(a)
in which E is the elastic Young modulus, Eep is the uniaxial elastoplastic modulus while surface a is the outermost active surface, and σ ˆa and εˆa are the stress-strain data points in which the uniaxial stress-strain curve has been discretized, see Figures 3 and 4. Then (a) E − Eep a ¯ H = (55) (a) EEep (N )
For the case a = N, the value of Eep , and hence HN , is given by the desired limiting hardening modulus. Remarkably, if we insert a hardening surface in a straight section of the stress(i+1) (i) ¯ i+1 = H ¯ i and 1/Hi+1 = 0. Then note that strain curve, we have Eep = Eep , so H γ˙ i+1 = 0 and the surface becomes irrelevant in terms of plastic work. Furthermore, if the uniaxial stress-strain curve happens to be bi-linear, naturally all 1/Hi = 0 except 1/H1 , and γ˙ i = 0 except for γ˙ 1 = γ. ˙ Therefore, the classical linearly-hardened 15
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 3: Example of a discretization of a stress-strain curve and the generation of the hardening √ surfaces; depicted in the 3τ − σ plane.
J2 -plasticity is recovered naturally without explicitly accounting for it, just as a consequence of the prescribed stress-strain curve. This is extremely important in order to obtain consistent predictions for different, arbitrary discretizations, a property not enjoyed, for example, by Mr´oz’s model [68]. 2.5. Cyclic isotropic hardening/softening There are many possible implementations for isotropic hardening/softening in the multisurface model. Since we are dealing with multiple surfaces to describe hardening, these surfaces may be used even to change the shape of the stress-strain curve during cyclic loading, for example from the monotonic backbone curve to a very different skeleton one. However, the purpose of the present paper is to show the capabilities of the model to describe multiaxial loading with cyclic hardening/softening. Therefore, since in the examples we will not account for a relevant change of shape in the skeleton curve, to keep the presentation simple, we consider an isotropic hardening function such that all surfaces remain proportional in size ri (p) = R (p) 0 ri
16
(56)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 4: Procedure for modification of the skeleton curve due to cyclic hardening.
where 0 ri = ri (0), R (p) is the isotropic hardening function and we made, as usual, an abuse of notation in ri (p) = ri (γ) to avoid proliferation of symbols. Then r dri 2 ′ 0 0 ′ r˙i (p) = p˙ = ri R (p) p˙ = R (p) ri γ˙ (57) dp 3 and obviously ρ˙ i (p) = r˙i (p) − r˙i−1 (p) =
0
ρi R′ (p) p˙
(58)
An scheme of the surface growth may be seen in Figure 4. The procedure that we use to determine the function R (p) is the follow-up of the peak stresses during subsequent cycles, which is a procedure useful in fatigue analysis. With the data from the change of peak stresses and the accumulated plastic strain during a cycle, we obtain the relation between R and p. With the knowledge of R (p), and hence of the change in the stresses σ ˆa , the derivative of the radius of any hardening surface dri /dp can be obtained. Hence, the relevant key function is R (p). A popular way of expressing the function R (p) is in relative terms, i.e. R/Rs ([27], [57]), where Rs is the saturated value. During cyclic loading, we can approximate R (p) measuring the peak stresses in succesive strain-controlled cycles. Assume that
17
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
when the strain range of a loop is known, the accumulative plastic strain of that first loop is approximately p0 . Then, during cyclic loading, the accumulated plastic strain is approximately pn = 0.25p0 + np0 , where n is the number of loops. The change of R (p) at different loop counting n, gives pairs {pn , Rn } , which can be immediately used to establish the relation R (p) in any desired format (piecewise linear, spline or adjusting a model function). Once the relation R (p), or more generally, the relation ri (p) is known, the expression for the derivative of the kinematic hardening moduli Hi can be obtained, see Figure 4. The inverse of the uniaxial elastoplastic tangent when surface a is active is, for a given p εa+1 (p) − εa (p) 1 1 1 = = + (59) (a) ¯ a (p) E H σ ˆa+1 (p) − σ ˆa (p) Eep (p) Since the elastic behavior is constant, the variations during a cycle are related by— see Figure 4 1 ∆ [εa+1 (p) − εa (p)] = (60) ∆ [ˆ σa+1 (p) − σ ˆa (p)] E Then, using the chain rule, taking the derivative of Eq. (59) respect to p—note that as it should be expected, this derivative vanishes if all surfaces harden equally in an isotropic manner, so dˆ σa+1 /dp = dˆ σa /dp— ¯ a (p) d 1/H dˆ σa+1 (p) dˆ 1 εa+1 (p) − εa (p) σa (p) (61) = − − dp E [ˆ dp dp σa+1 (p) − σ ˆa (p)]2
where dˆ σa = dp
r
2 dra 3 dp
(62)
Thus, from Eq. (51) we can compute in a recursive manner a−1 ¯ a (p) X d (1/Hi (p)) d 1/H d (1/Ha (p)) = − dp dp dp i=1
(63)
d (1/Hi (p)) dHi = −Hi2 (p) dp dp
(64)
and, obviously
In summary, from R (p) we can immediately compute ri (p), dri /dp and dHi /dp for all surfaces (the hardening of the last surface may be set constant). If a change of the shape of the stress-strain curve from monotonic to cyclic behavior is desired, 18
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
different functions for ri (p) can be selected to model such effect. 3. Fully implicit algorithmic implementation In this section we introduce a fully implicit integration algorithm based on the two-parameter Closest Point Projection (CPP) Algorithm [69]. 3.1. Incremental formulation The stress integration algorithm computes the stress t+∆t σ at “time” (or step) t + ∆t with the knowledge of: the elastic strains t εe at time t, the increment of total deformations ∆ε during the step, the history variables t αi , and either t p or t γ. The time-integration scheme, based on the classical radial return algorithm is depicted in Figure 2c. Using the trial stress rate Eq. (19), and defining the algorithmic trial stress tr σ = t σ + ∆tr σ, tr σ ≡ σ t εe + ∆ε = t σ + κ∆εv I + 2µ∆εd (65) Since during the absence of plastic flow the backstresses do not change (i.e. we can define the algorithmic backstresses tr
tr
α˙ i = 0),
αi := t αi , i = 1, ..., N
(66)
The corrector stress increment, using a Backward-Euler (radial-return) integration algorithm, is ∆ct σ = 2µ∆ctεe = −2µ∆γ t+∆t n ˆ (67) where
t+∆t
n ˆ is the normal to the yield surface at the end of the step. Then t+∆t
σ = tr σ − 2µ∆γ
t+∆t
n ˆ
(68)
Using the same integration rule, from Eq. (22) for surface 1 t+∆t
α = tα +
t+∆t
˜ 1 ∆γ1 H
t+∆t
n ˆ
(69)
where for convenience, we have absorbed the uniaxial equivalence factor 2/3 in t+∆t
˜ i := 2 H 3
19
t+∆t
Hi
(70)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Then, note that as in classical plasticity, we can obtain the flow direction directly from the deviatoric trial state—recall that the volumetric part is elastic t+∆t
tr
n ˆ=
t+∆t σd − t α n =: tr d t t+∆t k σ − αk k nk
(71)
However, to obtain t+∆t α in Eq. (69), which accounts for the correction to the dissipation (the energy not dissipated but microstructurally stored), we need to compute the consistency parameters which guarantee material compatibility. In order to do so, we need to develop an integration algorithm. We will be able to establish below a very efficient iterative algorithm based just on two scalars, namely ∆γ1 and ∆γ (i.e. a two-parameter CPP iterative algorithm), regardless of the number of surfaces being employed. The two functions to be minimized are the yield function f and the compatibility requirement g (∆γ, ∆γ1 ) = ∆γ − ∆γ a → 0
(72)
where we defined —note the different meaning of sub- and superindices ∆γ a (∆γ1 ) :=
a X
∆γi (∆γ1 ) = ∆γa (∆γ1 ) + ∆γ a−1 (∆γ1 )
(73)
i=1
The last recursive identity is given for future reference. In order to minimize the function, we need the derivatives respect to these two basic iterative variables, namely ∆γ and ∆γ1 . From Equation (37), using also Eq. (32), after some straightforward algebra, we can obtain dt+∆t αi−1 d∆γ a dt+∆t αi = t+∆t D1i : − t+∆t D2i (74) d∆γ1 d∆γ1 d∆γ1 where the algorithmic projectors and algorithmic radii vector changes are t+∆t t+∆t ρi ρi t+∆t t+∆t D1i = 1 − t+∆t P + t+∆t m ˆ i ⊗ t+∆t m ˆi k mi k k mi k r t+∆t 2 dt+∆t ρi t+∆t dt+∆t ρi R t+∆t 0 d D2i = m ˆ with = ρ i i t+∆t t+∆t t+∆t 3d p d p d p
20
(75) (76)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
in which we defined—note that t+∆t
mi :=
t+∆t
t+∆t
m ˆ i is also recursively defined
αi−1 − t αi and
t+∆t
t+∆t
m ˆ i :=
k
mi t+∆t m k i
(77)
and P is the fourth order deviatoric projector tensor; i.e. P = I− 31 I ⊗ I, where I is the fully symmetric 4th-order identity tensor. Then, taking the derivative of Eq. (43) respect to ∆γ1 , we can compute d∆γ a /d∆γ1 in a recursive manner using Eq. (73)2 . Consider ! b−1 X d∆γb d∆γi d∆γ b−1 d∆γb d∆γ b + = = + (78) d∆γ1 d∆γ1 d∆γ1 d∆γ1 d∆γ1 i=1 After some straightforward algebra we obtain t+∆t
dt+∆t αb−1 ∂∆γ b−1 + d∆γ b d∆γ1 ∂∆γ1 = q t+∆t t+∆t ˜b d∆γ1 n ˆ : t+∆t m ˆ b d t+∆t H Zb t+∆t 1+ n ˆ : t+∆t D2b + t+∆t Zb ∆γb 23 t+∆t H t+∆t H ˜b ˜b d t+∆t p (79) which is applied repeatedly for b = 1, ..., a. In this equation, to account in Eq. (43) for γ˙ i ≥ 0, we defined for convenience the Heaviside function t+∆t t+∆t 1 + sign m ˆ : n ˆ b t+∆t Zb = ∈ {0, 1} (80) 2 Zb t+∆t H ˜b
t+∆t
n ˆ:
t+∆t
D1b :
(b) (b−1) ¯b = H ¯ b−1 and H ˜ b → ∞; then ∆γ b = ∆γ b−1 , Remarkably, if Eep = Eep we have H d∆γ b /d∆γ1 = d∆γ b−1 /d∆γ1 and surface b has no effect in the formulation or the iterations.
3.2. Local iterative algorithm As mentioned, the implicit algorithm is based on the search for the values of ∆γ1 and ∆γ such that the following conditions are met t+∆t ˜ 1 ∆γ1 − t+∆t r1 → 0 f (∆γ, ∆γ1 ) := ktr σ − t αk − 2µ∆γ − t+∆t H (81) t+∆t g (∆γ, ∆γ1 ) := ∆γ a (∆γ1 ) − ∆γ → 0
21
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
The solution for [∆γ, ∆γ1 ]T is obtained by a Newton-Raphson algorithm until convergence is achieved in both f and g,
∆γ ∆γ1
(j+1)
=
∆γ ∆γ1
(j)
−
∂∆γ/∂f ∂∆γ/∂g ∂∆γ1 /∂f ∂∆γ1 /∂g
(j)
f g
(j)
(82)
where [·](j) implies values at iteration j. The derivatives in the tangent matrix in Eq. (82) are obtained as usual from its inverse matrix, whose components are ∂f ˜1 = − t+∆t H ∂∆γ 1 q ˜ 1 q d t+∆t r1 d t+∆t H ∂f = −2µ − 23 ∆γ1 t+∆t − 23 t+∆t ∂∆γ d p d p (83) a d∆γ ∂g = ∂∆γ1 d∆γ1 ∂g = −1 ∂∆γ which right-hand-side terms have already been computed above, see Eqs. (79) and (57). Remarkably, if linear hardening is considered, in which case only one surface ˜ i /dp = 0, dri /dp is constant and ∆γ a = ∆γ1 . Then, the function g is is relevant, dH identically satisfied and the algorithm gives the exact solution in just one iteration because all coefficients in the tangent are constants. As expected, this result is in line with the well-known classical radial return algorithm of Krieg and Key [58], so both the continuum and algorithmic classical plasticity are recovered. In the case when only isotropic hardening is nonlinear, g is also identically satisfied and the classical scalar-based algorithm is also recovered, see Ref. [19]. For the case when multiple surfaces are present, a guess (·)(0) is needed to start (0) the iterations. We always take ∆γ1 = 0. A simplest guess for the active surface is to assume a(0) = 1. If after the computation of the solution we obtain f¯2 > 0, then the index is increased a ← a+1, and the solution is recomputed; and so on. However, if a large number of surfaces is employed, a better guess to save computational time is to set a for the first iterations as the outermost surface such that f¯a = 0 and
t+∆t
n ˆ : tm ˆa >0
(84)
If after reaching a solution f¯a+1 > 0, then the active surface index is increased. On the contrary, if a solution with ∆γa > 0 is not obtained, the active surface index is 22
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
decreased. A layout of the computational algorithm is shown in Table 1 3.3. Consistent tangent for the global iterative algorithm To preserve the efficient asymptotic quadratic convergence rate of implicit Newton algorithms in global equilibrium iterations in finite element simulations, the algorithmic (consistent) tangent moduli must be used. These moduli relate the increment of stresses to the increment of strains according to the algorithmic implementation, i.e. t+∆t
dt+∆t σd dεd dt+∆t σ C = t+∆t = κI ⊗ I + t+∆t d : d ε d ε dε
(85)
where κI ⊗ I is the volumetric part (purely elastic and constant), dεd /dε ≡ P is the also constant deviatoric projector tensor and dt+∆t σd /dt+∆t εd is the relevant nontrivial algorithmic part. This tensor may be obtained taking the derivative of—see Figure 2c t+∆t d ˜ 1 ∆γ1 + t+∆t r1 ) t+∆t n σ = t α + (t+∆t H ˆ (86) as dt+∆t σ d dt+∆t n ˆ d∆γ1 t+∆t ˜ t+∆t ˜ 1 t+∆t n = ( H ∆γ + r ) + t+∆t H ˆ ⊗ t+∆t d 1 1 1 t+∆t d t+∆t d d ε d ε d ε t+∆t ˜ t+∆t d H1 d r1 + ∆γ1 t+∆t n ˆ ⊗ t+∆t d + t+∆t n ˆ ⊗ t+∆t d d ε d ε p Taking into account that ∆εp = ∆γ t+∆t n ˆ = 3/2∆p t+∆t n ˆ , we obtain dt+∆t n ˆ d∆γ dt+∆t εp = ∆γ + t+∆t d t+∆t d d ε d ε d∆γ1
and
dt+∆t p = dt+∆t εp
r
2 3
t+∆t
n ˆ , so
t+∆t
n ˆ⊗
d∆γ1 dt+∆t εd
dt+∆t p dt+∆t p dt+∆t εp : = dt+∆t εd dt+∆t εp dt+∆t εd
(87)
(88)
(89)
Using the chain rule ˜ 1 dt+∆t p ˜1 dt+∆t H dt+∆t r1 dt+∆t r1 dt+∆t p dt+∆t H = and = dt+∆t εd dt+∆t p dt+∆t εd dt+∆t εd dt+∆t p dt+∆t εd
23
(90)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Local iterative algorithm 1. Compute trial stress tr σ, Eq. (65) and both t+∆t n and t+∆t n ˆ , Eq. (71)
2. If tr f := t+∆t n − t r1 ≤ 0, step is elastic; accept trial solution and exit. Else, step is plastic: continue
3. Set ∆γ = ∆γi = 0 and guess a; either a = 1 or by Eq. (84) 4. For each iteration (j) Compute Set
t+∆t
t+∆t
˜i H
t+∆t
α1 = t α1 +
p and
t+∆t
t+∆t
˜ 1 ∆γ1 H
ri
t+∆t
t+∆t
p
with
t+∆t
p=
p
2/3 (t γ + ∆γ)
n ˆ and initialize ∆γ a = ∆γ1
For each hardening surface i = 2, ..., a Compute
t+∆t
mi and
t+∆t
m ˆ i via Eq. (77)
By Eq. (37) t+∆t αi = t+∆t αi−1 − t+∆t ρi t+∆t m ˆi ˜ i t+∆t αi − t αi : t+∆t n ∆γi = 1/t+∆t H ˆ
∆γ a ←− ∆γ a + ∆γi
Compute t+∆t f = t+∆t n − 2µ∆γ −
t+∆t
˜ 1 ∆γ1 − H
t+∆t
r1
Compute t+∆t g = ∆γ a − ∆γ If t+∆t f < tolf and t+∆t g < tolg accept the solution, go to step 5. Perform loop i = 2, ..., a to compute d∆γ a /d∆γ1 via Eq. (79) Compute the local tangent matrix, Eqs. (83) and (79) Update solution solving 2 × 2 linear system of Equations, e.g. Eq. (82) (j) 5. Solution: t+∆t εe = tr εe − ∆γ t+∆t n ˆ , t+∆t σ = σ t+∆t εe and t+∆t αi = t+∆t αi 6. If f¯a+1 > 0, a ←− a + 1 and go to step 4.
If solution needs ∆γ a < 0, a ←− a − 1 and go to step 4
7. Compute global tangent; see Section 3.3 Table 1: Local stress integration algorithm
24
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Taking into account that by Eq. (71)
d t+∆t n dt+∆t εd
= 2µ
2µ dt+∆t n ˆ I− = tr t+∆t d d ε k σ − t αk
t+∆t
n ˆ
t+∆t
(91)
n ˆ⊗
t+∆t
n ˆ
and that by Eq. (81) 0≡
dt+∆t f d∆γ = 2µ t+∆t n ˆ − 2µ t+∆t d − t+∆t d d ε d ε t+∆t ˜ t+∆t d H1 d r1 − t+∆t d ∆γ1 − t+∆t d d ε d ε
t+∆t
(92)
˜ 1 d∆γ1 H dt+∆t εd (93)
after some lengthy but straightforward algebra, we can factor-out the key derivative d∆γ1 = dt+∆t εd
2µ d∆γ + (2µ + t+∆t θ) d∆γ1
where t+∆t
θ=
r
2 3
t+∆t t+∆t H ˜
n ˆ
(94)
1
˜ 1 dt+∆t r1 dt+∆t H ∆γ1 t+∆t + t+∆t d p d p
!
(95)
Because at local convergence ∆γ a ≡ ∆γ, then d∆γ/d∆γ1 is given by Eq. (79). ˜ 1 = 2/3H ¯ and we Remarkably if only one surface is used we have d∆γ/d∆γ1 = 1, H naturally recover the usual expression for classical von Mises plasticity with mixed ¯ hardening if kinematic hardening is linear (dH/dp = 0). With the key Equation (94), all the quantities in the tangent Eq. (87) are known. The flow of the computational procedure for the tangent is Eq. (79) d∆γ Eq. (94) d∆γ1 −→ −→ t+∆t εp Eq. (89) Eq. (88) d d∆γ 1 dt+∆t εd −→ t+∆t d −→ t+∆t d n ˆ Eq. (92) d ε −→ t+∆t d d ε t+∆t ˜ 1 dt+∆t r1 Eq. (87) dt+∆t σd Eq. (85) p Eq. (90) dt+∆t H Eq. (89) d −→ t+∆t d −→ −→ t+∆t d −→ t+∆t C (96) , d ε dt+∆t εd dt+∆t εd d ε Obviously this diagram can also be employed to obtain the tangent in closed form. As it should be expected, it is straightforward to verify that if only one surface is 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
present, the classical tangent in Ref. [19] is recovered. 4. Large strains formulation As mentioned, the current formulation is a special case of classical von Mises plasticity with mixed hardening in which the effective kinematic hardening modulus is computed with the aid of internal variables. Therefore, as a remarkable advantage, the large strains implementation is immediate. To this end, we can use the algorithm of Eterovic and Bathe [65] based on logarithmic strains. This algorithm formulated in the full tensorial space is valid for kinematic hardening and uses a purely geometric pre- and post-processor, keeping the infinitesimal plastic stress integration unmodified. Using the multiplicative decomposition of the deformation gradient X = X e X p into elastic X e plastic X p parts, the elastic logarithmic strain in the intermediate configuration is obtained from the stretch tensor U e of the right polar decomposition X e = Re U e as E e = ln (U e )
(97)
The resulting incrementally additive framework using logarithmic strains mimics the infinitesimal theory, in which small elastic strains εe are interpreted as logarithmic elastic strains E e in the intermediate configuration, and Cauchy stresses σ as generalized Kirchhoff stresses T (or rotated Kirchhoff stresses in the elastic-isotropic case at hand). The hyperelastic stored energy is written in terms of these stress-strain measures as T = JU ′ (J) I + 2µE de (98) where JU ′ (J) and E de are respectively the volumetric and distorsional parts, the former obtained from the volumetric stored energy U ′ (J) as a function of the Jacobian determinant J. Algorithmically, the trial state is given by tr
Xe =
t+∆t t tX 0X e
(99)
and tr T d = 2µ tr E de with tr E de being the distorsional component of the trial elastic stretch tensor. Then, interpreting the backstress as a Kirchhoff-alike backstresses B, the small strains algorithm may be used without modification. Upon convergence, the final geometric update is performed from the converged elastic strains t+∆t t+∆t εe , as 0Ee ↔ tr t+∆t Re exp t+∆t0 E e (100) 0 Xe = 26
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
where we used the usual assumption that tr Re ≈ t+∆t0 Re [65], [63], [66], [59]. By systematic use of the chain rule, the tangent of the small strains algorithm may be converted to the typical tangent employed in finite element codes which relates second Piola-Kirchhoff stresses with Green-Lagrange strains in the reference configuration. Since the large strains pre- and post-processors are unchanged from the well-established classical plasticity and the experimental data used below are in the range of small strains, we refer for example to [65], [63] (Table 7.3.1), [66], [60], [59] for further details on the large strain numerical implementation using geometrical mappings. 5. Numerical examples In this section we include several numerical examples to demonstrate the predictive capabilities of the model and to demonstrate the performance of the fully implicit computational algorithm in actual finite element simulations. 5.1. Uniaxial example with cyclic isotropic hardening The purpose of this example is to show the capabilities to model cyclic hardening/softening during uniaxial cycles. The material for the simulation is steel BLY160. We used the material data from [57] and we predict the uniaxial cyclic experiments performed in Ref. [57]. In [57], the isotropic hardening is described with the relation between the normalized isotropic hardening parameter R/Rs and the accumulated plastic strain p, where Rs is the stabilized value. This function R/Rs is obtained in [57] by curve fitting from cyclic experiments shown in Fig. 18a of Ref. [57]. Since the saturated value Rs is also needed for our simulations, we take that data from Fig. 4 of Ref. [57], which contains experimental data of saturated stress values. Hence, from both figures in Ref. [57] we can obtain the needed relation R (p). However, we note that the shape of the skeleton curve is a saturated one with strain range 3% so the stress values have been adjusted by a ratio that is obtained by comparing the estimated yield stress of the saturated curve with the estimated yield stress of the first half loop curve. The resultant curve is the basic backbone which we used to obtain the initial kinematic hardening moduli field. With this data from Ref. [57], we performed simulations to predict the experimental observations shown in Fig. 3 of the same paper, which we reproduce in Figure 5a for the reader’s convenience. These experiments on three different specimens show the cyclic hardening behavior of the material at different deformation levels. In Figure 5b we show the predictions from our model. It can be seen that the predictions correlate quite well with the experimental observations. 27
a.
400
strain %
b.
300 200
stress (MPa)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
100 0 −100 −200 −300 −400 −0.04
−0.03
−0.02
−0.01
0
0.01
0.02
0.03
0.04
strain
Figure 5: (a) Experimental stress-strain cycles for prescribed strains of ∆ε = 1%, 2% and 3% for steel BLY160; adapted from [57], with permision form Elsevier. (b) Predictions of our model for the same experiments.
28
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
5.2. Examples under multiaxial, nonproportional loading. 5.2.1. Self-consistency of the formulation and recovery of classical bi-linear plasticity In order to show that classical von Mises plasticity with linear kinematic hardening is recovered in multiaxial cases regardless of the number of surfaces employed, we perform multiaxial simulations using as uniaxial stress-strain curve the bilinear one shown in Fig. 6a. In this figure we also show the locations of the different surfaces employed in the discretization and the slight modification to better appreciate the three cases in the simulations. We have prescribed the deformation path shown in Figure 6b. The result of the simulations using the three cases are shown in Figure 6c. It is seen that there is no relevant difference in the results. In all cases, the solution ˜ i /H ˜ i−1 is high enough to overcome the required is obtained in just one iteration if H tolerances in the local system of equations, except when a change of active surface is encountered. In such latter case, an additional iteration is needed just to update the position of the surface. 5.2.2. Experiments of Lamba and Sidebottom [61] We now show the predictions for the well-known nonproportional experiments of Lamba and Sidebottom [61] on oxygen-free high-conductivity (OFHC) copper. The nonproportional loading sequence (0 → 1 → · · · → 8) is shown in Figure 7. In this path, the positive and negative strain peaks are equal, and the angles between the different segments take very different values in order to test the response under change of direction, including a variety of angles [61]. To model these experiments we did not consider an isotropic hardening contribution. According to the loading sequence of Lamba and Sidebottom, the first segment 0 − 1 of the path is uniaxial shear loading, so the experimental stress-strain curve of this segment was used for generating the kinematic-hardening parameters of our model, see Eqs. (53) and (54), and the Poisson ratio to ν = 0.3. Figure 8 shows the comparison between our predictions (a) and the experimental results (b), the latter redrawn from Lamba and Sidebottom [61]. From the comparison of predicted and experimental results we can conclude that the model is capable of capturing most of the details of the loading history. Recall that the only data employed in the model is the proportional stress-strain curve for segment 0 − 1, which is shown in Figure 8 (1.a). A quantitative comparison of the von Mises stress in the characteristic points of the stress path is given in Table 2, where experimental values have been estimated from the plots in Ref. [61]. It is seen that the maximum differences are about 10%.
29
Uniaxial stress σ (GPa)
1
a
0.8 0.6
2 surfaces
0.4
4 surfaces
0.2
6 surfaces
0 0
0.2
0.4
0.6
0.8
Uniaxial strain ε (%) 1
1.5
1.2
1.4
1.6
c
b 1
Shear stress τ (GPa)
Shear strain γ/2 (%)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
0.5 0 −0.5 −1
0.5
0
−0.5 2 surfaces
−1
4 surfaces 6 surfaces
−1.5 −1
−0.5
0
0.5
−1.5
1
−1
−0.5
0
0.5
1
Axial stress σ (GPa)
Axial strain ε (%)
Figure 6: (a) Uniaxial stress-strain curve with different discretizations. Additional data: Poisson ˜ i = 103 H ˜ i−1 . (b) Prescribed strain path. (c) Predicted stress paths ratio ν = 0.3, H
30
Engineering shear strain γ %
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
1.5 1
6
1.0
5
0.5 0
0 −0.5
8 4 2
−1.0 −1.5 −0.4
7 3 −0.2
0
0.2
Axial strain ε %
0.4
a
Figure 7: Prescribed strain path in the tests from Lamba and Sidebottom [61]
Point 1 2 3 4 5 6 7 8
VM stress [MPa] (exp., approx.) 180.3 177.1 177.4 201.0 182.9 203.4 183.8 162.3
VM stress [MPa] (model) 185.5 185.1 166.1 178.5 164.5 179.3 198.8 164.8
% difference 2.88% 4.51% 6.37% 11.2% 10.1% 11.8% 8.16% 1.54%
Table 2: Approximate differences between experimental and computed von Mises stress at the turning points labelled in Fig. 7.
31
Shear Stress τ (MPa)
150 100
150 1(a)
5
50
4
0 −50
Axial Stress σa (MPa)
50
0
0
2
−150
100
100
−50
3 7
−100
1(b)
1 6
8
−100
Engrg.Shear Strain −0.01
0
γ
γ
Engrg.Shear Strain
0.01
−0.01
5 2
2(a)
0
−150
0.01
2(b)
100 0
0,1
0
7,8 3 6
−100
4
−200
100
3(a)
ε
0
0.004 −0.004 −0.002
0.002
Axial Strain
a
0
ε
−200
100
5
50
4
0
0
−50
−50
3
2
−100
7
−150
Axial Stress
−200
a
0.002 0.004
3(b)
1
8 6
50 0
−100
Axial Strain
−0.004 −0.002
Shear stress τ (MPa)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
−100
σ
a
0
−100 σ (MPa)
Axial Stress
(MPa) −200
100
−100
a
0
−150
100
Figure 8: Left column (a): Simulations of the Lamba-Sidebottom experiments using our model. The prescribed strain path is given in Figure 7. Right column (b): Experimental data from Lamba and Sidebottom, redrawn from [61]. The input data for the model is the monotonic curve of section 0-1
32
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
5.2.3. Experiments of Hamala et al [35] For further demonstration of the ability of the model for the description of multiaxial loading in other materials, we predicted the behavior of different specimens under four different strain-controlled tension-torsion loading paths given in Hamala et al [35], see also therein references. Then, the results of our simulations are compared with the experimental results in [35]. The uniaxial cyclic curve of the material (2124-T851 aluminum) given in Figure 8 of the same paper is used to determine our model, again using Eqs. (53) and (54). The four different multiaxial strain-loading paths (rhombus, circular, square and two-block) are simulated with a strain-driven procedure. The maximum amplitude in the prescribed √strain cycles shown in the small details in Fig. 9 is 1% (either axial ε or shear γ/ 3) in all cases. The comparison of the results from simulations and experiments along with prescribed strain paths are shown in Fig. 9. Since we did not consider isotropic hardening in this example, the first full cycle corresponds to a stabilized cycle in our model, and it is compared with the stabilized experimental results given in [35]. As shown in Fig. 9, the model predicts well the observed experimental stresspaths in all cases, in both the amplitudes and the different experimental shapes of the axial-shear stress√paths observed in the experiments. Since the plots given in the figure are in σ − 3τ axes, the distance to the origin gives a measure of the von Mises stress. The maximum differences in von Mises stress are found in the square and rhombus paths. In the corners of the square path, the maximum differences are about 15% at point labelled 2 in Fig. 9a, and about 10% at corner labelled 1. For the rhombus path, the maximum differences are found at corners labelled 1 and 3, with differences of about 10% and 18%, respectively. The maximum differences in the other two paths are about 10%. 5.2.4. Experiments of Tanaka et al [39] A series of experiments with nonproportional loadings were carried out also by Tanaka et al [36], [39]. Material parameters obtained from the experiments of Tanaka et al with loading paths of certain shapes have been used in different constitutive models for cyclic plasticity to derive a better description of the multiaxial behaviors of materials under those loading paths, see for example [39], [40], [41]. The prescribed loading path is given in Fig. 10a and the experimental result of the initial loading loop from [36] is redrawn in Fig 10b for comparison purposes. We note that the experimentally prescribed path given in [36], idealized in Fig. 10a, is given in terms of plastic strains, from machine-prescribed total strains with an iterative estimation therein done for the elastic part. We have performed our simulations using total strains, so an estimation of the maximum elastic contribution has been 33
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 9: Comparison of the predictions from our model against experimental cyclic data published in [35] for several prescribed multiaxial tension-torsion strain paths. (a) Square path. (b) Two-block path. (c) Rhombic path. (d) Circular path.
34
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
performed from stress values and the reported mean Young and shear modulus of the specimens, given in [36] as E = 203 MPa and G = 78 MPa. From Fig. 10b, it can be seen that the experimentally measured stresses in the initial strain-proportional path are not proportional, an effect also observed in other experiments in [36]. This can be attributed to experimental errors or to lack of compliance of the material state to the initial isotropic assumption. We show in Figure 10c the stress-strain curve employed in the simulations taken from Fig. 2 of Ref. [36], and in Fig. 10d the result of the simulations for the first cycle. Whereas there are several commented uncertainties that may affect the accuracy of the simulations, it is remarkable that the model captures well the shape of the stress path, including the special corner√effects in the path. The von Mises stress is the distance to the origin in the σ − 3τ plot. In Fig. 10d it is observed that the maximum difference is about 18% at the corner labelled 3, where isotropic hardening is already relevant. To include cyclic isotropic hardening in the simulations, from the peak stresses in the cycles in Fig 2 of Ref. [36] we have adjusted a function R(p), shown in Fig 10e. This curve, along with curve 10c was used for the simulations in Fig. 10d, which can be compared to those in Fig.3d of Ref. [36]. When including isotropic hardening, the error is about 10% at the mentioned corner. We note that we have not included path-dependent (directional) isotropic hardening. However, as commented, such effect could be included in the model allowing for a dependency r(γ, ξ). 5.2.5. Ramp-up loading experiments of Madrigal et al [71] In [71], a 90◦ tension-torsion out-of-phase experiment with continuously increasing strain amplitude is performed in AISI 4140 quenched and tempered steel, which is a material frequently used in the aerospace and automotive industries in elements that are subject to fatigue [71]. These experiments, along those from Lamba and Sidebottom, were used to test a previously proposed model [71], [72], [73], [74], similar to bounding surface models, but including a metric tensor which allows for the inclusion of different yield criteria, e.g. Tresca’s criterion [71]. This metric tensor also allows them to loosen the correlation between proportional axial and shear behaviors through additional material parameters in order to obtain better predictions for multiaxial loading. The 90o out-of-phase loading experiments with continuously increasing amplitude are challenging because they produce a large amount of additional multiaxial hardening and predictions are very sensitive to the uniaxial stress-strain curve and the size of the yield surface [71]. The experimental prescribed strain path of [71] is shown in Figure 11a. The experimental measured stress path obtained by Madrigal et al [71] is shown in Figure 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 10: Multiaxial cyclic experiments from [36]. (a) Prescribed strain path, 0 → 1 → 2 → 3 → 4. (b) Measured axial-shear stress path for the first cycle, redrawn from [36]. (c) Uniaxial monotonic stress-strain curve, obtained from Fig. 2 of [36] using a Young modulus of 203GP a and a shear modulus of 78GP a, which are average values given in Eq. (4) of [36]. (d) Predicted stress path for the first cycle. (e) Cyclic hardening curve R(p) − 1, obtained fitting approximate values of peak stress for subsequent cycles in Fig. 2 of [36]. (f) Predictions for several cycles including isotropic hardening (compare to Fig. 3d of [36]
36
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
11b. In order to perform the simulations, we have discretized the monotonic uniaxial stress-strain curve, given in Figure 5 of Ref. [71] and reproduced in Fig. 11c. We note that this material has a slightly different behavior in axial than in shear, which was the reason because of which Madrigal et al [71] developed a model including tensorial metrics. However, we performed the simulations using a von Mises criterion and taking the monotonic stress-strain curve for the model from axial loading; e.g. from Fig. 5 of Ref. [71]. In Figure 11d we show the results from our simulations for the ramp-up experiments with increasing amplitude. It is seen that the model reproduces the overall behavior, and quite accurately the axial stress values, having a difference precisely in the shear stress values, specially in the first loop. We attribute this difference in shear values to the commented slightly different behavior of the material in shear, which was the reason for the constitutive model proposal in Ref. [71]. 5.3. Finite element example One of the main purposes of multiaxial cyclic models is to compute stress redistribution during loading in a component subjected to multiaxial plastic cyclic loading. Then, an efficient finite element implementation of any proposed model, that allows for relatively large steps, still giving accurate predictions, is needed [5]. Accordingly, we have implemented the described implicit computational algorithm in our in-house finite element code Dulcinea, already used in simulations with other constitutive models [75], [76]. In order to verify the performance of the algorithm and, at the same time, to show a comparison with finite element simulations using the model of Ohno, we simulate a finite element example taken from the paper of Ohno et al [34], who used data from annealed OFHC copper. Ohno et al used this example to show the performance of the stress integration algorithm presented therein for the well-known model of Ohno. The simulation is performed on the finite element model of half of a singlehole plate fixed at one end and subjected to cyclic bending and a constant tension at the opposite end. The cyclic bending value uend y , which is varied from 0 mm to end 5 mm, is added after the axial tension σxx is increased from 0 MPa to 5 MPa. All the displacements are imposed in our model by means of a penalty formulation. The axial stresses are imposed by equivalent uniformly distributed loads at the nodes. Details of the example are given in Figure 12. In our simulations we obtain the hardening surfaces from the stress-strain curve shown in Figure 2 of Ref. [34] following the same procedure as in the previous examples. As in the simulations performed by Ohno, we include in this example 37
0.02 1000 Shear stress τ√3 (MPa)
Shear strain γ/√3
a 0.01
0
−0.01
b
500 0 −500
−1000 −0.02 −0.02
−0.01
0
0.01
−1000 −500
0.02
Axial strain ε
1000 Shear stress τ√3 (MPa)
1000
750
500
250
0
500
1000
Axial stress σ (MPa)
c Axial stress σ (MPa)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
d
500 0 −500
−1000 0 0
0.005 0.01 Axial strain ε
0.015
−1000 −500
0
Axial stress σ
500
1000
(MPa)
Figure 11: Experimental and simulated ramp-up 90◦ out-of-phase strain path with continuously increasing amplitude. (a) Prescribed experimental axial-shear strain path [71]. (b) Experimental axial-shear stress path [71]. (c) Axial monotonic stress-strain curve employed in the simulation, obtained from Fig. 5 of [71]. (d) Predicted stress path
38
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
a 100 Ø25
75
25
end xx
s
z
x
3.5
end
uy
y Time (step)
b end
uy
0 0.7
5
15
25
35
45
55
2.1 3.5
Figure 12: Plate with a hole, clamped at one end and subjected to an axial constant tension and a cyclic bending displacement at the other end; example taken from Ref. [34]. (a) Details of the geometry with dimensions in mm. (b) Prescribed cyclic bending displacement applied to the free end.
39
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
cyclic isotropic hardening obtained from Table 4 of Ref. [34], so all the presented capabilities are active and the efficiency of the consistent tangents may be verified in the general case. In Ref. [34] the authors made the implementation of their model in the commercial code Abaqus via user subroutine. The elements they used are the typical 20-node serendipity elements C3D20R with reduced integration to avoid volumetric locking. However, in our case, in order to include more integration points per element (important in stress concentration cases as the one at hand) and to use a mixed u/p formulation, we selected fully integrated 27/4 tri-quadratic elements with linear pressure interpolation. These are reliable elements which are known to pass the inf-sup condition and give optimal performance [77]. The finite element mesh is shown in Figure 13. We included a finer layer of elements as to compare our integration-point stress-strain results with the nodal ones in Ref. [34]. Since both the specimen and the loading pattern are symmetric, we used symmetry conditions in the finite element model, so only one half of the specimen has been modelled. In our simulations we have used the same step increment for cyclic bending as the one used in Ref. [34]. Even though this time step is too large for an accurate description of the behavior of the structure, this step size shows the robustness of the numerical implementation and facilitates the comparison with [34]. However, for the initial tension loading, to better describe the incursion from the elastic into the plastic regime and to show the effect in the convergence rate, we used 5 steps instead of the single step used by Ohno et al [34]. In Figure 13 we show the von Mises stress distribution in the upper face or the specimen (in MPa) at the end of the loading cycle. Note that the band plot results shown in the figure are unaveraged, nonsmoothed stresses extrapolated at nodes in order to show the quality of the mesh. The stress jumps at both ends of the specimen are due to the approximate approach used in prescribing the boundary conditions, and for our purposes those jumps are irrelevant in the simulation according to the Saint Venant principle. The stress-strain cyclic results shown in Figure 13c are obtained at the integration point closest to the middle point of the half circle in the direction of width and closest to the upper surface in thickness. As shown in Figure 13c, our predictions of the stress-strain hysteresis loops show the feature of cyclic hardening. These predictions are very similar to those given by Ohno et al in Figure 8b of Ref. [34] using their model. Apart from being different models (although both use a multilinear approximation obtained from stress-strain curves), the small discrepancies may also be attributed to the ratchetting effect included in Ref. [34], the different element type used, the finer mesh we have used and the actual location of the corresponding integration points, which is specially important in that zone due to the stress concentration promoted by the hole. 40
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Figure 13: Finite element mesh used for our simulations of the example given in [34] (a); elements are 27/4 u/p mixed elements (tri-quadratic interpolation in displacements and linear in pressure) with full 3 × 3 × 3 integration. Band plots correspond to nodal-unaveraged von-Mises stresses at the top layer (band-plot jumps are an estimate of the local error in the solution) (b). It is also indicated the location of node nearest to the integration point which corresponds to the stress-strain history plot (c)
41
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Iteration (1) (2) (3) (4)
Step 14 4.366E + 06 4.706E + 02 1.711E + 01 1.973E − 01
Iteration (1) (2) (3) (4)
Step 14 3.063E + 07 6.833E − 01 1.713E − 03 2.282E − 07
Residual norm Step 34 4.366E + 06 4.270E + 02 1.539E + 01 1.021E − 01 Energy norm Step 34 3.063E + 07 6.133E − 01 1.302E − 03 8.144E − 08
Step 54 4.366E + 06 4.839E + 02 1.663E + 01 2.542E − 01 Step 54 3.063E + 07 6.889E − 01 1.125E − 03 1.213E − 07
Table 3: Global equilibrium convergence values at steps 14, 34, 54 of the finite element simulation
The number of global equilibrium iterations employed for each step is shown in Figure 14. As usual, the first iteration corresponds to the imposition of the incremental loads of the step, giving the reference non-equilibrium values, whereas the remaining iterations are the actual iterative search for the equilibrium configuration. Since the first two steps are elastic, just two iterations (one for imposing load and one to regain equilibrium) give the solution to machine precision. The rest of the steps involve plastic strains, so more iterations are needed to reach the prescribed relative tolerances in forces of 10−7 and in energy of 10−14 . However note that usually only 4 iterations (3 equilibrium iterations) are needed to reach such tolerances, which is a performance similar to that of classical von Mises plasticity. Furthermore, in Table 3 we show force and energy convergence values for three typical steps. It can be observed that asymptotic quadratic convergence is obtained. 6. Conclusion In this paper we have presented a new model for multiaxial cyclic plasticity which follows the classical Prager’s formulation, including nonlinear kinematic hardening and cyclic hardening/softening. We also developed a fully implicit finite element implementation of the model. The model uses the idea from Mr´oz of discretizing the uniaxial stress-strain curve in several segments, resulting in a field of hardening surfaces. This approach is very simple in obtaining material parameters of the model. However, in contrast to Mr´oz’s proposal, our model is based on the traditional associative Prager rule and the yield surface is always the innermost surface. The 42
Number of global iterations per step
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
7 6 5 4 3 2 1 0
10
20
30
40
50
step number
Figure 14: Number of global iterations used for each loading step
outer surfaces are just hardening surfaces, a tool to describe directional and historydependent anisotropic hardening moduli. Our formulation and numerical algorithm naturally recover classical von Mises plasticity with mixed hardening in the case that kinematic hardening is linear in practice, regardless of the number of surfaces being employed. Hence, the predictions are self-consistent. Furthermore, since the model is in essence classical plasticity, the large strain implementation is straightforward and unmodified from the classical plasticity using the typical geometric pre- and post-processors in terms of logarithmic strains. We have validated the model with predictions of experimental results from both uniaxial and multiaxial tests taken from the literature. It is seen that the model is capable of giving good predictions in all cases, which include different materials and multiaxial loading paths. Finally, we have demonstrated the numerical efficiency of the fully implicit, Backward-Euler algorithmic formulation by the implementation in a finite element program and the simulation of a cyclic nonproportional loading example from the literature. Two effects are still not included in the formulation herein presented: the ratchetting effect and the path-dependent isotropic hardening. These effects can be accommodated in the present theory preserving Prager’s rule, respectively, by the inclusion of a ratchetting internal strain in the internal compatibility equation, and by allowing the full dependency of the isotropic hardening on ξ.
43
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Acknowledgments Partial funding from grant DPI2015-69801-R from the Direcci´on General de Proyectos de Investigaci´on of the Ministerio de Econom´ıa y Competitividad of Spain is acknowledged. We thank the authors of Ref. [71] for the access to the digital data of the experiment in their paper, shown in Fig. 11. FJM also acknowledges the support of the Department of Mechanical and Aerospace Engineering of University of Florida during the sabbatical period in which part of this work was done, and Ministerio de Educaci´on Cultura y Deporte of Spain for the financial support for that stay under grant PRX15/00065. References [1] Borja RI, Amies AP. Multiaxial cyclic plasticity model for clays. J Geotech Eng ASCE 1994; 120(6):1051-1070. 1 [2] Darwin D, Nmai CK. Energy dissipation in RC beams under cyclic load. J Struct Eng 1986; 112(8):1829-1846. 1 [3] Benavent-Climent, A. Seismic behavior of RC wide beam-column connections under dynamic loading. J Earthquake Eng 2007; 11(4):493-511. 1 [4] Gates NR, Fatemi A. A simplified cyclic plasticity model for calculating stressstrain response under multiaxial non-proportional loadings. Eur J Mech-A/Solid 2016; 59:344-355. 1 [5] Madrigal C, Navarro A, Chaves V. Numerical implementation of a multiaxial cyclic plasticity model for the local strain method in low cycle fatigue. Theor Appl Fract Mech 2015; 80:111-119. 1, 5.3 [6] Navarro A, Vallellano C, Chaves V, Madrigal C. A microstructural model for biaxial fatigue conditions. Int J Fatigue 2011;33(8):1048-1054. 1 [7] Teixeira P, Santos AD, Andrade Pires FM, C´esar de S´a JMA. Finite element prediction of ductile fracture in sheet metal forming processes. J Mater Proc Technol 2006;177:278-281. 1 [8] Khoei A. Computational plasticity in powder forming processes. Elsevier; 2010. 1
44
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[9] Cvitani´c V, Vlak F, Lozina Z. A finite element formulation based on nonassociated plasticity for sheet metal forming. Int J Plasticity 2008; 24(4):646687. 1 [10] Vladimirov IN, Pietryga MP, Reese S. Prediction of springback in sheet forming by a new finite strain model with nonlinear kinematic and isotropic hardening. J Mater Proc Technol 2009; 209:4062-4075 1 [11] Levkovitch V, Svendsen B. Accurate hardening modeling as basis for the realistic simulation of sheet forming processes with complex strain-path changes. AIP Conf Proc 2007; 908:1331-1336 1 [12] Milne I, Ritchie RO, Karihaloo BL (Eds). Comprehensive structural integrity: Cyclic loading and fatigue (vol 4). Elsevier; 2003. 1 [13] Li B, Reis L, De Freitas M. Simulation of cyclic stress/strain evolutions for multiaxial fatigue life prediction. Int J Fatigue, 2006;28(5):451-458. 1 [14] Ince A, Glinka G, Buczynski A. Computational modeling of multiaxial elastoplastic stress–strain response for notched components under non-proportional loading. Int J Fatigue 2014; 62:42-52. 1 [15] Rezaiee-Pajand M, Sharifian M. Accurate and approximate integrations of Drucker-Prager plasticity with linear isotropic and kinematic hardening. Eur J Mech A/Solid 2011; 30:345-361 1 [16] Moftakhar A, Buczynski A, Glinka G. Calculation of elasto-plastic strains and stresses in notches under multiaxial loading. Int J Fract 1994; 70(4):357-373. 1 [17] Zeng Z, Fatemi A. Elasto-plastic stress and strain behavior at notch roots under monotonic and cyclic loadings. J Str Anal Eng Des 2001; 36(3):287-300. 1 [18] Singh MNK, Glinka G, Dubey RN. Elastic-plastic stress-strain calculation in notched bodies subjected to non-proportional loading. Int J Fract 1986; 76(1):39-60. 1 [19] Simo JC, Hughes TJR. Computational Inelasticity. Springer; 1998. 1, 2.1, 2.1, 3.2, 3.3 [20] Lemaitre J, Chaboche JL. Mechanics of solid materials. Cambridge University Press; 1994. 1
45
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[21] De Angelis F, Taylor RL. An efficient return mapping algorithm for elastoplasticity with exact closed form solution of the local constitutive problem. Eng Comput 2015; 32(8): 2259-2291 1 [22] De Angelis F, Taylor RL. A nonlinear finite element plasticity formulation without matrix inversions. Fin Elem Anal Des 2016; 112:11-25. 1 [23] Valanis KC. Fundamental consequences of a new intrinsic time measure. Plasticity as a limit of the endochronic theory. Arch Mech 1978; 32:171. 1, 2.2 [24] Mr´oz Z. On the description of anisotropic workhardening. J Mech Phys Solid 1967; 15(3): 163-175. 1 [25] Garud YS. A new approach to the evaluation of fatigue under multiaxial loadings. J Eng Mater Technol ASME 1981; 103(2):118-125. 1 [26] Armstrong PJ, Frederick CO. A mathematical representation of the multiaxial Bauschinger effect. GEGB Report RD/B/N/731, Berkeley Nuclear Laboratories, UK; 1966. 1 [27] Chaboche JL. Time-independent constitutive theories for cyclic plasticity. Int J Plasticity 1986; 2(2):149-188 1, 2.5 [28] Dafalias YF, Popov EP. Plastic internal variables formalism of cyclic plasticity. J Appl Mech ASME 1976; 43: 645. 1 [29] Dafalias YF. The Concept and Application of the Bounding Surface in Plasticity Theory. In: Hult J, Lemaitre J (eds) Physical Non-Linearities in Structural Analysis. International Union of Theoretical and Applied Mechanics. Springer; 1981. 1 [30] Ohno N, Wang JD. Kinematic hardening rules with critical state of dynamic recovery. Part I: formulation and basic features for ratchetting behavior. Int J Plasticity 1993; 9:375–390. 1 [31] Mont´ans FJ. Implicit algorithms for multilayer J 2-plasticity. Comput Meth Appl Mech Eng 2000; 189(2):673-700. 1 [32] Caminero MA, Mont´ans FJ. An enhanced algorithm for nested surfaces plasticity using the implicit Mr´oz translation rule. Comput Struct 2006; 84(26):16841695. 1
46
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[33] Khoei AR, Jamali N. On the implementation of a multi-surface kinematic hardening plasticity and its applications. Int J Plasticity 2005; 21(9):1741-1770. 1 [34] Ohno N, Tsuda M, Kamei T. Elatoplastic implicit integration algorithm applicable to both plane stress and three-dimensional stress states. Finite Elem Anal Des 2013; 66:1-11 1, 5.3, 5.3, 12, 5.3, 13 [35] Halama R, Markopoulos A, Janco R, Bartecky M. Implementation of MAKOC cyclic plsticity model with memory. Adv Eng Soft 2017; in press, doi: 10.1016/j.advengsoft.2016.10.009 1, 5.2.3, 9 [36] Tanaka E, Murakami S, Ooka M. Effects of strain path shapes on nonproportional cyclic plasticity. J Mech Phys Solid 1985; 33(6):559-575 1, 2.1, 5.2.4, 5.2.4, 10 [37] Benallal A, LeGallo P, Marquis D. An experimental investigation of cyclic hardening of 316 stainless steel and of 2024 aluminum alloy under multiaxial loadings. Nucl Eng Des 1989; 114:345-353 1 [38] Calloch S , Marquis D. Triaxial tension–compression tests for multiaxial cyclic plasticity. Int J Plasticity. 1999; 5(5):521–49. 1 [39] Tanaka E. A nonproportionality parameter and a cyclic viscoplastic constitutive model taking into account amplitude dependences and memory effects of isotropic hardening. Eur J Mech/A Solid 1994; 13(2):155-173 1, 5.2.4 [40] Gates NR, Fatemi A. A simplified cyclic plasticity model for calculating stressstrain response under multiaxial non-proportional loadings. Eur J Mech A/Solid 2016; 59:344-355 1, 5.2.4 [41] Meggiolaro MA, Wu H, Pinho de Castro, JT. Non-proportional hardening models for predicting mean and peak stress evolution in multiaxial fatigue using Tanaka’s incremental plasticity concepts. Int J Fatigue 2016; 82:146-157 1, 5.2.4 [42] Taleb L, Cailletaud G. An updated version of the multimechanism model for cyclic plasticity. Int J Plasticity 2010; 26:859-874 1 [43] Cardoso RPR, Yoon JW. Stress integration method for a nonlinear kinematic/isotropic hardening model and its characterization based on polycrystal plasticity. Int J Plasticity 2009; 25:1684-1710 1 47
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[44] Mont´ans FJ. Implicit multilayer J2-plasticity using Prager’s translation rule. Int J Num Meth Eng 2001; 50:347-375 1 [45] Mont´ans FJ. Implicit plane stress algorithm for multilayer J2-plasticity using the Prager–Ziegler translation rule. Int J Num Meth Eng 2004; 59(3):409-418. 1 [46] Jiang Y, Sehitoglu H. Comments on the Mr´oz multiple surface type plasticity models. Int J Solid Struct 1996; 33(7):1053-1068. 1 [47] Mont´ans FJ, Caminero MA. On the consistency of nested surfaces models and their kinematic hardening rules. Int J Solid Struct 2007; 44(14):5027-5042. 1 [48] Chaboche JL. On some modifications of kinematic hardening to improve the description of ratchetting effects. Int J Plasticity 1991; 7:661-678. 1 [49] Abdel-Karim M. An extension for the Ohno-Wang kinematic hardening rules to incorporate isotropic hardening. Int J Pres Vessel Piping 2000; 87(4):170-176. 1 [50] Kossa A, Szabo L. Exact integration of the Von Mises elastoplasticity model with combined linear isotropic-kinematic hardening. Int J Plasticity 2009; 25:10831106 1 [51] Berisha B, Hora P, Wahlen A, Tong L. A combined isotropic-kinematic hardening model for the simulation of warm forming and subsequent loading at room temperature. Int J Plasticity 2009; 26:126-140 1 [52] Cao J, Lee W, Cheng HS, Seniw M, Wang HP, Chung K. Experimental and numerical investigation of combined isotropic-kinematic hardening behavior of sheet metals. Int J Plasticity 2009; 25:942-972 1 [53] Firat M. Cyclic plasticity modeling and finite element analyzes of a circumferentially notched round bar under combined axial and torsion loadings. Mater Design 2012; 34:842-852 1 [54] Hong PY, Pereira JM, Cui YJ, Tang AM, Collin F, Li XL. An elastoplastic model with combined isotropic-kinematic hardening to predict the cyclic behavior of stiff clays. Comput Geotech 2014; 62:193-202 1 [55] Aahnken R. Improved implementation of an algorithm for non-linear isotropic/kinematic hardening in elstoplasticity. Commun Num Meth Eng 1999; 15:745-754 1 48
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[56] Ohno N. A Constitutive model of cyclic plasticity with a non-hardening strain region. J Appl Mech. 1982; 49:721-727 1 [57] Xu L, Nie X, Fan J, Tao M, Ding R. Cyclic hardening and softening behavior of the low yield point steel BLY160: Experimental response and constitutive modeling. Int J Plasticity 2016; 78:44-63 1, 2.5, 5.1, 5 [58] Krieg RD, Key SW. Implementation of a time dependent plasticity theory into structural computer programs. In Constitutive Equations in Viscoplasticity: Computational and Engineering Aspects, JA Stricklin and KJ Saczalski eds, Appl Mech Div AMD-20, ASME, New York 1976. 1, 3.2 [59] Sanz MA, Mont´ans FJ, Latorre M. Computational anisotropic hardening multiplicative elstoplasticity based on the corrector elastic logarithmic strain rate. Comput Meth Appl Mech Eng 2017; 320:82-121 2.1, 2.1, 4 [60] Latorre M, Mont´ans FJ. A new class of plastic flow evolution equations for anisotropic multiplicative elastoplasticity based on the notion of a corrector elastic strain rate. Under review. arXiv:1701.00095 [cond-mat.soft] 2016 2.1, 2.1, 4 [61] Lamba HS, Sidebottom OM. Cyclic plasticity for nonproportional paths: Part 1 — Cyclic hardening, erasure of memory, and subsequent strain hardening experiments. Part 2 — Comparison with predictions of three incremental plasticity models. J Eng Mater Technol 1978; 100:93-111 5.2.2, 7, 8 [62] Zhang M, Ben´ıtez JM, Mont´ans. Capturing yield surface evolution with a multilinear anisotropic kinematic hardening model. Int J Solid Struct 2016; 81:329– 336 [63] Koji´c M, Bathe KJ. Inelastic Analysis of Solids and Structures. Springer; 2005. 2.1, 4 [64] Borja RI. Plasticity. Springer; 2013. 2.1 [65] Eterovic AL, Bathe KJ. A hyperelastic-based large strain elasto-plastic constitutive formulation with combined isotropic-kinematic hardening using the logarithmic stress and strain measures. Int J Num Meth Engrg 1990; 30(6):10991114. 4, 4
49
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
[66] Caminero MA, Mont´ans FJ, Bathe KJ. Modeling large strain anisotropic elastoplasticity with logarithmic strain and stress measures. Comput Struct 2011; 89(11): 826-843. 4 [67] Valanis KC, Wu HC. Endochronic representation of cyclic creep and relaxation of metals. J Appl Mech ASME 1975; 42:67-73. 2.2 [68] Mont´ans FJ, Caminero MA. On the consistency of nested surfaces models and their kinematic hardening rules. Int J Solid Struct 2007; 44(14): 5027-5042. 2.4 [69] Mi˜ nano M, Caminero MA, Mont´ans FJ. On the numerical implementation of the Closest Point Projection algorithm in anisotropic elasto-plasticity with nonlinear mixed hardening. Fin Elem Anal Des 2016; 121:1-17. 3 [70] Mont´ans FJ, Borja RI. Implicit J2-bounding surface plasticity using Prager’s translation rule. Int. J Num Meth Eng 2002; 55:1129:1166. [71] Madrigal C, Navarro A, Chaves V. Biaxial cyclic plasticity experiments and application of a constitutive model for cyclically stable material behaviour. Int J Fatigue 2016; 83(2):240-252 5.2.5, 11, 6 [72] Navarro A, Brown MW. A constitutive model for elastic-plastic deofrmation under cyclic multiaxial straining. Fatigue Fract Eng Mater Struct 1997; 20(5):747758 5.2.5 [73] Navarro A, Gir´aldez JM, Vallellano C. A constitutive model for elastoplastic deformation under variable amplitude multiaxial cyclic loading. Int J Fatigue 2005; 27:838-846 5.2.5 [74] Navarro A, Brown MW, Miller KJ. A multiaxial stress-strain analysis for proportional cyclic loading. J Strain Anal 1993; 28(2):125-133 5.2.5 [75] Latorre M, Mont´ans FJ. Fully anisotropic finite strain viscoelasticity based on a reverse multiplicative decomposition and logarithmic strains. Comput Struct 2016; 163:56-70 5.3 [76] Crespo J, Latorre M, Mont´ans FJ. WYPIWYG hiperelasticity for isotropic, compressible materials. Comput Mech 2017; 59:73-92. 5.3 [77] Bathe KJ. Finite Element Procedures, 2nd Ed. Klaus-J¨ urgen Bathe; 2014. 5.3
50