Geophysical Journal International Geophys. J. Int. (2017) 211, 252–262 Advance Access publication 2017 July 14 GJI Rock and mineral physics, rheology doi: 10.1093/gji/ggx297 A comparative study of Maxwell viscoelasticity at large strains and rotations Christoph E. Schrank,1,2 Ali Karrech,3 David A. Boutelier4 and Klaus Regenauer-Lieb5 1 School of Earth, Environmental and Biological Sciences, Queensland University of Technology, Brisbane, QLD, Australia. E-mail: christoph.schrank@qut.edu.au 2 School of Earth and Environment, University of Western Australia, Crawley, WA, Australia 3 School of Civil, Environmental and Mining Engineering, University of Western Australia, Crawley, WA, Australia 4 School of Environmental and Life Sciences, University of Newcastle, Callaghan, NSW, Australia 5 School of Petroleum Engineering, University of New South Wales, Sydney, NSW, Australia SUMMARY We present a new Eulerian large-strain model for Maxwell viscoelasticity using a logarithmic co-rotational stress rate and the Hencky strain tensor. This model is compared to the smallstrain model without co-rotational terms and a formulation using the Jaumann stress rate. Homogeneous isothermal simple shear is examined for Weissenberg numbers in the interval [0.1; 10]. Significant differences in shear stress and energy evolution occur at Weissenberg numbers >0.1 and shear strains >0.5. In this parameter range, the Maxwell–Jaumann model dissipates elastic energy erroneously and thus should not be used. The small-strain model ignores finite transformations, frame indifference and self-consistency. As a result, it overestimates shear stresses compared to the new model and entails significant errors in the energy budget. Our large-strain model provides an energetically consistent approach to simulating non-coaxial viscoelastic deformation at large strains and rotations. Key words: Creep and deformation; Fault zone rheology; Rheology: crust and lithosphere; Rheology and friction of fault zones. 1 I N T RO D U C T I O N The Maxwell body of a linear-elastic Hookean spring in series with a Newtonian dashpot is the simplest rheological model for geological deformation (e.g. Bailey 2006). It has been employed to describe many deformation processes and structures, for example mantle convection (Mühlhaus & Regenauer-Lieb 2005; Beuchert & Podladchikov 2010), Rayleigh–Taylor instabilities (Kaus & Becker 2007; Robin & Bailey 2009), lithospheric deformation (Kusznir & Bott 1977; Kaus & Podladchikov 2006; Korenaga 2007), postglacial rebound (Peltier 1976; Larsen et al. 2005), shear zones (Branlund et al. 2000; Mancktelow 2002; Sone & Uchide 2016) and folds (Mancktelow 1999; Zhang et al. 2000; Liu et al. 2016). In spite of their simplicity, current Maxwell models display well-known errors (e.g. Mühlhaus & Regenauer-Lieb 2005; Beuchert & Podladchikov 2010) when the associated strains and, importantly, rotations are large. Such conditions are met, for example, during plate motions (Tohver et al. 2010), in bending subducting plates (Conrad & Hager 1999), blocks of crust bounded by faults (Abels & Bischoff 1999), folds (Ramsay & Huber 1987), shear zones (Ramsay 1980), boudins (Mandal & Khan 1991), porphyroblasts (Johnson 2009) and–clasts (Passchier & Simpson 1986), and grains deforming by diffusion creep (Wheeler 2010). Large rotations pose a mathematical challenge when elasticity is considered in the rheology of highly 252  C transformed materials (Lehmann 1972a; Dienes 1979; Moss 1984; Xiao et al. 1997a, 1998a,b, 2006; Bruhns et al. 1999, 2001a,b; Colak 2004; Meyers et al. 2005; Karrech et al. 2011, 2012): one requires an objective formulation of the stress rate (time derivative of stress). If a volume element of a transforming body incurs a rigid-body rotation over an infinitesimal time increment, its state of stress and strain should not change relative to a coordinate system that corotates with this volume element. In other words, the notion of objectivity, or frame indifference, requires that the co-rotational stress and strain increments become zero during an instantaneous rigid-body rotation (see e.g. chap. 1.6 of Chakrabarty 2006). In an Eulerian (stationary) reference frame, the objectivity requirement is commonly addressed with a so-called co-rotational objective stress rate, which contains spins that account for elastic rotations (for an in-depth review, see Xiao et al. 2006). Here, we compare the self-consistent formulation of the Eulerian finite-transformation model for Maxwell viscoelasticity (called FT model hereafter) to two thermodynamically non-consistent Maxwell models, which are commonly used to model geologically relevant strains (e.g. Mancktelow 1999; Zhang et al. 2000; Moresi et al. 2002; Kaus & Becker 2007; Beuchert & Podladchikov 2010): the small-strain formulation without co-rotational terms, referred to as SS model, and the formulation using the Jaumann stress rate, called MJ model The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. Downloaded from https://academic.oup.com/gji/article/211/1/252/3965336 by guest on 30 January 2024 Accepted 2017 July 13. Received 2017 April 24; in original form 2017 January 31 Viscoelasticity at large transformations 2 M A X W E L L V I S C O E L A S T I C I T Y, C O - R O TAT I O N A L S T R E S S R AT E S A N D THE NEED FOR A FINITE MAXWELL MODEL The constitutive equation for Maxwell viscoelasticity of an incompressible material is Di j = σi′j σ̇i′j  1 L i j + L iTj = + , 2 2G 2η (1) where Di j is the stretching or deformation-rate tensor, defined as the symmetric part of the velocity-gradient tensor L i j = ∂vi /∂ x j with velocities vi and directions xi , the superscript T marks the matrix transpose, the overdot denotes a derivative with respect to time t, G is the shear modulus of an isotropic linear-elastic material, σi′j is the deviator (indicated by the dash) of the Cauchy stress tensor σi j , and η is the Newtonian viscosity. In the limit of infinitesimal deformations, the stretching tensor becomes the familiar strain-rate tensor for small strains. Eq. (1) states that elastic (first term) and viscous deformation rates (second term) are additive (serial coupling of an elastic spring with a viscous dashpot). The term σ̇i′j is the stress rate, which highlights the need for a co-rotational formulation in order to acknowledge rotations at finite strain. Formulations that include geometrically admissible co-rotational terms are called ‘objective’. Eq. (1) is solved by direct integration for the SS model while σ̇i′j is replaced with rates augmented by specific spins in the MJ and FT models (Section 3). It is obvious that the SS model ignores large rotations and strains and hence is not objective and not thermodynamically self-consistent. Errors are expected at large non-coaxial deformations. The magnitude of these errors can be estimated by comparing the SS model with other existing models such as MJ and now FT. The MJ model employs a classical objective co-rotational stress rate, the Jaumann rate (Jaumann 1911). Two important shortcomings of the MJ model under simple shear have been recognized in the geodynamic literature previously (Mühlhaus & Regenauer-Lieb 2005; Beuchert & Podladchikov 2010). To summarize these findings, we first recall the Maxwell relaxation time trel and the Weissenberg number W: trel = η G W = γ̇ trel, (2a) (2b) where γ̇ is a characteristic shear-strain rate, which in our study corresponds to the imposed simple-shear rate defined in eq. (4). Deformation observed at time scales well below trel is dominated by elastic contributions to the strain rate and stresses. After about five trel , deformation is essentially viscous. The Weissenberg number W measures to which degree the material behaviour is nonlinear and how much of the current strain is recoverable (White 1964; Dealy 2010). Mühlhaus & Regenauer-Lieb (2005) compared Maxwell models using the Jaumann and Green–Naghdi rates, two objective stress rates from the continuum-mechanics community, with the SS model frequently used in geological modelling. While these objective formulations contain the necessary geometric terms that ensure frame indifference, they do not satisfy the self-consistency condition of Bruhns et al. (1999). Consequently, for the same material properties, the Green–Naghdi and MJ models give different stress-strain curves at shear strains >1. The SS model (used in geology) and the Green–Naghdi model (used in some engineering applications, e.g. for car-crash simulations; Abaqus/Explicit 2015) produce qualitatively similar results but the MJ model (used in other engineering application, e.g. for calculations of the nonlinear creep of the steel carapace of a nuclear reactor under elevated temperatures; Abaqus/Standard 2009) exhibits spurious softening during numerical simulations of simple shear for W = 1. Mühlhaus & Regenauer-Lieb (2005) showed, however, that this effect vanishes if non-linear creep and/or plasticity are added to the constitutive model. Later, Beuchert & Podladchikov (2010) demonstrated numerically that the MJ model exhibits unphysical shear-stress oscillations at W = 10. They also showed that the choice, and even the neglect, of co-rotational terms become irrelevant at W ≤ 0.1. Therefore, two regimes of intermediate (0.1 < W < 1) and high (W ≥ 1) Weissenberg numbers exist in which one can expect problems with the SS and MJ models. The question arises: where in the Earth may the intermediate and high-W regimes occur? Beuchert & Podladchikov (2010) and Mühlhaus & RegenauerLieb (2005) aiming for modelling viscoelastic mantle convection argued that only the lithosphere exhibits sections with W ≫ 0.1. At the tectonic time and length scales of interest to these studies, the lithosphere behaves essentially elastic and hardly deforms. Hence, the problems of the MJ model, spurious softening and stress oscillations, were predicted to be irrelevant in viscoelastic mantle simulations. This argument is further supported by the notion that power-law flow and plasticity, both common in the lithosphere Downloaded from https://academic.oup.com/gji/article/211/1/252/3965336 by guest on 30 January 2024 hereafter. Self-consistency (i.e. adherence to thermodynamic principles) in the FT model is achieved by using the logarithmic spin introduced by Xiao et al. in the context of finite elastoplasticity (Xiao et al. 1997a, 1998a,b, 2006; Bruhns et al. 1999; Bruhns 2003). The choice of the logarithmic spin in the new FT model is motivated by two reasons. First, the logarithmic co-rotational rate of the Hencky strain was proven to be the only strain-rate measure, which is always identical with the stretching (Xiao et al. 1997a,b; Bruhns et al. 1999). The stretching, the symmetric part of the velocity gradient tensor, is the most commonly used and kinematically most appealing measure of deformation rates (cf. chapter 2 of Xiao et al. 1997b). Second, the logarithmic stress rate is the only choice, which delivers a self-consistent Eulerian constitutive model (Bruhns et al. 1999). Bruhns et al. (1999) prove self-consistency in the co-rotational Eulerian rate model through integrating the additive deformation-rate decomposition for the case of purely elastic stretching. If the rate formulation of the model is exactly integrable to deliver a truly elastic response, that is, fully reversible deformation without energy dissipation, then the Eulerian finite-strain model is defined as (thermodynamically) self-consistent. Only models that use the logarithmic spin are self-consistent. To analyse the relevance of the co-rotational logarithmic spin used in the FT model, homogeneous 2-D isothermal simple shear is examined because of its non-zero, skew-symmetric spin and high conceptual importance for understanding shear zones (e.g. Ramsay & Graham 1970; Ramsay 1980). The relative stress and energy states are evaluated for a shear-strain range of 0–10, providing the basis for comparison of the SS, MJ and FT models. The sections to come are organized as follows: first, the need for a proper finite Maxwell model is established. Second, analytical solutions for simple shear with the SS and MJ models are presented. Third, the FT model is derived and its numerical solution described. Fourth, all three models are compared. Finally, the results and their implications are discussed. 253 254 C.E. Schrank et al. G = 10 GPa W = G = 100 GPa W 1 = 1 FT only W FT only W = 0. = 1 SS 0. 1 FT (Pa.s) (Pa.s) SS FT (a) -1 Strain rate (s ) -1 Strain rate (s ) (b) Figure 1. Double-logarithmic contour plots of the decadic logarithm of W over those of viscosity and strain rate for the (a) lower- and (b) upper-bound estimate of the shear modulus. The white dotted lines denote parameter combinations with W = 0.1; the white dashed lines mark those with W = 1. Textboxes indicate the models, which are valid in the respective parameter fields separated by those lines. (e.g. Karato 2008), serve as stress limiters, which inhibit spurious softening and oscillations (Beuchert & Podladchikov 2010). However, we argue that deformation problems within the lithosphere, in particular at length scales ≤1 km, require the tackling of the problem of modelling Maxwell viscoelasticity in the intermediate and high-W regimes properly. The shear moduli of most rock-forming minerals are within the range of 10–100 GPa (Bass 1995). The effective viscosity of lithospheric rocks varies over several orders of magnitude due to its sensitivity to, for example, material composition, temperature, pressure, strain rate, stress, and grain size (e.g. Karato 2008). We adopt lower and upper bounds of 1018 and 1025 Pas, respectively. Shear-strain rates also cover several orders of magnitude. A typical lower bound for natural shear zones is 10−15 s−1 (Pfiffner & Ramsay 1982; Biermeier & Stüwe 2003). Estimates for the upper bound in nature are of order 10−10 s−1 for shear zones with widths ≤100 m (White & Mawer 1992). Computing W for the above parameter ranges shows that a significant subset (∼35 per cent to 50 per cent) of this parameter space is within the critical regimes of W (Fig. 1). Moreover, one can envisage geological scenarios during which the internal deformation of the lithosphere remains elastic, or transient viscoelastic, in large parts while the associated rotations become large, for example, in subducting plates and blocks of crust rotating about a vertical axis. In these cases, a proper objective treatment of large transformations, and elastic rotations in particular, seems desirable. This effort should also involve stress advection (e.g. Truesdell 1956), which is irrelevant in the simple-shear case at hand. 3 D E R I VAT I O N S : K I N E M AT I C S O F 2 - D SIMPLE SHEAR We consider 2-D plane-strain isothermal simple shear of a unit square, where x1 and x2 denote the horizontal and vertical direction, respectively. The lower-left corner of the square is centred on the origin of the coordinate system. Shear is parallel to x1 . The boundary conditions are v1 (x2 = 0) = v2 (x2 = 0) = 0 v1 (x2 = 1) = const. > 0 v2 (x2 = 1) = 0 ∂v1 (0 < x2 < 1) = const. > 0, ∂ x2 (3) where vi denote the velocities in the respective directions. Thus, the velocity-gradient tensor Lij has only one non-zero, constant entry: L 12 = ∂v1 = γ̇ , ∂ x2 (4) where γ̇ is the constant shear-strain rate applied to the unit square. Hence, the deformation-rate tensor Di j and the spin tensor Si j , the antisymmetric component of L i j , read ⎤ ⎡ γ̇ ⎢0 2⎥ ⎥ Di j = ⎢ ⎦ ⎣ γ̇ 0 2 ⎤ ⎡ (5a, b) γ̇ 0 ⎥ ⎢   1 2⎥ L i j − L iTj = ⎢ Si j = ⎦. ⎣ γ̇ 2 0 − 2 Eq. (5b) expresses the internal rotations inherent to simple shear. Conceptually, the simple-shear model embodies an idealized plane-strain shear zone with homogeneous strain distribution. In nature, shear zones generally display a heterogeneous strain pattern (Ramsay & Graham 1970; Ramsay 1980). However, this heterogeneity is a function of the scale of observation. One can identify length scales, at which deformation is, to first order, homogeneous. Hence, one may consider the following analysis as an examination of a homogeneously deforming representative elementary volume Downloaded from https://academic.oup.com/gji/article/211/1/252/3965336 by guest on 30 January 2024 SS = MJ = FT SS = MJ = FT Viscoelasticity at large transformations within an ideal shear zone, across which variations in parameters such as temperature, composition, viscosity, fluid content, etc., are negligible. 255 Inserting eq. (13) into the constitutive model (eq. 1) with the deformation-rate tensor for simple shear (eq. 5a) and substituting trel yield the system of coupled differential equations: ηγ̇ = trel (σ̇12 − γ̇ σ̇22 ) + σ12 3.1 Derivations: Maxwell model for small strain (SS model) 0 = trel (σ̇22 + γ̇ σ12 ) + σ22 . The SS model is solved by integrating eq. (1) after substituting simple-shear kinematics (eqs 3 and 5a), trel (eq. 2a), W (eq. 2b) and −1 t , which yields dimensionless time k = trel   −k (6a) σ12 (k) = W G 1 − e σ11 = σ22 = 0. (6b) Ė TOT (t) = σi′j : Di j = σ12 γ̇ . (7) Integrating eq. (7) over time with the initial condition E TOT (t = 0) = 0 gives the total energy per unit volume:   (8) E TOT (k) = W 2 G k + e−k − 1 . To decompose the total energy into elastic and viscous fractions, we use the additivity of the deformation rates to introduce the elastic and viscous fraction of strain rate (e.g. Kaus & Becker 2007) 1= γ̇vis γ̇el + = f el + f vis , γ̇ γ̇ (9) where fel and fvis denote the elastic (subscript el) and viscous (subscript vis) fraction of the strain rate, respectively. At infinite time, the ∞ = σ12 (k → ∞) = ηγ̇ = W G. shear stress is purely viscous: σvis Thus, the viscous and elastic stress fractions can be defined as σ12 (k) = 1 − e−k f vis = ∞ σvis f el = 1 − f vis = e−k . (10a) (10b) Inserting eqs (6), (9) and (10) into eq. (7) and integrating over time yield the stored and viscous energy per unit volume, again for zero initial energies: k f el σ12 dk = W 2 G E el = γ̇ 1 1 −2k + e − e−k 2 2 0 k 2 f vis σ12 dk = W G 2e E vis = γ̇ 0 −k  = 2 σ12 2G  1 −2k 3 . − e +k− 2 2 WG (15a) e−k [W sin (W k) − cos (W k)] + 1 1 + W2 WG σ22 (k) = (15b) e−k [W cos (W k) + sin (W k)] − W . 1 + W2 The first fraction in eq. (15a), b is the steady-state viscous shear ∞ stress for the MJ model: σvis = σ12 (k → ∞) = W G(1 + W 2 )−1 . It differs from the steady-state viscous stress of the SS model by the factor (1 + W 2 )−1 . This factor highlights an important effect of the stress-rotation terms included in the Jaumann stress rate: they decrease the maximum shear stress attained at steady state. Using eqs (7) and (15a), the total deformational energy per unit volume as function of k is  e−k ∞ ((W 2 − 1) cos(W k)+2W sin(W k)) E TOT (k) = σvis W k− 2 W +1  W2 − 1 + 2 . (16) W +1 σ12 (k) = Eq. (16) is also decomposed into elastic and viscous fractions. The elastic and viscous stress fractions read: f el = e−k (cos(W k) − W sin(W k)) (17a) f vis = e−k (W sin(W k) − cos(W k)) + 1. (17b) To conclude the energy decomposition, the time integrals of eq. (11) need to be solved again after substituting the appropriate expressions for the shear stress (eq. 15a) and elastic and viscous stress fractions (eq. 17). The results are cumbersome because of the lengthy expressions resulting from the integration of the trigonometric terms and are thus not shown here. Many mathematics software packages can conduct this integration, and the interested user may find this more convenient. (11a) 3.3 Derivations: Maxwell model at finite strain (FT model) (11b) 3.2 Derivations: Maxwell model with the Jaumann stress rate (MJ model) The elastic stress rate in eq. (1) (σ̇i′j ) is now replaced by the Jaumann stress rate (Jaumann 1911): τ̇i j = σ̇i j − Sik σk j + σik Sk j . The solution of eq. (14) in terms of W and dimensionless time k is (Appendix A1): (12) This stress rate is equivalent to a rigid-body rotation of the stress tensor neglecting second and third-order terms (Chakrabarty 2006). Since in the problem at hand σ11 = −σ22 ⇒ tr(σi j ) = 0 (e.g. Durban 1990), the dashes for deviatoric quantities are dropped in the following. Substituting the spin (eq. 5b), eq. (12) reads:   σ̇11 − γ̇ σ12 σ̇12 + γ̇ σ11 . (13) τ̇i j = σ̇21 + γ̇ σ11 σ̇22 + γ̇ σ12 The logarithmic finite-strain model requires slight modifications of eqs (1) and (12). The Cauchy stress σi j is replaced by the Kirchhoff stress i j = J σi j where J is the determinant of the deformationgradient tensor. We use the Kirchhoff stress with more general applications in mind, even though J = 1 for simple shear. The co-rotational Cauchy stress rate τ̇i j (eq. 12) is replaced by the co⌣ rotational rate of the Kirchhoff stress  i j ⌣ ˙ ij − i j =  ik k j + ik kj , (18) which involves a logarithmic spin of rotation discovered by Xiao and colleagues (Xiao et al. 1997a,b):  n  λ A +λ B 2 (A) (B) Pik Dkl Pl j , + i j =Wi j + λ A −λ B ln λ A − ln λ B A,B=1,A=B (19) where A and B are dummy integers denoting the principal directions, λ A is the Ath eigenvalue of the left Cauchy–Green strain Downloaded from https://academic.oup.com/gji/article/211/1/252/3965336 by guest on 30 January 2024 The total deformational power per unit volume is (14) 256 C.E. Schrank et al. (A) model computes lower energies than the FT model with a maximum difference at γ = 10, being about 1 per cent at W = 0.1 up to ca. 32 per cent at W = 0.9. At γ = 2, the difference is <5 per cent in all cases. In the evolution of stored energy, the SS model predicts larger elastic energy than the FT model, ranging from ca. 2 per cent at W = 0.1 and γ ≈ 0.5 to about 56 per cent at W = 0.9 and γ ≈ 2.3 (Fig. 3c). The curve shapes mirror those of the shear-stress evolution because elastic energy and shear stress are related by 2 G −1 (eq. 11a). The MJ model gives smaller energies E el = 0.5σ12 than the FT-model energies with the maximum difference appearing at γ = 10 for W ≤ 0.7, ranging from ca. 5 per cent at W = 0.1 to ∼ 95 per cent at W = 0.7. At larger W, negative ratios are obtained because the MJ theory computes negative stored energy there. In terms of dissipated viscous energy (Fig. 3d), the SS model first gives smaller and then larger values than the FT model. The difference of the smaller values does not exceed 10 per cent (at W = 0.9 and γ = 0.7). The difference of the lager values reaches ca. 35 per cent at the largest W at small γ . The MJ model exhibits the opposite behaviour with considerably larger magnitudes of discrepancy, attaining an excess of up to ∼ 260 per cent and a deficiency of ca. 30 per cent for the largest W. 4 R E S U LT S We present the evolution of normalized shear stress σ12 G −1 and the ratios of shear stress and energies up to a shear strain of γ = 10 for the intermediate- and high-W regimes. The former is represented by the discrete set of solutions with W ∈ {0.1; 0.2; . . . ; 0.9}, and the latter by the set W ∈ {1; 2; . . . ; 10}. 4.1 Shear-stress evolution The shear-stress curves start to diverge at γ ≈ 0.5 (Fig. 2). In general, shear stresses of the SS model and those of the MJ model bracket shear stresses of the FT model. Their discrepancy increases with increasing W. The spurious strain softening of the MJ model is significant for 0.5 ≤ W ≤ 1 while stress oscillations are visible for all greater W. The SS and FT models exhibit monotonously increasing shear stresses over the tested strain range. We continue with a detailed comparison of relative shear stresses and energies, considering the intermediate- and high-W regimes separately (Fig. 3). For the high-W regime, MJ results are omitted because their stress oscillations are highly problematic (see Fig. 2b and Appendix A3 for an in-depth analysis). 4.2 Intermediate-W regime The SS model predicts higher stresses than the FT model (Fig. 3a). The maximum magnitude of this discrepancy increases from ca. 1 per cent at W = 0.1 and γ ≈ 0.5 to ca. 25 per cent at W = 0.9 and γ ≈ 2.3. At γ ≤ 0.5, the difference never exceeds 5 per cent. The MJ model predicts lower stresses than the FT model (Fig. 3a). The maximum discrepancy increases nonlinearly with shear strain and W and reaches from <2 per cent at W = 0.1 to ca. 42 per cent at W = 0.9 and γ = 10. The difference does not exceed 3 per cent at γ ≤ 1. Since the total energy is the integral of the stress-strain curve, the evolutions of the relative total energies resemble those of the relative shear stresses (Fig. 3b). The SS model computes larger total energies than the FT model, ranging from ca. 1.5 per cent excess at W = 0.1 and γ ≈ 1.3 to ca. 20 per cent at W = 0.9 and γ ≈ 3.8. At γ = 1, the difference is <7 per cent. Again, the MJ 4.3 High-W regime The overall trends of shear-stress evolution remain similar to those of the intermediate-W regime. The SS model gives shear stresses higher than the FT model with a maximum discrepancy of ca. 30 per cent at W = 1 to ca. 350 per cent at W = 10 (Fig. 3e). This trend is reflected in the total-energy comparison (Fig. 3f), which reveals a maximum difference of ∼25 per cent at W = 1 up to ∼280 per cent at W = 10. A significant proportion of this energy error is due to problems with the computation of elastic energy (Fig. 3g). The SS model computes larger stored energy than the FT model by up to ca. 50 per cent at W = 1 up to more than one order of magnitude at W = 10. Accordingly, the differences in viscous-energy evolution are smaller (Fig. 3h). The largest negative difference ranges from ca. 10 per cent at W = 1 to about 70 per cent at W = 10; the largest positive difference covers a range of ∼13 per cent at W = 1 to ∼65 per cent at W = 10. 5 DISCUSSION The three examined models for viscoelastic deformation differ markedly. The SS model ignores large transformations, frame indifference, and self-consistency (sensu Bruhns et al. 1999) a priori. As a result, the FT and SS model differ most significantly during highly non-coaxial large deformations in the intermediateand high-W regimes. The MJ formulation performs worse in some cases than the SS model in spite of its use of a geometrically correct objective stress rate. However, it violates the self-consistency requirement, which incurs unphysical stress oscillations and negative stored energy (Fig. 3c, see also Appendix A3). This poses a serious problem for multi-physics formulations with energy feedbacks that rely on accurately tracking cross-scale energy fluxes through thermodynamic upscaling methods (Regenauer-Lieb et al. 2013, 2014). The FT model provides an energetically consistent large-transformation model that satisfies the objectivity and selfconsistency requirements. We now proceed to discuss the applicability of the examined Maxwell models in the light of our results. Downloaded from https://academic.oup.com/gji/article/211/1/252/3965336 by guest on 30 January 2024 tensor, Pik is its Ath eigenprojection (tensor product of an eigenvector with itself), and n is the number of distinct eigenvalues. The sum vanishes naturally if the eigenvalues are all equal. The spin of eq. (19) contains a correction term compared to that of the Jaumann formulation, which accounts for the stretching of the eigenprojections, not only their rotation. The logarithmic spin produces a smooth, stable response under simple shear, unlike the Jaumann spin, which produces aberrant oscillations (Lehmann 1972b; Dienes 1979). Moreover, since the co-rotational logarithmic rate of the Hencky strain tensor coincides with the deformation-rate tensor, the former constitutes an ideal energy conjugate of the Kirchhoff stress (Bruhns 2003). Therefore, the Maxwell model stated in eq. (1) involves another modification in the FT model, in addition to inserting the stress rate of eq. (18): the time derivative of the Eulerian logarithmic Hencky strain tensor replaces the stretching. The associated system of differential equations (eq. A9) and the numerical method for integrating them are detailed in Appendix A2. The particular qualities, which render logarithmic kinematics the most mathematically consistent tool for the constitutive behaviour of geomaterials, are further discussed in great depth by Xiao et al. (2006). Viscoelasticity at large transformations (a) 0.2 W = 0.2, k max = 50 257 W = 0.6, k max = 16.67 FT SS MJ 0.1 0 0.5 0 W = 0.7, k max = 14. 2 9 12 G -1 W = 0.3, k max = 33.33 0.5 0.2 0.1 0 0 W = 0.4, k max = 25 W = 0.8, k max = 12. 5 0.5 0.2 0 0 W = 0.9, k max = 11.11 0.4 0.5 0.2 0 0 0 2 (b) 4 6 8 10 0 2 W = 1, k max = 10 4 6 8 10 8 10 W = 6, k max = 1.67 4 0.5 FT SS 2 MJ 0 0 W = 2, k max = 5 W = 7, k max = 1. 4 3 4 1 2 0 0 W = 8, k max = 1.25 4 2 12 G -1 W = 3, k max = 3.33 2 1 0 0 W = 4, k max = 2. 5 W = 5, k max = 2 2 0 2 4 6 W = 10, k max = 1 6 4 2 0 4 0 W = 9, k max = 1. 11 6 4 2 0 3 2 1 0 8 10 0 2 4 6 Figure 2. Plots of normalized shear stress over shear strain for the SS, MJ and FT models in the (a) intermediate-W regime and (b) the high-W regime. Each plot title indicates the respective W and the dimensionless time kmax reached at final γ . Note that the results for W = 0.1 are omitted because the curves cannot be distinguished at this scale of observation. With engineering applications in mind, our comparison reveals a range of W where all discussed models deliver qualitatively similar results. At small strains and W < 0.1, the three tested models are identical. Under these conditions, the SS model provides the simplest choice. At intermediate W and γ < 1, the shear stresses of all three models are within ≤ ±10 per cent. If this accuracy is acceptable in this strain and W range, it does not matter practically which model is used. At γ > 1, the shear-stress deviations exceed 10 per cent at W > 0.3 in the MJ model, and at W > 0.4 in the SS model. Nevertheless, the shear-stress overestimation of the SS model does not exceed 25 per cent in the intermediate-W regime. One may perhaps envisage applications, where for the sake of minimizing the computational burden, this overestimate is acceptable. The MJ model exhibits spurious softening of the shear stresses at W ≥ 0.3 (Fig. 2a). This unphysical effect renders its use problematic at and beyond this W. In the high-W regime, shear stresses in the SS model exceed those of the FT model by ca. 10 per cent at γ ≈ 0.5 and keep on increasing quickly while the MJ model oscillates unrealistically (Fig. 2b). As it is clear that the formulation of the SS model is ill suited for large transformation, the rapid increase and high stress magnitude of the SS model in high-W regime cannot be considered acceptable. Therefore, we consider the FT model the best choice in the high-W regime. Downloaded from https://academic.oup.com/gji/article/211/1/252/3965336 by guest on 30 January 2024 W = 0.5, k max = 20 C.E. Schrank et al. Shear-stress ratio 1.25 1 SS FT 1.2 Total-energy ratio 258 0.9 1.15 0.8 1.1 0.7 1.05 MJ 0.6 FT 1.5 0.8 1.4 0.6 1.3 0.4 1.2 0.2 1.1 0 (c) 1 0 MJ -0.2 FT 1 2 3 4 5 6 7 8 9 10 0.9 1.1 0.8 1.05 MJ 0.7 FT (b) Viscous-energy ratio SS FT 1.15 0 1 2 3 4 5 6 7 8 9 10 MJ FT SS FT 1.3 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 2 1.2 1.5 1.1 1 (d) 1 0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10 10 SS 9 FT 8 7 6 5 4 3 2 1 3 2.5 2 1.5 Total-energy ratio Shear-stress ratio 3.5 (e) 2 1.5 (f) 12 SS FT 10 Viscous-energy ratio Elastic-energy ratio SS FT 2.5 8 6 4 SS FT 1.5 1 0.5 2 (g) 0 1 2 3 4 5 6 7 8 9 10 (h) 0 1 2 3 4 5 6 7 8 9 10 Figure 3. Top panel (light-grey canvas): comparison of model results in the intermediate-W regime. (a) Shear-stress ratio, (b) total-energy ratio, (c) storedenergy ratio and (d) viscous-energy ratio for the SS model (solid lines) and the MJ model (dashed lines). Each subfigure contains two plots: on the left, SS results are divided by FT results; on the right, MJ results are divided by FT results. Curve symbols denote the respective W, defined in the right panel of (d). Bottom panel (grey canvas): comparison of the SS to the FT model in the high-W regime. Results of the SS model are again divided by those of the FT model. The legend for the respective value of W is shown in (e). However, if energy consistency matters, the FT model is also better suited for most of the intermediate-W regime. As emphasized in a review by Hobbs et al. (2011), the time evolution of the stored energy of a coupled geological system holds the key to understanding dissipation and thus the formation of localized structures. This notion is deeply rooted in thermodynamic principles (Ziegler 1977; Coussy 2004). Therefore, if one wishes to tackle non-coaxial deformation problems, where the mechanical energy of the system is coupled to non-mechanical dissipation processes such as metamor- phic processes or heat transfer, energy consistency is important. The SS model overestimates stored mechanical energy in simple shear by well over 10 per cent at W > 0.3. Therefore, we recommend, as a pragmatic rule of thumb, to use the SS model in the intermediate-W regime only in cases where energy consistency is not required or W ≤ 0.3. The MJ model exhibits unphysical negative stored energy at W > 0.7, which obviously violates thermodynamic principles. In fact, all MJ models exhibit negative time gradients in stored Downloaded from https://academic.oup.com/gji/article/211/1/252/3965336 by guest on 30 January 2024 Elastic-energy ratio (a) 1 1 SS FT Viscoelasticity at large transformations 6 C O N C LU S I O N S We introduced a large-strain model for Maxwell viscoelasticity with a logarithmic co-rotational stress rate (the ‘FT model’). An analysis of homogeneous isothermal simple shear with the FT model compared to a small-strain formulation (the ‘SS model’) and a model using the Jaumann stress rate (the ‘MJ model’) leads to the following key conclusions: (1) At W ≤ 0.1, all models yield essentially identical results. (2) At larger W, the models show increasing differences for γ > 0.5. The SS model overestimates shear stresses compared to the FT model while the MJ model exhibits an oscillatory response underestimating the FT model. (3) The MJ model violates the self-consistency condition resulting in stress oscillations and should be disregarded. It does not deliver truly elastic behaviour. (4) In the intermediate-W regime, the shear-stress overestimates of the SS model may constitute acceptable errors if energy consistency is not important. If energy consistency is desirable, the SS model should not be used at W ≥ 0.3. (5) In the high-W regime, stresses in the SS model become unacceptably large. The FT model should be used in this domain. The FT model constitutes a physically consistent Maxwell model for large non-coaxial deformations, even at high Weissenberg numbers. It overcomes the conceptual limitations of the SS model, which is limited to small transformations, not objective and not self-consistent. It also solves the problem of the energetically aberrant oscillations of the MJ model. The FT model promises particular usefulness for simulations of bending plates, rotating blocks of crust, and lithospheric viscoelastic shear zones, especially those at the scale of hundreds of meters or smaller as well as those involving strong heterogeneity of material properties due to, for example, layering or inclusions. In such systems, high Weissenberg numbers and large rotations can be expected. AC K N OW L E D G E M E N T S The authors gratefully acknowledge funding by the Australian Research Council via Discovery project DP140103015. CS is grateful for support of the QUT High Performance Computing group. We are very thankful for the thorough and most helpful reviews of Greg Houseman, an anonymous reviewer, and the editor Jörg Renner. REFERENCES Abaqus/Standard, 2009. User’s Manual, Version 6.9, Hibbit, Karlsson, and Sorensen Inc., Pawtucket, Rhode Island, 2000. Abaqus/Explicit, 2015. User’s Manual, Version 6.14, Dassault Systèmes. Abels, A. & Bischoff, L., 1999. Clockwise block rotations in northern Chile: indications for a large-scale domino mechanism during the middle-late Eocene, Geology, 27, 751–754. Bailey, R.C., 2006. Large time step numerical modelling of the flow of Maxwell materials, Geophys. J. Int., 164, 460–466. Bass, J.D., 1995. Elasticity of minerals, glasses, and melts, in Mineral Physics and Crystallography: A Handbook of Physical Constants, pp. 45–63, ed. Ahrens, T.J., American Geophysical Union. Beuchert, M.J. & Podladchikov, Y.Y., 2010. Viscoelastic mantle convection and lithospheric stresses, Geophys. J. Int., 183, 35–63. Biermeier, C. & Stüwe, K., 2003. Strain rates from snowball garnet, J. Metamorphic Geol., 21, 253–268. Branlund, J.M., Kameyama, M.C., Yuen, D.A. & Kaneda, Y., 2000. Effects of temperature-dependent thermal diffusivity on shear instability in a viscoelastic zone: implications for faster ductile faulting and earthquakes in the spinel stability field, Earth planet. Sci. Lett., 182, 171–185. Bruhns, O., 2003. Objective rates in finite elastoplasticity, in IUTAM Symposium on Computational Mechanics of Solid Materials at Large Strains, pp. 151–160, ed Miehe, C., Springer. Bruhns, O., Xiao, H. & Meyers, A., 1999. Self-consistent Eulerian rate type elasto-plasticity models based upon the logarithmic stress rate, Int. J. Plast., 15, 479–520. Bruhns, O., Meyers, A. & Xiao, H., 2001a. Analytical perturbation solution for large simple shear problems in elastic-perfect plasticity with the logarithmic stress rate, Acta Mech., 151, 31–45. Bruhns, O., Xiao, H. & Meyers, A., 2001b. Large simple shear and torsion problems in kinematic hardening elasto-plasticity with logarithmic rate, Int. J. Solids Struct., 38, 8701–8722. Chakrabarty, J., 2006. Theory of Plasticity, 3rd edn, Elsevier ButterworthHeinemann. Colak, U.O., 2004. Modeling of large simple shear using a viscoplastic overstress model and classical plasticity model with different objective stress rates, Acta Mech., 167, 171–187. Conrad, C.P. & Hager, B.H., 1999. Effects of plate bending and fault strength at subduction zones on plate dynamics, J. geophys. Res., 104, 17 551– 17 571. Coussy, O., 2004. Poromechanics, John Wiley & Sons. Dealy, J.M., 2010. Weissenberg and Deborah numbers—their definition and use, Rheol. Bull., 79, 14–18. Dienes, J.K., 1979. On the analysis of rotation and stress rate in deforming bodies, Acta Mech., 32, 217–232. Durban, D., 1990. A comparative study of simple shear at finite strains of elastoplastic solids, Q. J. Mech. Appl. Math., 43, 449–465. Freeman, S. & Weissenberg, K., 1948. Some new rheological phenomena, Nature, 162, 320–323. Hobbs, B.E., Ord, A. & Regenauer-Lieb, K., 2011. The thermodynamics of deformed metamorphic rocks: a review, J. Struct. Geol., 33, 758–818. Jaumann, G., 1911. Geschlossenes System physikalischer und chemischer Differential-Gesetze, Sitzber. Akad. Wiss. Wien, Abt. IIa, 120, 385–530. Johnson, S.E., 2009. Porphyroblast rotation and strain localization: debate settled!, Geology, 37, 663–666. Karato, S., 2008. Deformation of Earth Materials: An Introduction to the Rheology of Solid Earth, Cambridge Univ. Press. Downloaded from https://academic.oup.com/gji/article/211/1/252/3965336 by guest on 30 January 2024 energy in the examined W range (Appendix A3, Fig. A1). Because the loading at the boundaries of the homogeneously deforming unit square is never removed and material properties remain constant, elastic strain and stored energy should not decrease. Therefore, the MJ model violates the self-consistency requirement and also the principle of energy conservation. The associated shear-stress decrease becomes noticeable at W = 0.3 (Fig. 2a). Hence, the MJ model should not be used at W ≥ 0.3. The use of the MJ model is thus restricted to the W interval ]0.1; 0.3[. This rather small parameter range and the fundamental conceptual shortcomings of the MJ model lead to the question if there is any practical merit in using it. Our analysis suggests that the FT model should be the tool of choice for highly non-coaxial finite viscoelastic deformation problems in the intermediate- and high-W regimes. This notion holds true if one considers the concepts of frame indifference and selfconsistency (sensu Bruhns et al. 1999) as meaningful physical constraints. Nevertheless, even though the FT model is objective and energetically consistent, it remains to be seen if it can be falsified in laboratory experiments. Experimental tests of the model with rotary rheometry and analogue materials are subject of ongoing work. In addition to its unique shear-stress response, the FT model predicts non-zero normal stresses in simple shear, which can also be observed and tested experimentally (Freeman & Weissenberg 1948). 259 260 C.E. Schrank et al. Sone, H. & Uchide, T., 2016. Spatiotemporal evolution of a fault shear stress patch due to viscoelastic interseismic fault zone rheology, Tectonophysics, 684, 63–75. Tohver, E., Trindade, R.I.F., Solum, J.G., Hall, C.M., Riccomini, C. & Nogueira, A.C., 2010. Closing the Clymene ocean and bending a Brasiliano belt: evidence for the Cambrian formation of Gondwana, southeast Amazon craton, Geology, 38, 267–270. Truesdell, C., 1956. Hypo-elastic shear, J. Appl. Phys., 27, 441–447. Wheeler, J., 2010. Anisotropic rheology during grain boundary diffusion creep and its relation to grain rotation, grain boundary sliding and superplasticity, Phil. Mag., 90, 2841–2864. White, J., 1964. Dynamics of viscoelastic fluids, melt fracture, and the rheology of fiber spinning, J. Appl. Polym. Sci., 8, 2339–2357. White, J. & Mawer, C., 1992. Deep-crustal deformation textures along megathrusts from Newfoundland and Ontario: implications for microstructural preservation, strain rates, and strength of the lithosphere, Can. J. Earth Sci., 29, 328–337. Xiao, H., Bruhns, O.T. & Meyers, A., 1997a. Hypo-elasticity model based upon the logarithmic stress rate, J. Elast., 47, 51–68. Xiao, H., Bruhns, O.T. & Meyers, A., 1997b. Logarithmic strain, logarithmic spin and logarithmic rate, Acta Mech., 124, 89–105. Xiao, H., Bruhns, O.T. & Meyers, A., 1998a. On objective corotational rates and their defining spin tensors, Int. J. Solids Struct., 35, 4001–4014. Xiao, H., Bruhns, O.T. & Meyers, A., 1998b. Strain rates and material spins, J. Elast., 52, 1–41. Xiao, H., Bruhns, O. & Meyers, A., 2006. Elastoplasticity beyond small deformations, Acta Mech., 182, 31–111. Zhang, Y., Mancktelow, N.S., Hobbs, B.E., Ord, A. & Mühlhaus, H.B., 2000. Numerical modelling of single-layer folding: clarification of an issue regarding the possible effect of computer codes and the influence of initial irregularities, J. Struct. Geol., 22, 1511–1522. Ziegler, H., 1977. An Introduction to Thermomechanics, Vol. 21, NorthHolland Publishing Co. APPENDIX A1 Analytical solution of the MJ model The MJ model for simple shear requires the solution of the following nonhomogeneous linear system of ordinary differential equations: ηγ̇ = trel (σ̇12 − γ̇ σ̇22 ) + σ12 0 = trel (σ̇22 + γ̇ σ12 ) + σ22 . (A1) For convenience, we insert the placeholders: −1 A = −trel (A2a) B = γ̇ (A2b) −1 , C = ηγ̇ trel (A2c) which leads, after rearrangement, to σ̇12 = Aσ12 + Bσ22 + C (A3a) σ̇22 = Aσ22 − Bσ12 . (A3b) The solution is σ12 (t) = e At [sin (Bt) Q 1 − cos (Bt) Q 2 ] − AC A2 + B 2 (A4a) σ22 (t) = e At [sin (Bt) Q 2 + cos (Bt) Q 1 ] − BC . A2 + B 2 (A4b) Downloaded from https://academic.oup.com/gji/article/211/1/252/3965336 by guest on 30 January 2024 Karrech, A., Regenauer-Lieb, K. & Poulet, T., 2011. Frame indifferent elastoplasticity of frictional materials at finite strain, Int. J. Solids Struct., 48, 397–407. Karrech, A., Poulet, T. & Regenauer-Lieb, K., 2012. Poromechanics of saturated media based on the logarithmic finite strain, Mech. Mater., 51, 118–136. Kaus, B.J.P. & Becker, T.W., 2007. Effects of elasticity on the Rayleigh– Taylor instability: implications for large-scale geodynamics, Geophys. J. Int., 168, 843–862. Kaus, B.J.P. & Podladchikov, Y.Y., 2006. Initiation of localized shear zones in viscoelastoplastic rocks, J. geophys. Res., 111, B04412, doi:10.1029/2005JB003652. Korenaga, J., 2007. Effective thermal expansivity of Maxwellian oceanic lithosphere, Earth planet. Sci. Lett., 257, 343–349. Kusznir, N.J. & Bott, M.H.P., 1977. Stress concentration in the upper lithosphere caused by underlying visco-elastic creep, Tectonophysics, 43, 247– 256. Larsen, C.F., Motyka, R.J., Freymueller, J.T., Echelmeyer, K.A. & Ivins, E.R., 2005. Rapid viscoelastic uplift in southeast Alaska caused by postLittle Ice Age glacial retreat, Earth planet. Sci. Lett., 237, 548–560. Lehmann, T., 1972a. Anisotropische plastische Formänderung, Rom. J. Techn. Sci. Appl. Mech., 17, 1077–1086. Lehmann, T., 1972b. Einige Bemerkungen zu einer allgemeinen Klasse von Stoffgesetzen für große elasto-plastische Formänderungen, Ing.-Arch., 41, 297–310. Liu, X., Eckert, A. & Connolly, P., 2016. Stress evolution during 3D singlelayer visco-elastic buckle folding: implications for the initiation of fractures, Tectonophysics, 679, 140–155. Mancktelow, N.S., 1999. Finite-element modelling of single-layer folding in elasto-viscous materials: the effect of initial perturbation geometry, J. Struct. Geol., 21, 161–177. Mancktelow, N.S., 2002. Finite-element modelling of shear zone development in viscoelastic materials and its implications for localisation of partial melting, J. Struct. Geol., 24, 1045–1053. Mandal, N. & Khan, D., 1991. Rotation, offset and separation of obliquefracture (rhombic) boudins: theory and experiments under layer-normal compression, J. Struct. Geol., 13, 349–356. Meyers, A., Bruhns, O. & Xiao, H., 2005. Objective stress rates in repeated elastic deformation cycles, PAMM, 5, 249–250. Moresi, L., Dufour, F. & Mühlhaus, H.B., 2002. Mantle convection modeling with viscoelastic/brittle lithosphere: numerical methodology and plate tectonic modeling, Pure appl. Geophys., 159, 2335–2356. Moss, W.C., 1984. On instabilities in large deformation simple shear loading, Comput. Methods Appl. Mech. Eng., 46, 329–338. Mühlhaus, H.-B. & Regenauer-Lieb, K., 2005. Towards a self-consistent plate mantle model that includes elasticity: simple benchmarks and application to basic modes of convection, Geophys. J. Int., 163, 788–800. Passchier, C.W. & Simpson, C., 1986. Porphyroclast systems as kinematic indicators, J. Struct. Geol., 8, 831–843. Peltier, W.R., 1976. Glacial-Isostatic Adjustment—II. The Inverse Problem, Geophys. J. Int., 46, 669–705. Pfiffner, O.A. & Ramsay, J.G., 1982. Constraints on geological strain rates: Arguments from finite strain states of naturally deformed rocks, J. geophys. Res., 87, 311–321. Ramsay, J.G., 1980. Shear zone geometry: a review, J. Struct. Geol., 2, 83–99. Ramsay, J.G. & Graham, R.H., 1970. Strain variation in shear belts, Can. J. Earth Sci., 7, 786–813. Ramsay, J.G. & Huber, M.I., 1987. The Techniques of Modern Structural Geology, Vol. 2: Folds and Fractures, Vol. 2, Academic Press. Regenauer-Lieb, K. et al., 2013. Multiscale coupling and multiphysics approaches in earth sciences: theory, J. Coupled Syst. Multiscale Dyn., 1, 49–73. Regenauer-Lieb, K. et al., 2014. Multiscale coupling and multiphysics approaches in earth sciences: applications, J. Coupled Syst. Multiscale Dyn., 1, 281–323. Robin, C.M.I. & Bailey, R.C., 2009. Simultaneous generation of Archean crust and subcratonic roots by vertical tectonics, Geology, 37, 523–526. Viscoelasticity at large transformations From the initial conditions σ12 (t = 0) = σ22 (t = 0) = 0, the constants Q1 and Q2 are calculated: Q1 = BC A2 + B 2 Q2 = − (A5a) AC . A2 + B 2 (A5b) Re-inserting all substitutions yields the final result:   ηγ̇ − t −1 −1 σ12 (t) = −1 e trel γ̇ sin (γ̇ t) − trel cos (γ̇ t) + trel trel + trel γ̇ 2 (A6a) σ22 (t) =   − t −1 e trel γ̇ cos (γ̇ t) + trel sin (γ̇ t) − γ̇ . (A6b) To compute the total energy, we integrate eq. (7) after substituting eq. (A6a) and W: E TOT (t) =  GW 2 GW γ̇ − t  t − e trel W 2 − 1 cos (γ̇ t) 2 2 W2 + 1 (W + 1) + 2W sin (γ̇ t) + Q 3 . Therefore, the Jacobian matrix describing the incremental constitutive behaviour before co-rotational correction is i j η = t (δik δ jl + δil δ jk ), (A12) + trel h̃ kl 2 where h̃ kl is the Hencky strain tensor in Voigt notation, and δi j is the Kronecker delta. Moreover, the stress stored at the end of the increment is i j (t = t0 + t) = i j − ik k j + ik kj . (A13) The Jacobian matrix and the Kirchhoff stress at the end of the time increment expressed by eqs (A12) and (A13), respectively, are sufficient to describe the constitutive behaviour in a finite-element code. We implemented the above formulation into the user routine UMAT (Abaqus/Standard 2009) to solve the FT model numerically. (A7) The integration constant Q3 is again determined from the initial condition E TOT (t = 0) = 0:   GW 2 W 2 − 1 Q3 = . (A8) (W 2 + 1)2 A2 Numerical solution of FT model From a numerical point of view, i j is completely defined at the beginning of each time increment, given its purely kinematic nature. Similarly, to calculate the quantity − ik k j + ik k j , it is customary to use the Kirchhoff stress of the beginning of the time increment in a material integration module, the choice we opt for here. The only remaining quantity needed to obtain the co-rotational ⌣ ˙ i j . For the rate of the Kirchhoff stress  i j is the time derivative  Maxwell model, the time increment is ˙ i j = 2ηḣ i j , i j + trel  where gi j is some second-order tensor, gi j (t = t0 ) is its value at the beginning of the increment, and gi j is the change in gi j over the time increment t. Hence, it can be seen that  t i j = 2η h i j − ti j (t = t0 ). + trel (A11) 2 (A9) where h i j = 0.5 ln bi j is Hencky’s strain tensor, and bi j is the left Cauchy–Green strain tensor (for derivations of these tensors, see e.g. Bruhns et al. 1999; Karrech et al. 2011). This equation is the modified version of the constitutive model in eq. (1), involving the substitution of the Kirchhoff stress, insertion of the Maxwell relaxation time, and replacement of the deformation-rate tensor with the time derivative of the Eulerian logarithmic Hencky strain tensor. We use the central-difference scheme to approximate the stress and strain increments as  t gi j = ġi j t = t0 + 2 t  t gi j = gi j (t = t0 ) + , (A10) gi j t = t0 + 2 2 A3 Analysis of oscillations in the MJ model The analytical solution in eq. (15) explains the spurious softening and stress oscillations in the MJ model (Fig. 2). Let us consider the expression for the shear stress in the following. It contains the trigonometric term T12 = W sin(W k) − cos(W k), which is the obvious cause of the spurious oscillations. The term√T12 has a halfperiod of π W −1 . Its amplitude magnitude is |A| = W 2 + 1. The latter expression shows why T12 has no notable influence on the shear-stress response at W ≤ 0.1. In addition, T12 is multiplied with the damping factor e−k . Therefore, the oscillations are damped at sufficiently large k ≥ 2π (Fig. A1). In the range of ten relaxation times, models with W < 1 exhibit spurious softening, while those at larger W show oscillations (Fig. A1). The gravity of this fundamental problem becomes perhaps clearer when one considers the evolution of stored energy in the MJ model. Eq. (15a) implies immediately that the elastic energy also contains trigonometric terms and therefore is subject to oscillations. This can be seen in Fig. 3(c), which reveals unphysical negative stored energy. Due to the kinematics of the examined simple-shear case, the stored energy in the system must not decrease but approach a constant maximum value ∞ 2 −1 ) G in the infinite-time limit. In other words, of E elmax = 0.5(σvis because constant-velocity boundary conditions are applied to a homogeneous body with constant material properties, it should not relax stored stresses. However, the elastic-energy solution for the MJ model exhibits negative gradients. The critical dimensionless time kcrit , after which the time derivative of the elastic energy first becomes negative, is plotted for a number of W in the right panel of Fig. A1. Beyond kcrit the MJ model is dissipative and violates the first law of thermodynamics. Since the principal of energy conservation is at the heart of continuum mechanics, we propose to abandon the MJ model. Downloaded from https://academic.oup.com/gji/article/211/1/252/3965336 by guest on 30 January 2024 ηγ̇ −1 trel + trel γ̇ 2 261 262 C.E. Schrank et al. Downloaded from https://academic.oup.com/gji/article/211/1/252/3965336 by guest on 30 January 2024 Figure A1. Left panel: contour plot of e−k [W sin(W k) − cos(W k)] in W–k space. Right panel: plot of kcrit over W. The time kcrit is the first point in time at which ∂ E el /∂k = 0. For discussion, see the text.