The role of viscous regularization in dynamical problems, strain localization and mesh dependency
TThe role of viscous regularization in dynamical problems, strainlocalization and mesh dependency
Alexandros Stathas a , Ioannis Stefanou a, ∗ a Institut de Recherche en G´enie Civil et M´ecanique (UMR CNRS 6183), Ecole Centrale de Nantes, Nantes, France
Abstract
In this paper we revisit the elasto-viscoplastic, strain-softening, strain-rate hardening, model as a meansto avoid strain localization on a mathematical plane in the case of a Cauchy continuum. In contrastto the work of previous researchers de Borst and Duretz (2020); Needleman (1988); Sluys and de Borst(1992); Wang et al. (1997), we assume that both the frequency ω and the wave number k belong to thecomplex plane, therefore a different expression for the dispersion relation is derived. We prove then thatunder these conditions strain localization on a mathematical plane is possible. The above theoreticalresults are corroborated by numerical analyses, where the total strain and plastic strain profiles exhibitmesh dependent behavior. Keywords: strain localization, Perzyna, Consistency elasto-viscoplasticity, traveling waves, meshdependency, bifurcation, finite elements, Lyapunov stability
1. Introduction
Strain localization is a phenomenon which is found throughout natural and man-made structures.It ischaracterized by nonlinearity as well as different characteristic time and spatial scales ranging from thenear instantaneous fracture of brittle materials to the vast geological time required for the formation ofintricate patterns in earth strata.On a mathematical level, strain localization is understood as a bifurcation from the initial homoge-neous deformation state of the structure to another equilibrium path. This automatically raises ques-tions concerning the uniqueness of the reference homogeneous solution and its stability. The question ofuniqueness of solution is decided from the conditions needed such that the determinant of the acousticcharacteristic tensor of the problem is equal to zero Rice (1976); Rudnicki and Rice (1975). As far asthe evaluation of stability is concerned, different criteria with varying degrees of applicability have beenproposed in the literature (see Bigoni (2012); Chambon et al. (2004)).In this paper the question of stability of the obtained equilibrium paths is decided based on the Lya-punov analysis Chambon et al. (2004); Lemaitre et al. (2020); Mawhin (2005); Rice (1976); Stefanouand Alevizos (2016). In a Classical Cauchy continuum it has been proven that in the context of a strainsoftening plasticity yield criterion, localization happens on a mathematical plane. This renders the so-lution obtained from numerical methods such as the Finite Element method mesh dependent. However,experimental evidence in geomaterials suggests that localization in nature does not occur on a mathe-matical plane, rather it involves a small zone of finite thickness that accommodates the majority of thedeformation (see Chambon et al. (2004); Chester and Chester (1998); Muhlhaus and Vardoulakis (1988);Sibson (2003); Vardoulakis and Sulem (1995) among others). To remedy this inconsistency betweenexperiments and analytical and numerical predictions two approaches are often found in the literature.The first approach seeks to incorporate a modified constitutive law including the effects of viscosityde Borst and Duretz (2020); Sluys and de Borst (1992); Sluys et al. (1993); Wang et al. (1996a, 1997),while the other starts from the introduction of micromorphic continua, which introduce characteristiclength scales to the mathematical problem De Borst (1991); de Borst and Sluys (1991); Germain (1973);Rattez et al. (2018c); Sluys et al. (1993); Sulem et al. (2011); Vardoulakis (2009) . ∗ Corresponding author
Email address: [email protected] (Ioannis Stefanou)
Preprint submitted to Elsevier February 23, 2021 a r X i v : . [ phy s i c s . a pp - ph ] F e b oing beyond and expanding on existing results, we explore, in this paper, the role of viscosity forthe regularization of the classical Cauchy continuum. Considerable amount of research has been madeon viscous regularization of localization under both quasi-static and dynamic conditions (see Loret andPrevost (1990); Needleman (1988); Sluys and de Borst (1992)). In particular, for the quasi-static case,it is noted that the elasto-viscoplastic Cauchy medium based on a power law for the viscoplasticity de-scription, will exhibit strain localization on a mathematical plane except if a particular procedure forthe time integration takes place Needleman (1988). Furthermore, it is mentioned in Sluys and de Borst(1992) that the regularizing properties of the elasto-viscoplastic medium are present only in the contextof dynamical analyses, where the higher order inertial terms are introduced. Moreover, for the dynamiccase, a conclusion taken from Needleman (1988); Sluys and de Borst (1992) is that the selection of aconsistency or Perzyna yield condition for the viscoplastic model is prefered over a power law based onthe fact that the third order terms of the Partial Differential Equation (PDE), gradually vanish as strainsoftening occurs. As it has been already presented in Abellan and de Borst (2006) and in particular inde Borst and Duretz (2020) the formulation of the elasto-viscoplastic material model, specifically theposition of the damper element into the idealized viscoplastic configuration is of great importance. Ithas been shown in de Borst and Duretz (2020) that a solution localizes into a mathematical plane oncethe configuration of the viscosity dashpot, the plasticity element and the elastic spring are in series. Thepresent paper addresses stability and localization questions for the parallel configuration as studies arenot conclusive yet. Hence the elasto-viscoplastic model with strain softening, which involves a parallelconnection between the plasticity element and the dashpot, will be examined at present using bifurcationand Lyapunov stability analysis. In particular, we examine in detail the stability of the reference solutionof uniform deformation of the elasto-viscoplastic problem.In previous works de Borst and Duretz (2020); Needleman (1988); Sluys and de Borst (1992); Wanget al. (1997) the question of localization of the deformation was addressed looking at the propagationvelocity of localization. In particular it is mentioned in Sluys and de Borst (1992) that the original ill-posed problem of strain-softening plasticity presents imaginary wave speeds corresponding to standingwaves, which cannot extend the localization zone. This is thought to be remedied by the introductionof viscosity. It is also mentioned that to properly test the conditions under which strain localization ispresent in the non-linear elasto-viscoplastic problem, we would require a closed form analytical expres-sion for the solution, which until now is impossible. Since no analytical known way exists in solving thesoftening, rate-dependent plasticity problem, the focus was shifted in the derivation of the dispersionrelationships. Namely, in the previous, it was assumed that if every mode has a positive velocity, thenthe corresponding part of the deformation will propagate and therefore, it will not concentrate in onlyone element of the model as it would be the case with a standing wave. An additional argument thatis presented in Sluys and de Borst (1992); Wang et al. (1996a), is the dispersive character of the par-tial differential equation, which is assumed to further regularize the problem as the deformation frontwidens due to the different velocity of the deformation modes. In all these works, the dispersion relationwas derived based on the assumption that the circular frequency ω is real, while the wavenumber k iscomplex: k = k r + k i i, k r , k i = α ∈ R , indicating spatial attenuation of the derived deformation modes,thus leading to characteristic localization length l = α . However, we find no reason to strictly assume ω ∈ R , which makes an important difference in the analysis.In this paper we depart from this main assumption by assuming both ω, k ∈ C , thus considering theproblem in its general form. Furthermore, due to application of the Lyapunov analysis we are interestedin the magnitude of the imaginary part of ω, ω i = Im[ ω ], which controls the evolution of the amplitudeof the oscillations and therefore, the stability of the reference homogeneous deformation state. We arguethat if the amplitude of the mode with zero wavelength increases faster than that of the others, then thismode will dominate the deformation profile and strain localization on the mathematical plane will bepossible. Theoretical proofs are presented showing this phenomenon and also, the emergence of travelingwaves of strain localization in some cases. The theoretical results are corroborated by numerical analy-ses where the total strain and plastic strain profiles exhibit mesh dependent behavior. These analysesprovide counter examples to viscous regularization in dynamical problems.This paper is organized as follows: In section 2 the linearized differential equation describing the 1Delasto-viscoplastic model equiped with a von-Mises viscoplasticity consistency criterion is presented basedon the works of de Borst and Duretz (2020); Sluys and de Borst (1992); Wang et al. (1996b,a, 1997). InPage 2ection 3 the dispersion relation of the linearized equation is derived between the wavenumber k and thefrequency ω , both assumed to belong in the complex plane C . Assuming k to be complex is indicative ofwaves exhibiting attenuation or amplification as they travel through space, fairly common in a varietyof applications. The assumption of a complex ω is crucial in order for the Lyapunov exponent s = − iω to acquire negative or positive real part, which controls the stability of the reference homogeneous de-formation state. Applications exhibiting waves with complex ω can be found in Bernard et al. (2001);Deschamps et al. (1997); Gerasik and Stastna (2010); Mainardi (1984, 1987); Marion (2013). The intro-duction of the Lyapunov exponent s in the context of the dispersion analysis is a fundamental differenceof our work compared to previous attempts. The dispersion relation provides a link between the Lya-punov exponent that controls the evolution of the bifurcation and the wave number, which characterizesthe localization width. Thus we can specify which bifurcation evolves the fastest determining the widthof localization. We establish that the localization width λ tends to zero (localization on a mathematicalplane) for values of ω at the poles of the dispersion function k ( ω ). Furthermore, we present the conceptof localization for complex ω, k . In section 4 we present a series of numerical analyses that validate theabove theoretical findings. Finally conclusions are presented in section 5.
2. The elasto-viscoplastic wave equation
Let us consider a body under homogeneous small deformation lying at rest. The equilibrium equation isgiven as: σ (cid:63)ij,j = 0 , (1)where σ (cid:63)ij is the developed stress field, i, j are indices indicating the directions of the components of thestress field ( i, j = 1 ... , j ) indicates the partial spatial derivative with respectto the j direction. The Einstein summation condition over repeated indices is implied. Consideringa perturbation ˜ u i to the reference displacement field u (cid:63)i of homogeneous deformation, we find a rela-tionship between the perturbed stress and displacement ˜ u fields, according to the conservation of linearmomentum:( σ (cid:63)ij + ˜ σ ij ) ,j = ρ ¨˜ u i , (2)˜ σ ij,j = ρ ¨˜ u i , (3)where ( · ) corresponds to the time derivative and ρ is the density of the material. Both displacementand stress variations are arbitrary respecting only the boundary and loading conditions such that ˜ u i =0 , ˜ σ i,j n j = 0 at the boundary of the body, where displacement and loading conditions are specified,respectively. In order to continue with the bifurcation analysis of the problem we need to look first atthe elasto-viscoplastic constitutive law we take into account. As mentioned in section 1 a variety of yield criteria and flow rules are available for modeling viscoplas-ticity. Here the von Mises yiled criterion with strain-hardening (softening) and strain-rate hardening isused together with the consistency approach Wang et al. (1997).
In an elasto-viscoplastic formulation the following relations hold as given in De Borst (1991); Wanget al. (1996b,a, 1997): F ( σ ij , ¯ (cid:15) vp , ˙¯ (cid:15) vp ) = 0 , (4)˙ ε ij = ˙ ε eij + ˙ ε vpij , (5)˙ σ ij = M eijkl ( ˙ ε kl − ˙ ε vpkl ) , (6)˙ ε vpij = ˙ λ ∂F∂σ ij , (7)¯ (cid:15) vp = (cid:90) t ˙¯ (cid:15) vp dt, (8)Page 3here F = F ( σ ij , ¯ (cid:15) vp , ˙¯ (cid:15) vp ) is the yield function incorporating the effects of strain and strain-rate hard-ening through the use of the accumulated deviatoric viscoplastic strain ¯ (cid:15) vp and its rate ˙¯ (cid:15) vp , respectively.The viscoplastic multiplier ˙ λ is given by the consistency condition: ˙ F = 0 , ˙ λF = 0. The time derivativeof the yield condition in this case is the following:˙ F = ∂F∂σ ij ˙ σ ij + ∂F∂ ¯ (cid:15) vp ˙¯ (cid:15) vp + ∂F∂ ˙¯ (cid:15) vp ¨¯ (cid:15) vp = 0 . (9)The von Mises yield criterion with strain-hardening (softening) and strain-rate hardening for the consis-tency approach reads Wang et al. (1996b,a, 1997): F ( σ ij , ¯ ε ij , ˙¯ ε ij ) = (cid:113) J ( σ ij ) − ( F + h ¯ (cid:15) vp + g ˙¯ (cid:15) vp ) , (10)where F is the initial yield strength of the material, h is a parameter indicating strain hardening of thematerial ( h < g is a param-eter indicating hardening of the material with increasing plastic strain-rate ( g < g >
0) is examined here.Alternatively one can define the classic von Mises yield condition (without the term g ˙¯ (cid:15) ij ) and assumethat after the material reaches the yield limit the viscoplastic strain-rate is given as ˙ ε vpij = FηF ∂F∂σ ij .This model is known as the Perzyna model. The results are the same as far as monotonic loading isapplied and non-holonomic behaviour of the material is excluded, provided that g = ηF . In the caseof stress reversal and subsequent unloading, however, the results between the consistency approach andthe Perzyna model will be different due to the elasto-viscoplastic component that the Perzyna modelpredicts during unloading (see Heeres et al. (2002)).In what follows the principal results of the elasto-viscoplastic bifurcation analysis are presented. Afull derivation of the elasto-viscoplastic constitutive relations is presented in the Appendix A. Applyingequation (10) to equations (6) and (A.12) of the Appendix A, we arrive at the relationship for the stressrate ˜ σ ij :˜ σ ij = M eijkl (cid:18) ˜ ε kl − C − h + C ˜ ε kl + hb kl − h + C ˙˜ λ (cid:19) , (11)where in order to simplify the notation we have replaced accordingly: ∂F∂ ˙ λ = h,C = ∂F∂σ ij M eijkl ∂F∂σ kl ,b kl = ∂F∂ ˙ λ ∂F∂σ kl . (12) We proceed now in deriving the general linearized perturbed equation of equilibrium for the given materiallaw under monotonic loading. Inserting equation (11) into equation (3) and taking into account thespatial derivative of equation (7), we obtain: M eijkl (cid:18) ˜ ε kl,j − C − h + C ˜ ε kl,j + h − h + C ˙˜ ε kl,j − M e − ijkl ρ ...˜ u i (cid:19) = ρ ¨˜ u i , (13)This equation describes the spatio-temporal evolution of perturbations from the reference solution ofhomogeneous deformation in 3D. We constrain our analysis to the study of 1D problems, since they constitute the simplest case to studylocalization and the regularization effects coming from the above material law. In this way direct parallelscan be drawn between our work and the main bulk of literature on the subject de Borst and DuretzPage 42020); Needleman (1988); Sluys and de Borst (1992); Wang et al. (1996b). For the shearing of an1D layer we assume that the shearing is coaxial to the direction of x and that the body deforms inthe direction x . Therefore ˜ u i = [0 , ˜ u ] T = [0 , ˜ u ] T . Since we are in a state of 1D deformation, only thederivatives along the 1D axis, x , survive, therefore we set ∂ ˜ u ∂x = ∂ ˜ u∂x . Taking into account the appropriatematerial constant M eijkl = D e = G , we proceed in deriving the perturbed linear momentum equationfor the shearing of an 1D elasto-viscoplastic layer: G ¯ h ∂ ˜ u∂x − ∂ ˜ u∂t (3 + ¯ h ) Gv s + ¯ η vp G (cid:18) ∂ ˜ u∂t∂x − v s ∂ ˜ u∂t (cid:19) = 0 , (14)where v s = (cid:113) Gρ and ¯ η vp G = ηF = g .This coincides with the elasto-viscoplastic equation derived by de Borst and Duretz (2020); Sluys andde Borst (1992); Wang et al. (1996b, 1997). However, the equation derived above describes the evolutionof a perturbation from the initial homogeneous deformation state. It will not be used as a descriptionof the total behavior of the material as it neglects the material behavior in unloading and we are notinterested in the solution of the elasto-plastic problem but only at the stability of the homogeneousdeformation state, in order to draw conclusions about strain localization Lemaitre et al. (2020). Herewe note that the equation (14) has time independent coefficients (autonomous system, see Brauer andNohel (1969)).As mentioned in Wang et al. (1997), this equation contains two components, a classical elastoplasticwave equation plus the higher order rate-dependent terms. The nature of this differential equation isdefined by the higher order derivatives. It is also stated in Sluys and de Borst (1992) that, in the limit ofhigh viscosity ¯ η vp → ∞ , only the rate terms contribute, since they travel with the elastic wave velocity.In this case, the implied deformation pulse will travel with the corresponding elastic wave velocity aspredicted in Loret and Prevost (1990) and Needleman (1988). We consider ¯ u = uu c , ¯ t = tt c , ¯ x = xx c , where u c , t c , x c are the characteristic displacement, time andlength, respectively. Applying these definitions to equation (14) we obtain: (cid:18) x c v s t c ∂ ¯ u∂ ¯ t − ∂ ¯ u∂ ¯ x ∂ ¯ t (cid:19) ¯ η vp t c ¯ h + x c v s t c h ¯ h ∂ ¯ u∂ ¯ t − ∂ ¯ u∂ ¯ x = 0 . (15)Introducing the characteristic velocity v c = x c t c , the result is written as: (cid:18) v c v s ∂ ¯ u∂ ¯ t − ∂ ¯ u∂ ¯ x ∂ ¯ t (cid:19) ¯ η vp t c ¯ h + v c v s h ¯ h ∂ ¯ u∂ ¯ t − ∂ ¯ u∂ ¯ x = 0 . (16)The above equation is linear and has solutions of the form:¯ u (¯ x, ¯ t ) = A exp [ i (¯ k ¯ x − ¯ ω ¯ t )] , (17)where ¯ k, ¯ ω ∈ C and A ∈ R is a constant indicating the wave amplitude. Finally, inserting thenon-dimensional solution (17) into the normalized equation (16) we arrive at:¯ k v s (¯ h ¯ t c − i ¯ η vp ¯ ω ) − v c ¯ ω [(3 + ¯ h )¯ t c − i ¯ η vp ¯ ω ] = 0 . (18)It is worth emphasizing that the choice of assuming both ¯ ω, ¯ k ∈ C is not studied extensively in theliterature. However, examples of the notion can be found in Mainardi (1984, 1987); Marion (2013).Replacing V = v c v s , T = ¯ η vp t c and solving the above for ¯ k we obtain:¯ k , = ± (cid:115) V ¯ ω (3 + ¯ h − iT ¯ ω )(¯ h − iT ¯ ω ) . (19)The above equation can be also derived directly by use of a Fourier transform on equation (16). Expand-ing ¯ k, ¯ ω in imaginary and real parts, ¯ k = ¯ k r + ¯ k i i, ¯ ω = ¯ ω r + ¯ ω i i as explained in Hayes (1970); Maugin(2007); Poeverlein (1962) indicates that in our analysis the dependence of the amplitude of the solutionPage 5as both a spatial and a temporal component. In particular the non-dimensional solution can be writtenas:¯ u (¯ x, ¯ t ) = exp ( − ¯ k i ¯ x + ¯ ω i ¯ t ) exp[ i (¯ k r ¯ x − ¯ ω r ¯ t )] , (20)where, without loss of generality, the amplitude constant in front of the exponential terms of solution(20) is set to unity. The first factor in the right hand side of equation (20) indicates a quantity thatincreases or decreases based on the relationship between (¯ ω i ¯ t and ¯ k i ¯ x ). In this paper, we define thatan observer moving along ¯ x with a velocity c i , such that the amplitude profile exp ( − ¯ k i ¯ x + ¯ ω i ¯ t ) remainsconstant, is moving with the amplitude velocity : c i = ¯ ω i ¯ k i . (21)Conversely the second factor of equation (20) indicates the classical wave solution. An observer movingwith a velocity c r such that the phase exp [ i (¯ k r ¯ x − ¯ ω r ¯ t )] remains constant is said to be moving with the phase velocity Marion (2013); Pain and Beyer (1993); Sluys and de Borst (1992): c r = ¯ ω r ¯ k r . (22)According to the definition of Lyapunov for continuous dynamical systems Rattez et al. (2018a); Mawhin(2005), for an equilibrium solution to be unstable, the amplitude of the initial perturbation must increasein time. According to Lyapunov stability analysis, a partial solution of the partial differential equation(16) is given by ¯ u (¯ x, ¯ t ) = exp (¯ s ¯ t + i ¯ k ¯ x ), where ¯ s is the Lyapunov exponent ¯ s = − i ¯ ω . From thiswe conclude that the important term whose sign determines the stability of the reference solution,corresponding to homogeneous deformation (17), is the imaginary part of ¯ ω , i.e. ¯ ω i . Therefore, for theperturbation to grow in amplitude, the term exp ( − ¯ k i ¯ x + ¯ ω i ¯ t ) must be increasing as the wave travels.Localization on a mathematical plane will happen if we can find appropriate ¯ ω i , ¯ k i terms for the amplitudeto be constantly increasing the fastest for the smallest possible wavelength ¯ λ → k r = π ¯ λ → ∞ ).
3. Dispersion analysis
Equation (19) consists of two multivalued functions ¯ k (¯ ω ) , ¯ k (¯ ω ) in the complex set ¯ ω ∈ C . Introductionof branch cuts along selected points of unambiguous value is needed for their study on the values of theirargument ¯ ω Arfken and Weber (1999).Noticing the square powers of V , ¯ ω inside the root, equation (19) can be simplified yielding:¯ k , (¯ ω ) = ± V (cid:18) h ¯ h (cid:19) ¯ ω (cid:18) hT i + ¯ ω (cid:19) (cid:18) ¯ hT + ¯ ω (cid:19) − . (23)Each of the two solutions contains right and left propagating waves based on the sign combinations of¯ ω r , ¯ k r . The second solution ¯ k ( ω ) presents the exact same points of interest as the first. The differencelies in the − e iπ = − k (¯ ω ) , ¯ k (¯ ω ). The colors on the right of Figure 1 are changed indicating a change ofthe argument from the upper half of the imaginary plane to the lower half of it, meaning that the twosolutions ¯ k , (¯ ω ) behave differently when it comes to the spatial amplification/attenuation coefficient ¯ k i .In particular, ¯ k (¯ ω ) predicts only attenuation waves, while ¯ k (¯ ω ) consists of amplification waves. Thischange indicates that the waves in Figure 1 travel in the same direction with opposite imaginary part of¯ k (¯ ω ). To us the propagation direction of the wave is not important because of spatial symmetry of thesolution (19).In the next section the points of interest of the ¯ k (¯ ω ) are presented and their behavior is explainedin the form of branch cuts and poles. As discussed previously, the same behavior is valid for ¯ k (¯ ω ).We choose to draw further conclusions in the form of plots over line from the combination of the twosolutions ¯ k , (¯ ω ). In particular, we focus on the positive real parts of the solutions ¯ k r (¯ ω ) , ¯ k r (¯ ω ) > k i (¯ ω ) , ¯ k i (¯ ω ) >
0. Thus, we investigate the function¯ k ( ω ) = | ¯ k r (¯ ω ) | + i | ¯ k i (¯ ω ) | as shown later in Figure 4.Page 6 .2. Poles and zeros Studying equation (23), the following points can be readily specified in the above form: • The third factor indicates a zero at the origin: ¯ ω O1 = 0. • The fourth factor becomes zero at position: ¯ ω O2 = − hT i . • The last factor indicates the presence of a pole at: ¯ ω P1 = − ¯ hT i . • The value of the function at complex infinity ¯ ω P2 → ∞ is found to be infinite in a complex sense,lim ¯ ω → ¯ ω P2 ¯ k (¯ ω ) → ∞ .For the purposes of our analysis the behavior of the dispersion function at the poles ¯ ω P1 , ¯ ω P2 is veryimportant as it will be shown to promote localization on a mathematical plane. Figure 1: On the left: Complex 3D plot of the ¯ k (¯ ω ) solution (19). On the right: Complex 3D plot of the ¯ k (¯ ω ) solution(19). Values on vertical axis indicate the solution’s magnitude, where the coloring indicates the argument of the function.Along the branch cut discontinuity, the difference in color indicates the jump in the argument of ¯ k (¯ ω ). Because of the fractional powers of the second and third term, equation (23) is a multivalued equation,since it is affected by the values of the argument. In order to remove the ambiguity from the functionwe need to constrain it in such a way that each value of the function corresponds to only one argument.For this we introduce branch cuts . A branch cut is a discontinuity in the function that is defined byarbitrarily joining the two points defined as branch points. The branch points are defined as points ofunambiguous value, where the argument of the function is exactly known for a particular value of thefunction and the values corresponding to other points in a region sufficiently close to the branch pointdepends on the argument of the complex number inserted in the function. Two such points for thecomplex function f (¯ ω ) = ¯ ω are ¯ ω = 0 , and | ¯ ω | → ∞ , ¯ ω ∈ C , because for these particular numbersthe value of the function is always zero and infinity respectively. However, around them the value of thefunction depends on the argument of the complex number (see below).We can translate this result to other points in the complex plane, namely to ¯ ω O2 , ¯ ω P1 . In a regionclose and around ¯ ω O2 = − i hT the complex number with starting point ¯ ω O2 that follows the curvesurrounding ¯ ω O2 changes its argument by 2 πi . However, the factor (cid:16) hiT − ¯ ω (cid:17) only changes by πi ,meaning that there is a sign difference between the starting and the end position along the closed curveat the same point. Similarly, the same happens in the region near the pole ¯ ω P1 , where for every 2 πi that the relative complex number starting at ¯ ω P1 changes following the surrounding curve, the factor (cid:16) ¯ hiT − ¯ ω (cid:17) − changes by − πi . However at complex infinity (¯ ω → ∞ ), both previous points are entailedby the curve at infinity. Therefore, the total change in the argument is πi − πi = 0. This means that, ¯ ω P2 it is not a branch point (it does not belong to the branch cut). The simplest cut is the one that followsthe line defined by the two branch points as shown in Figure 2. Since the point at complex infinity isnot a branch point then it can be expected to be an isolated singular point, namely a pole of n-order oran essential singularity Arfken and Weber (1999); Brown et al. (2009). In this case it can be proven tobe a simple pole as shown in Appendix B.1Introducing the mapping ¯ ω = z we notice ¯ ω → ∞ can be written as z when z →
0. The proper-ties of this mapping are explained in Brown et al. (2009) and in the Appendix B.1.Page 7 igure 2: Contours of the solution (19) indicating with red color the position of the branch points and with cyan the branchcut line that connects them.Figure 3: Complex 3D plot of the ¯ k ( z ) where z = ω . On the left part of the Figure the behavior in the region near ¯ ω infinity ¯ ω P2 and the pole at ¯ ω P1 is presented.We notice the two poles lying at positions ¯ ω P2 : z = 0 and ¯ ω P1 lying at ω P1 : z p = − i respectively. On the right the region close to the pole at infinity z = 0 is plotted. We notice the existenceof the zero ¯ ω O2 that due to the mapping now lies at z = 0 .
05 extremely close to z = 0, inside the unit circle. The plots showing the poles and zeros of the function ¯ k ( z ) with the mapping are shown in Figure 3.On the right part of the Figure the region around the poles and infinity is shown while on the left a de-tail is presented where the zero value ¯ ω O2 - that due to the mapping is found closer to the origin - is shown. Localization will happen when the amplitude of a particular perturbation mode as shown in equation 20is found to be continuously increasing faster than the rest. If this happens for the perturbation of thesmallest possible wavelength ¯ λ → k r → ∞ then localiazation on a mathematicalplane takes place.As stated in the previous paragraphs without loss of generality we focus on the positive real and imag-inary parts of the function | ¯ k (¯ ω ) | = | ¯ k r | + | ¯ k i | i . | ¯ k (¯ ω ) | has the same poles ¯ ω P and zeros ¯ ω O as theoriginal ¯ k (¯ ω ) with the added simplification that only the positive argument values of both functions¯ k (¯ ω ) , ¯ k (¯ ω ) are plotted. This simplifies our analysis with regards to the sign of ¯ k i , which contributes toPage 8 igure 4: Complex 3D plot of the ˜ k (¯ ω ) combination of the two solutions. This envelope incorporates all the waves thattravel to the positive part of the axis together with the highest spatial amplification coefficient. On the left: Complex 3DPlot of ˜¯ k (¯ ω ) = | ¯ k r | + i | ¯ k i | for values of ¯ ω r , ¯ ω i close to the pole value ¯ ω P1 = 0 . i . On the right: Complex 3D Plot of˜¯ k (cid:0) z (cid:1) = | ¯ k r | + i | ¯ k i | for values of ¯ ω r , ¯ ω i close to infinity, when z P2 = 0 , (cid:0) z (cid:1) → ∞ and the pole value ¯ ω P1 = 0 . i → z P = − i non-dimensional M1 M2material parameters stable unstable¯ V = v c v s T = gt c G h = hG Table 1: Non-dimensional material parameters used for the numerical analyses the exponential growth of the amplitude, but it is not crucial for the time evolution of the perturbationwhich is determined by ¯ ω .First we focus to the pole value at ¯ ω P1 which is shown in the left 3D plot of Figure 4. There thepole ¯ ω P and the first zero ¯ ω O are shown. The pole lies at the value ¯ ω P1 = − ¯ h ¯ T i with ¯ h <
0, cor-responding to a real and positive Lyapunov exponent ¯ s = i ¯ ω = ¯ h ¯ T >
0. Since for ¯ k r (¯ ω ) → ∞ when¯ ω = ¯ ω P1 , we conclude that localization on a mathematical plane is possible (¯ λ = π ¯ k r → s ≈ λ ). However, in both cases strain localization happens on a mathematical plane and, therefore,this analysis shows that viscoplasticity does not regularize this problem even in the presence of inertiaterms.The behavior of | ¯ k (¯ ω ) | at ¯ ω → ∞ meaning ¯ ω = ¯ ω P cannot be captured in this plot. For this reason weperform a change of variables in ¯ ω , replacing with ¯ ω = z . Now the value of ¯ ω at infinity corresponds to z = 0. The function ¯ k ( z ) is plotted on the right of Figure 4, where we capture the value at ¯ ω P2 → ∞ when z = 0 (¯ ω P ) and at the pole ¯ ω P1 = z P1 .The second zero value is also captured in the line plots ofFigures 6 and 7 with the help of the transformation ¯ k ( z ). In this case ¯ k ( z ) tends to infinity. Therefore ω P also constitutes a localization point. Having qualitatively described the behavior of the solutions of the dispersion equation we proceed withassigning specific values to the dimensionless parameters according to the material parameters takenfrom de Borst and Duretz (2020). These values are presented in Table 1 and describe the case of aviscoplastic material obeying the von Mises yield criterion with strain-hardening ( h >
0) and strain-ratehardening ( g > h <
0) and strain-rate hardening ( g > ω i . The contours ofPage 9he real | ¯ k r | and imaginary parts | ¯ k i | of the combination of solutions near the pole ω P are presented inFigure 5. We can define three cases for an 1D elasto-viscoplastic medium expanding to infinity in bothdirections around the origin: Figure 5: Contourplot of ¯ k r , | ¯ k i | for values of ¯ ω r , ¯ ω i close to the pole value ¯ ω P = − . i . The contour closer to ¯ k → We focus our attention on the line where Re[¯ ω ] = ¯ ω r = 0 (see also Figure 5, Case 1). In this case standingwaves are present in the medium. The amplitude of these standing waves is dependent on the valuesof ¯ k i , ¯ ω i . When ¯ k i >
0, while ¯ ω i = 0, the amplitude of the standing wave decreases with the distancefrom the origin as shown in equation (20). However, if ¯ ω i >
0, then the value of the amplitude of theoscillations at fixed positions, will grow exponentially with time. Thus strain localization will happeninside a length ¯ λ = π ¯ k r .Along the imaginary axis (¯ ω r = 0), there is only one position where ¯ k r → ∞ . This is the pole at¯ ω P1 which constitutes a branch point of the dispersion relation (18) (see also Figure 6-asymptote). Inthis point ¯ k r (¯ ω P1 ) → ±∞ as well as ¯ k i (¯ ω P1 ) → ±∞ . We notice that the value of ¯ ω i at the pole is acutoff value. There are no higher values of ¯ ω i for which standing waves are possible since ¯ k r = 0 for¯ ω i > ¯ ω P1 (see also Figure 6). For a perturbation from the initial homogeneous state to grow with timewe are interested only in ¯ ω i >
0. Therefore the behavior of ¯ k r for ¯ ω i < λ = π ¯ k r → ω i strain localization on a mathematical plane is inevitable.Next, we investigate the influence of ¯ k i to the evolution of the amplitude of the perturbation. Lookingat Figure 7 and considering the solution branch of equation (19) we notice that ¯ k i → ∞ for ¯ ω = ¯ ω P1 .Therefore, a standing wave with solution parameters (¯ ω, ¯ k ) defined the same as the pole ¯ ω P1 will exhibitlocalization on a mathematical plane as distance ¯ x from the origin increases. Another important case seen in bibliography Abellan and de Borst (2006); de Borst and Duretz (2020);Sluys and de Borst (1992); Wang et al. (1996a, 1997) is that of the traveling waves where the imaginarypart of angular frequency is zero ¯ ω i = 0 (see Figure 5, Case 2). Therefore the Lyapunov exponentis also zero ( s = − i ¯ ω i = 0). In this case the amplitude growth is dependent only on ¯ k i : ¯ u (¯ x, ¯ t ) =exp( − ¯ k i ¯ x ) exp[(¯ k r ¯ x − ¯ ω r ¯ t )]. The parameter ¯ k i corresponds to the parameter α in de Borst and Duretz(2020); Sluys and de Borst (1992); Wang et al. (1996a, 1997) and its inverse l = α − is thought toconstitute a critical length that is supposed to regularize the problem, damping the waves of higherwavenumber ¯ k r and therefore avoiding strain localization on a mathematical plane. Here we show thatin fact depending on the solution branch of equation (19) the contribution of ¯ k i to the solution’s behaviorcan instead be positive, indicating amplification of perturbations of higher wavenumber ¯ k r . We focusour attention on Figure 8. We follow the red line corresponding to ¯ ω i = 0. We notice that the dispersionPage 10 igure 6: Left: Evolution of ¯ k r with respect to ¯ ω i for parameter values ¯ ω r close to the pole value ¯ ω P1 = 0 . i indicatingtraveling waves around the pole (¯ ω r (cid:54) = 0). Right: | k r | along the lines of ¯ ω r = const for a range of values of ¯ ω i closeto infinity.For ¯ z r = 0 the imaginary axis ¯ ω i is parallel to the imaginary axis z i . Therefore The detail around z i = 0 isindicative of the behavior of function ¯ k (¯ ω ) as ¯ ω r = 0 , ¯ ω i → ∞ .Figure 7: Evolution of ¯ k i with respect to ¯ ω i for parameter values ¯ ω r close to the pole value ¯ ω P1 = 0 . i indicating travelingwaves around the pole (¯ ω r (cid:54) = 0). Left: | ¯ k i (¯ ω ) | along the lines of ¯ ω r = const for a range of values of ¯ ω i , detail around thepole region ¯ ω P1 . Right: | ¯ k i ( z ) | along lines of constant z r . For z r = 0 the imaginary axis ¯ ω i is parallel to the imaginaryaxis z i . Therefore The detail around z i = 0 is indicative of the behavior of function ¯ k (¯ ω ) as ¯ ω r = 0 , ¯ ω i → ∞ . relation predicts ¯ k r = 0 for (¯ ω r , ¯ ω i ) = (0 , ω r we notice a quasi-linear increase of the wavenumber ¯ k r . Figure 8 shows that | ¯ k r | increases monotonicallyand tends to infinity for | ¯ ω r | → ∞ . The latter can be proven mathematically (see Appendix B.1). Fromthis we establish that perturbations whose wavelength tends to zero ¯ λ → k i . - - -
10 10 20 30102030 - - Figure 8: Dispersion curves (¯ ω r , ¯ k r ) for different values of parameter ¯ ω i along the line of zero temporal coefficient ¯ ω i = 0(case 2). With red color and the value passing from the pole ¯ ω P1 = 0 . i purple color. Detail of the dispersion for lowvalues of ¯ ω r is shown on left. Figure 9 shows that for ¯ ω r tending to infinity, the value of | ¯ k i | increases monotonically and tends to aceiling value ¯ k i → c ∈ R . The latter can be proven mathematically (see Appendix B.1). Therefore,since when ¯ k r ( ω ) → ∞ , ¯ k i < λ is increasing the fastest as the perturbation travelsthrough the medium. Therefore, strain localization on a traveling mathematical plane will happen. Based on Figure 5 and the diagrams of Figures 6, 7, 8 and 9 very general cases of traveling waves canbe examined. In Figures 6 and 7 we examine the evolution of ¯ k r , ¯ k i with respect to the imaginary partof the frequency ¯ ω i , by considering the real part of the angular frequency ¯ ω r as a parameter.Page 11e already presented case 1 of standing waves ¯ ω r = 0 where the influence of the pole leads to strainlocalization due to ¯ ω i > k r (¯ ω ) → ∞ as seen on the left diagram of Figure 6. For the rest ofthe values of the parameter ¯ ω r , the wavenumber ¯ k r is bounded. Therefore, no strain localization on amathematical plane will take place in these cases.In Figures 8 and 9 temporal amplification ¯ ω i is introduced as a parameter, keeping ¯ ω r as the inde-pendent variable. Figure 8 is indicative of the dispersion relation of the medium. Away from ¯ ω r = 0the dispersion relation ¯ ω r , ¯ k r becomes linear for all values of ¯ ω i and the resulting traveling waves have acommon phase velocity. We notice here that as ¯ ω r → ∞ , ¯ k r (¯ ω r , ¯ ω i ) → ∞ . For the waves with the samevalue for the parameter ¯ ω i > ω r the spatial amplification coefficient | ¯ k i | presents a ceiling value. Thisresult is already proven for ¯ ω i = 0 (see Appendix B.1). The ceiling value depends on the parametervalue ω i , namely it increases as the parameter ¯ ω i increases. The ceiling value of the amplification coef-ficient ¯ k i corresponds to a wavenumber ¯ k r that tends to infinity ¯ k r → ∞ . This result shows that in thegeneral case of traveling waves with infinitesimal wavelength ¯ λ = π ¯ k r →
0, strain localization on a travel-ing mathematical plane will happen due to the combination of ¯ k i < ω i >
0, provided that ¯ k r → ∞ .In the above we focused mainly on limit cases related to strain localization in a elasto-viscoplasticstrain-softening ( h < g > ω P2 Based on the behavior close to infinity, ¯ ω → ∞ or z →
0, we get ¯ k (¯ ω ) = ¯ k ( z ) → ∞ . In contrast to real ∞ that can be either positive or negative or indeterminate based on whether we approach the value of z from above or below zero, complex ∞ is indeterminate as at the pole value ¯ ω P the limit along eachdirection surrounding the pole indicates differences in the real and imaginary parts of ¯ k (¯ ω ). Since arounda simple pole like ¯ ω P2 , where ¯ ω P2 r → ∞ , ¯ ω P2 i → ∞ the argument of a complex function changes by afull 2 π radians, the two limiting cases for the value of ¯ k (¯ ω P2 ) are Re[¯ k ] = ¯ k r → ∞ , Im[¯ k ] = ¯ k r → k ] = ¯ k r → , Im[¯ k ] = ¯ k r → ∞ . Since ¯ k r → ∞ when ¯ ω = ¯ ω P2 , we conclude again, that localizationon a mathematical plane is possible. Since, localization on a mathematical plane happens for values of¯ ω → ∞ , the rate of increase of the perturbation amplitude as given by the Lyapunov exponent ¯ s = − i ¯ ω i is unbounded. ¯ ω P By expanding the solution space allowing for complex ¯ ω i we allow a connection with the Lyapunov ex-ponent ¯ s used in stability analyses. The new solution space is richer regarding the perturbations we canintroduce in the visco-elastoplastic medium. Some key characteristics retained by the solution from thedefinitions already found in the literature, is the exclusion of standing waves of infinitesimal length thatgrow with an infinite Lyapunov coefficient as in the case of the pure ill-posed rate-independent plasticityproblem. - -
50 50 10078910111213 - - Figure 9: Evolution of ¯ k i for different values of temporal coefficient ¯ ω i including the case zero temporal attenua-tion/amplification ¯ ω i = 0 to the pole value ¯ ω P1 = 0 . i . For traveling waves around the pole (¯ ω r (cid:54) = 0) the spatial attenuationcoefficient ¯ k i is reaching an upper bound. Subfigure on left presents the curve of | k i | along the lines of ¯ ω i = const for arange of values of ¯ ω r while subfigure on right presents a detail around the pole region ¯ ω P1 . Page 12he introduction of the parameter T = ¯ η vp t c allows for the existence of a zero value on the imagi-nary axis. This zero in turn plays the role of the branch point, forcing the argument of the pole ¯ ω P atinfinity to turn by π/ k (¯ ω ) along the imaginary axis. If that zero wasnot there, then the point at infinity would be a branch point, therefore strain localization on a math-ematical plane would happen for infinite ¯ ω i as in the case of a strain-softening rate-independent material.The visco-elastoplastic medium discussed here under the expansion of its dispersion equation solutionnegates the instantaneous localization of deformation, however as dictated by the pole ¯ ω P amplificationof the infinitesimal wavelength perturbation is still possible. In other words, the value of the Lyapunovexponent ¯ s = − i ¯ ω i becomes bounded due to viscoplasticity but this is not true for the value of thewavenumber ¯ k r . For case M1 we note that ¯ h = 0 . >
0. In this case the points of interest of the dispersion relation(23) change leading to different behavior than the one previously presented. The relationships for thedetermination of zeros and poles of the function remain the same (see section 3.2).For these numerical values, the pole ¯ ω P1 is reflected due to the change of sign of ¯ h around the Real axis.This however is not true for the he zero ¯ ω O1 since the sign of ( − h ) does not change (provided that | ¯ h | < ω O1 . Sincenow the pole lies on ¯ ω i < λ → , ¯ k r → ∞ are attenuated with time, therefore strain localization on a mathematical plane is not possible in thiscase.
4. Numerical analysis
In this section, two sets of numerical analyses are performed in order to verify the stability and possiblestrain localization of the elasto-viscoplastic strain softening ( h < g > )0 waveequation. Two sets of analyses were performed. In the first set, a 1D model of an elasto-viscoplasticinfinite string was used on which, the linear differential equation of third order (equation (16)) is numer-ically solved for a variety of initial conditions. We comment on the displacement profiles that developwith time and we verify the theoretical findings of section 3. In the second set, we proceed in numericallysolving the fully non-linear problem. The difference between these two sets of analyses, lies in the factthat in the second case unloading is allowed to take place. Therefore, we can investigate its influence onthe strain localization profiles. The non linear numerical analyses show also strain localization and meshdependency.
To model the infinite string a large length and the Sommerfield open boundary conditions were used. Thelatter can be used since the partial differential equation in question is linear and by use of the Fouriertransform can be shown to have partial solutions in the form of A exp( i (¯ ω ¯ t − ¯ k ¯ x )). Three modes ofinducing the perturbation from the reference homogeneous state were examined making use of non zeroinitial conditions for the string displacement. This is achieved by varying the shape of the perturbationwith three different ways: an initial pinch of the string at the middle, a cosinus pulse centered in themiddle as well as a Gaussian distribution centered at the middle.Two sets of parameters were used for the analyses, where the sign of the hardening parameter ¯ h varies be-tween positive or negative in order to compare between strain-hardening and a strain-softening materialas shown in Table 1. The material parameters used in the unstable case are those provided by de Borstand Duretz (2020). The numerical analyses were performed using the method of Finite Differences. Inparticular a central difference scheme was used for the spatial discretization of the PDE problem, thedomain of length L = 15m is discretized into 250 segments, resulting in a coupled system of ODE’s whichwas solved by the algorithm IDA of the Mathematica T M software package Inc. (2020).
We present in Figure 10 the behavior of the string after an initial pinching -application of initial displace-ment conditions at the middle node of the discretized domain. Mesh convergence analysis has shownthat, in the unstable case M2, the localization width is equal to the mesh size. Therefore, the elasto-viscoplastic formulation does not regularize the underlying problem as presented in the introduction andPage 13he solution is mesh dependent. This is in accordance with the theoretical results presented in subsec-tions 3.4.1, 3.4.2, 3.4.3. We note also that the strain hardening material of the M1 case does not lead tostrain localization as expected (see section 3.4.6). x - u t = = = = = =
5s 2 4 6 8 10 12 14 x - u t = = = = = = Figure 10: Evolution of the pinching perturbation at different times. Left: strain hardening material case M1. On theright: strain softening material case M2.
The behavior of the string after application of cosine initial conditions for the two sets of parameters (seeTable 1) is presented in Figure 11. Again we notice a localization of deformation for the case of strainsoftening. x u t = = = = = = x - u t = = = = = = Figure 11: Evolution of the cosine perturbation at different times. Left: strain hardening material case M1. On the right:strain softening material case M2.
In order to verify the theoretical prediction, that the elasto-viscoplastic medium with strain softeninglocalizes on a mathematical plane under dynamic loading conditions, we superpose three different cosineperturbations in the medium by varying the width of each perturbation as shown in Figure 12. Theperturbation wavelengths are ¯ λ : 1, 5, 10, corresponding to perturbation widths of 0.5 ,2.5 , 5. In Figure12, the position of the perturbations from left to right is ¯ x = 3 . → ¯ λ = 10 , ¯ x = 7 . → ¯ λ = 1 and ¯ x =11 . → ¯ λ = 5 respectively. Figure 12 shows that localization is accumulating faster for the smallestperturbation length, verifying the theoretical findings of section 3. x u t = = = = = = Figure 12: Evolution of the different wavelength cosine perturbations. Faster strain localization is observed for the smallestwavelength.
Page 14 .1.3. Centered Gaussian initial condition
In Figure 13 we present the behavior of the string after a gaussian perturbation of the initial conditions-application of initial gaussian displacement conditions centered at the middle node of the discretizeddomain. In the strain softening case M2, we notice that the localization of the deformation is containedinto a narrow band of finite length, dependent on the width of the initial perturbation as shown on theright of Figure 13.We emphasize, that the lower bound of localization is a result of the mesh discretization. Furtherincrease of the mesh will lead to narrower bands as expected by theory mesh dependency. In Figure14 we compare among the perturbation profiles at different times for Gaussian perturbations of varyingwidth. As in the previous case, the narrowest perturbation localizes the fastest. x - u t = = = = = =
5s 2 4 6 8 10 12 14 x - u t = = = = = = Figure 13: Evolution of the Gaussian perturbation at different times. Left: strain hardening material case M1. On theright: strain softening material case M2. x u t = = = = = = Figure 14: Evolution of different Gaussian perturbations varying in their width. Faster strain localization is observed forthe smallest wavelength.
In order to model the effect of viscous reqularization on the strain localization and possible mesh depen-dence in fully nonlinear cases without neglecting the effect of unloading, a series of dynamic localizationanalyses was performed using the finite element analysis program ABAQUS Smith (2009). The fullyimplicit Newmark scheme was used. The parameters used for the Newmark scheme correspond to thetrapezoidal rule ( α = 0 , β = , γ = ) in order to avoid numerical damping.The use of a User Material Subroutine (UMAT) was favored in order to incorporate the Perzyna elasto-viscoplastic constitutive material law into ABAQUS. The material parameters leading to mesh indepen-dent solution were taken from de Borst and Duretz (2020), see Table 2. These values are quite low forreal physical applications but they are used following de Borst and Duretz (2020) in order to allow directcomparisons. Configuration D1 corresponds to a set of parameters leading to reqularization of strainlocalization and therefore to mesh independent results, while configuration D2 corresponds to a set ofparameters leading to strain localization and therefore to mesh dependent results. Care was taken toremove additional viscosity from the analysis except for the one strictly prescribed by the material. Theanalyses were performed using 2D, solid reduced integration elements CPE4R.Page 15nalyses D1 D2mesh independence Yes No Units τ
12 12 Pa ρ kgm G H -200 -1000 Pa c
20 20 Pa η vp = ηF = g
50 25 s
Table 2: Material parameters used for the ABAQUS numerical analyses, the first set of parameters is taken from de Borstand Duretz (2020). analyses ∆ t element number [s]25 0.00150 0.001100 0.001200 0.001400 0.0002800 0.0002 Table 3: Time increment ∆ t measured in seconds (s), used for the ABAQUS numerical analyses, such that the CFL yieldcriterion is satisfied for each mesh discretization. We study the pure shear of a 1D layer (see Figure 15). To avoid bending we block displacementsalong the length of the model. A shear traction τ is applied instantaneously on top of the model andit propagates towards the fixed end at the base of the layer. We extract our results after the pulse hasreturned to the free end of the model. For the duration of the analysis the time increment is kept smallerthan ∆ t = 0 .
001 s, which is smaller than the time ∆
CF L needed for the elastic wave to traverse thesmallest element dimension of the mesh as specified by the Courant-Friedrichs-Lewy (CFL) criterion.Considering the length of the model L = 1 m, and the material parameters of Table 2 leading to theelastic shear wave velocity c g = (cid:113) Gρ = 4 ms , for a mesh discretization of 200 elements, we obtain:∆ t CFL = ∗ = 0 . Figure 15: 2D model of a layer subjected to shear.
To investigate whether an analysis is mesh dependent or not we plot the profiles over the length of themodel of the engineering total shear strain γ and shear plastic strain rate ˙ λ for different number ofelements. Results are given at the end of the analysis for time t = 0 . γ and the plastic multiplier ˙ λ profiles are spreadover the model length and converge upon mesh refinement. From the analysis of de Borst and Duretz(2020) we can extract the material length scale of the problem which is equal to l = η vp c √ ρG = 0 . F = τ − c − hγ p , ˙ γ p = ˙ λ = Fη vp c respectively). This is in close agreement with the total strain and plastic strain profiles we present inFigures 16 and 17. It should be noted, however, that the above relation does not take into account thesoftening slope of the material. To this end the relation available in Wang et al. (1996b) can be used aswell. x [m]0.0010.0020.0030.004 Element number:
Element number: x [m]0.000.010.020.030.04 Element number:
Element number:
Figure 16: On the left: Total shear strain γ profiles at the end of the analysis for different mesh discretization of themodel for the material set D1. The response converges to a mesh independent solution as the number of elements increases.On the right: Detail near the base of the model of the total shear strain γ profiles at the end of the analysis for differentmesh discretization of the model for the material set D2. . The response does not converge to a mesh independent solutionas the number of elements increases. However, for the configuration D2 the response of the model is completely different. From the profiles of γ , ˙ λ in Figures 16, 17 we notice that strain, plastic strain-rate localize on one element, and the solutiondoes not show any signs of converging upon mesh refinement. The material length calculated from therelation above is 0.2m . This distance however is much larger than the variation of the strain, plasticstrain-rate profiles shown in Figure 17. The profiles show narrower localization for finer discretizations. x [m]0.0000.0020.0040.0060.0080.010 [ ] Element number:
Element number: x [m]0.00000.00020.00040.00060.00080.0010 Element number: Element number:
Figure 17: On the left: Plastic multiplier ˙ λ profiles at the end of the analysis for different mesh discretization of the modelfor the material parameters in set D1. The response is mesh independent. The noise near the end of the plastified region isdue to the automatic incrementation scheme used in ABAQUS. The results agree well with the material length defined byde Borst and Duretz (2020). On the right: Detail near the base of the model of the plastic multiplier profiles ˙ λ at the endof the analysis for different mesh discretization of the model, for the material parameters in set D2. The response localizesto a mesh dependent solution as the number of elements increases. In both plots a detail near the region of interest 0 . In order to further point out that the analyses diverge upon mesh refinement for the second set of mate-rial parameters D2, we present in Figure 18 the total strain profile for finer discretizations of the mesh.The results are taken at time t = 0 .
36 s after the pulse of initial stress 18 Pa, is reflected and has passedthe middle of the bar for a second time. Again the time increment of the analyses is kept smaller thanthe one specified by the CFL criterion. We notice that upon mesh refinement the profiles of total strain γ become narrower (strain localization in narrower bands), while the peak value at the base of themodel is constantly increasing without signs of convergence.Therefore, according to the above results we can conclude that the analyses with the material set pa-Page 17ameters D2 constitute a counterexample, about the beneficial role of viscous regularization in strainlocalization and mesh dependency. x [m]0.000.010.020.030.040.05 Element number:
Element number:
Figure 18: Total shear strain γ profiles at at time t = 0 .
36 s for different mesh discretization of the model for appliedload τ = 18 Pa. The response localizes to a mesh dependent solution as the number of elements increases. A detail nearthe region of interest 0 .
1m is plotted.
5. Conclusions
In this paper we investigated the regularization properties of elasto-viscoplasticity regarding strain lo-calization and mesh dependency under the presence of inertia. Even though for quasi-static cases it iswell-known that elasto-viscoplasticity of Perzyna or consistency type do not regularize strain localizationNeedleman (1988); Sluys and de Borst (1992) in the dynamical case the situation is not clear.Our approach is theoretical and numerical. After deriving the equilibrium equation of the model understrain-softening strain-rate hardening elasto-viscoplasticity, we study the Lyapunov stability of states ofuniform/homogeneous deformation. In order to avoid unnecessary complexity we focus on a 1D shearingexample. Our mathematical analysis differs from previous ones de Borst and Duretz (2020); Wang et al.(1997); Sluys and de Borst (1992); Needleman (1988) by considering the frequency ¯ ω to be a complexnumber in addition to the complex wavenumber ¯ k . This is an important point for investigating thestability of the homogeneous, reference state.Next, we proceed in finding the dispersion relationship between the complex wave number ¯ k and thecomplex frequency ¯ ω for an arbitrary perturbation. The dispersion equation presents a pole who isresponsible for strain localization on a mathematical plane and thus for mesh dependency. More specif-ically, the wavenumber becomes infinite, ¯ k r → ∞ , and the wavelength, ¯ λ →
0, for strain-softening,strain-rate hardening, which means that localization on a mathematical plane is indeed possible in thissystem (see section 3.4). In the same section we have also made an extensive discussion about thepossibility of traveling waves in the medium 3.4.2, 3.4.3 and their relation to strain localization andwave attenuation. Several qualitative observations of the poles of the dispersion equation are providedin sections 3.4.4, 3.4.5. The analysis is completed by some additional observations about the behaviorof propagating sinusoidal monochromatic pulses and their relation with strain localization and phase c r and amplitude c i velocities (see Appendix C).We validate our theoretical findings with 1D numerical analyses of an infinite string (see section 4.1).In particular, we investigate the effects of the perturbation mode on the localization mode. The per-turbation is introduced in different shapes via various initial conditions. The theoretical relationshipbetween the width of the perturbation and its rate of increase is confirmed. We also confirm that thesmallest perturbations propagate the fastest leading to strain localization and mesh dependency. Basedon these results we conclude that the elasto-viscoplastic model with strain softening is unable to restrictthe classical Cauchy continuum from localizing on a mathematical plane.However, our analysis up to this stage, is based on a linearized version of the problem. In this sense,we cannot claim an absolute proof that elasto-viscoplasticity and inertia prevent strain localization andmesh dependency. For this purpose, we perform fully nonlinear dynamic numerical analyses using thePage 18BAQUS commercial Finite Element software Smith (2009) with a strain-softening, strain-rate hard-ening, Perzyna elasto-viscoplastic user material (UMAT). An impicit Newmark scheme was employed.Special attention was given to avoid any artificial numerical damping in order to guarantee that the rightpartial differential equations are solved. The results show that the behavior strongly depends on the ma-terial parameters used for the analyses and that mesh dependent solutions are possible. These analysesconfirm the theoretical findings and provide a counter-example showing that viscosity and inertia do notregularize strain localization and mesh dependency. Acknowledgments
The authors would like to thank Professor Nicolas Mo¨es for his comments and suggestions.The two authors would like to acknowledge the support of the European Research Council (ERC) underthe European Union’s Horizon 2020 research and innovation program (Grant agreement no. 757848CoQuake).
ApendicesA. Derivation of the Elasto-viscoplastic wave equation.
In this Appendix the results presented in section 2 are obtained in detail.
A.1. Elasto-viscoplastic constitutive relationsA.1.1. Consistency model.
In an elasto-viscoplastic formulation the following relations hold: F ( σ ij , ¯ (cid:15) vp , ˙¯ (cid:15) vp ) = 0 , (A.1)˙ ε ij = ˙ ε eij + ˙ ε vpij , (A.2)˙ σ ij = M eijkl ( ˙ ε kl − ˙ ε vpkl ) , (A.3)˙ ε vpij = ˙ λ ∂F∂σ ij , (A.4)¯ (cid:15) vp = (cid:90) t ˙¯ (cid:15) vp dt. (A.5)The viscoplastic multiplier ˙ λ is given by the consistency condition˙ F = 0 , ˙ λF = 0 , (A.6)where F ( σ ij , ¯ (cid:15) vp , ˙¯ (cid:15) vp ) is the yield function incorporating the effects of strain and strain-rate hardeningthrough the use of the accumulated viscoplastic strain ¯ (cid:15) vp and its rate ˙¯ (cid:15) vp respectively. The timederivative of the yield condition in this case is given as:˙ F = ∂F∂σ ij ˙ σ ij + ∂F∂ ¯ (cid:15) vp ˙¯ (cid:15) vp + ∂F∂ ˙¯ (cid:15) vp ¨¯ (cid:15) vp = 0 , (A.7)starting from the von Mises yield criterion: F ( σ ij , ¯ ε ij , ˙¯ ε ij ) = (cid:113) J ( σ ij ) − F − h ¯ (cid:15) vp − g ˙¯ (cid:15) vp , (A.8)assuming dependence of yield on the deviatoric invariant of the stress tensor J ( s ij ) = (cid:113) s ij s ij , and˙¯ (cid:15) vp = (cid:113) ˙ ε vpij ˙ ε vpij where s ij = σ ij − σ ii , we obtain that ˙ λ = ˙¯ (cid:15) vpij therefore the yield criterion as well as theconsistency condition can be written as: F ( σ ij , λ, ˙ λ ) = 0 , (A.9)˙ F = ∂F∂σ ij ˙ σ ij + ∂F∂λ ˙ λ + ∂F∂ ˙ λ ¨ λ = 0 , (A.10)Page 19ultiplying (A.3) by ∂F∂σ ij . replacing ˙ ε vpij with help from (A.4) and replacing the term in the lefthandsidewith equation (A.10) we finally get: − ∂F∂λ ˙ λ − ∂F∂ ˙ λ ¨ λ = ∂F∂σ ij M eijkl (cid:18) ˙ ε kl − ˙ λ ∂F∂σ kl (cid:19) . (A.11)Grouping together the terms of ˙ λ and solving for ˙ λ we get:˙ λ = ∂F∂σ ij M eijkl − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˙ ε kl + ∂F∂ ˙ λ − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ¨ λ, (A.12)Inserting (A.12) into (A.3) we obtain:˙ σ ij = M eijkl (cid:32) ˙ ε kl − ∂F∂σ ij M eijkl ∂F∂σ kl − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˙ ε kl + ∂F∂ ˙ λ ∂F∂σ kl − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ¨ λ (cid:33) , (A.13)Replacing the time derivative with a variation taking advantage of the definition of variation we arriveat the constitutive equation describing the perturbed field of stress ˜ σ ij .˜ σ ij = M eijkl (cid:32) ˜ ε kl − ∂F∂σ ij M eijkl ∂F∂σ kl − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˜ ε kl + ∂F∂ ˙ λ ∂F∂σ kl − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˙˜ λ (cid:33) . (A.14) A.2. Derivation of the perturbed equation
Inserting (A.14) into (3) we arrive at: M eijkl (cid:32) ˜ ε kl,j − ∂F∂σ ij M eijkl ∂F∂σ kl − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˜ ε kl,j + ∂F∂ ˙ λ ∂F∂σ kl − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˙˜ λ ,j (cid:33) = ρ ¨˜ u i . (A.15)Substituting (A.4) into the last term of the right hand side of eq. (A.15) such that ˙ ε vpkl = ∂F∂σ kl ˙ λ we arriveat: M eijkl (cid:32) ˜ ε kl,j − ∂F∂σ ij M eijkl ∂F∂σ kl − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˜ ε kl,j + ∂F∂ ˙ λ − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˙˜ ε vpkl,j (cid:33) = ρ ¨˜ u i , (A.16)rewritting ˙˜ ε vpkl = ˙˜ ε kl − M e − ijkl ˙˜ σ ij : M eijkl (cid:32) ˜ ε kl,j − ∂F∂σ ij M eijkl ∂F∂σ kl − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˜ ε kl,j + ∂F∂ ˙ λ − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˙˜ ε kl,j − M e − ijkl ˙˜ σ ij,j (cid:33) = ρ ¨˜ u i , (A.17)inserting finally (3) we obtain: M eijkl (cid:32) ˜ ε kl,j − ∂F∂σ ij M eijkl ∂F∂σ kl − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˜ ε kl,j + ∂F∂ ˙ λ − ∂F∂λ + ∂F∂σ ij M eijkl ∂F∂σ kl ˙˜ ε kl,j − M e − ijkl ρ ...˜ u i (cid:33) = ρ ¨˜ u i . (A.18) A.2.1. 1D Example: shearing of a viscous Cauchy layer
We proceed in deriving the perturbed linear momentum equation for the shearing of a 1D elasto-visoplastic Cauchy layer. The basic kinematic equations and the constitutive relations with the vonMises yield criterion with linear strain and strain-rate hardening/softening, are given as follows: F ( τ, λ ) = √ τ − ¯ τ ( λ ) − ηF ˙ λ, (A.19)˙ γ = ˙ γ el + ˙ γ vp , (A.20)˙ γ vp = ˙ λ ∂F∂τ = FηF ∂F∂τ , (A.21)¯ τ ( λ ) = F + Hλ = F + H √ γ vp , (A.22)Page 20here τ is the shear stress, γ = 2 ε = 2 ε = u , is the engineering shear strain, F is the initial yieldstrength. h is a hardening/softening material parameters with units of pressure (MPa), while η is theviscosity parameter with units of time s .Applying the procedure described above and taking advantage of the viscosity potential F ( τ, λ, ˙ λ ) as-suming the consistency condition ˙ F = 0, we obtain the following:˙ F = ∂F∂τ ˙ τ + ∂F∂λ ˙ λ + ∂F∂ ˙ λ ¨ λ = ∂F∂τ ˙ τ − ∂ ¯ τ∂λ ˙ λ − ηF ¨ λ, (A.23)˙ τ = G ( ˙ γ − ˙ γ vp ) . (A.24)Multiplying the above equation (A.24) with ∂F∂τ and replacing ∂F∂τ ˙ τ we get:˙ τ √ G ( ˙ γ − ˙ γ vp ) √ , (A.25) ηF ¨ λ + ∂ ¯ τ∂λ ˙ λ = G ( ˙ γ − ˙ γ vp ) √ , (A.26)separating ˙ λ we obtain:˙ λ = G ˙ γ √ ∂ ¯ τ∂λ + 3 G + ηF ¨ λ ∂ ¯ τ∂λ + 3 G , (A.27)substituting ˙ λ to the original equation for the calculation of ˙ τ :˙ τ = G (cid:32) ˙ γ − G ∂ ¯ τ∂λ + 3 G (cid:33) + GηF √ λ ∂ ¯ τ∂λ + 3 , (A.28)substituting G ∂ ¯ τ∂λ = HG = ¯ h (linear mechanical softening) we arrive finally at:˙ τ = G (cid:18) ¯ h ¯ h + 3 (cid:19) ˙ γ + ηF ¯ h + 3 ¨ γ vp . (A.29)It should be noted that the results derived until now can be also obtained using the Perzyna modelinstead of the consistency approach in the case of monotonic loading (no stress reversal, so that therate dependent unloading overstress of the Perzyna model is not taken into account). Using the Perzynamaterial we assume the existence of the viscoplastic potential Ω( τ, λ, ˙ λ ) which constitutes a region outsidethe yield function that the stress vector τ is indeed applicable. The time derivative of the yield function˙ F is the given as:˙ F = ηF ¨ λ = ∂F∂τ ˙ τ + ∂F∂λ ˙ λ = ∂F∂τ ˙ τ − ∂ ¯ τ∂λ ˙ λ. (A.30)Using the same arguments as before, (multiplying (A.24) by ∂F∂τ and substituting (A.30)) we arrive atthe same expression for ˙ λ and ˙ τ .Perturbing and replacing (A.29) in equation (3) we obtain: G ¯ h h ∂ ˜ γ∂x + ηF h ∂ ˙˜ γ vp ∂x = ρ ¨˜ u, (A.31) G ¯ h h ∂ ˜ γ∂x + ηF h (cid:18) ∂ ˙˜ γ∂x − G ∂ ˙˜ τ∂x (cid:19) = ρ ¨˜ u, (A.32) G ¯ h h ∂ ˜ u∂x − ρ ∂ ˜ u∂t + ηF (3 + ¯ h ) (cid:20) ∂ ˜ u∂t∂x − v s ∂ ˜ u∂t (cid:21) = 0 , (A.33) G ¯ h ∂ ˜ u∂x − ∂ ˜ u∂t (3 + ¯ h ) Gv s + ¯ η vp G (cid:20) ∂ ˜ u∂t∂x − v s ∂ ˜ u∂t (cid:21) = 0 , (A.34)where v s = (cid:113) Gρ , ¯ η vp G = ηF . Page 21 .2.2. Normalizing the 1D elasto-viscoplastic wave equation. We consider ¯ u = uu c , ¯ t = tt c , ¯ x = xx c . Applying and differentiating into (A.34) we arrive at: x c v s t c ∂ ¯ u∂ ¯ t − ∂ ¯ u∂ ¯ x ∂ ¯ t ¯ η vp t c ¯ h + x c v s t c h ¯ h ∂ ¯ u∂ ¯ t − ∂ ¯ u∂ ¯ x = 0 . (A.35)Introducing the characteristic velocity v c = x c t c , the result is written as: v c v s ∂ ¯ u∂ ¯ t − ∂ ¯ u∂ ¯ x ∂ ¯ t ¯ η vp t c ¯ h + v c v s h ¯ h ∂ ¯ u∂ ¯ t − ∂ ¯ u∂ ¯ x = 0 . (A.36) A.2.3. Dispersion relationship
Inserting into the normalized equation (A.36) the nondimensional solution assuming both ¯ k, ¯ ω ∈ C :¯ u (¯ x, ¯ t ) = exp i (¯ k ¯ x − ¯ ω ¯ t ) , (A.37)we arrive at:¯ h ¯ k t c v s − i ¯ k v s ¯ η vp ¯ ω − (3 + ¯ h ) ¯ t c v c ¯ ω + iv c ¯ η vp ¯ ω = 0 . (A.38) B. Behavior of the dispersion relation near infinity.
B.1. Properties of the inverse transform ¯ ω = z We present here in more detail the inverse transform ¯ ω = z based on Brown et al. (2009). The inversetransform allows the analysis of the behavior of a function close to ∞ by inverting the independentvariable. Therefore, as the independent variable ¯ ω tends to ∞ , its inverse transform z tends to 0. Morespecifically the points lying outside the unit circle centered at the origin of the ¯ ω complex plane aremapped inside the circle, while the opposite is true for the points initially inside the unit circle. Thepoints lying on the unit circle remain on the circle.When a point ¯ ω = ¯ ω r + ¯ ω i i is the image of a nonzero point z = z r + z i i under the transformation ¯ ω = 1 /z ,the relationship between the real and imaginary parts in original ¯ ω and transformed z complex planesrespectively are given as: z r = ¯ ω r ¯ ω r + ¯ ω i and z i = − ¯ ω i ¯ ω r + ¯ ω i , (B.1)¯ ω r = z r z r + z i and ¯ ω i = − z i z r + z i . (B.2)When A,B,C,D ∈ R are real numbers satisfying the condition B + C > z r + z i ) + B z r + C z i + D = 0 (B.3)represents a circle or a line on the complex plane z. In particular, when A = 0 the equation of a line isreturned while when A (cid:54) = 0 by completing the squares we get: (cid:18) z r + B2A (cid:19) + (cid:18) z i + C2A (cid:19) = (cid:32) (cid:112) B + C − (cid:33) . (B.4)The above equation represents a circle under the condition mentioned previously. Similarly by substitu-tion of z r , z i we find:D(¯ ω r + ¯ ω i ) + B¯ ω r − C¯ ω i + A = 0 . (B.5)From the equations (B.3),(B.5) above it is clear that: • A circle in the z -plane (A (cid:54) = 0) not passing through the origin (D (cid:54) = 0) is transformed into a circlenot passing through the origin in the ¯ ω -plane. • A circle in the z -plane (A (cid:54) = 0) passing through the origin (D = 0) is transformed into a line notpassing through the origin in the ¯ ω -plane.Page 22 A line in the z -plane (A = 0) not passing through the origin (D (cid:54) = 0) is transformed into a circlepassing through the origin in the ¯ ω -plane. • A line in the z -plane (A (cid:54) = 0) passing through the origin (D = 0) is transformed into a line passingthrough the origin in the ¯ ω -plane.From the above remarks we conclude that in the two complex planes the directions of the real andimaginary axes coincide. Furthermore every line passing from the origin retains its direction.For our analyses we need to examine the behavior of ¯ k , (¯ ω ) along lines of constant ¯ ω r = c , ¯ ω i = c .The geometrical loci on the complex z -plane are given by equations: (cid:18) z r − c (cid:19) + z i = (cid:18) c (cid:19) , (B.6) z r + (cid:18) z i + 12 c (cid:19) = (cid:18) c (cid:19) . (B.7)We notice that due to the inverse transform properties, lines parallel to the real axis Re, lying in onehalf of the complex ¯ ω -plane are transformed into circles passing through the origin, whose center lies onthe opposite imaginary half of the z -plane. B.2. Application of the inverse mapping describing the point at complex infinity.
Applying this mapping in equation (23) yields:¯ k ( z ) = V z (cid:18) hiT − z (cid:19) (cid:18) ¯ hiT − z (cid:19) − , (B.8)where the following relations hold between the components of ¯ ω, z in their respective complex plane: z r = ¯ ω r ¯ ω r + ¯ ω i , z i = − ¯ ω i ¯ ω r + ¯ ω i . (B.9)Taking the limit as z → z → z → ¯ k (cid:18) z (cid:19) = V z (cid:32) −
12 3 + ¯ hiT z + 18 (cid:18) hiT (cid:19) z + ... (cid:33) − (cid:18)
12 3+¯ hiT z − (cid:16) hiT (cid:17) z + ... (cid:19) . (B.10)We notice the pattern in the denominator of the last term that we can replace with the Taylor series of − x around x → z → ¯ k ( z ) = V z (cid:32) −
12 3 + ¯ hiT z + 18 (cid:18) hiT (cid:19) z + ... (cid:33) (cid:32) (cid:32)
12 3 + ¯ hiT z + 18 (cid:18) hiT (cid:19) z + ... (cid:33) ++ (cid:32)
12 3 + ¯ hiT z + 18 (cid:18) hiT (cid:19) z + ... (cid:33) + ... . (B.11)From the above polynomial only the factor z tends to ∞ , therefore z = 0 ⇔ ¯ ω P → ∞ is a pole of firstorder Brown et al. (2009).For the real part of ¯ k at the pole ¯ ω P as z r tends to ∞ applying de l’ Hopital rule we can prove:lim ¯ z r → + Re (cid:2) ¯ k , (¯ z r + βi ) (cid:3) = lim ¯ z r → + ± V (cid:118)(cid:117)(cid:117)(cid:116) (3 + h ) + T z r h + T ¯ z r cos (cid:18) arg (cid:18) h − iT zr h − iT zr (cid:19)(cid:19) ¯ z r = ±∞ , (B.12)Page 23here β = 0 in the case of traveling waves.Furthermore, for the imaginary part of ¯ k ( z ) at the pole ¯ ω P as z r tends to ∞ we can prove that¯ k i → ± c ∈ R using the de l’ Hˆopital rule:lim ¯ z r → Im (cid:2) ¯ k , (¯ z r + 0 i ) (cid:3) = lim ¯ z r → ± V (cid:118)(cid:117)(cid:117)(cid:116) (3 + h ) + T z r h + T ¯ z r sin (cid:18) arg (cid:18) h − iT zr h − iT zr (cid:19)(cid:19) ¯ z r = ± c (B.13)The value of ¯ k i (¯ ω ) is independent of whether we move towards the left or the right of the real axis ¯ ω r .It only depends on the solution branch ¯ k (¯ ω ) , or ¯ k (¯ ω ) we follow. (see Figure 9). C. Stability and localization of a monochromatic sinusoidal propagating pulse
The stability of the homogeneous deformation for a strain-softening ( h <
0) strain-rate hardening ( g > c i , and the phase velocity, c f (see section3). Examining equation (20) we notice that the amplitude term can be described as an exponentialfunction in ¯ x traveling along ¯ x -axis with velocity c i . Similarly, the periodic part which travels with avelocity of c r . A monochromatic sinusoidal pulse, whose amplitude varies with time and distance, isgiven as:¯ p (¯ x, ¯ t ) = (cid:2) H (¯ k r ¯ x − ¯ ω r ¯ t ) − H (¯ k r ¯ x − ¯ ω r ¯ t − π ) (cid:3) ¯ u exp ( − ¯ k i ¯ x + ¯ ω i ¯ t ) exp ( i (¯ k r ¯ x − ¯ ω r ¯ t )) , (C.1)¯ p (¯ x, ¯ t ) = [ H (¯ x − c r ¯ t ) − H (¯ x − c r ¯ t − π )] ¯ u exp ( − (¯ x − c i ¯ t )) exp ( i (¯ x − c r ¯ t )) . (C.2)The Heaviside terms H ( · ) are multiplied to the original monochromatic solution to indicate the startand end of the monochromatic signal. Therefore, they travel with the velocity of the periodic wave. Inthis way we can describe the amplitude that corresponds to the wavelength of the pulse at a specifictime. Based on equation (C.2) the relationship between the velocities of the two exponentials comprisingthe pulse is indicative of the stability and possible strain localization of the solution. In particular thefollowing cases are possible. • c i < , c r > • c i > , c r > • c i > , c r < • c i < , c r < c i , c r refer to the wave moving opposite to the positive direction defined by thepositive ¯ x -axis. In the first case described above c i < , c r > x -axis while the amplitude towards the negative.Due to the construction of the amplitude function (negative exponential) this has as an effect that ev-ery perturbation is attenuating with time. Therefore stability of the reference solution of homogeneousdeformation is ensured and strain localization cannot take place.In the second case both amplitude and phase are moving towards the positive part of the ¯ x -axis asshown on the right part of Figure 19 and the left part of Figure 20. In this case the behavior of theperturbation is defined by the relative magnitudes of the velocities | c i | , | c r | . If | c i | < | c r | then the ve-locity of the negative exponential is lower than that of the perturbation. Therefore, the amplitude ofthe perturbation is attenuated and the reference solution is stable (see right part of Figure 19). In theopposite case, where the perturbation travels slower than the amplitude velocity, the perturbation grows,rendering the reference solution unstable (see left part of Figure 20). Since the amplitude is increasingthe fastest at the peak behind the pulse, displacement is localizing close to the tip and localization tothe smallest mesh dimension is inevitable.In the third case when c i < , c r <
0, again the amplitude function and the perturbation are travelingPage 24 igure 20: Conditions for the growth of the perturbation. Left: Evolution of the perturbation (blue curve) and its amplitude(red curve) at different times for c i > c r >
0. Right:Evolution of the perturbation (blue curve) and its amplitude (redcurve) at different times for 0 > c i > c r . The propagating pulse is the multiplication of the red and blue curves.Figure 21: Conditions for the decay and growth of the perturbation. Left: Evolution of the perturbation (blue curve) andits amplitude (red curve) at different times for 0 > c r > c i , decay of perturbation. Right: Evolution of the perturbation(blue curve) and its amplitude (red curve) at different times for c i > > c r , growth of perturbation. The propagatingpulse is the multiplication of the red and blue curves. towards the negative direction (see left part of Figure 20, right part of Figure 21). Again the questionof stability and localization is dependent on the relative magnitudes of the two velocities | c i | , | c r | . Inthis case, if the perturbation is traveling slower than the amplitude | c r | < | c i | then the amplitude of theperturbation is decreasing and no localization happens (see left part of Figure 21). When we considerthe case where | c r | > | c i | then the amplitude of the perturbation is increasing exponentially. The ampli-tude is increasing the fastest for the tip closer to the front of the pulse and localization to the smallestwavelength cannot be avoided (see right part of Figure 20.In the final case c i > , c r <
0, the perturbation is moving towards the negative part of the ¯ x -axiswhile the amplitude towards the positive (see right part of Figure 21). Due to the negative exponentialspatial profile of the amplitude function, the amplitude of the perturbation is always increasing thefastest at the tip in front of the pulse. Therefore, in this final case, the solution is unstable and strainlocalization is possible with the smallest possible wavelength. Figure 19: Conditions for decaying perturbation. Left: Evolution of the perturbation (blue curve, exp ( i (¯( x ) − c r ¯ t )) andits amplitude (red curve, exp ( − (¯ x − c i ¯ t ))) at different times for c i < < c r . Right: Evolution of the perturbation (bluecurve) and its amplitude (red curve) at different times for c r > c i >
0. The propagating pulse is the multiplication of thered and blue curves.
Page 25 eferences
Abellan, M.A., de Borst, R., 2006. Wave propagation and localisation in a softening two-phase medium.Computer Methods in Applied Mechanics and Engineering 195, 5011–5019. doi: .Arfken, G.B., Weber, H.J., 1999. Mathematical methods for physicists.Beltzer, A.I., 1988. Acoustics of solids. 1, Springer-Verlag Berlin Heidelberg. doi: .Benallal, A., 2005. On localization modes in coupled thermo-hydro-mechanical problems. ComptesRendus M´ecanique 333, 557–564.Bernard, A., Lowe, M.J.S., Deschamps, M., 2001. Guided waves energy velocity in absorbing and non-absorbing plates. The Journal of the Acoustical Society of America 110, 186–196. doi: .Bigoni, D., 2012. Nonlinear solid mechanics: bifurcation theory and material instability. CambridgeUniversity Press.Borcherdt, R.D., 1973. Energy and plane waves in linear viscoelastic media. Journal of GeophysicalResearch 78, 2442–2453. doi: .de Borst, R., Duretz, T., 2020. On viscoplastic regularisation of strain-softening rocks and soils. Inter-national Journal for Numerical and Analytical Methods in Geomechanics doi: .de Borst, R., Sluys, L.J., 1991. Localisation in a Cosserat continuum under static and dynamic loadingconditions. Computer Methods in Applied Mechanics and Engineering 90, 805–827. doi: .Brauer, F., Nohel, J., 1969. The Qualitative Theory of Ordinary Differential Equations: An Introduction.Dover Publications, New York.Brown, J.W., Churchill, R.V., et al., 2009. Complex variables and applications. Boston: McGraw-HillHigher Education,.Chambon, R., Caillerie, D., Viggiani, G., 2004. Loss of uniqueness and bifurcation vs instability: someremarks. Revue fran¸caise de g´enie civil 8, 517–535.Chester, F.M., Chester, J.S., 1998. Ultracataclasite structure and friction processes of the Punchbowlfault, San Andreas system, California. Tectonophysics 295, 199–221. doi: .Collins-Craft, N.A., Stefanou, I., Sulem, J., Einav, I., 2020. A cosserat breakage mechanics model forbrittle granular media. Journal of the Mechanics and Physics of Solids , 103975.Dafermos, C.M., Joseph, D.D., Leslie, F.M., 1986. The Breadth and Depth of Continuum Mechanics.January 1986. doi: .Datta, S.K., 1986. The Excitation and Propagation of Elastic Waves (J. A. Hudson). volume 28. doi: .De Borst, R., 1991. Simulation of strain localization: A reappraisal of the cosserat continuum. Engi-neering Computations 8, 317–332. doi: .Deschamps, M., Poir´ee, B., Poncelet, O., 1997. Energy velocity of complex harmonic plane waves inviscous fluids. Wave Motion 25, 51–60. doi: .Etse, G., Carosio, a., 1999. Constitutive equations and numerical approaches in rate dependent materialformulations. Mecom99 , 289–296.Farina, J.E.G., Main, I.G., 1985. Vibrations and Waves in Physics. The Mathematical Gazette 69, 242.doi: . Page 26erasik, V., Stastna, M., 2010. Complex group velocity and energy transport in absorbing media.Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 81. doi: .Germain, P., 1973. The method of virtual power in continuum mechanics. part 2: Microstructure. SIAMJournal on Applied Mathematics 25, 556–575.Hayes, W.D., 1970. Kinematic wave theory. Proceedings of the Royal Society of London. A. Mathematicaland Physical Sciences 320, 209–226.Heeres, O.M., Suiker, A.S., De Borst, R., 2002. A comparison between the Perzyna viscoplastic modeland the consistency viscoplastic model. European Journal of Mechanics, A/Solids 21, 1–12. doi: .Inc., W.R., 2020. Mathematica, Version 12.1. URL: . cham-paign, IL, 2020.Kounadis, 1995. Nonlinear Stability of Structures. volume 1. doi: , arXiv:arXiv:1011.1669v3 .Lemaitre, J., Chaboche, J.L., Benallal, A., Desmorat, R., 2020. M´ecanique des mat´eriaux solides-3e ´ed.Dunod.Loret, B., Prevost, J., 1990. Dynamic strain localization in elasto-(visco-)plastic solids, Part 1. Generalformulation and one-dimensional examples. Computer Methods in Applied Mechanics and Engineering83, 247–273.Lucarini, V., 2009. Evidence of dispersion relations for the nonlinear response of the lorenz 63 system.Journal of Statistical Physics 134, 381–400. doi: , arXiv:0809.0101 .Mainardi, F., 1984. On linear dispersive waves with dissipation, in: North-Holland Mathematics Studies.Elsevier. volume 97, pp. 307–317.Mainardi, F., 1987. Energy velocity for hyperbolic dispersive waves. Wave motion 9, 201–208.Marion, J.B., 2013. Classical dynamics of particles and systems. Academic Press.Maugin, G.A., 1992. The Thermomechanics of Plasticity and Fracture. Cambridge Texts in AppliedMathematics, Cambridge University Press. doi: .Maugin, G.A., 2007. Nonlinear kinematic wave mechanics of elastic solids. Wave Motion 44, 472–481.doi: .Mawhin, J., 2005. Alexandr mikhailovich lyapunov, thesis on the stability of motion (1892), in: LandmarkWritings in Western Mathematics 1640-1940. Elsevier, pp. 664–676.Muhlhaus, H.B., Vardoulakis, I., 1988. The thickness of shear hands in granular materials. G´eotechnique38, 331–331. doi: .Muschietti, L., Dum, C.T., 1993. Real group velocity in a medium with dissipation. Physics of Fluids B5, 1383–1397. doi: .Needleman, A., 1988. Material rate dependence and mesh sensitivity in localization problems. Computermethods in applied mechanics and engineering 67, 69–85.Outline, C., 2015. Viscoelasticity and Wave Propagation. doi: .Pain, H.J., Beyer, R.T., 1993. The Physics of Vibrations and Waves, Fourth Edition . The Journal ofthe Acoustical Society of America 94, 3529–3529. doi: .Parker, M., 2009. Group Velocity. Solid State and Quantum Theory for Optoelectronics , 789–795doi: .Peirce, D., Shih, C.F., Needleman, A., 1984. A tangent modulus method for rate dependent solids.Computers and Structures 18, 875–887. doi: .Page 27oeverlein, H., 1962. Sommerfeld-runge law in three and four dimensions. Physical Review 128, 956.Rattez, H., Stefanou, I., Sulem, J., 2018a. The importance of Thermo-Hydro-Mechanical couplingsand microstructure to strain localization in 3D continua with application to seismic faults. Part I:Theory and linear stability analysis. Journal of the Mechanics and Physics of Solids 115, 54–76. URL: https://doi.org/10.1016/j.jmps.2018.03.004 , doi: .Rattez, H., Stefanou, I., Sulem, J., Veveakis, M., Poulet, T., 2018b. Numerical analysis of strain local-ization in rocks with thermo-hydro-mechanical couplings using cosserat continuum. Rock Mechanicsand Rock Engineering 51, 3295–3311.Rattez, H., Stefanou, I., Sulem, J., Veveakis, M., Poulet, T., 2018c. The importance of Thermo-Hydro-Mechanical couplings and microstructure to strain localization in 3D continua with application toseismic faults. Part II: Numerical implementation and post-bifurcation analysis. Journal of the Me-chanics and Physics of Solids doi: .Rice, J.R., 1976. Localization of Plastic Deformation. Theoretical and Applied Mechanics 1, 207–220.doi: .Rudnicki, J.W., Rice, J.R., 1975. Conditions for the localization of deformation in pressure-sensitivedilatant materials. Journal of the Mechanics and Physics of Solids 23, 371–394. doi: .Sibson, R.H., 2003. Thickness of the Seismic Slip Zone. Bulletin of the Seismological Society of America93, 1169–1178. URL: https://doi.org/10.1785/0120020061 , doi: .Sluys, L.J., de Borst, R., 1992. Wave propagation and localization in a rate-dependent crackedmedium-model formulation and one-dimensional examples. International Journal of Solids and Struc-tures 29, 2945–2958. URL: http://dx.doi.org/10.1016/0020-7683(92)90151-I , doi: .Sluys, L.J., de Borst, R., M¨uhlhaus, H.B., 1993. Wave propagation, localization and dispersion in agradient-dependent medium. International Journal of Solids and Structures 30, 1153–1171. doi: .Smith, M., 2009. ABAQUS/Standard User’s Manual, Version 6.9. Dassault Syst`emes Simulia Corp,United States.Stefanou,, I., 2021. Numerical geolab. URL: . nantes, FR, 2021.Stefanou, I., Alevizos, S., 2016. Fundamentals of bifurcation theory and stability analysis. Modelling ofinstabilities and bifurcation in Geomechanics, 26th ALERT Doctoral School 725.Sulem, J., Stefanou, I., Veveakis, E., 2011. Stability analysis of undrained adiabatic shearing of a rocklayer with cosserat microstructure. Granular Matter 13, 261–268.Vardoulakis, I., 1996a. Deformation of water-saturated sand: I. uniform undrained deformation andshear banding. G´eotechnique 46, 441–456.Vardoulakis, I., 1996b. Deformation of water-saturated sand: Ii. effect of pore water flow and shearbanding. G´eotechnique 46, 457–472.Vardoulakis, I., 2009. Lecture notes on cosserat continuum mechanics with application to the mechanicsof granular media, in: 3rd National Meeting on Generalized Continuum Theories and Applications,Thessaloniki, Greece.Vardoulakis, I., 2018. Cosserat Continuum Mechanics: With Applications to Granular Media. volume 87.Springer.Vardoulakis, I., Sulem, J.., 1995. Bifurcation Analysis in Geomechanics (1st ed.). CRC Press. doi: .Wang, W.M., Sluys, L.J., De Borst, R., 1996a. Interaction between material length scale and imperfectionsize for localisation phenomena in viscoplastic media. European journal of mechanics. A. Solids 15,447–464. Page 28ang, W.M., Sluys, L.J., De Borst, R., 1996b. Interaction between material length scale and imperfectionsize for localization phenomena in viscoplastic media. European Journal of Mechanics, A/Solids 15,447–464.Wang, W.M., Sluys, L.J., De Borst, R., 1997. Viscoplasticity for instabilities due to strain softening andstrain-rate softening. International Journal for Numerical Methods in Engineering 40, 3839–3864.Wu, F.H., Freund, L.B., 1984. Deformation trapping due to thermoplastic instability in one-dimensionalwave propagation. Journal of the Mechanics and Physics of Solids 32, 119–132. doi: .Zbib, H.M., Aifantis, E.C., 1992. On the gradient-dependent theory of plasticity and shear banding.Acta Mechanica 92, 209–225. doi:10.1007/BF01174177