Tunnelling and Underground Space Technology 33 (2013) 171–185 Contents lists available at SciVerse ScienceDirect Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust Modeling time-dependent behavior of gas caverns in rock salt considering creep, dilatancy and failure S. Nazary Moghadam a, H. Mirzabozorg a,⇑, A. Noorzad b a b Civil Engineering Department, K.N. Toosi University of Technology, No. 1346, Valiasr Street, Mirdamad Intersection 19967, P.O. Box. 15875-4416, Tehran, Iran Faculty of Water and Environmental Engineering, Power and Water University of Technology, Tehran, Iran a r t i c l e i n f o Article history: Received 20 January 2012 Received in revised form 12 July 2012 Accepted 1 October 2012 Available online 16 November 2012 Keywords: Salt cavern Creep Non-linear analysis Finite element method Large deformation Viscoplasticity a b s t r a c t In this paper, a numerical model of the time-dependent behavior of underground caverns excavated in rock salt is presented. An elasto-viscoplastic constitutive model is utilized to describe dilatancy, shortterm failure as well as long-term failure during transient and steady state creep of rock salt. The constitutive model is employed by a Lagrangian finite element formulation to simulate the stress variation and ground movement during creep of rock salt around the cavern within the framework of large deformations. Finally, finite element analyses are provided in order to investigate the efficiency and applicability of the proposed computational model to represent time-dependent response of gas storage caverns in rock salt. Ó 2012 Elsevier Ltd. All rights reserved. 1. Introduction Rock salt has increasingly attracted the attention of petroleum industry as an ideal material providing valuable storage capacity for natural gas and other hydrocarbons. Due to very low permeability and high ductility of rock salt, there is also motivation from nuclear industry with particular interest related to the disposal of radioactive waste. Underground storages of petroleum products are more economical than surface tank storages which enable much higher volumes and pressures to be achieved. This becomes more evident taking into account the small amount of land required and security against external influences. The other advantage of underground storages in salt deposits is that salt caverns may be created by solution mining techniques instead of the more costly conventional excavation techniques. The above-mentioned advantages make salt deposits one of the prime alternatives for the construction of underground cavities for storage of petroleum products or radioactive wastes and hazardous chemical wastes. Time-dependent ground movement around an underground opening excavated in a rock salt formation can progress continuously leading to significant or even full closure of the opening, and also stress relaxation can occur around the opening up to small ultimate stresses. The consequence of this behavior is that time- ⇑ Corresponding author. Tel.: +98 21 88779623; fax: +98 21 88779475. E-mail addresses: saeed_nazarimoghadam@yahoo.com (S. Nazary Moghadam), mirzabozorg@kntu.ac.ir (H. Mirzabozorg), noorzad@pwut.ac.ir (A. Noorzad). 0886-7798/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.tust.2012.10.001 dependent inelastic behavior of rock salt has to be considered in the design and the stability analysis of underground excavations in rock salt. A crucial time-dependent property of rock salt is the time-dependent deformation behavior or creep behavior defined as the inelastic deformation in time without fracturing. The creep deformation can be decomposed into two components, namely transient and steady state creep deformation. Transient creep is the creep deformation component which continuously decreases in time until a stable state is reached, while steady state creep deformation continues indefinitely under constant or nearly constant stress state. In general, rock salt exhibits both transient and steady state creep and it is only the time interval elapsed from the application of the load that determines which of the two components of the creep deformation is dominant. In classic models, the creep deformation is generally described using viscoplastic theory based on the existing models for metals neglecting the inelastic volume changes caused by the closing or the opening and growth of microcracks. However, experimental evidence reveals that a typical manifestation of creep deformation of rock salt involves irreversible dilatancy and compressibility in a broad range of stress states. Hence, for a comprehensive description of the time-dependent behavior of rock salt it is important to choose a viscoplastic constitutive equation which considers irreversible compressibility–dilatancy, damage and failure as well. In order to investigate the time-dependent behavior of underground caverns excavated in rock salt, the constitutive modeling of rock salt is of vital importance. Constitutive models used to de- 172 S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 scribe the mechanical behavior of rock salt are usually based on viscoplasticity. Apparently, Perzyna (1966) was the first to suggest the theory of viscoplasticity describing the time-dependent inelastic behavior, in which combined treatment of rheological and plastic phenomena is essential, using the overstress concept. Viscoplastic problems can be solved using the numerical procedure formulated by Zienkiewicz and Cormeau (1972, 1974) based on the initial strain technique in the context of finite element discretization. In this time-marching procedure, after an elastic prediction of the stress change is made, if resulting stress state exceeds the yield criterion, viscoplastic strains are assumed to develop in time until the yield condition is satisfied within an acceptable tolerance. Several viscoplastic constitutive models have been developed and tested for the description of the time-dependent behavior of rock salt. In this regard, an empirical creep model of rock salt was developed for the Waste Isolation Pilot Plant known as the WIPP creep model in which the inelastic volume changes are neglected (Herrmann et al., 1980a,b). Cristescu and Hunsche (1992, 1993) proposed a viscoplastic constitutive model for rock salt based on triaxial experimental data in which transient creep, irreversible volumetric compressibility–dilatancy, damage and failure were considered. In that research, the effect of steady state creep and large deformation was neglected. In another research, in the same line of thought, Cristescu (1993) presented a procedure to determine a general viscoplastic constitutive equation which besides transient creep, irreversible compressibility–dilatancy, short term failure and damage, describes steady state creep. He introduced long-term creep failure utilizing an energetic criterion, and also a constitutive model for rock salt was given. Using the same procedure outlined in the previous researches conducted by Cristescu, Jin and Cristescu (1998) formulated a new viscoplastic constitutive model for transient creep of rock salt which was more appropriate than previous models to incorporate into a finite element code. The proposed constitutive model was employed by the authors to analyze the stress distribution around a plane strain borehole excavated in rock salt neglecting the effect of large deformation and steady state creep. Paraschiv-Munteanu and Cristescu (2001) derived a semi-analytical procedure based on finite element discretization for the analysis of stress variation during creep of rock salt around a plane strain borehole which describes compressibility– dilatancy during transient and steady state creep as well as damage and failure of the rock salt in the case of small deformations. Moreover, within the same framework, studies of Dawson and Munson (1983), Munson et al. (1990), Aubertin et al. (1991), Yahya et al. (2000), Peric and Crook (2004), Heusermann et al. (2003) and other researches (Chang and Zoback, 2010; Fuenkajorn and Phueakphum, 2010; Glamheden and Curtis, 2006; Li et al., 2005; Shao et al., 2003; Sharifzadeh et al., 2012; Wang et al., 2011; Xia et al., 2007; Zhu et al., 2008) are worth mentioning. In particular, Yahya et al. (2000) developed a unified viscoplastic model with a single set of equations and material constants for describing both short term and long term ductile behavior of rock salt. Furthermore, Peric and Crook (2004) carried out a research concerned with the aspects of computational strategies based on Lagrangian method for the simulation of large deformation problems encountered in the field of salt tectonics. Also, using a viscoplastic constitutive model for rock salt based on experimental data given by Yahya et al. (2000), a simulation of salt dome formation was presented by the authors. The objective of the present paper is to propose a computational algorithm to simulate the time-dependent response of underground openings excavated in rock salt. An elasto-viscoplastic constitutive model, which is comprised of both Peric’s viscoplastic power law representing transient creep and Cristescu’s steady state creep model, is employed by a Lagrangian finite element formulation to describe dilatancy, short-term failure as well as long- term failure during transient and steady state creep of rock salt in the framework of large deformations. The computational model is then implemented into a nonlinear finite element program, and afterwards the capability of the proposed computational algorithm to analyze the stress variation and ground movement during creep of rock salt around a storage cavern is investigated using an axisymmetric finite element representation. Finally, finite element analyses results obtained by the proposed model are presented and compared with those of WIPP steady state creep model. 2. Kinematics In a finite element representation, displacements can be expressed by the relationship u¼ 8 9 > > <u = > : v >¼ w ; X 8 9 > < pxi > = Ni pyi ¼ Np > : > ; pzi ð1Þ where the matrix N contains interpolation functions Ni and p is the nodal displacement vector. Based on Green’s strain, which is valid for large deflections, the nonlinear strain–displacement relation can be defined by the following tensor equation: Eij ¼   1 @ui @uj @uk @uk þ þ 2 @xj @xi @xi @xj ð2Þ where Eij is Green’s strain tensor. Also, Green’s strain vector E can be defined in terms of the linear strain vector El and the non-linear strain vector Enl as follows (Crisfield, 1991):   1 E ¼ fEx ; Ey ; Ez ; Exy ; Eyz ; Exz gT ¼ El þ Enl ¼ H þ AðhÞ h 2   @u @ v @w @u @ v @w @u @ v @w T ; ; ; ; ; ; ; ; h ¼ @x @x @x @y @y @y @z @z @z ð3Þ ð4Þ where h denotes the displacement gradient vector and A(h) is a suitably defined matrix operator which contains displacement derivatives. The matrixes A(h) and H are defined in Appendix A. The changes in the displacement gradient vector dh can be determined in terms of the interpolation functions Ni and nodal variables dp via dh ¼ Gdp ð5Þ where G contains the interpolation function derivatives. From the expression for A(h), it is clear that dA(h)h = A(h)dh, therefore, by taking the variation from Eq. (3), the incremental form of the Green’s strain can be expressed as dE ¼ Bnl ðpÞdp Bnl ðpÞ ¼ ½H þ AðhÞG ð6Þ ð7Þ where Bnl(p) is the nonlinear strain–displacement transformation matrix relating the increments of strain and displacement. 3. Elasto-viscoplastic constitutive relations 3.1. Elasto-viscoplastic model A comprehensive description of the time-dependent behavior of rock salt necessitates a viscoplastic constitutive model considering not only transient and steady state creep but also compressibility– dilatancy along with damage and failure. To this end, the viscoplastic constitutive equations proposed by Cristescu (1993), Jin and Cristescu (1998) and Paraschiv-Munteanu and Cristescu (2001) has been utilized to represent the mechanical behavior of rock salt. However, as the time intervals considered in this paper are in a range such that steady state creep is dominant, Cristescu’s viscoplastic transient model for rock salt may be cumbersome. There- S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 173 the basis of this argument, Eq. (13) can be rewritten in the rate form as follows: _ St E_ ¼ E_ e þ E_ v p ¼ E_ e þ E_ Tr v p þ Ev p ð14Þ where over-dot indicates the derivatives with respect to time and _ St E_ Tr v p and Ev p denotes transient and steady state creep strain rate, respectively. One of the main concepts incorporated in Cristescu’s constitutive model is the short-term failure surface which is the locus of maximum values of octahedral shear stresses determined with standard triaxial tests in relative short time intervals. The short-term failure surface bounds the domain where constitutive equations are valid. The following equation has been used to represent the short-term failure surface. sf ¼ c1  c2 exp c3 Fig. 1. Illustration of Cristescu’s steady state viscoplastic potential surface, shortterm failure surface and C–D boundary. fore, Peric’s viscoplastic power law (2004) has been employed to describe transient creep of rock salt. In this section it is assumed that compressive stresses and strains are positive while tensile stresses and strains are negative. The second Piola Kirchhof stress and Green’s strain components do not change when the material is subjected to rigid body motions. Therefore, any material description developed for infinitesimal displacement analysis using engineering stress and strain measures can directly be employed in large displacement and large rotation but small strain analysis, provided that second Piola Kirchhof stress and Green’s strain are used (Bathe, 1996). Therefore, in the following, constitutive equations are written in terms of second Piola Kirchhof stress and Green’s strain. Second Piola Kirchhof stress which energetically conjugates to the rate of Green’s strain is denoted by the vector S T ¼ fSx ; Sy ; Sz ; Sxy ; Syz ; Sxz g ð8Þ It should also be mentioned that, the constitutive functions used in this paper are dependent on stress invariants, namely mean stress 0 Sm and the second deviatoric stress invariant j2 defined explicitly as Sm ¼ Sx þ Sy þ Sz 3 n o S 0T ¼ S0x ; S0y ; S0z ; Sxy ; Syz ; Sxz ¼ fSx  Sm ; Sy  Sm ; Sz  Sm ; Sxy ; Syz ; Sxz g 0 j2  1  02 02 ¼ þ S2xy þ S2yz þ S2xz S þ S02 y þ Sz 2 x ð9Þ ð10Þ ð11Þ where S0 is the deviatoric stress vector. In addition, the octahedral shear stress s and the effective stress S are defined as rffiffiffiqffiffiffiffi pffiffiffi 2 0 2 j s¼ S¼ 3 2 3 ð12Þ Sm S ð15Þ In this equation, sf is the octahedral shear stress at failure, c1, c2 and c3 are material constants and S⁄ is the unit stress. Another crucial concept is the compressibility–dilatancy (C–D) boundary, i.e. X = 0, which separates the domain where a rock is compressible, i.e. X > 0, and the domain where the rock is dilatant, i.e. X < 0. In compressibility domain the inelastic volumetric strain rate is positive leading to an irreversible reduction of the absolute value of the volume as well as closure of microcracks and pores. In contrast, in dilatancy domain the inelastic volumetric strain rate is negative leading to an irreversible increase of the absolute value of the volume as well as opening of microcracks and pores. The C–D boundary has been defined by the following equation: X¼ s S þ f1 Sm S 2 þ f2 Sm S ¼0 ð16Þ where f1 and f2 are material parameters. In order to describe transient creep of rock salt the following viscoplastic flow rule is utilized: E_ Tr vp ¼ dETr @F vp ¼ K Tr hUðHÞi @S dt ð17Þ in which KTr is viscosity coefficient for transient creep, F is the viscoplastic potential function and the notation hi implies hUi ¼  U IF : U > 0 0 IF : U 6 0 ð18Þ The term U(H) is defined using Peric’s viscoplastic power law as follows: UðHÞ ¼ H H0 N 1 ð19Þ where H is the yield function, H0 is the yield stress and N is a material constant. As can be seen, the onset and continuance of transient creep deformation is governed by the stabilization boundary defined by the following equation: H  H0 ¼ 0 ð20Þ It is assumed that the total strain E can be decomposed into an elastic part Ee and an irreversible viscoplastic part Evp, and therefore the total strain can be expressed as The viscoplastic potential function for transient creep is defined assuming associated flow rule, i.e. F = H, and also the yield function is defined using the von Mises yield criterion as E ¼ E e þ Ev p H¼S ð13Þ In order to describe both transient and steady state creep, due to distinct volumetric behavior for the two kinds of creep, it is further assumed that the viscoplastic strain rate can be separated into transient creep strain rate and steady state creep strain rate. Indeed, unlike transient creep, the steady state creep must describe either dilatancy or incompressibility, otherwise the volume could shrink to zero by steady state creep (Cristescu and Hunsche, 1998). On ð21Þ It should also be mentioned that, the yield stress H0 can be a function of hardening parameter H0 defined by introducing a variable yield stress as a function of the equivalent viscoplastic strain Ev p as follows: H0 ¼ dH0 dEv p ð22Þ 174 S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185  1 4f 1 s 2 2 Sa f2 þ f2 þ S ¼ S 2f1  1 4f 1 s 2 2 Sb f2  f2 þ S ¼ S 2f1 ð27Þ ð28Þ Furthermore, the function Q is defined by the equation l @Q p s ¼ @ s S S ð29Þ where p and l are material parameters. As can be seen, in the compressibility region surfaces G = c represent horizontal straight lines in Sm  s plane indicating that in this region no steady state compressibility is possible, while in dilatancy region G = c depends on Sm which indicates that in dilatancy region steady state dilatancy is possible. It should also be mentioned that, classical failure criteria defined only based on stress invariants can account for short-term failure while failure of rock salt is strongly time-dependent. Therefore, an energetic long-term failure criterion based on inelastic strain work per unit volume which represents the progression in time of the damage has been used as follows: WV v p ðtÞ ¼ Z 0 t Sm E_ V v p dt _ _ _ E_ V v p ¼ Ev px þ Ev py þ Ev pz Fig. 2. Initial geometry and boundary conditions of the cavern. ð30Þ ð31Þ WV vp The equivalent viscoplastic strain is accumulated from the incremental equivalent viscoplastic strain defined by the following equation: dEv p E_ v p ¼ dt rffiffiffi 12 2 _2 1_2 Ev px þ E_ 2v py þ E_ 2v pz þ ¼ Ev pxy þ E_ 2v pyz þ E_ 2v pxz 3 2 in which denotes volumetric viscoplastic strain work per unit volume and E_ V v p termed volumetric viscoplastic strain rate. In compressibility domain W V v p increases which means that mechanical ð23Þ The equivalent viscoplastic strain can be considered as an internal variable which governs the history dependency of the material. It is worthy to note that the constitutive equation used for transient creep neglects any irreversible volumetric changes. Moreover, for the steady state creep of rock salt the following expression is employed: ESt vp ¼ dESt @G vp ¼ K St dt @S ð24Þ where KSt is a viscosity coefficient and G is the viscoplastic potential function for steady state creep defined by the equation 8  m  nþ1   > bS Sm s > þ Q Ss IF : ðs P sm Þ or ðs < sm and Sm 6 Sa Þ > nþ1 S S > > > > > >  m  nþ1   < bS Sa s þ Q Ss IF : ðs < sm and Sa < Sm < Sb Þ K St G ¼ nþ1 S S > > > > >    >         nþ1 nþ1 > > bS s m Sa nþ1 > þ Q Ss IF : ðs < sm and Sm P Sb Þ  SSb þ SSm : nþ1 S S ð25Þ in which b, m and n are material parameters and sm is the maximum value of octahedral shear stress of C–D boundary that can be taken as sm S ¼ f22 4f1 ð26Þ and also parameters Sa and Sb are obtained by intersecting the C–D boundary with a horizontal straight line in Sm  s plane for a fixed value of s < sm, according to Fig. 1, Fig. 3. Finite element configuration of the cavern analyzed by Owen and Hinton. 175 S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 Table 1 Material parameters for transient creep of the rock. Parameter H0 (MPa) N KTr (1/year) H (MPa) Value 10 1 0.075 0 energy is stored by closing microcracks, while in dilatancy domain WV v p decreases which shows that mechanical energy is released by opening and forming micro cracks. It can be concluded that, as stress state locates in the compressibility domain, damage decreases until the maximum value of W V v p is reached at the crossing of the C–D boundary, and afterwards damage begins to increase when dilatancy takes place. Consequently, the total decrease of Fig. 5. Finite element mesh of the cavern excavated at 800 m below the ground surface. Table 2 Material parameters for WIPP steady state creep model of the rock salt. Parameter D (MPaM/day) M Q (Cal/Mole) R (Cal/Mole °K) T (°K) Value 0.1257 4.9 12,000 1.987 304   V WV v p starting from its maximum value W v p max can be a measure of the damage of rock salt given by   dðtÞ ¼ W V vp max  WV v p ðtÞ ð32Þ where d(t) is the damage parameter at time t which is a moment of time during the period when dilatancy occurs. In this paper, as the proposed constitutive model neglects compressibility, the parame  ter W V vanishes. Damage parameter at failure df which is the vp max maximum energy release due to microcracking during dilatancy can be obtained as     Sm V df ¼ W V v p max  W v p failure ¼ d0 S  2 3 ð33Þ where d0 is a material parameter. Another steady state creep flow rule utilized in the present paper is the WIPP empirical creep law (Herrmann et al., 1980a,b) which is defined by the equation below: Fig. 4. Comparison of the results obtained by the finite element code with those obtained by Owen and Hinton; (a) variation of the radial displacement, (b) variation of the radial stress, and (c) variation of the tangential stress.   3 Q S0i;j DSM exp  E_ St v pi;j ¼ RT 2S ð34Þ 176 S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 Table 3 Material parameters for Cristescu’s steady state creep model of the rock salt. Parameter Value b (1/day) m n l p (MPa/day) f1 f2 S⁄ (MPa) S0 (MPa) c1 (MPa) c2 (MPa) 0.864  109 5 0.1 5 0.3677184  107 0.017 0.9 1 53 38 34.9 0.04 0.0727 c3 d0 (MPa) dEiv p ¼ E_ iv p dt i 6 T Ei rffiffiffi 12 2 2 1 2 Ex þ E2y þ E2z þ E¼ Exy þ E2yz þ E2xz 3 2 ð41Þ ð42Þ In addition, the change of time step length between successive intervals can be limited by introducing the factor K as follows: dt iþ1 6 Kdti ð43Þ As a result, the minimum value of time step lengths computed at integration points which satisfy restrictions represented by Eqs. (41) and (43) will be chosen for time integration process. For the explicit time marching scheme, stable and reasonably accurate solu- where D and M are material constants, R is the universal gas constant, Q is the activation energy and T is the temperature in Kelvin (K). The WIPP creep law has been used for comparison purpose. 3.2. Integration of constitutive equations For elasto-viscoplastic problems, the accuracy of the global solution is strongly affected by the accuracy of the computation of the stress increments. As the stress–strain relationship take an incremental form, a numerical integration procedure would be used at the integration point level. To this end, viscoplastic strain increment dEiv p occurring in a time interval dti can be determined using the time stepping scheme as h i dEiv p ¼ dt i ð1  hÞE_ iv p þ hE_ iþ1 vp ð35Þ in which for h = 0 the forward-Euler time integration scheme, which is also referred to as fully explicit method, is obtained while h = 1 gives the backward-Euler or fully implicit scheme. In this paper, the fully explicit method is employed, and therefore the change in viscoplastic strain is defined as dEiv p ¼ E_ iv p dti ð36Þ Since it is assumed that the total strain E can be decomposed into an elastic part Ee and a viscoplastic part Evp, and since the elastic strain Ee is related to the stress S by the elastic modular matrix C, stresses can be evaluated at the end of the time increment dti as follows:   S iþ1 ¼ C Eiþ1  Eiþ1 vp X E_ iv p dt i Eiv p ¼ ð37Þ ð38Þ Using this procedure, stresses are always updated from the end of the last converged equilibrium state. In addition, in order to calculate the tangent stiffness matrix, the incremental form of Eq. (37) is needed. Differentiation of Eq. (37) gives   dS i ¼ C dEi  dEiv p   dS i ¼ C Binl ðpÞdpi  E_ ivp dti ð39Þ ð40Þ It should also be mentioned that, although explicit scheme can simplify the time integration process considerably, accuracy of the solution can deteriorate severely with increase of time step length which may lead to unstable results. Therefore, the time step length should be limited in order to achieve a stable and accurate solution. To this end, the factor T is defined which limits the maximum effective viscoplastic strain increment dEiv p as a fraction of the total effective strain E as follows: Fig. 6. Comparison of the results obtained by Cristescu’s steady state creep model with those obtained by WIPP steady state creep model; (a) variation of the radial displacement, (b) variation of the radial stress, and (c) variation of the tangential stress. S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 177 virtual work in the framework of Lagrangian description, a known reference configuration is needed. In this regard, the total Lagrangian approach uses the initial unstressed configuration of the body as its reference configuration. As a result, the equilibrium equation of a deforming continuum body at time t + Dt based on total Lagrangian approach can be written as V ¼ V int  V ext ¼ g T dpv ¼ 0 Z V int ¼ S T dEv dV0 ¼ qTint dpv V ext ¼ qText dpv g ¼ qint  qext ¼ 0 ð44Þ ð45Þ ð46Þ ð47Þ where Vint is the virtual work due to internal forces, Vext is the virtual work due to external forces, g is the out-of-balance force vector, dpv is an arbitrary small virtual displacement, qint and qext are the internal and external force vectors, respectively, and S is second Piola Kirchhof stress corresponding to the time t + Dt but is measured using the initial unstressed configuration with volume V0 . Upon substitution of Eq. (6) into Eq. (45) the internal force vector can be obtained by the equation qint ¼ Z Bnl ðpÞT SdV0 ð48Þ Consequently, the equilibrium equation can be rewritten as g¼ Z Bnl ðpÞT SdV0  qext ¼ 0 ð49Þ In order to solve the nonlinear equilibrium equations represented by Eq. (49), a time stepping technique is used. The time stepping procedure linearizes the equilibrium equation with respect to time Fig. 7. Comparison of the results obtained by Cristescu’s steady state creep model with those obtained by WIPP steady state creep model; (a) variation with time of the radial convergence at point A, (b) variation with time of the vertical convergence at point C, and (c) variation with time of the vertical convergence at point D. tion can be obtained using a value of T in the range 0:1 6 T 6 0:2 but near the failure lower values in the range 0:01 6 T 6 0:05 is necessary, and also a value of K ¼ 1:5 is suitable (Zienkiewicz and Cormeau, 1974; Owen and Hinton, 1980). 4. Equilibrium equations Based on the principle of virtual work, a continuum body is in equilibrium if the virtual work of all forces acting on the body under a virtual displacement is zero. In order to apply the principle of Fig. 8. SFR contour for Cristescu’s steady state creep model in the case of p = 0. 178 S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 K t ¼ K t1 þ K r Z K t1 ¼ Bnl ðpÞT CBnl ðpÞdV0 Z K r ¼ GT b SGdV0 ð55Þ ð56Þ ð57Þ where Kt1 contains the linear stiffness matrix and the initial displacement stiffness matrix, Kr is the initial stress stiffness matrix, and b S is a 9  9 matrix containing stress components defined in Appendix A. Fig. 9. LFR contour for Cristescu’s steady state creep model in the case of p = 0. Table 4 Material parameters for transient creep of the rock salt. Parameter H0 (MPa) N KTr (1/day) H (MPa) Value 9.15 4 0.152064 0 using a truncated Taylor series expansion. Therefore, during the time step dti the incremental form of equilibrium equations are given as g iþ1 ¼ g i þ dg i ¼ 0 dg i ¼ Z h ii d Bnl ðpÞT S i dV0 þ ð50Þ Z h ii Bnl ðpÞT dS i dV0  dqiext h ii S i Gdpi d Bnl ðpÞT S i ¼ GT b ð51Þ ð52Þ where dSi is given by Eq. (40) and the time step counter i starts from a converged solution at the previous load increment. Using Eq. (40) and rearranging Eqs. (50)–(52), the displacement increment which occurs during time increment dti can be obtained by the following equations: K it dpi ¼ dF i Z dF ¼ Bnl ðpÞT C E_ v p dtdV0 þ dqext þ g ð53Þ ð54Þ in which dF is the incremental pseudo-load vector and Kt is the tangential stiffness matrix which is defined as Fig. 10. Results obtained by the present formulation in the case of p = 7 MPa; (a) variation of the radial displacement, (b) variation of the radial stress, and (c) variation of the tangential stress. S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 179 Fig. 11. Stress contours obtained by the present formulation in the case of p = 7 MPa; (a) radial stress contour, (b) vertical stress contour, (c) tangential stress contour, and (d) shear stress contour. 180 S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 sumed to be hydrostatic which means that the primary stresses are equal in all directions. Therefore, since at far enough distances from the excavation stresses must tend towards the primary values, the boundary conditions at the outer surface of the model are imposed such that the stresses at this surface become hydrostatic. This condition is simulated by applying distributed lateral edge loading, which is equal to the hydrostatic stress in the rock formation, at the outer surface of the model at time t = 0. Afterwards, due to the difference between the primary stress in the rock formation and the cavern internal pressure p, elastic stress redistribution occurs and thereafter creep deformation of the rock initiates. Moreover, the material properties for transient creep of the rock are presented in Table 1, and also elastic modulus, Poisson’s ratio and the density of the rock are E = 6.9 GPa, t = 0.4 and q = 2550 kg/m3, respectively. In Fig. 4a–c, the results obtained by the computational algorithm employed herein are compared with those achieved by Owen and Hinton, which is based on small displacement formulation, at stabilization of transient creep. Fig. 4a shows the variation of the radial displacement u with the radial coordinate x along the line AB, while Fig. 4b and c depict the distribution of the radial stress Sxx and the tangential stress Szz along the line AB, respectively. As can be seen from the figures, the results of Owen and Hinton are correctly reproduced by the finite element code. 5.2. Comparison between WIPP and Cristescu’s steady state creep model In the present analysis, using the proposed numerical procedure, Cristescu’s model is applied to simulate steady state creep of rock salt around a storage cavern and the results are compared Fig. 12. SFR contour obtained by the present formulation in the case of p = 7 MPa. 5. Numerical results Numerical simulation of the time-dependent behavior of storage caverns excavated in rock salt is provided in this section. The large finite element deformation formulation together with the elasto-viscoplastic constitutive equations derived in previous sections has been implemented in a nonlinear finite element program. Using this model, the efficiency and applicability of the proposed computational algorithm in modeling of the time-dependent behavior of salt caverns is investigated. In the following analyses, eight-nodded isoparametric elements with reduced, i.e. 2  2, Gaussian integration rule were employed. 5.1. Verification of the finite element code With the purpose of verifying the results from the finite element code, the problem of stresses induced in the vicinity of an excavated underground storage cavity solved by Owen and Hinton (1980) was recomputed. In this problem only the transient creep of a rock is considered, and also an axisymmetric representation is employed. The initial geometry, boundary conditions and the finite element mesh are shown in Figs. 2 and 3. As it can be seen from the figures, the bottom of the model is restrained against vertical displacements and the internal pressure of the cavern p is assumed to be zero. Also, as it is shown in the finite element mesh of the cavern, the center of the opening is located at the depth of 800 m below the ground surface. It is assumed that the cavern is instantaneously excavated at time t = 0, and also the gravitational body forces consistent with the density of the rock are applied to the model at time t = 0. For simplicity, in situ stress state before excavation of the cavity is as- Fig. 13. Displacement field obtained by the present formulation in the case of p = 7 MPa. S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 Fig. 14. Finite element mesh of the cavern excavated at 1400 m below the ground surface. with those of WIPP steady state creep model. The initial geometry, boundary conditions and initial conditions are the same as those represented in the first analysis while a finer finite element mesh shown in Fig. 5 is used. The displacements which are obtained in this condition are the summation of the displacements produced from geostatic condition and those elastically generated due to the excavation of the cavern together with the creep deformation of the rock salt. In reality, however, since the gravity forces exist before the excavation of the cavern, and since the relative displacements after excavation of the cavern are of concern, the geostatic displacements are removed from the total displacements resulted from the analyses. The material parameters for WIPP steady state creep model determined for New Mexico rock salt (Herrmann et al., 1980a,b) are presented in Table 2. In addition, the constitutive parameters for Cristescu’s steady state creep model (Cristescu, 1993; Jin and Cristescu, 1998; Paraschiv-Munteanu and Cristescu, 2001), determined based on the experimental data obtained for New Mexico rock salt (Herrmann et al., 1980a,b) as well as Avery Island rock salt (Hansen and Mellegard, 1980; Mellegard et al., 1981), are given in Table 3. According to the experimental data presented for Avery Island rock salt (Yahya et al., 2000) elastic modulus, Poisson’s ratio and the density of the rock salt are assumed to be E = 31 GPa, t = 0.38 and q = 2200 kg/m3, respectively. Fig. 6a–c show the variation of the radial displacement, radial stress and tangential stress with the radial coordinate along line AB at time t = 500 days, respectively. Also depicted in these figures are the instantaneous results obtained by elastic solution. In addition, the variation with time of the radial convergence u at point 181 A, the vertical convergence v at point C and the vertical convergence v at point D are plotted respectively in Fig. 7a–c in the both cases of Cristescu’s steady state creep model and WIPP steady state creep model. As shown, the absolute values of radial displacement obtained by Cristescu’s model are larger than those obtained by WIPP creep model near the cavern surface, although in farther distances from the cavern wall the difference between the two models diminishes gradually until the two models reach the elastic solution representing the primary conditions. In the neighborhood of the cavern surface, the absolute values of stress components obtained by Cristescu’s model as well as WIPP model are smaller than those obtained by elastic solution due to slow relaxation, while at greater distances from the cavity surface stress components approach asymptotically to the primary values representing the hydrostatic conditions. As observed, stress relaxation represented by Cristescu’s model is more considerable than the WIPP model which may bring about a sooner short-term failure for Cristescu’s model than WIPP model. It should also be mentioned that, due to stress relaxation around the cavern surface the rate of cavity surface convergence decreases steadily in the case of Cristescu’s model as well as WIPP model, and also Cristescu’s model shows faster cavity surface convergence compared with WIPP model. Finally, in order to investigate the stability of the cavern using Cristescu’s model, two parameters, namely short-term failure ratio (SFR) and long-term failure ratio (LFR), have been defined. Based on the Eq. (15), the parameter SFR is defined as the ratio of the current octahedral shear stress to the failure octahedral shear stress, i.e. SFR = s(t)/sf. Also, based on the Eqs. (32) and (33), the parameter LFR is defined as the ratio of the current damage parameter d(t), i.e. volumetric viscoplastic strain work per unit volume due to the dilatancy of rock salt, to the failure damage parameter, i.e. LFR = d(t)/df. The parameter SFR represents the allowable relative octahedral shear stress while the parameter LFR represents the allowable relative damage parameter. Therefore, the underground cavern is stable as long as the parameters SFR and LFR are smaller than one. For the SFR values greater than one, short-term failure of the rock salt occurs while for the LFR values greater than one, longterm failure of the rock salt is encountered. The SFR and LFR contours for the present simulation in the case of Cristescu’s model are illustrated in Figs. 8 and 9. According to the Fig. 9, as the LRF contour close to the cavern surface is not equal to zero, it can be concluded that the stress states in the neighborhood of the cavern surface are in the dilatancy domain. 5.3. Application to stability analysis of salt caverns In the last simulations, the capability of the proposed computational model to represent time-dependent response of gas caverns excavated in rock salt is studied. The initial geometry, boundary conditions, initial conditions and finite element mesh are the same as those represented in the previous section, except that here the opening is assumed to be filled with a gas under pressure, i.e. p – 0. The allowable range of cavern internal pressure mainly depends on the geostatic stress of the rock salt in the neighborhood of the cavern as well as the operational period of the cavern. In general, up to a specific value of p, a greater internal pressure contributes to the stability of the cavern. However, the maximum allowable internal pressure should be chosen such that the fracturing of rock salt and the increase in the permeability due to the microfracturing of rock salt can be prevented. It should also be mentioned that permeability of rock salt surrounding the cavern can increase considerably when the internal pressure is equal to the geostatic stress existing in the rock salt around the cavern (Rokahr et al., 1997). Therefore, in the present analyses, a cavern internal pressure which is appropriately lower than the geostatic stress 182 S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 Fig. 15. Comparison of the results obtained by the large displacements analysis with those obtained by small displacements analysis; (a) variation with time of the radial convergence at point A, (b) variation with time of the vertical convergence at point C, (c) variation with time of the vertical convergence at point D, and (d) variation with time of the cavern relative closure. of around 18 MPa is chosen. Also, as stress relaxation occurs due to the creep of rock salt, the minimum internal pressure of the cavern should be set such that the fracture strength of the rock salt cannot be exceeded in any location in the operational period of the cavern. As a result, in the present analyses, an internal pressure of 7 MPa is applied to the cavern surface at time t = 0 and then the cavern is permitted to creep for 30 years. The constitutive parameters for Cristescu’s steady state creep model of rock salt are given in Table 3 while the material properties for transient creep obtained for Avery Island rock salt (Yahya et al., 2000; Peric and Crook, 2004) are presented in Table 4. Also, elastic modulus, Poisson’s ratio and the density of the rock salt are the same as those given in the previous section. In Fig. 10a–c, the variation of the radial displacement, radial stress and tangential stress with the radial coordinate along line AB for p = 7 MPa at time t = 30 years are shown. Also plotted in these figures are the results obtained at time t = 0, i.e. those obtained by the elastic solution. In addition, radial stress Sxx, vertical stress Syy, tangential stress Szz, and shear stress Sxy for p = 7 MPa at time t = 30 years are illustrated in Fig. 11a–d. Also, the SFR contour and the displacement field for p = 7 MPa at time t = 30 years are presented in Figs. 12 and 13, respectively. As seen from the figures, stress components distribution and radial displacement distribution along line AB herein are qualitatively similar to those described in the previous section in the case of p = 0. It is worthy to note that, from the SFR contour, it can be observed that the highest values of the SFR are at the bottom of the cavern. However, the stress states shown for the case of p = 7 MPa do not cause the instability of the cavern. Based on the obtained results, as the LFR con- tour was equal to zero, it can be concluded that the stress states in the rock salt are in the compressibility domain in the case of p = 7 MPa. In another analysis, the potential effect of large deformations on the time-dependent response of salt caverns is investigated. With this aim, the center of the cavern is assumed to be at the depth of 1400 m below the ground surface while the cavern configuration is kept constant. The initial geometry, boundary conditions and initial conditions are the same as those represented in the previous analysis and the internal pressure is assumed to be p = 7 MPa. Also, the finite element mesh of the cavern is presented in Fig. 14. The simulation results comprising of the variation with time of the radial convergence at point A, the vertical convergence at point C, the vertical convergence at point D and the cavern relative closure, i.e. the ratio of the cavern volume reduction to the cavern initial volume, until the occurrence of the first short-term failure in one of the integration points are plotted respectively in Fig. 15a–d in the both cases of large displacements as well as small displacements. Also, the SFR and LFR contours, when the first short-term failure occurs at one of the integration points, are displayed respectively in Figs. 16a and b and 17a and b in the both cases of large displacements and small displacements. After excavation of the cavern, both Sm and s decrease with time due to the stress relaxation and if the stress state at an integration point is located in the dilatancy region, after some time a sudden failure at the integration point can be encountered by the stress state reaching the short-term failure surface. The amount of time which elapses until the occurrence of the short-term failure at an integration point depends on the octahedral shear overstress above S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 183 Fig. 16. SFR contours for the cavern excavated at the depth of 1400 m; (a) large displacements analysis and (b) small displacements analysis. the C–D boundary at the integration point. If the stress state is closer to the short-term failure surface at an integration point, the amount of time which elapses until the occurrence of the shortterm failure at the integration point will be shorter (Cristescu and Hunsche, 1998). According to the figures, in the present analysis, since the cavern is located at a greater depth in comparison with the previous analysis, the stress state after excavation of the cavern is closer to the short-term failure, and therefore the convergence of the cavern surface is faster here. Also, in the present analysis, as the LFR contour is not equal to zero in the neighborhood of the cavern surface, the stress state around the cavern surface is in the dilatancy region. From the SFR and LFR contours, it can be observed that the highest values of the SFR and LFR are located at the bottom of the cavern where the stress state is closer to the shortterm failure surface compared with other locations and the shortterm failure is encountered at this location. According to the figures, although, up to time around t = 500 days and for displacements up to around 1 meter, the values of cavity surface convergence and cavern relative closure obtained by the large displacements and small displacements analyses coincide, after time t = 500 days the results of these two analyses deviate and large displacements analysis shows faster convergence in comparison. Also, stress relaxation represented by large displacements analysis is more considerable than that of small displacements analysis which brings about an earlier shortterm failure for large displacements analysis than small displacements analysis. The values of SFR calculated by the large displacements analysis tend to localize in the cavern floor more than those of small displacements analysis. In addition, the time when the first short-term failure occurs at one of the integration points for large displacements analysis is equal to about 6 years while this value for small displacements analysis is equal to about 8 years. As it can be seen from the LFR contours, the calculated values of LFR around the cavern floor by the large displacements analysis is more considerable than those of small displacements analysis which means that in the case of large displacements analysis the volumetric viscoplastic strain work due to dilatancy and therefore the damage of the rock salt around the cavern floor is more noticeable than those of the small displacements analysis. However, in order to better estimate the time of the short-term failure occurrence, an implicit time integration scheme should be utilized. Also, after the occurrence of short-term failure in an integration point, the stress state at the integration point is located out of the domain where the Cristescu’s model is valid. Therefore, in order to better analyze the large displacements behavior of the salt cavern, the constitutive behavior of the rock salt beyond the occurrence of short-term failure, i.e. a fracture mechanic model of rock salt, should be added to the finite element model. 6. Conclusions The present paper has proposed a computational model to simulate the time-dependent behavior of underground storage caverns excavated in rock salt. An elasto-viscoplastic constitutive model describing dilatancy, short-term failure along with longterm failure during transient and steady state creep of rock salt, which is comprised of Peric’s viscoplastic power law representing transient creep and Cristescu’s steady state creep model, has been employed. The constitutive model is then implemented into a 184 S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 Fig. 17. LFR contours for the cavern excavated at the depth of 1400 m; (a) large displacements analysis and (b) small displacements analysis. Lagrangian finite element computer code, and afterwards in order to investigate the capability of the proposed computational algorithm, the stress variation and ground movement during creep of rock salt around a storage cavern has been analyzed numerically in the framework of large deformations and axisymmetric representation. According to the analyses, stress relaxation represented by Cristescu’s steady state creep model is more considerable in comparison with WIPP steady state creep model which may bring about a sooner short-term failure for Cristescu’s model than WIPP model, and also Cristescu’s model shows faster cavity surface convergence compared with WIPP model. To conclude, the computational model proposed herein represents a more conservative approach for stability analyses of the storage caverns than WIPP model which does not consider the irreversible volumetric changes. However, in order to deeper investigation of the presented computational model, further experimental data and measurements are required. Appendix A Matrixes A(h), H, and b S are defined by the following equations: 2 @u @x 60 6 6 60 6 AðhÞ ¼ 6 @u 6 @y 6 60 4 @u @z @v @x 0 3 @w @x 0 0 0 0 0 0 0 @u @y @v @y @w @y 0 0 07 7 7 @w 7 @z 7 7 07 7 @w 7 @y 5 0 0 0 0 0 @u @z @v @z @v @y @w @y @u @x @v @x @w @x 0 0 0 0 @u @z @v @z @w @z @u @y @v @y @v @z @w @z 0 @u @x @v @x 0 0 @w @x ðA1Þ 2 1 0 0 0 0 0 0 0 60 1 0 0 0 6 6 60 0 1 0 0 H¼6 60 1 0 1 0 6 6 40 0 0 0 0 0 2 0 1 0 3 07 7 7 0 0 0 17 7 0 0 0 07 7 7 1 0 1 05 0 0 0 ðA2Þ 0 0 0 0 1 3 Sx I Sxy I Sxz I 6 b S ¼ 4 Sxy I Sxz I Sy I 7 Syz I 5 Syz I 0 Sz I ðA3Þ where I is the 3  3 Identity matrix. References Aubertin, M., Gill, D.E., Ladanyi, B., 1991. A unified viscoplastic model for the inelastic flow of alkali halides. Mech. Mater. 11, 63–82. Bathe, K.J., 1996. Finite Element Procedures. Prentice-Hall, New Jersey. Chang, C., Zoback, M.D., 2010. Viscous creep in room-dried unconsolidated Gulf of Mexico shale (II): development of a viscoplasticity model. J. Petrol. Sci. Eng. 72, 50–55. Crisfield, M.A., 1991. Non-linear finite element analysis of solids and structures. Essentials, vol. 1. John Wiley & Sons, New York. Cristescu, N.D., Hunsche, U., 1992. Determination of a nonassociated constitutive equation for rock salt from experiments. In: Proceeding of IUTAM Symposium. Springer Verlag, Berlin, pp. 511–523. Cristescu, N.D., Hunsche, U., 1993. A constitutive equation for salt. In: Proceeding of 7th Int. Congr. on Rock Mech., Balkema, Rotterdam, pp. 1821–1830. Cristescu, N.D., 1993. A general constitutive equation for transient and stationary creep of rock salt. Int. J. Rock Mech. Min. Sci. 30, 125–140. Cristescu, N.D., Hunsche, U., 1998. Time Effects in Rock Mechanics. John Wiley & Sons, New York. S. Nazary Moghadam et al. / Tunnelling and Underground Space Technology 33 (2013) 171–185 Dawson, P.R., Munson, D.E., 1983. Numerical simulation of creep deformations around a room in a deep potash mine. Int. J. Rock Mech. Min. Sci. 20, 33–42. Fuenkajorn, K., Phueakphum, D., 2010. Effects of cyclic loading on mechanical properties of Maha Sarakham salt. Eng. Geol. 112, 43–52. Glamheden, R., Curtis, P., 2006. Excavation of a cavern for high-pressure storage of natural gas. Tun. Undergr. Sp. Technol. 21, 56–67. Hansen, F.D., Mellegard, K.D., 1980. Creep of 50-mm Diameter Specimens of Some Salt from Avery Island, Louisiana. RE/SPEC ONWI-104. Herrmann, W., Wawersik, W.R., Lauson, H.S., 1980a. Analysis of Steady State Creep of Southeastern New Mexico Bedded Salt. Sandia National Laboratories, SAND80-0558. Herrmann, W., Wawersik, W.R., Lauson, H.S., 1980b. Creep Curves and Fitting Parameters for Southeastern New Mexico Bedded Salt. Sandia National Laboratories, SAND80-0558. Heusermann, S., Rolfs, O., Schmidt, U., 2003. Nonlinear finite-element analysis of solution mined storage caverns in rock salt using the LUBBY2 constitutive model. Comput. Struct. 81, 629–638. Jin, J., Cristescu, N.D., 1998. An elastic viscoplastic model for transient creep of rock salt. Int. J. Plast. 14, 85–107. Li, Z., Liu, H., Dai, R., Su, X., 2005. Application of numerical analysis principles and key technology for high fidelity simulation to 3-D physical model tests for underground caverns. Tun. Undergr. Sp. Technol. 20, 390–399. Mellegard, K.D., Senseny, P.E., Hansen, F.D., 1981. Quasi-Static Strength and Creep Characteristics of 100-mm Diameter Specimens of Salt From Avery Island, Louisiana. RE/SPEC ONWI-250. Munson, D.E., Fossum, A.F., Senseny, P.E., 1990. Approach to first principles model prediction of measured WIPP (Waste Isolation Pilot Plant) in-situ room closure in salt. Tun. Undergr. Sp. Technol. 5, 135–139. Owen, D.R.J., Hinton, E., 1980. Finite Elements in Plasticity: Theory and Practice. Pineridge Press, Swansea. 185 Paraschiv-Munteanu, I., Cristescu, N.D., 2001. Stress relaxation during creep of rocks around deep boreholes. Int. J. Eng. Sci. 39, 737–754. Peric, D., Crook, A.J.L., 2004. Computational strategies for predictive geology with reference to salt tectonics. Comput. Methods Appl. Mech. Eng. 193, 5195–5222. Perzyna, P., 1966. Fundamental problems in viscoplasticity. Adv. Appl. Mech. 9, 243–377. Rokahr, R.B., Staudtmeister, K., Zander-Schiebenhofer, D., 1997. Development of a New Criterion for the Determination of the Maximum Permissible Internal Pressure for Gas Storage Caverns in Rock Salt. SMRI Research Project Report No. 1-94, IUB, Hannover. Shao, J.F., Zhu, Q.Z., Su, K., 2003. Modeling of creep in rock materials in terms of material degradation. Comput. Geotech. 30, 549–555. Sharifzadeh, M., Daraei, R., Sharifi Boroojerdi, M., 2012. Design of sequential excavation tunneling in weak rocks through findings obtained from displacements based back analysis. Tun. Undergr. Sp. Technol. 28, 10–17. Wang, G., Guo, K., Christianson, M., Konietzky, H., 2011. Deformation characteristics of rock salt with mudstone interbeds surrounding gas and oil storage cavern. Int. J. Rock Mech. Min. Sci. 48, 871–877. Xia, Y., Peng, S., Gu, Z., Ma, J., 2007. Stability analysis of an underground power cavern in a bedded rock formation. Tun. Undergr. Sp. Technol. 22, 161–165. Yahya, O.M.L., Aubertin, M., Julien, M.R., 2000. A unified representation of the plasticity, creep and relaxation behavior of rocksalt. Int. J. Rock Mech. Min. Sci. 37, 787–800. Zhu, W.S., Sui, B., Li, X.J., Li, S.C., Wang, W.T., 2008. A methodology for studying the high wall displacement of large scale underground cavern complexes and its applications. Tun. Undergr. Sp. Technol. 23, 651–664. Zienkiewicz, O.C., Cormeau, I.C., 1972. Visco-plastic solution by finite element process. Arch. Mech. 24, 873–889. Zienkiewicz, O.C., Cormeau, I.C., 1974. Visco-plasticity–plasticity and creep in elastic solids – a unified numerical solution approach. Int. J. Numer. Methods Eng. 8, 821–845.