A 3D Anisotropic Thermomechanical Model for Thermally Induced Woven-Fabric-Reinforced Shape Memory Polymer Composites

Soft robotic grippers offer great advantages over traditional rigid grippers with respect to grabbing objects with irregular or fragile shapes. Shape memory polymer composites are widely used as actuators and holding elements in soft robotic grippers owing to their finite strain, high specific strength, and high driving force. In this paper, a general 3D anisotropic thermomechanical model for woven fabric-reinforced shape memory polymer composites (SMPCs) is proposed based on Helmholtz free energy decomposition and the second law of thermodynamics. Furthermore, the rule of mixtures is modified to describe the stress distribution in the SMPCs, and stress concentration factors are introduced to account for the shearing interaction between the fabric and matrix and warp yarns and weft yarns. The developed model is implemented with a user material subroutine (UMAT) to simulate the shape memory behaivors of SMPCs. The good consistency between the simulation results and experimental validated the proposed model. Furthermore, a numerical investigation of the effects of yarn orientation on the shape memory behavior of the SMPC soft gripper was also performed.


Introduction
Soft robotic grippers can perform complex tasks in uncertain environments where it is difficult for rigid-bodied manipulators to perform tasks requiring flexibility. The main challenge faced by soft robotic grippers is that their stiffness can be variable during the grasping and transfer process [1]. During the grasping process, the stiffness of soft robotic grippers should be low to effectively buffer external impacts in complex environments and be able to adaptively grab irregularly shaped objects. During the object transfer process, the stiffness of soft robotic grippers should be high to maintain their configuration. Thermally induced shape memory polymer composites (SMPCs) are composed of thermally induced shape memory polymers (SMPs) and reinforcement. SMPs are a class of smart materials with the advantages of high stiffness, strength, and driving force, and SMPC-based soft actuators can be easily integrated with other adaptive functional components [2,3]. Furthermore, their stiffness varies with heat stimuli, which makes them a point of attraction in the field of soft robotic grippers [4][5][6][7].
The incorporation of a constitutive model is imperative to capture the thermomechanical and the shape memory behavior of SMPs and SMPCs. A great deal of research has been conducted on the thermomechanical modeling of SMPs. The modeling approaches of SMPs can be divided into two main types: the phase-transition-modeling approach and the viscoelastic modeling approach. The phase-transition-modeling approach assumes that SMPs are composed of a rubbery phase and a glass phase and that shape memory behavior can be realized by provoking the temperature induced transition between the glassy phase and the rubbery phase [8][9][10][11][12][13][14]. The viscoelastic modeling approach is based on the essential thermodynamic properties of SMPs and can be used to describe the entropy elasticity of SMPs above the glass transition temperature and the viscoelastic behavior below the glass transition temperature [15][16][17][18][19][20][21][22]. Moreover, there are also models that combine the concepts of phase transition and viscoelasticity theory [23][24][25][26][27].
Evidently, the constitutive models for SMPCs are more complicated than pure SMPs due to the introduction of reinforcements. Tan et al. [28] developed a constitutive model for unidirectional, continuous, carbon-fiber-reinforced SMPCs within a small strain range, and the effects of the inclination angle and the volume fraction of fiber on the thermomechanical and shape memory properties of SMPCs were investigated. Gu et al. [29] developed a finite strain viscoelastic model for unidirectional, continuous, fiber-reinforced SMPCs with internal state variables and considered the anisotropic thermal properties of SMPCs using a mesomechanics-based method. Li et al. [30] developed a thermomechanical model to describe the temperature-dependent elastic constants of unidirectional carbon-fiberreinforced SMPCs with various fiber volume fractions based on the phase transition theory. Wang et al. [31] developed a constitutive model for unidirectional carbon-fiber-reinforced SMPCs that accounted for interfacial bonding strength. Hong et al. [32] developed a constitutive model based on energy decomposition that accounted for thermal residual stress. Su et al. [25] developed an anisotropic thermomechanical constitutive model for woven-fabric-reinforced SMPCs and investigated the effects of fiber yarn orientation on the shape memory properties of SMPCs. Although the above studies have made great contributions to the modelling of SMPCs, more efforts are needed to improve the accuracy of the corresponding models and facilitate their further development In this paper, a novel 3D anisotropic thermomechanical model for thermally induced woven-fabric-reinforced SMPCs is developed based on Helmholtz free energy decomposition and the second law of thermodynamics. The total Helmholtz free energy is decomposed into the isotropic part of the matrix and the anisotropic part of the woven fabric. The energy of matrix part can be further decomposed into a hyperelastic part and a viscoelastic part, and the energy of the woven fabric part can be further decomposed into a hyperelastic part through the action of fiber stretching and fiber-fiber shearing. Furthermore, the rule of mixtures is modified to consider the stress distribution in the phases of the SMPCs, and stress concentration factors are introduced to consider fiber-matrix and fiber-fiber shearing interactions.
The paper is arranged as follows. Section 2 introduces a 3D anisotropic thermomechanical model for thermally induced woven-fabric-reinforced SMPCs. In Section 3, the model determination protocols for the material parameters' and model's verification are presented. Finally, the conclusions are drawn in Section 4.

Kinematics
In this section, a 3D anisotropic thermomechanical model for thermally induced woven-fabric-reinforced SMPCs is proposed based on Helmholtz free energy decomposition and the second law of thermodynamics along with the consideration of the temperaturedependent interfacial effects of the SMPCs during the shape memory cycle. To keep matters simple, the anisotropic thermal deformation and residual stress of the SMPCs were not considered. The corresponding rheological model is shown in Figure 1. The subscripts m, e, v, and f denote the single spring element of the matrix part, the spring in the Maxwell element, the dashpot in the Maxwell element, and the single spring element of the fabric part, respectively.  In this case, perfect bonding between the fabric and matrix is assumed, which leads to the following relation: where F is the total deformation gradient.
F can be further decomposed into an elastic part e F and a viscous part v F : ev  F F F (2) Then, the Cauchy-Green deformation tensors can be expressed as follows: The Green-Lagrange tensors are denoted as follows: where I is the second-order unit tensor.
The invariants of C and e C can be defined as follows: The yarn orientation unit vectors in the original configuration are denoted as 0 a and 0 b , as shown in Figure 2a, and the yarn orientation unit vectors in the current configuration are denoted as a , b , which can be formulated as follows: In this case, perfect bonding between the fabric and matrix is assumed, which leads to the following relation: where F is the total deformation gradient. F can be further decomposed into an elastic part F e and a viscous part F v : Then, the Cauchy-Green deformation tensors can be expressed as follows: The Green-Lagrange tensors are denoted as follows: where I is the second-order unit tensor. The invariants of C and C e can be defined as follows: The yarn orientation unit vectors in the original configuration are denoted as a 0 and b 0 , as shown in Figure 2a, and the yarn orientation unit vectors in the current configuration are denoted as a, b, which can be formulated as follows: The main modes of the deformation of the fabrics mainly include tension along the yarn orientation and shearing between the weft and warp yarn. To analyze the mechanical behavior of the fabrics, the invariants I 4 and I 7 are introduced [33]:  I are introduced [34]:

Constitutive Model for Matrix SMPs
It can be seen from Figure 1 Based on Equation (8) The Clausius-Duhem inequality based on the second law of thermodynamics can be expressed as follows [10]: where S is the second Piola-Kirchhoff stress tensor,  is entropy, q is the heat flux vector, and  represents the gradient operator.

Constitutive Model for Matrix SMPs
It can be seen from Figure 1 that the total Helmholtz free energy of SMPs ψ matrix can be decomposed into ψ m and ψ e : Based on Equation (8), . ψ matrix can be formulated as The Clausius-Duhem inequality based on the second law of thermodynamics can be expressed as follows [10]: where S is the second Piola-Kirchhoff stress tensor, η is entropy, q is the heat flux vector, and ∇ represents the gradient operator. Substituting Equations (4) and (9) into (10) leads to F v · F v −1 is the velocity gradient tensor. Equation (11) must account for arbitrary thermodynamic processes. Therefore, the second Piola-Kirchhoff stress of matrix S matrix and the velocity gradient tensor of the viscoelastic part can be expressed as follows: where ζ is a temperature-dependent viscosity parameter, which can be expressed as follows: where subscript l denotes parameters at T l , which is the lowest temperature below the glass transition temperature (T g ) in a thermomechanical test. Subscript h denotes the parameters at T h , which is the highest temperature above T g in a thermomechanical test. A is a material parameter, while φ is a weight function, which can be expressed as follows: Sensors 2023, 23, 6455 where g is a material parameter, and T r is the reference temperature.
Here, the Mooney-Rivlin model is adopted for ψ m and ψ e : where C 10 m , C 10 e , C 01 m , C 01 e , D m , and D e are temperature-dependent parameters, which can be expressed as follows: Here, the isochoric flow assumption (|F v | = 1) leads to the following relation: Based on Equations (12)- (17), the Cauchy stress of the matrix SMPs and the velocity gradient of the viscoelastic part can be expressed as follows:

Constitutive Model for Fabric
The woven fabric yarn can be stretched along the fiber orientation and sheared at a wider angle when a load is applied. Therefore, the total Helmholtz free energy of the fabric ψ f abric can be decomposed into a fiber-stretching part ψ ten and a fiber-shearing part ψ shr : Here, the polynomial function of I 4 a , I 4 b was adopted for ψ ten , and I 4 a , I 4 b are assumed to be equal to 1 under the following condition [33]: where k 1 , k 2 , k 3 are material parameters.
Here, the polynomial function of I 7 was adopted for ψ shr [33]: where k 4 , k 5 , k 6 are material parameters.  (10), the Cauchy stress of fabric can be expressed as follows: Generally, the total stress of SMPCs σ can be formulated using the rule of mixtures: where v m , v f are the volume fractions of the matrix SMPs and fabric, respectively, and v m + v f = 1. However, as demonstrated in previous experiments [25,31,34], the total stress of SMPCs cannot be accurately described by the rule of mixtures alone since a phase transition of the SMPs will occur with a temperature change, and the stress distribution in the phases of the SMPCs are also related to their microstructure. Therefore, an effective temperaturedependent fabric volume fraction of fabric is introduced based on phase transition theory, and Equation (23) is modified as follows: where v f _re f denotes the reference volume fractions of the fabric, and γ inter denotes the stress concentration factors to consider fiber-fiber and matrix-fiber shearing interactions, which are assumed to be equal to 1 when stretching along the yarn.

Methods for Determining Matrix SMPs Material Parameters
The material parameters in the matrix SMP part can be determined following the protocol outlined in our previous works [26].

Methods for Determining Remaining Material Parameters
With the parameters of the matrix SMP part determined, the parameters of k 1 , k 2 , k 3 in ψ ten and k 4 , k 5 , k 6 in ψ shr and v f _re f can be determined using the following procedure: 1.
Based on the Equations (18) and (22)- (24), k 1 , k 2 , k 3 can be determined by fitting the stress-strain curve of the 0 • SMPCs at T l with an initial assumption of v f _re f = v f .
v f _re f can be obtained by fitting the curve of the strain of the 0 • SMPCs in the loading step, cooling step, and unloading step in the shape memory cycle.

4.
Compare the fitting results with the experimental data. If good consistency has been achieved, then these parameters are determined. If not, modify the initial estimate of the constant v f _re f and return to (1).

5.
With the above parameter determined, γ inter can be obtained by fitting the curve of the strain of the 45 • SMPCs in the loading step, cooling step, and unloading step in the shape memory cycle.
According to the above parameter determination protocol, the parameters of k 1 , k 2 , k 3 in ψ ten and k 4 , k 5 , k 6 in ψ shr and v f _re f are shown). All determined material parameters are  Table 1, and the fitting results regarding k 1 , k 2 , k 3 , k 4 , k 5 , k 6 , v f _re f , and γ inter are shown in Figures 3-6. 3.
f_ref v can be obtained by fitting the curve of the strain of the 0° SMPCs in the loading step, cooling step, and unloading step in the shape memory cycle. 4. Compare the fitting results with the experimental data. If good consistency has been achieved, then these parameters are determined. If not, modify the initial estimate of the constant f_ref v and return to (1).

5.
With the above parameter determined, inter  can be obtained by fitting the curve of the strain of the 45° SMPCs in the loading step, cooling step, and unloading step in the shape memory cycle.
According to the above parameter determination protocol, the parameters of 1 2 3 , , k k k in ten  and 4 5 6 ,, k k k in shr  and f_ref v are shown). All determined material parameters are summarized in Table 1    f_ref v can be obtained by fitting the curve of the strain of the 0° SMPCs in the loading step, cooling step, and unloading step in the shape memory cycle. 4. Compare the fitting results with the experimental data. If good consistency has been achieved, then these parameters are determined. If not, modify the initial estimate of the constant f_ref v and return to (1).

5.
With the above parameter determined, inter  can be obtained by fitting the curve of the strain of the 45° SMPCs in the loading step, cooling step, and unloading step in the shape memory cycle.
According to the above parameter determination protocol, the parameters of 1 2 3 , , k k k in ten  and 4 5 6 ,, k k k in shr  and f_ref v are shown). All determined material parameters are summarized in Table 1

Model Verification
The proposed model was implemented in the commercial finite element software package ABAQUS/Standard via a user material subroutine (UMAT) to simulate the shape memory tests carried out in the study by Su et al. [25]. Eight-node linear hexahedron continuum elements (C3D8) are used. Loading and boundary conditions are set according to the design of the experiment. Heat transfer is ignored here, and the temperature is applied through the predefined field. First, the shape memory tests of the matrix SMPs were simulated. The temperature histories of the simulation in the cooling and reheating steps are consistent with the experimental data, as shown in Figure 7. During the shape memory cycle tests, the temperature of the experimental specimens cannot immediately reach the temperature of the experimental equipment because it takes time for the specimens to reach thermal equilibrium. Therefore, a temperature lag of 3.7 • C for SMP was assumed. The simulation results regarding the matrix SMPs are shown in Figure 8. It can be seen that the simulation results are in good agreement with the experimental data, thus verifying the effectiveness of the proposed model in predicting the shape memory behavior of matrix SMPs.
perimental data, as shown in Figure 7. During the shape memory cycle tests, the temperature of the experimental specimens cannot immediately reach the temperature of the experimental equipment because it takes time for the specimens to reach thermal equilibrium. Therefore, a temperature lag of 3.7 °C for SMP was assumed. The simulation results regarding the matrix SMPs are shown in Figure 8. It can be seen that the simulation results are in good agreement with the experimental data, thus verifying the effectiveness of the proposed model in predicting the shape memory behavior of matrix SMPs.  Second, the shape memory tests of the SMPCs with an initial yarn orientation of 0/90° and ±45° were simulated. The temperature histories of the simulation in the cooling and reheating steps are consistent with the experimental data, as shown in Figure 7. Temperature lag values of 3.2 °C and 3.7 °C for the 0/90° and ±45° woven-fabric-reinforced SMPCs perimental data, as shown in Figure 7. During the shape memory cycle tests, the temperature of the experimental specimens cannot immediately reach the temperature of the experimental equipment because it takes time for the specimens to reach thermal equilibrium. Therefore, a temperature lag of 3.7 °C for SMP was assumed. The simulation results regarding the matrix SMPs are shown in Figure 8. It can be seen that the simulation results are in good agreement with the experimental data, thus verifying the effectiveness of the proposed model in predicting the shape memory behavior of matrix SMPs.  Second, the shape memory tests of the SMPCs with an initial yarn orientation of 0/90° and ±45° were simulated. The temperature histories of the simulation in the cooling and reheating steps are consistent with the experimental data, as shown in Figure 7. Temperature lag values of 3.2 °C and 3.7 °C for the 0/90° and ±45° woven-fabric-reinforced SMPCs Second, the shape memory tests of the SMPCs with an initial yarn orientation of 0/90 • and ±45 • were simulated. The temperature histories of the simulation in the cooling and reheating steps are consistent with the experimental data, as shown in Figure 7. Temperature lag values of 3.2 • C and 3.7 • C for the 0/90 • and ±45 • woven-fabric-reinforced SMPCs were assumed. The simulation results are shown in Figure 9. It can be seen that the simulation accurately reproduces the shape memory recovery behavior of the SMPCs, which demonstrates the validity of the proposed model in predicting the shape memory behavior of SMPCs. The deviation between the simulation results and the experimental data in the final recovery stage might have been caused by interfacial failure and internal stress, and this deviation will be modified by introducing anisotropic thermal and stress internal stress and interfacial failure in our future research. Finally, the shape memory behavior of the SMPC gripper part was simulated to investigate the effects of yarn orientation. The main deformation mode of a soft robotic gripper is flexural deformation, as shown in Figure 10, since the human hand's ability to grasp objects is achieved through the movement of bones and joints. Here, the grasping of objects via hands is simulated through the flexural deformation of the individual gripper part with dimensions of 70 mm × 25 mm × 2 mm. Eight-node linear hexahedron continuum elements (C3D8) and six-node linear triangular prism elements (C3D6) are used. The applied position of the load and the boundary conditions are shown in Figure 11. The left area of the fixed boundary and the right area of the displacement boundary represent the two finger bones, and the middle part represents the flexible joint.
were assumed. The simulation results are shown in Figure 9. It can be seen that the simulation accurately reproduces the shape memory recovery behavior of the SMPCs, which demonstrates the validity of the proposed model in predicting the shape memory behavior of SMPCs. The deviation between the simulation results and the experimental data in the final recovery stage might have been caused by interfacial failure and internal stress, and this deviation will be modified by introducing anisotropic thermal and stress internal stress and interfacial failure in our future research. Figure 9. Comparison between experimental and simulation results regarding shape recovery for SMPCs.
Finally, the shape memory behavior of the SMPC gripper part was simulated to investigate the effects of yarn orientation. The main deformation mode of a soft robotic gripper is flexural deformation, as shown in Figure 10, since the human hand's ability to grasp objects is achieved through the movement of bones and joints. Here, the grasping of objects via hands is simulated through the flexural deformation of the individual gripper part with dimensions of 70 mm × 25 mm × 2 mm. Eight-node linear hexahedron continuum elements (C3D8) and six-node linear triangular prism elements (C3D6) are used. The applied position of the load and the boundary conditions are shown in Figure 11. The left area of the fixed boundary and the right area of the displacement boundary represent the two finger bones, and the middle part represents the flexible joint.    Finally, the shape memory behavior of the SMPC gripper part was simulated to investigate the effects of yarn orientation. The main deformation mode of a soft robotic gripper is flexural deformation, as shown in Figure 10, since the human hand's ability to grasp objects is achieved through the movement of bones and joints. Here, the grasping of objects via hands is simulated through the flexural deformation of the individual gripper part with dimensions of 70 mm × 25 mm × 2 mm. Eight-node linear hexahedron continuum elements (C3D8) and six-node linear triangular prism elements (C3D6) are used. The applied position of the load and the boundary conditions are shown in Figure 11. The left area of the fixed boundary and the right area of the displacement boundary represent the two finger bones, and the middle part represents the flexible joint.    Finally, the shape memory behavior of the SMPC gripper part was simulated to investigate the effects of yarn orientation. The main deformation mode of a soft robotic gripper is flexural deformation, as shown in Figure 10, since the human hand's ability to grasp objects is achieved through the movement of bones and joints. Here, the grasping of objects via hands is simulated through the flexural deformation of the individual gripper part with dimensions of 70 mm × 25 mm × 2 mm. Eight-node linear hexahedron continuum elements (C3D8) and six-node linear triangular prism elements (C3D6) are used. The applied position of the load and the boundary conditions are shown in Figure 11. The left area of the fixed boundary and the right area of the displacement boundary represent the two finger bones, and the middle part represents the flexible joint.   Initially, the temperature was set at 70 • C, and a rotation boundary condition of UR2 =−π/4 was applied. Then, the temperature was set to decrease from 70 • C to 20 • C at a cooling rate of −2 • C/min, which allowed for the shape to be maintained, followed by unloading when the temperature reached 20 • C. Finally, the temperature was set to increase from 20 • C to 70 • C at a heating rate of 2 • C/min in the free state. The simulation results are shown in Figures 12 and 13.
It can be seen from Figure 12 that the recovery rate of the 0/90 • SMPCs is faster than that of the ±45 • SMPCs owing to their larger stored recovery stress, as shown in Figure 13. However, the shape fixity ratio of the 0/90 • SMPCs is lower than that of the ±45 • SMPCs, as shown in Figure 14, because the fabric used in this paper was carbon fabric; additionally, as the fiber orientation angle decreased, the reinforcing effect of the carbon fiber was enhanced; thus, the carbon fiber was rendered elastic. When SMPCs are stretched at a high temperature and unloaded after cooling, the larger rebound stress of the fiber reduces the shape fixation rate of the 0/90 SMPCs compared to the ±45 SMPCs. Therefore, it is necessary to comprehensively consider the response rate and accuracy of an SMPC soft gripper to design an appropriate fiber orientation in practical applications.
Initially, the temperature was set at 70 °C , and a rotation boundary condition of UR2 =−π/4 was applied. Then, the temperature was set to decrease from 70 °C to 20 °C at a cooling rate of −2 °C /min, which allowed for the shape to be maintained, followed by unloading when the temperature reached 20 °C . Finally, the temperature was set to increase from 20 °C to 70 °C at a heating rate of 2 °C /min in the free state. The simulation results are shown in Figures 12 and 13.  It can be seen from Figure 12 that the recovery rate of the 0/90° SMPCs is faster than that of the ±45°S MPCs owing to their larger stored recovery stress, as shown in Figure 13. However, the shape fixity ratio of the 0/90° SMPCs is lower than that of the ±45° SMPCs, as shown in Figure 14, because the fabric used in this paper was carbon fabric; additionally, as the fiber orientation angle decreased, the reinforcing effect of the carbon fiber was enhanced; thus, the carbon fiber was rendered elastic. When SMPCs are stretched at a high Initially, the temperature was set at 70 °C , and a rotation boundary condition of UR2 =−π/4 was applied. Then, the temperature was set to decrease from 70 °C to 20 °C at a cooling rate of −2 °C /min, which allowed for the shape to be maintained, followed by unloading when the temperature reached 20 °C . Finally, the temperature was set to increase from 20 °C to 70 °C at a heating rate of 2 °C /min in the free state. The simulation results are shown in Figures 12 and 13.  It can be seen from Figure 12 that the recovery rate of the 0/90° SMPCs is faster than that of the ±45°S MPCs owing to their larger stored recovery stress, as shown in Figure 13. However, the shape fixity ratio of the 0/90° SMPCs is lower than that of the ±45° SMPCs, as shown in Figure 14, because the fabric used in this paper was carbon fabric; additionally, as the fiber orientation angle decreased, the reinforcing effect of the carbon fiber was enhanced; thus, the carbon fiber was rendered elastic. When SMPCs are stretched at a high temperature and unloaded after cooling, the larger rebound stress of the fiber reduces the shape fixation rate of the 0/90 SMPCs compared to the ±45 SMPCs. Therefore, it is necessary to comprehensively consider the response rate and accuracy of an SMPC soft gripper to design an appropriate fiber orientation in practical applications.

Conclusions
An anisotropic thermomechanical model for thermally induced woven-fabric-reinforced SMPCs was developed based on Helmholtz free energy decomposition and the second law of thermodynamics. The total Helmholtz free energy of the SMPCs was decomposed into an isotropic visco-hyperelastic matrix part and an anisotropic hyperelastic

Conclusions
An anisotropic thermomechanical model for thermally induced woven-fabric-reinforced SMPCs was developed based on Helmholtz free energy decomposition and the second law of thermodynamics. The total Helmholtz free energy of the SMPCs was decomposed into an isotropic visco-hyperelastic matrix part and an anisotropic hyperelastic fabric part. The stress distribution of the phases in the SMPCs was described using a modified rule of mixtures based on phase transition theory, and stress concentration factors were introduced to consider fiber-matrix and fiber-fiber shearing interactions. The shape memory tests presented in the study by Su et al. [25] were simulated, and a comparison between the simulation results and the experimental data verified the effectiveness of the proposed model. Finally, the effects of yarn orientation on the shape memory behavior of the soft robotic gripper were investigated using the proposed model.