Elastic–plastic and viscoplastic behaviors of pore-containing solids: a computational study on periodic structures

Porous materials are solids containing voids or pores ranging from nanometers to millimeters. The internal pores can either be engineered (pre-designed), caused inherently by the fabrication or natural forming processes or gradually developed during service (by gas diffusion, irradiation etc.). Porous materials are widely used across various engineering fields due to their low density, high surface area and tunable properties. For instance, porous materials serve in lightweight structural components because they offer a unique combination of low density, high specific strength/stiffness and energy absorption (Gibson and Ashby, 1997; Yeo et al., 2019). In biomedical applications, porous metals and ceramics enable implants and scaffolds that support tissue ingrowth (He et al., 2024; Čapek et al., 2016). Properties of open-cell porous metals/alloys can be tailored to match the real cancellous bones for orthopedic applications (Lewis, 2013). In addition, the high surface area and connectivity of porous materials make them ideal for catalysts, battery and fuel cell electrodes (Wittstock et al., 2010; Chen and Ruck, 2023).

Formation of porous materials can be categorized as follows: melt foaming, where gas bubbles are introduced into molten metals (Banhart, 2001; Ashby et al., 2000), dealloying, in which elements are selectively removed from an alloy (Erlebacher et al., 2001) and additive manufacturing, where lattices with porous features are designed via 3D printing (Lv et al., 2021). Pores in materials are typically classified into two types: open-cell and closed-cell. Open-cell materials contain interconnected pores that have been widely studied for their mechanical behaviors in ceramics and metals (Ishizaki et al., 2013; Wagh et al., 1991; Gibson and Ashby, 1997). In contrast, closed-cell porous materials consist of isolated pores enclosed within the solid matrix, resulting in lower permeability and generally higher stiffness and strength (Gibson and Ashby, 1997; Iqbal and Kamiński, 2024), compared with open-cell structures. Their nonlinear and hyperelastic properties have also been examined (Yang et al., 2022).

Mechanical properties of porous materials directly influence their performance under various loading conditions. One key factor that affects their elastic behavior is porosity (pore concentration or volume fraction). Stiffness of porous materials decreases as porosity increases, because there is less solid in the material to resist elastic deformation (Lu et al., 1999; Choren et al., 2013; Cui et al., 2017; Wagh et al., 1991). Compared to elastic modulus, Poisson’s ratio shows a weak dependence on porosity (Kováčik, 2006). In addition to porosity, microstructural morphology can notably affect elastic properties of porous materials (Kikuchi et al., 2011; Pirogova et al., 2024; Yeap et al., 2013; Lee et al., 2011). For a constant pore concentration, a square array of aligned pores exhibits higher stiffness compared to a staggered arrangement (Shen and Fathi, 2021).

Plastic properties of porous materials directly influence their strength, formability, ductility and toughness. Insights into the plastic response of porous materials are thus highly important. As porosity increases, yield strength, yield surface size and tensile strength generally decrease (Ji et al., 2006; Cui et al., 2017; Cao et al., 2018; Zohdi et al., 2002). Also, the spatial geometric arrangement of pores can influence the plastic responses of porous materials at fixed porosity, including plastic flow stress and internal plastic strain distribution under mechanical loading (Shen and Fathi, 2021). In addition, as in many dense metallic and polymeric materials, porous materials can display rate-dependent viscoplasticity, meaning their plastic response also depends on the loading rate. As a consequence, accounting for viscoplasticity is essential for accurately predicting their strength, deformation and failure across a range of strain rates. In general, a higher loading rate leads to higher yield strength (Wang et al., 2001; Védrine et al., 2025; Siegkas et al., 2012; Zhao et al., 2021), as faster loading suppresses the underlying crystal defect mechanisms or molecular movement, which in turn leads to higher yield stress (Yang and Li, 2025). Furthermore, at higher loading rates, necking of a tensile specimen is delayed due to the stabilizing effect of inertia, hence resulting in higher tensile strength (N’souglo et al., 2018).

Additive manufacturing enables the creation of periodic porous structures with desired functionalities for multidisciplinary applications (Maskery et al., 2018; Feng et al., 2022). The inherent periodicity allows the use of a unit-cell approach, where a single repeating cell can represent the mechanical behavior of the entire structure, facilitating efficient parametric analyses of the mechanical properties of the porous materials. Even for materials with random distributions of internal voids, using a unit cell as the representative volume element allows for easy extraction of general trends and systematically pinpointing the connection between microscopic and macroscopic features (Shen et al., 2022; Shen and Rose, 2024).

In the present work, a three-dimensional (3D) finite element unit-cell approach is employed to investigate the elastic–plastic and viscoplastic behaviors of three periodic closed-cell porous structures: simple cubic (SC), body-centered cubic (BCC) and face-centered cubic (FCC). Specifically, effective mechanical properties are evaluated as functions of porosity as well as pore arrangement. A primary objective is to study the effective rate-dependent viscoplastic behavior influenced by internal pores, which is rarely analyzed in the open literature. Before exploring the viscoplastic response, however, purely elastic and elastic–plastic simulations need to be conducted first to establish a baseline within the current modeling framework. These rate-independent results are also used for model verification, through comparisons with representative analytical expressions. Findings in this work provide insights into the porous material behavior. We also aim to generating practical information that can advance the design of 3D printed porous materials with tailored mechanical properties.

To investigate the effective elastic–plastic and viscoplastic behaviors of periodic pore-containing solids under uniaxial loading, 3D nonlinear finite element analysis (FEA) models are established. Figure 1(a) shows a $10×10×10mm3$ cubic domain containing spherical voids, with the FCC-type of spatial distribution. Due to the geometric symmetry and periodicity, only one octant is needed to simulate the infinite structure, as shown in Figure 1(b). Note that all its six faces are mirror symmetry planes. The same approach is used to simulate the SC and BCC void configurations. The applied loading and boundary conditions along with the definition of coordinate axes are schematically shown in Figure 2. Tensile loading is implemented through the prescribed displacement on the top face along the y-direction. The bottom face is constrained to prohibit displacement in the y-direction. For both the top and bottom faces, tangential movement of the nodal points is allowed. As for the vertical side boundary planes, the back and left faces are restricted in their out-of-plane directions, and the front and right faces are allowed to move but only in a parallel manner (so they remain vertical and planar during deformation). Macroscopic uniaxial loading can thus be ensured through the above implementation, and all the nodal points in each of the boundary planes move in a coordinated manner to maintain continuity across the boundary with the (imaginary) adjacent unit structure. An infinite FCC void configuration is thus preserved (same for the SC and BCC models). If the computational domain contains a fully dense homogeneous solid free of any void, the deformation field will be uniform with σ_yy being the only non-zero stress component. The existence of voids in a porous model will disturb the deformation field, which in turn affects the overall (effective) stress–strain behavior.

Numerical modeling is carried out by using the finite element software ABAQUS (version 2023). The solid part of the model is defined to be isotropic elastic–perfectly plastic, with Young’s modulus (⁠$E0$⁠) of 196 GPa, Poisson’s ratio (⁠$ν0$⁠) of 0.26 and yield strength (⁠$σ0$⁠) of 205 MPa. While this study does not target any specific material system, these parameters are representative of typical steel properties (Callister and Rethwisch, 2018). Plastic yielding follows the von Mises yield criterion and incremental flow theory. In the rate-dependent viscoplastic simulations, the above “static” yield strength is used as a reference state when the plastic strain rate is very low, namely 1 × 10^–6 s⁻¹ and below (Singh et al., 2011). The rate-dependent viscoplastic flow stress of the solid matrix follows

where $σe$ is the von Mises effective stress and R, the yield ratio, is a function of the equivalent plastic strain rate $ε˙p$⁠. Note that R defines the ratio of yield strength at higher strain rates to the static yield strength where R equals unity. In other words, the formulation utilizes the scaling parameter R to quantify the “strain-rate hardening” effect in viscoplasticity.

The yield ratios used in the current study are based on the Johnson–Cook constitutive law (Johnson and Cook, 1983). In its temperature independent form,

where $ε˙0$ is the reference strain rate (⁠$1×10−6s−1$ used here), and m is the strain-rate sensitivity parameter which dictates the viscoplastic response. Two different m values are considered in this study, 0.0316 and 0.1633, and these are termed the “intrinsic” strain-rate sensitivities since they represent the inherent property of the solid matrix. The first value follows that of AISI 1006 steel (Fang, 2005) after conversion to the current formulation. The second value is arbitrarily taken to be about five times greater, to examine the effect of intrinsic rate dependency on the effective plastic behavior of pore-containing solids. Note that the second value signifies a very high strain rate sensitivity among most engineering alloys at or near room temperature. The applied macroscopic strain rates used in the FEA simulations span the range from 0.0001 to 0.1 s⁻¹.

This study considers the closed-cell SC, BCC and FCC pore distributions with porosities up to 20 vol.%. Maintaining relatively low porosities ensures sufficiently thick solid segments between voids to facilitate systematic comparisons. Accommodating a wider range of porosity, including transitioning to open-cell structures, will be considered in future studies. Mesh convergence analyses were performed on the constructed models. A nominally uniform mesh with C3D10 elements (10-node, second-order tetrahedral solid elements) is employed in the modeling. The true stress and true strain measures are adopted to determine the overall stress–strain response of the porous models. The applied displacement on the top face is used to directly calculate the applied true strain. The reaction force, summed over the nodal points on the top face, is divided by the instantaneous top face area (including the solid and void portions), to obtain the overall true stress.

Figure 3 shows the simulated effective Young’s modulus as a function of porosity, for the cases of SC, BCC and FCC pore arrangements. In the figure the Young’s modulus values are normalized by that of the solid matrix so the results may be viewed as generic (independent of any particular material system). Young’s modulus decreases as porosity increases. An increase in pore concentration reduces the solid fraction of the porous structure, thereby decreasing its load-bearing capacity during elastic deformation and resulting in a lower Young’s modulus. The SC structure leads to a slightly greater modulus compared to BCC and FCC at the same pore fraction. Note that, for two-phase composite materials with stiffer particles embedded in a more compliant matrix, the case of aligned particles (like the SC structure) also leads to higher effective moduli (Shen, 2010). In the present scenario where the hard particles are replaced by voids, the same trend still holds true. The cases of BCC and FCC, both featuring a staggered pore array, generate nearly identical elastic response.

Also included in Figure 3 is an analytical curve that represents the dependency of Young’s modulus on porosity, which may be used to verify the numerical model. The relation is expressed as (Gibson and Ashby, 1997)

where E is the effective Young’s modulus, $E0$ is the Young’s modulus of the fully dense/non-porous material, $Vf$ is porosity and $C1$ and $C2$ represent two geometric constants. Based on typical ranges for metallic foams (⁠$C1∼0.1−0.3$⁠; $C2∼0.3−1.0$⁠), the coefficients are taken as $C1=C2=0.3$ in this study. As Figure 3 shows, the numerical results are in reasonable agreement with the analytical prediction. Note that Eq. (3) is not a strong function of C₁ and C₂ for their specified ranges. When other values of C₁ and C₂ are chosen, the analytical prediction is still close to the present numerical results.

Figure 4 shows the contour plots of the normal stress component σ_yy along the macroscopic loading direction, in the FCC models with (a) 5% and (b) 15% pore concentrations, when the overall applied tensile strain is 0.0005. It is evident that stress concentration occurs around the “equator” of each spherical pore. In regions close to the “north and south poles” of the voids, the σ_yy stress is approximately zero due to the traction-free condition. The general patterns in Figure 4(a) and (b) are similar, although the stress values are moderately lower in the case of 15% porosity. As porosity increases, the effective load-bearing section is reduced. The reduced overall stiffness is a direct consequence of the weaker load-bearing ability of the porous structure.

Figure 5 shows the simulated effective Poisson’s ratio as a function of porosity for the three porous structures considered. While a weaker dependence on porosity is apparent, there is still a decreasing trend of the effective Poisson’s ratio as the porosity increases. With an increasing porosity, the material’s ability to induce lateral deformation caused by axial loading is expected to decrease, due to the smaller solid fraction. Figure 5 also includes an analytical curve following this expression (Gibson and Ashby, 1997),

where $ν$ is the effective Poisson’s ratio of the porous material, and $ν0$ is the Poisson’s ratio of the solid matrix. $B$ and $C$ are two geometry-dependent constants that are typically assumed to be equal for metals, meaning Poisson’s ratio is insensitive to porosity. The analytical relation is therefore shown as a horizontal line in the figure.

The numerical results in Figure 5 illustrate the fact that the BCC and FCC arrangements lead to almost identical effective Poisson’s ratios which are higher than the SC counterpart. This tendency is consistent with the case of effective Young’s modulus, where the SC geometry displays a stiffer behavior. Figure 6(a) and (b) show the contour plots of lateral displacement in the x-direction in the SC structure with porosities of 5 and 15%, respectively, when the applied overall displacement (along the y-direction) is 0.005 mm. The magnitude of maximum lateral displacement is greater in Figure 6(a) than in (b), indicating a smaller effective Poisson’s ratio if the porosity is higher. Note that, for SC, the spherical voids are perfectly aligned in the x, y and z-directions and the distance between the edges of nearest-neighbor voids is the shortest compared with BCC and FCC under a fixed pore fraction. A high degree of mutual constraint between different regions of the solid matrix can thus be expected. Figure 7(a) and (b) show the contour plots of lateral displacement in the x-direction in the BCC and FCC structures, respectively, with 15% porosity for both, when the applied overall displacement in the y-direction is 0.005 mm. Their maximum displacements are nearly the same and are greater than that in Figure 6(b) of the same porosity (15%). With the staggered distribution of voids in BCC and FCC, the color bands with higher displacement magnitudes in Figure 7(a) and (b), for example dark blue and purple, can be seen to span a greater distance towards the left side of the images, compared to Figure 6(b). This is a manifestation that the BCC and FCC geometries allow for greater lateral movements and thus higher Poisson’s ratios.

Beyond the elastic regime, the deformation behavior is dominated by plasticity. Now, the focus shifts to the plastic response under uniaxial deformation. Different porous structures show distinct effective responses at various porosities, as demonstrated by the numerical results of true stress–true strain curves shown in Figure 8(a)–(c). For each void configuration (SC, BCC or FCC), the curve of 0% Vf obtained from the modeling is also the same as the constitutive properties for the solid matrix used as model input. The overall plastic flow stress decreases as the porosity increases. It can be seen that, with the existence of pores, perfect plasticity ensues once the entire solid matrix has yielded and the effective true stress stays constant. The constant plastic flow stress values can then be plotted as a function of porosity for the different void configurations, as shown in Figure 9.

In Figure 9 the three void arrangements display slightly different decreasing tendencies as the porosity increases, with the SC structure showing the highest plastic flow stress at a given porosity. A theoretical curve is also included in the figure, following the analytical expression (Gibson and Ashby, 1997):

where $σf$ is the effective plastic flow stress, and $C3$ and $C4$ represent geometric constants, typically ranging from 0.1 to 0.3 and from 0.3 to 0.7, respectively, for metallic foams. Herein, the coefficients are taken as $C3=C4=0.3$⁠. The numerical results are moderately lower than the analytical curve, similar to the case of elastic modulus (Figure 3). It is worth mentioning that Eq. (5) is also not a strong function of C₃ and C₄, so adjusting the constants leads to only slight shifting of the analytical curve.

Figure 10(a) and (b) show the contour plots of equivalent plastic strain inside the SC porous structure with the porosities of 5 and 20%, respectively. Here the overall applied tensile strain is at 0.05. As the pore fraction increases, the solid fraction of the material decreases, resulting in reduced resistance to plastic deformation. As observed in Figure 10(a), the porous structure with low porosity has more confined plasticity with concentrated high maximum plastic strains within a small domain adjacent to the void. As for Figure 10(b) with the higher void fraction, the deformation field exhibits a more diffuse nature, with less intense but more widespread plasticity in the solid matrix, indicating that the material is easier to deform plastically. A lower plastic flow stress can thus be expected. For the BCC and FCC structures, a similar trend can be found (not shown here). Their less aligned void configurations also promote shearing deformation in the solid matrix, thus leading to easier plastic deformation compared to SC.

The current study does not target any specific material system so a baseline perfectly plastic behavior is assumed. If the material displays post-yield strain hardening, the same trend is expected to hold but the plastic flow stress will no longer be a unique value and will depend on the applied strain. The evolution of local plastic strain and damage initiation will also be affected. These additional influences will be considered in future studies.

The previous section considers rate-independent plastic response. Attention is now turned to rate-dependent viscoplasticity. The modeling approach is the same as before, except that the applied strain rates (controlled by the speed of the pulling displacement) need to be specified. The plastic flow stresses pertaining to various applied strain rates can then be analyzed, with the objective of obtaining effective strain-rate sensitivity influenced by the porosity and pore configuration.

Figure 11 shows the simulated true stress–true strain curves with various pore concentrations subjected to loading at the applied strain rate of 0.01 s⁻¹. Figure 11(a)/(b), (c)/(d) and (e)/(f) correspond to the structures of SC, BCC and FCC, respectively. For each structure two intrinsic strain-rate sensitivities (m) of 0.0316 and 0.1633 used for the solid matrix are considered. In all cases, the true stress reaches a constant value after an applied true strain of about 0.005. The trend of decreasing flow stress caused by an increasing porosity is still apparent. For the same pore configuration, a higher strain-rate sensitivity (m = 0.1633) leads to higher flow stresses due to the greater strain-rate hardening effect. It also leads to a greater span of the flow stresses between the 0 and 20% porosities.

The results shown in Figure 11 are for a specific applied strain rate of 0.01 s⁻¹. We carried out simulations of applied strain rates from 0.0001 to 0.1 s⁻¹, and the results are compiled in Figure 12 using the effective yield ratio as a measure of the plastic flow stress. The yield ratio of the pure solid used as the model input, R, is defined in Section 2. For a given constitutive viscoplastic behavior (i.e. a given strain-rate sensitivity m), there is a set of R values describing the increasing yield strength with increasing plastic strain rate. In Figure 12, the effective yield ratios of the porous material are gathered, which is the effective plastic flow stress of the perfectly plastic porous structure normalized by the static yield strength of the solid matrix, $σ0$⁠. The effective yield ratio as a function of porosity is shown in Figure 12 for the cases of: (a) SC, $m$ = 0.0316, (b) SC, $m$ = 0.1633, (c) BCC, $m$ = 0.0316, (d) BCC, $m$ = 0.1633, (e) FCC, $m$ = 0.0316, (f) FCC, $m$ = 0.1633. It can be seen that an increase in applied strain rate leads to a stronger response. This strain-rate dependence of plastic flow stress is typical of experimental observations (Wang et al., 2001) and other numerical predictions (Li and Khraishi, 2024; Zhao et al., 2021). In actual materials, the underlying deformation mechanisms are less capable of keeping up with the increased strain rate and redistributing internal stress, therefore displaying a harder response.

All the curves in Figure 12 are linear. Note that in Figure 12(a), (c) and (e), the slope of the 0% porosity lines is calculated to be 0.0316 if the natural logarithmic strain rates are used, and it is 0.163 in Figure 12(b), (d) and (f). The m values used as model input are thus recovered for the pure solid model. With the existence of voids, the slope of the line appears to decrease with an increasing porosity. This effective strain-rate sensitivity value can be calculated for each line. Figure 13 plots the calculated effective strain-rate sensitivity as a function of porosity for the SC, BCC and FCC models, with two different input m values used for the solid matrix. The decreasing trend is clear, which is also in accord with experimental findings (Wang et al., 2001). As porosity increases, the amount of load-carrying solid material decreases, weakening the propensity of rate dependency and thereby reducing the effective strain-rate sensitivity of the porous materials. It is evident in Figure 13 that the case of higher intrinsic strain-rate sensitivity (m = 0.1633 for solid) results in a more severe decline of effective rate dependency as the porosity increases, compared to the case of m = 0.0316. While the void configuration has a moderate effect in the former case (with SC being the most rate-dependent), the difference diminishes in the latter case of m = 0.0316.

Figure 14(a) and (b) show the contour plots of equivalent plastic strain in the SC structure with 5% porosity and m of 0.0316, under the applied strain rates of 0.0001 and $0.1s−1$⁠, respectively, when the applied tensile strain is 0.05. Their strain patterns are nearly the same, indicating that varying the applied strain rate has a negligible effect on the resulting deformation field within the porous structure. Although a higher strain rate leads to an increase in plastic flow stress due to the strain-rate hardening effect, the local distribution of plastic deformation is governed by the geometric constraint defined by the porous geometry as well as the intrinsic strain-rate sensitivity. The plastic strain patterns in Figure 14(a) and (b) therefore remain essentially identical.

Figure 15(a) and (b) show the contour plots of equivalent plastic strain in the SC porous structure with 20% porosity, when the intrinsic strain-rate sensitivity of the solid matrix, m, is 0.0316 and 0.1633, respectively, under the applied tensile strain of 0.05. In both cases the applied strain rate is $0.01s−1$⁠. It is observed that the local maximum plastic strain is higher when m is lower, Figure 15(a). With a high m value in Figure 15(b), the deformation field is less concentrated due to the enhanced strain-rate hardening effect. When there is stronger strain-rate hardening, a highly strained region experiences significant strengthening where further localized deformation will be restricted. Redistribution of deformation to neighboring regions will take place, leading to a less concentrated plastic strain field. Similar observations can be made for the BCC and FCC pore configurations.

A 3D finite element unit-cell approach is employed to examine the elastic–plastic and viscoplastic behaviors of periodic closed-cell porous structures. Based on the material parameters of the solid matrix used as model input, the effective properties are simulated for the SC, BCC and FCC types of pore arrangements. Numerical results reveal that increasing porosity leads to significant reductions in Young’s modulus and plastic flow stress, and slight reduction in Poisson’s ratio. The effective properties are not very sensitive to the pore configuration, although the SC geometry consistently shows the stiffest response. For the rate-dependent viscoplastic response of the porous materials, an increase in applied strain rate leads to an enhanced strength with greater plastic flow stress. The rate dependency however is weaker than the solid counterpart, causing a decreasing trend of effective strain-rate sensitivity with porosity. This declining trend is stronger when the intrinsic strain-rate sensitivity for the solid matrix is greater. Examinations of the internal deformation fields reveal that the distribution of equivalent plastic strain is dictated by the pore geometry and intrinsic strain-rate sensitivity, essentially independent of the applied strain rate. A higher intrinsic strain-rate sensitivity leads to a less concentrated plastic strain field and greater overall strength of the porous material.