INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING Int. J. Numer. Meth. Engng. 47, 75–100 (2000) Consistent return mapping algorithm for chemoplastic constitutive laws with internal couplings J. Sercombe 1 , F.-J. Ulm 2 and H.-A. Mang1;∗ 2 1 Institute for Strength of Materials; Vienna University of Technology; Austria Concrete Department; Laboratoire Central des Ponts et Chaussees; Paris; France SUMMARY This paper discusses the algorithmic formulation of chemoplastic models with internal couplings. Developed within the framework of closed reactive porous continua, chemoplasticity allows for the introduction of the kinetics of a chemical reaction directly at the macro-level of material description. The resulting couplings (chemoelastic, chemoplastic) slightly complicate the integration of stresses at the material level, leading to chemoplastic hardening =softening. In this paper, the numerical integration of stresses for this type of constitutive laws with chemomechanical (internal) couplings is presented. Based on the denition of return mapping algorithms for standard plasticity, a closest-point projection algorithm is developed, which preserves the quadratic rate of convergence of the Newton–Raphson iterative scheme. Examples are presented with an increasing degree of complexity: perfectly plastic behaviour with chemical hardening for concrete at early ages, coupled creep-plastic behaviour, plastic behaviour with combined viscoplastic hardening =softening for concrete in high-rate dynamics. Copyright ? 2000 John Wiley & Sons, Ltd. KEY WORDS: chemoplasticity; return mapping; internal couplings; kinetic law PREFACE This paper is dedicated to the memory of the late Professor Richard H. Gallagher. The third author considers his scientic association with Professor Gallagher from 1975 till his death as a great privilege. He wishes to express his deep gratitude not only to a distinguished scientist and a genuine mentor but also to a great man and a dear friend. INTRODUCTION Modelling of granular materials such as concrete requires to account for several phenomena which may depend on the type of loading. These phenomena include but are not restricted to compaction, cracking, rate eects in dynamics and creep. Standard constitutive laws commonly used for the modelling of material behaviour (plasticity, viscoelasticity, viscoplasticity) have been applied with some success to the modelling of these phenomena (plasticity for cracking [1] or compaction [2], ∗ Correspondence to: H.-A. Mang, Institut fur Festigkeitslehre, Technical University of Vienna, Karlsplatz 13=202, A-1040 Vienna, Austria. E-mail: Herbert.Mang@tuwien.ac.at Contract=grant sponsor: Austrian Foundation for the Advancement of Scientic Research CCC 0029-5981/2000/010075–26$17.50 Copyright ? 2000 John Wiley & Sons, Ltd. Received 24 September 1998 76 J. SERCOMBE, F.-J. ULM AND H.-A. MANG visco-plasticity for rate eects in dynamics [3; 4], aging viscoelasticity for creep of concrete [5]). Nevertheless, for certain cases, phenomena of chemophysical nature have to be accounted for to model complex material behaviour. For a material such as concrete, this need mainly arises from its porous space in which dierent phenomena can take place: chemical reactions during the hydration process which lead to progressive hardening and shrinkage of the material, alkali-silicate reaction with the aggregates resulting in progressive softening and swelling of the material [6], leaching of the cement paste yielding progressive degradation of the material [7], viscous eects in the nanopores when concrete is subjected to high loading rates [8] resulting in an increase of strength, dehydration of concrete under high temperatures leading to progressive softening and damage of the material [9]. These limitations brought up the development of the theory of reactive porous media [10; 11] allowing for the introduction of the kinetics of a chemical reaction directly at the macrolevel of material description. As a particular case, it was used to model chemomechanical couplings within a closed reactive porous continuum [11; 12]. Within this specic framework, the complex chemomechanical behaviour of materials can be described by the cross-eects between the state variables of the model. By analogy with standard plastic hardening=softening, chemomechanical couplings have led to the denition of new concepts such as chemical or viscous hardening =softening, which basically allows one to describe the dependence of strength on a chemical variable. These new concepts slightly complicate the integration of stresses at the material level since generally, a kinetic constraint function (related to chemical evolutions) together with one yield criterion (related to plastic evolutions) are necessary and, furthermore, their positions in the stress space depend on both plastic variables (leading to standard plastic hardening=softening) and chemical variables (leading to chemical hardening=softening). In this paper, the numerical integration of stresses for this type of constitutive laws with chemoplastic hardening=softening is presented. Based on the denition of return mapping algorithms for standard plasticity [13], a closest-point projection algorithm is developed. The consistent tangent stiness matrix is derived in order to preserve the quadratic rate of convergence of the Newton– Raphson iteration scheme. In the rst part of this paper, the basic equations for chemoplastic constitutive laws with internal couplings are presented and the main dierences as compared to elastoplastic models are described. In the second part, the development of the closest-point projection algorithm is detailed, with particular emphasis on the denition of the trial stress state. In the third part of this paper, some applications to the modelling of concrete behaviour at early ages [12; 14; 15] and under high loading rates [16] are presented. In between, to demonstrate its generality, this algorithm is applied to the integration of stresses for a combined plasticity-creep model [17] which can be used for a great variety of materials. Finally, some numerical examples are shown which illustrate the coupled nature of these models. CONSTITUTIVE EQUATIONS Thermodynamic formulation of chemoplasticity with internal couplings Let us consider the framework of closed reactive porous media [11]. Closed refers to the elementary system which is considered as closed for the chemical constituents, i.e. there is no external Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS 77 supply. This allows to reduce the number of state variables related to the chemical process which therefore become internal state variables. Considering the free energy of the system, it reads as a function of the following variables: U the total strain tensor, Up , the plastic strain tensor and a hardening =softening variable (for clarity, we will restrict the derivation to one plastic hardening =softening mechanism). For granular materials, these variables usually account for the behaviour of the hardened part of the material (the skeleton deformation and cracking). On account of a chemical reaction occurring within the porous space of the material, an internal state variable (of scalar nature, which reects the isotropy of the internal chemical process) is introduced in the constitutive modelling. In a closed system, the variable is a measure of the reaction extent, i.e. of the progress of the chemical reaction considered [11]. Applying standard thermodynamics, the intrinsic dissipation is given by ˙ ¿0 ’ = b : U̇ − (1) where is the free energy of the material system which in this case reads = (U; Up ; ; ) = (U − Up ; ) + U (; ) (2) In (2), is the elastic strain energy, and U (; ) the ‘frozen’ energy. Use of (2) in (1) yields the state equations of the model: b= @ ; @(U − Up ) =− @U ; @ A=− @( + U ) @ (3) where b is the stress tensor, associated in dissipation (1) to plastic evolutions U̇p , is the hardening force, associated in dissipation (1) to plastic hardening =softening evolutions ˙ and A is the thermodynamic imbalance (chemical anity) which governs the evolution of the chemical reaction modelled by the internal variable . For the sake of clarity, it is useful to introduce the Legendre–Fenchel transform [10] of elastic strain energy , dened by ∗ = b : (U − Up ) − ; ∗ = ∗ (b; ) (4) which yields the state equations in the following form: U − Up = @ ∗ ; @b =− @U ; @ A= @( ∗ − U ) @ (5) In order to discuss and dene adequately the couplings between the state variables, the state equations are better formulated in dierential form on account of Maxwell symmetries [12; 14; 16]. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 78 J. SERCOMBE, F.-J. ULM AND H.-A. MANG Deriving (5), the following system, expressed in matrix form, is obtained: 2 ∗ @ @2 ∗ @2 ∗ @b2 @b@ @U@ db p dU − dU 2 2 2 @ @ @ U U U d d = − − 2 − @ @@ @@b d dA @2 ∗ ∗ 2 2 @ U @ ( − U) − @@b @@ @2 (6) Chemoelastic couplings Considering the rst equation of system (6), the following expression can be formulated: dU − dUp = : db + b d (7) wherein the following denitions were introduced: 1. = @2 ∗ =@b2 = C−1 is the inverse of the tensor of elastic moduli C linking the strain tensor dU − dUp to the stress tensor db. 2. b = @2 ∗ =@b@ is a second-order tensor which relates the increase of internal variable to the strains of chemical origin (dUc = b d). The tensor b is not necessarily a set of constant coecients. It can depend on the stress state b = b(b). 3. @2 ∗ =@b@ = 0 since the Legendre–Fenchel transform of the elastic strain energy, ∗ (b; ), does not depend on the hardening variable . Note that this term denes a plastic-damage coupling, which reects the possible dependence of the tensor of elastic moduli on the plastic state of the material C() [14; 18]. Integration and inversion of equation (7) yields the state equation of the chemoplastic model with internal couplings in the form b = C : (U − Up − Uc ) with dUc = b(b) d (8) Chemoplastic couplings The second equation of system (6) can be recast as follows: d = @ @ d + d @ @ (9) where the following denitions were introduced: (i) @=@ = − @2 U=@2 is the tangent modulus for plastic hardening =softening, relating the hardening force d to the increase of plastic hardening =softening variable d. (ii) @=@ = − @2 U=@@ is the tangent modulus for chemical hardening =softening (by analogy with plastic hardening=softening), relating the hardening force d to the extent of the internal chemical process d. (iii) @2 U=@@b is the counterpart of the plastic damage coupling previously considered (@2 ∗ = @b@ = 0). On account of Maxwell symmetries and because of the independence of the Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS 79 frozen energy U (; ) of the stress tensor b, the respective term in (6) does not appear in (9). Integration of (9) yields the hardening force = (; ). Hence, the hardening force which denes the elasticity domain of the material, does not only depend on the plastic hardening=softening variable but also on the reaction extent . The latter dependence leads to the concept of chemical hardening=softening, analogous to the classical concept of plastic hardening=softening. As in standard plasticity, the elastic domain is specied by a loading function f in the stress space: f = f(b; ) 6 0 (10) From equation (10), it is clear that the position of the yield surface (f = f(b; ) = 0) in the stress space depends on the reaction extent via the hardening force (; ). Therefore, the strength of any material (modelled by the plastic yield surface) can be related to any chemical process occurring in the material (modelled by the variable ). The instantaneous character of plastic evolutions (occurring only if f = 0 and df = 0, as in standard plasticity) is, however, preserved even when chemoplastic hardening =softening occurs. In this respect, standard evolution laws can still be used for the plastic strain tensor Up and the plastic hardening variable : dUp = d @g ; @b d = d @h @ (11) where g and h denote the plastic and hardening potentials, respectively, and d is the plastic multiplier. Macroscopic kinetic law On account of Maxwell symmetries, the third equation of system (6) can be recast in the following form: dA = b : db + @2 U @2 ( ∗ − U ) d d − 2 @ @@ (12) wherein: (a) b = @2 ∗ =@@b = @2 ∗ =@b@ is the counterpart of the chemoelastic coupling previously considered. It accounts for the eect of stress on the chemical process. (b) @2 ( ∗ − U )=@2 is a coecient of equilibrium for the chemical process, the evolution of which is governed by the chemical anity A. (c) @2 U=@@ is a coecient of equilibrium related to the chemical hardening=softening process. It accounts for the eect of plastic hardening=softening evolutions (d) on the chemical process (A). Integration of equation (12) leads to the following general expression of the chemical anity A: A = A(b; ; ) (13) Since the chemical anity A models, at the macroscopic level, the thermodynamic imbalance governing chemical evolutions, it has to be related to the rate of the chemical variable ˙ through Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 80 J. SERCOMBE, F.-J. ULM AND H.-A. MANG the denition of a suitable kinetic law. This relation is the key point of the constitutive modelling of any internal chemical process. It is usually obtained from chemo-physical considerations. This law depends on the nature of the process, whether it is inuenced by the stress state (b), the plastic state () or its own state (). Since the anity A depends on these quantities (see (13)), the chemical rate ˙ can simply be expressed as a function of chemical anity A: ˙ = f(A) = F(b; ; ) (14) CLOSEST POINT PROJECTION ALGORITHM Discretized form of the constitutive equations By analogy with Closest-Point Projection formulations for plasticity, the implicit backward-Euler dierence scheme will be applied to transform the continuous problem into a discrete constrained optimization problem [13]. The discretized version of the state equation (8), the ow rule (11) and the hardening rule (11) for the chemoplastic model with internal couplings are given hereafter: Un+1 = Un + ∇s (u) bn+1 = C : (Un+1 − Upn+1 = Upn + n+1 n+1 = n + n+1 (standard) Upn+1 Ucn+1 ) − @g @b n+1 @h @ n+1 Ucn+1 = Ucn + n+1 bn+1 (15) (16) (17) (18) (19) Classically, the discrete counterparts of the Kuhn–Tucker conditions have to be enforced: f(bn+1 ; n+1 ) 6 0 n+1 ¿ 0 (20) n+1 f(bn+1 ; n+1 ) = 0 with n+1 = (n+1 ; n+1 ). By analogy, it is necessary to dene limiting conditions for the evolution of the chemical variable. They are obtained from the discretized version of kinetic law (14) by considering a constant rate ˙ within the time step (rst-order integration scheme). Using the same backward-Euler scheme, the following condition can be obtained: F(bn+1 ; n+1 ; n+1 ) =−n+1 + tF(bn+1 ; n+1 ; n+1 ) = 0 (21) Contrary to standard plasticity, no conditions are imposed on the chemical increment n+1 which can either be positive or negative. Equation (21) can be considered, algorithmically, as the constraint condition for chemical evolutions: a ‘yield function’ with zero threshold. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS 81 Denition of the trial stress state Condition (21) implies in turn that for any stress state, even in the elastic domain, chemical evolutions will take place. Therefore, the standard denition of the trial elastic stress state, obtained by freezing plastic ow, and used as a starting point for the return mapping scheme, is not valid any more. Since chemical and plastic evolutions are not independent, it is necessary to freeze both plastic ow and chemical evolutions to obtain the trial elastic stress state: trial = C : (Un+1 − Upn − Ucn ) bn+1 (22) and the trial elastic state variables: trial n+1 = n trial n+1 = n (23) trial n+1 = n With respect to standard plasticity, this denition of the trial elastic stress state is not sucient to determine whether plastic evolutions should be considered or not. In fact, since chemical evolutions take place even if the stress state remains in the elastic domain, the trial elastic stress state is not the solution of the chemoelastic problem. Furthermore, since the plastic criterion depends on the chemical variable n+1 , its position in the stress space can evolve even if the loading point remains trial lies outside of the surface; chemoelastic. This in turn means that even if the trial elastic stress bn+1 chemical hardening=softening can bring it back inside the loading surface. Therefore, in order to check whether plastic evolutions should occur or not, it is necessary to introduce an intermediate step (see Figure 1 for a schematic representation of the successive steps in the b– plane) which consists of determining the chemoelastic solution assuming that the loading point will remain inside the plastic yield surface. This solution is obtained by calculating the following chemoelastic trial stress state: ce = C : (Un+1 − Upn − Uce bn+1 n+1 ) c ce with Uce n+1 = Un + bn+1 (n+1 − n ) (24) and the chemoelastic trial hardening force: ce = (n ; ce n+1 n+1 ) (25) which account only for the evolutions of the chemical variable . As will be shown later, this chemoelastic trial state is a particular solution of the coupled chemoplastic problem. Therefore, it will be explained in detail after the presentation of the general solution procedure. With this intermediate trial solution (24), (25), the plastic yield condition can be checked: ce ce ; n+1 )¿0 ⇔ n+1 ¿0 f(bn+1 (26) If this condition is fullled, plastic evolutions need to be taken into account. Since plastic and ce (24) cannot be chemical evolutions are not independent, the chemoelastic trial stress state bn+1 used as the starting point of the projection scheme for the chemoplastic problem. Instead, the trial trial (22) together with the trial elastic state variables (23) are used. elastic stress state bn+1 Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 82 J. SERCOMBE, F.-J. ULM AND H.-A. MANG Figure 1. Schematic representation of successive steps necessary for the updating of the stresses in a chemoplastic model with chemoplastic hardening =softening General solution for the chemoplastic problem Updating of the variables of the model is then performed within the framework of the closest point projection algorithm. Basically, the main features of the return mapping algorithm for chemoplastic models with chemoplastic hardening =softening can be explained in the simpler context of Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS 83 perfect plasticity ( = 0). In this case, the system to be solved at each iteration step (k) can be formulated as follows: (k) (k) (k) p(k) (k) @g (k) p − Ucn+1 + Ucn + (k) (27) Rn+1 =−Un+1 + Un + n+1 n+1 bn+1 = 0 @b n+1 with the following two constraint conditions: (k) (k) fn+1 = f bn+1 ; (k) =0 n+1 (k) (k) = F bn+1 ; (k) Fn+1 n+1 = 0 (28) Because of the coupled nature of the problem, equations (17) and (19) were combined into one set of equations (27). Therefore, the size of the set of equations remains the same as in standard (k) plasticity. It is only increased by one equation (28) resulting from the constraint condition Fn+1 =0 on chemical evolutions. Since Un+1 is xed during the return mapping stage, it follows that (k) (k) dUin+1 = dUpn+1 + dUcn+1 = −C−1 : dbn+1 (k) (k) (29) (k) in which dUin+1 stands for the innitesimal increment of inelastic strains. Linearization of (27) together with the use of equation (29) yields the following system: R(k) n+1 h + (k) n+1 i−1 (k) : dbn+1 (k) (k) @g @2 g @ (k) + d + b + d(k) = 0 @b n+1 n+1 @b@ @ n+1 n+1 (30) with (k) n+1 as the Hessian matrix of the system, given by " (k) n+1 = C−1 + (k) n+1 (k) #−1 (k) @2 g (k) @b + n+1 @b2 n+1 @b n+1 (31) Note that for complex plastic potentials g, the second-order derivative @2 g=@b@ appears in equation (k) can be expressed as a (30) in consequence of chemical hardening =softening. From (30), dbn+1 (k) (k) function of dn+1 and dn+1 , i.e. the two unknowns of the problem: ! (k) (k) h i @g @2 g @ (k) (k) (k) (k) (k) dbn+1 =− n+1 : Rn+1 + d + b + d (32) @b n+1 n+1 @b@ @ n+1 n+1 (k) To solve for d(k) n+1 and dn+1 , the two constraint conditions (28) must be linearized: (k) (k) @f @f @ (k) (k) fn+1 + : dbn+1 + d(k) = 0 @b n+1 @ @ n+1 n+1 (k) (k) @F @F (k) (k) Fn+1 + : db + d(k) = 0 n+1 @b n+1 @ n+1 n+1 Copyright ? 2000 John Wiley & Sons, Ltd. (33) (34) Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 84 J. SERCOMBE, F.-J. ULM AND H.-A. MANG Substituting (32) into (33) and (34) results in a system of two linear algebrac equations for d(k) n+1 and d(k) : n+1 ( ) ( ) 11 12 d(k) 1 n+1 (35) = 21 22 d(k) 2 n+1 with the following expressions for the scalars 11 ; 12 ; 21 ; 22 : (k) h i @g (k) @f (k) 11 = : n+1 : @b n+1 @b n+1 (k) h (k) (k) i @f @f @ @2 g @ (k) 12 = : n+1 : b + − @b n+1 @b@ @ n+1 @ @ n+1 (36) (k) h i @g (k) @F (k) 21 = : n+1 : @b n+1 @b n+1 (k) h (k) i @F @F (k) (k) 22 = : n+1 : bn+1 − @b n+1 @ n+1 and for the scalars 1 ; 2 : (k) 1 = fn+1 − (k) h i @f (k) : (k) n+1 : Rn+1 @b n+1 (k) h i @F (k) (k) (k) 2 = Fn+1 − : n+1 : Rn+1 @b n+1 (37) The system of algebrac equations (35) is solved at each iteration step (k) for (dkn+1 ; dkn+1 ). Then the stresses are updated using equation (32). Finally, the dierential increment of the inelastic i(k) strains dn+1 is obtained from equation (29). Note that plastic and chemical strains can only be calculated independently after convergence is achieved (when |R(k) n+1 | ≈ 0). Local iterations are (k) (k) | and criteria |fn+1 | and |Fn+1 | become performed in this way until the norm of the residuals |R(k) n+1 less than prescribed tolerances. Therefore, the application of a closest-point projection algorithm to the integration of the stresses for a chemoplastic model with internal couplings does not lead to additional diculties with respect to its application to standard plasticity. The main modications lie in: (1) the derivation of the Hessian matrix (31) which includes the gradients of the second-order tensor b, k ; dkn+1 ) at each iteration by means of solving a system of two (2) the solution for (dn+1 algebrac equations. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS 85 General solution for the chemoelastic problem As a special case of the system of algebrac equations (35), the chemoelastic solution d(k) n+1 is (k) obtained by solving the second equation under consideration of the condition dn+1 = 0: d(k) n+1 = 2 22 (38) Once the increment of the internal variable d(k) n+1 is determined, the stresses can be updated from (k) = 0 and = 0. This particular solution leads the chemoelastic trial equation (32) with d(k) n+1 n+1 ce dened by equation (24). It corresponds to the solution for the case that the stress stress state bn+1 point remains in the elastic domain of the material, as illustrated in Figure 1. General solution for the chemoplastic problem with chemoplastic hardening =softening Finally, the details of the general case with chemoplastic isotropic hardening=softening are described hereafter and illustrated in Box 1. (i) The residuals (27) include h an iadditional equation related to the hardening law (18). (ii) The generalized moduli A(k) n+1 (39) remains the same as for the case of the integration of elasto-plastic constitutive laws with plastic hardening =softening only (the Hessian matrix i h (k) n+1 being dened however, by equation (31)). This is a consequense of the fact that the ow rule and the hardening rule still hold without any modications in the formulation of chemoplastic models with internal couplings. The only dierence with respect to (k) which cannot be directly dened standard plasticity originates from the expression of dn+1 (k) since chemical hardening =softening leads to some evolutions of the as a function of dn+1 hardening force : ! (k) (k) @ @ (k) (k) (k) −1 d (40) with Dn+1 = − dn+1 =−Dn+1 dn+1 − @ n+1 n+1 @ n+1 (iii) A similar but more elaborated system with respect to (35) is obtained for the calculation (k) of d(k) n+1 and dn+1 . The complete expression is given in Appendix I. Derivation of the consistent tangent moduli Rewriting the stress–strain relations (16), the discrete ow rule (17) and the discrete evolution law (19) in dierential form, with attention restricted to the case of perfect plasticity, the following system is obtained: dbn+1 = C : (dUn+1 − dUpn+1 − dUcn+1 ) with dUpn+1 = n+1 dUcn+1 = n+1 @g @2 g @2 g @ : dbn+1 + dn+1 + n+1 dn+1 @b2 n+1 @b n+1 @b@ @ n+1 @b : dbn+1 + b(bn+1 ) dn+1 @b (41) (42) n+1 Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 86 J. SERCOMBE, F.-J. ULM AND H.-A. MANG Box 1. General closest-point projection iteration (0) (0) (0) 1: Initialize: k = 0; Upn+1 = Upn ; Ucn+1 = Ucn ; n+1 = 0; (0) n+1 = 0; n+1 = 0 (0) 2. Check yield conditions and evaluate ow rule=hardening law residuals: (k) (k) (k) fn+1 = f(bn+1 ; n+1 ) (k) (k) (k) Fn+1 = F(bn+1 ; n+1 ; (k) n+1 ) ( R(k) n+1 = (k) −Uin+1 + Uin (k) p −n+1 + np ) (k) @g ( (k) ) @b bn+1 n+1 (k) (k) + n+1 + n+1 (k) @h 0 @ (27) n+1 (k) ¡TOL1 ; If: fn+1 (k) Fn+1 ¡TOL2 and |R(k) n+1 |¡TOL3 3. Compute consistent tangent moduli: h i−1 (k) n+1 i−1 h = A(k) (k) n+1 @2 h (k) n+1 @@b n+1 then EXIT (k) (k) @2 g @b@ n+1 # (k) 2 (k) @ h −1 Dn+1 + n+1 2 @ n+1 (k) n+1 4. Obtain increment of plastic multiplier and reaction extent: #−1 ( ) ( (k) ) " 11 12 dn+1 1 = (k) 2 21 22 dn+1 (39) (35) 5. Obtain the incremental stress tensor and hardening force (k) @g (k) ( (k) ) bn+1 @b h i dbn+1 n+1 (k) (k) (k) (k) (k) =− An+1 Rn+1 + dn+1 (k) + dn+1 @ (k) −1 dn+1 −D @h @ n+1 @ n+1 6. Update state variables and associated forces (k+1) (k) (k) bn+1 = bn+1 + dbn+1 (k+1) (k) (k) = + d n+1 n+1 n+1 (k+1) = (k) + d(k) n+1 n+1 n+1 (k+1) (k) (k) n+1 = n+1 + dn+1 ! (k) (k+1) (k+1) @ (k) (k) (k) −1 (k) Ui −1 i(k) n+1 = n+1 + Dn+1 dn+1 − @ dn+1 n+1 = Un+1 −Cn+1 : dbn+1 n+1 Set k = k + 1 and GOTO 2. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 87 CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS Using equations (42) in equation (41) leads to the algorithmic relation @g @2 g @ d − b + d dbn+1 = n+1 : dUn+1 − n+1 n+1 @b n+1 @b@ @ n+1 (43) where n+1 is the algorithmic moduli, previously dened by equation (31). Substituting (43) into the discrete consistency condition: @f @f @ dfn+1 = : dbn+1 + dn+1 = 0 (44) @b n+1 @ @ n+1 and the analogous discrete condition for chemical evolutions: @F @F dFn+1 = : db + dn+1 = 0 n+1 @b n+1 @ n+1 (45) yields a system of two algebrac equations for dn+1 and dn+1 : @f ( )( ) : [n+1 ] : dUn+1 @b n+1 11 12 dn+1 = @F 21 22 dn+1 : [n+1 ] : dUn+1 @b n+1 (46) wherein (11 ; 12 ; 21 ; 22 ) are still given by (36). From this system, dn+1 and dn+1 can be expressed as functions of dUn+1 : 1 @f @F dn+1 = : n+1 : dUn+1 (47) 22 −12 n+1 @b @b n+1 1 @F @f : n+1 : dUn+1 (48) 11 −21 dn+1 = n+1 @b @b n+1 with n+1 ; the determinant of the (11 ; 12 ; 21 ; 22 ) matrix. Finally, substitution of (47) and (48) into (43) yields the expression for the algorithmic tangent moduli for chemoplastic models with internal couplings: 1 @g @f @F db = n+1 − : n+1 ⊗ n+1 : 22 −12 dU n+1 @b @b @b n+1 n+1 − 1 n+1 @F @f 11 −21 @b @b n+1 n+1 : n+1 ⊗ n+1 : @2 g @ b + @b@ @ (49) n+1 APPLICATIONS In this part, constitutive models with internal couplings are presented. The number of internal couplings is progressively increased. A simple example, consisting of a square plate (of dimensions 1 × 1) subjected to a state of plane stress is considered for each of the three described applications. Is is discretized by means of four linear nite elements and analysed with various boundary and loading conditions, as illustrated in Figure 2. Attention will be focused on the couplings considered and their eect on the response of the material. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 88 J. SERCOMBE, F.-J. ULM AND H.-A. MANG Figure 2. Boundary conditions for the three examples: (1) chemoplasticity; (2) combined creep-plasticity; (3) viscous hardening plasticity Chemoplasticity for concrete at early ages Formulation. The rst application concerns the modelling of concrete behaviour at early ages [12; 14], that is to say, soon after casting. In this particular case, hydration of the cement particles in concrete is the important internal process (a chemical reaction) which governs the mechanical response of the material, via its mechanical characteristics (elastic modulus, strength) which depend on the amount of hydrated cement particles. In order for the chemical reaction to proceed, water in the concrete pores has to reach the unhydrated cement particules. Therefore, the hydration rate is governed by this micro-diusion of water which has to by-pass the already formed Cement Silicate Hydrates (CSH) to reach the unhydrated particles. The mechanical characteristics of concrete depend on the amount of already formed CSH (solid phase of the material). Hence, they depend on the mass of water (internal variable m) already combined in the solid phase. In the model, this dependence is dened through chemical hardening, that is to say the dependence of the hardening force on the mass of water already combined, m. For convenience, the internal state variable m can be replaced by the hydration degree = m=m∞ ; which varies between 0 (no water combined in the solid phase, m = 0) and 1 (all the water is combined in the solid phase m = m∞ , where m∞ is the total mass of water available for hydration). For the given case, the system of equations (6), describing the couplings between the state variables of the chemoplastic model, takes the following form: 1 p dU−dU E() G d = 0 dA 0 0 0 0 s @E()=@ ∞ 1− b (1 − s ) E()2 ft∞ 1 − 0 −A0 db d d (50) This system of equations is based on the following assumptions: (1) Poisson’s ratio does not depend on the hydration state . Only the elastic modulus E() depends on the hydration process (G = E is independent of elastic modulus E). s =(1−s )1+(1−0 )=(E∞ h−0 i2 )b) accounting on the (2) Linear chemoelastic couplings (b = ∞ s s =(1 − s )1 d with ∞ ¡0 as the one hand for chemical shrinkage (chemical strains dUc = ∞ volumetric shrinkage strain at complete hydration, s , the shrinkage threshold, and 1 the unit tensor) and, on the other hand, for the chemical stiening of the material (E() = E∞ h−0 i= (1 − 0 ) with E() as the elastic modulus depending on the hydration degree ; E∞ as the elastic modulus at complete hydration, and 0 as the percolation threshold). Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 89 CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS Table I. Material parameters for the chemoplastic material E∞ (GPa) 40 ft∞ (MPa) Ea =R (K) T (K) 0 s ∞ s A0 (1 =s) 0·2 2 4000 293 0·05 −0·0003 0·3 5 (3) Linear chemo-hardening (@=@ = ft∞ =(1−0 ) with ft∞ as the tensile strength at complete hydration) associated with a tension-cut-o criterion: I1 −() = 0 (51) wherein I1 = tr(b) is the rst invariant of the stress tensor. No plastic hardening =softening is considered in this case. Associated plasticity is assumed. (4) First-order reaction kinetics (@A=@ = −A0 with A0 ¿0 as a constant parameter). The couplings @A=@b (inuence of the stress state on the hydration reaction) and @A=@ (inuence of cracking on the hydration reaction) are neglected relative to @A=@; accounting for the inuence of the hydration state on the hydration reaction [12]. Integration of the system of equations (50) on the basis of these asumptions leads to the following state equations: U − Up = 1 s h−s i G : b + ∞ 1 E() (1−s ) () = ft∞ h−0 i (1−0 ) (52) A() = A0 (1−) To complete the system of equations (52), a macroscopic kinetic law (14) needs to be chosen. For the hydration reaction, it is well known that an Arrhenius-type law can be adopted [14], reading Ea Ea ˙ = A() exp − = A0 (1−) exp − (53) RT RT where Ea is the hydration activation energy, R the universal constant for ideal gas and T the absolute temperature. As can be seen from (53), the kinetic law for the hydration reaction does not depend on the solution of the mechanical problem (stress state and plastic variables). Therefore, in this application, the chemomechanical problem is uncoupled and the chemical solution can be determined prior to the mechanical analysis [15; 19]. Complete expressions for the system (35) and the consistent stiness matrix (49) for this application are given in Appendix II. Chemoplastic analysis of the plate. In this example, the square plate is clamped at the four edges (see case (1) in Figure 2). No external loading is applied to the plate to show that, depending on the boundary conditions, chemical shrinkage of the material (resulting from the hydration reaction) can lead to damage (plasticity) of this simple structure. The material parameters used in the analysis are given in Table I. Discretization in time of the problem is directly linked to the hydration process. The time step has to be suciently small to ensure a realistic evaluation of the hydration degree. In this analysis, Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 90 J. SERCOMBE, F.-J. ULM AND H.-A. MANG Figure 3. Evolution of the hydration degree and the tensile stress y (left) and of the chemical strain yc and the plastic p strain y (right) a constant time step of 1·0 h was applied during 200 h. The results, showing the evolution of plastic and chemical strains (yp ; yc ), together with the hydration degree and the tensile stress y ; are given in Figure 3. Chemical shrinkage (yc ) starts at t = 18 h. Since the plate is clamped on four edges, this shrinkage leads to tensile stresses (y ). Once the plastic yield stress ft () is exceeded, plasticizing occurs (yp ). These plastic strains, however, do not inuence the chemical process (see the evolution of the hydration degree or chemical shrinkage yc ) which leads to a progressive strengthening of the material (see the increase of tensile stress y as compared to the evolution of the hydration degree ). This is the result of the decoupling between the chemical evolutions and the mechanical behaviour. Structural applications of chemoplasticity can be found in [20; 21] (chemoplasticity for steelbre reinforced concrete applied to the structural optimization of tunnel-platforms for the French national railway company SNCF) and [19; 22] (chemoplasticity for shotcrete applied to the 2-D structural analysis of tunnels driven according to the New Austrian Tunnelling Method). Combined plasticity-creep model Formulation. This application deals with the integration of stresses for a combined plasticitycreep model applicable to a variety of material response (metals under high temperatures, polymers, : : :). In the framework of closed reactive porous media, creep can be modelled as a chemomechanical coupling [11]. It is this approach that will be considered in the following. Chemical hardening =softening will, however, not be considered to show that the algorithm developed in this paper is equivalent to previously dened return mapping schemes for this type of coupled problems [17; 23; 24]. For the modelling of coupled creep-plastic evolutions, system (6) takes the following form: ( p dU−dU dA Copyright ? 2000 John Wiley & Sons, Ltd. ) = s |s| s ( ) | s | db d − (54) Int. J. Numer. Meth. Engng. 47, 75–100 (2000) CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS 91 The derivation of this system is based on the following assumptions: c (a) A non-linear chemo-elastic coupling leading to deviatoric √ creep strains dU = (s=| s |)d (with s = b−I1 =31 as the deviatoric stress tensor and | s | = s : s as the second invariant of this tensor). (b) The creep evolutions are associated with a chemical reaction of extent and related to the application of deviatoric stresses @A=@b = s=| s | and to the current state of the chemical process @A=@ = −, with ¿0 as a constant parameter. (c) No hardening (of plastic or chemical origin) is considered (d = 0). Therefore, the system of equations (6) can be reduced to a 2 × 2 system (54). Integration of the system (54) leads to the following state equations for the combined creepplasticity model: U−Up = : b + Uc with dUc = s d |s| A(b; ) = | s | − (55) As a rst approximation, the kinetic law for creep can be chosen as linear such that the standard ‘strain-hardening creep’ model [17] is recovered: 1 1 ˙ = A(b; ) = (| s | − ) (56) where and are parameters related to the time-scale and magnitude of creep. They can also depend on the temperature. Finally, for the plastic evolutions, associated ow is formulated based on the von Mises criterion assuming a perfectly plastic behaviour: r f = | s| − 3 ft = 0 2 and dUp = s d |s| (57) In (57), ft is the tensile strength of the material. With regard to the previous application, the chemomechanical problem for this combined plasticity-creep model is coupled since the solution of the chemical problem (56) depends on the mechanical problem (namely on the stress state) and vice versa (55). Note, however, the symmetry of the system (54) resulting from the choice of the chemoelastic couplings. This also leads a symmetric consistent tangent stiness matrix (49). Complete expressions for the system (35) and the consistent tangent stiness matrix (49) for this application are given in Appendix III. Combined plasticity-creep analysis of the plate. In this analysis, the plate is only clamped on two edges (see case (2) in Figure 2). Displacements of equal magnitude (increments of u imp = v imp = 0·001 with a time step of 0·5 h) in the x- and y-direction are prescribed on the remaining two edges such that it results in a biaxial tensile stress state in the plate. The material parameters are given in Table II. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 92 J. SERCOMBE, F.-J. ULM AND H.-A. MANG Table II. Material parameters for the combined creep-plasticity model E (GPa) 200 ft (MPa) (GPa.h) (GPa) 0·2 500 40 20 p Figure 4. Evolutions of plastic and creep strains, y ; yc (left), and of the tensile stress y (right) The numerical results illustrating the combination of creep and plasticizing of the plate are shown in Figure 4 where the evolutions of the creep strain yc and of the plastic strain yp are illustrated, together with the tensile stress y . Creep of the plate occurs as soon as the plate is loaded (see the creep strain yc ). It results in a pronounced non-linear evolution of the tensile stress y . Once the plastic yield stress is reached, the tensile stress y remains equal to the yield stress ft , because perfectly plastic behaviour was assumed. This in turn leads to a stabilization of the evolution of creep strains (see yc ) since they are, as the plastic evolutions, mainly related to the deviatoric stresses (see equation (56)). The coupling between the chemical problem (creep) and the mechanical one (stress state and plastic evolutions) is clearly illustrated by this simple example. Viscous hardening plasticity for concrete in high rate dynamics Formulation. The last application is concerned with the modelling of concrete behaviour in high rate dynamics. It was shown experimentally that the observed increase of strength with an increasing loading rate (up to 100 =s) is related to the presence of water in the nanopores of concrete [26] (pores of calcium silicate hydrates). Therefore, a viscous phenomenon [8] can be thought to be at the origin of rate eects for strain rates up to 100=s. At the macroscopic level, this microscopic viscous phenomenon can be accounted for by introducing a scalar viscous variable [16], denoted as x (internal variable). Induced rate eects under dynamic loading (increase of strength) are then described in the model by viscous hardening, which mainly leads to express the hardening force (modelling the increase of strength) as a function of the variable x (modelling the phenomenon at the origin of this increase of strength). This application is used to demonstrate that the framework of closed reactive porous media can be extended to the modelling of phenomena involving internal water movements within the porous space of materials. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS 93 Considering again the system of equations (6) describing the possible couplings between the state variables of the model, it takes the following form when applied to concrete modelling in high rate dynamics: s 0 |s| db dU − dUp a @f () x 1 c d d (58) 1 + a log fc () 0 = @ x∞ x x∞ dAv dx s 0 −K |s| The system of equations (58) results from the following assumptions: 1. The viscous phenomenon at the origin of rate eects mainly depends on the deviatoric stress components (@Av =@b = s=| s |) and on its own state (@Av =@x =−K with K¿0 as a constant parameter). Hence, for purely hydrostatic loading, there are no viscous rate eects [16]. The inuence of plastic hardening/softening is negligible (@Av =@ ≈ 0). 2. The visco-elastic coupling is non-linear leading to deviatoric viscous strains (dUv = (s=| s |) dx). 3. Viscous and plastic hardening=softening are not independent (@=@ = @fc ()=@[1 + a log (x=x∞ )] depends on viscous strain x and @=@x = afc ()=(xx∞ ) depends on the plastic hardening =softening variable ). For the plastic criterion, a Drucker–Prager yield surface can be chosen: (59) f(b; ) = 0 I1 + | s | + k0 ∗ (; x) = 0 p p wherein 0 = (1 − r) 2=3=(1 − 2r) and k0 = r 2=3=(1 − 2r) are two material parameters depending on the ratio r = fbc =fc with fbc and fc as the biaxial and uniaxial compressive strength, respectively. For plastic hardening, a hyperbolic relation was adopted [15]: ( − u )2 (60) fc () = ! + (1 + !) 1 − u2 where w = fy =fc is the ratio between the elastic limit in compression fy and the compressive strength fc . u is the value of the hardening variable when the compressive strength is reached. Associated plasticity (g = f) is considered again together with a non-associated hardening rule (h = f(b; ) = 0 I1 + | s | + (; x)). 4. The viscous hardening relation is linear in the log-scale (@=@x = afc ()=(xx∞ ) = f(1=x) with a as a material parameter and x∞ as the value of the viscous strain x when the tensile strength is reached under static loading) which reects the experimental dependency of the strength on the strain rate [25; 26]. Integration of the system of equations (58) yields the state equations of the model as follows: s U − Up = : b + Uv with dUv = dx |s| x (61) (; x) = fc () 1 + a log x∞ A(b; x) = | s | − Kx Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 94 J. SERCOMBE, F.-J. ULM AND H.-A. MANG Table III. Material parameters for viscous hardening plasticity E (GPa) 30 fc (MPa) fbc (MPa) N (GPa.s) K (GPa) w a u 0·2 30 34·5 1000 200 0·25 −0·08 0·0022 To complete the system (61), the kinetic law for the viscous process at the origin of rate eects needs to be chosen. In a rst approximation, a linear relationship was used [16]: ẋ = 1 (| s | − Kx) N (62) wherein N¿0 is the viscosity associated with the viscous phenomenon occurring in high-rate dynamics within the nanopores of concrete. With respect to the two previous applications, the system is fully coupled. The solution of the viscous problem (62) depends on the mechanical solution (stress state) and vice versa (61). Furthermore, the mechanical solution (stress state) depends on the plastic solution which cannot be obtained without solving for the coupled viscoplastic hardening=softening force (; x). Note that neglecting the eect of plastic hardening=softening on the viscous process (@Av =@ ≈ 0) leads to an unsymmetric system of equations (58) for the coupled problem which results in a loss of symmetry of the consistent tangent stiness matrix. Complete expressions for the latter and the system (35) are given in the Appendix IV. Viscoplastic analysis of the plate. The same boundary conditions (xed boundaries and prescribed displacements at the remaining boundaries, respectively) as in the previous example are used for this analysis. However, the prescribed displacements are chosen such (increment of u imp = v imp = − 0·0002) that a biaxial compressive stress state is achieved in the plate (see case (3) in Figure 2). Three calculations are performed with a decreasing time step of 0·1; 0·001 and 0·00001 s. This results in loading conditions varying from quasi-static to dynamic (˙x = ˙y is equal to 0·002; 0·2, and 20=s, respectively). The material is chosen as elastoplastic with viscoplastic hardening for the purpose of consideration of the sensitivity with respect to the rate of local straining. Its parameters are given in Table III. In Figure 5 the stress–strain curves and the evolution of the hardening force (; x) are shown for the three loading rates. In the left part of Figure 5, the (biaxial) stress–strain response is mainly governed by the combination of viscoplastic hardening=softening. This is reected by the shape of the stress–strain curves showing a decrease of non-linearity with increasing strain rate (consistent with experimental observations in uniaxial compression [25]). This results from the eect of viscous hardening (x) (increasing the initial plastic threshold and the biaxial yield stress) on the plastic hardening relation (; x) = fc ()(x) (see equation (61)). This coupling is better observed from the curves describing the evolution of the hardening force for the three calculations. These curves are obviously not parallel. They have dierent maximum values for the hardening force. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS 95 Figure 5. Stress (y )-strain (y ) curves (left) and evolution of the hardening force (; x) (right) for the three strain rates A structural application of viscous hardening plasticity can be found in Reference [27] (viscous hardening plasticity for ber-reinforced concrete applied to 3-D structural analysis of the impact resistance of containers to drop tests from a height of 5 m). CONCLUSIONS In this paper, the formulation of a return mapping algorithm for the integration of stresses and the updating of state variables for chemoplastic constitutive laws with internal couplings was presented. Details of the application of this algorithm to the modelling of the behaviour of concrete at early ages and under high loading rates were discussed with an emphasis on the numerical treatment of the problem. To show the generality of the present algorithm, it was also applied to the integration of stresses for a combined plasticity-creep constitutive law. Numerical simulations were then presented for each one of these applications with special emphasis on the couplings of the models. From a numerical standpoint, each one of these tests showed the usefulness of the algorithm for increasingly complex cases: perfectly plastic behavior with chemical hardening, coupled creep-plastic behaviour, plastic behaviour with combined visco-plastic hardening =softening. The following conclusions can be drawn from these analyses: (i) the numerical integration of the stresses for a chemoplastic model with internal couplings requires solution of a 2 × 2 system of algebrac equations at each time step; (ii) this algorithm is unconditionally stable, i.e. its stability is independent of the chosen time step; the algorithm preserves the quadratic rate of asymptotic convergence which is typical for elasto-plastic constitutive laws; (iii) it can be used to update the state variables for any type of plastic constitutive law combined with a chemophysical phenomenon (creep, chemical reaction, viscous eect, . . .), which in turn has an inuence on the mechanical properties of the material; (iv) the extension of the presented algorithm to multi-surface plasticity remains straightforward, provided the internal phenomenon is described by a single variable; (v) analysis with combined chemoplastic hardening =softening can be performed without further diculties. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 96 J. SERCOMBE, F.-J. ULM AND H.-A. MANG APPENDIX I Expressions of (11 ; 12 ; 21 ; 22 ) in Box 1: (k) @g ( (k) )T h (k) i @b n+1 @f @f (k) (k) 11 = An+1 ; @b n+1 @ n+1 @h @ n+1 (k) bn+1 (k) (k) )T h i @f @f (k) (k) A ; 12 = @ n+1 @b n+1 @ n+1 −D−1 @ n+1 (k) @g ( )T (k) (k) h i @b n+1 @F @F (k) −1 ; − D A 21 = (k) n+1 n+1 @b n+1 @ n+1 @h @ n+1 ( ( 22 = (k) bn+1 )T (k) (k) (k) h i (k) @F @ (k) @F @F @F (k) −1 ; − D − (63) A − n+1 −1 @ @b n+1 @ n+1 n+1 @ @ n+1 @ n+1 −Dn+1 @ n+1 Expressions of (1 ; 2 ) in Box 1: (k) 1 = fn+1 (k) 2 = Fn+1 ( (k) (k) )T h i @f @f (k) ; A(k) n+1 Rn+1 @b n+1 @ n+1 ( (k) (k) )T h i @F @F (k) −1 A(k) ; −Dn+1 n+1 Rn+1 @b n+1 @ n+1 − − (64) APPENDIX II Since for the chemoplastic model for concrete at early ages, the chemomechanical problem is uncoupled, the chemical solution can be found prior to the mechanical one. Considering the discretized version of kinetic law (53), the following condition has to be satised: Ea =0 (65) Fn+1 = − n+1 + t A0 (1 − n+1 ) exp − RT The solution of this linear equation leads the chemical variable n+1 (or the increment n+1 ). Replacing I1; n+1 in the tension cut-o (51) by its expression in function of I1;trial n+1 , n+1 and Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS 97 n+1 , leads to the following equation: s ∞ hn+1 − 0 i n+1 −ft∞ =0 − 9K( ) − 9K( ) fn+1 = I1;trial n+1 n+1 n+1 n+1 (1 − s ) (1 − 0 ) | {z } (66) I1; n+1 Since n+1 is known, linearization of (66) with respect to leads the plastic solution in the form n+1 = 1 fn+1 = 9K(n+1 ) 11 (67) s Note that (66) is only valid if ¿s . The corresponding term (9K(n+1 )∞ =(1 − s )n+1 ) does not appear in case this condition is not fullled (6s ). K(n+1 ), the bulk modulus, is given by K(n+1 ) = E∞ hn+1 − 0 i E(n+1 ) = 1 − 2 1 − 2 (1 − 0 ) (68) with the Poisson’s ratio. The expression of the consistent tangent stiness matrix for the chemoplastic model for concrete at early ages is db = C (n+1 ) − K(n+1 )1 ⊗ 1 (69) dU n+1 APPENDIX III Replacing |sn+1 | in the discretized version of the von Mises criterion (57) and of the kinetic law trial |, n+1 and n+1 leads to the system: (56) by its expression as function of |sn+1 r 3 trial fn+1 = | s | − 2Gn+1 − 2Gn+1 − ft | {z } 2 | sn+1 | (70) t trial (| sn+1 | − 2Gn+1 − 2Gn+1 − n+1 ) Fn+1 = − n+1 + Linearization of (70) with respect to and leads to the expressions of (11 ; 12 ; 21 ; 22 ) for the combined plasticity-creep model: 11 = 2G; 12 = 2G t t (2G + ) 21 = 2G ; 22 = 1 + (71) 1 and 2 are respectively equal to fn+1 and Fn+1 in (70). System (35) is then solved to obtain the plastic and chemical solutions (n+1 , n+1 ). G, the shear modulus is given by G= Copyright ? 2000 John Wiley & Sons, Ltd. E 2(1 + ) (72) Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 98 J. SERCOMBE, F.-J. ULM AND H.-A. MANG The expression of the consistent tangent stiness matrix for the combined plasticity-creep model is 4G 2 db trial trial trial trial = C − 2Gnn+1 ⊗ nn+1 − trial (n+1 + n+1 )(Idev − nn+1 ⊗ nn+1 ) (73) dU n+1 | sn+1 | trial trial trial with nn+1 = sn+1 =| sn+1 |, and Idev = I − 1=31 ⊗ 1 the deviatoric operator with I standing for the fourth-order unit tensor. Note that = 2G(1 + t=)¿0, the determinant of system (71), is always strictly positive. APPENDIX IV Replacing I1; n+1 and | sn+1 | in the discretized version of the Drucker–Prager criterion (59) and of trial the kinetic law (62) by their expressions as function of I1;trial n+1 , |sn+1 |, n+1 and xn+1 leads to the system xn+1 trial trial fn+1 = 0 (I1; n+1 − 9K0 n+1 ) + | sn+1 | − 2G(n+1 + xn+1 ) +k0 fc 1 + a log x∞ {z } {z } | | I1; n+1 | sn+1 | (74) t trial (| s | − 2Gn+1 − 2Gxn+1 − Kn+1 ) Fn+1 = −xn+1 + N n+1 Linearization of (74) with respect to and leads to the expressions of (11 ; 12 ; 21 ; 22 ) for viscous hardening plasticity: xn+1 @fc 1 + a log 11 = 9K02 + 2G − k0 @ n+1 x∞ a 12 = 2G − k0 fc xn+1 x∞ (75) t 21 = 2G N t (2G + K) 22 = 1 + N 1 and 2 are respectively equal to fn+1 and Fn+1 in (74). System (35) is then solved at each (k) (k) |¡TOL and |Fn+1 |¡TOL. Expressions iteration step (k) until convergence is achieved, i.e. |fn+1 (74) and (75) are only valid if xn+1 ¿0 (see [16]). @fc =@n+1 is given by (1 + !) @fc @fc = −2 (n+1 − u ) if (n+1 ¡u ) and = 0 if (n+1 ¿u ) (76) @ n+1 u2 @ n+1 The expression of the consistent tangent stiness matrix for viscous hardening plasticity is: 1 db trial trial = C − (3K0 1 + 2Gnn+1 ) ⊗ [3K0 22 1 + 2Gnn+1 (22 − 12 t=)] dU n+1 1 trial trial ⊗ [2Gnn+1 (11 t= − 21 ) − 3K0 21 1] − 2Gnn+1 − 4G 2 trial trial ⊗ nn+1 ) (n+1 + xn+1 )(Idev − nn+1 trial | sn+1 | Copyright ? 2000 John Wiley & Sons, Ltd. (77) Int. J. Numer. Meth. Engng. 47, 75–100 (2000) CHEMOPLASTIC LAWS WITH INTERNAL COUPLINGS 99 with = 11 22 − 12 21 , the determinant of system (75), which is always strictly positive since k0 ¡0, @fc =@n+1 ¿0 (plastic hardening), a¡0 (viscous hardening), and xn+1 ¡x∞ (see [16]). ACKNOWLEDGEMENTS The rst and the third writer would like to acknowledge the support of this research by the Austrian Foundation for the Advancement of Scientic Research within the framework of the Joint Research Initiative ‘Numerical Simulation in Tunnelling’. REFERENCES 1. Chen WF. Plasticity in Reinforced Concrete. McGraw-Hill: New York, 1982. 2. Mestat P. Constitutive laws for geomaterials and nite element modelling. Technical Report, Laboratoires des Ponts et Chaussees, Paris, 1993 (in French). 3. Mould JC, Levine H, Tennant D. Evaluation of a rate-dependent three invariant softening model for concrete. In Fracture and Damage of Quasi-Brittle Structures, International Conference, Bazant ZP, Bittnar Z, Jirasek M, Mazars J (eds). Madrid, Spain, E & FN SPON, London, UK, 1994; 231–242. 4. Bicanic N, Zienkiewicz OC. Constitutive model for concrete under dynamic loading. Earthquake Engineering and Structural Dynamics 1983; 11:689–710. 5. Bazant ZP. (ed.). Mathematical Modelling of Creep and Shrinkage in Concrete. Wiley: Chichester, England, 1988. 6. Larive C. Apport combines de l’experimentation et de la modelisation a la comprehension de l’alcali-reaction et de ses eets mecaniques. Ph.D. Thesis, Laboratoire Central des Ponts et Chaussees, Paris, France, 1997 (in French). 7. Torrenti JM, Mainguy M, Adenot F, Tognazzi C. Modelling of leaching in concrete. In Computational Modelling of Concrete Structures, de Borst R, Bicanic N, Mang H, Meschke G (eds). vol. 1, Badgastein, Austria. Balkema: Rotterdam, The Netherlands, 1998; 531– 538. 8. Rossi P. A physical phenomenon which can explain the mechanical behaviour of concrete under high strain rates. Materials and Structures 1991; 24:422– 424. 9. Ulm F-J, Coussy O, Bazant ZP. The ‘chunnel’ re: thermal spalling due to chemoplastic softening in rapidly heated concrete. I: constitutive modelling. Journal of Engineering Mechanics, ASCE 1999; 125(3):272–282. 10. Coussy O. Mechanics of Porous Continua. Wiley: Chichester, West Sussex, England, 1995. 11. Coussy O, Ulm F-J. Creep and plasticity due to chemo-mechanical couplings. Archive of Applied Mechanics 1996; 66:523–535. 12. Ulm F-J, Coussy O. Modeling of thermochemomechanical couplings of concrete at early ages. Journal of Engineering Mechanics, ASCE 1995; 121(7):785 –794. 13. Simo JC, Hughes TJR. Computational Inelasticity. Springer: Berlin, 1998. 14. Ulm F-J, Coussy O. Strength growth as chemo-plastic hardening in early age concrete. Journal of Engineering Mechanics, ASCE 1996; 122(12):1123–1132. 15. Hellmich Ch, Ulm F-J, Mang HA. Multisurface chemoplasticity for shotcrete. Part I: Material modeling and numerical aspects. Journal of Engineering Mechanics; ASCE 1999; 125(6):692–701. 16. Sercombe J, Ulm F-J, Toutlemonde F. Viscous hardening plasticity for concrete in high-rate dynamics. ASCE Journal of Engineering Mechanics 1998; 124(9):1050 –1057. 17. Lyons LPR, Li X, Duxbury PG, Grook AJL. Integration of isotropic =anisotropic combined plasticity and creep. In Third International Conference on Computational Plasticity, Fundamental and Applications, Owen DRJ, Oñate E, Hinton E (eds). Barcelona, Spain, Pineridge Press: Swansea, UK, 1992; 249–262. 18. Ortiz M, Simo JC. An analysis of a new class of integration algorithms for elastoplastic constitutive relations. International Journal for Numerical Methods in Engineering 1986; 23:353 –366. 19. Hellmich Ch, Ulm F-J, Mang HA. Multisurface chemoplasticity for shotcrete. Part II: 2D structural analysis of tunnels driven according to the NATM. Journal of Engineering Mechanics; ASCE 1999; 125(6):702–714. 20. Ulm F-J, Coussy O. Couplings in early-age concrete: from material modelling to structural design. International Journal for Solids and Structures 1998; 35(31–32):4295– 4312. 21. Ulm F-J, Coussy O, Hellmich Ch. Chemoplasticity: a review of evidence. In Proceedings of the International Conference on Computational Modelling of Concrete Structures, de Borst R, Bicanic N, Mang H, Meschke G (eds). Badgastein, Austria. Balkema: Rotterdam, The Netherlands, 1998; 421– 440. 22. Hellmich Ch, Ulm F-J, Mang HA. Elucidation of the soil-shotcrete interaction problem governing the NATM, based on visco- and chemoplastic material laws. In Proceedings of the International Conference on Computational Modelling of Concrete Structures, de Borst R, Bicanic N, Mang H, Meschke G (eds). Badgastein; Austria, vol 1. Balkema: Rotterdam, The Netherlands, 1998; 855–866. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000) 100 J. SERCOMBE, F.-J. ULM AND H.-A. MANG 23. Kojic M, Bathe K-J. The ‘eective-stress-function’ algorithm for thermo-elasto-plasticity and creep. International Journal for Numerical Methods in Engineering 1987; 24:1509–1532. 24. Li X, Duxbury PG, Lyons P. Coupled creep-elastoplastic-damage analysis for isotropic and anisotropic nonlinear materials. International Journal of Solids and Structures 1994; 31(9):1181–1206. 25. Bischo PH. Perry SH. Compressive behaviour of concrete at high strain rates. Materials and Structures 1991; 24:425–450. 26. Toutlemonde F. Shock strength of concrete structures. From material behavior to structural design. Ph.D. Thesis, Laboratoire Central des Ponts et Chaussees, Paris, France, July 1995 (in French). 27. Sercombe J, Toutlemonde F, Ulm F-J. From material behaviour to structural design for concrete in high rate dynamics. In Computational Modelling of Concrete Structures, de Borst R, Bicanic N, Mang H, Meschke G (eds), vol. 2, Badgastein, Austria, Balkema: Rotterdam, The Netherlands, 1998; 633– 642. Copyright ? 2000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng. 47, 75–100 (2000)

1/--страниц