The present paper concerns a computational study of a three-dimensional (3D) unit cells with identical spherical voids in a von Mises matrix. The objective is to estimate the effective plastic flow surface of 3D microstructures. The work originality is to deal with identical spherical overlapping voids, covering a wide range of stress triaxiality ratios. The effective plastic flow surface is computed for nine distinct loadings on four distinct microstructures. The result indicates that the classical Gurson–Tvergaard–Needleman (GTN) model obtained using the Fritzen et al. [40] parameters, matches with our numerical simulations.

The mathematical development of yield criteria for plastic porous solids has been widely investigated, see Rice and Tracey [1], Gurson [2], Tvergaard [3], Needleman and Tvergaard [4], Becker et al. [5], Koplik and Needleman [6], Sun and Wang [7], Ponte Castaneda [8], Michel and Suquet [9], Gologanu et al. [10,11], Garajeu [12], Zuo et al. [13], Garajeu and Suquet [14], Gologanu et al. [15], Faleskog et al. [16], Ma and Kishimoto [17], Corigliano et al. [18], Pardoen and Hutchinson [19], Zhang et al. [20], Gologanu et al. [21], Negre et al. [22], Kim et al. [23], Wen et al. [24], Zaïri et al. [25], McElwain et al. [26], Monchiet et al. [27], Zaïri et al. [28], Besson [29], Laiarinandrasana et al. [30], Li and Karr [31], Nielsen and Tvergaard [32], Vadillo and Fernandez Saez [33], Zadpoor et al. [34], Lin et al. [35], Dunand and Mohr [36], Li et al. [37], Mroginski et al. [38], Fei et al. [39], Fritzen et al. [40], Madou and Leblond [41,42], Yan et al. [43], Benhizia et al. [44], Khdir et al. [45] and Khdir et al. [46].

Essentially, because of the role of porosities regarding the ductile fracture process, these voids are the consequence of manufacturing processes.

The mathematical derivations of these criteria are generally based on the continuum-based micromechanical framework, for which the starting point is the microstructural representation of the porous medium. The non triviality of the theoretical problem leads to define a basic unit cell containing one centered void for the material volume used to represent the microstructure.

Gurson [2] proposed the most widely used micromechanics-based yield criterion to analyze plastic porous solids containing spherical voids. The Gurson model is based on the following assumptions: isotropy, incompressibility of the surrounding matrix and rigid plasticity for the local yielding of the surrounding matrix material, which obeys to the von Mises criterion. The Gurson resulting macroscopic yield criterion for porous medium is hydrostatic pressure dependent. It integrates the porosities volume fraction as a model parameter and accounts for a possible void growth driven by the local plastic deformation of the surrounding matrix material. As pointed out by Tvergaard [3], the Gurson model gives an upper bound of the macroscopic yield stress as a function of the mean stress for voids periodic arrangement.

In order to get closer to two-dimensional finite element simulation results on a periodic unit cell, Tvergaard [3] proposed to introduce heuristic parameters in the Gurson yield criterion. These adjustable parameters have no direct physical meaning, but may be correlated to interaction effects between voids.

The extension of the Gurson model by Tvergaard [3], known as the Gurson–Tvergaard–Needleman model (GTN), was widely used by researchers to check its capability to highlight the poro-plastic behavior of many porous materials. Benzerga and Leblond [47] and Besson [48] reviewed the various extensions of the Gurson model based on enhanced micromechanical approaches or upon phenomenological generalizations to take into consideration the void shape or the matrix material features such as isotropy, kinematic hardening, visco-plasticity, compressibility and anisotropy.

Using micromechanical approaches, Ponte Castaneda [8] and Sun and Wang [7] proposed, respectively, upper and lower bounds for the overall yield surface of porous media. Using the variational technique, introduced by Ponte Castaneda [8], Garajeu and Suquet [14] proposed another upper bound which overcomes the well known basic drawbacks of the Gurson criterion at low stress triaxiality values.

The effect of void shape on the macroscopic yield response of porous materials was investigated by several authors. For more details, see Gologanu et al. [10,11], Garajeu [12], Yee and Mear [49], Gologanu et al. [15,21], Son and Kim [50], Siruguet and Leblond [51], Flandi and Leblond [52], Li and Huang [53], Li and Steinmann [54], Monchiet et al. [27], Gao et al. [55], Keralavarma and Benzerga [56], Lin et al. [35], Lecarme et al. [57], Scheyvaerts et al. [58], Zaïri et al. [59], Danas and Aravas [60], Madou and Leblond [41], Khdir et al. [45] and Khdir et al. [46]. Some studies concerned with the influence of the 3rd invariant of the macroscopic stress tensor have been presented by Hsu et al. [61] and Gao et al. [55,62].

The mathematical developments have reached a high degree of sophistication. The resulting yield criteria generally involve a certain number of parameters with no physical significance. This may be explained by the fact that these micromechanics based models consider, as a material volume element, an elementary volume element containing a single void. As the voids are diluted in the matrix material, the interactions between voids are neglected.

In order to be statistically representative, the material volume element should contain sufficient information about the porous microstructure, in particular the void distribution.

The material response of porous media was also investigated using computational micromechanics. This approach is emerging as a powerful tool to bring a better understanding of void distribution effects and interaction phenomena on the mechanical behavior of random porous media.

The main advantage of the computational homogenization is its ability to directly compute the mechanical fields on the random porous media by representing explicitly the microstructure features such as shape, orientation and distribution of voids.

Although many studies were dedicated to the development of yield criteria for plastic porous media, only a few works are devoted to three-dimensional (3D) computational homogenization involving multiple voids. To the best of our knowledge, only Bilger et al. [63,64], Fritzen et al. [40], Khdir et al. [45], Vincent et al. [65] and Khdir et al. [46] used this approach to estimate the overall yield surface of porous materials. Their computations were limited to spherical voids and ellipsoidal shape in Vincent et al. [65] and Khdir et al. [46]. The calculations of Bilger et al. [63,64] and Moulinec and Suquet [66,67] were performed on the basis of 3D Fast Fourier Transform. The pore clustering effect on the overall material response was the key point of their investigations.

Fritzen et al. [40] assumed the random porous media as a volume of porous material, which is periodically arranged. The results highlighted by Fritzen et al. [40] led to the extension of the GTN yield criterion in order to overcome the analytical/numerical discrepancies. In Vincent et al. [65], numerical simulations are performed on 3D cells containing numerous spherical or ellipsoidal (oblate) voids in a von Mises or GTN matrix. The cavities, in this study, cannot overlap (each void is surrounded by a matrix spherical shell to avoid overlapping) and the numerical method is based on Fast Fourier Transform. Khdir et al. [45] focused their investigations on the porous materials containing two populations of voids. Their results showed that, there is no significant difference between a double and a single population of voids for an identical fraction of porosities. A computational homogenization of random porous media, including spherical identical overlapping voids, is used in order to determine their overall yield surface. Vincent et al. [65,68] have derived a simple analytical expression of the effective flow surface obtained by generalizing the Gologanu et al. [11] approach to compressible materials. The material under consideration exhibits two populations of voids, where the smaller voids are spherical voids and the larger are spheroidal and randomly oriented inside the material.

The computational investigations performed in this study account for the complex coupling existing between void distribution, void shape (overlapping or non overlapping) and external loading mode, in order to compare the simulation results with some Gurson type yield criteria and verify the extension of the GTN yield criterion provided by Fritzen et al. [40] in the case of random porous media containing overlapping spherical voids.

2Numerical homogenization2.1Porous microstructuresThe porous media considered in the computations are made of perfectly plastic matrix, obeying to the commonly used isotropic von Mises yield criterion. The matrix yield stress is constant for all studied microstructures. All the results (the effective stresses) can be normalized with respect to this. In order to disregard hardening effects in the investigation and to compare the simulation results with the most common analytical models, the plastic flow is assumed perfect. The consideration of the path-dependent plasticity and the material elastic strain is a noteworthy difference with respect to existing works where only few studies have accounted for these effects. The matrix material is stiff enough to overcome any yield strain effects. The voids are assumed to be pressure free. The porous media are represented by 3D cubic cells containing a large number of identical, overlapped and randomly distributed, spherical pores, see Fig. 1, in order to assure that the studied material volume element is sufficiently large if compared to porosities.

Four microstructures are studied with different porosities, p=0.13=13%, 0.23, 0.4 and 0.5. Contrarily to existing literature, see Bilger et al. [63,64], Fritzen et al. [40] and Khdir et al. [45], the overlapping voids effects are studied in this work

2.2Morphology and finite element meshThe computation method of a 3D elastic plastic porous material, containing random spherical overlapping voids, is used. To generate our microstructures, we chose at first, randomly points M1, M2,..., Mi,..., MN according to Poisson process, see Jeulin and Ostoja Starzewski [69], Kanit et al. [70], Kanit et al. [71], El Moumen et al. [72], El Moumen et al. [73,74] and Kaddouri et al. [75]. Then, we construct each ith spherical void (i=1, 2,..., N) for each center Mi regardless of the repulsion distance in the construction of materials with overlapping voids. Fig. 2 shows examples of used porous material with N=200 random overlapping spherical voids.

A standard small strain approximation was used for the simulations. A geometrically linear description is employed to allow comparison with existing analytical models which are valid in the linear setting only. The infinitesimal strain tensor is defined as the symmetric part of the gradient of the displacement field. The pores are assumed to be pressure free. More specifically, the surface of the pores is assumed as free boundary and the traction vector is zero on this surface. Note that in the given geometrically linear setting the evolution of the porosity due to applied loading cannot be attributed for.

The finite element method with multi-phase element technique was chosen for the numerical computations. The finite element mesh, associated with the image of the microstructure, is obtained, using the so-called multi-phase element technique. This method was developed by Lippmann et al. [76] and extensively used for the homogenization of virtual and real images by several authors, in linear elasticity and in thermal conductivity by Kanit et al. [70], Kanit et al. [77], Kanit et al. [71], El Moumen et al. [72], El Moumen et al. [78], El Moumen et al. [73,74], Kaddouri et al. [75], Djebara et al. [79] and in the case of plasticity, these methods are extended to approximate the elastic–plastic behavior of particle-reinforced composites as shown in Khdir et al. [45,80], Benhizia et al. [44], Khdir et al. [46] or for polycrystalline aggregates behavior in Barbe et al. [81]. In the multi-phase element technique, an image of the microstructure is used to attribute the proper phase property to each integration point of a regular mesh, according to the color of the underlying voxel. It consists of superimposing a regular finite element grid, see Fig. 3a, on the image of the microstructure, see Fig. 3b. The obtained meshed microstructure is presented in Fig. 3c.

The main drawback of this simple technique is that in the same finite element two different phases can be present. The element edges do not necessarily follow the interfaces of phases (matrix and voids) in the microstructure. For a better representativity of the microstructure, and for more precision, we used a quadratic finite element of 20 nodes with a full integration of 27 integration points. First of all, a mesh size test is performed, on a Representative Volume Element (RVE) with N=200 voids with porosity p=0.13, see Fig. 2, to determine the optimum mesh grid. For that purpose, we mesh the same microstructure, with different grids by increasing their size from 125 to 42875 ﬁnite elements, see Fig. 4.

The time of calculation and the memory required in terms of ﬁnite elements number are presented in Fig. 5. To perform this test, the sixth loading from Table 2 characterized by the two loading parameters (α=1, β=1/2) is used.

The convergence of the macroscopic property as a function of ﬁnite elements number is presented in Fig. 5. The retained mesh density is the one which allows the determination of the macroscopic property with a good precision in a minimum time. It appears that the stability of the curve starts from 8000 ﬁnite elements for N=200. This mesh size estimates the macroscopic property with a precision of approximately 3% in 240h. The mesh size of 27000 ﬁnite elements estimates the macroscopic property with a precision of approximately 1% in 810h. As a result of convergence, a mesh density of 27000 ﬁnite elements was adopted in this study for all simulations. As result of convergence, the mesh density of 27000 ﬁnite elements (729000 integration points) was enough to represent accurately the geometry of the porosity and has insured the overall response convergence. This mesh density was adopted in this study for all simulations.

2.3Number of realizations and RVE sizesIn the statistical approach used in our study, the mesh of n different realizations for each volume N of porous material is hard and long to achieve with a free mesh, whereas this process can be automated with a regular mesh. The technique, developed by Kanit et al. [70], is used to determine the RVE size of different microstructures with overlapping voids. For the case of non overlapped voids the RVE size is evaluated by Khdir et al. [80]. The RVE is the volume that allows the estimation of the effective property with one realization. This technique uses the finite elements computations on different realizations with the same size. Note that the realizations in the same conﬁguration have the same number of voids N. The number of realizations n required for each number of voids N in each volume, is defined according to the statistical approach given by Kanit et al. [70], and are summarized in Table 1. The number n of different realizations considered for each volume size N is given in Table 1. This number is chosen so that the obtained mean value of the flow surface, of our studied porous material, and variance do not vary up to a given precision (less than 2% here). As a result, the overall flow surface in a porous material can be determined either by a few number of measurements on large volumes, or by many realizations for small volumes of material, see Kanit et al. [70] and Khdir et al. [80].

Fig. 6 presents an example of three realizations containing 200 voids, 8 realizations containing 50 voids and 27 realizations containing 5 voids.

The convergence of apparent properties to the effective properties allowing the determination of the RVE size is deducted from the effective yield surface study for the two extreme volume fractions p=0.13 and p=0.5. Different realizations are used to evaluate the mean value of the apparent property for each volume. In this test, the sixth loading characterized by the two loading parameters (α=1, β=1/2) are also used. It should be mentioned here that the apparent property is the results given by volumes smaller than the RVE size. Then, the mean values and the confidence intervals are plotted in terms of voids number N. Fig. 7 shows that the results dispersion decreases when the volume size increases.

The error bars decreases when the volume size increases and tends to zero for the RVE size. The error bars, for each voids number N, are determined with n realizations as defined in Table 1. In the case of p=0.13, the RVE is in the range of 50 voids but for p=0.5 the RVE has a size in the order of 100 voids. For better precision, we adopt for all simulations a RVE of 200 voids.

2.4Boundary conditions and loadingSince the porous media is hydrostatic pressure dependent, the boundary conditions imposed to the RVE should involve a wide range of stress triaxiality ratios. The stress triaxiality parameter T=Σm/Σeq is defined as the ratio of the overall hydrostatic stress Σm and the overall von Mises equivalent stress Σeq, given by Σm=trace(Σ)/3 and Σeq=(2/3)(Σ′:Σ′). Σ is the applied macroscopic stress tensor and Σ′ denotes its deviatory part.

Due to the computational robustness, the following mixed boundary conditions are imposed: E11(t)=tε˙0(β+α), E22(t)=tε˙0(β−α), E33(t)=tε˙0β and E12=E13=E23=0. The shear components assigned values of the overall stress tensor are zero. The terms α and β, introduced to control the diagonal components of the overall strain tensor E, are two loading parameters, ε˙0>0 is a prescribed deformation rate, and t is the simulation time. These conditions are numerically realized by imposing only the normal mean displacement on the unit cell surface. Under the given loadings, the mean shear tractions on the surfaces vanish and the shear components of the macroscopic stress tensor are always zero. Thus, the von Mises and the hydrostatic stresses depend only on the macroscopic stress diagonal components. The stress triaxiality is indirectly assigned by the stress two measures, Σm and Σeq, which are implicitly defined by the mixed boundary conditions through the two loading parameters α and β. Due to the chosen mixed loading, the two stress measures contain only three independent parameters. The different values of α and β used to obtain different stress triaxiality ratios are listed in Table 2.

3Brief survey of used analytical models3.1Gurson criterionGurson [2] considered a matrix perfectly plastic with a spherical cavity, obeying to plasticity criterion of von Mises yield strength σ0, subjected to conditions of uniform and homogeneous strain rate applied to the outer edge. The macroscopic criterion, represents the plastic potential with the function of the flow surface depending on the macroscopic stress and the voids volume fraction, verify: a(Σeq/σ0)2+2pcosh(3Σm/2σ0)−1−p2=0.

3.2Gurson–Tvergaard criterionThe modifications of the original Gurson material model performed by Tvergaard [3], for porous media involving the introducing of parameters qi to describe plastic properties of the material taking into account the interaction between cavities, are given by the following threshold function: (Σeq/σ0)2+2q1pcosh(3q2Σm/2σ0)−1−q3p2=0.

Several values of the parameters qi, to approximate the structures real behavior, were determined and proposed by many authors, see Table 3.

Values of GTN model parameters, q3=q12.

References | q1 | q2 |
---|---|---|

Gurson [2] | 1 | 1 |

Tvergaard [3] | 1.5 | 1 |

Koplik and Needleman [6] | 1.25 | 1 |

Zuo et al. [13] | 1.4 | 1 |

Faleskog et al. [16] | 1.46 | 0.93 |

Ma and Kishimoto [17] | 1.35 | 0.95 |

Corigliano et al. [18] | 1.08 | 0.99 |

Zhang et al. [20] | 1.25 | 1 |

Negre et al. [22] | 1.5 | 1.2 |

Kim et al. [23] | 1.58 | 0.91 |

McElwain et al. [26] | 1.31 | 1.16 |

Nielsen and Tvergaard [32] | 2 | 1 |

Vadillo and Fernandez Saez [33] | 1.46 | 0.93 |

Dunand and Mohr [36] | 1 | 0.7 |

Fei et al. [39] | 1.8 | 1 |

Yan et al. [43] | 1.55 | 0.9 |

In this section, we present and discuss the numerical study we conducted. We first start with the asymptotic stress response, then we pass to a subsection on the representativity, and end with a subsection on the local plastic strain fields.

4.1Asymptotic stress responseThe asymptotic response of the ideally plastic porous microstructures was systematically examined by plotting the overall normalized von Mises equivalent and hydrostatic stresses as a function of the overall von Mises equivalent strain. For each considered microstructure presented in Fig. 2, these two stress measures are plotted for each volume fraction: p=0.13, 0.23, 0.4 and 0.5, see Fig. 8. The arrows illustrate the evolution of the two stresses measures for the nine different loading cases given in Table 2, showing that the macroscopic von Mises equivalent stresses and the hydrostatic stresses vary inversely for the nine different loading cases.

Macroscopic normalized von Mises equivalent stress Σeq/σ0 and normalized macroscopic hydrostatic stress Σm/σ0 as a function of the macroscopic von Mises equivalent strain Eeq. The curves correspond to the different loading cases given in Table 2.

Fig. 8 shows that the porous microstructures are subjected to a stationary response beyond a certain strain. Similar observations on large volume computations have been pointed out by Fritzen et al. [40], Khdir et al. [45,46]. This result is obvious (there is an asymptotic limit), because we suppose that the geometry is fixed. In reality, this one can evolve during the calculations, so the cavities can grow during the loading, and, in this case, there is no more stationary response. In order to define the numerical yield points, the overall stresses at the end of the simulation are considered.

4.2RepresentativityThe size of the RVE is conditioned by the number of voids which should be chosen large enough to ensure that the volume element is representative. This representativity was investigated in terms of the mechanical responses by Huet [82], Drugan and Willis [83], Kanit et al. [70,77], Jiang et al. [84], Kanit et al. [71], El Moumen et al. [73,74,78], Djebara et al. [79] and Kaddouri et al. [75]. These authors have studied the effects of the volume element size on the elastic stiffness. Khdir et al. [80] have investigated these effects on the elastic plastic response. Vincent et al. [68] proposed a new measure for the deviation from isotropy with respect to the effective plastic flow surface. In the latter case, made of two phases with highly contrasted properties, Khdir et al. [80] have shown that the minimum size of the volume element in the yield and post yield region must be greater than the minimum size required in the elastic domain. This question, which arises in 3D computational homogenization, has to be systematically accounted for several volume elements with different sizes (i.e. containing different number of pores) are simulated for several porosities.

4.3Local plastic strain fieldsExamples of the distribution of local plastic strain fields for the studied porosities: p=0.13, 0.23, 0.4 and 0.5 are illustrated in Fig. 9.

The cases (α, β)=(1, 0) and (α, β)=(0, 1) corresponding to Σm=0 and Σeq=0, are the lowest and highest triaxiality ratios respectively and the case (α, β)=(1, 1/4) corresponds to the intermediate case. The observations are presented at the end of the loadings. The main observation is the local fields differences between different porosities. We can see that plastic zones are more important for low volume fractions, decreasing progressively as the volume fraction increases. The maximum accumulated plastic strain values are observed regardless of the volume fraction. Therefore, we have the same amount of accumulated plastic strain which concentrates in the matrix independently of the volume fractions. The stress triaxiality which refers to the ratio of equivalent stress and mean stress is not sufficient to represent the actual material behavior. For isotropic viscoplastic polycrystals the 3rd invariant can influence the effective constitutive response, see Böhlke [85]. For porous media, Kim et al. [23], Danas et al. [86] and Cazacu et al. [87] investigated the subject. Usually the 3rd invariant is used in the context of the Lode angle. In our work, the influence of the 3rd invariant should be eliminated from the presented results. We evaluated the Lode angle for different realizations undergoing the same loading. The small fluctuations found confirm the choice of mixed boundary conditions.

4.4Comparison between numerical results and analytical criteriaThe common representation of the overall yield surface, plotting the overall von Mises equivalent stress as a function of the overall hydrostatic stress, is adopted to illustrate the computational data. In Fig. 10, an example of the relationship between these two stress measures is represented, for real porosity values: p=0.13, 0.23, 0.4 and 0.5.

The normalization is performed with respect to the matrix yield stress σ0. The hydrostatic pressure dependency of the porous material macroscopic yield response is clearly highlighted in Fig. 10. Although an asymptotic response is reached by the perfectly plastic cubic cell, see Fig. 8. The relationship between the two measures of stress does not follow a straight line but a rather complex non proportional path, also observed by Fritzen et al. [40]. The yield points are obtained at the end of the simulation. The computed yield data strongly point out the convexity of the overall yield surface.

4.5GTN model for random porous mediaTo improve its agreement with the computational results, the GTN model can be compared and calibrated using our computed data. The GTN model parameters are as follows, q1=1.6−p, q2=0.9 and q3=q12. The model represents all computed data in a very satisfactory manner as shown in Fig. 11.

We get the same expressions as those found by Fritzen et al. [40]. This calibration can be compared with the values reported in Table 3, usually obtained by calibration on two-dimensional finite element simulation results using plane stress, plane strain or axis-symmetric periodic unit cell models.

Fig. 11 shows that there is no significant effect of the void shape, overlapping or non-overlapping spherical voids, on the overall yield response. This result does not match the micromechanics-based analytical yield criteria developed in the literature; see e.g. Gologanu et al. [10,11,15,21] and Madou and Leblond [41]. This difference is due to the fact that the cited author's results deal with ellipsoidal non overlapping voids contrarily to our case where we considered overlapping spherical voids.

The theoretical developments consider unit cell representations, periodicity is assumed and void shape dependence of the overall yield response is thus expected. From a micromechanics-based analytical point of view, the overall response of a random porous medium is evaluated from that of a unidirectional unit cell averaged over all orientations; see Zaïri et al. [88] and Vincent et al. [89]. Shape dependence is thus preserved in the micromechanics-based analytical models.

The exact number of calculations is: 216+27×5+2×20+8×50+8×100+3×200 (from Table 1)+4× (on RVE)=80876. If the large volume of computations (80876 calculations) performed in this study show no significant effect of the void shape on the volume average behavior, this could be a consequence of the cubic cell microstructure in which the pores are simultaneously randomly distributed and oriented in space. However, this statement must be verified for a larger range of shape ratios.

5Concluding remarksThe macroscopic mechanical behavior of porous materials with perfectly plastic von Mises type plasticity was studied on a computational basis. The overall yield surface of plastic porous media, with overlapping identical spherical voids, was investigated via computational micromechanics. The random microstructures consisting of spherical overlapping voids are used in finite element simulations and the deviatoric and hydrostatic limit stresses are captured for all computations.

For the class of porous materials under consideration, it is found that the analytical predictions obtained in the exhaustive literature on the topic are all rather close to the numerical prediction for small pore volume fractions. Notably, the analytical models are rather close to the numerical results for pore fractions up to 2.5%. For higher pore fractions, significant deviations between the various models are found. Therefore, our numerical results are compared with a modification of the computationally attractive GTN model given in Fritzen et al. [40]. This extension is appealing due to a single additional parameter and its ability to fit all calculated points to an excellent extent.

Based on the statistical results and results obtained for different types of boundary conditions and numbers of pores, the computational results can be considered representative under the stated constitutive assumptions. The computational results were investigated in terms of representativity and were related to the Gurson-type yield criteria.

Our numerical results for overlapping porous materials are close to results obtained by the Gurson-type yield criteria given in Fritzen et al. [40]. The overall yield surfaces were found nearly the same for overlapping or non-overlapping spherical voids provided that they are randomly distributed in a large volume element. So we deduce that there is no difference between overlapping or non-overlapping porous material if the porosity is the same.

Future research can be dedicated to the consideration of more general loading conditions and of the role played by the 3rd invariant of the stress tensor. Such studies could consider explicit variation of the Lode angle which is almost constant in our study.

Conflicts of interestThe authors declare no conflicts of interest.