A hygro-thermo-mechanical constitutive model for shape memory polymers filled with nano-carbon powder

ABSTRACT The nano-carbon powders are often used as fillers to endow the shape memory polymers (SMPs) with electroconductivity. It has been found that the shape memory effects (SMEs) of SMPs filled with nano-carbon powder can be triggered both by temperature and by water. To reveal the driving mechanism of SMEs, a constitutive model for describing the thermally activated and moisture activated SMEs of these shape memory polymer composite (SMPCs) is developed here. Because both of the SMEs share the same driving mechanism, the variable moisture is incorporated into the framework of a thermo-mechanical modeling approach to disclose the effect of moisture on the thermoviscoelastic properties. The SMPCs are regarded as isotropic materials and the effect of carbon powder on the mechanical properties of the matrix is also considered in the paper. Because the complete recovery may not be reached even they are exposed to the stimulus environment long enough, the blocking mechanism is also considered here. This is the mainly new contribution compared to the early work. Using the method of parameter determination presented here, the effectiveness of the proposed hygro-thermo-mechanical constitutive model is confirmed by comparing the model results with the test data of uniaxial deformation from the literature. GRAPHICAL ABSTRACT


Introduction
By virtue of the integration of stimulus perception and instruction execution, smart materials have gained significant attention in many fields. The successful application of smart materials has greatly promoted the development of relevant industrial areas. As a class of smart materials, shape memory polymers (SMPs) and shape memory polymer composites (SMPCs) can be used in intelligent medical devices, self-driven actuators and 4D printing industry [1][2][3] by utilizing the shape memory effects (SMEs) triggered by the environmental stimuli [4][5][6][7][8][9]. In comparison with other shape memory materials, SMPs and SMPCs have the advantages of light weight, low cost, and ease of manufacturing. Therefore, the long-term prospects of the materials are promising.
In the application of SMPs and SMPCs, it is vital to predict their SMEs for the smart structures. Therefore, it is crucial to develop corresponding constitutive models for SMPs and SMPCs with different driving mechanisms. To date, the most common SMPs are triggered by the temperature because the thermally activated SMPs can be easily prepared and their SMEs can be well controlled. The modeling approaches for the thermally activated SMPs can be classified into the phase transition approach and thermoviscoelastic approach [10][11][12][13][14]. For the semicrystalline SMPs, the temporary shape can be locked by the crystalline phase during cooling, and the shape recovery can be realized by melting in the heating process. Based on the theory of multiple natural configurations, Moon et al. [15] modeled the transition of the crystalline phase to describe the recovery to the initial state. The major advantage of the phase transition is convenience of application. Therefore, the modeling approach is also utilized to capture the shape memory mechanism of the amorphous SMPs, though the phase transition does not naturally exist in the amorphous polymers. Using the phase transition approach, Liu et al. [16] firstly modeled the SMEs of the amorphous shape memory epoxy resin under the small deformation condition. In their work, the material domain was phenomenally divided into the glassy and rubbery phases to account for the storage and release of the deformation. Qi et al. [17], Wang et al. [18], and Zhao et al. [19] further improved the phase transition approach to make it suitable for describing the SMEs of the amorphous SMPs under large deformations. With the thermoviscoelastic approach, Diani et al. [20] firstly modeled the time and temperature dependent viscoelastic behavior of amorphous polymers. To capture the multitudes of relaxation processes, Nguyen et al. [21], Westbrook et al. [22], and Yu et al. [23] further improved the model by incorporating the multiple discrete structural and stress relaxation chains. To make it with a more apparent physical significance, the thermoviscoelastic model developed by Su et al. [24] is based on the second law of thermodynamics. Recently, a series of SMPCs have been made to enhance or enable new stimulus approaches and novel SMEs. According to reinforcements, they can be classified into fiber reinforced and particle reinforced SMPCs. Generally, the micromechanical method is often used to describe the effects of reinforcements on the thermomechanical behaviors of the SMP matrix.
Similar to the effect of temperature on the material property, Hu et al. [25] found that the stiffness and strength of polymers can also be weakened by moisture. Furthermore, Yang et al. [26] found that the glass transition temperature (θ g ) of SMPs can be notably reduced by moisture. Based on the theories of polymer solution and relaxation, Lu et al. [27] studied the solution-driven SME of thermally activated SMPs in depth. It was found that the imbibed solvent molecules have a plasticizing effect on the polymer network in the form of diffusion. It releases the interactive forces among tangled polymer chains and causes shape recovery. Xiao and Nguyen [28] found that absorption of solvent at a low concentration can increase the mobility of the molecular chains and result in the decreasing of θ g . Assuming that both the thermally and solvent activated SMPs share the same deformation mechanism, Xiao [29] simulated solvent activated SMEs in the framework of the thermoviscoelastic model. It has been demonstrated that the approach is effective. Similarly, the moisture diffusion theory is incorporated into the internal state variable modeling approach in our previous work to describe the SMEs that can be triggered both by temperature and moisture [30,31]. It was found that the model results fit well with the test data for both thermally activated and moisture activated SMEs. Besides, we also found that the modeling approach based on thermodynamics with internal state variables is also effective for SMPCs [32]. Therefore, a similar modeling approach incorporated with the moisture diffusion theory is still used here to characterize the hygro-thermomechanical shape memory behavior of the SMPs filled with nano-carbon powder.
This paper is arranged in the following order. Fundamental relations of the finite deformation constitutive model for the hygrothermally activated SMPCs are introduced in Section 2. Subsequently, the method used for parameter identification is firstly introduced in Section 3. Then the comparison between the test data and the model prediction is also presented in this section. The last section reviews our work and forecasts future work.

Constitutive relations
As demonstrated previously, structural relaxation can characterize the process of the polymer chains from the non-equilibrium state to an equilibrium configuration for the thermally activated SMPs [21]. It has been proved that the fictive temperature [33], order parameters [34], or other internal variables [35] can well describe the deviation of the current state from the equilibrium configuration. Compared with the classical theories, the constitutive models based on thermodynamics with internal state variables have the following advantages. Firstly, the thermal and mechanical loadings are accounted for in a single consistent modeling approach. Secondly, the basic thermodynamic potential depends on the temperature, the stress and, a set of internal variables [36]. In our previous work, it has been proved that the modeling approach based on thermodynamics with internal state variables can be effectively extended to describe the hygro-thermomechanical SMEs of the hygrothermally activated SMPs [30]. Therefore, a similar modeling approach is still used here to characterize the SMEs of the hygrothermally activated SMPs filled with nano-carbon powder as shown in Figure 1.

Deformation and stress relations
The deformation gradient of the material configuration can be firstly divided into the mechanical part F m and hygrothermal part F wθ as follows: Then, F m is further separated into the elastic part F e and viscous part F v in the following form: By applying the polar decomposition, F e can also be split into a stretch (V e ) and a rotation (R e ) as F e = V e R e . Here, the nano-carbon powder reinforced SMPCs are considered as isotropic materials. Therefore, F wθ can be written as where J wθ is the volume ratio. The total stress can be divided into a nonequilibrium part T e and two equilibrium parts T hi as follows: Here, T e can be given by where J e ¼ det F e ð Þ and E e ¼ ln V e . L e is the isotropic elasticity tensor which can be written as where λ and G are Lamé constants. Because a fast viscous flow rate can be triggered by heating, the deformation of the nonequilibrium part can be driven by a very weak stress at high temperatures. Therefore, hyperelastic part 1 is used to describe the mechanical response at high temperatures. Here, the eight-chain network model is utilized to characterize the stress-deformation response of equilibrium part 1 as follows: where u r1 and K b are the shear modulus and bulk modulus, respectively. N is the number of Kuhn segments. The effective stretch ratio λ chain is given by λ chain ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi trðBÞ=3 q , along  the hyperelastic part and J h is the volumetric deformation. B 0 is the deviatoric part of B which can be written as B 0 ¼ B À 1 3 trðBÞI. L is the Langevin function which is defined as L β ð Þ ¼ À 1=β þ coth β ð Þ. Traditionally, one hyperelastic process is common in the constitutive modeling of SMPs with good SMEs. Noting that not all SMPs can fully recover to their initial shapes, even they are exposed to the stimulus environment long enough. The reason should be that the microstructure of the materials during recovery may differ from the initial state. It can be considered that the recovery process after cooling would be restrained by part of deformation orientated molecular chains during the pre-deformation process. Besides, the change of interfacial adhesion in nanocomposites may also hinder the recovery process. Therefore, a second hyperelastic process is used here to denote the blocking mechanism. Because the blocking effect only exists in the recovery process, the effect of the hyperelastic part 2 can be ignored in the pre-deformation. Note that this method is phenomenological. The hyperelastic part 2 takes effect only during the recovery process for the SMPs that can not notably reach their initial shapes. At the beginning of the recovery, the pre-deformation of hyperelastic part 2 is ignored and the stress generated by this part is equal to 0. Therefore, the reversed stress generated by this part in the following process impedes the recovery driven by the hyperelastic part 1. It is assumed that the stress-deformation response of hyperelastic part 2 is similar to part 1. Therefore, the eight-chain network model is also used to simulate the mechanical behavior of hyperelastic part 2 as follows.
where u r2 is the shear modulus of the hyperelastic part 2. The definitions of other parameters are the same as those in Equation (7a).

Structural relaxation and viscous flow
Based on thermodynamics with internal state variables, Lion and Yagimli [35] successfully described the material thermomechanical behavior in the vicinity of the glass transition. It has been demonstrated that the nonequilibrium response of the hygrothermally activated SMPs during the glass transition can be well described by thermodynamics with internal state variables. Assuming that the nano-carbon powders are uniformly distributed, the SMPCs can be regarded as isotropic materials. Therefore, it is also utilized to present the structural relaxation of the hygrothermally activated SMPCs. In continuum mechanics, the Gibbs free energy g can be chosen as the thermodynamic potential to characterize the thermoviscoelastic properties of glass-forming materials. It is assumed that g depends on the stress tensor T, temperature θ, moisture w, volume fraction of the filler φ and internal variables α. Therefore, the time derivative of g(T, θ, w, φ, α) can be written as Note that φ is a constant for a certain SMPC, the above equation is rewritten as In our early work, it has been proved that the second law of thermodynamics can be met by the formulation with the Gibbs free energy based on thermodynamics with internal state variables for hygrothermally activated SMPs. Therefore, the modeling approach for structural relaxation is also used here for the hygrothermally activated SMPCs. For simplicity, the modeling approach for structural relaxation will be briefly introduced here. The detailed deduction procedure can be found in Gu et al. [30].
In a reference state of thermodynamic equilibrium, all the values of the variables can be considered as reference values. To model the nonequilibrium behavior of the SMPCs, a set of time-dependent perturbation functions ΓðtÞ, �ðtÞ, #ðtÞ, ωðtÞ, η(t) and δ(t) is introduced. Therefore, the vicinity of the reference configuration can be obtained by the sufficiently small fluctuations of the variables, as follows: and α ref are the reference values of the strain tensor E, stress tensor T, temperature θ, moisture w, specific entropy s per unit mass, and internal variable tensor α. For convenience, the symmetric tensors are written in Voigt notation, e.g. E→ E. In our early work, the evolution rules for the variation of the perturbation functions are given by where D 0 is the compliance matrix. k 0 and j 0 are the thermal and hygral expansion coefficients, respectively. e, w 0 , and w 1 are the parameters determining the coupling effects between the internal variable and other variables. τ is the structural relaxation time which can be presented by the Adam-Gibbs model, as follows: where τ 0 is the reference values, B 0 is the activation parameter. It is assumed that moisture w has the similar effect with temperature θ, Equation (19) is modified as where p is the equivalent coefficient. For the materials in the glassy state, the Eyring model can be applied to the stressactivated viscoplastic flow, with the following equation: where k B and η ref are Boltzmann constant and reference shear viscosity, respectively. Q and E a are activation parameters. s y and � τ are the yield stress and equivalent shear stress, respectively. To cover the vicinity of the glass transition under the hygrothermal environment, Equation (21) is modified accordingly as follows:

Moisture diffusion
Carter and Kibler [37] found that the moisture absorbed can be split into the free and bound phases. The free phase molecules diffuse with a diffusion coefficient D γ and become bound with a probability per unit time γ. With a probability per unit time β, molecules can be emitted from the bound phase and become free in the meantime. The following equations show the evolution rules for the number of free molecules (n) and the number of bound molecules (N) in the direction of thickness (2Δ).
where n ∞ is the equilibrium value for the free molecules per unit volume. r � can be given by ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where The moisture content at time t can be calculated as where w ∞ is the equilibrium value of the weight fraction of moisture. Therefore the moisture contents at time t for the free molecules and bound molecules can be respectively calculated as

Results and discussions
The model is implemented in a C++ program to carry out the numerical calculation. To verify the validity of the model proposed above, a comparison of the simulation results and test data found in Yang et [38]. and Yang [39] is presented in this section.

Parameter determination
Yang et [38]. and Yang [39] systematically studied the effects of temperature and moisture on the polyurethane SMPs filled with nano-carbon powder. It is found that θ g of the SMPCs can be effectively reduced by raising the carbon powder content and increasing the immersing time in the water. Zhou et al. [40] found that the fractional free volume of the system f depends on the temperature and moisture. In the theory of free volume, the volume of materials can be divided into two parts. One part is occupied by the molecules and the other unoccupied part is the free volume which is distributed over the whole substance in the form of voids. The mobility of the molecular chains is due to the existence of the free volume. Evidently, both the moisture and filler added into the polymer matrix can enlarge the free volume of the material. Assuming that the effect of filler on the free volume is the same as the moisture, the following equations can be obtained.
where w 0 and φ 0 are the moisture and filler volume fraction of the reference state, respectively. Their values can be set to 0. k 0 is the thermal expansion coefficient. θ g0 is the glass transition temperature for the pure dry SMPs. Parameter ζ denotes the contribution of the moisture to the free volume. Because the water can be divided into the free and bound phases, their contributions should be distinguished by where ζ n , ζ N , and ξ are the proportional coefficients of the free water, bound water, and filler volume fraction, respectively. In the condition of θ = θ g0 , the variation of θ g can be calculated from Equations (31) and (32), as follows: Yang et al. [26] experimentally studied the effect of the moisture on the glass transition temperature for the polyurethane SMPs. They found that θ g can be observably affected by the bound water, while the effect of free water can be ignored. Therefore, Equation (33) can be simplified as In the work of Yang et al. [38], differential scanning calorimetry tests were carried out to measure the glass transition temperature of the polyurethane SMPCs. It was found that θ g decreases with the increase of the carbon powder in volume fraction. As shown in Figure 2, the value of ξ/k 0 can be obtained by the method of ordinary least squares.
In the work of Yang [39], the absorbed water was quantitatively divided into the free and bound water by the relevant tests. Therefore, the parameters used in the rule of moisture diffusion can be determined from the corresponding tests. Specifically, w ∞ is chosen as the saturation value during immersion. Based on a nonlinear regularization program, the other parameters, i.e. β, γ, and D γ , can be determined. In the paper, the abbreviation CBX is used to denote the volume fraction of X% carbon powder. The parameter fitting for the moisture diffusion of CB0 is shown in Figure 3. Note that the model results (curves) fit well with the test data (points). Along with the evolution law of the bound water, the parameter ξ N /k 0 used in Equation (34) can be determined by fitting to the test data as shown in Figure 4.
Under the assumption that the phase transition parameters β and γ remain unchanged after the carbon powder is added to the SMP matrix, D γ can be determined from the process of water absorption. As shown in Figure 5, the diffusion coefficients D γ of CB4, CB7, CB10, CB13, and CB15 are determined by fitting to the test data. In the meantime, the free and bound water absorbed in these SMPCs is predicted. Therefore, the parameter ξ N /k 0 used for these SMPCs can be obtained by fitting to the test data as shown in Figure 6. It  is shown that the results are in good agreement with the test data, which also demonstrates the validity of the model.
To investigate the mechanical properties of the SMPCs, uniaxial tensile tests were carried out by Yang [39] at room temperature. It is found that the stiffness increases    with the increase of the carbon powder. Note that the distribution form of carbon powders affects the interfacial adhesion in nanocomposites, which may lead to the variation of the stiffness of the SMPCs. A precise method for estimating the variation rule is rather complicated and the parameters involved also need to be tested by additional experiments. For convenience, the following mixture rule is used to predict the effect of carbon powder on the stiffness for the SMP matrix in the glassy state: where E m and E f are Young's moduli of the SMP matrix and carbon powder, respectively. As shown in Figure 7, the Young's modulus of the carbon powder is obtained by fitting to the test data at room temperature. Similarly, the shear modulus u r in the eight-chain network model is also presented by using the mixture rule, as follows: where u rm and u rf are the shear moduli of the SMP matrix and carbon powder, respectively. As indicated by Zhou et al. [40], the effects of temperature and moisture on the modulus can be given as where E r is Young's modulus of the dry materials in the rubbery state. E u is Young's modulus of the dry materials in the glassy state which can be calculated by Equation (35) . m 1 and m 2 are the coefficients of moisture and temperature, respectively. θ refs is the reference temperature of the material stiffness. The Poisson's ratio μ is considered as temperature dependent and can be given by where ϕ g is the volume fraction of the glassy phase. μ g and μ r are the Poisson's ratios of the glassy and rubbery phases, respectively. As indicated by Qi et al. [17], ϕ g can be obtained by where θ m is the reference temperature and Z is the parameter that restricts the width of the glass transition zone.  Thermomechanical tests were carried out to investigate the shape memory properties by Yang [39]. The following 4 steps were involved in the shape memory cycle. At θ g +10°C, the sample was firstly uniaxially stretched to the pre-strain of 20% at a constant strain rate of 5 × 10 −3 /s. Then the sample was cooled to room temperature in 4 minutes, while the constraint was maintained. After that, a temporary shape was obtained after unloading at a constant strain rate of 5 × 10 −3 /s. Finally, the free recovery was triggered either by heating or by moisture absorption. Note that our study aims to predict the recovery process. Therefore, the other parameters used in our model can be determined from the first 3 steps. The comparison of the model results and test data for the first 3 steps is plotted in Figure 8. It can be found that the stress-strain response of CB0 can be well described by the model. A notable deviation is found in the first step of CB4. The reason should be that the Young's modulus of CB4 predicted by the mixture rule is less than the value tested, as shown in Figure 7.
In summary, the parameters used in the model can be determined by the methods mentioned above. Their values are listed in Table 1.

Model verification
To validate the predictability of the model proposed in the paper, the comparison of the model results with the test data is presented in this section. In the work of Yang [39], the free recovery of CB0 and CB4 was activated either by heating or by moisture absorption. Specifically, the pre-stretched specimen was heated to a high temperature that is well above θ g at a rate of 2°C/min without any constraint to achieve the thermally activated shape recovery. To obtain the moisture-driven SME, the specimen was immersed completely in room temperature water for recovery without any heating process. Note that our model is suitable for describing both the temperature and moisture activated SMEs. Therefore, the heating or moisture absorption triggered free recovery is simulated to contrast with the test data.
As shown in Figure 9, the comparison of the model results with the test data for the recovery ratios of CB0 and CB4 during heating is presented. It should be noted that the test data of the recovery ratios for both CB0 and CB4 are close to 100%. It also means that the blocking effect for the recovery can be ignored during heating. Therefore the equilibrium part 2 shown in Figure 1 is set to invalid in the simulation. It is found that the model results fit well with the test data. The main inconsistency exists in the initial slope of the recovery ratio. The reason is that the viscous flow triggered by temperature is a little faster in our model.
Note that the specimen can not return to its initial shape in the work of Yang [39] even after a long enough immersion time. The final recovery ratio distinctly decreases with the increase of the filler, which is also shown in Figure 10. It should be noted that the final recovery ratio of CB4 is only equal to 90%, after a long enough immersion time (54 hours). It can be regarded as the blocking mechanism that takes effect for the moisture activated SME. Therefore, the effect of hyperelastic part 2 is considered here to denote the blocking mechanism. The pre-deformation of part 2 is ignored at the beginning of the recovery. Once the recovery is started, the hyperelastic part 2 is under the state of compression which is the inverse of hyperelastic part 1. It is assumed that the parameters of part 2 are the same as that of part 1, except for the shear modulus u r . The determination of shear modulus u r2 for part 2 is phenomenological. It is calculated that u r2 = 0.15 u r1 for CB0 and u r2 = 0.2 u r1 for CB4. As shown in Figure 10, the comparison of the model results with the test data for the recovery ratios of CB0 and CB4 in the process of moisture absorption is presented. It can be found that the model results fit well with the test data of CB0. The main inconsistency still exists in the initial slope of the recovery ratio of CB4. Note that the viscous flow triggered by moisture is a little slower in our model. This should be the main reason for this phenomenon.

Conclusions
Concerning the application of thermally activated SMPs in smart structures of aeronautics and astronautics, the electrically conductive fillers are often used to endow the SMPs matrix with electroconductivity. Although the hygro-thermo-mechanical properties and the SMEs of the SMPs filled with nano-carbon powders have already been experimentally investigated intensively, few constitutive models can be used to well describe those phenomena. In our work, a hygro-thermo-mechanical constitutive model for the SMPCs is developed by incorporating the variable moisture into the framework of a thermo-mechanical modeling approach for the first time. Based on thermodynamics with internal state variables, the structural relaxation and viscous flow are well represented in the paper. For convenience, the SMPCs are deemed as isotropic materials in the modeling work. The effect of the nano-carbon powders on the mechanical properties of the SMPCs is considered by using the mixture rule. In the process of immersion, the absorbed moisture is divided into the free phase and bound phase. A moisture diffusion model is proposed to disclose the evolution rule of the moisture diffusion and phase transition. The respective effects of the filler and water on the glass transition temperature are represented based on the theory of free volume. Besides, a method of parameter determination is also presented in the paper. It can be found that the comparisons between the model results and the test data for thermally triggered and moisture triggered SMEs present a good agreement. Thus the effectiveness of the proposed constitutive model is confirmed. Note that the SMEs of the SMPCs are triggered by the ambient temperature or moisture in this study. In practice, the electricity activated SMEs of the SMPs filled with nano-carbon powders are commonly adopted to achieve the driving ability in smart structures. Therefore, the correlation of electric field with Joule heat should be disclosed and the constitutive model for the electricity activated SMEs will be proposed in our future work. Besides, the distribution form of carbon powders in the SMPCs should also be considered.