Multilayer modeling framework for analyzing thermo-mechanical properties and responses of twisted and coiled polymer actuators

The twisted and coiled polymer actuator (TCPA) has a complex multi-scale structure consisting of crystalline micro-fibrils and an amorphous matrix at the micro-scale, which are organized into a macro-scale fiber. When the polymer fiber undergoes twisting and coiling, its mechanical and thermal properties become variable. In this study, we developed a multi-layer modeling framework capable of accurately predicting the effective mechanical and thermal properties, as well as the thermo-mechanical responses of the TCPA. Our numerical results demonstrate that the effective mechanical and thermal properties of the TCPA are influenced by the radius and twisting angle of the polymer fiber. By analyzing the precise mechanical and thermal properties, the numerical calculated driving responses exhibit good agreement with experimental data. We also examined the influence of initial helical radius, helical pitch and fiber radius on the driving responses of the TCPA. The proposed numerical model can be further utilized to optimize the driving responses of the TCPA by adjusting geometric parameters and the twisting angle of the polymer fiber.


Introduction
Artificial muscles have emerged as a type of soft actuation that imitate the compliant behavior of biological muscles.They serve to address the limitations of conventional actuators in bio-inspired applications [1][2][3][4][5][6].These artificial muscles rely on external stimuli, such as voltage, current, pressure, temperature, or light, to achieve reversible contraction, expansion, and rotation.As a result, they find applications in various fields, including smart exoskeletons, prosthetics [7,8], wearables [9], intelligent surgical tools [10], and even humanoid robots [11,12].Depending on the driving materials and characteristics, artificial muscles can be categorized as piezoelectric materials [13], shape memory alloy [14], shape memory polymer [15], dielectric elasticity (DE) [16], electrochemical drive [17], among others.In this study, we mainly focus on fibrous artificial muscles.
Fibrous artificial muscles consist of polymer fibers, which are available in various forms, including artificial polymer fiber (e.g.Polyethylene and Nylon 66 [18,19]), artificial inorganic fibers (e.g.carbon nanotube and graphene fibers [20,21]), natural fibers (e.g.silk and cotton [22,23]), and composite fibers [24].Semi-crystalline fibers, which exhibit highly oriented structures along their length and exhibit anisotropic thermal expansion behavior, are commonly used as the constituent material for fibrous artificial muscles.In recent years, a specific type of thermally driving fibrous artificial muscle called the twisted and coiled polymer actuator (TCPA) has gained increasing attention.TCPA offers numerous advantages, including a high power-to-weight ratio, large contractions, lightweight design, and cost effectiveness.The TCPA is made from widely available low-cost polymer fibers, such as nylon fishing lines and sewing threads [25].The TCPA can contract by up to 49%, lift loads over 100 times heavier than human muscle of the same length and weight, and generate 5.3 kilowatts of mechanical work per kilogram of muscle weight [25].Nylon 66 used as the TCPA precursor fiber is an anisotropic material with an axial thermal shrinkage rate of 4% when affected by temperature [25].The effect of the twisting angle on the shrinkage rate can reach 18% [26].This development has generated substantial interest in the field due to its potential applications and advantages.
In recent years, there have been several experimental studies focused on improving the driving performance of the TCPAs, such as the work by Higueras-Ruiz et al [27] and references therein.Temperature and twisting angle can improve the driving response of the TCPA [26,28].The helical pitch of the TCPA has a large effect on the drive performance.Pyo et al [29] experimentally studied the effect of helical pitch on driving performance, and the results showed that when the helical pitch increased, TCPA driving performance increased first and then decreased.The spring index, the ratio of mean coil diameter to the fiber diameter, affects the driving strain of TCPA by 15% when the TCPA spring index increases from 1.1 to 1.7 [25].In addition, TCPA's heating and cooling mode and speed, heating voltage, current and power will all have an impact on its driving performance [30][31][32][33].
However, based on complex geometric configurations of the TCPA, there are relatively few theoretical and numerical works to study the effective mechanical and thermal properties of the TCPA.Yang and Li [34] established a theoretical framework based on a multi-scale approach to obtain the thermomechanical driving response of artificial muscles.They incorporated Castigliano's second theorem to account for the combined effect of the axial applied force and recovered torque on the macro-scale artificial muscle's thermally triggered tensile actuation.Hunt et al [28] expanded upon Yang's model by conducting a finite element analysis using COMSOL software.In Hunt's study, a local cylindrical coordinate system was established in COMSOL's Curvilinear Coordinates module to incorporate the material and thermal properties of the TCPA.Wu and Zheng [35] proposed an elastic-rod-theorybased model to capture the quantitative relationship between tensile actuation and fabrication load for twisted and coiled polymers.Tang [36] presented a numerical model to evaluate the thermal stress in twisted fibers.Their model successfully predicted the deformation of the TCPA under varying temperatures, and its accuracy was validated against experimental results.Liu et al [37] analyzed the effects of preparation conditions, actuation load, and operating temperature on the tensile stroke of twisted and coiled artificial muscles.Although these models present a wide range of driving responses of the TCPA, the limitation is that these models cannot determine the effective mechanical and thermal properties which vary with the radius and twisting angle of the polymer fiber.
In this paper, we develop a fully three-dimensional multilayer model for the characterization of the thermo-mechanical behaviors of the TCPA.The model is based on the multi-scale theory and is implemented numerically using a combination of Abaqus software and Python code.The Python code is utilized to enhance the capabilities of finite element analysis and tackle complex tasks that are difficult or impossible to accomplish manually in Abaqus.The primary objective of this model is to evaluate the effective mechanical and thermal properties of the TCPA and predict its thermo-mechanical responses.The paper is organized as follows: In section 2, we provide an overview of the multi-scale theoretical model for TCPA.Section 3 presents the multi-layer model that is specifically designed to analyze the effective mechanical and thermal properties of the TCPA.Furthermore, in section 4, we discuss the thermomechanical responses of the TCPA for a range of parameter variations.Finally, in section 5, we summarize the findings of this research and present our conclusions.Overall, this paper aims to contribute to the understanding and characterization of the thermo-mechanical behaviors of the TCPA through the development and implementation of a multi-layer modeling framework.

Multi-scale theoretical model for TCPA
To analyze the thermo-mechanical properties and responses of the TCPA, we numerically solve a multi-scale theoretical model of the TCPA which provided a top-down strategy, spanning from macro-scale helical spring analysis down to molecular level chain interaction [34].In the multi-scale theoretical framework, the mechanical analysis of twisted and coiled fibers can be divided into four levels: macro-scale, meso-scale, micro-scale and nano-scale.The nano-and micro-scales can be utilized with direct input the material properties from existing literatures [38][39][40][41].We will focus on the theories of macroscale and meso-scale.For completeness, a summary of the formulation is provided below.

Nano-and micro-scale modeling for TCPA
The multi-scale theoretical model begins by determining the material properties of the fiber at the nano-and micro-scale levels.Previous research has provided information on the material properties of a composite plate consisting of crystalline micro-fibrils and an amorphous matrix [38,39].Based on the shear lag theory developed by Cox [40], the longitudinal effective axial Young's modulus E 1 and the effective thermal expansion coefficient α 1 at the micro-scale can be determined by where ) , l f is the length of the discontinuous reinforcing phase, and r f is the radius of microfibrils.
, G f 12 = C c 66 and α f 1 represent Young's modulus and thermal expansion coefficient of the fiber phase in the micro-scale.V f is the volume fraction.E m 1 and α m 1 represent Young's modulus and thermal expansion coefficient of the matrix phase in the micro-scale.K is a constant which is a measurement of the reinforcing efficiency.C c ij are elastic moduli of the nano-scale crystals which can be obtained from theoretical calculations and experimental measurements [34].The relationship of the axial direction modulus between nano-scale and micro-scale is related by a scale factor λ. The transverse direction modulus almost unchanged.
is the shear modulus of the matrix phase in the micro-scale, and where x m g is the glassy phase volume fraction within the amorphous matrix.S is a constant and T g represents the glass transition temperature.
Referring to the micromechanical relationship proposed by Bogetti and Gillespie [41], the transverse material parameters can be expressed as where k T is the effective bulk modulus.v f 12 and α f 2 represent Poisson's ratio and thermal expansion coefficient of the microfiber phase in the micro-scale, v m 12 and α m 2 represent Poisson's ratio and thermal expansion coefficient of the amorphous matrix phase in the micro -scale.G 23 is the material parameters related to G f 23 and G m 23 .v 12 is the material parameters related to v f 23 and v m 23 .

Meso-scale modeling
When the polymer fibers are twisted, the fibers are assumed to be a layered cylinder with helical anisotropy, as shown in figure 1. Considering the concentric hollow cylinders, radial displacement u, axial displacement w and tangential displace- The corresponding strain components are Based on the axial symmetry of the structure, the equilibrium differential equation can be expressed as The constitutive relation for a helical anisotropy cylinder is The solution of the radial displacement u can be expressed as where coefficients C 1 and C 2 are determined by the continuity and boundary conditions.Figure 2 illustrates the continuity and boundary conditions of the multi-layer concentric hollow cylinders.The superscript 'k' represents the mechanical quantities of the kth layer cylinder.The radius of the kth layer is given by r k = r 0 + k rf−r 0 N , where r 0 is the inner radius, which is a very small value.The relationship between the twisting angle and the layer numbers can be expressed as θ k f = arctan( r k rf tan(θ f )), as shown in figure 3. Here, θ f is the twisting angle on the outermost cylinder.For simplicity, the superscript 'k' is omitted in equation ( 10) and subsequent expressions.The continuity and boundary conditions of the multi-layer concentric hollow cylinders hold From the given r 0 and r k , all C 1 and C 2 coefficients of the kth layer in the displacement expression can be obtained using equation (A.1) (refer to appendix A for more details).Subsequently, the stress components can be determined based on the literature [38].
The theoretical analysis of twisted and coiled polymer fiber (figure 3(a)) involves decomposing the twisted polymer fiber into concentric hollow cylinders (figure 3(b)).Each cylinder is then hypothetically cut and opened into a layered structure.Existing literature can provide input parameters for nanoand micro-scale analyses.The analysis of each lamina on the meso-scale is based on the classical laminate theory [42].When the polymer fiber is twisted, the Cartesian coordinate system {1, 2, 3} is used, where 1 denotes the principal directions of the off-axis behavior.θ f represents the twisting angle on the outermost surface of the cylinder, as depicted in figure 3(b).When θ f = 0, namely the polymer fiber is not twisted which represents an on-axis lamina, the direction of the cylindrical coordinate system is the same as the Cartesian coordinate system, as shown in figure 3(c).When θ f ̸ = 0, namely the polymer fiber is twisted, the angle between the cylindrical coordinate and the Cartesian coordinate is the twisting angle θ f , as shown in figure 3(d).Therefore, the stressstrain relationship of the on-axis lamina can be expressed as where [K] is the stiffness matrix of the composite lamina.
According to the composite laminate theory, the stiffness matrix of the off-axis lamina can be obtained by transforming the stiffness matrix of the on-axis lamina where C ij is the stiffness matrix of the on-axis lamina, and m = cos θ f , n = sin θ f .The relationship of thermal expansion coefficient between the on-axis and off-axis lamina can be expresses as

Effective mechanical and thermal properties of macro-scale TCPA
Based on the analysis of a helically anisotropic cylinder in section 2.2, we now present the relaxation method to calculate the effective mechanical and thermal properties of the TCPA [38].Given an axial extensional strain ω 0 of the cylinder, the axial resultant force can be determined, and the effective Young's modulus and Poisson's ratio of each layer cylinder are obtained.Similarly, by applying shear strain ν 0 to the cylinder, the resultant torque is obtained, and the effective shear modulus of each layer cylinder is determined.The resultant forces and torques are represented by the matrix below where τ = 2 ´rf r 0 σ θz π r 2 dr and f = 2 ´rf r 0 σ z π rdr.The coefficients A, B, C and D can be expressed in terms of ω 0 and ν 0 Specially, When resultant torque is zero, the effective Young's modulus and Poisson's ratio can be obtained.When the resultant force is zero, the effective shear modulus can be obtained.The effective Young's modulus, Poisson's ratio, and shear modulus of each cylinder can be expressed as follows From the relations in equation ( 21), the effective mechanical properties of the TCPA can be expressed as For effective thermal expansion coefficient, let ν 0 = ω 0 = 0, and the thermal torque τ T and force f T are determined for a given increment ∆T.The corresponding thermal strain components are determined by the following equation Finally, the effective thermal expansion coefficient can be expressed as

Numerical implementation
In this section, we conduct multi-layer modeling framework to solve numerically the above multi-scale model of the TCPA and to investigate the effective mechanical and thermal properties, as well as the thermo-mechanical responses of the TCPA.
Our aim is to demonstrate that these simulations can accurately capture the actuation responses of the TCPA.To perform these simulations, we utilize the commercial finite element software Abaqus and Python code within the Abaqus environment [43].Figure 4 schematically depicts the multi-layer numerical model of the TCPA, which consists of many concentric hollow cylinders.Each layer cylinder exhibits different mechanical and thermal properties when the polymer fibers are twisted.Building upon the previous section, we use Python code to build a multi-layer numerical model of twisted fiber by assigning material parameters in each layer cylinder.
To model the driving response of the TCPA, we utilize three-dimensional coupled temperature-displacement elements (C3D8T) in the Abaqus.The outside surface of the TCPA is subjected to a thermal load, which can take the form of uniform loading [18].Displacement boundary conditions are imposed on the outside surface of the fiber, allowing the actuator to only swell in the length direction.Once the material properties of each layer cylinder are determined, the equilibrium equations of the system are solved at each step using a nonlinear solver based on the Newton-Raphson method.The flowchart outlining the multi-layer modeling framework is shown in figure 5.
The multi-layer model predicts the effective mechanical and thermal properties as a function of fiber radius, as shown in figure 6.In figure 6(a), it can be observed that the effective Young's modulus increases with the fiber radius.The concentric hollow cylinders exhibit a sequential increase in the effective Young's modulus from the inside to the outside.Additionally, figure 6(a) demonstrates the influences of the numbers of concentric hollow cylinder layers on the effective Young's modulus.The radius of polymer fiber is 0.571 mm, helical radius is 1.083 mm, θ f = 37.82 • , helical pitch is 1.23 mm, and the winding angle of the TCPA is 8.45 • .The samples were heated from 30 • C to 80 • C. The results in the figure indicate that the effective Young's modulus of the twisted cylinder increases with an increase in the number of concentric hollow cylinder layers.When the number of layers exceeds 30, the value of the effective Young's modulus no longer changes with the number of layers.Figure 6(b) illustrates that the effective Poisson's ratio increases with an increase in the fiber radius.The number of layers has no effect on the effective Poisson's ratio as the layer numbers exceeds 30, which follows a similar trend to the effective Young's modulus.Figure 6(c) shows the effect of the radius on the variation of the effective thermal expansion coefficient with the layer numbers.It can be seen that the change in the effective thermal expansion coefficient is small with the layer numbers.
The effective mechanical and thermal properties of the TCPA as a function of the twisting angle are predicted using a multi-layer finite element method.These properties are plotted in figure 7 for several representative layer numbers.Figure 7(a) shows the variation of the effective Young's modulus with the twisting angle.Among all effective Young's modulus values, there exists an optimal value that maximizes Young's modulus of the TCPA.For example, we found the maximum effective Young's modulus at a twisting angle of 40 • .Figure 7(b) depicts the variations of the effective Poisson's ratio with the twisting angle.Similar to the effective Young's modulus, the effective Poisson's ratio also reaches a maximum value at a twisting angle of 40 • .Figure 7(c) shows the effect of twisting angle on the variation of the effective thermal expansion coefficient with the layer numbers.The variation trend is consistent with the effect of radius.
From the numerical prediction results in figure 7(a), it can be seen that the largest effective Young's modulus occurs as the twisting angle of the TCPA is 40 • .The largest Young's modulus corresponds to the largest driving response.Its mechanism is based on the macro-scale energy equation of Castigliano's theorem [44].When a TCPA is subjected to the axial force F, the recovery moment M rec can be generated.More details can be found in section 3 of the work by Yang and Li [34].The macro-scale energy can be expressed as   f is the cross-sectional area of the fiber and r f is the fiber radius.f 11 , f 12 and f 22 are the collected constants, and are expressed as where L c is the final length of the TCPA.α c is the helical angle, and R is the helical radius.The normal force

and the bending moment
The driving displacement δ of a TCPA is directly found by applying Castigliano's theorem [34,44]: When the effective moduli of TCPA is obtained numerically, the driving displacement can be obtained by using equation (33) as shown in figure 7(d).From figure 7(d), it is seen that the driving response reaches optimal value when the twisting angle is 40 • .The key point of the above energy method is to obtain the effective moduli Ē and Ḡ, which can be accurately calculated from the proposed multi-layer model.To ensure accurate numerical solutions, it is necessary to study the convergence behavior of the driving forces.Figure 8 presents the convergence behavior of the driving forces for different layer numbers (5, 10, 20, 30, 40, and 50).It is observed that when the layer numbers of concentric hollow cylinders reach 30, the driving forces converge, ensuring accuracy in the calculations.

Thermo-mechanical responses of the TCPA and discussions
In this section, we employ the multi-layer finite element model as described in sections 3 to analyze and understand the thermo-mechanical responses of the TCPA.Whenever possible, we compare the results of our numerical simulations with those obtained from experimental data.

Comparison with experimental results
In this section, we present four sets of experimental data to validate and evaluate the reliability of the numerical implementation.Figure 9 illustrates a comparison between the model predictions and the experimental data of Yang and Li [34] for uniaxial thermal actuation at several temperatures.The driving displacements are plotted as a function of the temperature.Our model results are in good agreement with the experimental observations reported in the literature.It is evident from figure 9 that the driving strain increases consistently as the temperature rises.
To further assess the performance of the proposed numerical model, we analyze the driving forces of the TCPA in relation to different heating wire powers and heating times, as presented in figures 10 and 11, respectively.The corresponding parameters used on these simulations are shown in table 1.The thermally driven test setup and sample picture  Comparison of the variations of the driving strain with temperature between the experiments and numerical results for TCPA [27].used for studying the driving response of the TCPA is shown in figure 12.The predicted driving forces align reasonably well with the experimental results.Our numerical model effectively captures the driving responses of the TCPA for various heating wire powers.Moreover, the driving forces increase as the heating powers are escalated.
Additionally, our proposed model can also be applied to predict the driving strain under specific temperature conditions, as shown in figure 13.We compare the model predictions with the corresponding experimental data.As can be seen from the figure, the driving strain of the TCPA escalates with increasing temperature, and our model predictions align with the experimental results reported by [26].

Effects of geometrical parameters on driving responses
After validating the accuracy the numerical model, we proceed to investigate the influence of various geometrical parameters on the driving forces.Initially, we analyze the effects of the radius of polymer fibers and the helical radius of the TCPA on the driving forces.Figure 14 illustrates a comparison of the predicted driving forces for different radius of polymer fibers.The driving forces are plotted as a function of the helical radius, with a fixed helical pitch h = 1.23 mm.As the helical radius increases, the driving forces increase until reaching a maximal value, and then decline with  Comparison of the variations of the driving strain with temperature between the experiments and numerical results for TCPA [39].
further increments in the helical radius.Evidently, the driving forces reach their maximum value when the helical radius is 0.8 mm.In figure 15, we depict the driving forces varying with the helical pitch for different fiber radii, while keeping the helical radius constant at R = 1.083 mm.It is worth mentioning that the driving forces attain their maximum value at the helical  pitch h = 2 mm.Subsequently, as the helical pitch increases, the driving forces begin to decrease.This conclusion is consistent with the experimental results [29].
Furthermore, figure 16 showcases the effects of the helical radius on the driving forces for different helical pitches with a fixed polymer fiber radius of r f = 0.571 mm.As the helical pitch increases, the driving forces escalate until reaching a maximum value, and then decline with further increments in the helical pitch.In addition, a larger helical radius indicates a higher corresponding optimal helical pitch.

Conclusions
In this paper, we have successfully introduced a numerical scheme to simulate the thermo-mechanical properties and responses of the TCPA based on the theoretical model proposed by Yang and Li [34].Our study has accomplished the following tasks: (1) We have proposed an efficient multi-layer model to solve the governing equation system of the model.Through this numerical method, we were able to determine the effective mechanical and thermal properties of the TCPA.(2) We have comprehensively simulated the thermo-mechanical responses of the TCPA under heating loads.The obtained numerical results exhibit a high degree of consistency with the corresponding experimental data.These results confirm accuracy of the current model and the efficiency of the numerical algorithm.(3) Based on the numerical simulation results, we have analyzed the optimal geometric features of the TCPA.These findings strongly support the improvement of the driving performance of the TCPA.
Compared to previous modeling works, our current model offers several advantages: We have considered the entire thermo-mechanical system of the TCPA, allowing us to predict all driving responses.We have successfully simulated the thermo-mechanical properties within the TCPA, enabling the proper optimization of certain geometrical features.We have proposed a parameterized numerical model based on Python code running in the Abaqus environment.The utilization of Python code has facilitated the execution of challenging tasks that would have been difficult or impossible to establish manually in Abaqus.
In summary, our present multi-layer model can be effectively applied to multi-stage helical materials and artificial muscles with twist and coiled structures.This study is highly significant in the development of techniques to simulate the driving responses of the TCPA under complex driving environments.

Figure 1 .
Figure 1.Schematic drawings of the concentric hollow cylinder.

Figure 2 .
Figure 2. The continuity and boundary conditions.

Figure 3 .
Figure 3. (a) Schematic drawings of the twisted and coiled polymer actuator.(b) The twisted fiber is decomposed into concentric hollow cylinders, and cut and separated into a lamina.(c) The corresponding on-axis lamina.(d) The corresponding off-axis lamina.

Figure 4 .
Figure 4. Multi-layer modeling framework of the TCPA.

Figure 5 .
Figure 5. Flowchart of the multi-layer modeling framework to solve the thermo-mechanical responses of the TCPA.

Figure 6 .r 4 f 2 is the polar moment of area, I = π r 4 f 4
Figure 6.(a) Effective Young's modulus, (b) Poisson's ratio and (c) thermal expansion coefficient as a function of the fiber radius under several representative layer numbers.

Figure 7 .
Figure 7. (a) Effective Young's modulus, (b) Poisson's ratio and (c) Thermal expansion coefficient as a function of the twisting angle under several representative layer numbers.(d) The driving displacement δ of a TCPA is directly obtained by applying Castigliano's theorem.

Figure 8 .
Figure 8. Driving force response of the TCPA with representative layer numbers.

Figure 9 .
Figure 9.Comparison of the variations of the driving strain with temperature between the experiments and numerical results for TCPA[27].

Figure 10 .
Figure 10.Comparison of the driving force between the experiments and multi-layer finite element results for different powers of heating wire.

Figure 11 .
Figure 11.Comparison of the driving force between the experiments and multi-layer finite element results for different heating times.

Table 1 . 6 Figure 12 .
Figure 12.(a) Photograph of a thermally driven test setup used for studying the driving response of the TCPA.(b) TCPA sample.

Figure 13 .
Figure 13.Comparison of the variations of the driving strain with temperature between the experiments and numerical results for TCPA[39].

Figure 14 .
Figure 14.Variations of the driving force with helical radius and radius of precursor fibers, and fixed helical pitch h = 1.23 mm.

Figure 15 .
Figure 15.Variations of the driving force with helical pitch and radius of precursor fibers, and fixed helical radius R = 1.083 mm.

Figure 16 .
Figure 16.Variations of the driving force with helical radius and helical pitch, and fixed the radius of precursor fibers r f = 0.571 mm.