In this paper, we extend the dynamic spherical cavity expansion model for rate-independent materials developed inrefs. [1,2,3]toviscoplasticmedia. Forthatpurpose, wedescribethematerialbehaviorwithanisotropicPerzyna- type overstress formulation [4, 5] in which the material rate-dependence is controlled by the viscosity parameter η. The theoretical predictions of the cavity expansion model, which assumes that the cavity expands at constant velocity, are compared with ﬁnite element simulations performed in ABAQUS/Explicit [6]. The agreement between theory and numerical simulations is excellent for the whole range of cavitation velocities investigated, and for diﬀerent values of the parameter η. We show that, as opposed to the steady-state self-similar solutions obtained for rate-independent materials [1, 2, 3], the material viscosity leads to time-dependent cavitation ﬁelds and stress relaxation as the cavity enlarges. In addition, we also show that the material viscosity facilitates to model the shock waves that emerge at the highest cavitation velocities investigated, controlling the amplitude and the width of the shock front. Keywords: Dynamic cavity expansion, Shock waves, Viscoplasticity, Stress relaxation 1. Introduction Beginning with the pioneering works of Hopkins [7] and Goodier [8], dynamic cavity expansion theories have been widely applied to model high-velocity penetration of metallic targets. The cavity expansion models use the solution of the pressure for expansion of a cavity in an inﬁnite medium as the contact pressure experienced by the projectile during the penetration process. Over the years, several works have shown that penetration models based on cavity expansion theories provide fair predictions for the resisting force, ballistic limit and penetration depth [9, 10, 1, 11, 12, 13]. Mostofthecavityexpansiontheoriesassumethattheelastoplasticﬁeldsthatdevelopnearthecavityaresteady- state and self-similar [10, 14, 15, 1], such that the solution of the problem depends on only one space variable, and it ∗Corresponding author. E-mail address: jarmarti@ing.uc3m.es Preprint submitted to Acta Mechanica February 4, 2020 2 isindependentoftime. Steady-stateself-similarapproacheshavebeenappliedtomodelpenetrationoftargetswhich exhibit a range of material response including compressibility, anisotropy, strain hardening, pressure sensitivity and porosity [10, 1, 13, 16, 17, 18]. However, since viscosity leads to time-dependent cavitation ﬁelds [15, 19, 20], steady- state similarity solutions are not possible for rate-dependent materials, imposing a limitation on the application of cavity expansion theories to model the penetration of targets with strain-rate dependent response. An important exception is the work of Warren and Forrestal [19], in which the cavity expansion model derived by Hopkins [7] was extended to consider the strain-hardening and strain-rate eﬀects on the penetration of 6061-T651 aluminum targets impacted by spherical-nosed rods. Warren and Forrestal [19] showed that, for strain-rate sensitive materials, the cavitation ﬁelds depend on the current cavity radius. Moreover, the problem of dynamic expansion of pressurized cavities in viscoplastic compressible solids, with application to geological materials, was addressed by Cleja-Tigoiu et al. [21] to study the crater formation by hypervelocity projectile impact in concrete structures. The authors provided a complete description of the distribution of stresses and strains as functions of time and distance from the center of the cavity, and for a wide range of applied cavity pressures. Moreover, Buchely and Marañon [22, 23] compared the cavitation ﬁelds obtained with the steady-state self- similar model developed by Masri and Durban [2] for rate-independent Mises solids with ﬁnite element simulations performed for strain-hardening and strain-rate dependent J2 materials. The computations of Buchely and Marañon [22, 23] conﬁrmed previous theoretical considerations reported in refs. [15, 19, 20], and showed that the cavitation ﬁelds for rate-dependent materials are functions of time, and cannot be modeled with steady-state self-similar approaches. However, Buchely and Marañon [22, 23] limited their attention to relatively low cavity expansion velocities (m≤0.1, see Eq. (27)) and the eﬀect of material viscosity on shock waves formation at high cavitation velocities was not addressed. Buchely and Marañon [22, 23] also did not pay attention to the eﬀect of stress relaxation that, for viscoplastic materials, plays an important role in the cavitation ﬁelds that develop near the cavity. In an eﬀort to provide further understanding of the eﬀect of material rate-dependence on dynamic cavitation, in this paper we extend the cavity expansion theory developed in refs. [1, 2, 3] to viscoplastic media. The material behavior is described with an isotropic Perzyna-type overstress formulation [4, 5], with associated ﬂow rule, in which the material rate-dependence is controlled by a single viscosity parameter η. Unlike the self-similar cavity expansion theories, the proposed model has two independent variables: the radial coordinate and the time. The theoretical results, which show the eﬀects of strain-rate dependence and stress relaxation in the stress ﬁelds that develop near the cavity, are compared with ﬁnite element simulations performed in ABAQUS/Explicit [6]. Excellent agreement between theoretical predictions and numerical results is obtained for the whole range of cavitation velocities investigated ( 0.1≤m≤0.6, see Eq. (27)). To be highlighted that, in agreement with the results reported by Molinari and Ravichandran [24], both theoretical and ﬁnite element results show that material viscosity controls the amplitude and the width of the shock waves that emerge at the highest cavitation velocities investigated. 3 2. Constitutive framework We consider an isotropic von Mises material described with an overstress Perzyna-type constitutive model [4, 5] that can be represented by the Bingham rheological model shown in Fig. 1, see ref. [25]. The yield function takes the form: Φ =τe−τy (1) whereτyis the yield stress, which is assumed to be constant, and τeis the eﬀective Kirchhoﬀ stress deﬁned by: τe=/radicalbigg3 2s:s,s=τ−τh1, τh=1 3τ:1 (2) whereτis the Kirchhoﬀ stress tensor, sis its deviatoric part, τhis the hydrostatic stress, and 1is the unit second order tensor. We have followed the work of Holzapfel [26] and used the Kirchhoﬀ stress in the formulation of the constitutive equations. Kirchhoﬀ stress appears to be an adequate stress measure to model ﬁnite strain metal plasticity, where the observed material hardening mainly results from isochoric dislocation accumulation. The rate of deformation tensor is taken to be the sum of an elastic, de, and a viscoplastic part, dvp: d=de+dvp(3) The elastic part is related to the rate of the stress by the following hypo-elastic law: Oτ=C:de(4) withOτbeing an objective stress rate and Cthe isotropic elastic tensor given by: C=E 1 +νI/prime+E 3 (1−2ν)1⊗1 (5) whereEis the Young’s modulus, νis the Poisson’s ratio, and I/primeis the unit deviatoric fourth order tensor. Moreover, assuming an associated ﬂow rule, the viscoplastic part of the rate of deformation tensor is: dvp=˙λ∂τe ∂τ(6) where ˙λis a non-negative viscoplastic multiplier given by the following kinetic equation [4]: 4 Symbol Property and units Value ρ0Initial density (kg/m3) 7800 EYoung’s modulus ( GPa), Eq. (5) 210 νPoisson’s ratio, Eq. (5) 0.3 τyYield stress ( MPa), Eq. (1) 240 ηViscosity parameter ( Pa·s), Eq. (7) 1.5, 150, 1500 Table 1: Material parameters. ˙λ=1 η/angbracketleftΦ/angbracketright (7) whereηis a viscosity parameter and /angbracketleftx/angbracketright=1 2(x+|x|)is the Macaulay operator. Moreover, the dissipation consistency: τ:dvp=τe˙¯εvp(8) leads to the identity: ˙λ=˙¯εvp(9) where ˙¯εvpis the eﬀective viscoplastic strain rate. The material parameters related to the elastic properties, Eq. (5), the yield stress, Eq. (1), and the viscoplastic kinetic relation, Eq. (7), are given in Table 1. Note that, for the sake of simplicity, we neglect the (potential) strain hardening and thermal softening of the material. Moreover, in Section 5 we carry out a parametric analysis to investigate the eﬀect of viscosity in the cavitation ﬁelds. In the parametric analysis diﬀerent values of ηare used, see Table 1. We are aware that most metals and alloys present a nonlinear viscoplastic behavior and display values of ηsmaller than 1500 Pa·s, see for instance ref. [27], where η≈280 Pa·sfor the highest strain rate sensitive branch /parenleftbig˙¯εvp≥1110 s−1/parenrightbigof a carbon steel. However, considering a wide range of values of ηhelps to bring to light the eﬀect of material viscosity in the cavitation ﬁelds and in the formation of shock waves at high cavitation velocities. In addition, assuming that ηis constant facilitates the interpretation of results. Note that actual materials presenting a nonlinear viscous behavior may be modeled considering a generalization of the viscoplastic framework used in this work. For this purpose, the function /angbracketleftΦ/angbracketrightshould be modiﬁed accordingly (e.g. refs. [28, 29, 30]). Nevertheless, including nonlinear viscosity in the theoretical and ﬁnite element models is left for a future work. 3. Theoretical model In this section we extend the dynamic spherical cavity expansion model developed in refs. [1, 2, 3] to viscoplastic materials. We also incorporate into the model the artiﬁcial viscosity approach used in refs. [31, 18] to capture the 5 formation of plastic shock waves at high cavitation velocities in rate-independent materials. The goal is to study whether this artiﬁcial viscosity is required to model shocks in viscoplastic media or if, on the contrary, the material viscosity is suﬃcient to regularize the problem and describe the emergence of shocks. Consider a pressurized spherical cavity of instantaneous radius aexpanding at constant velocity in an inﬁnite medium. We denote (r,θ,φ )the radial, circumferential and azimuthal Eulerian coordinates of a material point with reference to a spherical system with the origin located at the center of the cavity O, see Fig. 1. We introduce the dimensionless variable ξ=r/asuch that the derivative with respect to the radial coordinate is: ∂() ∂r=∂() ∂ξ∂ξ ∂r= ()/prime1 a(10) where the prime superscript denotes diﬀerentiation with respect to ξ. Moreover, the time derivative is [14]: ˙() =∂() ∂ξ˙ξ= ()/prime˙a a(v−ξ) (11) wherev=˙r/˙ais the dimensionless radial velocity. Since ˙ais constant (constant-velocity expansion), the dimension- less acceleration is ˙v=¨r/˙a. ξ Strain and stress-free undisturbed zone with ρ ρ0 Elastic zoneElasto-viscoplastic zoneξw ξi ξ=1 OξpE, ν τy τ τeMaterial behavior Figure 1: Scheme of the cavitation ﬁeld induced by a pressurized spherical cavity expanding dynamically in an elasto-viscoplastic medium. The internal pressure pis applied at the cavity wall where the dimensionless radial coordinate is ξ=r/a= 1. The remote ﬁeld, atξ > ξw, is undisturbed. Behind the elastic wave front, located at ξ=ξw, an elastic range develops. The transition between the elastic and elasto-viscoplastic regions is denoted by ξi, and the location of the viscoplastic shock wave by ξp(a viscoplastic shock wave emerges at high cavitation velocities). The material behaviour can be represented by the Bingham rheological model shown in the upper right part of the illustration. This ﬁgure is adapted from Durban and Masri [1]. 6 For spherical deformations, and isotropic material behavior, the non-zero components of the Kirchhoﬀ stress tensor areτrrandτθθ=τφφ, which leads to the following expression for the eﬀective Kirchhoﬀ stress: τe=τθθ−τrr (12) Similarly, the non-zero components of the rate of deformation tensor are: drr=∂˙r ∂r(13) dθθ=dφφ=˙r r(14) The elastic strain rates, Eq. (4), can be rewritten as: de rr=1 E˙τrr−2ν E˙τθθ (15) de θθ=de φφ=1−ν E˙τθθ−ν E˙τrr (16) In previous expressions, since the material spin is zero, the objective derivativeO ()(see Eq. (4)) has been replaced by a time derivative˙(). Moreover, the viscoplastic strain rates, Eq. (6), can be rewritten as: dvp rr=−˙¯εvp(17) dvp θθ=dvp φφ=1 2˙¯εvp(18) Using equations (10)-(18), the additive decomposition of the rate of deformation tensor, Eq. (3), leads to the following two expressions for the radial and circumferential directions: v/prime= (v−ξ) [(1−2ν)Trr−2νTe−¯εvp]/prime(19) v ξ= (v−ξ)/bracketleftbigg (1−2ν)Trr+ (1−ν)Te+1 2(¯εvp)/bracketrightbigg/prime (20) 7 whereTrr=τrr/EandTe=τe/Eare the dimensionless radial and eﬀective Kirchhoﬀ stresses, respectively, and ¯εvp=/integraltextt 0˙¯εvpdγis the eﬀective viscoplastic strain, where trefers to time. Moreover, during plastic ﬂow, we have that Φ≥0, which allows to rewrite the kinetic equation (7) as: (v−ξ) (¯εvp)/prime=a t∗˙a(Te−Ty) (21) where equations (1) and (11) have been used. For a constant-velocity expansion, i.e. ˙a=constant , the current cavityradiuscanbecalculatedas a=a0+ ˙at, wherea0isthecavityradiusat t= 0(initialcavityradius). Therefore, the dynamic cavity expansion model developed in this paper is time-dependent, thus requiring to specify the time tat which the corresponding cavitation ﬁelds are evaluated [15, 19]. Moreover, in previous expression t∗=η/E is the so-called relaxation time [5] and Ty=τy/Eis the dimensionless yield stress. Equation (21) shows that the rate-dependence of this problem is inversely proportional to the cavity radius a, leading to stress relaxation as the cavity enlarges (see Section 5). For the limiting cases in which a→∞(i.e.t→∞) ort∗→0(i.e.η→0, rate-independent material), the cavitation ﬁelds become steady-state and self-similar. Moreover, theratiobetweenthecurrentandinitialmaterialdensitiescanbeobtainedfromthemassconservation equation which, in Eulerian description, yields: ˙ρ+ρtr (l) = 0 (22) wherelis the spatial velocity gradient. In absence of material spin l=d, and, using equations (11), (13) and (14), previous expression can be rewritten as: (v−ξ)ρ/prime ρ+v/prime+ 2v ξ= 0 (23) Inserting (19) and (20) into equation (23), the following expression is obtained: ρ/prime ρ=−3 (1−2ν)T/prime h (24) whereTh=τh/Eis the dimensionless hydrostatic stress. Integrating equation (24) from a given coordinate ξto the elastic wave front ξw, whereρ=ρ0andTh= 0(see Fig. 1), the following relation for the density ratio is obtained: J−1= exp [−3 (1−2ν)Th] (25) whereJ=ρ0/ρis the material Jacobian (i.e. the determinant of the deformation gradient tensor). 8 Moreover, the balance of linear momentum along the radial direction, in Eulerian description, takes the form: ∂/parenleftbigJ−1¯τrr/parenrightbig ∂r+2J−1 r(¯τrr−¯τθθ) =ρ¨r (26) where ¯τrrand¯τθθarethetotalradialandcircumferentialKirchhoﬀstressesgivenbythesumofthe elasto-viscoplastic stresses and the viscousstresses, i.e. ¯τrr=τrr+τv rrand¯τθθ=τθθ+τv θθ. On the one hand, τrrandτθθare related to the elastic and viscoplastic deformation rates by the constitutive equations presented in Section 2 (i.e. equations (15)-(18) in this Section). On the other hand, τv rrandτv θθare related to the total deformation rates by means of an artiﬁcial viscosity [32] to be presented in the sequel. Using equations (10) and (11), and the relation J=ρ0/ρ, the linear momentum equation (26) can be rewritten as: ¯T/prime rr+J¯Trr/parenleftBig J−1/parenrightBig/prime−2 ξ¯Te=m2(v−ξ)v/prime(27) where ¯Trr=¯τrr Eand ¯Te=¯τe Eare dimensionless total Kirchhoﬀ stresses, and m=˙a/√ E/ρ0is the dimensionless cavitation velocity deﬁned as the ratio between the cavity expansion velocity ˙aand the elastic wave velocity in a long rod/radicalbig E/ρ0. Note that the artiﬁcial viscosity approach is only used in the deﬁnition of the stresses that enters into the balance of linear momentum (see refs. [31, 18]). Following the development presented in refs. [33, 32, 31, 18], the dimensionless viscous Kirchhoﬀ stresses, Tv rr=τv rr/EandTv θθ=τv θθ/E, are calculated as: Tv rr=  −ηrrm|v/prime|ifv/prime<0 0 if v/prime≥0(28) and Tv θθ=  −ηθθm|v|ifv<0 0 if v≥0(29) being ηrr=c2 0(∆ξ)2ρ ρ0m/vextendsingle/vextendsinglev/prime/vextendsingle/vextendsingle+cL∆ξ/radicalbiggρ ρ0(30) and ηθθ=c2 0/parenleftbigg∆ξ ξ/parenrightbigg2ρ ρ0m|v|+cL∆ξ ξ/radicalbiggρ ρ0(31) the radial and circumferential artiﬁcial viscosity parameters, respectively. In equations (30) and (31), c0andcLare 9 constant parameters (their values are given in Section 5.2), and ∆ξis a given increment along the dimensionless radial coordinate ξ. Recall that in refs. [31, 18] the artiﬁcial viscosity was essential to model plastic shocks in rate-independent materials. In this paper, as mentioned before, the idea is to investigate if the artiﬁcial viscosity is necessary to model shock waves in viscoplastic media. Moreover, to identify the emergence of a shock we use the approach developed in ref. [3] for rate-independent materials. Namely, we assume that a shock wave is formed if the following condition is met: m2(v−ξ)2+Trr−1 3 (1−2ν)= 0 (32) All the details of the derivation of expression (32) are given in ref. [3]. In the dynamic cavity expansion model developed in this paper, equations (19), (20), (21), (25), (27) and (28), provide a system of ﬁve diﬀerential equations and one algebraic relation with derivatives of six unknowns (v,¯εvp,J,Trr,Tv rr,Te), withTv θθbeing calculated from the algebraic relation (29). Integration is performed over the dimensionless radial coordinate ξ, from the elastic wave front ξwup to the cavity wall ξ= 1. The following boundary conditions have to be satisﬁed: v= 1, Trr=−JP, atξ= 1 (33) and v= 0, Trr=Te= 0, J = 1,atξ=ξw (34) whereP=p/Eis the dimensionless applied pressure (assumed constant, see Fig. 1). The compatibility of the ﬁeld variables at the interface between elastic and elasto-viscoplastic zones (see Fig. 1), ξ=ξi, has been enforced using the closed-form elastic solution derived in ref. [1]. The problem is solved using the numerical integration scheme proposed in refs. [1, 2, 3]. 4. Finite element model ThissectiondescribesthemainfeaturesoftheaxisymmetricﬁniteelementmodeldevelopedinABAQUS/Explicit [6] to simulate dynamic spherical cavitation. This same model has been used in refs. [34, 31]. The specimen is a sphere with radius Rs= 300mm which has a small cavity of radius Rc= 0.5mm in its center. To reduce the com- putational cost of the simulations, only the θ≥0half of the specimen has been analyzed (see Fig. 2). The sample is initially at rest and undeformed, while a constant internal pressure pis applied at the cavity wall (consistently with the theoretical model, see Fig. 1). Note that the dimensions of the sphere and the cavity are taken such that, for all the numerical simulations, the stress waves generated by the application of the cavity pressure reach the 10 Figure 2: Axisymmetric ﬁnite element model developed in ABAQUS/Explicit [6] to simulate dynamic spherical cavitation. This ﬁgure is adapted from refs. [34, 31]. external free boundary of the sphere after the cavity has started to expand at constant speed (i.e. cavitation has occurred [34]). The model has been meshed with 240,000four-node axisymmetric reduced integration elements, CAX4R in ABAQUS notation. The mesh shows radial symmetry, in order to retain the symmetry of the problem, with 200 elements along the circumferential direction and 1200 along the radial direction. The elements size is constant along the circumferential direction whereas it decreases along the radial direction as the cavity is approached with a bias ratio of 100. The dimensions of the elements located near the cavity are 12µm×4µm. Large gradients of stress and strain arise near the cavity, and within the shock wave that emerges at high cavitation velocities, which may lead to severe mesh distortion. In order to circumvent this problem, the Arbitrary Lagrangian Eulerian (ALE) adaptive meshing available in ABAQUS/Explicit [6] has been used, see refs. [34, 31]. Moreover, note that ABAQUS/Explicit [6] introduces artiﬁcial damping in the calculations in order to attenuate the numerical solution and ensure stability, see ref. [35]. The code generates bulk viscosity pressures, which are linear and quadratic, respectively, in the volumetric strain rate. This artiﬁcial viscosity is controlled by two parameters /pi1 (linear viscosity) and χ(quadratic viscosity). In the calculations shown in Section 5, the default values /pi1= 0.06 andχ= 1.2have been used. For the implementation of the constitutive model presented in Section 2 in ABAQUS/Explicit [6] a user material subroutine VUMAT was developed. The integration of the constitutive equations has been performed using the 11 radial return mapping algorithm, see refs. [36, 37]. 5. Data analysis and discussion This section is split into two parts. In Section 5.1, the ﬁnite element simulations are used to identify the time required to reach the constant-velocity expansion of the cavity, for a wide range of dimensionless cavitation velocities 0.1≤m≤0.6, and the three values of the viscosity parameter given in Table 1. Similar range of values ofmwas investigated in refs. [16, 31, 18]. In Section 5.2, the elasto-viscoplastic cavitation ﬁelds predicted by the theoretical model are compared with the ﬁnite element results for diﬀerent loading times. 5.1. Key ﬁnite element results Figure 3 shows the radial velocity of a node located at the cavity wall ˙aversus the loading time tfor ﬁnite element simulations performed using three values of the viscosity parameter η, namely: 1.5,150and1500 Pa·s (see Table 1). The node taken for the velocity measurements is indicated in Fig. 2. Results are reported for four diﬀerent applied pressures, p= 3696,10479,36645and78351MPa, such that the corresponding cavitation velocities are ˙a≈520,1040,2080and3120m/s, respectively (Figs. 3(a), 3(b), 3(c) and 3(d)). These cavitation velocities correspond to m= 0.1,0.2,0.4and0.6(recall from Section 3 that m=˙a/√ E/ρ0). Note that, for a given applied pressure, the cavitation velocity is virtually independent of the material viscosity. This behaviour is, most likely, due to the fact that for low values of mthe inﬂuence of viscosity in the cavitation ﬁelds is small, and for large values of mthe stress state is nearly hydrostatic, see Section 5.2. Note also that the relationship between the applied pressure and the cavitation velocity is nonlinear, in agreement with the results reported in refs. [16, 31, 18] for rate-independent materials. Moreover, the ˙a×tcurves display a concave-downward shape, for the four values ofm, and the three viscosity parameters considered. The rate of increase of the expansion velocity slows down astincreases such that ˙aeventually becomes constant, indicating the onset of cavitation. The time to reach the constant-velocity expansion, hereinafter denoted by tcav, and referred to as cavitation time, is taken as the loading time for which the condition ¨a≤1×106m/s2is met. The condition ¨a≤1×106m/s2is chosen because for this value of ¨athe results for the three viscosity parameters investigated virtually overlap each other. Fig. 4 shows the cavitation time tcavversus the dimensionless cavitation velocity m, for the three values of η investigated. The cavitation time tcavdecreases nonlinearly with m(see also refs. [34, 31]), such that the tcav×m curves display a concave-upward shape, with faster decrease of tcavfor short values of m. Moreover, the cavitation time is greater as the material viscosity increases, i.e. the tcav×mcurve is shifted upwards as the value of η increases for the whole range of cavitation velocities investigated. The inﬂuence of viscosity in the cavitation time is greater for small values of m. The cavitation times obtained for η= 1500 Pa·sare from 2to7times greater than those obtained for η= 1.5 Pa·sandη= 150 Pa·s. It becomes apparent that material viscosity delays cavitation, extending the transient behavior that precedes the constant-velocity expansion of the cavity. Moreover, note that 12 0100200300400500600 0 10 20 30 40 50 601.5 150 1500 Loading time, t ( μs)Finite element results η=1.5 Pa s η=150 Pa s η=1500 Pa sm=0.1 . . .Radial velocity of the cavity wall, a (m/s). (a) 020040060080010001200 0 5 10 15 20 25 30B D FFinite element results m=0.2 Loading time, t ( μs)η=1.5 Pa s η=150 Pa s η=1500 Pa s. . .Radial velocity of the cavity wall, a (m/s). (b) 04008001200160020002400 0 5 10 15B D F Loading time, t ( μs)Radial velocity of the cavity wall, a (m/s)Finite element results m=0.4. η=1.5 Pa s η=150 Pa s η=1500 Pa s. . . (c) 060012001800240030003600 0 1 2 3 4 5B D F Loading time, t ( μs)Finite element results m=0.6 η=1.5 Pa s η=150 Pa s η=1500 Pa s. . .Radial velocity of the cavity wall, a (m/s). (d) Figure 3: Finite element results. Radial velocity of a node located at the cavity wall ˙aversus the loading time tfor three values of the viscosity parameter η= 1.5,150and1500 Pa ·s. Results are shown for four diﬀerent dimensionless cavitation velocities: (a) m= 0.1, (b)m= 0.2, (c)m= 0.4, (d)m= 0.6. 13 01020304050 0 0.1 0.2 0.3 0.4 0.5 0.6 0.71.5 150 1500Finite element resultsη=1.5 Pa s η=150 Pa s η=1500 Pa s. . .Cavitation time, tcav (μs) Dimensionless cavitation velocity, m Figure 4: Finite element results. Cavitation time tcavversus dimensionless cavitation velocity mfor three values of the viscosity parameterη= 1.5,150and1500 Pa ·s. in ref. [34] it was stated that, if the penetration process is much faster than the time required for constant-velocity cavity expansion to appear, then constant-velocity cavitation ﬁelds are not suﬃcient to describe the underlying physical phenomena. Therefore, the results presented in this paper, which reveal that viscosity delays cavitation, suggest that the increase of material rate-dependence may hinderthe application of constant-velocity cavitation ﬁelds to model high-velocity penetration problems. Nevertheless, this is still an open issue that requires further research. 5.2. Veriﬁcation of the theoretical model The veriﬁcation is performed comparing the stress ﬁelds predicted by the theoretical model with the ﬁnite element results. The comparison is carried out for loading times greater than tcav, for which the constant-velocity expansion solution developed in Section 3 is valid. Moreover, despite that the Kirchhoﬀ stress tensor τwas used to formulate the constitutive framework and the theoretical model (see Sections 2 and 3), with the aim of facilitating the interpretation of the results, since the most commonly used measure of stress is the Cauchy stress tensor σ, the stress ﬁelds shown in the sequel are given in terms of the components of the dimensionless Cauchy stress tensor Σ=σ/E. The tensor Σis readily calculated as Σ=J−1T, whereT=τ/Eis the dimensionless Kirchhoﬀ stress tensor (see Section 3). Moreover, note that all the theoretical stress ﬁelds correspond to the elasto-viscoplastic stresses (see Section 3). The viscous stresses lack of physical meaning, and their role in the theoretical model is limited to facilitate the integration of the governing equations when shock waves appear at high cavitation velocities and the viscosity of the material is low [31, 18]. If the material viscosity is high, the artiﬁcial viscosity is no longer needed to model the formation of shocks, as will be discussed further below. 14 Figure 5 compares theoretical (TM) and ﬁnite element results (FEM) for m= 0.1. Thislowcavitation velocity does not lead to the emergence of shock waves, the artiﬁcial viscosity was not necessary to obtain the theoretical predictions, and the parameters c0andcLwere set to zero (see Section 3). Recall from Section 3 that we consider a shock wave to form when the condition (32) is met. Figures 5(a)-(b) show the evolution of the dimensionless radial Cauchy stress Σrrversus the dimensionless radial coordinate ξforη= 1.5 Pa·sandη= 1500 Pa·s, respectively. Note that the material viscosity hardly inﬂuences the radial stress ﬁelds. The Σrr×ξcurves, as they are displayed in this paper with increasing negative values of the radial stress along the y-axis, show a concave-upward shape, with the maximum value of the radial stress Σrr=−0.0177(compression) located at the cavity wall ξ= 1, and the elastic wave front located at ξw= 11.60. The agreement between theoretical predictions and ﬁnite element results is very good. In addition, the results reported for diﬀerent loading times, for both viscosity parameter values, virtually overlap each other, which indicates that the radial stress ﬁelds are nearlysteady-state and self-similar. Namely, for η= 1.5 Pa·sresults are shown for t= 10µs and 20µs, and forη= 1500 Pa·sthe loading times considered are t= 40µs,50µs andt→∞ (only theoretical predictions). Recall that the limiting case of t→∞corresponds to the rate-independent material response. Moreover, note that, as anticipated, the loading times are equal or greater than the cavitation times tcavreported in Fig. 4. Namely, the smaller loading times considered correspond, approximately, to the cavitation time, see Fig. 4. Figures 5(c)-(d) show the variation of the dimensionless eﬀective Cauchy stress Σewith the dimensionless radial coordinate ξfor the same two values of the viscosity parameter. Note that, near the cavity, the eﬀective stress is one order of magnitude smaller than the absolute value of the radial stress. On the other hand, as in the case of the radial stress ﬁelds, the agreement between theoretical predictions and ﬁnite element results is very good. Forη= 1.5 Pa·s, Fig. 5(c), the eﬀective stress slightly decreases as we move away from the cavity wall until the elastic/elasto-viscoplastic interface is reached at ξi= 6.04. Forξi<ξ <ξw, which corresponds to the elastic zone, the eﬀective stress Σeshows a steeper decline, featuring a concave-upward shape. Note that the results obtained for the two loading times considered, t= 10µs andt= 20µs, are very similar, showing that the eﬀective stress ﬁelds, for this value of the viscosity parameter, are nearlysteady-state and self-similar. In contrast, for η= 1500 Pa·s, Fig. 5(d), there are meaningful diﬀerences between the results obtained near the cavity for diﬀerent loading times. The eﬀective stress decreases as the loading time increases, i.e. for this value of the viscosity parameter the relaxation of the material as the cavity increases inﬂuences the eﬀective stress ﬁelds. The decrease of the eﬀective Cauchy stress with the radial coordinate within the elasto-viscoplastic zone for t→∞(rate-independent limit) is due to the inﬂuence of material compressibility in Σe. Recall that the yield criterion is formulated in terms of the Kirchhoﬀ stress, see Eq. (1). Figure 6 compares the theoretical and ﬁnite element stress ﬁelds for a greater cavity expansion velocity m= 0.2. 15 -0.02 -0.015 -0.01 -0.005 0 1 3 5 7 9 11 13 15B D F HDimensionless Cauchy radial stress, Σrr Dimensionless radial coordinate, ξTM - t=10 μs TM - t=20 μs FEM - t=10 μs FEM - t=20 μs m=0.1 η=1.5 Pa s .ξw=11.60 (a) -0.02 -0.015 -0.01 -0.005 0 1 3 5 7 9 11 13 15TM40 TM50 TMINF NM40 NM50Dimensionless Cauchy radial stress, Σrr Dimensionless radial coordinate, ξm=0.1 η=1500 Pa s .ξw=11.60TM - t=40 μs TM - t=50 μs TM - t FEM - t=40 μs FEM - t=50 μs8 (b) 00.0010.002 1 3 5 7 9 11 13 15B D F HTM - t=10 μs TM - t=20 μs FEM - t=10 μs FEM - t=20 μsDimensionless Cauchy effective stress, Σe Dimensionless radial coordinate, ξm=0.1 η=1.5 Pa sξw=11.60 ξi=6.04 . (c) 00.0010.002 1 3 5 7 9 11 13 15TM40 TM50 TMINF NM40 NM50 Dimensionless radial coordinate, ξξw=11.60 ξi=6.04Dimensionless Cauchy effective stress, Σe TM - t=40 μs TM - t=50 μs TM - t FEM - t=40 μs FEM - t=50 μs m=0.1 η=1500 Pa s .zoom8 0.0010.00120.00140.0016 1 1.5 2 2.5 (d) Figure 5: Comparison between theoretical results (TM) and ﬁnite element results (FEM) for m= 0.1. Dimensionless radial Cauchy stress Σrrversus dimensionless radial coordinate ξfor (a)η= 1.5 Pa ·sand (b)η= 1500 Pa ·s. Dimensionless eﬀective Cauchy stress Σeversus dimensionless radial coordinate ξfor (c)η= 1.5 Pa ·sand (d)η= 1500 Pa ·s. The theoretical results were obtained without artiﬁcial viscosity, c0= 0andcL= 0. 16 This cavitation velocity neither leads to the emergence of shocks (the condition (32) is not met), so that we set c0= 0andcL= 0in the theoretical model. As in Fig. 5, results are shown for η= 1.5 Pa·sandη= 1500 Pa·s. The agreement between theoretical and ﬁnite element calculations is very good. Figures 6(a)-(b) show the evolution of the radial Cauchy stress Σrrversus the dimensionless radial coordinate ξ. The inﬂuence of the material viscosity on the Σrr×ξcurves is very small. Moreover, the radial stress ﬁelds obtained for diﬀerent loading times practically overlap each other. The diﬀerence with the Σrr×ξcurves shown in Figs. 5(a)-(b) is that the absolute value of the radial stress near the cavity is ≈3times greater. The evolution of the eﬀective Cauchy stress Σewith the dimensionless radial coordinate ξis displayed in Figs. 6(c)-(d). Compared to the results reported in Figs. 5(c)-(d), the Σe×ξcurves form= 0.2show increased gradients of eﬀective stress near the cavity. The increase of the cavitation velocity leads to greater strain rates and greater material compressibility near the cavity wall, which boosts the eﬀective Cauchy stress. The greater the cavitation velocity and the material viscosity, the more important the eﬀective Cauchy stress gradients near the cavity are. Moreover, note that in the case of η= 1500 Pa·s, Fig. 6(d), the stress relaxation leads to signiﬁcant diﬀerences in the eﬀective stress near the cavity for the loading times considered, which illustrates that the material viscosity leads to time-dependent cavitation ﬁelds [38]. As the loading time increases, the gradients of stress in the elasto- viscoplastic zone get gradually reduced. As mentioned before, the decrease of the eﬀective Cauchy stress near the cavity fort→∞(rate-independent limit) is due to the inﬂuence of material compressibility in Σe. Figures 7 and 8 show the theoretical and ﬁnite element stress ﬁelds for m= 0.4and three values of the viscosity parameterη= 1.5,150and1500 Pa·s. This cavitation velocity leads to the formation of a shock wave. The minimum value of mfor which the shock condition (32) is met is 0.34for the three viscosity parameters (recall that the shock condition was developed for rate-independent materials). Moreover, in the case of η= 1.5 Pa·s, the theoretical results were obtained using the artiﬁcial viscosity with c0= 0.55andcL= 0.25. In contrast, for η= 150 Pa·sandη= 1500 Pa·s, the theoretical solution could be obtained without artiﬁcial viscosity because thelargematerial viscosity facilitated the integration of the governing equations for the loading times investigated. However, for larger loading times, t>22µs andt>226µs, due to stress relaxation, the artiﬁcial viscosity becomes also necessary to capture the shock wave for η= 150 Pa·sandη= 1500 Pa·s, respectively. Moreover, note that for this cavitation velocity, the theoretical model also reproduces very well the ﬁnite element results. The evolution of the radial Cauchy stress Σrrversus the dimensionless radial coordinate ξforη= 1.5,150and 1500 Pa·sis shown in Figs. 7(a)-(b)-(c), respectively. As for the lower cavitation velocities investigated before, the material viscosity and the loading time hardly aﬀect the radial stress ﬁelds. On the other hand, an important diﬀerence with the results reported for m= 0.1andm= 0.2is that the Σrr×ξcurves, as they are displayed in this paper with increasing negative values of the radial stress along the y-axis, do not show a concave-upward shape. Form= 0.4, the radial stress decreases almost linearly as we move away from the cavity wall, with a quick drop 17 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 1 2 3 4 5 6 7 8TM5 TM15 NM5 NM15 Dimensionless radial coordinate, ξξw=5.80Dimensionless Cauchy radial stress, Σrr TM - t=5 μs TM - t=15 μs FEM - t=5 μs FEM - t=15 μs m=0.2 η=1.5 Pa s .ξi=4.18 (a) -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0 1 2 3 4 5 6 7 8TM20 TM40 TMinf NM20 NM40 Dimensionless radial coordinate, ξξw=5.80Dimensionless Cauchy radial stress, Σrr TM - t=20 μs TM - t=40 μs TM - t FEM - t=20 μs FEM - t=40 μs m=0.2 η=1500 Pa s .ξi=4.18 8 (b) 00.00050.0010.00150.0020.0025 1 2 3 4 5 6 7 8TM5 TM15 NM5 NM15Dimensionless Cauchy effective stress, Σe ξw=5.80 Dimensionless radial coordinate, ξm=0.2 η=1.5 Pa s .TM - t=5 μs TM - t=15 μs FEM - t=5 μs FEM - t=15 μs ξi=4.18 (c) 00.00050.0010.00150.0020.0025 1 2 3 4 5 6 7 8TM20 TM40 K NM20 NM40Dimensionless Cauchy effective stress, Σe TM - t=20 μs TM - t=40 μs TM - t FEM - t=20 μs FEM - t=40 μs Dimensionless radial coordinate, ξξw=5.80 ξi=4.18m=0.2 η=1500 Pa s . 8 (d) Figure 6: Comparison between theoretical results (TM) and ﬁnite element results (FEM) for m= 0.2. Dimensionless radial Cauchy stress Σrrversus dimensionless radial coordinate ξfor (a)η= 1.5 Pa ·sand (b)η= 1500 Pa ·s. Dimensionless eﬀective Cauchy stress Σeversus dimensionless radial coordinate ξfor (c)η= 1.5 Pa ·sand (d)η= 1500 Pa ·s. The theoretical results were obtained without artiﬁcial viscosity, c0= 0andcL= 0. 18 atξp= 2.29which corresponds to the shock wave. The shock wave location virtually coincides with the elastic / elasto-viscoplastic interface ξp≈ξi(see also refs. [31, 18, 3]). Figures 8(a)-(b)-(c) display the variation of the eﬀective Cauchy stress Σewith the dimensionless radial coor- dinateξ. Forη= 1.5 Pa·s, Fig. 8(a), the eﬀective stress ﬁeld is hardly dependent on the loading time considered. However, the time dependence of the eﬀective stress ﬁelds becomes more important as the viscosity parameter increases. For η= 150 Pa·s, Fig. 8(b), the value of Σenear the cavity is slightly smaller for t= 10µsthan fort= 5µs. Moreover, for the greatest material viscosity considered η= 1500 Pa·s, Fig. 8(c), the value of the eﬀective stress at the cavity wall for 7µs is≈1.75times greater than for 20µs. Note that, for the latter loading time, near the cavity, the eﬀective stress is approximately 80times smaller than the absolute value of the radial stress, i.e. the triaxiality is roughly −80. Note that for high cavitation velocities and applied pressures, the radial, circumferential and azimuthal stresses are similar, and much greater in absolute value than the equivalent stress (see Eq. (12)). Therefore, the ratio between the radial and the eﬀective stress is roughly the stress triaxiality. The large hydrostatic pressure explains that, while the eﬀective stress ﬁelds are strongly dependent on the loading time, the radial stress ﬁelds are not (compare Figs. 7 and 8). Moreover, the peakof eﬀective stress located ξp= 2.29 corresponds to the shock wave. The formation of a shock leads to a local increase of the strain rate [39, 31], which boosts the eﬀective stress. As the loading time increases, the peak of stress becomes smaller (as further discussed below). Figures 9 and 10 show the radial and eﬀective Cauchy stress ﬁelds predicted by the theoretical model and the ﬁnite element simulations for the greatest cavitation velocity investigated in this work m= 0.6. Note that the intensity of the shock increases with the cavitation velocity. The artiﬁcial viscosity with c0= 0.55andcL= 0.25 was used to obtain the results for η= 1.5and150 Pa·s(note that in Figs. 7 and 8 the results for m= 0.4 andη= 150 Pa·swere obtained without artiﬁcial viscosity). In contrast, in the case of η= 1500 Pa·s, for the loading times considered, the artiﬁcial viscosity was not needed to integrate the theoretical model. Nevertheless, fort >12µs the artiﬁcial viscosity becomes also necessary to capture the shock wave for η= 1500 Pa·s. As for the lower cavitation velocities, the agreement between theoretical predictions and ﬁnite element calculations is excellent. Figures 9(a)-(b)-(c) show the radial Cauchy stress Σrrversus the dimensionless radial coordinate ξforη= 1.5, 150and1500 Pa·s, respectively. The Σrr×ξcurves, as they are displayed in this paper with increasing negative values of the radial stress along the y-axis, feature a concave-downward shape with a drastic drop at ξp= 1.68 which corresponds to the shock wave location. Note that, as in the case of m= 0.4shown in Figs. 7 and 8, we have thatξp≈ξi. Both in the theoretical model and the ﬁnite element simulations, the change in the radial stress induced by the shock, while abrupt, is not a discontinuity, i.e. the shock has a ﬁnite width due to the regularizing eﬀect of artiﬁcial/material viscosity. In the case of η= 1500 Pa·s, Fig. 9(c), it becomes apparent that the width 19 ξp=2.29-0.2 -0.15 -0.1 -0.05 0 1 2 3 4TM3 TM10 NM3 NM10 Dimensionless radial coordinate, ξDimensionless Cauchy radial stress, Σrr m=0.4 η=1.5 Pa s .TM - t=3 μs TM - t=10 μs FEM - t=3 μs FEM - t=10 μs ξw=2.90 (a) -0.2 -0.15 -0.1 -0.05 0 1 1.5 2 2.5 3 3.5 4TM5 TM10 NM5 NM10 Dimensionless radial coordinate, ξDimensionless Cauchy radial stress, Σrr m=0.4 η=150 Pa s .ξp=2.29 ξw=2.90TM - t=5 μs TM - t=10 μs FEM - t=5 μs FEM - t=10 μs (b) -0.2 -0.15 -0.1 -0.05 0 1 2 3 4TM7 TM20 NM7 NM20TM - t=7 μs TM - t=20 μs FEM - t=7 μs FEM - t=20 μs Dimensionless radial coordinate, ξDimensionless Cauchy radial stress, Σrr m=0.4 η=1500 Pa s .ξp=2.29 ξw=2.90 (c) Figure 7: Comparison between theoretical results (TM) and ﬁnite element results (FEM) for m= 0.4. Dimensionless radial Cauchy stress Σrrversus dimensionless radial coordinate ξfor (a)η= 1.5 Pa ·s, (b)η= 150 Pa ·sand (c)η= 1500 Pa ·s. The theoretical results for η= 1.5 Pa ·swere obtained with c0= 0.55andcL= 0.25. In contrast, the theoretical results for η= 150 Pa ·sand η= 1500 Pa ·swere obtained without artiﬁcial viscosity, c0= 0andcL= 0. 20 00.0010.0020.0030.0040.005 1 2 3 4TM3 TM10 NM3 NM10 Dimensionless radial coordinate, ξm=0.4 η=1.5 Pa s .ξw=2.90TM - t=3 μs TM - t=10 μs FEM - t=3 μs FEM - t=10 μsDimensionless Cauchy effective stress, Σe ξp=2.29 (a) 00.0010.0020.0030.0040.005 1 1.5 2 2.5 3 3.5 4TM5 TM10 NM5 NM10 Dimensionless radial coordinate, ξDimensionless Cauchy effective stress, Σe m=0.4 η=150 Pa s .TM - t=5 μs TM - t=10 μs FEM - t=5 μs FEM - t=10 μs ξp=2.29 ξw=2.90 (b) 00.0010.0020.0030.0040.005 1 2 3 4TM7 TM20 NM7 NM20 Dimensionless radial coordinate, ξDimensionless Cauchy effective stress, Σe m=0.4 η=1500 Pa s .TM - t=7 μs TM - t=20 μs FEM - t=7 μs FEM - t=20 μs ξp=2.29 ξw=2.90 (c) Figure 8: Comparison between theoretical results (TM) and ﬁnite element results (FEM) for m= 0.4. Dimensionless eﬀective Cauchy stress Σrrversus dimensionless radial coordinate ξfor (a)η= 1.5 Pa ·s, (b)η= 150 Pa ·sand (c)η= 1500 Pa ·s. The theoretical results for η= 1.5 Pa ·swere obtained with c0= 0.55andcL= 0.25. In contrast, the theoretical results for η= 150 Pa ·sand η= 1500 Pa ·swere obtained without artiﬁcial viscosity, c0= 0andcL= 0. 21 of the shock is smaller as the loading times increases, i.e. as the material relaxes, which shows the inﬂuence of material viscosity in the shock width. Moreover, the comparison with the results reported in Figs. 7 and 8 shows that the shock wave is located closer to the cavity as mincreases (see also refs. [3, 31, 18]). Outside the shock front, as for the lower cavitation velocities investigated, the material viscosity hardly aﬀects the radial stress ﬁelds, which virtually overlap each other for diﬀerent loading times and the three viscosity parameters considered (note that the hydrostatic stress is very large). The evolution of the eﬀective Cauchy stress Σewith the radial coordinate is shown in Figs. 10(a)-(b)-(c). For the smaller material viscosity, Fig. 10(a), the eﬀective stress is largely independent of the loading time. The decrease of Σewith the radial coordinate within the elasto-viscoplastic zone is, mostly, due to the eﬀect of material compressibility on the eﬀective Cauchy stress. The kink in the Σe×ξcurves corresponds to the elastic / elasto- viscoplastic interface ξi, which virtually coincides with the shock wave location (as mentioned before). The decrease of the eﬀective stress for ξ>ξicorresponds to the elastic zone. For greater values of the viscosity parameter, Figs. 10(b)-(c), the Σe×ξcurves become time dependent, as in the case of m= 0.4shown in Figs. 8(b)-(c). For both η= 150 Pa·sandη= 1500 Pa·s, within the elasto-viscoplastic zone, the eﬀective stress is greater for 5µs than for 10µs, i.e. the material relaxes with the loading time. The drastic increase of eﬀective stress located at ξp= 1.68 corresponds to the shock wave. Both ﬁnite element simulations and theoretical predictions show that the width of the shock wave and the maximum eﬀective stress attained within the shock wave decrease with the loading time. This result is consistent with the conclusions obtained in ref. [24], where it was shown that viscosity plays a fundamental role in the shock structure. 6. Concluding remarks and perspectives The main contribution of this paper is to extend the dynamic spherical cavity expansion theory developed in refs. [1, 2, 3] to viscoplastic Perzyna-type materials. The theoretical model has been compared with ﬁnite element simulations for a wide range of cavitation velocities. We have shown that the material viscosity: (i) increases the time needed for the cavity to reach the constant-velocity expansion, (ii) leads to time-dependent cavitation ﬁelds and(iii)facilitatestomodeltheshockwavesthatemergeathighcavitationvelocities. Moreover, therate-dependent cavity expansion model presented in this work could be readily extended to include, for instance, plastic anisotropy [18], porosity [16, 31], strain hardening [2] and thermal softening [40]. The application of the model to study the inﬂuence of material viscosity in the penetration resistance of rate-dependent materials is also left as future work. In this regard, note that the application of constant-velocity cavitation ﬁelds to model high-velocity penetration phenomena requires the transient behavior that precedes cavitation to be short which, attending to the results presented in this paper, may not be the case if the viscosity of the material is suﬃciently high. 22 -0.4 -0.3 -0.2 -0.1 0 1 1.2 1.4 1.6 1.8 2TM2.5 TM10 NM2.5 NM10Dimensionless Cauchy radial stress, Σrr Dimensionless radial coordinate, ξTM - t=2.5 μs TM - t=10 μs FEM - t=2.5 μs FEM - t=10 μs ξw=1.93m=0.6 η=1.5 Pa s . -0.2 0 1.6 1.8zoom ξp=1.68 (a) -0.4 -0.3 -0.2 -0.1 0 1 1.2 1.4 1.6 1.8 2TM5 TM10 NM5 NM10Dimensionless Cauchy radial stress, Σrr Dimensionless radial coordinate, ξm=0.6 η=150 Pa s .TM - t=5 μs TM - t=10 μs FEM - t=5 μs FEM - t=10 μs ξw=1.93 ξp=1.68-0.2 0 1.6 1.8zoom (b) -0.4 -0.3 -0.2 -0.1 0 1 1.2 1.4 1.6 1.8 2TM5 TM10 NM5 NM10TM - t=5 μs TM - t=10 μs FEM - t=5 μs FEM - t=10 μsDimensionless Cauchy radial stress, Σrr Dimensionless radial coordinate, ξm=0.6 η=1500 Pa s . ξw=1.93-0.2 0 1.6 1.8zoom ξp=1.68 (c) Figure 9: Comparison between theoretical results (TM) and ﬁnite element results (FEM) for m= 0.4. Dimensionless radial Cauchy stress Σrrversus dimensionless radial coordinate ξfor (a)η= 1.5 Pa ·s, (b)η= 150 Pa ·sand (c)η= 1500 Pa ·s. The theoretical results for η= 1.5 Pa ·sandη= 150 Pa ·swere obtained with c0= 0.55andcL= 0.25. In contrast, the theoretical results for η= 1500 Pa ·swere obtained without artiﬁcial viscosity, c0= 0andcL= 0. 23 00.0040.0080.012 1 1.2 1.4 1.6 1.8 2TM2.5 TM10 NM2.5 NM10TM - t=2.5 μs TM - t=10 μs FEM - t=2.5 μs FEM - t=10 μs ξw=1.93m=0.6 η=1.5 Pa s .Dimensionless Cauchy effective stress, Σe Dimensionless radial coordinate, ξ (a) 00.0040.0080.012 1 1.2 1.4 1.6 1.8 2TM5 TM10 NM5 NM10 ξp=1.68 ξw=1.93Dimensionless Cauchy effective stress, Σe Dimensionless radial coordinate, ξm=0.6 η=150 Pa s .TM- t=5 μs TM - t=10 μs FEM - t=5 μs FEM - t=10 μszoom 00.004 1.6 1.8 (b) 00.0040.0080.012 1 1.2 1.4 1.6 1.8 2TM5 TM10 NM5 NM10Dimensionless Cauchy effective stress, Σe Dimensionless radial coordinate, ξm=0.6 η=1500 Pa s .TM- t=5 μs TM - t=10 μs FEM - t=5 μs FEM - t=10 μs ξp=1.6800.0040.0080.012 1.6 1.8zoom ξw=1.93 (c) Figure 10: Comparison between theoretical results (TM) and ﬁnite element results (FEM) for m= 0.4. Dimensionless eﬀective Cauchy stress Σrrversus dimensionless radial coordinate ξfor (a)η= 1.5 Pa ·s, (b)η= 150 Pa ·sand (c)η= 1500 Pa ·s. The theoretical results for η= 1.5 Pa ·sandη= 150 Pa ·swere obtained with c0= 0.55andcL= 0.25. In contrast, the theoretical results for η= 1500 Pa ·swere obtained without artiﬁcial viscosity, c0= 0andcL= 0. 24 Acknowledgements The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme. Project PURPOSE, grant agreement 758056. RRwishestoacknowledgethesupportofCNPq,ConselhoNacionaldeDesenvolvimentoCientíﬁcoeTecnológico, grant number 306058/2018-9.