*2.4. Converting to Smoothed Solid CAD Models*

The STL formatted files consist of only triangle surfaces and the triangle edges are sharp. When dividing STL files to finite elements directly, the triangle surfaces are divided into further small elements that result in a considerable number of nodes and elements. Therefore, the vertices of the triangle surfaces should be interpolated by mathematically smooth surfaces. This smoothing can be performed using Recap® and Fusion 360® software, both products of Autodesk, Inc. in California, CA, USA. Firstly, the STL files with the triangle meshing were converted to surface models with quad meshing (Figure 4b). The quad meshed surfaces were then interpolated and smoothed by T-spline surfaces (Figure 4c). Finally, boundary representation solid models were generated based on the T-spline surface models (Figure 4d). The resulting solid models were then capable of being analysed in commercial FEA software.

#### *2.5. Hexahedron Dominant Meshing*

Although the geometries were smoothed by the T-spline interpolation, they were still too complex for hexahedron meshing to be applied. Therefore, mixed hexahedron and tetrahedron meshing was employed. These two kinds of elements were joined by the pyramid mesh elements. The mesh divisions were performed using Ansys® Academic Research Meshing, Release 19.2 [49]. In this study, three representative geometric models were analysed. Figure 5 shows the mesh divisions of these models and Table 1 summarises the numbers of the nodes and the elements in each. Where possible, the models were meshed with hexahedron or pyramid elements.

**Figure 5.** Mesh divisions for the models.

**Table 1.** Numbers of nodes and elements for the models.


#### *2.6. Finite Element Analyses*

The deformation behaviour of three different specimen models was calculated with the commercial FEA software Ansys® Academic Research Mechanical, Release 19.2 [50]. The large deflection was taken into account to analyse the deformations up to the plateau region. As this study focused on the static mechanical properties of polyurethane foams, the static implicit method was employed and the damping or the dynamic characteristics were neglected.

#### *2.7. Strut Material Model*

A specimen without pores is needed in order to measure the stress–strain relationship of the matrix material. The diameters of the struts are less than 0.1 mm and form a complex microstructure. Foam was compressed between plates that were heated to 150 ◦C in order to obtain a parent material specimen without pores. The original thickness of the foam was 50 mm and the compressed specimen had a thickness of 0.7 mm. The measured density of the specimen was 1200 kg/m<sup>3</sup> .

Tensile testing was performed to obtain the tensile stress–strain relationship. The test equipment was a universal testing machine AGS-X 10 kN with a 500 N load cell, products of SHIMADZU CORPORATION. The specimen was cut into 50 × 5 mm<sup>2</sup> rectangular shape specimens and then a tensile test was performed under the strain rate 0.01 s −1 . The difference between the grippers was regarded as the elongation of the specimen.

Figure 6 shows the measured nominal stress–strain curve. The experimental result is approximated by the neo-Hookean (Equation (1)) and Mooney–Rivlin (Equation (2)) hyper elastic models, respectively.

$$\mathcal{W} = \mathbb{C}\_{10}(\mathbb{I}\_1 - \mathfrak{Z}) + \frac{1}{d}(J - 1)^2 \tag{1}$$

$$\mathcal{W} = \mathbb{C}\_{10}(\mathbb{I}\_1 - \mathfrak{Z}) + \mathbb{C}\_{01}(\mathbb{I}\_2 - \mathfrak{Z}) + \frac{1}{d}(f - 1)^2 \tag{2}$$

*W* is the strain energy density, ¯*I*<sup>1</sup> and ¯*I*<sup>2</sup> are the first and second deviatoric strain invariants, and *J* is the determinant of the deformation gradient. Table 2 shows the material constants *C*10, *C*01, and *d*. Because the matrix material is thought to be incompressive, *d* was calculated to let the initial Poisson's ratio *ν* equal to 0.48. The Mooney–Rivlin model was employed in this study, as it shows better agreement with the experimental result than the neo-Hookean model.

**Figure 6.** The result of the tensile test for the matrix material and its approximations by hyperelastic models.



#### *2.8. Boundary Conditions*

The foam model specimen was uniaxially compressed between two rigid shell plates, as in Figure 7. The lower plate was fixed preventing any translational or rotational displacements. Translational displacement was applied to the upper plate, whilst all other degrees of freedom were constraint. Frictionless contacts between the foam model and the rigid walls were defined while using the penalty method with a stiffness factor of 0.01. Self-contacts between the struts were not considered, as this study focuses on the buckling behaviour in the transitions to the plateau regions. Finally, remote displacements were used to constraint the specimen lateral boundaries from rigid translational and rotational movement. The average values of the displacements of the nodes on the boundaries that correspond to these directions were fixed. This would allow for deformation, but not rigid body movement.

**Figure 7.** Boundary condition for the uniaxial compression analyses of the foam models.

#### *2.9. Experimental Measurement for the Macroscopic Stress Strain Relationships*

The uniaxial compression tests for the actual foam specimens were performed to compare with the FEA results. The testing method was similar to ISO3386-1 [51]. The 25 × 25 × 10 (mm<sup>3</sup> ) sized specimens were cut from the centre parts of the moulded foams. The specimens were set into the same equipment as seen in Section 2.7 with the compression plates. The lower plate was perforated by 6 mm holes arranged in a latticed pattern with 20 mm distances, so that the air in the foam could be ventilated. Firstly, the foams were compressed to achieve 75% nominal strain with the speed 50 mm/s as the pre-compression. Afterwards, the load was taken off with the same speed and the foams were left for 60 s. After that, the foams were compressed again with the same speed and compressive strain to measure the load and the displacement.

#### **3. Results**

#### *3.1. The Deformed Shapes of the Models*

Figure 8 shows the deformed shapes of the different specimen models at the macroscopic nominal compressive strains *ε <sup>c</sup>* = 0.05, 0.25 and 0.50 respectively. The coloured contour represents the Von–Mises equivalent strains *ε eq*. The struts bend in the linear elastic region (*ε <sup>c</sup>* = 0.05), as mentioned in Section 1. After that, some struts start to buckle, which indicates a transition to the plateau region (*ε <sup>c</sup>* = 0.25). Finally, the models gradually become denser and transfer into the densification region (*ε <sup>c</sup>* = 0.50). Because self-contacts were not applied in the foam models, the struts did not touch, but instead overlapped. The results of the analyses enable the microscopic behaviour of the struts to be carefully observed.

**Figure 8.** The deformed shapes of the models with the distributions of the equivalent strain *ε eq* .

### *3.2. Macroscopic Stress-Strain Relationships*

The FEA results were compared with the experiment results to validate the accuracy of the presented analysis method. Figure 9 shows both the experimental and FEA results of the relations between the nominal compressive stress and strain. The slopes of the stress-strain curves for the FEA results start decreasing in the strain region around 0.05 as compared to the smaller strain region. It is thought to mean the transition from the linear elastic regions to the plateau regions.

The models appear to be in good agreement with experiments in the linear elastic and the plateau regions, and up to the strain of 0.30. Differences of the stresses between the experimental and FEA results at the nominal compressive strain of 0.25 were 0.1%, 16.5%, and 6.6% for the models A, B, and C, respectively. In contrast, the FEA results are stiffer than the experimental results in larger strain regions than 0.30. After reaching the strain of 0.30, the slopes of the stress strain curves start increasing again. This behaviour looks similar to the transition to the densification regions; however, self-contacts were not

enabled within the model and the stress increase occurs far too early in the strain regions. The presented method should be modified when applied for the densification region.

**Figure 9.** The experimental and FEA results in the relations between the macroscopic compressive stress and strain.

#### **4. Discussions**

The compressive response of Polyurethane foam geometries was simulated while using FE methods and then compared with experiments. Foam specimens were scanned using X-ray CT and analysed to obtain geometries for FE simulations. The simulation results were in good agreement with experiments up to 0.3 strain. The finite element model over-predicted stresses, beyond that strain. Three different specimens were scanned and modelled to ensure the repeatability of results.

Elastic buckling appears to be one of the dominant deformation mechanisms. The finite element simulation results seem to have captured the strut deformation behaviour in agreement to relevant literature [1]. Similar deformation mechanisms have been captured with virtually generated cell structures, such as Kelvin's cells [2,4,8,9,11,12,23,29] or Voronoi polyhedrons [24,27,30,31]. Previous models based on X-ray CT scanned foam structures were mostly limited to small strains (up to 5.31%) [43–45].

Hexahedron dominant meshing was used for the large deformation analyses. Struts in foam materials can be long and narrow. Euler–Bernoulli beams have been widely employed for analytical calculations [2,4] and numerical simulations [8,9,24,27,29–31]. However whilst beam models might be beneficial in reducing complexity and calculation time, they might also add stiffness to the structure and result in higher stress predictions in comparison to the experimental values. Hexahedron meshes in large deformation problems have been used for simplified geometries [11,12]. The presented smoothing method and hexahedron dominant meshing are recommended for the complex X-ray scanned geometries.

The foam struts at the lateral specimen boundary were unconstraint. Similarly to other studies [31–33,36,43], compressive loads were applied in the model by using rigid plates. Contact was defined between the foam specimen and rigid plates. The modelled specimens were smaller than those that were used for experiments. However the boundary conditions seem to have been sufficient in capturing the strut behaviour for strains up to 0.3. The effect of surrounding material at the boundary might have been effectively negligible for up to the strain of interest due to the high porosity of the foam. However, more sophisticated boundary conditions might be required for achieving better accuracy beyond 0.3 strain, or for lower porosity foams. The surrounding cell structures could affect the computed region with bending moments, forces, or contacts between the struts, particularly as the foam densifies. These effects could potentially be taken into account by considering periodic boundary conditions [11,12]. However, this type of boundary condition requires the geometry in the model to be periodic and, therefore, might be more difficult to apply in models of stochastic foam geometries.

The finite element model over-predicted stresses, beyond 0.3 strain. Figure 10 shows an example of the deformed modelled specimen at 0.3 strain. As the specimen is compressed, struts that were initially away from the boundary, might then deform and come into contact with the loading plates at the boundary of the specimen. This could cause an increase in the stress response. This could arguably also occur during experiments; however, the model size is considerably smaller than the specimen size in experiments; therefore, the effect of these interactions would be more pronounced in the finite element model simulation. A mitigating approach could be to selectively enable contacts between the loading platens and parts of the foam, i.e., only applying contacts to the nodes on the boundary of the foam rather than the whole specimen. Increasing the model domain size could also improve results. However, a larger model would also increase the computational cost. A damage model was not included in this study. The inclusion of a damage model could potentially improve the accuracy at higher stains.

Analysing the models up to the densification region using implicit FE methods, with the periodic boundary conditions or with larger domains, remains a challenge. Additionally, investigating the effect of strut length and cross-section on the buckling behaviour of struts and the effect of the cell size variation on the linearity of the stress-strain response could inform manufacturing processes for future products. These would be the topics of future work.

**Figure 10.** Deformed shape of the model A at the strain of 0.30, when a strut contacts the upper rigid wall.

#### **5. Conclusions**

Polyurethane foam specimens, which were intended for automotive seat pads, were scanned using X-ray computed tomography. The scans were converted to 3D CAD models and used to simulate uniaxial compression test using the finite element method. The methodology for the scanning and analyses was described, and the analysis results were compared with the experiments. All three numerical models sufficiently captured the material behaviour in the linear elastic and plateau region of the stress-strain curve. The conclusions for this study are summarised below:


The presented method was successfully used to analyse foam structures and it provided a tool in understanding the mechanism of compressive deformations in polyurethane foams. Commercial CAD products and open source software were used for creating a solid mesh for FE analysis from X-ray scans. The chosen approach was perhaps more efficient by comparison to alternative specialised software at a higher cost or in-house development of custom tools. The dependence of the foam macroscopic mechanical behaviour on the microstructural features can now be further investigated to inform manufacturing processes for future polyurethane foam products.

**Author Contributions:** Conceptualization, M.I.; methodology, M.I. and R.G.; validation, M.I., P.S., B.S. and N.M.; formal analysis, M.I.; investigation, M.I.; writing—original draft preparation, M.I.; writing—review and editing, M.I., R.G., P.S., B.S. and N.M.; supervision, P.S., B.S. and N.M.; project administration, N.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing is not applicable to this article.

**Acknowledgments:** The authors kindly acknowledge Bridgestone Corporation, which supplied the specimens for the research.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**

