NOMENCLATURE
- A
cross-sectional area of SMA wires
- As
austenite starting temperature
- Af
austenite finishing temperature
- B
bulk modulus
- CA
austenite stress influence coefficient
- CM
martensite stress influence coefficient
- E
Young’s modulus
- EA
Young’s modulus of austenite phase
- EM
Young’s modulus of martensite phase
- F
applied vertical force
- G
shear modulus
- I
second moment of area
- k
spring stiffness
- l
length of SMA wires
- L
length of cantilever beam
- Ms
martensite starting temperature
- Mf
martensite finishing temperature
- T
temperature
$\delta_{\text{tip}}$tip deflection
- $\varepsilon$
SMA strain
- $\varepsilon$L
SMA maximum recoverable strain
- $\nu$
Poisson’s ratio
- $\phi$
diameter of SMA wires
- $\Theta$
thermoelastic coefficient
- $\Omega $
phase transformation coefficient
- $\rho $
mass density
- $\sigma $
SMA stress
- $\xi $
martensite volume fraction
1.0 INTRODUCTION
Shape-memory alloys (SMAs) based on nickel-titanium (NiTi) have attracted enormous attention from researchers since their discovery in the 1960s by the Naval Ordnance Laboratory, who referred to such alloys as “Nitinol”.(Reference Buehler and Wiley1) This is due to their well-known and unique properties, namely the shape-memory effect (SME) and their pseudoelasticity (PE). These properties are related to the changes in the SMA crystal structure from an austenite phase to a martensite phase (forward transformation) and vice versa (reverse transformation) upon variation of the applied temperature. These SME properties have been utilised in many industrial applications, such as circular stents to enlarge blood vessels,(Reference Frotscher, Nortershauser, Somsen, Neuking, Bockmann and Eggeler2,Reference Petrini and Migliavacca3) heat engines(Reference Tobushi, Date and Miyamoto4,Reference Banks5) and robotic arms.(Reference Kim, Han, Song and Ahn6)
A great deal of research has been conducted over the past few decades to study the use of the SME properties of SMAs in morphing aerospace applications. In engineering, morphing refers to continuous shape changes where one entity deforms upon actuation without the movement of discrete parts relative to each other.(Reference Thill, Etches, Bond, Potter and Weaver7) Morphing is a short form of metamorphose, and shape morphing of an aircraft refers to geometrical changes, such as the adjustment of aerofoils.(Reference Barbarino, Bilgen, Ajaj, Friswell and Inman8) Using morphing technology, a smooth change in aerofoil shape can be achieved when compared with conventional hinged control surfaces.(Reference Iannucci, Evans, Irvine, Patoor and Osmont9)
The unique SME characteristic of SMAs has been exploited for the design of actuators to deform various aerospace structures, such as aircraft wings. Jardine et al.(Reference Jardine, Bartley-Cho and Flanagan10) utilised SMA torque tubes that were capable of producing approximately 3,500in-lb blocking torque to morph a military aircraft wing. This resulted in a 5° spanwise twist and an improvement in the aerodynamic performance by 8–12%. Researching a composite morphing aerofoil section, Dong et al.(Reference Dong, Boming and Jun11) exploited SMA springs to deform the skins between the leading and trailing edges, achieving a maximum vertical deflection of 5mm.
Some researchers have employed SMA wires to morph aerofoil sections and wings, yielding promising deformation results. Iannucci and Fontanazza(Reference Iannucci12) investigated a shape-memory wire concept, in addition to other smart morphing concepts, where NiTi SMA wires were embedded into the surface of a cambered wing. The expansion and contraction of the SMA wires were modelled based on equivalent temperature changes in the wires (positive and negative thermal expansion). A wing twist with approximately 5.27mm out-of-plane displacement was achieved upon SMA actuation. Kang et al.(Reference Kang, Kim, Jeong, Lee and Ahn13) and Rim et al.(Reference Rim, Kim, Kang and Lee14) morphed a balsa wing using resistive heating of six SMA wires, achieving a maximum flap deflection of 21° when supplying a current of 3.3A. Almeida et al.(Reference Almeida, Santos and Otubo15) exploited NiTi SMA wires to actuate a morphing wing rib made of polylactic acid (PLA) material, manufactured using a 3D printer. Their experimental results showed that the 1.05–1.06% strain recovery of appropriately located SMA wires could produce trailing-edge deflections of 11.6–14.0mm. Aerodynamically, the trailing-edge deflections resulted in an improvement of the maximum lift-to-drag ratio by 10–13% at small angles of attack when compared with a hinged-flap aerofoil.
Due to their structural flexibility and anisotropic properties, honeycomb and corrugated structures have been exploited by several researchers to morph aerofoils/wings. Iannucci(Reference Iannucci16) patented a camber morphing aerofoil concept where the aerofoil consisted of a honeycomb core structure and a polyurethane-reinforced glass/carbon fibre composite skin. SMAs have been proposed as smart materials for actuation mechanisms. Ruangjirakit(Reference Ruangjirakit17) introduced a composite morphing aerofoil section consisting of a lower corrugated carbon fibre composite skin. In experimental tests, a linear actuator was used to deform the aerofoil section, resulting in a deflection angle of approximately 4–6°. Thill et al.(Reference Thill, Etches, Bond, Potter and Weaver18) experimentally investigated a composite corrugated sandwich structure for the skin panels of a morphing trailing edge. The corrugated structure was manually actuated using two huge cross-bars with a torque wrench, achieving trailing-edge deflection angles of up to 12°. The corrugated outer surfaces of aerofoils/wings should not significantly affect the aerodynamic flow at low Reynolds number.(Reference Terwagne, Brojan and Reis19) To achieve a trade-off between aerofoil stiffness and deformability/flexibility, Kudva(Reference Kudva20) exploited two layers of flexible honeycomb core with a centre laminate, silicone skins and an aluminium tip for the trailing-edge section of an uninhabited combat air vehicle (UCAV). Experimentally, the trailing-edge section of the UCAV was actuated using a series of eccentuators, achieving a maximum deflection angle of 20° in wind-tunnel tests.
Numerically, several attempts have been made to model the SME of SMAs using finite-element methods, by implementing different constitutive models into FEA software. Gao et al.(Reference Gao, Qiao and Brinson21) claimed that Birman(Reference Birman22) applied the Tanaka model in his study of sandwich panels reinforced with NiTi SMA fibres, but that study focused more on an analysis of their buckling stability for different (sinusoidal and uniform) SMA distributions by applying previously reported experimental stress–temperature behaviour(Reference Cross, Kariotis and Stimler23). They developed a one-dimensional (1D) SMA model in a user subroutine for ABAQUS finite-element analysis software by implementing a Brinson constitutive model. The 1D model was tested for cases of isobaric, isothermal and constrained recovery, and subsequently applied to a hybrid composite with a self-healing function to close open cracks. More recent research on modelling the SMA for the actuation of structures includes a UMAT in ABAQUS to actuate smart landing gear for an airship,(Reference Alipour, Kadkhodaei and Ghaei24) a user element (UEL) subroutine in ABAQUS to actuate an adaptive strip,(Reference Solomou, Machairas and Saravanov25) a subroutine in ABAQUS to investigate the latent heat and performance of the SMA,(Reference Tabesh, Lester, Hartl and Lagoudas26) an exploitation of an existing superelastic SMA material model with certain boundary conditions to deform a lightweight actuator structure(Reference Okabe, Sugiyama and Inayoshi27) and SMA models based on the SMA thermal expansion to actuate a chevron, a cambered surface and a wing trailing edge.(Reference Iannucci12,Reference Turner, Cabell, Cano and Silcox28,Reference Barbarino, Pecora, Lecce, Concilio, Ameduri and Calvi29) Sellitto and Riccio(Reference Sellitto and Riccio30) adopted a simplified SMA material model based on the Coefficient of Thermal Expansion (CTE) in ABAQUS to actuate a hinged spoiler trailing edge and a rear upper panel of a vehicle. Sufficient target displacements of 10mm were achieved in both cases. Those authors also presented extensive reviews on SMAs and their applications in morphing structures, and forecast that a total of 70,000 scientific articles and 50,000 U.S. patents related to SMAs would be produced during the next 10 years, due to their promising potential.
However, most FEA research in this field has been implemented implicitly,(Reference Alipour, Kadkhodaei and Ghaei24,Reference Solomou, Machairas and Saravanov25,Reference Okabe, Sugiyama and Inayoshi27,Reference Turner, Cabell, Cano and Silcox28) which is the main distinction from the current presented research. Moreover, most of the actuated structures had simple forms such as a beam for landing gear(Reference Alipour, Kadkhodaei and Ghaei24) or an adaptive strip,(Reference Solomou, Machairas and Saravanov25) with the results being limited to displacement curves without any illustration of the actuated structures. In addition, some of this research has focused more on the superelastic properties of the SMA (without including temperature changes), with the SME properties only being investigated for a constant load case and without any examples of actuated structures.(Reference Tabesh, Lester, Hartl and Lagoudas26) Hardly any work has been carried out on modelling the SME of SMAs explicitly in commercial finite element analysis software, and explicitly for LS-DYNA, especially for applications to real and complex structures such as composite morphing aerofoils.
This paper presents the development of a user-defined material (UMAT) model for the SME of an SMA wire. The SMA wire was modelled using a beam element with one integration point, mimicking a truss element. The Tanaka one-dimensional SMA constitutive model(Reference Tanaka, Hayashi and Itoh31) was implemented within an explicit LS-DYNA UMAT subroutine written in FORTRAN to simulate the thermomechanical behaviour of the SMA wire. For the case of an SMA wire connected in series with a linear spring, a Newton–Raphson iteration method was applied to evaluate the mechanical responses for each change in the applied temperature. The stiffness of the spring could represent the stiffness of the actuated geometrical structure. The simulation results were compared with analytical results obtained using MATLAB. Another comparison of the UMAT model with an analytical solution was then conducted for an SMA-actuated cantilever beam. The UMAT was then applied in a mesh sensitivity study, a key parameters study, and a constant load case. Finally, the UMAT was extended to multiple SMA wires for the actuation of aluminium and composite aerofoil sections.
The originalities of this research are the new UMAT for the SMA in the explicit LS-DYNA software, and the new concept of SMA-actuated aerofoils based on the exploitation of corrugated surfaces. Even though there are various studies on 1D SMAs, hardly any research has focussed on an explicit formulation for the SMA, especially in the explicit LS-DYNA software. Robust verification of this is found in the fact that no explicit SMA material models have been successfully developed for LS-DYNA, particularly including SME properties. The new, explicit SMA material model developed in this paper offers several advantages, including a significant saving in simulation time due to the avoidance of matrix inversion steps, and the great potential for application of the UMAT of the SMA in dynamic (fast actuation) and impact (on SMA-actuated morphing aerostructures) simulations in future research. Regarding the concept of morphing aerofoils, although much research have utilised corrugated surfaces as outer skins of aerofoils, none of them have utilised SMA actuators between the corrugated cells. The morphing concept presented in this paper, where the deformation of the aerofoils can be smoothly controlled by the SMA wires between the corrugated cells, is therefore considered to be original.
This work is of considerable importance in the field of morphing structures, as the UMAT developed for the SMA can serve as a design tool for the virtual design of realistic composite morphing aerofoils/wings. Moreover, the explicit LS-DYNA software is capable of multiphysics-based incompressible fluid dynamics (ICFD) simulations, which will be used in future research to simulate and analyse the movement of actuated structures under realistic airflow conditions with the inclusion of Fluid–Structure Interaction (FSI).
2.0 FINITE-ELEMENT MODELLING OF SMA WIRES USING A USER-DEFINED MATERIAL (UMAT) MODEL IN EXPLICIT LS-DYNA
The commercial explicit finite element analysis (FEA) software LS-DYNA does not include a material model that is capable of simulating the SME behaviour of SMAs. The most relevant material model is MAT_030 i.e. MAT_SHAPE_MEMORY.(32) However, this material model is ideal to capture the superelastic behaviour of SMAs at a defined, constant temperature above the austenite finishing temperature, A f. It is not suitable for actuation processes of an SMA-actuated structure, where temperature changes take place to activate the SME. Therefore, in this research, a new user-defined material model for SMA wires was developed by implementing a one-dimensional SMA constitutive model in FORTRAN.
The Tanaka SMA constitutive model was implemented.(Reference Tanaka, Hayashi and Itoh31)
\begin{equation}\sigma - {\sigma _o} = E\left( \xi \right)\left( {\varepsilon - {\varepsilon _o}} \right) + \Theta \left( {T - {T_o}} \right) + \Omega \left( \xi \right)\left( {\xi - {\xi _o}} \right)\!,\end{equation}
where σ, ɛ, T and 
${\rm{\xi }}$
are the stress, strain, temperature and martensite volume fraction of the SMA, respectively. The subscript ‘o’ denotes the initial state of each variable. E, Θ and 
$\Omega$
are the Young’s modulus, the thermoelastic coefficient and the phase transformation coefficient of the SMA, respectively.
In this paper, the thermal expansion of the SMA is neglected, which is a common assumption when modelling SMAs because of their much smaller strain (10−5) compared with the SMA recovery strain (10−2). The Young’s modulus and the phase transformation coefficient are functions of the martensite volume fraction.
\begin{equation}E\left( \xi \right) = {E_A} - \xi \left( {{E_A} - {E_M}} \right)\end{equation}
\begin{equation}\Omega \left( \xi \right) = - {\varepsilon _L}E\left( \xi \right)\!,\end{equation}
where E A and E M are the Young’s modulus of the SMA in the austenite phase and martensite phase, respectively, while ɛL is the SMA maximum recoverable strain. The SMA constitutive Equation (1) therefore simplifies to
\begin{equation}\sigma - {\sigma _o} = E\left( \xi \right)\left( {\varepsilon - {\varepsilon _o}} \right) - {\varepsilon _L}E\left( \xi \right)\left( {\xi - {\xi _o}} \right)\end{equation}
The martensite volume fractions for a forward martensitic transformation and a reverse transformation have an exponential form.(Reference Tanaka, Hayashi and Itoh31) They can be obtained from the transformation kinetics derived from the 1st law (the conservation of energy principle) and the 2nd law of thermodynamics. For a forward martensitic transformation (cooling stage), the martensite volume fraction is defined as
\begin{equation}{\xi _{A \to M}}=1 - {e^{{a_M}\left( {{M_s} - T} \right) + {b_M}\sigma }}\end{equation}
For a reverse transformation (heating stage), it is defined as
\begin{equation}{\xi _{M \to A}} = {e^{{a_A}\left( {{A_s} - T} \right) + {b_A}\sigma }}\end{equation}
where A s and M s are the austenite and martensite starting temperatures, respectively. a A, b A, a M and b M are material constants defined in Equations (7) and (8), where A f and M f are the austenite and martensite finishing temperatures, respectively. C A and C M are the stress influence coefficients for the austenite and martensite phase, respectively.
\begin{equation}{a_A} = \frac{{\ln \left( {0.01} \right)}}{{{A_s} - {A_f}}},\;{b_A} = \frac{{\ln \left( {0.01} \right)}}{{{C_A}\left( {{A_s} - {A_f}} \right)}} = \frac{{{a_A}}}{{{C_A}}}\;\end{equation}
\begin{equation}{a_M} = \frac{{\ln \left( {0.01} \right)}}{{{M_s} - {M_f}}},\;{b_M} = \frac{{\ln \left( {0.01} \right)}}{{{C_M}\left( {{M_s} - {M_f}} \right)}} = \frac{{{a_M}}}{{{C_M}}}\end{equation}
The simplified SMA constitutive Equation (4) was solved for the SMA stress and strain by implementing the Newton–Raphson iteration method for a complete heating–cooling cycle. A stress function was obtained by setting Equation (4) to zero. For the reverse transformation (heating stage), with zero initial stress and strain and initial martensite volume fraction of unity (fully martensite, 
${{\rm{\xi }}_{\rm{o}}}$
= 1), the stress function and its derivative with respect to stress are
\begin{equation}f\left( \sigma \right) = \sigma - E\left( \xi \right)\varepsilon + {\varepsilon _L}E\left( \xi \right)\left( {\xi - 1} \right)\end{equation}
\begin{equation}\frac{{\partial f\left( \sigma \right)}}{{\partial \sigma }}=1 - \frac{{\partial E\left( \xi \right)}}{{\partial \sigma }}\varepsilon + {\varepsilon _L}\left( {E\left( \xi \right)\frac{{\partial \xi }}{{\partial \sigma }} + \xi \frac{{\partial E\left( \xi \right)}}{{\partial \sigma }} - \frac{{\partial E\left( \xi \right)}}{{\partial \sigma }}} \right)\end{equation}
For the martensitic transformation (cooling stage), with zero initial martensite volume fraction (fully austenite, 
${{\rm{\xi }}_{\rm{o}}}$
= 0), the stress function and its derivative with respect to stress are
\begin{equation}f\left( \sigma \right) = \sigma - {\sigma _o} - E\left( \xi \right)\left( {\varepsilon - {\varepsilon _o}} \right) + {\varepsilon _L}E\left( \xi \right)\xi \end{equation}
\begin{equation}\frac{{\partial f\left( \sigma \right)}}{{\partial \sigma }}=1 - \left( {\varepsilon - {\varepsilon _o}} \right)\frac{{\partial E\left( \xi \right)}}{{\partial \sigma }} + {\varepsilon _L}\left( {E\left( \xi \right)\frac{{\partial \xi }}{{\partial \sigma }} + \xi \frac{{\partial E\left( \xi \right)}}{{\partial \sigma }}} \right)\end{equation}
For each increment of temperature on the heating stage, or each decrement of temperature on the cooling stage, the current stress function and its differentiation were used with the current state of stress of the SMA wire to calculate the next state of stress of the SMA wire, by using Equation (13). The iteration was repeated until the equation converged, or simply the left-hand side (LHS) was equal to the right-hand side (RHS).
\begin{equation}{\sigma _{n+1}} = {\sigma _n} - \frac{{f\left( {{\sigma _n}} \right)}}{{f'\left( {{\sigma _n}} \right)}}\end{equation}
The Newton–Raphson iteration method was applied because there are two main unknowns that must be solved after each temperature change, viz. the martensite volume fraction in Equation (5) or (6) and the SMA stress in Equation (4) (including the stress functions and their derivatives in Equations (9-12)). The martensite volume fraction (Equation (5) or (6)) was a function of the SMA stress, and the SMA stress (Equation (4)) was a function of the martensite volume fraction. Neither of them could be solved directly by substituting the value of the temperature, hence the Newton–Raphson iteration method was required to solve for the SMA stress.
Here, the SMA strain ɛ was obtained from LS-DYNA for each time step and treated as a constant during the differentiation from Equation (9) to (10) and from Equation (11) to (12). The results for the SMA stress and strain were compared with an analytical solution obtained using MATLAB, where the SMA strain was expressed as a function of the SMA stress, which had to be taken into account in the differentiation of the stress functions.
\begin{equation}\varepsilon = - \frac{{\sigma A}}{{kl}},\;\;\frac{{\partial \varepsilon }}{{\partial \sigma }} = - \frac{A}{{kl}},\end{equation}
where k, l and A are the spring stiffness, and the SMA wire length and cross-sectional area, respectively.
The main differences between the user-defined material (UMAT) model and the analytical model are the derivatives in Equations (10) and (12), where the UMAT does not require the spring stiffness k and the SMA wire length, l. This demonstrates that the newly developed and generalised UMAT can be used for any SMA wire having arbitrary length and connected to any actuated structure having any stiffness. In other words, the LS-DYNA user does not have to work out and define the stiffness of the actuated structure or the length of the SMA wires in the keyword input stage, which simplifies the modelling stage.
In the explicit LS-DYNA software, nodal forces associated with mass/inertia and damping were included. Another requirement for the explicit analysis was that the time step must be less than the Courant time step (the time for a sound wave to travel across an element). These were some of the reasons why the explicit solver was preferred for dynamic analysis. Furthermore, explicit analysis did not require the numerical solver to invert the stiffness matrix, which is an expensive operation, especially for large models. The UMAT developed in the explicit LS-DYNA software requires the definition of the shear modulus and bulk modulus of the SMA material, which change during its phase transformation. Because the SMA Young’s modulus changes during the martensite-to-austenite phase transformation (heating stage) and the austenite-to-martensite phase transformation (cooling stage), the shear modulus and the bulk modulus of the SMA also changed, as they are functions of the SMA Young’s modulus. Their relations with the SMA Young’s modulus are shown in Equations (15) and (16), respectively.
\begin{equation}G\left( \xi \right) = \frac{{E\left( \xi \right)}}{{2\left( {1 + \nu } \right)}}\end{equation}
\begin{equation}B\left( \xi \right) = \frac{{E\left( \xi \right)}}{{3\left( {1 - 2\nu } \right)}}\end{equation}
Specific values of the material constants were assigned to the bulk modulus and the shear modulus. In the Newton–Raphson iteration, they were assigned to the respective historical values of the variables. After solving at each time step, the historical values of the variables were passed to the material constants. This was another important difference between the programmed UMAT in the explicit LS-DYNA and the MATLAB analytical solution used for comparison. Neither the shear nor the bulk modulus was required in the MATLAB code. This is also one of the crucial differences between the current development of the UMAT explicitly in LS-DYNA and previous research on UMATs implicitly in ABAQUS.
The SMA material properties that were used in Equations (1-16) are listed in Table 1 with other additional parameters. For a clearer illustration, the structure of the UMAT is shown in Fig. 1. The SMA wire was modelled as a beam element with one integration point, mimicking a truss element. The beam element was used instead of a truss element because the number of integration points must be defined in a ‘database extend binary’ keyword file, otherwise no results for the historical values of the variables would be available in the output file. This integration point could not be assigned to a truss (or cable) element that used a resultant formulation. A circular cross section was assigned to the beam, which was appropriate for a wire.
At a given time step, the beam element of the SMA processed the applied temperature as an input, and output the SMA stress and strain. Due to the change in the SMA strain, the nodes connected to the ends of the SMA wire were subjected to acceleration (unless the node was subjected to a certain boundary condition e.g. a clamped end of the SMA wire as shown in Fig. 3). Then the velocity and the displacement of the nodes were computed. From the displacement, the strain of the actuated element (i.e. a shell element) to which the nodes were connected was calculated. Finally, from the strain of the element, the stress in the element was solved. Consequently, all neighbouring elements of the actuated structure were accelerated, displaced and stressed. The cycle was repeated for the next time step, with a new incremental temperature input for the SMA wire. This process is illustrated in the higher-level flow chart shown in Fig. 2.
Table 1 NiTi SMA material properties and additional parameters
Figure 1. UMAT structure.
Figure 2. UMAT of SMA for structure actuation in explicit LS-DYNA.
Figure 3. SMA wire connected in series with a linear spring.
The results of the SMA material model were first compared with the MATLAB analytical solution for two cases: an SMA wire connected in series with a linear spring, and an SMA-actuated cantilever beam. For the first case, the spring stiffness may represent the stiffness of actual geometrical structures. For the SMA-actuated cantilever beam, the LS-DYNA results were compared with the analytical solution of the SMA–spring model, with the spring stiffness defined by Equation (17).
\begin{equation}k = \frac{F}{{{\delta _{tip}}}} = \frac{{3EI}}{{{L^3}}}\end{equation}
The variables F, 
$\delta$
tip, E, I and L are the vertical force acting on the tip, the beam tip deflection, the Young’s modulus, the second moment of area and the length of the cantilever beam, respectively.
Once the comparisons were completed, several design parameters such as the spring stiffness, and the SMA wire length and cross-sectional area were varied to study their effects on the SMA thermomechanical behaviour. A mesh sensitivity study and a constant load case were also analysed. The SMA material model was finally utilised as a design tool to design and actuate aluminium and composite morphing aerofoils.
3.0 COMPARISON BETWEEN USER-DEFINED MATERIAL (UMAT) MODEL OF THE SMA AND MATLAB ANALYTICAL SOLUTION
The SMA material model was initially compared with analytical solutions for two cases: (1) an SMA wire connected in series with a linear spring, and (2) an SMA-actuated cantilever beam.
3.1. An SMA wire–linear spring model
The results of the SMA material model were first compared with the MATLAB analytical solution for the case of an SMA wire connected in series with a linear spring (Fig. 3). The SMA thermomechanical behaviour predicted by the LS-DYNA simulation was compared with the SMA stress and strain given by the analytical solution in MATLAB. The differences between the numerical and analytical models were described in the previous section 2.
The SMA wire was modelled as a beam element having one integration point, mimicking a truss element, while the linear spring was modelled using a discrete element. The diameter and length of the SMA wire were 0.5 and 100mm, respectively. Other SMA material properties were listed in Table 1. The stiffness of the linear spring was 3.5N/mm. One complete heating–cooling cycle was applied on the SMA wire, and the SMA stress and strain were evaluated.
The results for the SMA stress and strain given by both the UMAT model and the MATLAB analytical solution are depicted in Fig. 4. The red and blue lines represent the SMA stress and strain during the heating and cooling stage, respectively.
Figure 4. Comparison of the user-defined SMA material model (UMAT) with the MATLAB analytical solution: SMA stress and strain for an SMA wire connected in series with a linear spring.
The SMA wire was first heated from 25°C to 85°C. As the temperature exceeded the austenite starting temperature, A s, the wire contracted as it recovered its initial memorised length. Consequently, the wire extended the linear spring, and hence the SMA stress increased until the martensite-to-austenite transformation completed. On cooling, as the temperature was decreased to below the martensite starting temperature, M s (added to the stress-influenced part), a martensitic transformation took place where the SMA crystal structure changed from the austenite phase to a martensite phase. During this forward transformation, the martensite volume fraction, 
$\xi$
, increased from zero to one, and the elastic modulus of the SMA wire reduced from 75 to 40GPa. As a result, the energy stored in the linear spring from the previous reverse transformation was converted to extend the SMA wire.
It can be observed from Fig. 4 that the austenite finishing temperature, where the SMA stress and strain started to remain constant at the end of the heating stage, was higher than 65°C (Table 1). This was because of the effect of the stress influence coefficient on the transformation temperature. A similar explanation applied to the higher martensite starting temperature at the beginning of the cooling stage.
The LS-DYNA simulation results (solid lines) were compared with the analytical solutions (dotted lines), showing extremely good agreement. At the end of the heating stage, the final SMA stress of the numerical simulation was slightly lower than the analytical solution.
3.2. An SMA-actuated cantilever beam
Secondly, a comparison between the FE simulation (the UMAT of the SMA) and the MATLAB analytical solution for the SMA-actuated cantilever beam (Fig. 5) is presented. The SMA wire was modelled with similar dimensions and properties as described in section 3.1. The cantilever beam was modelled using an elastic aluminium material having a density of 2.81 × 10−6kg/mm3, a Young’s modulus of 71.7GPa, and a Poisson’s ratio of 0.33. The cantilever beam dimensions were 225 × 1 × 10mm3, and it was meshed with 90 3D shell elements. One end of the beam was fixed in all degrees of freedom, while the other end was actuated by an SMA wire in the vertical direction. A complete heating–cooling cycle was applied on the SMA wire. The SMA stress and strain, as well as the beam tip deflection, were evaluated and compared.
Figure 5. SMA-actuated cantilever beam.
The SMA stress and strain are shown in Fig. 6. When the SMA wire was heated and the temperature exceeded the austenite starting temperature, the SMA stress increased and the SMA strain decreased (the SMA wire contracted), as shown by the red lines. The SMA stress and strain reached constant magnitudes, which were 12.6 MPa and 1.58%, respectively, when the martensite-to-austenite phase transformation completed. When cooled, as the temperature dropped to below the martensite starting temperature (added to the stress-influenced part), the SMA stress decreased while the SMA strain increased, as shown by the blue lines. They reached constant magnitudes when the martensitic transformation completed.
Figure 6. SMA stress and strain of the SMA-actuated cantilever beam.
The LS-DYNA simulation results, shown by the solid lines, were compared with the analytical results (shown by the dotted lines) with a spring stiffness of k = 3EI/L 3 (refer to Equation (17)). The tip deflection of the cantilever beam is depicted as a function of the SMA temperature in Fig. 7 in comparison with the analytical result. The numerical and analytical maximum tip deflections were −1.57mm and −1.58mm, respectively, showing a 0.6% difference. The tip deflections correspond to the SMA maximum recoverable strain of 1.6%.
Figure 7. Tip deflection of the SMA-actuated cantilever beam.
4.0 FEA SIMULATION RESULTS: MESH SENSITIVITY AND KEY PARAMETER STUDIES, AND A CONSTANT LOAD CASE
4.1. Mesh sensitivity study
A mesh sensitivity study was conducted on the SMA material model. The SMA wire–linear spring model with the properties and dimensions described in section 3.1 was used in this study. The length of the SMA wire was kept constant at 100mm, while the number of beam elements along the wire length was increased from 1 to 2, 5, 10 and 100. One complete heating–cooling cycle was applied on the SMA wire, and the resulting SMA stress and strain are depicted in Fig. 8. Excellent agreement of the SMA thermomechanical behaviour was obtained regardless of the number of beam elements used. This finding indicated that one element was adequate to accurately predict the SMA behaviour. This study was quite important, as it was found that, for the UMAT of the SMA in the explicit LS-DYNA software, a total strain formulation should be used instead of an incremental strain formulation for the SMA.
Figure 8. Mesh sensitivity: the SMA stress and strain.
4.2. Key parameter study
The influence of several key parameters on the SMA thermomechanical behaviour was investigated. The parameters were the stiffness of the linear spring (k), and the length (l) and the cross-sectional area (A) of the SMA wire. The SMA wire–linear spring model was simulated with the range of k = 1, 10 and 100N/mm, l = 100 and 250 mm and diameter, D = 0.1, 0.2, 0.4 and 0.5mm. The combinations of these parameters resulted in normalised stiffness (
$\pi$
kl/4 A or kl/D 2) values of 6,250, 10,000, 25,000, 62,500, 100,000 and 250,000N/mm2. Two combinations of the parameters were simulated for each normalised stiffness. The SMA stress and strain were evaluated for one complete heating–cooling cycle.
The SMA stress and strain results obtained for the above-mentioned normalised stiffness (
$\pi$
kl/4 A) values are depicted in Fig. 9. It could be observed that the recovery stress, the residual stress and the transformation temperatures (A f, M s and M f) increased with an increase in the spring stiffness or the SMA length, and with a decrease in the SMA cross-sectional area. The recovery strain decreased in magnitude, but the residual strain first increased and then decreased, as the normalised stiffness was increased. The results for the SMA residual stress and strain need further fundamental understanding of the SMA constitutive model. A similar tendency in the residual stress and strain that remained at the end of the cooling phase was also observed in the actuation of helicopter rotor blades, experimentally and numerically, even though the Brinson SMA constitutive model was applied.(Reference Epps and Chopra33)
Figure 9. SMA stress and strain for several normalised stiffness (
$\pi$
kl/4 A) values.
The changes in the SMA stress, SMA strain and transformation temperatures were greatly influenced by the changes of the spring stiffness and the SMA cross-sectional area. These effects could be observed if the SMA stress and strain were compared separately for each value of the spring stiffness and the length or cross-sectional area of the SMA wires.
The results in Fig. 9 showed that, to achieve a high actuation displacement or a high recovery strain, the normalised stiffness (
$\pi$
kl/4 A) should be minimised. This could be achieved by minimising the stiffness of the actuated structure or maximising the SMA cross-sectional area.
4.3. A constant load case
The SMA material model was then tested for a constant load case. The properties and dimensions of the SMA wire were similar to the wire described in section 3.1. One end of the SMA wire was constrained in all degrees of freedom, while the other end was subjected to an axial tensile load. The load resulted in an SMA stress of approximately 10MPa. One complete heating–cooling cycle was applied. Because the SMA stress was constant, only the SMA strain was obtained from the finite element simulation, as depicted in Fig. 10. The initial SMA strain was positive before the heating cycle, because of the elongation of the SMA wire due to the applied tensile load. When heated, the SMA strain decreased non-linearly to almost 1.6% as the SMA wire contracted. When cooled, it increased back to the initial strain.
Figure 10. SMA strain for a constant load case.
5.0 FEA SIMULATION RESULTS: MORPHING AEROFOIL SECTIONS
5.1. Aluminium morphing aerofoils
Having compared the results of the SMA material model with the analytical solution and testing it in various cases, it was then applied to actuate aluminium morphing aerofoils. To achieve the morphing objective of a sufficiently high trailing-edge deflection, aerofoil sections consisting of a rigid D-nose spar, a corrugated section/surface and a trailing-edge section were investigated. The aerofoils had a chord length of 300mm, a width/span of 10mm, a rigid D-nose spar at 25% (75mm) in front of the chord and a flexible section at 75% (225mm) aft of the chord. As depicted in Fig. 11, two different positions of the corrugated section/surface were considered, firstly as a lower corrugated skin (Fig. 11a) and secondly as a flexible corrugated cantilever beam (Fig. 11b). The corrugated section/surface consisted of 12 cells, with a side length of approximately 5mm. The aerofoil sections would be more practical if a highly flexible structure, such as low-density flexible foam or a flexible honeycomb structure, filled the empty space between the D-nose spar and the trailing-edge section. One of the possible configurations was two layers of honeycomb core sandwiched between the corrugated cantilever beam and the silicone skins, such as the morphing method implemented by Kudva et al.(Reference Kudva20) However, this was omitted from the current analysis because of the assumption of low stiffness in the chordwise direction, which could be neglected.
Figure 11. Morphing corrugated aerofoil sections.
The NACA 0012 aerofoil sections were modelled with an elastic aluminium (7075) material, having a density of 2.81 × 10−6kg/mm3, an elastic modulus of 71.7GPa, and a Poisson’s ratio of 0.33. They were modelled with 3D shell elements, having thickness of 1mm for the D-nose spars and 0.5mm for the upper skin, corrugated sections/surfaces and trailing-edge sections. The width of the aerofoil sections was 10mm, with ten shell elements across the width. The SMA wires were modelled with beam elements having one integration point, with a diameter of 0.5mm and a maximum recoverable strain of 1.6%. One SMA wire was positioned in each cell of the corrugated surface, hence a total of 12 SMA wires were used for both cases.
The D-nose spar was fixed in all degrees of freedom, while a symmetrical boundary condition was applied to the side edges of the skin, the corrugated surface and the trailing-edge section. One complete heating–cooling cycle was applied on the SMA wires, and the deformations of the aluminium morphing aerofoils were evaluated.
The deformations of the aluminium aerofoils after SMA actuation are depicted in Fig. 12. For the aerofoil section with a corrugated lower skin, the resulted trailing-edge deflection was 7.8mm. The deflection of the trailing edge was improved to 65.9mm (eight times greater) when a corrugated section was positioned as a flexible cantilever beam in between the D-nose spar and the trailing-edge section. This larger deflection was expected, as there was no upper skin to provide additional stiffness to the aerofoil section in the chordwise direction. If two layers of honeycomb core with silicone skins were incorporated between the D-nose spar and the trailing-edge section, the trailing-edge deflection could possibly be reduced to one-third or one-half of the current value, depending on the type of core and the skin stiffness. This would require another study with greater resources in term of computational time and capability.
Figure 12. Trailing-edge deflections of aluminium morphing aerofoil sections.
5.2. Composite morphing aerofoils
Finally, the SMA material model was applied to actuate composite morphing aerofoils, rather than the isotropic aluminium aerofoils. A second aerofoil configuration consisting of a corrugated cantilever section/surface in between the D-nose spar and the trailing-edge section was used in the finite-element simulations.
The composite aerofoils had material properties of CFRP T300/914. The density was 1.53 × 10−6kg/mm3, the Young’s moduli were 129GPa (E11) and 8.4GPa (E22), the shear modulus was 4.2GPa and the Poisson’s ratio was 0.34 (v 12). The corrugated section had two different lay-ups: (a) a symmetrical lay-up [90°/45°/−45°/90°]s and (b) an anti-symmetrical lay-up (90/45/−45/90/90/45/−45/90)°, where 90° was the spanwise direction. The D-nose and trailing-edge section had the anti-symmetrical lay-up. The thickness of all parts was 1mm. A complete heating–cooling cycle was applied on the SMA wires, and the deformations of the composite aerofoils were analysed.
As depicted in Fig. 13, the trailing-edge deflections of the composite aerofoils with a corrugated section/surface having a symmetrical lay-up and an anti-symmetrical lay-up were 52.5mm and 52.0mm, respectively. As expected, there was a small amount of twist for the composite aerofoil with a corrugated section/surface having a symmetrical layup, because the components of the D matrix, which were D13 and D23, were not zero. The trailing-edge deflections were smaller than the aluminium aerofoil because of the higher stiffness of the composite aerofoils, due to the greater thickness and modulus.
Figure 13. Trailing-edge deflections of composite morphing aerofoil sections: corrugated section with (a) a symmetrical lay-up and (b) an anti-symmetrical lay-up.
6.0 CONCLUSIONS
A new UMAT for NiTi SMA wires was successfully developed and implemented into the commercial explicit LS-DYNA code. The results of the SMA material model were initially compared with MATLAB analytical solutions for the actuation of an SMA wire–linear spring structure and an SMA-actuated cantilever beam. The model was then tested to determine its mesh sensitivity and the effect of key parameters, and applied in a constant load case. It was further applied to actuate aluminium and composite aerofoil sections consisting of corrugated surfaces. Trailing-edge deflections of up to 65.9mm and approximately 52mm were achieved, respectively. The novelties of this research, including the new explicit formulation for the SMA wires and the new concept of morphing aerofoils with SMA wires arranged between the corrugated cells, contributed to the achievement of these results.
The composite corrugated section is currently manufactured using 3D printing technology, and the prototype will be tested and compared using another series of finite-element simulations to validate the SMA-actuated morphing aerofoils. Because of the low fibre volume fraction of the 3D-printed composite, only bench-top tests will be carried out. That study will be presented in another research paper. In future research, the SMA material model could be improved or extended to 2D and 3D forms of SMAs (i.e. SMA strips and tubes). The SMA model may also be applied for the actuation of a composite morphing wing under flight aerodynamic pressure loads, which can be represented experimentally by using a series of weights located along the wing span. The research may also be extended to FSI simulations of the composite morphing wing in the explicit LS-DYNA.
ACKNOWLEDGEMENTS
The research is funded by the Ministry of Education Malaysia and IIUM, which is hereby gratefully acknowledged.
