Accelerated Aging on the Compression Properties of a Green Polyurethane Foam: Experimental and Numerical Analysis

The aim of this work is to evaluate the changes in compression properties of a bio-based polyurethane foam after exposure to 90 °C for different periods of time, and to propose a method to extrapolate these results and use a numerical approach to predict the compression behaviour after degradation for untested conditions at different degradation times and temperatures. Bio-based polymers are an important sustainable alternative to oil-based materials. This is explained by the foaming process and the density along the material as it was possible to see in a digital image correlation analysis. After 60 days, stiffness was approximately decreased by half in both directions. The decrease in yield stress due to thermo-oxidative degradation had a minor effect in the foaming directions, changing from 352 kPa to 220 kPa after 60 days, and the transverse property was harshly impacted changing from 530 kPa to 265 kPa. The energy absorption efficiency was slightly affected by degradation. The simulation of the compression stress-strain curves were in accordance to the experimental data and made it possible to predict the changes in mechanical properties for intermediate periods of degradation time. The plateau stress for the unaged foam transverse to the foaming direction presented experimental and numerical values of 450 kPa and 470 kPa, respectively. In addition, the plateau stresses in specimens degraded for 40 days present very similar experimental and numerical results in the same direction, at 310 kPa and 300 kPa, respectively. Therefore, this paper presents important information regarding the life-span and degradation of a green PUF. It provides insights into how compression properties vary along degradation time as function of material operation temperature, according to the Arrhenius degradation equation.


Introduction
The price of the human desire to constantly improve its quality of life has been paid over the last centuries at the expense of our planet; if no actions are taken, the damages could be irreversible as the global population is now estimated to be 7.7 billion and it is expected to reach 9.5 billion by 2050 [1]. That means we must provide more food, crops and housing for all those people which would increase the damage to our land. However, by the means of science, we are now able to process more waste, to produce more food with less land and to develop bio-based and biodegradable materials, which is a promising path toward sustainability [2,3]. Plastics are responsible for around 10% of all generated waste in the world and they comprise 60-90% of marine litter. An alternative and more sustainable way to reduce the ecological damage that polymers cause is the development and application of bio-based materials. They have been the focus of many scientific studies represents the foam's modulus of elasticity. The long plateau region extends from the yield stress point, where the linear elastic region ends, up to the start of an exponential growth in stress. In this plateau region, the stress is almost constant as the strain increases, and it represents the failure of the cells being crushed. Moreover, since the size of the cells may be different throughout the foam, some cells will fail before others, which may cause a non-constant plateau (strain hardening). Furthermore, a stress hardening behavior may appear before fracture. It appears as an exponential stress growth after the cells are completely smashed and the foam's density tends towards the polymer's [23].
The plateau region is also the main factor responsible for the high energy absorption characteristic of the foams, since it is basically pure plastic deformation. The specific energy (work per unit of volume) is the area below the stress-strain curve, and it is calculated as follows: furthermore, the efficiency of the energy absorption can be described as the ratio between the specific energy and stress. Efficiency of energy absorption is more suitable to represent the ability of a material to absorb energy, since it takes in consideration the amount of stress being applied. The energy absorption efficiency can be calculated using Equation (1.2): where U is the deformation energy, η is the energy absorption efficiency, and σ(ε) is the stress at its respective strain [23].
Marvi-Mashhadi et al. [24] carried out a high fidelity simulation of foam's mechanical properties. The authors reported a divergence between numerical and experimental results due to the gap of two years between unloading/reloading tests and the monotonic ones as they have not considered the properties varying with the aging time. Lou et al. [25] studied the accelerated ageing of rubber foams regarding compression properties. However, they did not elaborate a numerical tool to simulate other results. Hence, there is a gap in the literature regarding the prediction of mechanical properties after ageing of cellular materials, especially towards bio-based and sustainable materials.
Therefore, the aim of this work is to report the temperature degradation of a bio-based polyurethane foam regarding their compression properties. Additionally, it will provide a numerical approach to simulate those properties for untested periods of degradation. Thus, the present work pushes the boundaries of the literature towards bio-based materials and their applications.

Materials
The raw materials of the PUF used in this work were bought from Kehl Company and are cataloged by the manufacturer as IC200 and KT1106-R for the isocyanate and polyol, respectively. The isocyanate is a Methylene Diphenyl Diisocyanate (MDI) with a density of 1.22 g/cm 3 and viscosity between 170 and 250 mPa·s at room temperature. The polyol is a bio-based material from a blend of vegetable oils (mainly castor oil). Its density is 1.0 g/cm 3 and its viscosity is around 3000 mPa·s.

Manufacturing Process
The density of the foam chosen for this work was 100 kg/m 3 (0.1 g/cm 3 ) because it was the smallest density in closed expansion that showed a good dimension standardization by filling every corner with the mold. Additionally, Linul et al. (2017) [26] reported that the best density regarding energy absorption is 100 kg/m 3 . A density around this value was also studied in the works of Liu et al. (2019) [27], Mazzuca et al. (2021) [28], and Iqbal et al. (2022) [29].
The raw materials were homogenized by hand with a mass proportion of 1:1 with 8 g of each material. In order to standardize the samples, a manufacturing optimal process was applied and it is described in the flowchart in Figure 1.  Specimens that did not fit the standardization criteria were discarded. A ground steel mold closed with screws was used to contain the expansion. Figure 2 shows the manufacturing process of the foam from both components to the finished cube. The expanding direction was marked and the specimens were kept away from light, moisture and heat until the necessary amount of samples were manufactured, then, they were all together exposed to accelerated degradation.

Ageing Process
The Arrhenius equation describes the rate of a chemical reaction based on the temperature and is written as shown in Equation (2.1), [30,31]: where K is the Arrhenius degradation rate constant of the desired phenomenon; A is a pre-exponential and experimental factor; E a is the activation energy for the respective phenomenon; R the universal gas constant (8.3145 J/K/mol); and T is the temperature in Kelvin [K]. Furthermore, the E a can be obtained from a modified Arrhenius equation that relates two different degradation rates with their respective temperatures as shown in Equation (2.2): ln where T is the temperature in the accelerated environment (elevated temperature) and T is the operation temperature of the material. The E a represents the minimum amount of energy a sample must posses in order to undergo a transformation. Hence, by increasing the thermal degradation activation energy, the material thermal stability also increases. [32,33]. Consequently, using the obtained activation energy it is possible to estimate the acceleration factor by modifying the Arrhenius equation as shown in Equation (2.3) [34]: where AF is the acceleration factor that represents the time ratio of the reaction under different temperatures. In the work conducted by Berardi (2019) [35], a polyurethane foam was exposed for 4.5 months at 70°C. The operation temperature and activation energy was assumed to be 20°C and 60 kJ/mol, respectively. The author's acceleration factor was around 29 which allowed the simulation of a degradation at room temperature over more than a decade [36]. Moreover, Sandia National Laboratories [37] elaborated a report where an Arrhenius activation energy was around 78 kJ/mol. Additionally, Lee et al. (2018) [38] reported an activation energy value of 76 kJ/mol for a compression spring constant. Thus, in order to simulate degradation for long periods of time at room temperature, in the present work a bio-based PUF has been kept at 90°C for 10, 20, 30, 40, 50, and 60 days, which considering an average activation energy of 70 kJ/mol will represent an acceleration factor as a function of the operation temperature as shown in Figure 3. Here, it is shown how the AF changes according to the desired operation temperature (smaller than the accelerated degradation temperature of 90°C), since if it gets closer to 90°C, the acceleration is going to be smaller up to 1 if the temperatures are the same (that means no acceleration). Thus, the degradation in this work can simulate the degradation of the PUF for several decades depending on the operation temperature. For instance, if the temperature the material will be exposed to in real life is 40°C, the acceleration factor will be approximately 62. Therefore, as the time of exposure at elevated temperature in this work was 60 days, the degradation of compression properties of those 60 days at 90°C will represent a degradation of around 10 years (60 · 62/365) at 40°C.

Methods
The ageing process was carried out in a computer-controlled convection oven by Cienlab that kept the temperature constant for the whole time. After degradation, specimens presented a little dimensional changing, losing the adequate parallelism. Thus, they were sanded back to their original shape. The mass change (Equation (3.1)) in the aging process considered the mass before the sanding whereas the mass loss was only due to degradation: where ml% is the mass loss percentage, m 0 is initial mass, m d is the mass after the degradation process [31,39].

Compression Test
The D1621-16 standard of the American Society for Testing and Materials (ASTM) [40] specifies the test method for compressive properties of rigid cellular plastics. Therefore, it is suitable for the application in this work and widely applied for polymeric foam characterizations. It has been used in the works of Liu (2019) [27], Liu (2020) [41], and Li (2021) [42]. They performed compression tests in polyurethane foam, poly(vinyl chloride) foam and polystyrene foam, respectively.
According to the standard, the specimens should have a cross section area of at least 25.8 cm 2 , and the height should not exceed the width and depth. Therefore, considering the cross section as a square, the width and depth of the specimens should be at least 50.8 mm.
In that way, the dimensions adopted in this work were 50.8 × 50.8 × 50.8 mm 3 .
Considering the foam to be a transversely isotropic cellular material as the results presented by Tita (2012) [43] and Tagarielli (2005) [44], the compression tests were carried out in two directions for each ageing time, the first one being the expansion direction (upwards) and the second being orthogonal to the expansion. Five test specimens were created for each direction for a total of ten test specimens for each ageing time.
The tests were carried out in a universal testing machine with test speed of 5 mm/min, which represents 10% of the specimen initial height per minute. The strain values were obtained with the Digital Image Correlation (DIC) technique using the software GOM Correlate. In order to do so, one side of each specimen was painted with a black background with a pattern of white dots over it, so the software can recognize the deformation by comparing the position of the white dots. The pictures were taken every 5 s with a resolution of 6000x4000 pixels. Room lighting was good enough to provide a good pattern recognition on GOM Correlate. The test setup is shown in Figure 4. Four transverse (two for each plane orthogonal to the load-normal plane) extensometers, 45 mm long, were placed on the image of the specimens for the acquisition of the true transverse strain. The longitudinal strain was calculated from the displacement (acquired from the DIC) of the machine's cross head.
The strength data were obtained every 0.2s and a 12 degree polynomial fitting using Python Scipy Optmize Curve_fit was used to extract 100 values of strength and time that perfectly fit the raw data. Then, the same procedure was performed with the Digital Image Correlation (DIC) strain data (raw values every 5 s) to obtain 100 data sources of strain/displacement and time. Finally, strength and strain/displacement data were matched as they represent the same test time.

Numerical Implementation
The prediction of the degradation of the mechanical behavior was implemented using an Abaqus UMAT user subroutine. The constitutive model used in this work was based on the transversely isotropic model proposed by Tagarielli [44] that was derived from the generalized Hill's effective stress [45]: where E 1 is the stiffness in the transversely isotropic directions and E 3 is the reference stiffness value; ν 12 and ν 13 are the two Poisson's ratio for uniaxial loading. The shear properties for the transversely isotropic model were applied as calculated by Xu (2022) [11]. The stress was calculated using the Jaumann objective stress rate in Equation (3.3) [46]: where ∇ σ is the Jaumann stress rate; G is the shear modulus of the foam; D e is the elastic strain rate matrix; λ is the Lamé's first parameter; and I is an identity matrix. The elastic strain rate matrix is equal to the total strain rate matrix (D e = D) if the yield criterion has not been satisfied. Otherwise, the strain rate decomposition followed the rule (D = D e + D p ). Then, the material stress rate can be calculated by Equation (3.4): whereσ is the stress rate and W is the spin matrix. Both spin and total strain rate matrix were calculated by the deformation gradient given by Abaqus as shown in the following equations: where F is the deformation gradient;Ḟ is the deformation gradient rate; and L is the velocity gradient. Furthermore, in order to calculate the effective stress for the yield criterion, a transversely isotropic adaptation from the von Mises criterion was adopted as shown in Equation (3.8):σ where Q is the yield transversely isotropy matrix that represents the ratio of yield strength relative to a reference direction (dir. 3 in this case) as shown in Equation (3.9): the five constants are defined by: if after the elastic trial the yield function returns f > 0, the plasticity calculation begins in order to bring f back to 0. Moreover, if the plastic strain rate is assumed to be normal to the yield surface (associated flow), the model results in a plastic strain rate as described in Equation (3.11) where the d f /dσ represents the mapping of the effective strain rate.
Furthermore, the plastic strain rate matrix (D p ) can be obtained as shown in Equation (3.12): applying the strain decomposition it is possible to calculate the Jaumann stress rate (Equation (3.3)) and then the stress rate for plasticity (Equation (3.4)). The interaction between Abaqus and the calculation performed by the UMAT is represented in Figure 5.  The degradation of the compression properties was implemented considering the stress to be a function of time exposure in the convection oven. Thus, as shown in Equation (2.3) it will represent a degradation over longer periods of time for any application temperature under 90 • C. Therefore, the stiffness matrix C shown in Equation (3.2) obtains the shape of C(t), as well as the hardening function for the plasticity implementation H(ε p , t).
The element used in this simulation was a C3D20 20-node quadratic brick with 1000 elements. Two parallel rigid plates were positioned at the top and bottom bounders of the cube with a friction coefficient of 0.001, one plate was clamped and the other was displaced 35 mm in the desired direction to simulate the same procedure carried out in the experimental tests. This load mechanism is quite similar to the one used by Linul (2017) [26]. This macro approach represents with certain precision what happens throughout the whole material. Micro behavior approaches are more indicated when aiming to represent the stress gradient along the material, similarly to Voronoi diagram seeded cells and Laguerre tessellation, as shown in the works of Hössinger-Kalteis (2022) [47] and Shiravand (2022) [20], respectively. However, the micro failure behavior of the foam is not the emphasis of this work, instead, it aims to predict the degradation of a bio-based material over time and for this reason the macro approach was considered more reasonable.

Results
After the ageing process, the color of the foam became darker as the ageing time increased. Figure 6 shows the inner part of the specimens after each accelerated degradation time. The mass loss for ageing times up to 60 days was studied. Figure 7 shows the progression in mass change. A logarithmic mass loss fitting was adopted in the shape of ml = 2ln(0.084t + 1), as t is the number of days in the convection chamber.

Compression Test Results
The compression tests were carried out in two directions, where the direction of expansion was set as the main one and named Direction 3 (Dir3). The directions transverse to the expansion were called Direction 1 (Dir1) and Direction 2 (Dir2). Moreover, Dir1 and Dir2 were considered to be equal, since there is no reason for them to be different based on the manufacturing process and many cellular material models in the literature have made the same consideration [44,48,49]. The convention adopted in this work follows the classical transverse isotropic formulation, where E 1 = E 2 = E 3 . Figure 8 shows the compression stress-strain results for the degraded foams in both directions studied in this work. The decrement in stiffness and compression strength were well pronounced as the degradation time increased. The three compression regions are explicit in Figure 8, where the elastic region reached approximately the same strain for all the ageing times in each direction. While for Dir3 the elastic region reached a strain of around 6%, for Dir1 it extended to about 7%. The plateau region for both directions ended at around 50% of strain. Figure 9 shows the beginning of the stress-strain curve in order to better represent the elastic region of the materials. The decreasing stiffness with the increasing ageing time as well as the mass loss presented by the specimens show that the damage caused by temperature is most likely to be due to sisions of the polyurethane molecule than to the formation of cross-links between molecular chains that would increase the stiffness of the material. In the work of Mohammadi et al. [50] it was shown that during ageing at high temperatures, there is both molecular sision along the macromolecules and further cross-links between the molecular chain which is responsible for PUF's embrittlement. Regarding the plastic region, the relationship between the plateau and the foam's energy absorption efficiency is shown in Figure 10 as the densification (end of plateau) decreases the efficiency of energy absorption, therefore, this bio-based PUF as plastic energy absorption material works better for strains up to 50∼60%.
The higher compression and plateau stress of the unaged foam for both directions could represent a lower energy absorption efficiency as it is inversely proportional to the stress as shown in Equation (1.2). However, there were no noteworthy differences among the samples for each ageing time. Thus, it means that the ageing process did not harm nor increase the energy absorption efficiency, where the higher strength of the unaged foam also showed higher toughness. The difference in the mechanical properties between both directions can be easily explained by their DIC strain fields shown in Figure 11. The strain field in Dir3 shows a concentration of strain in the top of the specimens, which represents the last region of the mold filled with polymeric material. Furthermore, this top region is expected to have a smaller density than the rest of the specimen. Thus, it matches the strain field shown in Figure 11a. Since the material is being treated as transversely isotropic, Dir1 should exhibit a strain field more homogenized which is exactly as can be seen in Figure 11b. The hardening plateau presented in Figure 8a can also be explained by the DIC images, since the size of the cells is expected to be different along Dir3, the cells on the top of the specimens will fail before the ones below them. Thus, while the lower cells have not failed, the stress necessary to change the configuration of the sample slightly increases with the strain, which is in accordance to the behaviour of a typical cellular material strain hardening [23].

Numerical Prediction Results
By applying the Tagarielli yield criterion as function of the time of exposure in the convection chamber, it was possible to estimate the compression behavior of the bio-based PUF under untested ageing periods of time. The simulated stress-strain curves can properly fit the experimental ones depending on the amount of fitting parameters used in the σ(t) and E(t). In order to simplify and reduce the number of parameters and present a typical Arrhenius degradation behavior of a negative radical function of time [34], two functions were adopted in the form of Equations (4.1) and (4.2): additionally, the plasticity tangent modulus (H) was obtained by the derivative of the two variable functions represented by Equation (4.3), which was the optimal shape for thē σ(ε p , t) function based on the fitting of the experimental data.
It is important to realize that for large deformation theory, Abaqus uses nonlinear geometry configuration Thus, the strain values calculated by the software are always logarithmic strains. Hence, the strain values used in the stress-strain curves for the fitting parameters must be logarithmic, unlike the most common way to present compression values (engineering strain) since compression logarithmic strains can go beyond 100%. The fitting tool was the Python Scipy Optimize Curve_fit using the non-linear least square method. The fitting curves for the yield stress and compression elasticity modulus are shown in Figure 12, where the yield stress was obtained at the end of linearity. Furthermore, the fitting parameters obtained in this work are shown in Table A1 of Appendix A.  Figure 12a,b show that the fitting functions were adequate for the experimental compression test results. Thus, with their application in the subroutine it was possible to predict the compression properties over intermediary untested periods of time for both Dir3 and Dir1. As shown in Figure 13, which compares the experimental results with the numerical ones in Figure 13a,b for Dir3 and Dir1, respectively. The difference between the experimental and numerical values are very close given to the stiffness fitting function, as the Dir3 fitted compression elastic modulus represented by the curve in Figure 12b presents higher values for the 50 and 60 days aged foam than the experimental values for those ageing periods. Thus, it was expected that the linear elastic inclination for those numerical curves would be higher than their respective experimental curves and it is exactly what is shown in Figure 13a. However, those differences were quite small and considering the simplicity of the fitting functions the numerical results were rather accurate. The same explanation can be attributed to the yield difference for the Dir3 10 days aged foam. The plasticity H modulus were able to describe most of the plasticity behavior of the foams, especially up to the maximum energy absorption region. Furthermore, Figure 13c-f present the Dir3 simulation for untested ageing periods showing that the parameters work well for any time between 0 and 60 days. Moreover, since the tangent modulus difference is smooth between the experimental curves for Dir3, the 5 days curves practically represents a mean curve between the 0 and 10 days aged foam and the same behavior can be observed for the following ones in Dir3. Figure 13g-j also show simulations for untested aging times but for Dir1. However, as seen in Figure 13g, which represents 5 days of accelerated aging, the curve does not represent an average of 0 and 10 days, instead, it grows quite as expected until the yield point, but then, in plasticity, the stress reaches and basically follows the 10 days curve. This can be explained by the difference in the tangent modulus behavior between the unaged and the 10 days aged foam. For the unaged foam there is a noteworthy decrement in stress before the beginning of the plateau while for the 10 days aged foam the decrement was almost negligible, thus, the 5 days aged simulation inherited an average behavior between the 0 and 10 days tangent modulus, presenting a significant, yet smaller than the unaged, decrement in stress before the plateau, which was enough to place it just above the 10 days aged foam stress-strain curve. This difference in H was high enough to influence just the 5 days simulation. For the others, the behavior was basically an average stress, which means that the parameters for Dir1 were all adequate and represented a plausible simulation of the accelerated ageing process.
The maximum stress before the beginning of the plateau region (σ max ) can be considered almost the same as the plateau stress (σ pl ) for all the Dir1 results except the unaged foam, where there is a significant stress drop before the plateau begins and it was well described by the unaged simulation. As expected, Dir3 numerical results showed a hardening plateau, as discussed in Section 4.1, matching the experimental behavior. Such characteristic can be simulated for a large range of ageing times as shown by Figure 13c-f. Table 1 summarizes the numerical and experimental plateau results. The numerical σ pl for all the Dir1 ageing times were very close to the respective experimental results, except for the 50 days aged foam where plateau hardening inclination (H pl ) and maximum energy absorption strain occurred (ε η max ), since the σ pl is not constant. The H pl represents the hardening plateau more characteristic of the Dir3 behavior and shows that the model properly describes this phenomenon since small differences in the numbers do not interfere at the ε η max value. This numerical simulation is important due to the Arrhenius acceleration factor shown in Figure 3 which allows the simulation of the degradation under practical applications. Thus, the years of operation exhibited in Figure 3 changes depending on the time of exposure in the accelerated ageing environment and the simulation of a specific ageing time can represent with good accuracy the degradation in the compression properties of the bio-based foam for the desired operation time and temperature.

Conclusions
This paper investigated the behavior of a bio-based PUF when exposed to a temperature of 90°C for periods of 10, 20, 30, 40, 50 and 60 days, focusing on the compression properties. Additionally, this work aims to provide a numerical tool to predict the changes in those properties for untested periods of time. In order to do so, compression tests were carried out and numerical simulations were made using an Abaqus UMAT subroutine. The exposure of the foam to elevated temperature caused a variation in the foam's color, geometry and mass. In addition, after comparing the mechanical compression behavior of the unaged foam with the ones exposed to 90°C, a decrement was detected in stiffness and yield stress. Furthermore, this decrement was noticed after each interval of time in a behavior much like a negative radical function. This function was fitted in the degradation curve and provided degradation parameters to be used in the simulations. The implementation of these properties in UMAT proved to be promising, as the model was able to represent the experimental results with a good precision. That means that the fitting functions applied in this work are important information that can represent how the compression properties change over time. In addition, an estimation for the mechanical properties of the PUF in different temperatures can be made by calculating the acceleration factor for a desired working temperature by the Arrhenius equation. Therefore, the goals of this work have been achieved as it provides experimental data and a numerical tool to simulate the thermo-oxidative degradation of the foam, which is a key piece of information for the growing range of application of bio-based materials.

Conflicts of Interest:
The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: PUF Polyurethane foam DIC Digital image correlation Dir3 Foaming (expansion) direction Dir1 Direction transverse to the expansion

Appendix A
Fitting parameters from the Python scipy.optmize.curve_fit function used in the Abaqus UMAT subroutine for the yield strength, elastic and tangent modulus.