Research paper

# Analysis of the mechanical behavior of porous materials containing two populations of voids under dynamic spherical loading

Mohammed El Ansī, Cédric Sartori*, Christophe Czarnota

Université de Lorraine, CNRS, Arts et Métiers Sciences et Technologie, LEM3, F-57000, Metz, France

## ARTICLE INFO

Handling editor: J Molinari

Keywords:

Micro-inertia effects

Porous microstructures

Numerical homogenization

Void interaction

## ABSTRACT

A computational homogenization analysis is performed on three-dimensional representative volume elements (RVE) that contain two distinct populations of voids, to demonstrate the influence of the interaction between cavities. RVEs are constructed as cubic elastic perfectly-plastic matrices embedding two families of spherical voids, and subjected to dynamic loading under homogeneous kinematic boundary conditions. Multiple micro-structure models are considered by varying the number, position, and size of voids to evaluate the micro-inertia contribution to the overall macroscopic stress. Velocity fields within numerical RVEs are investigated to reveal the interaction of voids and their key role in the dynamic macroscopic response of the porous material. Based on numerical simulation results, a homogenization analytical model considering void interaction is proposed to describe the mechanical behavior under dynamic loadings. This model relies on adjusting the strain rate level at the boundary of unit cells composing the porous material. Comparison between our numerical and analytical results with those obtained using the classical Taylor homogenization scheme highlights the limitations of the Taylor model in the case of porous materials under dynamic loading.

## 1. Introduction

The design of materials with desired mechanical properties requires a profound understanding of their mechanical behavior and the ability to predict fracture phenomena. This is particularly pertinent in the case of ductile materials, where fracture occurs through three distinct stages: void nucleation, void growth, and void coalescence, see Benzerga and Leblond (2010). In this context, ductile materials undergoing damage are often modeled as porous materials. Rice and Tracey (1969) were among the first to consider the response of a spherical void embedded in an infinite medium. The matrix is assumed to be rigid perfectly plastic with a uniform loading in the remote boundary. Using a variational approach in a static analysis, the authors succeeded to identify the void radius expansion. In 1977, Gurson (1977) has developed one of the most classical models of ductile damage. His work, based on a micro-mechanical approach, provides an approximate yield criterion for ductile porous material. The analysis was performed on a hollow sphere made of rigid perfectly plastic material obeying the von Mises criterion, and subjected to uniform boundary strain rate. This model was improved by incorporating additional coefficients and parameters to match with numerical analysis (Tvergaard (1982); Needleman and

Tvergaard (1984); Koplik and Needleman (1988)). Ponte Castañeda (1991) and Michel and Suquet (1992) developed a bounding model to obtain accurate bounds of the effective potential for porous materials. The shape of voids is an additional parameter which influences the porous materials behavior. Following the approach of Gurson (1977), Gologanu et al. (1993, 1994) developed models for spheroidal cavities embedded in a finite volume. Morin et al. (2016) explored the effect of shape and voids orientation on the response of porous materials using numerical simulations. Reboul and Vadillo (2018) modified the Gurson's (1977) model to include strain rate sensitivity to analyze the expansion of a long cylindrical void.

Most of works cited above have used a unit cell containing a single void to represent a porous medium. Several studies have explored interactions between voids considering RVE containing multiple voids. Vincent et al. (2008, 2009) developed an approximation of the effective flow surface, similar to the Gurson-Tvergaard-Needleman criterion, of a porous material containing two populations of voids with different sizes and shapes. Julien et al. (2011) developed a semi-analytical model for the behavior of saturated viscoplastic materials containing two populations of spherical cavities of different sizes.

In addition to analytical approaches, many studies use three-

<footnote>* Corresponding author.
E-mail address: cedric.sartori@univ-lorraine.fr (C. Sartori).</footnote>
dimensional numerical VERs to propose micromechanical models. For instance, Garajeu et al. (2000) developed a micro-mechanical approach to ductile damage using numerical simulations. Their study focused on two types of unit cells: cylindrical and spherical. The primary objective of these simulations was to validate the consistency of the analytical model, which was designed to simulate ductile damage through void growth in a nonlinear viscous matrix. In the work of Vishnu et al. (2023), the authors use cubic numerical RVEs containing multiple microstructural configurations (number and size of voids, porosities, ...) to evaluate the mechanical characteristics of these RVEs compared to simple unit cells and to study the growth of voids. It is interesting to note that among the many conclusions of this work, it is highlighted that the interactions between pores directly affect the growth of voids. Kim et al. (2003) conducted finite element calculations on porous representative volume elements to explore the relationship between void shape and material anisotropy. They focused on determining the effective yield criterion for a porous, anisotropic material. The methodology developed by the authors is based on the Hill-Mandel homogenization theory, considering that the RVE is subjected to uniform strain rate boundary conditions.

Using three-dimensional Fast Fourier Transformation (FFT) calculations, Bilger et al. (2007) studied the effect of non-uniform distribution of voids on the macroscopic response of porous media with rigid-perfectly plastic matrix. Fritzen et al. (2012) built their studies on finite element simulations to investigate the effective response of ductile metals containing spherical pores with volume fraction less than 30%. To investigate the effects of different populations of voids, Khdir et al. (2014) conducted numerous computations on porous media containing two populations of spherical voids. They concluded that the response of a material with two populations of voids can be represented by an equivalent medium, of same overall porosity, containing a random distribution of voids of the same size. From a similar viewpoint, Espeseth et al. (2021) conducted numerical simulations on a unit cell containing a single large spherical void (primary void). The matrix is modeled using a porous plasticity model that considers the effects of void (secondary voids) size. The simulations showed a close relationship between the secondary void size and an intrinsic length scale (a measure of the material's characteristic size). More recently, Kowalczyk-Gajewska et al. (2024) explored the elastic properties of porous materials produced by additive manufacturing, focusing on three types of pore distributions with cubic symmetry. Their experimental results were compared with mean-field estimates and numerical calculations. This approach demonstrated excellent consistency between analytical estimates and numerical predictions across a wide range of porosities and for all examined microstructures, which included various pore distributions. The study highlighted the substantial impact of morphological characteristics on the anisotropy degree in the effective elastic stiffness of porous materials. This was notable even when the solid phase was isotropic and the pores were spherical.

As illustrated above, extensive studies have been carried out for porous materials under quasi-static loading conditions. In case of high strain rate phenomena, e.g. during planar impact tests, materials particles located in the vicinity of voids sustain large accelerations which develop at the local scale. This micro-inertia may have a significant effect that could not be neglected when porous materials are facing high strain rate loadings. Carroll and Holt (1972) carried out one of the pioneering works in this context. They studied the behavior of a hollow sphere under dynamic compression. The material was perfectly plastic, micro inertia effects were found scaled by the square of the void radius. Johnson (1981) revisited the model of Carroll and Holt (1972) and studied spallation problems and void growth in copper materials. Cortes (1992) enriched the Carroll and Holt (1972) work by including the influence of material viscosity, strain hardening and thermal softening. These works studied the dynamic response of porous materials under spherical loading conditions.

Intending to enhance micro-inertia analytical models for application

across various loading conditions, Wang (1997), Wang and Jiang (1997) developed an approximate dynamic yield criterion based on Gurson's (1977) model. From a different point of view, Molinari and Mercier (2001) proposed an explicit relationship for the dynamic macroscopic stress tensor for porous materials. This approach has led to various micro-inertia dependent models used in several problems, like spall fracture (Czarnota et al. (2008), Jacques et al. (2010), Versino et al. (2018)), dynamic cracking (Jacques et al. (2012, 2015)) or shock propagation in cellular and porous solids (Czarnota et al. (2017); Barthélémy et al. (2016); Czarnota et al. (2020); Massarwa et al. (2024)).

Following the dynamic homogenization approach proposed in Molinari and Mercier (2001), Sartori et al. (2015, 2016) considered the effect of shape of voids under dynamic loading and developed a micro-mechanical model for porous materials with spheroidal voids. More recently, Subramani et al. (2020) investigated the response of material containing cylindrical parallel voids, and showed that the dynamic macroscopic response is influenced by both the void length and the radius of the circular cross section.

To our knowledge, only few authors were interested to three-dimensional computations of voided media under dynamic loadings. Marvi-Mashhadi et al. (2021) carried out a finite element analysis on real porous microstructures and elucidated the effect of porosity on the creation of necks and fragments in ductile rings under dynamic expansion. Micro inertia effects were accounted for through the explicit description of the porous microstructure. Vishnu et al. (2022) conducted 3D finite element simulations to investigate the impact of porous microstructures, from additively manufactured materials, on dynamic shear localisation. Results showed that voids serve as initial sites for shear band nucleation, influencing the speed and direction of their development. An increase in the number of voids, given a constant void volume fraction, leads to more shear bands, while larger pores accelerate the inception and development of these bands. Recently, in a study of dynamic formability of ductile porous sheet specimens, Nieto-Fuentes et al. (2022) showed that a homogenized porosity numerical model produces the same overall results and general trends to those obtained from a finite element model incorporating spatial and size distributions of explicitly resolved voids. In the configurations tested by Nieto-Fuentes et al. (2022), micro-inertia effects have a slight impact on the dynamic formability prior to necking, possibly elucidating the similarity in results obtained from both models.

The present study aims to explore the impact of void distribution in porous materials under dynamic loading conditions, utilizing 3-D finite element simulations. Section 2 provides an in-depth description of the porous material model, focusing on the mesh generation process and the implementation of finite element analysis. Sections 3 and 4 present the results obtained from the finite element simulations, examining both the macroscopic response and the velocity fields within RVEs. In Section 5, an analytical method is introduced to predict the macroscopic behavior of porous materials with two distinct void populations, and these predictions are compared with the results from the finite element analysis. Finally, Section 6 summarizes the key findings and conclusions of this research.

## 2. Main framework of the problem

### 2.1. Modeling the porous material with spherical voids

Several works represent a porous material as a hollow sphere (spherical matrix containing a concentric spherical void), such as the standard damage model of Gurson (1977). Under static loading, the macroscopic response depends only on the porosity of the unit cell. By contrast, under high strain rate loading, micro-inertia effects induce a size effect and the dynamic behavior of the hollow sphere depends on both, the porosity and the overall size of the RVE, or alternatively, on the void size (see Carroll and Holt (1972); Molinari and Mercier (2001)).
Therefore, the use of a hollow sphere as a porous material RVE is appropriate if voids share identical sizes and when they are isotropically distributed within the matrix material, which is rarely observed in real materials. Under static conditions, the hollow sphere model also allows the description of materials composed of different sizes of spherical cells, where each cell has the same void volume fraction. This is not the case under dynamic loading where the void size enters into play, see Jacques et al. (2012) for a definition of an equivalent void size for this type of voided solids.

Here a Representative Volume Element describing a porous material containing spherical voids of multiple sizes is considered, see Fig. 1a. At the macroscopic scale, we denote $\underline{\Sigma}$ the macroscopic stress tensor and $\hat{D}$ the macroscopic strain rate tensor. At the scale of a given unit cell, defined by its porosity $f$ and internal void radius $a$, we denote by $\underline{\Sigma}$ the stress tensor at the outer boundary of the hollow sphere, and $\hat{D}$ the strain rate tensor.

The behavior of each unit cell is described by using the dynamic homogenization scheme summarized in the next section.

## 2.2. Unit cell dynamic behavior

Considering a given spherical unit cell of inner radius $a$ and porosity $f$, see Fig. 1b. Denoting by $b$ the external boundary of the unit cell, the porosity $f$ is expressed as:

$$f = \frac{a^3}{b^3} \qquad (1)$$

Homogeneous boundary conditions of the following type are assumed:

$$\nu = \underline{D} \cdot x \quad \text{on } \partial V \qquad (2)$$

where $\partial V$ represents the boundary of the unit cell and $\underline{D}$ is supposed uniform on $\partial V$. Elasticity is neglected and the matrix material is assumed incompressible. Based on the work of Molinari and Mercier (2001), the stress at the level of the unit cell, $\underline{\Sigma}$, is the sum of a static contribution (micro-inertia independent term), $\underline{\Sigma}^{\text{sta}}$, and a dynamic contribution (micro-inertia dependent term), $\underline{\Sigma}^{\text{dyn}}$:

$$\underline{\Sigma} = \underline{\Sigma}^{\text{sta}} + \underline{\Sigma}^{\text{dyn}} \qquad (3)$$

In Eq. (3), $\underline{\Sigma}^{\text{sta}}$ is related to the constitutive response of the matrix material and can be evaluated from any potential or flow surface developed for porous material containing spherical voids (Gurson (1977);

Tvergaard (1982); Needleman and Tvergaard (1984); Koplik and Needleman (1988); Ponte Castagneda (1991); Michel and Suquet (1992); Gologanu et al. (1993, 1994); Gârâjeu and Suquet (1997); Rousselier (2001); Benzerga and Leblond (2010)). It may thus include viscous rate dependence (Flandi and Leblond (2005)), strain hardening, or incorporate plastic anisotropy (Cazacu and Revil-Baudard (2021), Hosseini et al. (2022), Ferreira et al. (2023)). The second term of the right-hand side of Eq. (3), $\underline{\Sigma}^{\text{dyn}}$, is dependent on the porosity $f$ and found to be scaled by the square of the void radius, and the matrix mass density, $\rho$, see Molinari and Mercier (2001).

In the case of spherical loading, which is the focus of the present paper, the dynamic stress is expressed as follows:

$$\Sigma_m^{\text{dyn}} = \rho a^2 \left[ \dot{D}_m \left( f^{-1} - f^{-2/3} \right) + D_m^2 \left( 3f^{-1} - \frac{5}{2}f^{-2/3} - \frac{1}{2}f^{-2} \right) \right] \qquad (4)$$

where $\underline{\Sigma}_m = \frac{\text{tr}(\underline{D})}{3}$ is the spherical part of the second-order tensor $\underline{\Sigma}$, and $\dot{D}_m$ is the time derivative of the spherical part of the strain rate tensor. The theory developed by Carroll and Holt (1972) for the dynamic compaction of powders is then retrieved. Except for the mass matrix density $\rho$, note that no matrix material parameters are involved in the dynamic contribution.

By combining Eqs. (1) and (4) one gets:

$$\Sigma_m^{\text{dyn}} = \rho b^2 \left[ \dot{D}_m \left( f^{-1/3} - 1 \right) + D_m^2 \left( 3f^{-1/3} - \frac{5}{2} - \frac{1}{2}f^{-4/3} \right) \right] \qquad (5)$$

which highlights, as an alternative, that micro-inertia effects are scaled by the square of the cell size.

For a given loading state defined by the set $(\dot{D}_m, D_m)$, Eq. (4) demonstrates that the difference of micro-inertia effects experienced by two homothetic cells depends solely on the variation in the size of the embedded void. Furthermore, under the same loading condition, Eq. (5) highlights that micro-inertia effects in two spherical cells with identical external dimensions but varying void radii (i.e., distinct porosities) differ by an extent regulated by porosity.

In the following, we restrict our attention to the case of time-independent strain rate loading, i.e. $\dot{D}_m = 0$. Eq. (4) thus reduces to:

$$\Sigma_m^{\text{dyn}} = \rho a^2 D_m^2 \left( 3f^{-1} - \frac{5}{2}f^{-2/3} - \frac{1}{2}f^{-2} \right) \qquad (6)$$

and Eq. (5) can be written as:

$$\Sigma_m^{\text{dyn}} = \rho b^2 D_m^2 \left( 3f^{-1/3} - \frac{5}{2} - \frac{1}{2}f^{-4/3} \right) \qquad (7)$$

<figure><img src="image_1.png" /><figcaption>Fig. 1. Representative Volume Element of a porous material containing spherical voids, subjected to homogeneous strain-rate condition at the macroscopic scale (a), at the level of the unit cell (b).</figcaption></figure>
Note that the dynamic stress contribution in Eqs. (6) and (7) exhibit negative values within the range [0,0.479] which encompasses the values of porosity considered in this work.

## 3. Finite element modeling

### 3.1. Description of the unit cells

Dynamic micromechanical computations have been carried out with the finite element code ABAQUS/Explicit (2022). The matrix material is elastic perfectly-plastic, obeying the J2 flow theory. The Young's modulus value is taken significantly high, $E = 10^{11}$ GPa, and the Poisson's ratio is taken as $\nu = 0.35$. So, a rigid perfectly-plastic response is almost ensured. The tensile yield stress and mass matrix density are representative of steel with $\sigma_0 = 500$ MPa and $\rho = 7800$ kg/m$^3$.

In the present work, several simulations have been performed using 3D models representing two cubic unit cells, each containing a single spherical void at the center. These unit cells are denoted by UC1 and UC2 in the following. Such family of unit cells may inform on the behavior of a porous material embedding an isotropic distribution of spherical voids, all of the same size.

UC1 and UC2 both share the same side length denoted by $l$, with $l = 790\mu m$, but have distinct void radius: $a_1 = 49\mu m$ and $a_2 = 70\mu m$. The corresponding porosities are obtained from the following relationships:

$$f_1 = \frac{4\pi a_1^3}{3l^3}, f_2 = \frac{4\pi a_2^3}{3l^3} \qquad (8)$$

which implies $f_1 = 0.001$ and $f_2 = 0.002915$ for UC1 and UC2 respectively.

Fig. 2 presents the mesh of the 3D unit cells. The unit cell UC1 (resp. UC2) has been discretized using 269048 (resp. 329696) three-dimensional elements with 8-nodes and reduced integration with hourglass control (C3D8R in ABAQUS notation). In order to achieve convergence of the results while limiting computation costs, each cell was meshed using a spherical partition positioned at a distance from the center of the void, approximately twice the pore radius. The inner region contains elements having an aspect ratio close to one and a size of 2 µm. The outer region of the domain was meshed with a bias in node spacing extending from the partition boundary to the external boundary.

The unit cells UC1 and UC2 serve a dual purpose in this study. Firstly, they are utilized to conduct dynamic calculations on individual cells, enabling a detailed examination and comparison of the discrepancies between finite element analysis results and theoretical predictions. Secondly, unit cells UC1 and UC2 play a key role in building the numerical RVEs. These RVEs are specifically designed to investigate the interactions of voids under dynamic loading conditions. This aspect of their application is further elucidated in the subsequent section.

### 3.2. Numerical RVE containing two populations of spherical voids

The dynamic behavior of porous material containing two populations of voids is analyzed using RVEs defined from a stack of eight unit cells, all out of UC1 or UC2 defined in section 3.1. All RVEs have an identical volume with $L$, the side length, verifying:

$$L = 2l \qquad (9)$$

The overall porosity $\hat{f}$ of a given RVE is obtained from:

$$\hat{f} = \theta_1 f_1 + \theta_2 f_2 \qquad (10)$$

$\theta_1$ and $\theta_2$ being the volume fractions of unit cells UC1 and UC2 inside the RVE:

$$\theta_1 = \frac{V_1}{V} \quad \text{and} \quad \theta_2 = \frac{V_2}{V} \qquad (11)$$

where $V_1$ (resp. $V_2$) is the volume occupied by all unit cells UC1 (resp. UC2) inside the RVE, and $V = V_1 + V_2$ is the total volume of the RVE. Note that unit cells UC1 and UC2 have identical volume. The RVE formed by combining $N_1$ cells of type UC1 and $N_2$ cells of type UC2 is denoted as $B_{N_1} S_{N_2}$, with the total number of cells verifying $N_1 + N_2 = 8$.

As an example, Fig. 3 illustrates the Representative Volume Element (RVE) denoted as $B_4 S_4$, created by stacking four UC1 and four UC2.

Therefore, a total of five RVE families were created by varying number of unit cells UC1 and UC2, as outlined in Fig. 4. In fact, multiple RVE configurations could be generated. However, given the spherical nature of the applied load and the arrangement of the eight unit cells within the RVE, we have restricted our analysis to RVE configurations that fulfill the following condition: A RVE is considered well-balanced when the center of mass coincides with the central point (0,0,0) and

<figure><img src="image_2.png" /><figcaption>Fig. 2. Finite element model. Geometry and mesh of one eight of each cell is shown for (a) UC1 with internal void radius $a_1 = 49\mu m$ ($f_1 = 0.001$) and for (b) UC2 with internal void radius $a_2 = 70\mu m$ ($f_2 = 0.002915$).</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Fig. 3. Example of the RVE construction, labeled as B₄S₄, formed by stacking four cells of type UC1 and four cells of type UC2. The center of the RVE coincides with the reference frame (O; x₁; x₂; x₃). The RVE is subjected to homogeneous strain rate with the following boundary conditions: v₁(±l, x₂, x₃) = v₂(x₁, ±l, x₃) = v₃(x₁, x₂, ±l) = ±lD̃m, imposed at the external faces.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 4. RVE made of eight unit cells taken out of UC1 and UC2 with corresponding overall porosity, f̂, and volume fraction of UC1 inside the RVE, ϑ₁.</figcaption></figure>

when the assembly displays cubic symmetry. Therefore, we have limited our examination to RVEs where unit cells of the same family experience the same type of loading in all three principal directions. Consequently, voids belonging to the same family are surrounded by an identical distribution of voids relative to the three reference directions. This approach ensures a coherent and systematic investigation of the material behavior under specified loading conditions.

### 3.3. Boundary conditions, solving and post processing

In this work, we restrict our attention to the case of high strain rate

spherical expansion. This state of loading is representative of the large stress triaxiality state observed at the spall plane, which develops as a result of the intersection of stress-relief waves caused by impacts or explosions (Antoun et al., 2003). In the following, the macroscopic strain rate tensor is written as:

$$ \hat{D} = \hat{D}_m \underline{I}, \quad \hat{D}_m > 0 \qquad (12) $$

where ■m = tr(■) is the spherical part of the second-order tensor ■, and I stands for the second order identity tensor.

The macroscopic response of all cells (including UC1 and UC2) was
analyzed considering a homogeneous strain rate field applied at the external boundaries. A cartesian coordinate system is established within a reference frame (O; x₁; x₂; x₃) corresponding to the principal directions of the macroscopic loadings. The following boundary condition holds:

$$v_1(\pm h, x_2, x_3) = v_2(x_1, \pm h, x_3) = v_3(x_1, x_2, \pm h) = \pm h \hat{D}_m \quad (13)$$

with h = l in the case of RVEs containing two populations of voids, see Fig. 3, and h = l/2 for unit cells UC1 and UC2. The boundary condition (Eq. (13)) is assigned to a unique master node, from which resulting force components are calculated.

Following the approach of Sartori et al. (2015), we adopted a two-stage simulation procedure where two distinct numerical simulations are conducted using Abaqus/Explicit. First, a numerical simulation under quasi-static conditions, in which geometric change of the body configuration is ignored (small-strain theory), is performed, followed by a non-linear dynamic computation. The aim of the first calculation is twofold: (i) to determine the static part of the macroscopic stress, neglecting the effects of $\hat{D}_m$ on the convective part of particle acceleration and (ii) to provide initial conditions for the non-linear dynamic computation conducted under a given constant $\hat{D}_m$. This second point prevents strong disturbances that can occur when a large velocity is applied instantaneously to a cell initially at rest. Fig. 5 illustrates a schematic of the history of the velocity applied at external faces, following the boundary condition Eq. (13), within the two-stage simulation: up to time $t^+$ during the quasi-static loading, over a time period $t^* - t^+$ for the dynamic loading (see TAB. A1 and A2 in Appendix A for the values of $t^* - t^+$ used in our calculations).

We utilized an explicit solver in both calculation stages to facilitate the transfer of results from the first simulation to the second one. Note that we conducted verification tests and confirmed that employing an implicit solver (Abaqus/Standard) for the static phase yields identical results (results not shown here). As a matter of fact, the choice of the value for the time $t^+$ has no impact on the results, and an arbitrary value was adopted.

Please note that the velocity field inside the RVE, obtained at the end of the first computation stage (static loading) shows negligible evolution over the very short time period $t^* - t^+$ of the second stage (dynamic loading). Additionally, it's worth noting that, except for regions close to the cell corners, the matrix material in well-balanced RVEs, at the start of the dynamic calculation (or equivalently, at the end of the static calculation), is nearly entirely in a plastic state.

The macroscopic stress components are calculated at the end of each computation stage as follows:

$$\hat{\Sigma}_i = \frac{\hat{F}_i}{A_0} \quad (14)$$

where $A_0 = (2h)^2$ denotes the initial area of the cubic unit cell surface, and $\hat{F}_i$ the associated resultant force component. Note that the time period $t^* - t^+$ (see TAB. A1 and A2 in Appendix A), over which the analysis is conducted in the dynamic case, is very small with the porosity increasing by 2% of the initial value, and as such, the initial area is used in Eq. (14). The static stress components are obtained from the first numerical simulation as:

$$\hat{\Sigma}_i^{\text{sta}} = \hat{\Sigma}_i(t^+) \quad (15)$$

while the dynamic stress components are determined from the following difference:

$$\hat{\Sigma}_i^{\text{dyn}} = \hat{\Sigma}_i(t^+) - \hat{\Sigma}_i^{\text{sta}} \quad (16)$$

where $\hat{\Sigma}_i(t^+)$ is the macroscopic stress (sum of the static and the dynamic contributions) obtained at the end of the second stage of the numerical procedure. The macroscopic mean stresses are computed according to:

$$\hat{\Sigma}_m^{\text{sta}} = \frac{\hat{\Sigma}_1^{\text{sta}} + \hat{\Sigma}_2^{\text{sta}} + \hat{\Sigma}_3^{\text{sta}}}{3} \\ \hat{\Sigma}_m^{\text{dyn}} = \frac{\hat{\Sigma}_1^{\text{dyn}} + \hat{\Sigma}_2^{\text{dyn}} + \hat{\Sigma}_3^{\text{dyn}}}{3} \quad (17)$$

$$\hat{\Sigma}_m = \hat{\Sigma}_m^{\text{sta}} + \hat{\Sigma}_m^{\text{dyn}}$$

RVEs studied in this work have cubic symmetry. The macroscopic behaviors are identical along the three axial directions resulting in vanishing macroscopic equivalent stresses.

# 4. Simulation results

## 4.1. Unit cell results

We first investigate the response of the configurations UC1 and UC2 under spherical expansion. In this section, macroscopic quantities will be noted without the symbol ‘.’ above them. The boundary condition

<figure><img src="image_5.png" /><figcaption>Fig. 5. Illustration of the applied velocity $v_i$ at the external face $x_i = \pm h$ of the RVE, see Eq. (13). A two-stage simulation is performed with Abaqus/Explicit: a quasi-static calculation (without geometrical nonlinearities) up to time $t^+$, followed by a dynamic calculation (including geometrical nonlinearities) over a duration $t^* - t^+$. At the end of the quasi-static loading, the unit cell is nearly fully plastic. The resulting velocity and stress fields obtained from the first simulation are used as initial conditions for the non-linear dynamic second computation.</figcaption></figure>
(Eq. (13)) is considered substituting $\hat{D}_m$ by $D_m$ and assigning $h = l/2$. The static and dynamic contributions of the macroscopic stress are addressed and the velocity field inside each RVE is analyzed.

### 4.1.1. Static macroscopic stress

We compare the numerical results with the analytical expression provided by Tvergaard (1988) and inherited from the Gurson-Tvergaard (GT) flow surface:

$$ \Omega = \left(\frac{\Sigma_{\text{eq}}^{\text{sta}}}{\sigma_0}\right)^2 + 2q_1 f \cosh\left(\frac{3}{2} \frac{\Sigma_m^{\text{sta}}}{\sigma_0}\right) - 1 - (q_1 f)^2 = 0 \quad (18) $$

The parameter $q_1$ has been introduced by Tvergaard (1981) to ensure closer agreement, at low void volume fractions, between numerical results of periodic arrays of voids and predictions of the model. Under spherical expansion, the yield function Eq. (18) reduces to:

$$ \Sigma_m^{\text{sta}} = \frac{2}{3} \sigma_0 \ln\left(\frac{1}{q_1 f}\right) \quad (19) $$

which represents the critical stress at which the void grows without bound.

With $\sigma_0 = 500\text{MPa}$, numerical computations provide a mean stress ratio $\Sigma_m^{\text{sta}}/\sigma_0$ of 4.506 and 3.793 for UC1 and UC2 respectively, which are accurately compared to analytical values predicted by Eq. (19) with $q_1 = 1.16$. This value is adopted in the following to estimate the static stress of cubic unit cells under spherical loading. This calibration can be compared against values found in the literature (e.g. $q_1 = 1.08, 1.25, 1.46, 1.6, 2$ in Corigliano et al. (2000); Tvergaard (1981); Vadillo and Fernández-Sáez (2009); Kim et al. (2004); Nielsen and Tvergaard (2009), respectively). Note that under static loading conditions with $D_m > 0$ (expansion), the mean stress $\Sigma_m^{\text{sta}}$ is positive.

### 4.1.2. Dynamic macroscopic stress

Fig. 6 displays, in a dynamic stress ratio $\Sigma_m^{\text{dyn}}/\sigma_0$ versus strain rate $D_m$ diagram, the results obtained from finite element calculations (with $\sigma_0 = 500\text{MPa}$) along with the results predicted by the analytical model given by Eq. (7) with $b = 490\mu\text{m}$, and $a = 49\mu\text{m}$ for UC1, $a = 70\mu\text{m}$ for UC2. Recall that the dynamic stress $\Sigma_m^{\text{dyn}}$ is obtained from the difference

<figure><img src="image_6.png" /><figcaption>Fig. 6. Dynamic stress ratio $\Sigma_m^{\text{dyn}}/\sigma_0$ vs. strain rate for two unit cells of identical size. The matrix material yield strength is $\sigma_0 = 500 \text{ MPa}$. Finite Element results (symbols) were determined for UC1 and UC2 of identical size length ($l = 790 \mu\text{m}$), considering four values of the strain rate. Analytical results were obtained from the theory developed for spherical hollow spheres, Eq. (7) with $b = l(3/4\pi)^{1/3} = 490 \mu\text{m}$, and $a = 49 \mu\text{m}$ (dashed line), $a = 70 \mu\text{m}$ (solid line).</figcaption></figure>

between the macroscopic stress obtained at the end of dynamic calculation (second stage) and the stress resulted from the static calculation (end of the first stage). The value of $b$, representing the external spherical shell radius in the analytical approach, is determined as $b = l(3/4\pi)^{1/3}$, with $l = 790 \mu\text{m}$ as previously mentioned. Recall that the dynamic stress is independent of the yield stress of the matrix material, such that a dimensionless quantity is employed solely for the sake of representation.

In Fig. 6, it is shown that the Finite Element results are in very good agreement with the results obtained from the analytical model, considering Eq. (6) or, equivalently, Eq. (7). This observation demonstrates the quadratic influence of the mean strain rate $D_m$ on the dynamic stress of cubic unit cells under spherical loading. In addition, as elucidated in the analysis following Eq. (7), when comparing the dynamic stress of same size unit cells UC1 and UC2, one can observe the more pronounced micro-inertia effects experienced by UC1.

As an illustration of precise values for UC1, when $D_m = 13\ 856\text{s}^{-1}$, one has $\Sigma_m^{\text{dyn}}/\sigma_0 = -3.576$ from the theory which is compared to $-3.532$ obtained from numerical calculations. The same comparative analysis with UC2 provides $\Sigma_m^{\text{dyn}}/\sigma_0 = -0.847$ from the theory and $-0.850$ obtained from the finite element model. TAB. 1 gathers FEM results between the dynamic stress ratio and analytical results deduced from Eq. (6), for values of $D_m$ used for numerical simulation results reported in Fig. 6.

An important finding emerges from this comparison: for dynamic stress, there is no need to introduce a coefficient, such as the parameter $q_1$, which is used to adjust the static response of cubic RVEs to match that of spherical RVEs. Given this observation, the analytical hollow sphere model (Eq. (6) or equivalently Eq. (7)) can be used directly to predict the dynamic stress of a cubic unit cell, under spherical loading at constant strain rate ($\dot{D}_m = 0$).

Besides, a peculiar trend is observed when comparing static and dynamic loadings. The macroscopic stress, defined as the sum of a static and a dynamic part, is always positive under static expansion loading while under dynamic loading, the macroscopic stress may reveal negative values. This is a clear manifestation of micro-inertia, as also highlighted in Sartori et al. (2015). Specifically, using data from TAB. 1 and the results reported in section 4.1.1, one obtains $\Sigma_m^{\text{sta}}/\sigma_0 = \Sigma_m^{\text{dyn}}/\sigma_0 = 4.096$ for UC1 and $\Sigma_m^{\text{sta}}/\sigma_0 = 3.703$ for UC2 when $D_m = 4\ 619\text{s}^{-1}$, and $\Sigma_m^{\text{dyn}}/\sigma_0 = -1.76$ for UC1 and $\Sigma_m^{\text{dyn}}/\sigma_0 = 2.297$ for UC2 when $D_m = 18\ 475\text{s}^{-1}$.

### 4.1.3. Velocity field inside the unit cell

In this section, we conduct a comparison between the velocity field obtained from finite element calculations on UC1 and UC2 and the incompressible velocity field proposed by Gurson (1977) for the hollow sphere model. The goal is to clarify whether the velocity field adopted for quasi static loading is still valid under dynamic loading. For a general loading state, the trial velocity field proposed by Gurson (1977) is of the form:

$$ \underline{v} = \underline{D}'\underline{x} + D_m \left(\frac{b}{r}\right)^3 \underline{x} \quad (20) $$

where $\underline{D}'$ is the deviatoric part of the strain rate tensor, $b$ the external radius of the hollow sphere, and $r$ is the norm of the position vector $\underline{x}$. In the case of spherical loading, Eq. (20) simplifies to the following expression for the radial component:

$$ v_r = \frac{D_m}{f} \frac{a^3}{r^2} \quad (21) $$

where $a$ is the void radius, and $f = a^3/b^3$ is the porosity of the hollow sphere. When the loading is a spherical expansion ($D_m > 0$), the radial component coincides with the velocity norm. Note that this last relationship
provides the exact solution of a hollow sphere made of a fully incompressible matrix and subjected to a uniform remote strain rate $D_m$. It's important to recall that, except in the vicinity of the cubic unit cell's external faces, the matrix material is fully plastic at the time of our analysis. Thus, with a large elastic modulus, the elastic strain remains limited compared to the plastic contribution.

Fig. 7 shows a comparison between the norm of the velocity field obtained from Eq. (21), with $a \le r \le af^{-1/3}$, and the corresponding value calculated from FEM simulations with $D_m = 13\ 856s^{-1}$. The analysis is conducted at the end of the first stage of the computation procedure, knowing that the evolution of the velocity field from $t^+$ to $t^+$ is negligible. The velocity field around the void, as obtained from numerical computations, is displayed along three distinct paths all originating from the center of the cell (i.e., the center of the void). These paths are designated as Path 1 (from point O to point P1), Path 2 (from O to P2), and Path 3 (from O to P3), see Fig. 2. Along these paths, the position of a point varies as follows: from 0 to $l\sqrt{3}/2$ for Path 1, $l\sqrt{2}/2$ for Path 2, and $l/2$ for Path 3.

We observe that all responses overlap in the vicinity of the void up to a distance of approximately $l/4$. Discrepancies become apparent as we move further away from the void, which is attributed to the cubic shape of the unit cell. Consequently, a radial velocity field develops around the void. Furthermore, from this, by determining the velocity magnitude at the void boundary, one can estimate the amplitude of the strain rate at the remote boundary of a cubic cell, using Eq. (21).

## 4.2. RVE results: reference case$B_4S_4$

The reference RVE, denoted by $B_4S_4$ in Fig. 4, has been selected to illustrate the main outcomes of our investigation. This RVE contains four UC1 cells and four UC2 cells, arranged to ensure the RVE is well balanced, see Fig. 8. The overall porosity of the RVE, as deduced from Eq. (10) is $\hat{f} = 1.958\ 10^{-3}$. In this section, the remote strain rate is $\hat{D}_m = 13\ 856s^{-1}$.

Other configurations have been subjected to the same analysis. Results of all RVEs described in Fig. 4 are discussed in section 5.3.

### 4.2.1. Static macroscopic stress

Under quasi-static loading, the macroscopic stress is not depending

<figure><img src="image_7.png" /><figcaption>Fig. 8. Reference configuration denoted by $B_4S_4$. This RVE, made of an assembly of four UC1 and four UC2 unit cells, is well balanced (the center of mass coincides with the center of the cell and the RVE has cubic symmetry). The two arrows define paths from which the velocity will be extracted.</figcaption></figure>

on the intensity of the macroscopic strain rate. In addition, numerical simulations conducted here show that the macroscopic response of the RVE $B_4S_4$ under static loading align well with the response of a unit cell with porosity $\hat{f}$. In fact, the macroscopic hydrostatic stress computed from numerical simulation under static loading provided the value of 4.04 which is well compared to $\frac{\hat{\Sigma}_m^{\text{sta}}(\hat{f})}{\hat{\sigma}_0} = 4.06$ calculated using Eq. (19) with $q_1 = 1.16$.

### 4.2.2. Dynamic macroscopic stress

The dynamic macroscopic stress obtained from the FEM calculation

<figure><img src="image_8.png" /><figcaption>Fig. 7. Velocity at time $t^+$ vs. position for configurations (a) UC1 ($a = 49\ \mu\text{m}$, $f = 0.001$), and (b) UC2 ($a = 70\ \mu\text{m}$, $f = 0.002915$), subjected to a uniform remote strain rate $D_m = 13\ 856s^{-1}$. Comparison of the norm of the velocities along Path 1, Path 2, and Path 3 (see Fig. 2), obtained through finite element calculations, with the analytical velocity norm obtained from Eq. (21).</figcaption></figure>
is $\frac{\Sigma^{dyn}}{\sigma_0} = -1.560$ which lies between the values obtained for UC1 ($-3.532$) and UC2 ($-0.847$) under the same imposed strain rate, see section 4.1.2, Table 1. It appears that the overall dynamic response is closer to the value obtained in UC2, seemingly reflecting the concept of the weakest link, the response of the aggregate being strongly influenced by the cell with the highest porosity. Note that the volume average calculated from the two last values overestimates by 40% the results for $B_4S_4$, indicating that a rudimentary approach, based on isolated calculations conducted on unit cells, cannot provide a sound estimate for an aggregate.

The outcome $\frac{\Sigma^{dyn}}{\sigma_0} = -1.560$ determined for $B_4S_4$ is subjected to a comparative analysis using two distinct unit cell models. The first model considers a cubic unit cell, mirroring the external volume and porosity of the RVE $B_4S_4$. From the results of section 4.1.2, the analytical approach can be used to provide a very good estimate of the dynamic stress for this cubic unit cell. Using Eq. (7) with $f = \hat{f}$ and $b = 2l(3/4\pi)^{1/3}$, one gets $\frac{\Sigma^{dyn}}{\sigma_0} = -5.82$ which is almost four times the value obtained with the RVE $B_4S_4$.

The second model utilizes a unit cell with identical porosity $\hat{f}$ containing a void whose volume coincides with the average void volume of $B_4S_4$. This model corresponds to a unit cell of side length $l = 790$ µm embedding a void of radius 61.30 µm. In this case, the calculation yields $\frac{\Sigma^{dyn}}{\sigma_0} = -1.45$ providing a closer match to the numerical simulation result. Both scenarios illustrate the impact of void radius on the dynamic response of the material. More important, this analysis highlights how critical is the definition of an equivalent unit cell that can effectively replace an array of voids of different sizes within an RVE under dynamic loading.

### 4.2.3. Velocity field

An in-depth analysis of the velocity field generated within the RVE $B_4S_4$ is crucial for understanding the complex interactions between voids in this RVE. Fig. 9 shows color V-contour maps at the end of the first stage of the computation procedure, representing the amplitude of the particle velocity in the plane $x_3 = -l/2$ passing through the centers of four voids, see Fig. 8. The results show that the velocity field adopts a nearly spherical configuration in the vicinity of each void. As the distance from the void increases, the influence of neighboring voids or cell boundaries becomes more visible. Furthermore, one can notice that the magnitude of the velocity around big voids are larger than around small ones. Thus, the dynamic contribution that takes place inside regions delineated by unit cells UC1 is less significant, while unit cells UC2 contributes largely to the dynamic response of the RVE $B_4S_4$. This result is consistent with results of section 4.2.2.

To supplement this analysis, we define two node paths crossing the entire RVE, and passing through the centers of two successive voids from the small to the large void, see Fig. 8. Fig. 10a shows the evolution along Path 4 of $|v_{path4}|$ and $||v||$ respectively representing the absolute value of the velocity magnitude component along $x_1$, and the norm of the velocity vector. Data are extracted at time $t^+$. It appears that both evolutions coincide over almost the entire range $[-l, +l]$. The same trend is revealed in Fig. 10b, which shows evolutions along Path 5 of $||v||$, and

<caption>**Table 1**</caption>

FEM results compared to analytical calculations of $\Sigma^{dyn}/\sigma_0$ for UC1 ($a = 49\mu m$, $f = 0.001$), and UC2 ($a = 70\mu m$, $f = 0.002915$), subjected to a uniform remote strain rate $D_m$.

<table><thead><tr><th>Dm [s<sup>-1</sup>]</th><th>FEM results (UC1)</th><th>Eq. (6)</th><th>FEM results (UC2)</th><th>Eq. (6)</th></tr></thead><tbody><tr><td>4619</td><td>-0.410</td><td>-0.397</td><td>-0.090</td><td>-0.094</td></tr><tr><td>9238</td><td>-1.576</td><td>-1.589</td><td>-0.378</td><td>-0.377</td></tr><tr><td>13856</td><td>-3.532</td><td>-3.576</td><td>-0.847</td><td>-0.850</td></tr><tr><td>18475</td><td>-6.266</td><td>-6.357</td><td>-1.497</td><td>-1.511</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Distribution of the norm of the velocity (in $m/s$) inside the reference RVE $B_4S_4$ made of an assembly of four UC1 and four UC2 cells, with an overall porosity $\hat{f} = 1.958 10^{-3}$. The cutting plane $x_3 = -l/2$ is used, see Fig. 8. The results are depicted at the end of the first stage of the computation procedure. The remote strain rate is $\hat{D}_m = 13 856s^{-1}$.</figcaption></figure>

$|v_{paths}|$ the absolute value of the velocity vector projection over the direction (1,1,1). Some differences emerge at $x_1 = -l$ in Fig. 10a, where $||v|| = 18.727m/s$ is larger than the imposed velocity $|v_1| = \hat{D}_m = 10.946m/s$. This results from interactions with the external boundary and the presence of large voids positioned orthogonally to Path 4, inducing non-vanishing components $v_2$ and $v_3$. Another effect of void interactions is put in evidence in Figs. 10a and b where the velocity magnitude at the void boundaries of the small void fluctuates greatly. Depending on the location of the nodes, $||v||$ adopts the following values at the small void boundary: 229.94 $m/s$ and 284.11 $m/s$ (Fig. 10a), 194.75 $m/s$ and 302.48 $m/s$ (Fig. 10b). These absolute differences are more limited at the larger void boundaries where the velocity varies from 521.36 $m/s$ and 537.24 $m/s$ in Fig. 10a, and 518.69 $m/s$ and 556.67 $m/s$ in Fig. 10b.

It's important to note that, due to the positioning of the voids within $B_4S_4$, other node paths passing from a small to a large void and aligned with any principal axis or diagonal axis of the cube would yield identical results. As a matter of fact, Fig. 10 is representative of the velocity field for the two families of voids embedded inside the RVE.

### 4.2.4. Void-void interaction parameters

It is admitted that the velocity distribution around voids characterizes interactions inside the RVE. To gain insight into these interactions at the void level inside $B_4S_4$, particle velocities were extracted from all nodes on the surfaces of both a small void (radius $a_1$) and a large void (radius $a_2$). Subsequently, the following dimensionless parameter was calculated for each node $i$:

$$ \xi_i = \frac{v_i}{\hat{D}_m} \qquad (22) $$

Here, $v_i$ stands for the radial velocity at a specific node $i$ on the void surface and $\hat{D}_m$ represents the velocity imposed at the external boundary of the RVE. $\xi_i$ characterizes the effect of several interactions. To provide reference values for the interaction parameter in the case of UC1

$$ \text{and UC2, the following values are obtained: } \left( \frac{3}{4\pi f_1^2} \right)^{1/3} \approx 30.4 \text{ for UC1} $$

$$ \text{and } \left( \frac{3}{4\pi f_2^2} \right)^{1/3} \approx 62.0 \text{ for UC2. These results are deduced taking advantage of the very good predictive capability of the analytical model} $$
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Evolution (a) along Path 4 of the norm of the velocity vector $||v||$, and $|v_{path4}|$, the absolute value of the velocity component along $x_1$, and (b) along Path 5 of $||v||$ and $|v_{path5}|$, the absolute value of the velocity vector projection over the direction (1,1,1). Path 4 and Path 5 cross the entire RVE, and pass through the centers of two successive voids from the small to the large void, see Fig. 8. The macroscopic hydrostatic strain rate is $\hat{D}_m = 13\,856.4s^{-1}$. The results are depicted at the end of the first stage of the computation procedure within the reference RVE $B_4S_4$ consisting of an assembly of four UC1 and four UC2 cells, with an overall porosity $\hat{f} = 1.958\,10^{-3}$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 11.** (a) and (b) – Distribution of the normalized radial velocity $\epsilon_i = \frac{y_i}{|D_m|}$ under imposed strain rate $\hat{D}_m = 13856.4s^{-1}$, at nodes i located on the surface of the top two voids of Fig. 8: (a) the void of radius $a = 70\,\mu m$ and (b) the void of radius $a = 49\,\mu m$. (c) and (d) - Location of nodes of interest defined by its polar angle $\theta$ and azimuthal angle $\phi$ in a spherical coordinate system centered at the void center. The results correspond to the reference RVE $B_4S_4$ consisting of an assembly of four UC1 and four UC2 cells, with an overall porosity $\hat{f} = 1.958\,10^{-3}$.</figcaption></figure>
(developed for spherical unit cells) to reproduce the response of a cubic unit cell.

Fig. 11a and b illustrate the distribution of $\xi|_i$ for the top two voids of Fig. 9 which are representative of the two families of voids embedded inside the RVE. In this figure, the position of each node is determined by its polar angle $\theta$ ranging from 0 to $\pi$ and its azimuthal angle $\phi$ ranging from $-\pi$ to $\pi$ in a spherical coordinate system centered at each void center, as shown in Fig. 11c and 11d. Some particular nodes, located at the surface of each void, are defined in these figures to highlight specific interactions in distinct directions.

Values of $\xi|_A$ at point A and $\xi|_B$ at point B demonstrate how small voids are affected by the large voids in the diagonal and longitudinal directions, respectively. Correspondingly, $\xi|_A$ and $\xi|_B$ represent the effects of small voids on the large void in the aforementioned directions. Points C, D, and D' are associated with interactions between the external boundary of the RVE and voids.

From Fig. 11a and b, one observes that the heterogeneity in the distribution of $\xi|_i$ around the small void is more pronounced, with an average value of $\xi_1 = 23.38$ and a standard deviation of $\Delta\xi_1 = 2.74$, while the corresponding values for the large void are $\xi_2 = 48.46$ and $\Delta\xi_2 = 1.18$, respectively. Results obtained from other configurations, as shown in section 4.3, actually indicate a more complex scenario for cases where the number of small and large voids embedded in the RVE is not the same.

In order to better represent void-void interactions for each void, $\xi_1^*$ and $\xi_2^*$ are introduced for the small and large void, respectively, as the averages of $\xi|_i$ calculated over 1/8th of the surface of each void identified by the shaded portion in Fig. 11c and d. These areas are delineated by white rectangles in Fig. 11a and b, outlining the respective range of variation for $\theta$ and $\phi$. The following values are obtained: $\xi_1^* = 27.16$ with a standard deviation of 0.36, $\xi_2^* = 47.36$ with a standard deviation of 0.08.

### 4.3. RVE results

Similarly to the case $B_4S_4$, in the remaining RVEs, the distribution of $\xi|_i$ around one large and one small void is indicative of the velocity field for the two families of embedded voids. The analysis is conducted with $\hat{D}_m = 13\ 856s^{-1}$, focusing on the two top voids shown in Fig. 4, positioned with their centers on the plane $x_3 = -1/2$. Fig. 12 shows the distribution of $\xi|_i$ around the two aforementioned voids for all RVEs except for $B_0S_8$ and $B_8S_0$ which contain a unique family of voids.

The corresponding average and standard values of the interaction parameters $\xi_1$ and $\xi_1^*$ for the small void, and $\xi_2$ and $\xi_2^*$ for the large void, as defined in the previous section, are reported in Table 2. The dynamic stress ratios, obtained from finite element calculations, are also presented in Table 2.

As seen in Table 2 and it is clear that for RVEs out of $B_0S_8$ and $B_8S_0$ which contain a unique family of voids, $\xi_2$ is larger than $\xi_1$. Thus, the larger voids deform more significantly than the smaller ones. Moreover, as their number increases, the deformation becomes less localized around these larger voids. In fact, the values of $\xi_1$ and $\xi_2$ both decrease as the overall porosity increases. Given that $\xi_1$ and $\xi_2$ are associated with the average velocity around small and large voids respectively (e.g. see Eq. (22)), even though $\xi_2 > \xi_1$, when one replaces a large void by a small one, the overall porosity decreases, and the inertial contribution becomes more significant. This is confirmed by the values of $\frac{\xi_m^{dyn}}{\sigma_0}$ presented in Table 2.

It is interesting to consider the cases of $B_8S_0$ and $B_0S_8$, obtained from a stacking of eight identical unit cells UC1, or UC2. In fact, the response of each of those two stackings can be approximated based on the response of the corresponding unit cell, for which we derived the interaction parameters as $(\frac{3}{4\pi f_1^2})^{1/3} \approx 30.4$ for UC1 and $(\frac{3}{4\pi f_2^2})^{1/3} \approx$

62.0 for UC2. These values, which consider the interaction between voids of identical size through the unit cell porosity, closely match the average values of $\xi_1$ and $\xi_2$ as reported in Table 2 for $B_8S_0$ and $B_0S_8$, respectively.

Most important, our results enable the suggestion of a theoretical approach to characterize the dynamic behavior of RVEs containing two populations of voids, based on averaging the dynamic response of isolated unit cells. To illustrate this point, let us note that the value of $\xi_2$ obtained for $B_8S_0$ is lower than the corresponding values obtained for any cases where two populations of different voids are embedded in the RVE. This implies that higher values than $\hat{D}_m$ should be imposed at the level of UC2, to get similar amplitudes of the velocity at the boundary of the small void inside the RVE. Similarly, based on the same comparative analysis and reasoning applied to $\xi_1$, the strain rate to be imposed at the level of UC1 should be lower than $\hat{D}_m$.

The following section aims at proposing an alternative to a Taylor type averaging approach, which assumes that all unit cells experience a strain rate of same amplitude.

## 5. Proposed theoretical approach for porous materials with two populations of voids

Based on the numerical investigation detailed previously, an analytical model for the dynamic response of a porous material containing spherical voids of various sizes is derived in this section. More specifically, the analysis is restricted to the case of two spherical unit cells, denoted by $SC_i$ ($i=1, 2$). Each cell is characterized by its inner radius $a_i$ and porosity $f_i$, the outer radius being given by $b_i = a_i f_i^{-1/3}$.

### 5.1. Macroscopic response at the scale of the RVE and averaging

We consider a two-scale model, and the porous material is viewed as an aggregate of spherical unit cells of identical size. Elasticity is neglected and the matrix material is assumed incompressible. From Fig. 1, we recall that $(\hat{\Sigma}; \hat{D})$ are the macroscopic quantities at the RVE level. We adopt the following notation: $(\Sigma_i; D_i)$ are the corresponding quantities at the level of the unit cells $SC_i$. Obviously, the link between macroscopic quantities, i.e. $(\hat{\Sigma}; \hat{D})$ defined at the level of the RVE, and $(\Sigma_i; D_i)$ defined at the unit cell level must be characterized.

The RVE is subjected to a homogeneous strain rate boundary condition, $\nu = \hat{D} \cdot x$ on the outer boundary, with:

$$ \hat{D} = \langle D_i \rangle \qquad (23) $$

where $\langle \cdot \rangle = \frac{1}{V} \int dV$, with V the volume of the overall domain. With voids of different sizes, the porous material naturally exhibits a complex microstructure. Despite this complexity, we assume that all unit cells experience uniform loading conditions. In addition, the analysis of section 4.3 shows that the classical Taylor-type scheme where $D_i = \hat{D}$ cannot be considered as a valid option to reproduce finite element results. Here, we propose to generalize the approach by assuming that a linear relationship holds for a given unit cell $SC_i$:

$$ D_i = \alpha_i \hat{D} \qquad (24) $$

where $\alpha_i$ is a scalar parameter which may vary from one unit cell to another.

The following equation, obtained from combining Eq. (23) with Eq. (24), yields a relationship between the parameters $\alpha_i$:

$$ \langle \alpha_i \rangle = 1 \qquad (25) $$

Eq. (24) provides more flexibility in distributing the deformation in each sub-region. In agreement with Finite Element results, the
<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Distribution of the normalized radial velocity $\xi|_{i} = \frac{y_{i}}{l_{D_{m}}}$ around the two top right voids of Fig. 4, having their centers located on the plane $x_{3} = -l/2$; (a–d) the large void, (e–h) the small void embedded in the RVEs $B_{0}S_{8}$, $B_{2}S_{6}$, $B_{4}S_{4}$, $B_{6}S_{2}$, and $B_{8}S_{0}$ under imposed strain rate $\hat{D}_{m} = 13\,856.4s^{-1}$. $\theta$ and $\phi$ are defined similarly to Fig. 11c and d.</figcaption></figure>
<caption>**Table 2**</caption>

Average and standard values of the interaction parameters $\xi_1$ (small void of radius $a = 49 \ \mu m$) and $\xi_2$ (large void of radius $a = 70 \ \mu m$) over the entire domain, and corresponding $\xi_1^*$ and $\xi_2^*$ calculated over 1/8th of the surface of each void corresponding to white rectangles in Fig. 12.

<table><thead><tr><th></th><th>f̂ (x 10<sup>-3</sup>)</th><th>Σ̂<sub>m</sub><sup>dyn</sup> /σ<sub>0</sub></th><th>ξ<sub>1</sub></th><th>Δξ<sub>1</sub></th><th>ξ<sub>2</sub></th><th>Δξ<sub>2</sub></th><th>ξ<sub>1</sub><sup>*</sup></th><th>Δξ<sub>1</sub><sup>*</sup></th><th>ξ<sub>2</sub><sup>*</sup></th><th>Δξ<sub>2</sub><sup>*</sup></th></tr></thead><tbody><tr><td>B<sub>0</sub>S<sub>8</sub></td><td>1.000</td><td>-3.565</td><td>60.25</td><td>0.51</td><td>/</td><td>/</td><td>59.98</td><td>0.29</td><td>/</td><td>/</td></tr><tr><td>B<sub>2</sub>S<sub>6</sub></td><td>1.479</td><td>-2.203</td><td>37.30</td><td>1.39</td><td>64.40</td><td>1.21</td><td>37.62</td><td>0.82</td><td>63.31</td><td>0.13</td></tr><tr><td>B<sub>4</sub>S<sub>4</sub></td><td>1.958</td><td>-1.560</td><td>23.38</td><td>2.73</td><td>48.46</td><td>1.18</td><td>27.16</td><td>0.36</td><td>47.36</td><td>0.08</td></tr><tr><td>B<sub>6</sub>S<sub>2</sub></td><td>2.436</td><td>-1.320</td><td>2.53</td><td>0.22</td><td>39.69</td><td>3.47</td><td>2.18</td><td>0.09</td><td>39.36</td><td>2.06</td></tr><tr><td>B<sub>8</sub>S<sub>0</sub></td><td>2.915</td><td>-0.856</td><td>/</td><td>/</td><td>30.12</td><td>0.50</td><td>/</td><td>/</td><td>29.35</td><td>0.05</td></tr></tbody></table>

parameter $\alpha_i$ is introduced at the level of unit cells to account for void-void interactions. Such interactions have been characterized at a local scale in section 4.2.4 by $\xi_i^*$, the radial velocity at the void surface normalized by the outer expansion velocity. Note that the classical Taylor-type scheme, as adopted in Czarnota et al. (2008), Jacques et al. (2010, 2012), is retrieved when $\alpha_i = 1$.

Similar to the dynamic response of an individual unit cell, see Eq. (3), the following stress decomposition holds at the level of the RVE:

$$ \hat{\Sigma} = \hat{\Sigma}^{\text{sta}} + \hat{\Sigma}^{\text{dyn}} \quad (26) $$

where $\hat{\Sigma}^{\text{sta}}$ is the static part of the macroscopic stress considering the overall porosity $\hat{f}$ of the RVE being analyzed. $\hat{\Sigma}^{\text{dyn}}$ represents the dynamic part obtained through volume averaging over the RVE:

$$ \hat{\Sigma}^{\text{dyn}} = \theta_1 \hat{\Sigma}_1^{\text{dyn}} + \theta_2 \hat{\Sigma}_2^{\text{dyn}} \quad (27) $$

where $\theta_1$ and $\theta_2 = 1 - \theta_1$ have been defined in Eq. (11).

The analysis now is restricted to the case of hydrostatic loading. Therefore, Eqs. 23 and 24 yield:

$$ \hat{D}_m = \langle D_{im} \rangle \quad (28a) $$

$$ D_{im} = \alpha_i \hat{D}_m \quad (28b) $$

In addition, we consider that all unit cells are of the same size, i.e. $b = b_1 = b_2$. Recall that $\hat{D}_m$ is constant, so $\hat{D}_m$ is neglected in Eq. (4). Making use of Eqs. (7, 28b), Eq. (27) thus reduces to:

$$ \hat{\Sigma}_{m}^{\text{dyn}} = \rho b^2 \hat{D}_{m}^{2} [\theta_{1} \alpha_{1}^{2} G(f_{1}) + \theta_{2} \alpha_{2}^{2} G(f_{2})], \text{ with } G(f_{i}) = 3f_{i}^{-1/3} - \frac{5}{2} - \frac{1}{2} f_{i}^{-\frac{4}{3}} \quad (29) $$

A quadratic dependence on $\alpha_i$ is proposed in Eq. (29). This emphasizes the necessity of accurately evaluating void interactions when micro inertia is involved in the mechanical behavior of the porous material.

### 5.2. $\alpha_i$ parameters identification:

The parameters $\alpha_1$ and $\alpha_2$ are related through Eq. (25). Therefore, an additional equation is needed to derive an analytical model in the case of a material containing two families of spherical voids.

Recall that the interaction parameters $\alpha_i$ are introduced to adjust the strain rate at the level of each unit cell, see Eq. (28b), with the aim to reflect void-void interactions. It appeared, from numerical simulation results, that void-void interaction induces a heterogeneous velocity field characterized by the parameters $\xi_1$ and $\xi_2$ in section 4.2.4. In addition, the velocity field inside the RVE only slightly evolves over the time period $t^* - t^+$ of the dynamic computation stage. Therefore, one can assume that considering the case of a porous material under static loading informs about void-void interactions at high strain rates, at least for the considered loading conditions and material (constant strain rate, perfectly plastic matrix material). This is consistent with analytical approaches founded on a dynamic homogenization scheme for which static velocity fields are employed to derive the macroscopic response of porous materials under dynamic loading (Molinari and Mercier, 2001; Sartori et al., 2015, 2016; Subramani et al., 2020).

In the context of spherical loading, the rate of plastic work within the

RVE under static conditions is obtained from volume average of plastic power per spherical shells,

$$ \hat{\Sigma}_{m}^{\text{sta}}(\hat{f}) \hat{D}_{m} = \theta_{1} \Sigma_{m}^{\text{sta}}(f_{1}) D_{1m} + \theta_{2} \Sigma_{m}^{\text{sta}}(f_{2}) D_{2m} \quad (30) $$

Using Eq. (28b), one obtains:

$$ \hat{\Sigma}_{m}^{\text{sta}}(\hat{f}) = \theta_{1} \alpha_{1} \Sigma_{m}^{\text{sta}}(f_{1}) + \theta_{2} \alpha_{2} \Sigma_{m}^{\text{sta}}(f_{2}) \quad (31) $$

$\Sigma_m^{\text{sta}}(f_1)$, $\Sigma_m^{\text{sta}}(f_2)$, and $\hat{\Sigma}_m^{\text{sta}}(\hat{f})$ are determined from the analytical solution (19) derived for the hollow sphere model under spherical loading. Note that numerical simulation results of section 4.2.1 showed that $\hat{\Sigma}_m^{\text{sta}}(\hat{f})$ provides a good estimate of the macroscopic static response of the RVE made of two populations of voids.

Finally, from a combination of Eqs. (19), (25) and (31), one gets the theoretical expressions of $\alpha_1$ and $\alpha_2$, when $\theta_1 \in ]0; 1[$:

$$ \alpha_1 = \frac{1}{\theta_1} \left[ \frac{\ln(\hat{f}) - \ln(f_2)}{\ln(f_1) - \ln(f_2)} \right], \quad \alpha_2 = \frac{1 - \alpha_1 \theta_1}{1 - \theta_1} \quad (32) $$

with $\hat{f}$ the overall porosity of the RVE, $f_1$ the porosity of the unit cell SC1, $f_2$ the porosity of the unit cell SC2, $\theta_1$ the volume fraction of SC1 and $\theta_2$ the volume fraction of SC2. Once the parameters $\alpha_i$ are obtained, the dynamic macroscopic stress can be determined from Eq. (29). Obviously, Eq. (32) is valid when the RVE contains two populations of voids. In the case where the RVE is made of a single population of voids, one has $D_{im} = \hat{D}_m$. One observes from Eq. (32) that in the limiting case of a minimal presence of SC1, $\lim_{\theta_1 \to 0} \alpha_1 = (f_1/f_2 - 1) / \ln(f_1/f_2)$ is non-zero and $\lim_{\theta_1 \to 0} \alpha_2 = 1$. Likewise, when the proportion of SC2 is negligible, one has $\lim_{\theta_1 \to 1} \alpha_1 = 1$ and $\lim_{\theta_1 \to 1} \alpha_2 = (f_2/f_1 - 1) / \ln(f_2/f_1)$.

We propose to compare results from the proposed model (denoted $\hat{\Sigma}_m^{\text{dyn}}$ in the following) and estimates from the approach based on a Taylor type scheme (denoted $\hat{\Sigma}_m^{\text{dyn}}|_T$ in the following) obtained using Eq. (29) with $\alpha_1 = \alpha_2 = 1$. The difference between the two predictions $\Delta_1 = \frac{1}{\sigma_0} (\hat{\Sigma}_m^{\text{dyn}} - \hat{\Sigma}_m^{\text{dyn}}|_T)$, as well as the relative difference $\Delta_2 = \frac{\hat{\Sigma}_m^{\text{dyn}} - \hat{\Sigma}_m^{\text{dyn}}|_T}{|\hat{\Sigma}_m^{\text{dyn}}|}$ are investigated. The following material parameters are adopted $\sigma_0 = 500 \text{ MP a}$, $\rho = 7800 \text{ kg/m}^3$, and the macroscopic strain rate is $\hat{D}_m = 13 \ 856 \text{ s}^{-1}$. Various cases are considered for a fixed value of the void radius in SC1, $a_1 = 49 \ \mu m$ and three different values of the void radius characterizing the second population of voids: $a_2 = 70, 80$ and $90 \ \mu m$. Note that the relative difference $\Delta_2$ depends solely on the porosities of SC1 and SC2, as well as on $\theta_1$ while $\Delta_1$ additionally depends on the material parameters and strain rate level.

In Fig. 13a and c, the outer radius of the spherical cells SC1 and SC2 is fixed at $b = 490 \ \mu m$ resulting in a fixed porosity of $f_1 = 1.10^{-3}$ for SC1. The porosity for SC2 is $f_2 = 2.91 \ 10^{-3}$, $4.35 \ 10^{-3}$ and $6.19 \ 10^{-3}$ when $a_2 = 70, 80$ and $90 \ \mu m$ respectively. The overall porosity is given by $\hat{f} = v_1 f_1 + (1 - v_1) f_2$ which decreases with $v_1$. In Fig. 13b and d, where $v_1 =$
<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Evolution of $\Delta_1$, the difference in the dynamic stress ratio between the proposed model and the Taylor approach, in (a) with respect to $\vartheta_1$ using $b = 490 \ \mu m$, in (b) with respect to the cell radius $b$ for $\vartheta_1 = 0.5$. The corresponding relative difference $\Delta_2$ is portrayed in (c) and (d). The comparison is made for three assemblies with a fixed radius of the void in SC1, $a_1 = 40 \ \mu m$, and $a_2 = 70 \ \mu m$, $80 \ \mu m$, $90 \ \mu m$ for the larger void embedded in SC2. An alternative representation is proposed in Fig. B1 of Appendix B where the evolutions of $\Delta_1$ and $\Delta_2$ are shown with respect to the overall porosity $\hat{f}$.</figcaption></figure>

0.5, all three porosities $f_1$, $f_2$, and $\hat{f}$ decrease as $b$ increases. As an alternative representation of the results, Fig. B1 in Appendix B shows the evolution of $\Delta_1$ and $\Delta_2$ with respect to the overall porosity $\hat{f}$.

Fig. 13 demonstrates that the gap between the two approaches increases with the size of large voids $a_2$. The relative difference between the two estimates can be significant, exceeding 100%, as illustrated by the case when $a_2 = 90 \ \mu m$ at $\vartheta_1 = 0.4$ in Fig. 13c. This translates to about 1 GPa in $\Delta_1$ as shown in Fig. 13a. In the three scenarios analyzed in Fig. 13a and b, as the void radius $a_2$ increases, the value of $\vartheta_1$ at which $\Delta_1$ reaches its peak value increases too. Conversely, with the increase in $a_2$, the maximum values of $\Delta_2$ are obtained with a reduction in $\vartheta_1$, as shown in Fig. 13b.

The case $\vartheta_1 = 0.5$ is now selected to examine the influence of the external shell radius $b$. Fig. 13b shows that the difference between the two approaches, $\Delta_1$, becomes quite important when the outer radius of the unit cells is larger than $300 \ \mu m$ (i.e. $f_1 < 0.004$). This finding emphasizes the value of the proposed model, particularly when dealing with porous materials featuring small voids and when the disparity between the radii of the two populations is large.

In the next section, the theoretical model is tested for configurations examined previously by means of numerical simulations.

## 5.2. Comparison of analytical model and numerical results

The different cases of Fig. 4 are considered for a comparative analysis. Table 3 summarizes the values of $\Sigma_m^{dyn}/\sigma_0$ obtained in finite element calculations for the considered assemblies, as well as the corresponding values calculated from the new proposed approach. As in FEM calculations, the tensile yield stress is $\sigma_0 = 500$ MPa and the matrix mass density is $\rho = 7800 \text{kg/m}^3$. The macroscopic strain rate in FEM and in the analytical model is $\hat{D}_m = 13 \ 856 \text{s}^{-1}$.

The variation of $\Sigma_m^{dyn}/\sigma_0$ with respect to $\vartheta_1$, as predicted by the analytical model in Eq. (29), is illustrated in Fig. 14. Additionally, the numerical simulation results from Table 3 are reported using symbols. A close agreement between the FEM results and the theoretical estimates is observed from Table 3. The most significant difference occurs in the case of $B_6S_2$ ($\vartheta_1 = 0.25$), with the theoretical value underestimating by approximately 11% the value obtained from the numerical simulation. Table 3 shows that $\alpha_1$ and $\alpha_2$ both increase with the overall porosity $\hat{f}$. It is interesting to note that similar trends were revealed for the interaction parameters $\xi_1$ and $\xi_2$ which were defined from the velocity distribution around the voids in finite element simulations (see Table 2).

The model presented in Eq. (29) with $\alpha_1$ and $\alpha_2$ given by Eq. (32) accurately captures the non-linearity observed in FEM numerical results,
<figure><img src="image_14.png" /><figcaption>**Fig. 14.** Comparison between finite elements results (FEM) for the dynamic contribution $\hat{\Sigma}_m^{dyn}$ normalized by $\sigma_0$ (tensile yield stress of the matrix material) and analytical results (Eq. (29)) considering various void arrangements $v_1$. Results obtained from the Taylor homogenization scheme (Eq. (29)) with $\alpha_1 = \alpha_2 = 1$) and the unit cell of equivalent radius $(a^3)^{1/3}$ are also shown for comparison. The porosity $\hat{f}$ decreases when $v_1$ increases according to $\hat{f} = v_1f_1 + (1 - v_1)f_2$, with $f_1 = 1.10^{-3}$, $f_2 = 2.915 \ 10^{-3}$.</figcaption></figure>

<caption>**Table 3**</caption>

FEM results for tested assemblies composed of various arrangements of UC1 ($a_1 = 49\mu m, f_1 = 0.001$) and UC2 ($a_2 = 70\mu m, f_2 = 0.002915$) compared to analytical calculations of $\hat{\Sigma}_m^{dyn}/\sigma_0$ given by Eq. (29) with $f_1 = 0.001, f_2 = 0.002915, b = 490\mu m$. The case $B_8S_0$ (resp. $B_0S_8$) is compared to the analytical model using $\theta_1 \sim 0$ (resp. $\theta_1 \sim 1$). The matrix material properties are: $\sigma_0 = 500$ MPa and $\rho = 7800$ kg/m$^3$. The macroscopic strain rate is $\hat{D}_m = 13 \ 856s^{-1}$.

<table><thead><tr><th rowspan="2">RVE</th><th rowspan="2">f&#x0302; (x 10<sup>-3</sup>)</th><th colspan="2">FEM</th><th colspan="4">THEORY</th></tr><tr><th>&Sigma;&#x0302;<sub>m</sub><sup>dyn</sup>/&sigma;<sub>0</sub></th><th>&theta;<sub>1</sub></th><th>&theta;<sub>2</sub></th><th>&alpha;<sub>1</sub><br>Eq. (32)</th><th>&alpha;<sub>2</sub><br>Eq. (32)</th><th>&Sigma;&#x0302;<sub>m</sub><sup>dyn</sup>/&sigma;<sub>0</sub><br>Eq. (29)</th></tr></thead><tbody><tr><td>B<sub>0</sub>S<sub>8</sub></td><td>1.000</td><td>-3.565</td><td>1</td><td>0</td><td>1</td><td>/</td><td>-3.576</td></tr><tr><td>B<sub>2</sub>S<sub>6</sub></td><td>1.479</td><td>-2.203</td><td>0.75</td><td>0.25</td><td>0.846</td><td>1.462</td><td>-2.375</td></tr><tr><td>B<sub>4</sub>S<sub>4</sub></td><td>1.958</td><td>-1.560</td><td>0.5</td><td>0.5</td><td>0.744</td><td>1.255</td><td>-1.662</td></tr><tr><td>B<sub>6</sub>S<sub>2</sub></td><td>2.436</td><td>-1.320</td><td>0.25</td><td>0.75</td><td>0.670</td><td>1.110</td><td>-1.188</td></tr><tr><td>B<sub>8</sub>S<sub>0</sub></td><td>2.915</td><td>-0.856</td><td>0</td><td>1</td><td>/</td><td>1</td><td>-0.850</td></tr></tbody></table>

whereas the Taylor model exhibits linear dependence to $\theta_1$, as predicted by Eq. (29) using $\alpha_1 = \alpha_2 = 1$. Therefore, the Taylor model deviates significantly from the FEM results, with a relative difference about 42% observed for $B_4S_4$. In contrast, the proposed model achieves results closer to the FEM findings, with a relative difference of 6% for the same RVE.

In dynamic applications that involve strong micro inertia effects, the $\alpha_i$ parameters are key factors to adjust the inertial contribution at the level of each void. The inherent limitations of the Taylor model become apparent when attempting to account for this aspect. Consequently, employing the Taylor model for materials with varying void sizes under dynamic loads poses a significant challenge. As illustrated in the previous section, this concern becomes even more critical when the size of the small voids within a family is small and the difference with the radius of the bigger voids is large.

Another noteworthy aspect is the comparison between FEM results and an equivalent hollow sphere model. In this model, the void radius is obtained by averaging void volumes while preserving the overall porosity. Specifically, this implies considering a hollow sphere containing a void with radius $a = (\theta_1 a_1^3 + \theta_2 a_2^3)^{1/3}$.

Even if the equivalent hollow sphere model seems to provide with close agreement finite element results of several configurations, the

value of $\hat{\Sigma}_m^{dyn}/\sigma_0$ for the case $B_6S_2$ is far from the actual finite element result. Moreover, it is important to highlight that it does not explicitly account for the interactions elucidated in the finite element calculations. In contrast, our model, by introducing parameters $\alpha_i$, enables the consideration of interaction effects.

## 6. Conclusion

In this paper, a computational homogenization analysis has been carried out on three-dimensional representative volume elements (RVEs) to explore the effects of void interactions in porous materials containing two populations of voids under dynamic loading. The matrix material has elastic perfectly plastic properties. The spatial distribution of voids inside the RVE has been chosen to ensure equilibrium and cubic symmetry. We examined five microstructure models, composed of eight cubic unit cells containing a central spherical void, to assess the impact of void distribution on macroscopic stress responses.

Following the approach of Sartori et al. (2015), a two-stage simulation procedure was implemented to prevent strong perturbations that may arise when a large velocity is applied instantaneously. The velocity field inside the RVEs and the dynamic macroscopic stress were obtained considering a spherical strain rate loading with constant strain rate applied at the outer boundary of the RVEs.

The numerical analysis of unit cells confirmed that under static loading, an adjusting parameter is needed to correct the discrepancy between the analytical hollow sphere model and numerical results. It is noteworthy that our results show that there is no need to introduce such an adjusting parameter for the dynamic macroscopic stress.

A meticulous analysis of the velocity within the assemblies was conducted, aiming to understand the interaction between voids. The findings revealed that the velocity fields around each void was nearly spherical. In particular, it is demonstrated that material particles located near the larger voids experience more intense velocities compared to those situated around smaller voids. As a consequence, bigger voids largely contribute to the dynamic response of the porous material. A dimensionless parameter, defined from the radial velocity at the void surface normalized by the outer expansion velocity was introduced to highlight the effect of void interaction at the local scale, paving the way to propose an analytical approach.

Based on numerical results, the proposed analytical model involves an adjustment of the strain rate at the level of unit cells. This is achieved by interaction parameters defined at the larger scale of a spherical unit cell. The proposed model is validated with comparison against finite element results. It is also highlighted that the Taylor scheme is not suitable to describe a porous material containing two families of voids when micro inertia effect cannot be neglected. We showed that this issue is more critical when small voids are embedded in the porous material (small local porosity) and when the difference between the pore sizes of the two families is large.

In this work, we restricted our attention to the case of purely hydrostatic loading of well-balanced RVEs made of two types of unit cells. This type of loading has revealed, for porous materials containing spherical voids, the most significant contribution of micro-inertia effects (see Sartori et al., 2015). However, the Finite Element model proposed in this paper could be used to analyze other loading configurations by adjusting the form of the macroscopic strain rate tensor (Eq. (12)) and subsequently the boundary conditions (Eq. (13)). This will be the subject of future investigations. Additionally, considering only two scalar parameters to account for interactions are well adapted for the examined configurations. In fact, a unique void is representative of all voids belonging to the same population. Further development of the model would be to consider more complex arrangements of two families of voids, wherein the interaction parameter should account not only for the size, the porosity and the volume fraction of void populations, but also for variations in the voids' surroundings within a given population. One
potential enhancement of the model could involve expressing the interaction parameters $\alpha_i$ in a tensorial form.

### CRediT authorship contribution statement

**Mohammed El Ansı:** Writing – review & editing, Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis, Conceptualization. **Cédric Sartori:** Writing – review & editing, Visualization, Validation, Supervision, Software, Project administration, Methodology, Investigation, Formal analysis, Conceptualization.

**Christophe Czarnota:** Writing – review & editing, Visualization, Validation, Supervision, Software, Project administration, Methodology,

Investigation, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

## Appendix A: time periods adopted for nonlinear dynamic calculations

This appendix compiles the time periods utilized in our finite element procedure. A two-stage simulation is performed with Abaqus/Explicit: a quasi-static calculation (without geometrical nonlinearities) up to time $t^+$, followed by a dynamic calculation (including geometrical nonlinearities) over a duration $t^* - t^+$. The choice of the value for the time $t^+$ has no consequence on the results, and an arbitrary value was adopted. During the very short time period $t^* - t^+$, the body configuration has slightly evolved, with the porosity increasing by 2% of the initial value. TAB. A1 presents the time periods used for calculations conducted at various strain rates for unit cells UC1 and UC2 described in Fig. 2. TAB. A2 lists the time periods for RVEs depicted in Fig. 4 under an imposed strain rate of $\hat{D}_m = 13\ 856\ s^{-1}$.

### TAB. A1

Values of the time period $t^* - t^+$ adopted for nonlinear dynamic calculations (second stage of the numerical procedure, see Fig. 5) on cubic unit cells UC1 and UC2 described in Fig. 2. The body configuration slightly evolves during this very short time period, with the porosity increasing by 2% of the initial value.

<table><thead><tr><th></th><th>UC1<br>$f_1 = 0.001$<br>$a_1 = 49\mu m$</th><th>UC2<br>$f_2 = 0.002915$<br>$a_2 = 70\mu m$</th></tr><tr><th>$Dm$ [s<sup>-1</sup>]</th><th>$t^* - t^+$ [ns]</th><th>$t^* - t^+$ [ns]</th></tr></thead><tbody><tr><td>4619</td><td>1.44</td><td>4.22</td></tr><tr><td>9238</td><td>0.72</td><td>2.11</td></tr><tr><td>13856</td><td>0.48</td><td>1.40</td></tr><tr><td>18475</td><td>0.36</td><td>1.05</td></tr></tbody></table>

### TAB. A2

Values of the time period $t^* - t^+$ adopted for nonlinear dynamic calculations (second stage of the numerical procedure, see Fig. 5) on RVEs $B_0S_8$

$B_2S_6$, $B_4S_4$, $B_6S_2$, and $B_8S_0$ described in Fig. 4, under an imposed strain rate of $\hat{D}_m = 13\ 856\ s^{-1}$. The body configuration slightly evolves during this very short time period, with the porosity increasing by 2% of the initial value.

<table><thead><tr><th></th><th>$\hat{f}$<br>x 10<sup>-3</sup></th><th>$t^* - t^+$ [ns]</th></tr></thead><tbody><tr><td>$B_0S_8$</td><td>1.000</td><td>0.48</td></tr><tr><td>$B_2S_6$</td><td>1.479</td><td>0.71</td></tr><tr><td>$B_4S_4$</td><td>1.958</td><td>0.94</td></tr><tr><td>$B_6S_2$</td><td>2.436</td><td>1.17</td></tr><tr><td>$B_8S_0$</td><td>2.915</td><td>1.40</td></tr></tbody></table>

## Appendix B. Alternative representation of the differences between our proposed model and the Taylor approach

We propose in this appendix an alternative representation of the comparative analysis conducted in section 5.2 between the proposed model and estimates from the approach based on a Taylor type scheme. Two variables are investigated in section 5.2: the difference between the two predictions,

$$ \Delta_1 = \frac{1}{\alpha_0} (\hat{\Sigma}_m^{dyn} - \hat{\Sigma}_m^{dyn}|_T), \text{ and the relative difference } \Delta_2 = \frac{\hat{\Sigma}_m^{dyn} - \hat{\Sigma}_m^{dyn}|_T}{|\hat{\Sigma}_m^{dyn}|} $$

Fig. B1 provides an alternate representation of the results reported in Fig. 13. Here the evolutions of $\Delta_1$ and $\Delta_2$ are shown with respect to the overall porosity. Fig. 13a and c, the porosity in SC1 is fixed $f_1 = 1.10^{-3}$ and the overall porosity $\hat{f} = v_1f_1 + (1 - v_1)f_2$ with $f_2 = 2.91\ 10^{-3}$, 4.35 $10^{-3}$ and 6.19 $10^{-3}$ when $a_2 = 70$, 80 and 90 $\mu m$ respectively. In Fig. 13b and d, where $v_1 = 0.5$, both $f_1 = a_1^3/b^3$ and $f_2 = a_2^3/b^3$ decrease as $b$ increases. The same applies to the overall porosity given in this case by $\hat{f} = (a_1^3 + a_2^3)/2b^3$.
<figure><img src="image_15.png" /><figcaption>**Fig. B1.** Evolution of $\Delta_1$, the difference in the dynamic stress ratio between the proposed model and the Taylor approach with respect to the overall porosity $\hat{f} = v_1f_1 + (1 - v_1)f_2$ in (a) at fixed $b = 490 \ \mu m$, in (b) at fixed $\vartheta_1 = 0.5$. The corresponding relative differences $\Delta_2$ are portrayed in (c) and (d). The comparison is made for three assemblies with a fixed radius of the void in SC1, $a_1 = 40 \ \mu m$, and $a_2 = 70 \ \mu m$, $80 \ \mu m$, $90 \ \mu m$ for the larger void embedded in SC2.</figcaption></figure>

References

Antoun, T., Seaman, L., Curran, D., Kanel, G.I., Razorenov, S.V., Utkin, A.V., 2003. Spall Fracture. Springer-Verlag.

Barthélémy, R., Jacques, N., Kerampran, S., Vermeersch, F., 2016. Modelling of micro-inertia effects in closed-cell foams with application to acoustic and shock wave propagation. Int. J. Solids Struct. (97-98), 445-457.

Benzerga, A.A., Leblond, J.-B., 2010. Ductile fracture by void growth to coalescence. Adv. Appl. Mech. 44, 169-305.

Bilger, N., Auslender, F., Bornert, M., Moulinec, H., Zaoui, A., 2007. Bounds and estimates for the effective yield surface of porous media with a uniform or a nonuniform distribution of voids. Eur. J. Mech. A Solids. 26 (5), 810-836.

Carroll, M.M., Holt, A., 1972. Static and dynamic pore-collapse relations for ductile porous materials. J. Appl. Phys. 43 (4), 1626-1636.

Castañeda, P.P., 1991. The effective mechanical properties of nonlinear isotropic composites. J. Mech. Phys. Solids. 39 (1), 45-71.

Cazacu, O., Revil-Baudard, B., 2021. Tension-compression asymmetry effects on the plastic response in bending: new theoretical and numerical results. Mech. Res. Commun. 114, 103596.

Corigliano, A., Mariani, S., Orsatti, B., 2000. Identification of gurson-tvergaard material model parameters via kalman filtering technique. i. theory. Int. J. Fract. 104, 349-373.

Cortes, R., 1992. The growth of microvoids under intense dynamic loading. Int. J. Solids Struct. 29 (11), 1339-1350.

Czarnota, C., Molinari, A., Mercier, S., 2017. The structure of steady shock waves in porous metals. J. Mech. Phys. Solids. 107, 204-228.

Czarnota, C., Jacques, N., Mercier, S., Molinari, A., 2008. Modelling of dynamic ductile fracture and application to the simulation of plate impact tests on tantalum. J. Mech. Phys. Solids. 56, 1624-1650.

Czarnota, C., Molinari, A., Mercier, S., 2020. Steady shock waves in porous metals: viscosity and micro-inertia effects. Int. J. Plast. 135, 102816.

Espeseth, V., Morin, D., Faleskog, J., Børvik, T., Hopperstad, O.S., 2021. A numerical study of a size-dependent finite-element based unit cell with primary and secondary voids. J. Mech. Phys. Solids. 157, 104493.

Ferreira, A.R., Proença, S.P.B., Benallal, A., 2023. Yield criteria for voided materials with anisotropic matrix behavior. Eur. J. Mech. A Solids., 105079

Flandi, L., Leblond, J.-B., 2005. A new model for porous nonlinear viscous solids incorporating void shape effects-i: theory. Eur. J. Mech. A Solids. 24 (4), 537-551.

Fritzen, F., Forest, S., Böhlke, T., Kondo, D., Kanit, T., 2012. Computational homogenization of elasto-plastic porous metals. Int. J. Plast. 29, 102-119.

Garajeu, M., Michel, J., Suquet, P., 2000. A micromechanical approach of damage in viscoplastic materials by evolution in size, shape and distribution of voids. Comput. Methods Appl. Mech. Eng. 183 (3-4), 223-246.

Gäräjeu, M., Suquet, P., 1997. Effective properties of porous ideally plastic or viscoplastic materials containing rigid particles. J. Mech. Phys. Solids 45, 873-902.

Gologanu, M., Leblond, J.-B., Devaux, J., 1993. Approximate models for ductile metals containing non-spherical voids—case of axisymmetric prolate ellipsoidal cavities. J. Mech. Phys. Solids. 41 (11), 1723-1754.

Gologanu, M., Leblond, J.-B., Devaux, J., 1994. Approximate Models for Ductile Metals Containing Nonspherical Voids—Case of Axisymmetric Oblate Ellipsoidal Cavities. J. Eng. Mater. Technol. 116, 290-297.

Gurson, A.L., 1977. Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I—Yield Criteria and Flow Rules for Porous Ductile Media. J. Eng. Mater. Technol. 99 (1), 2-15.

Hosseini, N., Fuentes, J.N., Dakshinamurthy, M., Martinez, J.R., Vadillo, G., 2022. The effect of material orientation on void growth. Int. J. Plast. 148, 103149.

Jacques, N., Czarnota, C., Mercier, S., Molinari, A., 2010. A micromechanical constitutive model for the dynamic damage and fracture of ductile materials. Int. J. Fract. 162 (1), 159-175.

Jacques, N., Mercier, S., Molinari, A., 2012. Effects of microscale inertia on dynamic ductile crack growth. J. Mech. Phys. Solids. 60, 665-690.

Jacques, N., Mercier, S., Molinari, A., 2015. A constitutive model for porous solids taking into account microscale inertia and progressive void nucleation. Mech. Mater. 80, 311-323.

Johnson, J.N., 1981. Dynamic fracture and spallation in ductile solids. J. Appl. Phys. 52 (4), 2812-2825.

Julien, J., Garajeu, M., Michel, J.-C., 2011. A semi-analytical model for the behavior of saturated viscoplastic materials containing two populations of voids of different sizes. Int. J. Solids Struct. 48 (10), 1485-1498.

Kdir, Y.-K., Kanit, T., Zaïri, F., Naït-Abdelaziz, M., 2014. Computational homogenization of plastic porous media with two populations of voids. Materials Science and Engineering: A 597, 324-330.

Kim, J., Gao, X., Srivatsan, T., 2003. Modeling of crack growth in ductile solids: a three-dimensional analysis. Int. J. Solids Struct. 40 (26), 7357-7374.

Kim, J., Gao, X., Srivatsan, T.S., 2004. Modeling of void growth in ductile solids: effects of stress triaxiality and initial porosity. Eng. Fract. Mech. 71 (3), 379-400.
Koplik, J., Needleman, A., 1988. Void growth and coalescence in porous plastic solids. Int. J. Solids Struct. 24 (8), 835-853.

Kowalczyk-Gajewska, K., Maj, M., Bieniek, K., Majewski, M., Opiela, K.C., Zieliński, T. G., 2024. Cubic elasticity of porous materials produced by additive manufacturing: experimental analyses, numerical and mean-field modelling. Arch. Civ. Mech. Eng. 24 (1), 34.

Marvi-Mashhadi, M., Vaz-Romero, A., Sket, F., Rodriguez-Martinez, J., 2021. Finite element analysis to determine the role of porosity in dynamic localization and fragmentation: application to porous microstructures obtained from additively manufactured materials. Int. J. Plast. 143, 102999.

Massarwa, E., Czarnota, C., Molinari, A., 2024. Finite Element Modeling of steady plastic shockwaves in porous metals: role of size, shape, and spatial distribution of voids. Int. J. Impact Eng. 184, 104817.

Michel, J., Suquet, P., 1992. The constitutive law of nonlinear viscous and porous materials. J. Mech. Phys. Solids. 40 (4), 783-812.

Molinari, A., Mercier, S., 2001. Micromechanical modelling of porous materials under dynamic loading. J. Mech. Phys. Solids. 49 (7), 1497-1516.

Morin, L., Leblond, J.-B., Tvergaard, V., 2016. Application of a model of plastic porous materials including void shape effects to the prediction of ductile failure under shear-dominated loadings. J. Mech. Phys. Solids. 94, 148-166.

Needleman, A., Tvergaard, V., 1984. An analysis of ductile rupture in notched bars. J. Mech. Phys. Solids. 32 (6), 461-490.

Nielsen, K.L., Tvergaard, V., 2009. Effect of a shear modified gurson model on damage development in a fsw tensile specimen. Int. J. Solids Struct. 46 (3-4), 587-601.

Nieto-Fuentes, J.C., Jacques, N., Marvi-Mashadi, M., N'souglo, K.E., Rodriguez-Martinez, J.A., 2022. Modeling dynamic formability of porous ductile sheets subjected to biaxial stretching: actual porosity versus homogenized porosity. Int. J. Plast. 158, 103418.

Ponte Castañeda, P., 1991. The effective mechanical properties of nonlinear isotropic composites. J. Mech. Phys. Solids 39, 45-71.

Rousselier, G., 2001. The Rousselier model for porous metal plasticity and ductile fracture. Handbook of Materials Behavior Models 2, 436-445.

Reboul, J., Vadillo, G., 2018. Homogenized gurson-type behavior equations for strain rate sensitive materials. Acta Mech. 229 (8), 3517-3536.

Rice, J.R., Tracey, D.M., 1969. On the ductile enlargement of voids in triaxial stress fields. J. Mech. Phys. Solids. 17 (3), 201-217.

Sartori, C., Mercier, S., Jacques, N., Molinari, A., 2015. Constitutive behavior of porous ductile materials accounting for micro-inertia and void shape. Mech. Mater. 80, 324-339.

Sartori, C., Mercier, S., Jacques, N., Molinari, A., 2016. On the dynamic behavior of porous ductile solids containing spheroidal voids. Int. J. Solids Struct. 97, 150-167.

Subramani, M., Czarnota, C., Mercier, S., Molinari, A., 2020. Dynamic response of ductile materials containing cylindrical voids. Int. J. Fract. 222, 197-218.

Tvergaard, V., 1981. Influence of voids on shear band instabilities under plane strain conditions. Int. J. Fract. 17, 389-407.

Tvergaard, V., 1982. On localization in ductile materials containing spherical voids. Int. J. Fract. 18, 237-252.

Tvergaard, V., 1988. 3d-analysis of localization failure in a ductile material containing two size-scales of spherical particles. Eng. Fract. Mech. 31 (3), 421-436.

Vadillo, G., Fernández-Sáez, J., 2009. An analysis of Gurson model with parameters dependent on triaxiality based on unitary cells. Eur. J. Mech. Solids 28 (3), 417-427.

Versino, D., Bronkhorst, C., 2018. A computationally efficient ductile damage model accounting for nucleation and micro-inertia at high triaxialities. Comput. Methods Appl. Mech. Eng. 333, 395-420.

Vincent, P.-G., Monerie, Y., Suquet, P., 2008. Ductile damage of porous materials with two populations of voids. Comptes Rendus Mécanique. 336 (1-2), 245-259.

Vincent, P.-G., Monerie, Y., Suquet, P., 2009. Porous materials with two populations of voids under internal pressure: I. instantaneous constitutive relations. Int. J. Solids Struct. 46 (3-4), 480-506.

Vishnu, A.R., Marvi-Mashhadi, M., Nieto-Fuentes, J., Rodriguez-Martinez, J., 2022. New insights into the role of porous microstructure on dynamic shear localization. Int. J. Plast. 148, 103150.

Vishnu, A.R., Vadillo, G., Rodríguez-Martínez, J.A., 2023. Void growth in ductile materials with realistic porous microstructures. Int. J. Plast. 167, 103655.

Wang, Z.-P., 1997. Void-containing nonlinear materials subject to high-rate loading. J. Appl. Phys. 81 (11), 7213-7227.

Wang, Z.-P., Jiang, Q., 1997. A yield criterion for porous ductile media at high strain rate. J. Appl. Mech. 64 (3), 503-509.