6 In this paper, we have developed a linear stability analysis to predict the formation of necking instabilities 7 in porous ductile plates subjected to dynamic biaxial stretching. The mechanical behavior of the material is 8 described with the Gurson-Tvergaard-Needleman constitutive relation for progressively cavitating solids (Gurson, 9 1977; Tvergaard, 1981, 1982; Tvergaard and Needleman, 1984) which considers the voids to be spherical and 10 the matrix material isotropic with yielding defined by the von Mises (1928) criterion. The analytical model is 11 formulated in a two-dimensional framework in which the multiaxial stress state that develops inside the necked 12 region is approximated with the Bridgman (1952) correction factor, superimposing a hydrostatic stress state to the 13 uniform stress field that develops in the plate before localization. As opposed to the linear stability models published 14 so far to model dynamic necking in ductile plates, which consider the material to be fully dense and incompressible, 15 the approach developed in this paper provides new insights into the interplay between porosity and inertia on plastic 16 localization. In addition, the predictions of the theoretical model for the critical strain leading to necking formation 17 have been compared with unit-cell finite element calculations performed in ABAQUS/Explicit (2019). Satisfactory 18 quantitative and qualitative agreement has been found between the theoretical and computational approach for 19 loading paths ranging from plane strain tension to nearly equibiaxial tension, loading rates varying from 100 s−120 to 10000 s−1, and different values of the initial void volume fraction ranging from 0 .01 to 0 .1. Both analytical 21 and finite element results suggest that the influence of porosity on necking localization increases, due to early 22 voids coalescence, as the loading rate increases and the loading path approaches equibiaxial tension. The original 23 formulation developed in this paper serves as a basis for analytically modeling the dynamic formability of porous 24 ductile plates, and it can be readily extended to consider more complex porous plasticity theories, e.g. constitutive 25 models which consider the anisotropy of the material (Benzerga and Besson, 2001) and/or voids with different 26 shapes (Gologanu et al., 1993; Monchiet et al., 2008). 27 Keywords: 28 ∗Corresponding author. E-mail address: knsouglo@ing.uc3m.es Preprint submitted to European Journal of Mechanics A/Solids July 3, 2021 2 Dynamic formability, Porous plasticity, Linear stability analysis, Necking instabilities 29 1. Introduction 30 Linear stability theories have been widely accepted as an effective approach for identifying the mechanisms 31 which control necking localization in metallic materials subjected to dynamic loading. This technique is based 32 on the superposition of a small perturbation on the homogeneous solution of the problem considered, so that if 33 the perturbation grows faster than the background solution, the plastic flow is unstable and a non-homogeneous, 34 neck-like deformation field can develop. 35 36 Chou and Carleone (1977) formulated a one-dimensional perturbation analysis to investigate the stability of 37 plastic shaped-charge jets subjected to dynamic loading and thus provide insights into the mechanisms which 38 control the formation of multiple necks and fragments in the jets at large strains. The jet was modelled as a 39 slender rod with circular cross-section deforming at constant stretch rate. Moreover, the material was considered 40 to be incompressible, rate-independent and with non-linear strain hardening. Consistent with the Consid` ere (1885) 41 criterion, it was shown that the perturbation was stable –i.e. the necking condition was not met– if the stress 42 was smaller than the slope of the stress-strain curve. On the other hand, it was shown that for unstable jets, the 43 smaller the flow stress, and the larger the material density and the stretching rate, the slower the growth rate 44 of the perturbation, and so the development of the necking pattern. Moreover, the one-dimensional analysis did 45 not yield a finite critical wavelength which causes the fastest growth, so that the perturbation growth increased 46 monotonically with the decrease of the perturbation wavelength. Chou and Carleone (1977) pointed out that it is 47 necessary to account for the multiaxial stress state which develops in a necked section to demonstrate the existence 48 of a critical wavelength for which the perturbation growth is the fastest. This critical perturbation wavelength 49 –also referred to as critical necking wavelength– was presumed to determine the average spacing between necks 50 in the multiple localization pattern that forms in the shaped-charge jets. Shortly after, Walsh (1984) extended 51 the stability analysis of Chou and Carleone (1977) considering Bridgman (1952) approximation to model within a 52 one-dimensional framework the stress enhancement that occurs in a necked section. The material was considered 53 to be incompressible with constant yield strength. According to Chou and Carleone (1977) predictions, it was 54 shown that the local hydrostatic stress that develops in the neck slows down the growth of short wavelength 55 perturbations, while inertia refrains the growth of long wavelength perturbations, leading to the promotion of a 56 critical wavelength for which the growth rate of the disturbance is maximum (and thus the mechanical energy 57 3 required to trigger a neck is minimum). In addition, Walsh (1984) applied dimensional considerations to the 58 perturbed stretching jet and showed that the joint effects of strain rate, material density, flow stress and cross- 59 section radius on the perturbation growth –as earlier identified by Chou and Carleone (1977)– can be collected into 60 a single dimensionless parameter which represents the inertial resistance to motion. A later study of Fressengeas 61 and Molinari (1985) on the stability of slender bars modelled with a thermo-viscoplastic constitutive relation and 62 submitted to uniaxial tension demonstrated that inertia is responsible for the increased ductility that metallic 63 materials generally display at high strain rates (e.g. see Wood (1967) and Hu and Daehn (1996)), while thermal 64 softening was shown to favour early necking formation. Consistently with Walsh (1984) analysis, Fressengeas and 65 Molinari (1985) pointed out that consideration of multiaxial stress effects in the necking zone using Bridgman 66 (1952) correction suppresses short wavelengths and delays plastic localization. 67 68 Following the works of Walsh (1984) and Fressengeas and Molinari (1985), one-dimensional stability models with 69 Bridgman (1952) correction for the axial stress –which stand out for their simplicity– have been used by several 70 authors to investigate the effect of different aspects of the mechanical behaviour of ductile metals on necking 71 formation. For instance, Zhou et al. (2006) modelled the material behaviour using von Mises (1928) plasticity 72 with yield stress dependent on strain and strain rate, and showed that the strain rate sensitivity slows down the 73 growth of perturbations, so that the necks grow more slowly, which delays plastic localization. Finite difference 74 simulations of a round bar subjected to constant stretching velocity confirmed the predictions of the stability 75 analysis such that increasing the strain rate sensitivity increased the material ductility. Zaera et al. (2014) and 76 Rodr´ ıguez-Mart´ ınez et al. (2015) investigated the effect of strain induced martensitic transformation and dynamic 77 recrystallization on the inception of multiple necks in steel and titanium bars subjected to dynamic tensile loading. 78 A key feature of the works of Zaera et al. (2014) and Rodr´ ıguez-Mart´ ınez et al. (2015) was to incorporate into 79 the one-dimensional stability analysis the thermomechanical coupling due to plastic dissipation and the energy 80 release due to the exothermic character of the microstructural transformations. The theoretical predictions were 81 compared with finite element calculations, which confirmed the stabilizing effect of the martensitic transformation 82 on the flow stress, and the destabilizing effect of the dynamic recrystallization, which favours the formation of 83 early necks. Recently, N’souglo et al. (2018) extended the linear stability analysis of Zhou et al. (2006) to consider 84 materials with yielding described with the porous plasticity criterion developed by Gurson (1977). The theoretical 85 predictions for the evolution of the necking strain with strain rate found reasonable agreement with finite element 86 results of long rods subjected to dynamic stretching at strain rates varying from 1000 s−1to 50000 s−1. 87 4 88 Guduru and Freund (2002) derived a three-dimensional stability analysis for rapidly extending cylindrical 89 rods modelled with the constitutive law developed by St¨ oren and Rice (1975) which considers the material to 90 be incompressible, rate-independent with strain hardening. The results obtained for the critical perturbation 91 wavelength were compared with the average neck spacing (number of necks) obtained in the ring expansion tests 92 performed by Grady and Benson (1983) on 1100-0 aluminum and OFHC copper specimens, and satisfactory 93 agreement was obtained for strain rates of the order of 104s−1. The stability analysis of Guduru and Freund 94 (2002) also predicted that the number of necks and the necking strain increase non-linearly with the applied 95 strain rate, consistently with the experimental data obtained by Altynova et al. (1996) for 6061 aluminum and 96 OFHC copper rings expanded at velocities ranging from 50 to 300 m /s. Shortly after, Mercier and Molinari (2003) 97 formulated a three-dimensional linear stability analysis for a cylindrical bar subjected to dynamic extension. As in 98 Guduru and Freund (2002), the material behaviour was modelled with the constitutive relation of St¨ oren and Rice 99 (1975). Mercier and Molinari (2003) showed that inertia slows down the growth of long wavelength perturbations 100 while multidimensional effects (stress multiaxiality) and strain hardening extinct short wavelength perturbations, 101 so that the critical perturbation wavelength depends on the specimen deformation, with the average neck spacing 102 predicted by the stability analysis decreasing with the increase of the strain hardening and the stretch in the bar. 103 104 Linear stability models have also been developed to study the formation of dynamic necks in plates and shells. 105 For instance, Fressengeas and Molinari (1994) extended the analysis performed by Hutchinson et al. (1978) for 106 a plate subjected to plane strain tension to consider inertia effects. A non-linear viscous material was employed 107 in the analysis. Fressengeas and Molinari (1994) showed that viscosity preferentially slows down the growth of 108 short wavelength perturbations –consistent with the results obtained by Zhou et al. (2006) using a 1D-stability 109 analysis, see discussion above– and compared the theoretical predictions for the critical wavelength with the average 110 fragments size obtained in experiments and non-linear numerical simulations reported by Karpp and Simon (1976), 111 so that satisfactory agreement was obtained for copper shaped-charged jets deforming at high strain rates of the 112 order of 8000 s−1. Moreover, Shenoy and Freund (1999) generalized the quasi-static analysis of the plane tension 113 test developed by Hill and Hutchinson (1975) to include dynamic effects. The critical perturbation wavelength 114 was shown to decrease with increasing strain rate, so that the number of necks predicted by the theoretical 115 model increased with the applied loading velocity, in qualitative agreement with the ring expansion experiments 116 of Grady and Benson (1983) and Altynova et al. (1996). Years later, Mercier et al. (2010) developed a linear 117 5 stability analysis to model necking localization in experiments performed with tantalum and copper hemispherical 118 shells subjected to dynamic expansion at strain rates varying between 103s−1and 104s−1. The motion of the 119 hemispheres was not spherically divergent, resulting in the formation of multiple necking bands parallel to each 120 other in an axisymmetric layer located near the equatorial plane where the material stretched under plane strain 121 conditions. The plane strain model of Mercier et al. (2010) was an extension of earlier formulations developed 122 by Fressengeas and Molinari (1994) and Mercier and Molinari (2003) to consider von Mises materials with flow 123 stress dependent on strain, strain rate and temperature. Based on a so-called effective strain rate sensitivity 124 parameter, which merged the thermal and the visco-plastic material characteristics, Mercier et al. (2010) obtained 125 a criterion to determine using the stability analysis the critical strain leading to necking formation and the average 126 spacing between necks. The theoretical predictions were in agreement with the experiments for both tantalum 127 and copper specimens. Shortly after, Jouve (2013, 2015) formulated three-dimensional stability models to predict 128 the formation of necking instabilities in plates of arbitrary thickness subjected to plane strain tension and biaxial 129 tension, respectively. The mechanical behaviour of the material was modelled with von Mises (1928) plasticity, 130 associated flow rule and flow strength dependent on plastic strain, plastic strain rate and temperature. Adiabatic 131 conditions of deformation were considered. The analysis was performed for moderate strain rates within the range 132 10 s−1−20 s−1. Consistent with Fressengeas and Molinari (1994), Jouve (2013, 2015) showed that viscosity 133 essentially refrains the growth of the shortest wavelength instabilities, and identified a loading-path-dependent 134 critical value of material viscosity above which inertia effects on necking formation are negligible (for the range of 135 strain rates considered). Moreover, Zaera et al. (2015) developed a 2D linear stability analysis to model necking 136 formation in thin-plates plates subjected to in-plane biaxial loading for strain rates of the order of 1000 s−1, i.e. 137 above the range of strain rates considered by Jouve (2013, 2015). The model of Zaera et al. (2015) was essentially 138 an extension of earlier formulation of Dudzinski and Molinari (1991) to account for inertia effects, considering 139 Bridgman (1952) correction for the axial stress (in-plane normal stress perpendicular to the necking band). As 140 in Jouve (2013, 2015), the material was modelled with von Mises (1928) plasticity, and a thermo-viscoplastic 141 constitutive relation for the evolution of the yield stress. The predictions of the stability analysis were used to 142 construct forming limit diagrams that were in qualitative agreement with results obtained from finite element 143 calculations, and showed that the critical perturbation wavelength increases as the loading path moves away from 144 uniaxial tension and approaches equibiaxial tension, and so does, therefore, the spacing between necks in the 145 multiple necking patterns predicted in the numerical calculations (see also Rodr´ ıguez-Mart´ ınez et al. (2017)). 146 147 6 The stability analysis formulation of Zaera et al. (2015) has been later extended to consider different plasticity 148 theories with the aim of studying the effect of various aspects of the mechanical behaviour of ductile metals on 149 dynamic necking localization in thin plates. For instance, N’souglo et al. (2020) employed the isotropic form of 150 the Cazacu et al. (2006) yield criterion, which accounts for the effect of tension-compression asymmetry in flow 151 stresses, and showed that the necking strain is significantly greater for a material that displays a larger flow stress 152 in uniaxial compression than in uniaxial tension. The stability analysis calculations performed for plates subjected 153 to dynamic plane strain tension, at strain rates ranging from 4000 s−1to 80000 s−1, were compared with unit-cell 154 finite element simulations assuming that a perturbation mode turns into a necking mode when a critical value of the 155 perturbation growth rate is met. Satisfactory qualitative and quantitative agreement between stability analysis and 156 finite element results was obtained for all the perturbation wavelengths and strain rates investigated. Moreover, 157 N’souglo et al. (2021) considered plastic orthotropic materials modelled with the yield criterion developed by Hill 158 (1948) and performed the stability analysis for loading paths ranging from plane strain tension to equibiaxial 159 tension, and for strain rates varying between 100 s−1and 50000 s−1. Assuming that a neck is formed when the 160 perturbation growth rate reaches a critical value –as in N’souglo et al. (2020)–, the stability analysis was employed 161 to construct forming limit diagrams that found good correlation with the results obtained from the non-linear 162 two zone model developed by Jacques (2020) and with unit-cell finite element calculations. N’souglo et al. (2021) 163 concluded that the effect of anisotropy on dynamic formability increases as the loading path moves away from 164 plane strain and approaches equibiaxial tension, but it is not significantly affected by the prescribed strain rate. 165 Very recently, Rodr´ ıguez-Mart´ ınez et al. (2021) developed a linear stability analysis for plates subjected to plane 166 strain tension modelled with the isotropic Drucker (1949) yield criterion, which depends on both the second and 167 third invariant of the stress deviator. The analysis was performed for effective strain rates ranging from 10−4s−1to 168 8·104s−1, and the theoretical results were compared with uni-cell finite element calculations, as in N’souglo et al. 169 (2020, 2021). Both stability analysis and finite element results showed that for quasi-static loadings the specimen 170 elongation when necking occurs is virtually insensitive to the third invariant of the stress deviator. In contrast, for 171 large values of the prescribed effective strain rate, the influence of the third invariant becomes important, notably 172 for long wavelengths, so that the specimen elongation when necking occurs increases as the ratio between the yield 173 stresses in uniaxial tension and in pure shear increases, while the necking strain decreases. 174 175 However, to the authors’ knowledge, all the linear stability models developed so far to study the formation 176 of dynamic necks in plates and shells –most of them cited in previous two paragraphs– consider the material to 177 7 be fully dense and incompressible. As of today, there is no specific formulation for porous materials. In order 178 to fill this gap, in this paper we have developed a two-dimensional perturbation analysis for thin plates modelled 179 with the Gurson-Tvergaard-Needleman constitutive relation (Gurson, 1977; Tvergaard, 1981, 1982; Tvergaard and 180 Needleman, 1984) and subjected to in-plane dynamic loading. The stress multiaxiality effects in the necked region 181 have been approximated using Bridgman (1952) correction. Moreover, we have investigated loading paths ranging 182 from plane strain tension up to equibiaxial tension, and strain rates varying from 100 s−1to 10000 s−1. An original 183 contribution of this work is that we have enhanced the calibration procedure developed by N’souglo et al. (2020) 184 for the stability analysis, so that in this paper a perturbation mode is considered to turn into a necking mode 185 based on a function of the perturbation growth rate which takes into account the combined effect that loading 186 path and loading rate have on the necking strain. The stability analysis predictions obtained with this novel 187 necking criterion have been compared with unit-cell finite element calculations, and satisfactory agreement has 188 been obtained for all the strain rates and loading paths considered. An important outcome of this work is that 189 both analytical and finite element results suggest that the influence of porosity on necking localization increases, 190 due to early voids coalescence, as the loading rate increases and the loading path approaches equibiaxial tension. 191 2. Analytical model 192 In this section, we develop a 2D perturbation analysis to model the formation of necking instabilities in porous 193 plates subjected to dynamic biaxial stretching. We consider a thin plate of initial thickness h0and edges of initial 194 length L0subjected to constant in-plane velocities, VxandVy, on opposed sides, see Fig. 1. Cartesian coordinates 195 (X, Y, Z ) define the location of a material point in the undeformed state, while the location of the same material 196 point in the deformed configuration is ( x, y, z ). The origin of both Lagrangian and Eulerian coordinate systems 197 is located at the center of mass of the plate. Note that the out-of-plane direction Zis assumed of plane stress. 198 The elastic strains are neglected, and the mechanical behavior of the plate is described with the porous plasticity 199 constitutive relation developed by Gurson (1977), Tvergaard (1981, 1982) and Tvergaard and Needleman (1984) 200 which accounts for void growth and coalescence. The multiaxial stress state that develops inside a necked region is 201 modeled with Bridgman (1952) approximation superimposing a hydrostatic stress state to the uniform stress field 202 that develops in the plate during homogeneous deformation (before necking localization, see the works of Zhou 203 et al. (2006), Zaera et al. (2015) and Xavier et al. (2021), among others). 204 2.1. Governing equations 205 The strain rates in the loading directions are: 206 8 Figure 1: Schematic representation of the geometry and boundary conditions of the problem addressed: a thin plate of initial thickness h0and edges of initial length L0is subjected to constant and opposed stretching velocities, VxandVy, on opposed sides. The Lagrangian and Eulerian Cartesian coordinate systems, ( X, Y, Z ) and ( x, y, z ), respectively, with origin located at the center of mass of the plate, are associated to the applied velocity field. dxx=∂Vx ∂x dyy=∂Vy ∂y(1) which are related to the corresponding stretches, FxxandFyy, and stretch rates, ˙Fxxand ˙Fyy, as: 207 Fxx=˙Fxx dxx Fyy=˙Fyy dyy(2) Moreover, the through-thickness stretch is: 208 Fzz=h h0(3) 9 where his the current thickness of the plate. 209 210 The balance of mass provides a relationship between the initial void volume fraction f0, the current void volume 211 fraction f, and the stretches: 212 FxxFyyFzz=1−f0 1−f(4) which leads to the following expression for the density ratio: 213 ρ ρ0=1−f 1−f0 (5) where ρandρ0are the current and initial material densities. 214 215 The Gurson-Tvergaard-Needleman flow potential is expressed as follows (Gurson, 1977; Tvergaard, 1981, 1982; 216 Tvergaard and Needleman, 1984): 217 Φ = σ2 xx+σ2 yy−σxxσyy ¯σ2! + 2q1f∗coshq2(σavg xx+σavg yy+σavg zz) 2¯σ −1−(q1f∗)2= 0 (6) with q1andq2being material parameters, ¯ σthe flow strength of the matrix material (see equation (11)), and f∗218 the effective porosity, defined as: 219 f∗=  f if f < f c fc+(fu−fc) (f−fc) (ff−fc)if fc⩽f⩽ff fu if f > f f(7) where fcis the void volume fraction at which voids coalesce, ffis the void volume fraction at final fracture of the 220 material and fu= 1/q1is the ultimate void volume fraction. 221 222 10 Moreover, σxxandσyyare the uniform stresses in the loading directions (stress field before necking localization), 223 while σavg xx,σavg yyandσavg zzare the average stresses defined as: 224 σavg xx=σxx+σh=B(φ)σxx σavg yy=σyy+σh=σyy+σavg xx−σxx σavg zz=σh=σavg xx−σxx(8) where σhdenotes the hydrostatic stress that develops in a necked section and B(φ) is the correction factor 225 introduced by Bridgman (1952): 226 B(φ) =p 1 +φ−1ln 1 + 2 φ+ 2p φ(1 +φ) −1 (9) with: 227 φ=h 8∂2h ∂x2(10) While Bridgman (1952) correction for necking was developed for elastic-perfectly plastic materials under quasi- 228 static loading conditions, in this work we assume that it can be used to approximate the axial stress in the necked 229 section of porous materials with flow strength dependent on strain and strain rate. Moreover, in future works, we 230 shall consider including in the linear stability analysis the one-dimensional model recently developed by Audoly 231 and Hutchinson (2016, 2019) to determine the stress state inside the neck, and we shall compare the results with 232 the approach used in this paper based on the classical theory of Bridgman (1952). 233 234 The flow strength of the matrix material is defined as: 235 ¯σ=σ0 ε0+ ¯εn˙¯ε ˙ε0m (11) where σ0,ε0,n, ˙ε0andmare material parameters, and ¯ εis the effective plastic strain: 236 11 ¯ε=Zt 0˙¯ε(τ) dτ (12) with ˙¯εbeing the effective plastic strain rate. We are aware that at high strain rates, due to local adiabatic heating, 237 the effect of thermal softening on the mechanical behavior of ductile materials is important. Nevertheless, the 238 main objective of this work is to investigate the effect of material porosity on dynamic necking localization, so 239 that neglecting the temperature sensitivity of the material facilitates the interpretation of results. The reader is 240 referred to Section 6 .4 of the recent paper of N’souglo et al. (2021) to obtain insights into the effect of adiabatic 241 heating on necking localization in anisotropic ductile plates subjected to dynamic biaxial loading for loading paths 242 varying from plane strain to equibiaxial tension. 243 244 Assuming that the rate of macroscopic plastic work is equal to the rate of equivalent plastic work in the matrix 245 material, it follows that: 246 σavg xxdxx+σavg yydyy+σavg zzdzz= (1−f) ¯σ˙¯ε (13) 247 Moreover, considering an associated plastic flow rule, the following relationships between strain rates and 248 stresses are obtained: 249 dxx=˙λ2σxx−σyy ¯σ2+q1q2f∗ ¯σsinhq2(σavg xx+σavg yy+σavg zz) 2¯σ dyy=˙λ2σyy−σxx ¯σ2+q1q2f∗ ¯σsinhq2(σavg xx+σavg yy+σavg zz) 2¯σ dzz=˙λ −σxx+σyy ¯σ2+q1q2f∗ ¯σsinhq2(σavg xx+σavg yy+σavg zz) 2¯σ(14) where ˙λis the plastic flow proportionality factor. 250 251 The evolution of the void volume fraction is defined as: 252 12 ˙f= (1−f)˙λ3q1q2f∗ ¯σsinhq2(σavg xx+σavg yy+σavg zz) 2¯σ (15) where the matrix material has been assumed to be incompressible. 253 254 The momentum balance equations in the loading directions take the form: 255 Fyy∂(hσavg xx) ∂X=ρh0∂Vx ∂t Fxx∂(hσyy) ∂Y=ρh0∂Vy ∂t(16) Note that the momentum balance equation in the out of plane direction has been disregarded, see Zaera et al. 256 (2015) and N’souglo et al. (2021). 257 258 The solution of the problem at any given loading time, before necking localization, is called homogeneous or fun- 259 damental solution S1= (V1 x, V1 y, d1 xx, d1 yy, d1 zz, F1 xx, F1 yy, F1 zz, h1,˙¯ε1,¯ε1,¯σ1,˙λ1, σ1 xx, σ1 yy, σavg, 1 xx, σavg, 1 yy, σavg, 1 zz, ρ1, f1, φ1)T, 260 which is obtained by solving the system of equations (1)-(6), (8) and (10)-(16) under the following initial and 261 boundary conditions: 262 Vx(X, Y, 0) = ˙ ε0 xxX; Vy(X, Y, 0) = ˙ ε0 yyY; ¯ ε(X, Y, 0) = 0 h(X, Y, 0) = h0; f(X, Y, 0) = f0; ρ(X, Y, 0) = ρ0 Vx L0/2, Y, t =−Vx −L0/2, Y, t = ˙ε0 xxL0/2 Vy X, L0/2, t =−Vy X,−L0/2, t = ˙ε0 yyL0/2 where ˙ ε0 xxand ˙ε0 yydefine the initial strain rates in XandYdirections, respectively. Hereinafter, ˙ ε0 xxwill be also 263 referred to as imposed initial major strain rate. Note that the applied loading is selected through input values of 264 ˙ε0 xxandχ= ˙ε0 yy/˙ε0 xx, where the ratio χwill be hereinafter referred to as loading path ( χ= 0 corresponds to plane 265 strain tension, and χ= 1 to equibiaxial tension). 266 13 2.2. Linear stability analysis 267 At some time t1, we superimpose a small perturbation δS(X, t)t1=δS1eiξX+η(t−t1)to the fundamental solution, 268 where δS1= (δVx, δVy, δdxx, δdyy, δdzz, δFxx, δFyy, δFzz, δh, δ ˙¯ε, δ¯ε, δ¯σ, δ˙λ, δσ xx, δσyy, δσavg xx, δσavg yy, δσavg zz, δρ, δf, δφ )T269 is the perturbation amplitude, ξis the wavenumber and ηis the instantaneous growth rate of the perturbation at 270 time t1. While ηandξare considered time independent (frozen coefficients method) for the sake of simplicity, it 271 has to be noted that a time dependent numerical solution for the linear perturbation analysis has been recently 272 developed by Xavier et al. (2020, 2021) for von Mises materials subjected to plane strain tension. Moreover, note 273 that the perturbation is inserted perpendicular to the main loading direction, i.e. on lines X=constant , because 274 this is the preferential orientation for an isotropic material to trigger a neck within the range of loading paths 275 considered 0 ≤χ < 1, see Zaera et al. (2015) and Jacques (2020). The perturbed solution of the problem is 276 S=S1+δS, with |δS| ≪ S1 , so that keeping only the first-order terms in the increments δS1, equations (1)-(6), 277 (8) and (10)-(16) are linearized are follows: 278 ˆF1 xxδˆdxx+ˆd1 xxδˆFxx−iˆξδˆVx= 0 ˆF1 yyδˆdyy+ˆd1 yyδˆFyy= 0(17)  ˆη−ˆd1 xx δˆFxx−ˆF1 xxδˆdxx= 0  ˆη−ˆd1 yy δˆFyy−ˆF1 yyδˆdyy= 0(18) δˆh−δˆFzz= 0 (19) 1 ˆF1xxδˆFxx+1 ˆF1yyδˆFyy+1 ˆF1zzδˆFzz−1 1−ˆf1δˆf= 0 (20) δˆρ+1 1−f0δˆf= 0 (21) 14 ˆQ1δˆσavg xx+ˆQ2δˆσavg yy+ˆQ3δˆσavg zz+ˆQ4δˆσxx+ˆQ5δˆσyy+ˆQ6δˆ¯σ+ˆQ7δˆf= 0 (22) δˆσavg xx−δˆσxx−2 3ˆσ1 xxδˆφ= 0 δˆσavg yy−δˆσyy−δˆσavg xx+δˆσxx= 0 δˆσavg zz−δˆσavg xx+δˆσxx= 0(23) 8 ˆF1 xx2 δˆφ+ˆh1 ˆξ2 δˆh= 0 (24) δˆ¯σ−ˆP1δˆ¯ε−ˆP2δˆ˙¯ε= 0 (25) δˆ˙¯ε−ˆηδˆ¯ε= 0 (26) ˆσ1 xxδˆdxx+ˆd1 xxδˆσavg xx+ ˆσ1 yyδˆdyy+ˆd1 yyδˆσavg yy+ˆd1 zzδˆσavg zz− 1−ˆf1 δˆ¯σ− 1−ˆf1 ˆ¯σ1δˆ˙¯ε+ˆ¯σ1δˆf= 0 (27) δˆdxx−ˆd1 xx ˆ˙λ1δˆ˙λ−ˆM1δˆσavg xx−ˆM2δˆσavg yy−ˆM3δˆσavg zz−ˆM4δˆσxx−ˆM5δˆσyy−ˆM6δˆ¯σ−ˆM7δˆf= 0 δˆdyy−ˆd1 yy ˆ˙λ1δˆ˙λ−ˆR1δˆσavg xx−ˆR2δˆσavg yy−ˆR3δˆσavg zz−ˆR4δˆσxx−ˆR5δˆσyy−ˆR6δˆ¯σ−ˆR7δˆf= 0 δˆdzz−ˆd1 zz ˆ˙λ1δˆ˙λ−ˆS1δˆσavg xx−ˆS2δˆσavg yy−ˆS3δˆσavg zz−ˆS4δˆσxx−ˆS5δˆσyy−ˆS6δˆ¯σ−ˆS7δˆf= 0(28) 15  ˆη−ˆT5 δˆf−ˆ˙f1 ˆ˙λ1δˆ˙λ−ˆT1δˆσavg xx−ˆT2δˆσavg yy−ˆT3δˆσavg zz−ˆT4δˆ¯σ= 0 (29) iˆξˆF1 yyˆσ1 xxδˆh+iˆξˆh1ˆF1 yyδˆσavg xx−1 ˜L2ˆρ1ˆηδˆVx= 0 δˆVy= 0(30) This system of equations has been adimensionalized by introducing the following dimensionless groups: 279 ˆVx=Vx h0˙¯ε1;ˆVy=Vy h0˙¯ε1;ˆdxx=dxx ˙¯ε1;ˆdyy=dyy ˙¯ε1;ˆdzz=dzz ˙¯ε1;ˆFxx=Fxx;ˆFyy=Fyy ˆFzz=Fzz;ˆh=h h0;ˆ˙¯ε=˙¯ε ˙¯ε1;ˆ¯ε= ¯ε;ˆ¯σ=¯σ σ0;ˆ˙λ=˙λ σ0˙¯ε1; ˆσxx=σxx σ0 ˆσyy=σyy σ0; ˆσavg xx=σavg xx σ0; ˆσavg yy=σavg yy σ0; ˆσavg zz=σavg zz σ0; ˆρ=ρ ρ0;ˆf=f; ˆφ=φ ˆη=η ˙¯ε1;ˆξ=ξh0;ˆP1=1 σ0∂¯σ ∂¯ε t1;ˆP2=˙¯ε1 σ0∂¯σ ∂˙¯ε t1;ˆQ1=σ0∂Φ ∂σavg xx t1=ˆQ2=ˆQ3 ˆQ4=σ0∂Φ ∂σxx t1;ˆQ5=σ0∂Φ ∂σyy t1;ˆQ6=σ0∂Φ ∂¯σ t1;ˆQ7=∂Φ ∂f t1 ˆM1=σ0 ˙¯ε1∂dxx ∂σavg xx t1=ˆM2=ˆM3=ˆR1=ˆR2=ˆR3=ˆS1=ˆS2=ˆS3;ˆM4=σ0 ˙¯ε1∂dxx ∂σxx t1=ˆR5 ˆM5=σ0 ˙¯ε1∂dxx ∂σyy t1=ˆR4=ˆS4=ˆS5;ˆM6=σ0 ˙¯ε1∂dxx ∂¯σ t1;ˆM7=1 ˙¯ε1∂dxx ∂f t1=ˆR7=ˆS7 ˆR6=σ0 ˙¯ε1∂dyy ∂¯σ t1;ˆS6=σ0 ˙¯ε1∂dzz ∂¯σ t1;ˆ˙f1=˙f1 ˙¯ε1;ˆT1=σ0 ˙¯ε1∂˙f ∂σavg xx t1=ˆT2=ˆT3 ˆT4=σ0 ˙¯ε1∂˙f ∂¯σ t1;ˆT5=1 ˙¯ε1∂˙f ∂f t1;1 ˜L2=ρ0 h0˙¯ε12 σ0 280 A non-trivial solution for δˆS1is obtained by imposing that the determinant of the matrix of coefficients of 281 the system of linear algebraic equations (17)-(30) is equal to zero, which leads to a fourth order polynomial for 282 ˆηwhich depends on the adimensionalized fundamental solution ˆS1and on the adimensionalized wavenumber ˆξ. 283 The adimensionalized wavenumber is related to the normalized perturbation wavelength as L0/h0= 2πˆξ−1(see 284 N’souglo et al. (2020, 2021)). Note that the normalized perturbation wavelength is usually referred to as necking 285 wavelength (N’souglo et al., 2020; Rodr´ ıguez-Mart´ ınez et al., 2021). The root of the polynomial with physical 286 16 meaning is the one with the greatest positive real part (Dudzinski and Molinari, 1991), hereinafter denoted by 287 ˆη+. According to the notation used by El Ma¨ ı et al. (2014) and Rodr´ ıguez-Mart´ ınez et al. (2021), among others, 288 ˆη+will be referred to as instantaneous instability index, and it will be used to compute the so-called cumulative 289 instability index I=Rt 0ˆη+˙¯ε1dτ, which tracks the history of the instantaneous grow rate of all the growing modes. 290 The cumulative instability index generally provides better predictions than the instantaneous instability index for 291 the critical strain leading to necking formation at high strain rates (see the recent paper of Rodr´ ıguez-Mart´ ınez 292 et al. (2021)). Assuming that a perturbation mode turns into a necking mode when the cumulative instability index 293 reaches a critical value which depends on the loading path and the strain rate, see Section 3.1, the linear stability 294 analysis is used in Section 3 to construct forming limit diagrams which are compared with unit-cell finite element 295 simulations performed in ABAQUS/Explicit (2019) to demonstrate the ability of the analytical model to capture 296 the effect of porosity on necking localization (this is same methodology used by N’souglo et al. (2021) to assess 297 the ability of the linear stability analysis to capture the effect of anisotropy on necking localization in materials 298 described with Hill (1948) yield criterion). Details of the unit-cell finite element model –which is not shown here 299 for the sake of brevity– can be found in Sections 4 .1 and 5 of Rodr´ ıguez-Mart´ ınez et al. (2017) and N’souglo et al. 300 (2021), respectively. In the finite element calculations, the flow strength of the matrix material, equation (11), has 301 been used along with the Gurson-Tvergaard-Needleman model pre-implemented in ABAQUS/Explicit through a 302 user-defined subroutine VUHARD (ABAQUS/Explicit, 2019). 303 304 The linear stability analysis has been implemented in the software Wolfram Mathematica. In Supplementary 305 material , we include the links to download all the codes developed, which come with specific instructions for 306 reproducing the stability analysis results reported in Section 3. 307 3. Results 308 This section of the paper is split into three parts. In Section 3.1, we show the procedure for calibration of 309 the linear stability analysis using unit-cell finite element calculations. Moreover, the effect of initial void volume 310 fraction on the critical strain leading to necking formation is investigated in Section 3.2, while the effect of loading 311 rate is explored in Section 3.3. Results are presented for loading paths ranging from plane strain tension ( χ= 0) 312 to nearly equibiaxial tension, for imposed initial major strain rates varying between 100 s−1and 10000 s−1, and 313 for different initial void volume fractions spanning from 0 .01 to 0 .1. For all the calculations performed, the initial 314 thickness of the plate is h0= 2 mm. In the unit-cell finite element simulations, the normalized imperfection 315 17 amplitude is ∆ = 0 .2% (see Section 5 of N’souglo et al. (2021)). The baseline numerical values of the material 316 parameters employed in the stability analysis and in the finite element simulations are given in Table 1. We use 317 standard values representative of general purpose steel. Note that in some calculations presented in Sections 3.2 318 and 3.3 we disregard voids coalescence and fracture considering that fc→ ∞ andff→ ∞ . Note also that, 319 unlike the stability analysis, the finite element calculations performed for comparison with the analytical model 320 predictions consider the elastic response of the material, which is assumed to be isotropic and linear. 321 Symbol Property and units Value ρ0Initial density (kg /m3) 7850 G Elastic shear modulus (GPa) 75.19 K Bulk modulus (GPa) 196.1 q1 Material parameter, Eq. (6) 1.25 q2 Material parameter, Eq. (6) 1 fc Material parameter, Eq. (7) red 0.12 ff Material parameter, Eq. (7) 0.5 σ0Material parameter (MPa), Eq. (11) 300 ε0Material parameter, Eq. (11) 0.023 n Material parameter, Eq. (11) 0.12 m Material parameter, Eq. (11) 0.01 ˙ε0Material parameter (s−1), Eq. (11) 1000 Table 1: Baseline numerical values of the material parameters employed in the stability analysis and in the finite element simulations. The values of the parameters corresponding to the flow strength of the matrix material are taken from N’souglo et al. (2021). 3.1. Calibration of the linear stability analysis 322 We have enhanced the procedure developed by N’souglo et al. (2020) to calibrate the linear stability analysis 323 so that in this paper a perturbation mode is assumed to turn into a necking mode when the cumulative instability 324 index reaches a critical value which depends on the imposed initial major strain rate and on the loading path 325 (recall that N’souglo et al. (2020, 2021) assumed Icto be constant): 326 Ic ˙ε0 xx, χ =a+bχ+c˙ε0 xx+dχ2+e( ˙ε0 xx)2(31) where a,b,c,dandeare calibration coefficients to be determined. For that purpose, we have devised a two-steps 327 methodology in which the results of the stability analysis for the evolution of the major necking strain εneck xxwith 328 the necking wavelength L0/h0are fitted to unit-cell calculations for different cell lengths in which the material 329 behaviour is modelled with von Mises plasticity (i.e. imposing f0= 0 on the Gurson-Tvergaard-Needleman 330 18 constitutive relation): 331 1. For plane strain tension χ= 0 and initial major strain rates ˙ ε0 xx= 1000, 5000 and 10000 s−1, the best fitting 332 between stability analysis and finite element results is obtained for Ic= 3.8, 4.33 and 4 .9, respectively, so 333 thata= 3.65722, c= 0.000144833 s and e=−2.05556 ·10−9s2(see Fig. 2(a)). 334 335 2. For imposed initial major strain rate ˙ ε0 xx= 10000 s−1and loading paths χ= 0.5 and 0 .75, the best fitting 336 between stability analysis and finite element results is obtained for Ic= 5.78 and 6 .19, respectively, so that 337 b= 1.84 and d=−0.16 (see Fig. 2(b)). 338 The major necking strain is the logarithmic strain in the major loading direction when necking occurs and 339 it is determined both in the stability analysis and in the finite element simulations as in N’souglo et al. (2021), 340 so that the reader is referred to Section 6 .1 therein for the procedure to obtain εneck xx, and no additional details 341 are provided here for the sake of brevity. The εneck xx−L0/h0curves display a nonlinear concave-upward shape, so 342 that the major necking strain first decreases, reaches a minimum for an intermediate value of L0/h0, and then 343 increases (the increase being less as the loading rate decreases, so that for 100 s−1the minimum of the εneck xx−L0/h0344 curve is shallow). The larger values of the major necking strain for small and large values of L0/h0are due to 345 the stabilizing effects of stress multiaxiality and inertia, respectively (Fressengeas and Molinari, 1994; Mercier and 346 Molinari, 2003). The minimum value of εneck xx, also referred to as critical major necking strain and hereinafter 347 denoted by εc xx, determines the size of the necks L0/h0cwhich form earlier, i.e. the necks which require less 348 energy to develop (N’souglo et al., 2020). Note that the critical major necking strain εc xxincreases as the loading 349 path moves away from plane strain and approaches equibiaxial tension, and the value of L0/h0cdecreases with 350 the imposed initial major strain rate. Moreover, following the procedure detailed in Section 6 .1 of N’souglo et al. 351 (2021), the critical major necking strain εc xxand the corresponding critical minor necking strain εc yy(i.e. the 352 logarithmic strain in the minor loading direction when necking occurs, i.e. when εc xxis measured), obtained from 353 the finite element calculations and the linear stability analysis for different loading paths and loading rates, are 354 used in the following sections of this paper to construct forming limit diagrams. Namely, the predictions of the 355 linear stability analysis obtained with the calibrated necking criterion (31) will be compared with finite element 356 calculations in which the material is modeled with the Gurson-Tvergaard-Needleman constitutive relation (i.e. 357 f0̸= 0). Note that, while the results of the stability analysis for the von Mises material are fitted (see Fig. 2), the 358 results for the porous material will be predicted. 359 19 00.20.40.60.811.21.4 0 1 2 3 4 5 6FEM: ε0 xx=1000 s-1 FEM: ε0 xx=5000 s-1 FEM: ε0 xx=10000 s-1 LSA: ε0 xx=1000 s-1 (Ic=3.8) LSA: ε0 xx=5000 s-1 (Ic=4.33) LSA: ε0 xx=10000 s-1 (Ic=4.9) L0/h0χ=0 f0=0 (von Mises)Major necking strain, εxxneck (L0/h0)cεxxcFEM: ε0 xx=1000 s-1 FEM: ε0 xx=5000 s-1 FEM: ε0 xx=10000 s-1 LSA: ε0 xx=1000 s-1 (Ic=3.8) LSA: ε0 xx=5000 s-1 (Ic=4.33) LSA: ε0 xx=10000 s-1 (Ic=4.9) . . . . . . (a) 00.40.81.21.62 0 1 2 3 4 5 6FEM: χ=0.5 FEM: χ=0.75 LSA: χ=0.5 (Ic=5.78) LSA: χ=0.75 (Ic=6.19) L0/h0f0=0 (von Mises)Major necking strain, εxx ε0 =10000 s-1 xx.neck εxxc (L0/h0)cFEM: χ=0.5 FEM: χ=0.75 LSA: χ=0.5 (Ic=5.78) LSA: χ=0.75 (Ic=6.19) (b) Figure 2: Calibration of the linear stability analysis (LSA) with the finite element calculations (FEM). Major necking strain εneck xx versus L0/h0(dimensionless cell length and necking wavelength in the finite element calculations and in the linear stability analysis, respectively), for f0= 0 (von Mises). (a) Results are shown for χ= 0 (plane strain tension) and three different imposed initial major strain rates ˙ ε0 xx= 1000, 5000 and 10000 s−1. (b) Results are shown for imposed initial major strain rate ˙ ε0 xx= 10000 s−1and two different loading paths χ= 0.5 and 0 .75. 3.2. The effect of porosity 360 The calculations presented in this section are performed for two different imposed initial major strain rates 361 ˙ε0 xx= 500 s−1and 10000 s−1, for which the value of the inertia parameter is ˜L−1≈0.005 and ≈0.1, respectively. 362 This range of strain rates lies within the loading conditions attained in dynamic forming processes (e.g. see Dariani 363 et al. (2009), Golovashchenko et al. (2013) and Li et al. (2017), among others). Moreover, the numerical results 364 reported by N’souglo et al. (2020) and Zheng et al. (2020) indicated that the effect of inertia on neck retardation 365 starts to be important for ˜L−1in the range 0 .06−0.1. Initial void volume fractions ranging from 0 .01 to 0 .1 366 are considered. This range of values of f0is chosen because sintered and additively manufactured metals contain 367 residual porosities that can be as large as 10% (e.g. see Aboulkhair et al. (2014), Heidari et al. (2020) and Marvi- 368 Mashhadi et al. (2021)). In addition, exploring a wide range of initial void volume fractions helps to enlighten the 369 effect of porosity on necking formation for different loading paths. We have also performed calculations for f0= 0, 370 which corresponds to von Mises (1928) plasticity (fully dense material). 371 372 Fig. 3 compares forming limit diagrams (critical major necking strain εc xxversus critical minor necking strain 373 εc yy) predicted by the stability analysis with finite element calculations for two values of the initial void volume 374 fraction: f0= 0 (von Mises) and f0= 0.1. For the porous material, void coalescence and fracture have been 375 disregarded, i.e. f∗=fconsidering fc→ ∞ andff→ ∞ , see equation (7). The imposed initial major strain rate 376 20 is ˙ε0 xx= 10000 s−1. There is a monotonic increase of εc xxas the value of χincreases so that the minimum critical 377 major necking strain corresponds to plane strain tension ( χ= 0). This trend is consistent with experimental, 378 numerical and analytical results published elsewhere (Kuroda and Tvergaard, 2000; Brunet and Morestin, 2001; 379 Banabic et al., 2010). Moreover, as compared to the results obtained with the von Mises material, increasing the 380 initial porosity up to 0 .1 shifts the εc xx−εc yycurve downwards. The drop of εc xxis maximum for χ= 0 (≈15%), 381 and decreases as the loading path approaches equibiaxial tension. Moreover, note that the agreement between the 382 theoretical predictions and the numerical simulations is excellent for the whole range of loading paths considered. 383 In this regard, recall that the calibration of the stability analysis, equation (31), was performed using the numerical 384 results obtained for ˙ ε0 xx= 10000 s−1for the von Mises material, which explains the perfect matching between finite 385 element calculations and stability analysis for f0= 0. On the other hand, the results of the analytical model for 386 f0= 0.1 are fully predictive , i.e. none of the finite element results for this initial void volume fraction was used for 387 the calibration of the stability analysis, so that the agreement between both approaches makes apparent the ability 388 of the linear perturbation model to capture the effect of porosity on the necking formability of porous metals. 389 00.20.40.60.811.2 0 0.2 0.4 0.6 0.8 1 1.2FEM: f0=0 (von Mises) FEM: f0=0.1 LSA: f0=0 (von Mises) LSA: f0=0.1 Critical minor necking strain, εyyCritical major necking strain, εxxc cε xx =10000 s-1 0. WITHOUT COALESCENCE ε xx= ε yyc cFEM: f0=0 (von Mises) FEM: f0=0.1 LSA: f0=0 (von Mises) LSA: f0=0.1 Figure 3: Comparison between finite element results (FEM) and linear stability analysis predictions (LSA). Forming limit diagrams, critical major necking strain εc xxversus critical minor necking strain εc yy, for loading conditions ranging from plane strain stretching (χ= 0) to nearly equibiaxial tension. Results are shown for initial void volume fraction f0= 0 (von Mises) and f0= 0.1 disregarding void coalescence, i.e. f∗=f, see equation (7). The imposed initial major strain rate is ˙ ε0 xx= 10000 s−1. 390 21 Fig. 4 illustrates the influence of void coalescence on forming limit diagrams calculated using finite element 391 simulations and predicted with the stability analysis for three different values of the initial void volume fraction 392 f0= 0.01, 0 .02 and 0 .03. Namely, Fig. 4a shows results for fc→ ∞ andff→ ∞ (i.e. f∗=f), so that 393 the constitutive model accounts for the growth of voids while neglecting coalescence and fracture, and Fig. 4b 394 considers that fc= 0.12 and ff= 0.5 (see Table 1). We have checked that in all the calculations included in 395 Fig. 4b necking occurs before fracture (recall that the stability analysis only models necking). As in Fig. 3, the 396 imposed initial major strain rate is ˙ ε0 xx= 10000 s−1. The quantitative agreement between the predictions of the 397 analytical model and the numerical simulations is excellent for all the loading paths investigated, and for the three 398 values of the initial void volume fraction considered. Note that, if coalescence and fracture are disregarded, Fig. 399 4a, the influence of the initial porosity on the forming limit diagrams, for the three values of f0considered, is 400 small. For plane strain tension the critical major necking strain εc xxis slightly smaller as the initial void volume 401 fraction increases, while the trend is reversed near equibiaxial tension. In contrast, accounting for voids coalescence 402 boosts the effect of initial void volume fraction on the forming limit diagrams, see Fig. 4b. Specifically, the critical 403 major necking strain εc xxdecreases with the value of f0for all the loading paths investigated, and the difference 404 between the results obtained for different values of the initial porosity increases as the loading path approaches 405 equibiaxial tension. Note that, as the initial void volume fraction increases, less plastic deformation is needed to 406 meet the coalescence criterion, leading to earlier drop of the flow stress, which favors necking formation, shifting 407 downwards the εc xx−εc yycurve. For instance, the stability analysis predicts that coalescence occurs for εc yy= 0.35, 408 0.19 and 0 .07 in the cases of f0= 0.01, 0.02 and 0 .03, respectively, while the values obtained with the finite element 409 calculations for the same initial void volume fractions are εc yy= 0.35, 0.19 and 0 .08. 410 Fig. 5 shows forming limit diagrams for lower imposed initial major strain rate ˙ ε0 xx= 500 s−1(20 times 411 lower than in Fig. 4). The results are qualitatively the same as those obtained for 10000 s−1, and the agreement 412 between finite element simulations and analytical predictions is excellent for all the loading paths. If coalescence 413 and fracture are disregarded, fc→ ∞ andff→ ∞ , the results obtained for f0= 0.01, 0 .02 and 0 .03 are very 414 similar, see Fig. 5a. Near plane strain tension the value of εc xxis slightly smaller as f0increases, while this trend 415 is reversed as the loading path approaches equibiaxial tension, as in the case of ˙ ε0 xx= 10000 s−1shown in Fig. 416 4a. On the other hand, including coalescence and fracture in the calculations, fc= 0.12 and ff= 0.5, see Fig. 417 5b, yields differences in the results obtained for the three values of f0considered which increase as the loading 418 path moves away from plane strain, as in the case of ˙ ε0 xx= 10000 s−1shown in Fig. 4b. Moreover, note that 419 the main quantitative difference between the forming limit diagrams obtained for 10000 s−1and 500 s−1is that, 420 22 00.20.40.60.811.2 0 0.2 0.4 0.6 0.8 1 1.2FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 Critical minor necking strain, εyyCritical major necking strain, εxxc cε xx= ε yyc cWITHOUT COALESCENCEε xx =10000 s-1 0. FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 (a) 00.20.40.60.811.2 0 0.2 0.4 0.6 0.8 1 1.2FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 Critical minor necking strain, εyyCritical major necking strain, εxxc cWITH COALESCENCEε xx =10000 s-1 0. ε xx= ε yyc cFEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 (b) Figure 4: Comparison between finite element results (FEM) and linear stability analysis predictions (LSA). Forming limit diagrams, critical major necking strain εc xxversus critical minor necking strain εc yy, for loading conditions ranging from plane strain stretching (χ= 0) to nearly equibiaxial tension. Results are shown for three values of the initial void volume fraction f0= 0.01, 0.02 and 0 .03. The imposed initial major strain rate is ˙ ε0 xx= 10000 s−1. Results corresponding to: (a) fc→ ∞ andff→ ∞ (i.e.f∗=f), so that the constitutive model accounts for the growth of voids while neglecting coalescence and fracture, and (b) fc= 0.12 and ff= 0.5 (see Table 1). 23 independently of the value of the initial void volume fraction, increasing the loading rate shifts the εc xx−εc yycurves 421 upwards due to the stabilizing effect of inertia which delays necking formation (we further elaborate on this issue 422 in Section 3.3). 423 These results highlight the influence of void coalescence in the formability limits of porous metals, and reinforce 424 the idea that the stability analysis captures the mechanisms controlling the formation of dynamic necks for different 425 loading paths, loading rates and void volume fractions. 426 00.20.40.60.811.2 0 0.2 0.4 0.6 0.8 1 1.2FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 Critical minor necking strain, εyyCritical major necking strain, εxxc cε xx= ε yyc cWITHOUT COALESCENCEε xx =500 s-1 0. FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 (a) 00.20.40.60.811.2 0 0.2 0.4 0.6 0.8 1 1.2FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 Critical minor necking strain, εyyCritical major necking strain, εxxc cε xx= ε yyc cWITH COALESCENCEε xx =500 s-1 0. FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 (b) Figure 5: Comparison between finite element results (FEM) and linear stability analysis predictions (LSA). Forming limit diagrams, critical major necking strain εc xxversus critical minor necking strain εc yy, for loading conditions ranging from plane strain stretching (χ= 0) to nearly equibiaxial tension. Results are shown for three values of the initial void volume fraction f0= 0.01, 0.02 and 0 .03. The imposed initial major strain rate is ˙ ε0 xx= 500 s−1. Results corresponding to: (a) fc→ ∞ andff→ ∞ (i.e. f∗=f), so that the constitutive model accounts for the growth of voids while neglecting coalescence and fracture, and (b) fc= 0.12 and ff= 0.5 (see Table 1). 3.3. The effect of loading rate 427 The calculations presented in this section are performed for imposed initial major strain rates ranging from 428 100 s−1to 10000 s−1, and for three different values of the initial void volume fraction f0= 0.01, 0.02 and 0 .03. 429 430 Fig. 6 compares forming limit diagrams calculated using finite element simulations and predicted by the 431 stability analysis for four different values of the imposed initial major strain rate: 500 s−1, 2000 s−1, 5000 s−1and 432 10000 s−1. The initial void volume fraction is f0= 0.01. Results disregarding voids coalescence and fracture are 433 24 shown in Fig. 6a. The theoretical predictions are in excellent accord with the numerical results. As anticipated 434 in previous section, the stabilizing effect of inertia enhances the formability of the material as the strain rate 435 increases. Notice that the gap between the curves obtained for different strain rates is roughly constant for 436 all the loading paths considered, showing that the relative contribution of inertia to neck retardation decreases 437 as the loading path moves away from plane strain and approaches equibiaxial tension, in agreement with the 438 finite element simulations performed by N’souglo et al. (2021) for anisotropic materials modelled with Hill (1948) 439 criterion. Results considering voids coalescence and fracture are shown in Fig. 6b. We have checked that necking 440 occurs before fracture (recall that the stability analysis only models necking). As in Fig. 6a, the gap between the 441 curves obtained for different strain rates is very similar for all loading paths. Note that the difference with the 442 results reported in Fig. 6a for the simulations with f∗=fis that the formability of the material is reduced near 443 equibiaxial tension (therefore the relative influence of coalescence on the formability is more as equibiaxial tension 444 is approached). The stability analysis predicts that coalescence occurs for εc yy= 0.35, 0.38, 0.43 and 0 .51, in the 445 cases of ˙ ε0 xx= 500 s−1, 2000 s−1, 5000 s−1and 10000 s−1, respectively, while the values obtained with the finite 446 element simulations are εc yy= 0.38, 0.40, 0.44 and 0 .50. These results show that increasing the loading rate delays 447 the onset of coalescence. 448 449 The influence of the loading rate on necking formation is further investigated in Fig. 7, which pictures the 450 comparison between finite element simulations and linear stability analysis for the evolution of the critical major 451 necking strain εc xxwith imposed initial major strain rate ˙ ε0 xxfor three different values of the initial void volume 452 fraction f0= 0.01, 0.02 and 0 .03. Note that the range of values of ˙ ε0 xxconsidered spans two decades of strain rate. 453 Results are shown for two different loading paths χ= 0 and 0 .5. The εc xx−˙ε0 xxcurves display a concave-upward 454 shape, showing that the increase of the necking strain with the loading rate is nonlinear, in agreement with the 455 finite element results shown in Rodr´ ıguez-Mart´ ınez et al. (2017) and N’souglo et al. (2021) for fully dense materials 456 modelled with von Mises (1928) and Hill (1948) yield criteria, respectively. 457 Results corresponding to fc→ ∞ andff→ ∞ (i.e. f∗=f) are shown in Figs. 7a and 7b. The differences 458 between the εc xx−˙ε0 xxcurves obtained for different initial void volume fractions are very small, consistently with 459 the results previously shown in Figs. 4a and 5a. The agreement between stability analysis predictions and finite 460 element calculations is excellent for the whole range of strain rates considered in the case of χ= 0, see Fig. 7a. 461 Moreover, for the loading path χ= 0.5, the stability analysis underestimates the critical major necking strain 462 calculated with the finite element simulations for the lower strain rates investigated (underestimation of ≈6% for 463 25 00.20.40.60.811.2 0 0.2 0.4 0.6 0.8 1 1.2FEM: ε0 xx=500 s-1 FEM: ε0 xx=2000 s-1 FEM: ε0 xx=5000 s-1 FEM: ε0 xx=10000 s-1 LSA: ε0 xx=500 s-1 LSA: ε0 xx=2000 s-1 LSA: ε0 xx=5000 s-1 LSA: ε0 xx=10000 s-1 Critical minor necking strain, εyyCritical major necking strain, εxxc cε xx= ε yyc cf0=0.01 WITHOUT COALESCENCE FEM: εxx=500 s-1 FEM: εxx=2000 s-1 FEM: εxx=5000 s-1 FEM: εxx=10000 s-1 LSA: εxx=500 s-1 LSA: εxx=2000 s-1 LSA: εxx=5000 s-1 LSA: εxx=10000 s-1. . . . . . . .0 0 0 0 0 0 0 0 (a) 00.20.40.60.811.2 0 0.2 0.4 0.6 0.8 1 1.2FEM: ε0 xx=500 s-1 FEM: ε0 xx=2000 s-1 FEM: ε0 xx=5000 s-1 FEM: ε0 xx=10000 s-1 LSA: ε0 xx=500 s-1 LSA: ε0 xx=2000 s-1 LSA: ε0 xx=5000 s-1 LSA: ε0 xx=10000 s-1 Critical minor necking strain, εyyCritical major necking strain, εxxc cε xx= ε yyc cf0=0.01 WITH COALESCENCE FEM: εxx=500 s-1 FEM: εxx=2000 s-1 FEM: εxx=5000 s-1 FEM: εxx=10000 s-1 LSA: εxx=500 s-1 LSA: εxx=2000 s-1 LSA: εxx=5000 s-1 LSA: εxx=10000 s-1. . . . . . . .0 0 0 0 0 0 0 0 (b) Figure 6: Comparison between finite element results (FEM) and linear stability analysis predictions (LSA). Forming limit diagrams, critical major necking strain εc xxversus critical minor necking strain εc yy, for loading conditions ranging from plane strain stretching (χ= 0) to nearly equibiaxial tension. Results are shown for four values of the imposed initial major strain rate ˙ ε0 xx= 500 s−1, ˙ε0 xx= 2000 s−1, ˙ε0 xx= 5000 s−1and ˙ε0 xx= 10000 s−1. The initial void volume fraction is f0= 0.01. Results corresponding to: (a) fc→ ∞ andff→ ∞ (i.e.f∗=f), so that the constitutive model accounts for the growth of voids while neglecting coalescence and fracture, and (b) fc= 0.12 and ff= 0.5 (see Table 1). 26 ˙ε0 xx= 100 s−1). 464 Results corresponding to fc= 0.12 and ff= 0.5 are shown in Figs. 7c and 7d. The εc xx−˙ε0 xxcurves obtained 465 for plane strain tension χ= 0 are very similar to those obtained in Fig. 7a, showing that necking occurs before the 466 coalescence criterion is met, for this loading path, for the whole range of strain rates considered. In contrast, for the 467 loading path χ= 0.5, there are differences between the results obtained for the three initial void volume fractions 468 considered, so that increasing the value of f0shifts the εc xx−˙ε0 xxcurves downwards. In fact, the comparison with 469 the results obtained in Fig. 7b brings out that for f0= 0.02 and 0 .03 the critical major necking strain is less 470 when coalescence is taken into account, specially, for strain rates greater than ≈2000 s−1. Notice that the results 471 obtained with both finite element calculations and linear stability analysis show the same overall trends and there 472 is also quantitative agreement for most of the loading rates considered (the accord is particularly good at high 473 strain rates for the three void volume fractions considered). 474 4. Concluding remarks 475 In this paper, we have developed an analytical model to predict the formation of dynamic necks in ductile 476 plates subjected to in-plane biaxial stretching for loading paths ranging from plane strain to nearly equibiaxial 477 tension, and loading rates varying from 100 s−1to 10000 s−1. The novelty with respect to other stability analysis 478 currently available in the literature is that it takes into account the material porosity using the Gurson-Tvergaard- 479 Needleman constitutive relation for progressively cavitating solids (Gurson, 1977; Tvergaard, 1981, 1982; Tvergaard 480 and Needleman, 1984). The calibration of the analytical model has been performed assuming that a perturbation 481 mode turns into a necking mode when the cumulative instability index meets a critical value that depends on the 482 loading path and the strain rate. Excellent qualitative and quantitative agreement has been obtained between 483 formability diagrams predicted with the stability analysis and calculated with finite element simulations for initial 484 void volume fractions spanning from 0 .01 to 0 .1. Both analytical and numerical results show that the influence 485 of porosity on necking localization increases, due to early voids coalescence, as the loading rate increases and the 486 loading path approaches equibiaxial tension. The stability analysis derived in this paper, which stands out for its 487 simplicity, is intended to provide practical indications on the effect of porosity on the dynamic necking formability 488 of sheet metal products. 489 27 00.20.40.60.811.2 100 1000 10000FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 Imposed initial major strain rate, εxxCritical major necking strain, εxxc 0χ=0 WITHOUT COALESCENCE .FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 (a) 00.20.40.60.811.2 100 1000 10000FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 Imposed initial major strain rate, εxxCritical major necking strain, εxxc 0χ=0.5 WITHOUT COALESCENCE .FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 (b) (c) 00.20.40.60.811.2 100 1000 10000FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 Imposed initial major strain rate, εxxCritical major necking strain, εxxc 0χ=0.5 WITH COALESCENCE .FEM: f0=0.01 FEM: f0=0.02 FEM: f0=0.03 LSA: f0=0.01 LSA: f0=0.02 LSA: f0=0.03 (d) Figure 7: Comparison between finite element results (FEM) and linear stability analysis predictions (LSA). Evolution of the critical major necking strain εc xxwith imposed initial major strain rate ˙ ε0 xxfor three different values of the initial void volume fraction f0= 0.01, 0.02 and 0 .03. Results are shown for two different loading paths: χ= 0 for (a) fc→ ∞ andff→ ∞ (i.e.f∗=f) and (c) fc= 0.12 andff= 0.5, and χ= 0.5 for (b) fc→ ∞ andff→ ∞ (i.e.f∗=f) and (d) fc= 0.12 and ff= 0.5 (see Table 1). 28 Funding 490 The research leading to these results has received funding from the European Research Council (ERC) under 491 the European Union’s Horizon 2020 research and innovation programme. Project PURPOSE, grant agreement 492 758056. 493 Conflict of interest 494 The authors declare that they have no conflict of interest. 495 Supplementary material 496 In the following link we provide access the codes developed in Wolfram Mathematica to compute the stability 497 analysis results presented in Section 3: https://doi.org/10.21950/2XVATI. 498 Author contributions 499 M. Anil Kumar : Conceptualization; Data curation; Formal analysis; Investigation; Software; Validation; 500 Writing - original draft. K. E. N’souglo : Conceptualization; Data curation; Formal analysis; Investigation; 501 Methodology; Software; Supervision; Validation; Writing - original draft; Writing - review & editing. J. A. 502 Rodr´ ıguez-Mart´ ınez : Conceptualization; Formal analysis; Funding acquisition; Investigation; Methodology; 503 Project administration; Resources; Supervision; Validation; Visualization; Roles/Writing - original draft; Writing 504 - review & editing. 505