Nonlinear Finite Element Modelling of Thermo-Visco-Plastic Styrene and Polyurethane Shape Memory Polymer Foams

: This paper presents nonlinear ﬁnite element (FE) models to predict time- and temperature-dependent responses of shape memory polymer (SMP) foams in the large deformation regime. For the ﬁrst time, an A SMP foam constitutive model is implemented in the ABAQUS FE package with the aid of a VUMAT subroutine to predict thermo-visco-plastic behaviors. A phenomenological constitutive model is reformulated adopting a multiplicative decomposition of the deformation gradient into thermal and mechanical parts considering visco-plastic SMP matrix and glass microsphere inclusions. The stress split scheme is considered by a Maxwell element in parallel with a hyper-elastic rubbery spring. The Eyring dashpot is used for modelling the isotropic resistance to the local molecular rearrangement such as chain rotation. A viscous ﬂow rule is adopted to prescribe shear viscosity and stress. An evolution rule is also considered for the athermal shear strengths to simulate macroscopic post-yield strain-softening behavior. In order to validate the accuracy of the model as well as the solution procedure, the numerical results are compared to experimental responses of Styrene and Polyurethane SMP foams at different temperatures and under different strain rates. The results show that the introduced FE modelling procedure is capable of capturing the major phenomena observed in experiments such as elastic and elastic-plastic behaviors, softening plateau regime, and densiﬁcation.


Introduction
Shape memory polymers (SMPs) are a new class of smart materials that have promising properties for different applications such as biomedical devices, actuators, and microelectromechanical systems. SMPs can respond to external stimuli and change their configuration to the original shape they "remember". Non-mutagenic, non-toxic, and biocompatible SMPs have found an increasingly key role in clinical applications in recent years. Their potential medical applications have also been exponentially increasing recently due to new developments in the manufacturing of Polyurethane (PU)-based SMP foams by means of cold hibernated elastic memory (CHEM) processing. The glass transition temperature (T g ) of CHEM foams can be adjusted for shape restoration/self-deployment of clinical devices when they are used in the human body. While inserted into the human body via small catheters, they can be miniaturized, deformed, and consequently, regain a larger pre-set shape in the right position [1]. In addition, soft robotics has benefited significantly from using SMP foams [2][3][4]. Applications of two-way and multiple-way SMPs in soft robotics had been reviewed by Scalet [5]. SMPs have been recently used in four-dimensional (4D) printing technology. For examples, Bodaghi et al. [6,7] presented self-bending/morphing/rolling structures fabricated by 4D printing and simulated their thermomechanical behaviors by using a newly-emerged simple computational tool. The shape memory behavior can be magnified in foams due to their intrinsic property of high compressibility. Several research works have been conducted within the last decade on the preparation and characterization of SMP foams for various applications.
The first applied research on SMP foams with the aim of developing buoyancy aid materials for marine applications is dated back to the 1960s. Since then, other characteristics of these materials such as their high specific strength/stiffness, high energy absorption, high biocompatibility, and flame retardancy have been further realized in different industrial areas including transportation, biomedicine, construction, and aerospace [8,9]. The popularity and versatility of these composite materials have thus attracted a lot of attention and has led to numerous studies on the improvement of the above-mentioned characteristics.
While extensive experimental studies have been performed on the thermomechanical properties of SMPs [10][11][12][13][14][15], the thermomechanical properties of SMP foams are yet to be explored. Although earlier efforts [16,17] in using rheological models could clarify the characteristic thermomechanical behavior of SMPs, the lack of strain storage and release mechanisms decrease the accuracy of the predictions. Other methods such as the mesoscale model [18,19] and molecular dynamic simulation [20] were later proposed for a detailed understanding of the underlying mechanisms. For instance, Hedayati et al. [21] proposed a multi-scale modelling approach to solve the problem of crack propagation in additively manufactured porous biomaterials. To macroscopically elucidate the thermomechanical behavior of SMPs, a variety of phenomenological approaches [22][23][24][25][26][27][28] have been recently presented. These approaches which rely on different hypotheses can be classified into two groups: mesoscale and macroscale. Chen and Lagoudas [27] provided a time-independent model for SMPs that considered a modulation of frozen and active phase effects. At temperatures below the glass transition temperature (T g ), there is a frozen phase, and at temperatures above T g , the active phase is predominant. The combination of these two phases could be useful at a temperature of T g . The stored strain (ε s ) was also suggested to determine the strain storage and release mechanisms. In this model, the parameters could be simply identified by macroscopic testing and the simulation could fairly represent the essential shape memory responses. Boatti et al. [29] introduced a 3D finite-strain phenomenological model for shape memory polymers. Nguyen et al. [28] examined the effects of stress relaxation on the structure of SMPs. They showed that the shape memory phenomena of SMPs mainly resulted from the remarkable change of molecular chain mobility caused by the glass transition [30]. Li and Xu [31] programmed SMPs at T g and stated that this programming should be done in such a way that the pre-strain is less than the strain yielding. The latest constitutive models of SMP foam are visco-plasticity models [28,32,33] which are the most accurate and complex models used in polymer modelling. Some of these models are impressively effective in predicting nonlinear, timeand temperature-dependent responses of various polymers. Recently, Walter et al. [34] presented a simple model according to the neo-Hookean model in commercial finite element (FE) software package for shape memory polyester urethane urea (PEUU) foam. Among the weak points of these models are their need for developing additional user-defined codes in the commercial software packages, a large number of experimental parameters required for the proper prediction of the material behavior, as well as being computationally costly. Due to the recent rapid evolution of the visco-plastic constitutive modelling field, new material models have been constantly presented by different researchers [35].
The literature review indicates the lack of FE modelling of SMP foams that is crucial for the design and analysis of advanced smart structures. The present study aims at developing nonlinear FE models to simulate thermo-visco-plastic behaviors of SMP foams in the large deformation regime. With the aid of this subroutine, the behavior of SMP foams under various conditions such as different temperatures and strain rates can be predicted, and the effect of applied stress on foam under loading can be assessed. A phenomenological constitutive model is reformulated adopting the multiplicative decomposition of the deformation gradient into thermal and mechanical parts considering visco-plastic SMP matrix and damaged/undamaged glass microsphere inclusions. The time-discrete counterpart of the constitutive model is developed for integrating the evolution equation by the midpoint Runge-Kutta method. The time-discrete model is implemented into the FE package with the aid of a VUMAT subroutine written using the FORTRAN programming language to predict the SMP thermomechanical behaviors. Capabilities of the model and proposed solution are demonstrated through several numerical simulations and comparisons with experimental results on Styrene and PU SMP foams in cubic and cylindrical shapes. Compressive thermomechanical behaviors of SMP foams are investigated at different temperatures and under low and high strain rates. The model and solution procedures are expected to be useful computational tools for the design and analysis of SMP foam structures under different thermomechanical loadings with nonlinear time-and temperature-dependent features.

Constitutive Modelling
The first material model for SMP-based syntactic foams was presented by Li and John [36]. In their model, they considered the visco-elastic and visco-plastic behavioral aspects of SMP foams along with the creep effects, and they proposed a rheological model for it. By performing a series of experiments, Li and Nettles [37] improved the model and presented a one-dimensional (1D) model to predict the behavior of SMP foams under compressive stress. Several theories have been presented to predict the mechanical and physical behavior of SMPs. Although rheological models [10,16,17] were proposed earlier to predict thermomechanical behaviors, they were not accurate enough because they did not take into account the effects of creep and stress relaxation. Pascon and Coda [38] proposed a model based on visco-hyper-elastic behavior and implemented the neo-Hookean hyperelastic law for their model.
The thermo-visco-elastic constitutive model previously developed by Boyce et al. [32] was reformulated and improved by Xu and Li [30] to predict the thermo-visco-plastic response of SMP foams. This model, which is also used in this study, is an expansion of the features of the rheological thermo-visco-plastic Arruda-Boyce model illustrated in Figure 1. In this model, the general deformation gradient tensor (F) is decomposed into two constituents, namely thermal (F T ) and mechanical (F M ) deformation tensors. The mechanical deformation tensor is then divided again into elastic (F e P ) and viscous (F v P ) parts. Figure 2 schematically shows these decompositions for deformation gradient tensors and indicates how they relate to each other. The constitutive relations for the mechanical stress response and thermal deformation are further divided into equilibrium and non-equilibrium parts. In addition, the time effects on the viscous deformations and non-equilibrium viscous deformations are also considered in this model.
Glass micro-balloons are used in the structure of some novel foams to improve their mechanical properties [39]. The schematic of a foam structure with its matrix embedded with micro-balloons is shown in Figure 3. As can be seen, the distribution of micro-balloons in the foam can be nonuniform. According to Xu and Li [33], the addition of these microspheres to the SMP matrix improves the properties of the foam.  Glass micro-balloons are used in the structure of some novel foams to improve their mechanical properties [39]. The schematic of a foam structure with its matrix embedded   Glass micro-balloons are used in the structure of some novel foams to improve their mechanical properties [39]. The schematic of a foam structure with its matrix embedded with micro-balloons is shown in Figure 3. As can be seen, the distribution of micro-balloons in the foam can be nonuniform. According to Xu and Li [33], the addition of these microspheres to the SMP matrix improves the properties of the foam. Glass micro-balloons are also considered in the governing equations considered for the foam structure. As shown in Figure 2, any desired thermomechanical path can be originated from an initial undeformed and unheated configuration position ( ), eventually leading to a spatial configuration ( ) that undergoes thermal and mechanical deformations. The decomposition of the deformation gradient tensor can be represented by the multiplication of its mechanical and thermal parts: where FM and/or FT contribute to taking the material configuration from Ω0 to Ω. Since the material is isotropic, the thermal deformation gradient tensor can be expressed as follows: where JT= det(FT) is the determinant of the thermal deformation gradient which indicates The volumetric thermal deformation and I denotes the second-order identity tensor. The rule of mixtures is applied to examine the composition of the syntactic foam: where indicates the deformation of the SMP matrix, represents the deformation of the glass microsphere inclusions and refers to the volume fraction of the polymer matrix. The condition of the microspheres (damaged/undamaged) is also considered in Fi as stated in Equation (S1) in the Supplementary Materials. After some manipulation on relationships and considering thermo-visco-plastic terms (the procedure of which can be found in the Supplementary Materials), the Cauchy stress can be found in terms of initial elastic response and nonlinear hyper-elastic behavior. Governing equation for stress and deformation gradient are presented in: Glass micro-balloons are also considered in the governing equations considered for the foam structure. As shown in Figure 2, any desired thermomechanical path can be originated from an initial undeformed and unheated configuration position (Ω 0 ), eventually leading to a spatial configuration (Ω) that undergoes thermal and mechanical deformations. The decomposition of the deformation gradient tensor can be represented by the multiplication of its mechanical and thermal parts: where F M and/or F T contribute to taking the material configuration from Ω 0 to Ω. Since the material is isotropic, the thermal deformation gradient tensor can be expressed as follows: where J T = det(F T ) is the determinant of the thermal deformation gradient which indicates The volumetric thermal deformation and I denotes the second-order identity tensor. The rule of mixtures is applied to examine the composition of the syntactic foam: where F p indicates the deformation of the SMP matrix, F i represents the deformation of the glass microsphere inclusions and φ p refers to the volume fraction of the polymer matrix. The condition of the microspheres (damaged/undamaged) is also considered in F i as stated in Equation (S1) in the Supplementary Materials. After some manipulation on relationships and considering thermo-visco-plastic terms (the procedure of which can be found in the Supplementary Materials), the Cauchy stress can be found in terms of initial elastic response and nonlinear hyper-elastic behavior. Governing equation for stress and deformation gradient are presented in: where σ is the Cauchy stress which is a combination of elastic σ e p and hyperelastic σ he p components. Moreover, σ ve p and σ v p represent the visco-elastic and viscous stresses. As for deformation gradients, we have: Actuators 2021, 10, 46 where F ve p and F he p are the visco-elastic and hyperelastic gradient deformation tensors, respectively. The elastic stress can be written as: where α r and α g denote the thermal expansion coefficients of the material in respectively rubbery and glassy states, and T and T 0 are current and initial temperatures, respectively.
J e p =det F e p is the elastic volume change. Furthermore, L e p = 2G P I 4 + λ p I 2 × I 2 denotes the elasticity tensor, where G P and λ p are shear modulus and Lamé constant of foam in the glassy state, and I 4 and I 2 are the fourth-and second-order identity tensors, respectively.
By applying Arruda-Boyce eight chain model [40] to the nonlinear rubbery spring, we have: where J he p is the volume change calculated as J he is the isochoric left Cauchy-Green tensor which considers the vastly different volumetric and deviatoric behaviors displayed by most amorphous polymers. Moreover, is the derivative of B [41]. λ chain = I n1 /3 denotes the effective stretch, where I n1 = tr B represents the first invariant. Moreover, k b , µ r and λ L are respectively the bulk modulus, rubbery modulus, and locking stretch of SMP, and λ L refers to the locking stretch indicating the rigidity between two elements; The Langevin function L is defined as L(β) = cot(β) − 1 β . To consider the plastic gradient deformation, viscous stretch (D v ) should be calculated. Viscous flow rule is used for describing D v as: where .
γ v is the plastic shear strain rate, and normalized vector n is defined as: The shear strain rate states the isotropic resistance to the local molecular rearrangement such as chain rotation. Nguyen et al. [28] developed the following equation for the shear strain rate based on the Eyring model [42]: .
where τ = σ e p /2 represents the equivalent shear stress, σ e p being the deviatoric part of σ e p , C 1 and C 2 indicate Williams-Landel-Ferry (WLF) constants, Q refers to the activation parameter, η g denotes the shear viscosity at T g , T f is the fictive temperature that shows the nonlinearity of structural relaxation [43], and s is the athermal shear strength.

Materials
In this study, two types of SMP foam, namely Polyurethane (PU) SMP and Styrene SMP (Veriflex) are considered for numerical simulations. The mechanical properties of SMP foams are highly dependent on the temperature at which their properties are measured, and they usually change drastically at the glass transition temperature T g .

PU SMP Foams
Shape memory Polyurethane foams (PU SMPs) are used in a wide range of applications including biomechanics [44], aerospace [45], seating parts [46], transportation [47], and furniture [47]. The most important advantage of this foam is that, in addition to having excellent biocompatibility properties, it can be moulded into the desired intricate shape without any need to be cut at a later stage [46]. The thermomechanical characterization of this type of foam has been investigated under various conditions (strain levels, temperatures, and lateral constraints) [48]. The microstructure of PU foam consists of several spherical voids with dimensions around 100-300 µm, see Figure 4a [49]. Since these voids reduce the mechanical properties of the foam, the volume fraction of the matrix materials should be considered in constitutive models. In this study, three types of PU SMP foams, namely PU-I, PU-II, and PU-III, are investigated. Table 1 lists the geometrical and physical properties of the noted foams. Shape memory Polyurethane foams (PU SMPs) are used in a wide range of applications including biomechanics [44], aerospace [45], seating parts [46], transportation [47], and furniture [47]. The most important advantage of this foam is that, in addition to having excellent biocompatibility properties, it can be moulded into the desired intricate shape without any need to be cut at a later stage [46]. The thermomechanical characterization of this type of foam has been investigated under various conditions (strain levels, temperatures, and lateral constraints) [48]. The microstructure of PU foam consists of several spherical voids with dimensions around 100-300 µm, see Figure 4a [49]. Since these voids reduce the mechanical properties of the foam, the volume fraction of the matrix materials should be considered in constitutive models. In this study, three types of PU SMP foams, namely PU-I, PU-II, and PU-III, are investigated. Table 1 lists the geometrical and physical properties of the noted foams.   The glass transition temperature of this type of foam is not constant and depends on the volume fraction of the matrix. Nevertheless, the foam material model is capable of varying the module depending on T g . Under compressive loading, the stress-strain curve of PU foam consists of three main stages: linear elastic, plateau, and densification. PU foams are capable of storing pre-strains for a long period of time (up to two months), and after this time the strain would be completely recovered [53]. The main difference between Styrene foam and PU foam is the presence of glass microspheres dispersed in the SMP matrix of Styrene foam. Figure 4b [50] shows a microscale image of a syntactic foam obtained by scanning electron microscope (SEM) imaging technique [46]. The matrix polymer is a Styrene-based thermoset SMP resin T g = 62 • C which is commercially available by Cornerstone Research Group (CRG) Industries (OH, United States) [54] under the name of Veriflex. The glass micro-balloons (Q-CEL6014) have an 85 mm average outer diameter, an effective density of 0.14 g/cm 3 , and a wall thickness of 0.8 mm. Figure 4b shows that microspheres occupy the matrix only partially, which prevents the excessive presence of voids in the matrix and thus improves its mechanical properties as compared to PU foams. Different shape memory foam types have been modelled in the study the geometrical and material properties of which are summarized in Table 2.  50 20 As can be seen in Table 2, PU foams have identical transition temperatures (T g ) while of Veriflex is more than that in PU foams. All the foams have identical low temperatures (T l ), while high temperature (T h ) is different for different foams. A significant parameter for the investigation behaviour of foams is shear modulus G P . It is obvious that this figure for Veriflex foam is higher than that in other foams, as glass micro-balloons were used in the structure of this foam boosting shear modulus. In PU foams, shear modulus depends on the relative density of foam. The denser the foams, the higher the shear modulus would be. It is obvious that the density (ρ) for PU-II and PU-III is less than that in PU-I. Therefore, the shear modulus G P for PU-II and PU-III is less than that in PU-I.

Experimental Setup
The specimens were compressed under a constant strain rate. The tests were carried out according to UNI 6350-68 in wet conditions (deionized water) at different temperatures. The test devices were equipped with a thermoregulated chamber to heat the foams. Usually, for each test, three specimens were prepared, and the average result of these samples were reported. The details of the experimental tests can be found in [11,50,52].

FE Modelling
There exist several built-in material models in ABAQUS (Dassault Systèmes, France) FE package. ABAQUS material model library is known to possess highly stable, fast, and accurate models. However, the library still does not include material models for polymers undergoing large strain deformation and having visco-plastic characteristics [55]. Therefore, the behaviors of these materials have to be defined by additional coding. ABAQUS supports both explicit and implicit user material models by means of VUMAT and UMAT which implement constitutive models via an appropriate user subroutine [56]. An explicit algorithm (VUMAT) is always recommended to be written before an implicit algorithm (UMAT) is made, due to the fact that VUMAT does not require the calculation of the tangent stiffness matrix. The VUMAT subroutines work in this way: at the beginning of each increment, deformation gradient changes resulting from mechanical and thermal loading are calculated, and then, the total deformation gradient is quantified [57]. After calculating the total deformation, the total stress is measured by Equation (9) and the obtained value is reported in this increment by ABAQUS. The development of VUMAT implemented in this study is briefly described in the following paragraphs.
The softening factor s t , total deformation gradient F t , temperature T t , and time step dt are the parameters that vary during each increment and are used as inputs for the sequential increment. The VUMAT is employed to update the stress state F t+dt at the end of the increment. It is also needed to store the updated deformation gradient, F t+dt , strain-softening parameter s t+dt , and temperature T t+dt , which all will be applied at the beginning of the next increment [35]. The procedure of update in the parameters from the beginning to the end of each increment can be expressed as The midpoint Runge-Kutta method is an integration method employed in this VU-MAT [55]. If the midpoint scheme's accuracy is not satisfied at any point, the increment must be divided into several smaller sub-increments. Then, the accuracy of the sub-increments is continuously tested and in case of errors that are too large, the sub-incrementation process is resumed with even smaller increments. Figure 5 illustrates the steps defined for this VUMAT subroutine and shows how stress is calculated in each increment.
In order to have the elastic deformation gradient at each increment, the total deformation gradient should be found at the midpoint. To do this, an average of the deformation gradient at the beginning and the end of the increment is taken.
For validation of the subroutine code, the numerical results were compared for both foam types. A simple compression was performed on the cellular structures to study the Arruda-Boyce model applied in the VUMAT user material. An SMP foam was then created with multiple QUAD element (C3D8R) which is a cubic element with eight nodes. As the Arruda-Boyce material model is a visco-plastic strain rate-dependent material model, deformation rates were taken constantly for the analyses.
The midpoint Runge-Kutta method is an integration method employed in VUMAT [55]. If the midpoint scheme's accuracy is not satisfied at any point, the incre must be divided into several smaller sub-increments. Then, the accuracy of the sub-i ments is continuously tested and in case of errors that are too large, the sub-increm tion process is resumed with even smaller increments. Figure 5 illustrates the step fined for this VUMAT subroutine and shows how stress is calculated in each increm In order to have the elastic deformation gradient at each increment, the total d mation gradient should be found at the midpoint. To do this, an average of the d mation gradient at the beginning and the end of the increment is taken.
For validation of the subroutine code, the numerical results were compared for foam types. A simple compression was performed on the cellular structures to stud Arruda-Boyce model applied in the VUMAT user material. An SMP foam was then ated with multiple QUAD element (C3D8R) which is a cubic element with eight n

Results and Discussion
In this section, the numerical stress-strain curves are compared to their experimental counterparts for two SMP foams. Various temperatures (T) and strain rates (SR) were applied to the foam models. Figure 6 compares the numerical and experimental [33] stress-strain curves under compression for Styrene SMP foam. It can be observed that the material shows different behaviors at different temperatures. While the Styrene foam shows a monotonic linear plastic hardening at the high temperature of 100 • C, a hardeningsoftening-hardening elastic-plastic behavior is observed at low and mid temperatures of 20 • C and 60 • C. There is no residual strain in the initial elastic part where the total deformation gradient F is equal to the elastic deformation gradient due to elastic bending and stretching of cell structures and compression of the void within the closed-cell. The softening section, also known as the plateau region, occurs due to the plastic deformation, buckling, and rupture of the cell structures. The end of the plateau regime can be identified by the start of the densification regime when the stress increases sharply. As can be seen, the linear elastic region decreases drastically when the temperature increases. As structural relaxation time (T s ) shortens, the molecular mobility of the polymer increases, and according to Equation (11), shear strength (s) decreases. Figure 6 reveals that the main aspects of SMP foams such as elastic and elastic-plastic behaviors, plastic growth, and softening plateau are captured well by the simulation.
ing-hardening elastic-plastic behavior is observed at low and mid temperatures of 20 °C and 60 °C. There is no residual strain in the initial elastic part where the total deformation gradient F is equal to the elastic deformation gradient due to elastic bending and stretching of cell structures and compression of the void within the closed-cell. The softening section, also known as the plateau region, occurs due to the plastic deformation, buckling, and rupture of the cell structures. The end of the plateau regime can be identified by the start of the densification regime when the stress increases sharply. As can be seen, the linear elastic region decreases drastically when the temperature increases. As structural relaxation time ( ) shortens, the molecular mobility of the polymer increases, and according to Equation (11), shear strength ( ) decreases. Figure 6 reveals that the main aspects of SMP foams such as elastic and elastic-plastic behaviors, plastic growth, and softening plateau are captured well by the simulation.  Figure 7. As can be seen, the stress level generated in the foam is reduced by increasing the temperature. It is also observed that the stress distribution in the Styrene foam is nonuniform with its maximum magnitude at the top and bottom edges of the model. In other words, stress concentration occurs in the corners due to the fact that the foam models are constrained at the noted edges.  Figure 7. As can be seen, the stress level generated in the foam is reduced by increasing the temperature. It is also observed that the stress distribution in the Styrene foam is nonuniform with its maximum magnitude at the top and bottom edges of the model. In other words, stress concentration occurs in the corners due to the fact that the foam models are constrained at the noted edges.
Finally, the thermomechanical behavior of SMP PU foams at different temperatures and loaded with different strain rates are simulated for different specimen types. Three PU specimen types with different geometrical and material characteristics obtained from three different experimental studies in the literature have been considered. The chemical compositions of PU SMP foams reported in [50,52] are slightly different which has led to different mechanical properties and material model constants in Table 2.
First, we investigate the sole effect of temperature variations on the mechanical performance of PU SMP foams. Figure 8 compares the strain-stress diagrams obtained from computational modelling with experimental stress-strain curves [50,52] loaded at different temperatures. The 3D contour of stress distribution of the relevant models at the end of loading is also depicted in Figure 9. Figure 8 reveals that under compressive loading, all SMP foam models have a hardening-softening (elastic fully-plastic) behavior at all temperatures (35,55, and 70 • C). The foam experiences an initial hardening elastic deformation followed by a plastic softening deformation. Figure 8 shows that increasing the temperature reduces the elastic strain range, the stress generated in the material, and the structural relaxation time (τ s ) as increasing the temperature enhances the molecular mobility of the polymer. A similar conclusion was achieved for Styrene SMP foams as presented in Figures 6 and 7. The stress distribution of the same (cylindrical) models is illustrated in Figure 9. Again, the maximum stress level occurs at the top and bottom edges of the model. Finally, Figure 8 implies that the FE model has been capable of accurately predicting the elastic behavior, the elastic-plastic behavior, and the softening plateau regime. Finally, the thermomechanical behavior of SMP PU foams at different temperatures and loaded with different strain rates are simulated for different specimen types. Three PU specimen types with different geometrical and material characteristics obtained from three different experimental studies in the literature have been considered. The chemical compositions of PU SMP foams reported in [50,52] are slightly different which has led to different mechanical properties and material model constants in Table 2.
First, we investigate the sole effect of temperature variations on the mechanical performance of PU SMP foams. Figure 8 compares the strain-stress diagrams obtained from computational modelling with experimental stress-strain curves [50,52] loaded at different temperatures. The 3D contour of stress distribution of the relevant models at the end of loading is also depicted in Figure 9. Figure 8 reveals that under compressive loading, all SMP foam models have a hardening-softening (elastic fully-plastic) behavior at all temperatures (35,55, and 70 °C). The foam experiences an initial hardening elastic deformation followed by a plastic softening deformation. Figure 8 shows that increasing the temperature reduces the elastic strain range, the stress generated in the material, and the structural relaxation time ( ) as increasing the temperature enhances the molecular mobility of the polymer. A similar conclusion was achieved for Styrene SMP foams as presented in Figures 6 and 7. The stress distribution of the same (cylindrical) models is illustrated in Figure 9. Again, the maximum stress level occurs at the top and bottom edges of the model. Finally, Figure 8 implies that the FE model has been capable of accurately predicting the elastic behavior, the elastic-plastic behavior, and the softening plateau regime. Next, the results of the model constructed based on [50] are presented in Figures 10  and 11. As can be seen, the stiffness of the PU foam decreases drastically by increasing the Next, the results of the model constructed based on [50] are presented in Figures 10  and 11. As can be seen, the stiffness of the PU foam decreases drastically by increasing the temperature from 20 • C to 70 • C. The stress-strain behavior changes from hardeningsoftening-hardening at 20 • C to softening-hardening at 70 • C. It can be observed that even at low temperatures, the foam model hardens at large strains due to the densification. The stress distribution of the foam model at 20 • C is also different from that at 70 • C ( Figure 11). Again, the results show that nonlinear FE modelling is an efficient tool to accurately predict thermomechanical SMP behaviors. temperature from 20 °C to 70 °C. The stress-strain behavior changes from hardening-softening-hardening at 20 °C to softening-hardening at 70 °C. It can be observed that even at low temperatures, the foam model hardens at large strains due to the densification. The stress distribution of the foam model at 20 °C is also different from that at 70 °C ( Figure  11). Again, the results show that nonlinear FE modelling is an efficient tool to accurately predict thermomechanical SMP behaviors.  The effect of change in both the strain rate and temperature on the compressive behavior of PU foams is illustrated in Figures 12 and 13. The numerical and experimental [11] results have been compared in terms of strain-stress diagrams in Figure 12. As shown in Figure 12, the foam's stress level declines with a decrease in the strain rate from 250%/min to 2.5%/min. In fact, any increase in loading duration (∆t) leads to an increase in dimensionless relaxation quantity (Δ ) (see Equation (S12)), a decrease in the viscous Figure 11. Stress distribution of PU SMP foams at different temperatures at the end of loading. The stress contours are related to the stress-strain results of the models illustrated in Figure 10  The effect of change in both the strain rate and temperature on the compressive behavior of PU foams is illustrated in Figures 12 and 13. The numerical and experimental [11] results have been compared in terms of strain-stress diagrams in Figure 12. As shown in Figure 12, the foam's stress level declines with a decrease in the strain rate from 250%/min to 2.5%/min. In fact, any increase in loading duration (∆t) leads to an increase in dimensionless relaxation quantity (∆ζ) (see Equation (S12)), a decrease in the viscous flow, and thus a reduction in the stress. Comparing the results at different temperatures as presented in Figures 12 and 13 shows that the effect of the strain rate remains significant when the temperature is increased. At each temperature, the strain rate does not have a great influence on the stress distribution in the PU SMP models. It should be mentioned that a few assumptions were considered in this study for the models. One of the assumptions was that the volume fraction remains the same during the loading process. For more accurate modelling, the dependence of volume fraction on pressure and temperature should also be considered. In addition, the VUMAT subroutine It should be mentioned that a few assumptions were considered in this study for the models. One of the assumptions was that the volume fraction remains the same during the loading process. For more accurate modelling, the dependence of volume fraction on pressure and temperature should also be considered. In addition, the VUMAT subroutine presented can be developed for other filler types used in SMP foam structures [58]. One of the most applicable fillers is carbon nanotubes. In future works, the VUMAT subroutine could be adjusted for SMP foams filled by nanotubes. Furthermore, the interaction between the particles could be considered similar to what has previously been used in relevant computational Lagrange techniques such as smoothed particle hydrodynamics (SPH) [59] by taking distance, particle mass, and particle interaction mechanism into account. Therefore, by implementing these considerations in the FE model, the numerical results would have better agreement with experimental data.

Conclusions
The main goal of this study was to simulate the large-deformation thermomechanical behavior of SMP foams by implementing a thermo-visco-plastic constitutive model in the ABAQUS FE package with the aid of an appropriate VUMAT subroutine. The model was developed based on the multiplicative decomposition of the deformation gradient into thermal and mechanical parts considering visco-elastic SMP matrix and glass microsphere inclusions. The stress split scheme was also assumed by a Maxwell element in parallel with a hyper-elastic rubbery spring. A phenomenological evolution rule for the athermal shear strengths was also applied to simulate the macroscopic post-yield strain-softening behavior. A solution algorithm for the time-discrete constitutive SMP model was developed and integrated by the midpoint Runge-Kutta method. Numerical simulations of the thermomechanical behaviors of Styrene and PU SMP foams in cubic and cylindrical shapes at different temperatures subjected to compressive loading with low and high strain rates were performed, and their results were compared with experimental results. It was seen that the thermo-visco-plastic constitutive model is capable of replicating the behavior of SMP foams under simultaneous thermal and mechanical loadings. A very good qualitative and quantitative correlation was observed between numerical and experimental results which verified the high capability of the developed FE model and the solution procedure in predicting the behavior of SMP materials. The results showed that, as compared to temperature, the strain rate has much less influence on the foam's behavior. The procedures presented in this paper provide a computationally efficient tool for designing and analysis of SMP foams under thermomechanical loads. This model could be improved by applying some additional considerations such as volume fraction variations. : fictive temperature, Equation S10: fictive temperature, Equation S11: function of fictive temperture, Equation S12: dimensionless quantity for measuring the ratio of the relaxation process, Equation S13: structural relaxation time, Equation S14: global isobaric volumetric thermal deformation, Equation S15: athermal shear strengths rule. More data that support the findings of this study are available from the author Hamid Reza Jarrah upon reasonable request.

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