A Review on Mechanical Models for Cellular Media: Investigation on Material Characterization and Numerical Simulation

Cellular media materials are used for automobiles, aircrafts, energy-efficient buildings, transportation, and other fields due to their light weight, designability, and good impact resistance. To devise a buffer structure reasonably and avoid resource and economic loss, it is necessary to completely comprehend the constitutive relationship of the buffer structure. This paper introduces the progress on research of the mechanical properties characterization, constitutive equations, and numerical simulation of porous structures. Currently, various methods can be used to construct cellular media mechanical models including simplified phenomenological constitutive models, homogenization algorithm models, single cell models, and multi-cell models. This paper reviews current key mechanical models for cellular media, attempting to track their evolution from their inception to their latest development. These models are categorized in terms of their mechanical modeling methods. This paper focuses on the importance of constitutive relationships and microstructure models in studying mechanical properties and optimizing structural design. The key issues concerning this topic and future directions for research are also discussed.


Introduction
A cellular medium is a porous structure with a solid material (polymer/metal/ceramic) as the matrix and a large number of bubbles. It is a very common material in nature, especially in the biological field. Cellular media provide a variety of protective measures and structures for organisms, and have high specific strength (e.g., turtle shells and wood) [1][2][3][4][5][6][7]. Inspired by nature, researchers began to have a strong interest in the study of the diversity and mechanical properties of porous media, and began to design and analyze the mechanical behavior of different porous media materials and structures. Cellular media materials have been widely used in packaging boxes, transportation, for daily necessities, industrial manufacturing, national defense, military industries, and in biomechanical research (substrate material of elastic matrix sensor) [8], due to their many advantages including their efficient heat and sound insulation, good compression, mechanical properties and ductility, and low cost [9][10][11][12][13][14].The mechanical properties of a cellular medium are determined by both matrix and cell structure. Investigating the constitutive relationship between these can help us better understand the mechanical properties of a material and have an irreplaceable effect on the optimization of its properties [15][16][17].
The internal structure deformation and evolution of cellular media is complex, and there are many problems in traditional experimental methods. Therefore, numerical simulation based on fine analysis model has become a valid research method. There are two numerical simulation methods to investigate the mechanical properties of cellular medium. One method ignores the microstructure and explores the simplified constitutive relationship. The other method builds the relationship between the mechanical properties and microstructure. Then, the geometric modeling of the cellular medium becomes critical.
In this paper, we first introduce the characterization of a cellular medium's mechanical properties. Then, we provide the latest comments of key mechanical models, including classification of these models based on modeling scale and methodology, as well as a discussion of their advantages and disadvantages. We summarize and analyze these models in detail in Table 1. Finally, we propose prospects for these mechanical models. This work mainly focuses on the research progress of cellular medium materials in recent years, provides a reference and a compilation of important documents for researchers in the same research direction, and points out the focus of the research and potential application directions. Table 1. Performance comparison of various cellular medium mechanical models.

Models
Advantages Disadvantages

Phenomenological constitutive model
• Can obtain the mathematical expression between stress and strain by fitting experimental data; • Can study the effect of void porosity and strain rate on the stress-strain curve of foam material; • Can obtain important parameters in the compressive stress-strain curve of foam: linear elasticity (elastic modulus), yield plateau (yield stress), and densification region (densification strain).
Cannot consider the effect of the microstructure (cell size, cell shape and cell distribution) on the stress-strain curve of the foam material; Cannot describe the stress-strain relationship at higher strain rates (10 5 s −1~) or under explosive loading; The densification region fits poorly at high relative density (low void porosity).

Homogenization algorithm model
• Considered as a multiphase composite of matrix and pores; • Can quickly calculate the stress, strain, and stiffness.
Cannot give the detailed strain and stress field in each phase (matrix or pores); Cells are limited to the ellipsoidal topological structure.
Single cell model • Use the deformation mechanism of a single cell under load to represent the deformation mechanism of the entire material; • Can obtain relative modulus and relative strength of the foam.
Fit poorly when the relative density is higher than 0.3; The cell shape is too simple and far from the realistic structure.

Multi-cell model
• Can observe deformation process of the microstructure, local stress and strain; • Modeling with algorithms.
Complicated modeling algorithm; Need large-scale commercial software; Affected by random factors, and the algorithm stability is not high; High cost.
A cellular medium has its own unique characteristics. Compared with thin-walled tube and laminate composite materials, cellular media materials are isotropic and have a relatively large compressive strain, stable stress plateau, and lower density. The void  Artificial cellular medium (a) square honeycomb structure, (b) hexagonal honeycomb [2,29], (c) foaming materials [30], and stress-strain curves of cellular medium under (d) quasi-static compression [30] and (e) dynamic compression [31].
A cellular medium has its own unique characteristics. Compared with thin-walled tube and laminate composite materials, cellular media materials are isotropic and have a relatively large compressive strain, stable stress plateau, and lower density. The void  [2,29], (c) foaming materials [30], and stress-strain curves of cellular medium under (d) quasi-static compression [30] and (e) dynamic compression [31]. A cellular medium has its own unique characteristics. Compared with thin-walled tube and laminate composite materials, cellular media materials are isotropic and have a relatively large compressive strain, stable stress plateau, and lower density. The void porosity can exceed 90%, and they can be made into ultralight components. It is precisely due to the fact that cellular media materials have these excellent properties that they have received increasing attention in material research.
Previous research on cellular media has mainly focused on the influence of structural unit morphology, relative density, cell size, and defects on the material modulus, yield strength, and plateau stress. Research has also been conducted on the cell deformation characteristics of cellular media under quasi-static and dynamic loading conditions. The schematic diagrams of quasi-static loading and dynamic loading are shown in Figure 3. However, research on the energy absorption characteristics of cellular media under dynamic loads, especially the response characteristics under explosive loads, has not been fully investigated. due to the fact that cellular media materials have these excellent properties that they have received increasing attention in material research. Previous research on cellular media has mainly focused on the influence of structural unit morphology, relative density, cell size, and defects on the material modulus, yield strength, and plateau stress. Research has also been conducted on the cell deformation characteristics of cellular media under quasi-static and dynamic loading conditions. The schematic diagrams of quasi-static loading and dynamic loading are shown in Figure 3. However, research on the energy absorption characteristics of cellular media under dynamic loads, especially the response characteristics under explosive loads, has not been fully investigated.

Quasi-Static Mechanical Properties
The mechanical properties of cellular media are significantly influenced by relative density, ambient temperature, and compressive strain rate. The stress-strain curve can be divided into three stages when the cellular medium is compressed, namely linear elasticity, plateau, and densification, as shown in Figure 2d. When the strain is small, the stress-strain curve shows linear elasticity, and the stress causes the cell struts to bend. Then, a plateau appears, and the stress is almost constant as the strain increases. The cell walls and struts yield, corresponding to the plastic deformation of the cellular medium. Finally, as the cell walls are squeezed together, the cellular medium is compacted, and the stress increases rapidly. The mechanical properties of cellular media are affected by their matrix properties, relative density (void porosity), and cell size, as shown in detail in Table 2.

Quasi-Static Mechanical Properties
The mechanical properties of cellular media are significantly influenced by relative density, ambient temperature, and compressive strain rate. The stress-strain curve can be divided into three stages when the cellular medium is compressed, namely linear elasticity, plateau, and densification, as shown in Figure 2d. When the strain is small, the stress-strain curve shows linear elasticity, and the stress causes the cell struts to bend. Then, a plateau appears, and the stress is almost constant as the strain increases. The cell walls and struts yield, corresponding to the plastic deformation of the cellular medium. Finally, as the cell walls are squeezed together, the cellular medium is compacted, and the stress increases rapidly. The mechanical properties of cellular media are affected by their matrix properties, relative density (void porosity), and cell size, as shown in detail in Table 2. Some researchers found that under the same porosity, the mechanical properties did not change with the cell size [35][36][37]49,50]. Some researchers hold different views. Smits [38] found that when the cell size was large, the compressive strength of PU foam decreased as the cell diameter decreased. However, when the cell diameter decreased to a certain extent, the compressive strength increased as the cell diameter decreased. The increase in compressive strength came from the increase in the polymer fraction in the cell walls [51][52][53]. Chen [39] investigated the mechanical behavior of 3D layered nanofoam, and believed that the solid surface had fewer atoms than the inside (so it had excess surface energy). They found as the cell cross-sectional diameter decreased, the influence of surface effects on modulus and strength increased. Notario [40] prepared PMMA foams with porosity of 0.5 and cell sizes of 200 nm to 11 µm. Compared with microcellular foams, the thermomechanical properties and Shore hardness of nanocellular foams were significantly improved by 25% and 15%. Zhu [30] investigated the effect of void porosity, cell size, and cell shape on the static-compressive properties of PMMA foams. They found that porosity had greatest influence on the mechanical properties [54], cell size and cell shape also affect the compressive properties.
According to the above examples from the literature, cell sizes have different effects on the mechanical properties of cellular medium. We agree with Smits [38] that while the cell size has a turning point, different relative densities and matrix materials correspond to different turning points. Therefore, the correlation between mechanical properties and cell structure requires further research.

Dynamic Mechanical Properties
As a cushioning material, cellular media are subjected to dynamic loads in many practical applications, and the loading process is usually accompanied by a high loading strain rate. The protective foam in the side door panels of automobiles has a strain rate of up to 1500 s −1 when the velocity of side collision is 90 km/h. The strain rate of the foam energy-absorbing box of the automobile reaches 200 s −1 at 40 km/h [47]. When the cellular medium is subjected to larger strain and higher strain rate under dynamic loading, the static test result and model cannot exactly evaluate the mechanical properties of the cellular medium. Therefore, for high strain rate loading scenarios, it is important to investigate their mechanical behavior under dynamic loads [34,[55][56][57].
There are different opinions on whether cellular media have strain rate effects. Some researchers found that the strain rate had little effect on cellular medium under dynamic compression [41][42][43]58]. In addition, some studies have shown that cellular media have a strong strain rate effect [31,[44][45][46][47][48]. The deformation modes of cellular media under dynamic load are significantly different from that under quasi-static load. Cellular media have obvious deformation localization and stress enhancement under dynamic load [59][60][61][62].
Ramachandra [45] used the pressure head test to find that the effect of energy absorption on closed-cell foams with a relative density of 8% at a velocity below 10 m/s basically did not change, but at a velocity above 10 m/s, its energy absorption effect increased significantly. Sun [46] proposed that the strain rate effect of cellular media was closely related to the matrix material and cell structure. Different microstructures of cellular media experience different deformation modes and stress and strain distributions during the compression. Bouix [47] investigated the stress-strain behavior of PP foams in the strain rate range of 0.01-1500 s −1 and found that low-porosity foams had strain rate hardening due to thicker cell walls, resulting in high yield stress. However, the high-porosity foams were not sensitive to strain rate hardening. Luong [31] found that the yield strength of PVC foam material increased with the increase of strain-rate, as shown in Figure 2e. Chen [48] found that when the strain-rate exceeded 113 s −1 , the compressive strength of EPS foam increased rapidly with the increase of strain-rate, and the densification strain decreased slightly with the increase of strain-rate. In the strain-rate range of 1-280 s −1 , Young's modulus had no obvious strain rate dependence. There are also different opinions on whether cellular media can truly attenuate stress waves. Researchers used shock wave generators, explosives, and light gas guns to load metal foams to study the attenuation effect of shock waves [63]. They found that foams had better attenuation characteristics for shock waves, and the shock wave pressure decayed exponentially. Cruttmann [64] used finite element simulations to find that cellular media mainly absorbed impact loads through internal cell deformation and collapse to attenuate stress waves. Branch [65,66] found that the periodically arranged cell structure [67] can adjust the elastoplasticity of a cellular medium. Cellular media materials with simple cubic structures were prone to stress concentration, which induced emission flow and attenuated impact load energy. Cellular media materials with face-centered cubic structures were not prone to stress concentration (due to layer-to-layer symmetry) and did not induce emission flow. Radford [68] used aluminum foam as a bullet to impact a dense material, and investigated the propagation of shock waves in aluminum foam. They found that the energy carried by the shock wave was severely dissipated in the deformation of the cell structure, and the shock wave stress was related to the density of the aluminum foam.
Most of the research on the dynamic mechanical properties of cellular media are concentrated in the range of low and medium strain rates. We believe that there is now a turning point in cellular media. Above this point it has a strain rate effect, and the turning point corresponding to different relative density, cell size, and matrix material is different. Regarding the high strain rate loading of cell media, especially the strain rate higher than 10 5 s −1 and the mechanical response behavior under explosive loading, it is necessary to thoroughly study the complex wave system effects inside.

The Constitutive Relationship of Cellular Medium
The constitutive relationship can be written as the evolution equation of generalized force and generalized displacement with time. For the mechanical properties of solid materials, the constitutive relationship can be expressed as a stress-strain curve. The constitutive relationship of a cellular medium is greatly affected by relative density, temperature, and strain rate. A macroscopic theoretical cellular medium model can be constructed by experimentation to characterize the foam cushioning performance under specific conditions. Therefore, knowledge of mechanical properties is a prerequisite for establishing the constitutive relationship.
Through investigating of constitutive equations, we can fully understand the mechanical properties of materials and optimize their material properties. Many researchers put forward simplified constitutive models, which can be import into the finite element model to facilitate the calculation of complex mechanical properties [55,[69][70][71][72][73][74][75][76].

Phenomenological Constitutive Model
The phenomenological constitutive model is the mathematical expression between stress and strain obtained by fitting experimental data. It can be used to describe physical phenomena such as linear elastic region, plateau region, and densification region. However, it cannot explain the microstructure deformation and failure mechanism of the foam material, such as the weak part of the material, the stress concentration, and the local deformation [77].
Rusch [78] first raised a universal model: where E f is the apparent Young's modulus, and ψ(ε) = Aε m + Bε n is a dimensionless function of the compressive strain ε. Parameters A, B, m, and n are curve-fitting constant from a typical logarithmic plot of ψ(ε) against ε. This model was not accurate in describing the densification region. Based on the uniaxial stress-strain test results, Liu and Subhash [79] proposed a sixparameter constitutive model of polyurethane foam to fit the linear elasticity, yield plateau, and densification region: where parameters A, B, α, and β are constant for a given initial density and strain rate. A is the scaling factor for yield stress; B plays a role of shifting the lower asymptote of the function. And parameters C and γ describe the densification region, varying C implies stretching or shrinking of the curve on the strain axis. Avalle [80] improved the curve knee at the elastic and plateau region connection: where A and E are the plateau stress and the elastic modulus, exponent m represents the connection between the elastic region and the plateau region. And parameters B and n describe the densification region.
Goga [77] proposed Maxwell model: where k is equal to the elastic modulus; c is equal to the plateau stress; k p is the slope of the plateau region; and k D (ε) = γ(1 − e ε ) n , which describes the densification region Cousins [81] proposed a constitutive equation with foam strain rate effect: where K is the plateau stress, n and r are constants used to simulate strain dependence and strain rate dependence respectively, typically in the range of 0 ≤ n ≤ 4 and 0 ≤ r ≤ 0.3. The constitutive equation can describe the mechanical behavior well in the strain range of 0.1-0.7, but it cannot fit the linear elasticity and densification region well. Sherwood [82] considered the influence of the strain rate and proposed a stress-strain framework of cellular medium, as shown in the following formula: A i ε n is the strain function, A i 's are determined using a tenth-order polynomial fit. M ε, . ε 0 is a quasi-static strain rate. n(ε) = a + bε is a linear function of strain, which is the slope of the relationship between the stress and strain rate [83]. a and b are material constants in the ranges of 0-0.0442 and 0.0343-0.145 [83].
Zhang [84] considered the influencing factors of temperature and strain rate to establish the uniaxial constitutive relationships of polypropylene foam, polyurethane foam and polystyrene foam according to the work of Nagy [83]. They extended the constitutive relationship to triaxial loading and successfully used it in a finite element program. The one-dimensional constitutive model is: where σ 0 (ε) represents the nominal stress response at a constant quasi-static strain rate, T r is the glass transition temperature. The meaning of n(ε) = a + bε is reflected in Equation (6).
Polymers 2021, 13, 3283 9 of 29 C 1 and C 2 are Williams-Landel-Ferry (WLF) function coefficients to be determined from temperature effect experiments. Jeong [85,86] added a modulus function (strain rate dependence) to the Avalle's model [80], and proposed a constitutive model for the stress-strain relationship of cellular medium at low strain rates: where the meanings of A, B, E, m and n are reflected in Equation (3) [80] and these parameters can be obtained at a quasi-static strain rate of 0.001 s −1 . They thought the modulus function M ε, . ε worked as a scale factor of the shape function linearly varying along the strain axis between a reference strain rate (quasi-static strain rate 0.001 s −1 ) and another strain rate (0.1 s −1 ). And they proposed two modulus functions M ε, . ε , as shown in Equations (9) [85] and (10) [86].
where a and b are obtained by linearly approximating of the experimental stress-strain data.
M ε, where a 1 , a 2 , b 1 and b 2 , are obtained by fitting the experimental stress-strain data at strain rate of 0.1 s −1 . c 1 , c 2 , c 3 and c 4 are determined by solving the linear simultaneous Equation (12).
In the dynamic load, cellular media will inevitably encounter the impact of shock waves. Therefore, different shock wave models of cellular media are proposed to fit different types of stress-strain curves [46,87], as shown in Figure 4a. Reid and Peng [88] proposed a one-dimensional shock model. They assumed that the constitutive relationship of cellular media was independent on strain rate and obeyed the rigid-perfectly plasticlocking (R-PP-L) relationship. There were only yield stress and densification strain in the model. Lopatnikov [89][90][91] considered the effect of the elastic section of a cellular medium on the mechanical properties at low velocities on the basis of R-PP-L model, and proposed an elastic-perfectly plastic-rigid (E-PP-R) model, which took into account the elastic modulus. Zheng [92] proposed a rigid-linear hardening plastic-locking (R-LHP-L) model. A shock model under high-velocity impact and a transition model under medium-speed impact were established. In addition, Harrigan [93] used the rigid-softeninghardening (RSH) model and elastic-softening-hardening (ESH) model to fit the impact load on oak tree. Pattofatto [94,95] and Zheng [96] proposed a simple rigid-power-law harding (R-PLH) model to better evaluate the stress-strengthening level of impact enhancement.
This model gave a closed or semi-closed solution of the physical quantities on the shock wave front.

Homogenization Algorithm Model
A cellular medium can be considered as a two-phase composite with matrix and pores. The homogenization algorithm model [97][98][99][100][101][102][103][104] combines a macro-scale and a micro-scale to predict the effect of microstructure on the overall properties of the material. The two scales are connected by a representative volume element (RVE). In a micromacro-approach, at each macro-point ̅ , we know the macro-strain ̅ and need to compute the macro-stress . At micro-level, the macro-point is viewed as the center of a RVE with domain and boundary ∂ , as shown in Figure 4b. The Hill-Mandell condition, expressing the equality between energies at both scales, transforms the relation between macro-strains ̅ and macro-stresses into a relation between average strains 〈 〉 and average stresses 〈 〉 over the RVE, with ̅ = 〈 〉 and = 〈 〉. The matrix (phase 0) has volume and volume fraction = ⁄ , where is the volume of the RVE. The pores (phase 1) have a total volume and a volume fraction = ⁄ = 1 − (subscript 0 refers to the matrix and 1 to the pores). The volume averages are defined as: Although the phenomenological constitutive models can fit the compressive stressstrain curves of cellular media well, they cannot explain the mechanism of the compression process.

Homogenization Algorithm Model
A cellular medium can be considered as a two-phase composite with matrix and pores. The homogenization algorithm model [97][98][99][100][101][102][103][104] combines a macro-scale and a micro-scale to predict the effect of microstructure on the overall properties of the material. The two scales are connected by a representative volume element (RVE). In a micro-macro-approach, at each macro-point x, we know the macro-strain ε and need to compute the macro-stress σ. At micro-level, the macro-point is viewed as the center of a RVE with domain ω and boundary ∂ω, as shown in Figure 4b. The Hill-Mandell condition, expressing the equality between energies at both scales, transforms the relation between macro-strains ε and macrostresses σ into a relation between average strains ε and average stresses σ over the RVE, with ε = ε and σ = σ . The matrix (phase 0) has volume V 0 and volume fraction v 0 = V 0 /V, where V is the volume of the RVE. The pores (phase 1) have a total volume V 1 and a volume fraction v 1 = V 1 /V = 1 − v 0 (subscript 0 refers to the matrix and 1 to the pores). The volume averages are defined as: where x is the macro-point, and x is the micro-point. The averages over ω (the entire RVE), ω 0 (the matrix phase) and ω 1 (the pores phase) are related by The strain averages and stress averages can be written as [97,101]: Under linear boundary conditions, the strain averages per phase are related by a strain concentration tensor B ε as follow: To express this strain concentration tensor B ε , assumptions should be made on the micro-mechanics. This model is based on the Eshelby tensor S(I, C 0 ) and assumes that the strain field inside the inclusion (pore) is uniform and related to the macro-strain [103].
where Eshelby tensor S(I, C 0 ) depends on the geometry of the pore (I) and the stiffness tensor of the matrix C 0 . The per-phase average strains are related to the macro-strain ε by: We consider two-phase composites where the pores (phase 1) have the same shape, orientation and stiffness tensor C 1 . Thus the macroscopic material behavior can be written in the form [30,97,[101][102][103] where C is the homogenized macro-stiffness tensor. The homogenization model can predict heterogeneity by considering the volume composition of composite materials [97,[100][101][102][103][104]. Compared with the finite element (FE) calculation which takes up much CPU computing time, the homogenization model can quickly calculate the stress, strain, and stiffness of composite materials. However, it cannot give the detailed strain and stress field in each phase, and the inclusions are limited to the ellipsoidal surface.
The homogenization model is only available when the matrix phase and the inclusion phase are both isothermal linear elastic materials. However, a cellular medium is not just a linear elastic material. To study the plasticity of material, it is necessary to adopt a progressive damage model. The progressive damage model of composite material belongs to the category of continuum damage mechanism. Different from the damage model that studies the evolution of a single macroscopic crack inside a composite material, the progressive damage model introduces damage variables to express the microdefects in the entire material [105][106][107][108]. These micro-defects can be considered to be continuously distributed inside the material. They will cause the actual internal stress of the material to be higher than the characteristic stress obtained by external measurement. It can be expressed that the stiffness of the material with defects is lower than the stiffness of the nondamaged material, as shown in Figure 4b. Based on the Kachanov-Rabotnov theory, the damage variable D is introduced to describe the constitutive relationship of the damaged material. The one-dimensional characterization stress σ and the real stress (macro-stress) σ can be simplified as: Based on the equivalent strain, the constitutive relationship of the three-dimensional damaged material can be obtained. The three-dimensional characterization stress σ is: where S is the stiffness with no damage, and S is the characteristic damage stiffness related to the damage variable. D is a second-order tensor, and S is a fourth-order tensor. The stiffness matrix of the damaged material S(D) is where D ij is the element of the damage variable matrix D. S ijkl is the element of the nondamaged material matrix S.The progressive damage model describes the generation and evolution of damage variable D. Damage evolution is a function of plasticity: damage initiates at the onset of plasticity. The damage variable D in damage evolution is the key issue and foundation. It directly reflects the damage mechanism of the material. In damage mechanics, the damage variable D is used to describe the degradation process of material.
where ϕ(f) is the damage evolution function, which reflects the law of the internal damage of the material. Generally, D = 0 in the non-destructive state, and D = 1 in the complete damage state. The Hashin damage model describes damage factor D as a function of the Hashin failure criterion coefficient f. The model assumes that each phase of material obeys Weibull distribution. The damage evolution process is actually the process of strain energy release. In the process of strain energy release, the material will soften, and the macroscopic manifestation is the degradation of elastic modulus and the decrease of load-bearing capacity. The functional relationship between damage factor D and failure factor f can be defined as a variety of functional forms, so that the occurrence and evolution of damage can be described flexibly. The damage evolution function ϕ(f) includes: the step function damage evolution model Equation (24), power function damage evolution model Equation (25), linear softening damage evolution model Equation (26) and exponential function damage evolution model Equation (27) [97,[102][103][104].
where f is the failure factor of the material, f min is the minimum value, and f max is the maximum value. D max is the maximum value of material damage coefficient. α and β are the material response parameters. Although the above models can describe the stress-strain relationship of cellular media well, they cannot explain the mechanism during stretching or compression. The constitutive models of high strain rate, cell size and doped particles cannot be obtained. Therefore, it is necessary to deeply investigate the constitutive models of loading conditions and cell structure of cellular medium.

Microstructure Mechanical Models of Cellular Medium
The deformation and evolution of a cellular medium's internal structure is complex, and there are many problems in traditional experimental methods. Numerical simulation based on fine analysis model has become a valid research method.

Single Cell Model
The single cell model is based on analyzing the deformation mechanism of a tiny element (a cell) under load, as shown in Figure 5. Gibson and Ashby [28] first proposed a GA model suitable for open-cell and closed-cell foam, which simplified the cellular medium into countless cells composed of cell edges and cell walls. Under pressure, the cell edges bend, and the cell walls stretch. The relationship between the relative modulus E r , and relative strength σ r , and relative density ρ r can be obtained by: where φ is the volume fraction of cell edge, which is a correction term of 0-1. Since the value of φ is unrelated to density only in an ideal state, it is difficult to determine the value for actual closed-cell materials, especially when the density is high. Therefore, there are errors in the Gibson-Ashby model. Polymers 2021, 13, x FOR PEER REVIEW 15 of 31 Figure 5. Deformation mechanism (a) elastic collapse and (b) plastic collapse [28], and single cell model (c) open-cell unit [109], (d) closed-cell unit [110], and (e) cell model under different uniaxial loadings [111].
Choonghee [109,112] investigated the mechanical strength of open-cell PMMA foam. They assumed that the open-cell foam was composed of many equal cubic frames, and the simplified structural unit is shown in Figure 5c. Choonghee obtained the relationship between relative elastic modulus and relative density: where C is a constant fitted by experimental data. Chen [110] A can be calculated as: When the relative density is high, the closed-cell constitutive model prediction is slightly lower than the measured result. They pointed out that when the thickness of cell wall is different from the thickness of cell edge, the phenomenon of "too low" prediction will appear, and the cell structure tends to be round rather than square. Chen believed that when the relative density was lower than 0.3, the relative strength of cellular medium could be predicted well. Triawan [111] proposed an improved GA model and derived elastic modulus empirical equations under different uniaxial loadings according to experimental results.  [28], and single cell model (c) open-cell unit [109], (d) closed-cell unit [110], and (e) cell model under different uniaxial loadings [111].
Choonghee [109,112] investigated the mechanical strength of open-cell PMMA foam. They assumed that the open-cell foam was composed of many equal cubic frames, and the simplified structural unit is shown in Figure 5c. Choonghee obtained the relationship between relative elastic modulus and relative density: where C 0 is a constant fitted by experimental data. Chen [110] explored the constitutive model of closed-cell foam. The basic unit is similar to the open-cell structure unit, as shown in Figure 5d. The cell wall is added to the open-cell model. The thickness of cell wall t f is equivalent to the thickness of cell edge t e both is t/2, and the length of cell edge is also l.
According to the Gibson-Ashby model, E r of closed-cell constitutive model can be calculated as: A r can be calculated as: When the relative density is high, the closed-cell constitutive model prediction is slightly lower than the measured result. They pointed out that when the thickness of cell wall is different from the thickness of cell edge, the phenomenon of "too low" prediction will appear, and the cell structure tends to be round rather than square. Chen believed that when the relative density was lower than 0.3, the relative strength of cellular medium could be predicted well. Triawan [111] proposed an improved GA model and derived elastic modulus empirical equations under different uniaxial loadings according to experimental results.
On this basis, researchers developed tetradecahedron models (Kelvin model) [113][114][115] and octahedron models [116][117][118], as shown in Figure 6. These single cell models provide an easy method to analyze the mechanical properties, but the simple shape and regular permutation of them are too idealistic and far from the realistic structure. They cannot represent the mechanical properties of the entire cellular medium. On this basis, researchers developed tetradecahedron models (Kelvin model) [113][114][115] and octahedron models [116][117][118], as shown in Figure 6. These single cell models provide an easy method to analyze the mechanical properties, but the simple shape and regular permutation of them are too idealistic and far from the realistic structure. They cannot represent the mechanical properties of the entire cellular medium.

Multi-Cell Model
Many scholars have conducted research on the finite element (FE) method. For cellular media, modeling is important. Only when a model more closely describes real conditions, can the calculation agree with experimental values.
Some researchers [119][120][121][122][123][124] established 2D multi-cell models to investigate the mechanical properties of cellular media, as shown in Figure 7. Sun [120] thought that entrapped gas hardly affects plateau stress, but it did affect the deformation modes. In the densification region, the entrapped gas significantly enhanced crushing stress, as shown in Figure 8a,b. Generally, the effect of entrapped gas on the hexagonal cells is more important than that on the circular cells, indicating that it is dependent on the cell morphology.

Multi-Cell Model
Many scholars have conducted research on the finite element (FE) method. For cellular media, modeling is important. Only when a model more closely describes real conditions, can the calculation agree with experimental values.
Some researchers [119][120][121][122][123][124] established 2D multi-cell models to investigate the mechanical properties of cellular media, as shown in Figure 7. Sun [120] thought that entrapped gas hardly affects plateau stress, but it did affect the deformation modes. In the densification region, the entrapped gas significantly enhanced crushing stress, as shown in Figure 8a,b. Generally, the effect of entrapped gas on the hexagonal cells is more important than that on the circular cells, indicating that it is dependent on the cell morphology. Zhang [119] established a circle arc (c-a) model, as shown in Figure 7c. Under quasi-static conditions, the cellular medium deformed into an X-shape, and the broken zone spread until it was completely densified. At moderate impact velocities, V-shaped deformation can be found. At high velocities, the crushing zone is positioned high on the impact side, which is an I-shaped deformation, and then spreads along the crushing direction, as shown in Figure 8c.
Mukhopadhyay [122][123][124] investigated the mechanical properties (Young's modulus, Poisson's ratio, and shear modulus) of 2D honeycomb material based on a regular hexagon-cell model and irregular hexagon-cell model. They considered the viscoelasticity and the irregularity of the model, and found that the Young's modulus and shear modulus depend on the viscoelastic parameters, while the Poisson's ratio was independent on the viscoelastic parameters and frequency. They put forward the expressions of elastic modulus under static loading and steady-state vibration conditions. Zhang [119] established a circle arc (c-a) model, as shown in Figure 7c. Under quasistatic conditions, the cellular medium deformed into an X-shape, and the broken zone spread until it was completely densified. At moderate impact velocities, V-shaped deformation can be found. At high velocities, the crushing zone is positioned high on the impact side, which is an I-shaped deformation, and then spreads along the crushing direction, as shown in Figure 8c.
Mukhopadhyay [122][123][124] investigated the mechanical properties (Young's modulus, Poisson's ratio, and shear modulus) of 2D honeycomb material based on a regular hexagoncell model and irregular hexagon-cell model. They considered the viscoelasticity and the irregularity of the model, and found that the Young's modulus and shear modulus depend on the viscoelastic parameters, while the Poisson's ratio was independent on the viscoelastic parameters and frequency. They put forward the expressions of elastic modulus under static loading and steady-state vibration conditions. Polymers 2021, 13, x FOR PEER REVIEW 18 of 31 Some researchers established 3D multi-cell models to investigate mechanical properties of cellular medium, as shown in Figure 9. Song [125] used the finite element software ABAQUS to simulate the dynamic mechanical behavior of closed-cell foam with the Voronoi models [61,[126][127][128][129][130] and tetradecahedron models [131][132][133]. At low impact velocities, during the initial compression process, plastic deformation first begins to appear in the middle of the tetradecahedron model. At high impact velocities, cell crushing propagates downwards from the impacted surface, which exhibits an "I" shaped deformation mode, as shown in Figure 10a. Some researchers established 3D multi-cell models to investigate mechanical properties of cellular medium, as shown in Figure 9. Song [125] used the finite element software ABAQUS to simulate the dynamic mechanical behavior of closed-cell foam with the Voronoi models [61,[126][127][128][129][130] and tetradecahedron models [131][132][133]. At low impact velocities, during the initial compression process, plastic deformation first begins to appear in the middle of the tetradecahedron model. At high impact velocities, cell crushing propagates downwards from the impacted surface, which exhibits an "I" shaped deformation mode, as shown in Figure 10a. Polymers 2021, 13, x FOR PEER REVIEW 19 of 31 Figure 9. 3D multi-cell models: (a) Voronoi model [33]; (b) tetradecahedron model [33]; (c) ellipsoidal-cell model [113]; (d) body-centered cubic (BCC) model [134]; (e) face-centered cubic (FCC) model [134]; (f) "Ad-hoc" model [134].
Li [135] and Xue [33] concluded that cellular media showed different deformation modes under different impact velocities. They observed that at low impact velocities, the deformation zones were randomly distributed, while at medium velocities there were transitional deformation modes. When the impact velocity was high enough, the broken zones were locally concentrated in the impacted area and then spread along the direction of propagation. Figure 10e compares the strain rate enhancement and inertia enhancement at different impact velocities. In the random mode, the stain rate enhancement decides the strength enhancement. In the band front mode, the inertia enhancement decides the strength enhancement. Figure 9. 3D multi-cell models: (a) Voronoi model [33]; (b) tetradecahedron model [33]; (c) ellipsoidal-cell model [113]; (d) body-centered cubic (BCC) model [134]; (e) face-centered cubic (FCC) model [134]; (f) "Ad-hoc" model [134].
Li [135] and Xue [33] concluded that cellular media showed different deformation modes under different impact velocities. They observed that at low impact velocities, the deformation zones were randomly distributed, while at medium velocities there were transitional deformation modes. When the impact velocity was high enough, the broken zones were locally concentrated in the impacted area and then spread along the direction of propagation. Figure 10e compares the strain rate enhancement and inertia enhancement at different impact velocities. In the random mode, the stain rate enhancement decides the strength enhancement. In the band front mode, the inertia enhancement decides the strength enhancement.
Dabo [134] established an ad-hoc model to generate closed-cell foam microstructures, as shown in Figure 9f. Figure 10c showed that the ad-hoc model was the closest to the real microstructure and compressive behavior of the cellular medium material. The adhoc model adopted a dynamic method that considered the effect of temperature on the material characteristics. It can simultaneously consider cell nucleation, growth through gas expansion, and the interactions between growing cells.
De Giorgi [113,136] developed a randomly distributed ellipsoidal cell model based on the Kelvin cell model, as shown in Figure 9c. Fang [137,138] presented a three-dimensional (3D) mesoscopic model of closed-cell foam material [139].
These models assume equal thickness of cell wall, which is inconsistent with reality. And they are affected by random factors, and the algorithm stability is not high. There is unsatisfactory practical application of these models, and they can only provide a simple analytical method of cellular media.
To truly describe the microstructure of cellular medium, some researchers use computed tomography (CT) or X-ray scanning technology to characterize the cellular medium, as shown in Figure 11. This technology can display the material microstructure layer by layer intuitively and nondestructively. The 3D microstructure of cellular medium can be further acquired by image reconstruction technology, which restores the topological structure to the greatest extent [60,[130][131][132][133][134][135][136][137][138][139][140][141][142][143][144][145][146][147][148][149]. It has authenticity unmatched by other methods. Dabo [134] established an ad-hoc model to generate closed-cell foam microstructures, as shown in Figure 9f. Figure 10c showed that the ad-hoc model was the closest to the real microstructure and compressive behavior of the cellular medium material. The ad-hoc model adopted a dynamic method that considered the effect of temperature on the material characteristics. It can simultaneously consider cell nucleation, growth through gas expansion, and the interactions between growing cells.  [125], (b) deformation process of Voronoi model at different impact velocities [33], (c) stress-strain curves of Ad-hoc model, BCC model, FCC model and Kelvin model [134], (d) stress-strain curves of metal foam at 90 m/s [33], and (e) crushing strength under various impact velocities [33].
To truly describe the microstructure of cellular medium, some researchers use computed tomography (CT) or X-ray scanning technology to characterize the cellular medium, as shown in Figure 11. This technology can display the material microstructure layer by layer intuitively and nondestructively. The 3D microstructure of cellular medium can be further acquired by image reconstruction technology, which restores the topological structure to the greatest extent [60,[130][131][132][133][134][135][136][137][138][139][140][141][142][143][144][145][146][147][148][149]. It has authenticity unmatched by other methods. Figure 11. Flow chart of foam geometry development based on micro-CT: a flat panel detector and micro-focus X-ray source of 100 kV were used for imaging. The foam specimen was placed on a motor controlled rotating stage and radioscopic projections were taken after each degree of rotation. And the grey projections were transferred to binary images and processed to reconstruct full 3D digital foam geometry. Tetrahedral elements were used to mesh the geometry due to complex geometrical structure [34].
Huang [141] established the finite element model of cenosphere epoxy syntactic foams in ABAQUS/Explicit software using X-ray scanning technology. As shown in Figure 12a, the stress is concentrated near the top and bottom of any single cenosphere as well as the junction between adjacent cenospheres. Microcracks tend to initiate in the cenosphere connection area of the matrix. Then, the microcracks spread and merge with other microcracks and voids left by the broken cenospheres to form macrocracks. Figure 11. Flow chart of foam geometry development based on micro-CT: a flat panel detector and micro-focus X-ray source of 100 kV were used for imaging. The foam specimen was placed on a motor controlled rotating stage and radioscopic projections were taken after each degree of rotation. And the grey projections were transferred to binary images and processed to reconstruct full 3D digital foam geometry. Tetrahedral elements were used to mesh the geometry due to complex geometrical structure [34].
Huang [141] established the finite element model of cenosphere epoxy syntactic foams in ABAQUS/Explicit software using X-ray scanning technology. As shown in Figure 12a, the stress is concentrated near the top and bottom of any single cenosphere as well as the junction between adjacent cenospheres. Microcracks tend to initiate in the cenosphere connection area of the matrix. Then, the microcracks spread and merge with other microcracks and voids left by the broken cenospheres to form macrocracks.
Wang [150] used micro-CT technology to numerically model the yield characteristics of aluminum (Al) foam [151][152][153]. Figure 12b presents deformation modes under compression and tension. During the compressive process, the plastic hinges are formed first, and then the plastic collapse arises in the weak area [143]. During the tension process, crack first occurs at the cell walls in the weak area, then spreads the load, and then completely fractured.
Kader [34] developed a FE model based on micro-CT to study the elastoplastic behavior of shock wave propagation and the cell collapse mechanisms of aluminum foams. As shown in Figure 12c, a sudden rise of free surface speed was discovered at 845 m/s, indicating the development of the shock wave. It describes elastic and plastic wave propagation characteristics of the impact load in the loading direction. The elastic wave front tagged with a red line is migrating in the direction of shock. The result shows that the cell collapses in an uneven manner and continue to collapse as the waves move.
However, the modeling algorithm based on X-ray and CT is complicated, with many intermediate processing steps, and often need large-scale commercial software. The cost of using X-ray or CT scanning technology is high, which limits its wide application. Polymers 2021, 13, x FOR PEER REVIEW 22 of 31 Figure 12. (a) Predicted von Mises stress distribution and failure process during uniaxial compression [141], (b) deformation patterns under uniaxial compressive loading and tension loading [150], (c) comparison of free surface velocitytime curve between experiment and FE simulations for 845 m/s flyer plate velocity, and the visual wave propagation is also correlated with velocity-time curve [34], (d) cell collapse mechanism with time for 845 m/s flyer plate velocity, and the red line indicates the periphery of a single pore [34], and (e) elastic wave propagation characteristic through a 8 mm thickness foam at 845 m/s [34].
Wang [150] used micro-CT technology to numerically model the yield characteristics of aluminum (Al) foam [151][152][153]. Figure 12b presents deformation modes under com-  [141], (b) deformation patterns under uniaxial compressive loading and tension loading [150], (c) comparison of free surface velocity-time curve between experiment and FE simulations for 845 m/s flyer plate velocity, and the visual wave propagation is also correlated with velocity-time curve [34], (d) cell collapse mechanism with time for 845 m/s flyer plate velocity, and the red line indicates the periphery of a single pore [34], and (e) elastic wave propagation characteristic through a 8 mm thickness foam at 845 m/s [34].

Discussion
At present, researchers have performed much work characterizing mechanical properties, establishing the constitutive relationship and microstructure model, and optimizing current cellular media designs. This paper reviewed the mechanical properties characterization of cellular media, and provided the latest comments of key mechanical models, including classifying these models based on modeling scale and methodology, tracking the historical progress of these models, and discussing their scope of application, advantages and disadvantages.

•
The effects of cell size, matrix, relative density, strain rate and testing methods on the mechanical properties of cellular media were systematically investigated through quasi-static mechanical testing and dynamic mechanical testing. The investigation on mechanical properties of cellular media under explosive loading and ultra-high strain rates involves complex stress waves, and systematic analysis methods and influence laws have not yet been obtained. • The constitutive model can quickly and accurately describe the stress-strain relationship of a cellular medium, and obtain mechanical parameters (elastic modulus, plateau stress, densification strain, etc.). However, it cannot explain the microstructure deformation mechanism during stretching or compressing. The constitutive model with cell parameters (cell size and cell shape), high strain rate (greater than 10 4 s −1 ) and explosive loading needs further study. • The microstructure model can not only obtain the stress-strain curve, but also explain the microstructure deformation mechanism of cellular media, such as cell deformation, weak parts of the material, stress concentration and local deformation. However, most of the models are based on various algorithms and describe the local microstructure (RVE) of cellular media, which is greatly affected by random factors and is not stable. It is a great challenge to establish a stable, low-cost, real and effective microstructure model of cellular media. • Molecular dynamics, phase field method and peridynamics can be used to establish nanoscopic models to investigate the structural evolution and properties of molecules and atoms in cellular media. Due to these methods involve nanoscale, simulation results and experimental characterization results are not comparable, and the scale range is small and the calculation amount is large, so they are not introduced in this paper.
The research of cellular media is growing rapidly and has been widely used in fields of packaging boxes, transportation, daily necessities, industrial manufacturing, and aerospace. There will be more applications in the future:

•
Expand the mechanical constitutive models of cellular media to suitable for various complex loads, and propose more constitutive models to understand the multifunctional performance of cellular media. • Novel design of cellular metamaterials with unprecedented properties and advanced manufacturing technology [164][165][166][167]. • Molecular dynamics [168][169][170], phase field method [171][172][173], peridynamics [174] combined with multi-physics coupling calculations can obtain the nanoscale structural evolution and multifunctional performance of cellular media.

Conclusions
This paper mainly focuses on the research progress of cellular medium materials in recent years, introduces the mechanical properties characterization of cellular media, and provides mechanical models to describe their stress-strain relationship. The phenomenological constitutive models and the homogenization algorithm models can quickly obtain the stress-strain relationship of cellular media, but they cannot investigate the mechanism of the microstructure deformation. The microstructure models can not only obtain the stress-strain curve, but can also observe the microstructure deformation of a cellular medium. However, these models are not stable and require a large amount of calculation. Molecular dynamics, the phase field method, and peridynamics are the research focus of cellular media at the molecular scale and atomic scale. This paper provides a reference and a compilation of important documents for researchers pursuing the same research direction, and points out the focus of the research and potential application directions.