A Multiscale Analysis on the Superelasticity Behavior of Architected Shape Memory Alloy Materials

In this paper, the superelasticity effects of architected shape memory alloys (SMAs) are focused on by using a multiscale approach. Firstly, a parametric analysis at the cellular level with a series of representative volume elements (RVEs) is carried out to predict the relations between the void fraction, the total stiffness, the hysteresis effect and the mass of the SMAs. The superelasticity effects of the architected SMAs are modeled by the thermomechanical constitutive model proposed by Chemisky et al. 2011. Secondly, the structural responses of the architected SMAs are studied by the multilevel finite element method (FE2), which uses the effective constitutive behavior of the RVE to represent the behavior of the macroscopic structure. This approach can truly couple the responses of both the RVE level and structural level by the real-time information interactions between two levels. Through a three point bending test, it is observed that the structure inherits the strong nonlinear responses—both the hysteresis effect and the superelasticity—of the architected SMAs at the cellular level. Furthermore, the influence of the void fraction at the RVE level to the materials’ structural responses can be more specifically and directly described, instead of using an RVE to predict at the microscopic level. Thus, this work could be referred to for optimizing the stiffness, the hysteresis effect and the mass of architected SMA structures and extended for possible advanced applications.


Introduction
Cellular materials are widely used for their high strength-to-weight ratio and high energy absorption performance (Gibson and Ashby [1]; Ashby et al. [2]). For instance, honeycomb, folded cellular materials and foam are usually used as the core of the sandwich structures for dissipating the kinetic energy, damping or reducing the weight of the structure (Ashby et al. [2]; Yazdani et al. [3]; Garcia-Moreno [4]; Hangai et al. [5]; Strano et al. [6]). However, honeycomb and folded cellular materials have high manufacturing costs and moisture problem, as well as buckling problems (Rashed et al. [7]). The mechanical behaviors of foams are too difficult to be accurately measured for their stochastic cells, which always results in the excessive use of materials to satisfy the safety factor. To overcome these shortcomings, partially-ordered foams allowing limited structural control of the pore and spatial distribution of pore levels, such as metal syntactic foams (Taherishargh et al. [8][9][10]; by the effective behaviors of homogenized RVEs, and the strain states of the RVEs are given by the associated integration points. The FE 2 method shows good performance on multiscale modeling of SMA-based materials. Xu et al. [39] proposed a 3D FE 2 model for simulating composites with stiff SMA fibers embedded in a soft matrix. This model was validated by the literature and showed a good ability in modeling the superelasticity and the shape memory effects of SMA composites. Following previous works, it is therefore necessary and feasible to investigate the behavior of architected SMAs with a multilevel finite element model. Towards a better understanding of the studied architected SMA, unit cells (RVEs) with different void fractions are introduced to study the superelasticity effect of the materials and structures. In Section 2, a constitutive model for SMA is introduced briefly. Then, a parametric analysis with different void fractions at the cellular level is performed to predict the relations between the void fraction, the total stiffness, the hysteresis effect and the mass of the architected cellular SMAs. In Section 3, after a short presentation for the multilevel finite element model, multiscale modelings on architected SMA structures are carried out to simulate the structural and the cellular responses simultaneously.

Cellular Response
To investigate the architected SMA structures, RVEs with different geometrical parameters are studied firstly, then the structural response is investigated further by the computational homogenization method. First of all, a short presentation of the SMA constitutive model is given in the next subsection.

SMA Constitutive Model
The SMA model, proposed by Chemisky et al. [40], is adopted for the thermomechanical behavior modeling of the architected SMA. This model follows the work of Peultier et al. [41], who proposed a macroscopic phenomenological SMA approach based on the Gibbs free energy. This model was implemented on ABAQUS via user-defined materials (UMAT ) and validated by experiments. It was later improved by Chemisky and Duval; see Chemisky et al. [40]; Duval et al. [42]. Here, only a brief introduction is presented for the readers due to the limited length of the paper. The SMA constitutive model used here is able to describe four different strain mechanisms: the elastic strain ε e , the thermal expansion strain ε th , the martensitic transformation strain ε tr and the twin accommodation strain ε tw . To this end, the total strain is decomposed in the following form: where the above strains are formulated as: where S denotes the isotropic fourth order compliance tensor, α represents the isotropic thermal expansion tensor, T re f gives the temperature without expansion strain, f gives the martensitic volume fraction in SMA, the transformation strainε tr describes the mean strain over the martensite related to martensite reorientation and the mean strainε tw and self-accommodated martensitic volume fraction f FA are introduced to describe the twin accommodation over the martensite. The potential energy of the SMA model is based on the following Gibbs free energy expression: where U A and U M are the austenitic and the martensitic internal energy and S A and S M are the austenitic and the martensitic entropy. Several terms, such as ∆U = U M − U A , ∆S = S M − S A and ∆T = T − T re f , are also introduced. Considering T 0 = ∆U ∆S the equilibrium temperature of transformation, a linear variation of entropy around T 0 is defined as B = −∆S. The terms H tr , H f and H tw comprise, respectively, a set of material parameters characterizing interactions between grains of the microscopic RVE, between variants inside grains and between twins. C v is the transformation latent heat coefficient. The control equation of the thermodynamic system is established by introducing the Clausius-Duhem inequality: Substituting Equation (3) to Equation (4) and considering the thermo-elastic balance conditions, the Clausius-Duhem inequality is reduced to: Here, driving forces related to this dissipation expression are introduced: 1. Transformation driving force related to f : where the definition ofḟ FA = ζ FAḟ is introduced here.

2.
Orientation force related toε tr : where σ denotes the deviatoric part of the stress tensor σ.

3.
Twin accommodation related toε tw : By introducing a series of predefined critical thresholds related to these driving forces, the evolution of phase transformation is controlled. For more clarity of the relations between the main control equations, the constitutive equations for this model are summarized in Table 1.
The implementation of this model in a finite element code is realized; for more details about the implementation process, see Chemisky et al. [40] and Duval et al. [42]. Moreover, the evolution of tension-compression asymmetry and internal loops during the partial loadings can also be well simulated with this adopted model. In the following subsection, architected SMA RVEs with different geometric parameters are studied to see the cellular response of the material.

Convergence Analysis for the RVE Mesh
An architected SMA RVE is introduced to describe the microscopic structure; see Figure 1. It is formed by excavating cylindrical holes through the center of each face of a cube SMA. The size of the cube is given by 1 mm × 1 mm × 1 mm, while the radius of the cylindrical holes is given as 0.38 mm. According to computational homogenization theory, periodic boundary conditions are introduced into the RVE by the multi-point constraints (MPCs) on ABAQUS: where u is the displacement vector, x is the coordinates of nodes andε is the strain load applied on the RVE. The notations + and − denote the nodes on opposite boundaries; the notation ∆ represents the incremental case. Thus, a set of reference points (RPs) for three different directions is introduced in MPCs and the meshes on the opposite boundaries must be the same in order to create periodic boundaries; see Figure 1. Both the continuum 3D solid element with full integration (labeled C3D8 on ABAQUS) and continuum 3D solid element with reduced integration (labeled C3D8R on ABAQUS) for the mesh in Figure 1a are studied in order to optimize the computation efficiency on the RVE. Following the constitutive model introduced in the last subsection, the material parameters of a conventional NiTi alloy are given in Table 2, which are identified with the experimental data of Sittner et al. [43].
To check the effective response of the RVE, a tensile strain up to 10% in the Y direction is applied on the RVE at a constant temperature of 30°C.  To describe the effective behavior of RVE, the averaged values, such as stressσ, strainε and martensitic volume fractionf t , are introduced. The averaged values are computed by volume averaging over the whole 1 mm 3 cube from Element 1 to N, where N denotes the total element number. Note the averaged strain is equal to the prescribed strain due to the periodic boundary conditions, Equation (9). Figure 2a shows the stress-strain relations simulated with different meshes and element types. The hysteresis effect is observed during the loading/unloading cycle. This nonlinear response is caused by the inner forward phase transformation of SMA. The averaged martensitic volume fractionf t increases along with the loads until reaching a maximum value; see Figure 2b. When the unloading begins,f t decreases immediately until recovering to zero. From the comparison of the curves, it is observed that the mesh with 720 elements is fine enough to model this architected SMA RVE. Furthermore, the RVEs simulated by the C3D8R and C3D8 elements have exactly the same responses; see Figure 2.
Generally speaking, the mesh with 720 C3D8R elements has relatively enough accuracy and lower computational cost than other meshes. Thus, meshes like density and C3D8R elements are used for the RVE in the following work.

Cells with Different Geometries
In this subsection, five types of RVEs with different void volume factions ξ are studied, as illustrated in Figure 3. The size of the cube is given by 1 mm × 1 mm × 1 mm, while the radius of the cylindrical hole varying along with ξ is given in Table 3. The material parameters of SMA remain the same as those in the last subsection.  A tensile strain up to 10% and a compressive strain up to −10% in the Y direction are respectively applied on the RVEs to simulate effective responses for different RVEs. Figure 4 gives the curves of averaged stress versus the averaged stain along the loading direction, simulated by the above RVEs respectively. The absolute value of stress in different RVEs at a certain strain level increases along with the decreasing of the void volume fraction, as depicted in Figure 4. In more detail, Figure 5 gives the relation of maximum stress and material volume fraction (1-ξ) when the absolute value of strain is up to 10%. The trend of the curve shows that the higher the material volume fraction is, the faster the stress increases. This trend could be explained by exploring the stress distribution on the RVE at the maximum strain level 10%, as shown in Figures 6 and 7. The high stresses are mainly distributed over the middle of the pillars along the loading direction. Lower material volume fraction results in a stronger stress concentration effect. Thus, the stiffness of the RVE decreases more rapidly than the material volume fraction. This test could be a reference for designers to balance the stiffness and the mass of the architected structure.   Figure 8 gives the martensitic volume fraction averaged over SMA versus the strain averaged over the cube. The martensitic volume fractionf t increases along with the increasing of the strain level, but stops before it reaches one. More specifically, each RVE has a certain maximumf t value in this loading case. The maximum martensitic volume fraction versus the material volume fraction at the maximum absolute strain level of 10% is depicted in Figure 9. It shows that the lower the material volume fraction is, the more rapidly the martensitic volume fraction decreases. This also means the usage rate of SMA is at a very low level when the material volume fraction is low, and it could result in a waste of SMA. For example, the hysteresis effect of the SMA architected structure is not proportional to the SMA mass, but to the product of the SMA mass and the SMA usage rate. Thus, the hysteresis effect could be very weak when the SMA usage rate is very low. Like the stiffness discussed before, this phenomenon could also be explained by the stress concentration where only the high stress area transforms into martensite; see Figures 10 and 11. This also means SMA cannot completely transform from austenite into martensite. Thus, as long as the unloading begins, the reverse transformation starts immediately, resulting in the reduction off t .       To design a light weight architected SMA structure, it is required to minimize the material volume fraction in order to decrease the weight and cost of the structure. However, the necessary stiffness of the structure should be firstly satisfied. It should be noticed that the waste of SMA should also be avoided as much as possible due to SMA's high price. Therefore, the results given in Figures 5 and 9 can be a reference for the designers to balance the stiffness, mass and SMA usage according to their specific requirements.
Following the tests at the RVE level, the structural response of the architected SMAs is studied in the next section.  Figure 11. Distribution of martensitic volume fractionf t (named SDV7 in colorbar) on the RVEs with different void volume fractions ξ in compression loading at strain level −10%.

Structural Response
In this section, the structural responses of the architected SMAs are focused on. The constitutive behavior of the macroscopic structure is represented by the computational homogenized RVEs in the last subsection. To simulate the structural response, a short introduction of the multilevel finite element method (FE 2 ) is presented firstly.

FE 2 Formulation
Generally speaking, an architected SMA with a periodic microscopic structure can be divided into two scales and be simulated by the multiscale homogenization method. As depicted in Figure 12, each point at the macroscopic level is represented by a periodic RVE after homogenization. The constitutive behavior and the stress of each macroscopic point are transferred from the RVE by the computational homogenization technique, while the macroscopic strain is applied on the RVE with periodic boundary conditions. Both levels are simulated by FEM, which can capture their mechanical fields accurately. This approach is implemented on ABAQUS via the user subroutine UMAT. The real-time interaction of the two levels is realized by the iteration of the Newton-Raphson method. This method is able to compute both the macroscopic and microscopic responses of the structures simultaneously. The authors have developed this FE 2 approach on the ABAQUS platform for SMA-based composites; see Xu et al. [39] for more detailed formulations. Related valuable works on multiscale modeling of SMA composite could be also referred to, such as Kohlhaas and Klinkel [37], Chatzigeorgiou et al. [38], Chatzigeorgiou et al. [44] and Fatemi Dehaghani et al. [45].

Beam with Three Kinds of Cells Subjected to Three-Point Bending
The multiscale finite element method mentioned above has been validated by the reference and shows good ability in modeling the superelasticity and the shape memory effect of the SMA composites [39]. Therefore, this multiscale model is used to simulate the multiscale response of the architected SMA structure herein. A 3D beam subjected to three-point bending load is shown in Figure 13.
This beam is composed of architected SMA. The width, length and height of the beam are 5 mm, 20 mm and 5 mm, respectively. Edge Y = 5 mm, Z = 10 mm is fixed in the Y and Z directions. A displacement load up to 0.5 mm in the Y direction is applied on the edges Y = Z = 0 mm and Z = 20 mm, Y = 0 mm. Node X = 0 mm, Y = 5 mm, Z = 10 mm is fixed in the X direction in order to eliminate the rigid body displacement in the X direction. Considering the symmetry of the structure and the boundaries, only the left half of the structure is simulated in order to reduce the computation cost. To do this, an additional displacement constraint in the Z direction is given on face Z = 10 mm. Since the deformation in the X direction is not obvious in the three-point bending test, one element is used in this direction. Each edge in the Y direction is meshed by two elements and in the Z direction by four elements; see Figure 14. The continuum 3D solid element with incompatible modes (labeled C3D8I in ABAQUS) is adopted for the modeling of the macroscopic beam since it is enhanced by incompatible modes to improve its bending behavior. The RVEs studied in the last section are used herein with void volume fractions ξ of 40.7%, 72.5% and 82.2%, respectively.
As the macroscopic constitutive behavior on each integration point is not clear, it has to be represented by the effective behavior of the associated RVE at each macroscopic increment. A brief flow diagram, showing how this multiscale problem is solved, is illustrated in Figure 15. Specifically, the effective behavior of the RVE is computed by seeking the relations between the averaged stresses and averaged strains over the RVE via a series of loading tests. Once the effective behavior is obtained, the macroscopic problem is to be solved. Considering the nonlinear response of the RVE during loading, the macroscopic convergence has to be checked in each macroscopic iteration of the Newton-Raphson framework. In an iteration, the strain states of RVEs are updated with the macroscopic strains, and in return, the macroscopic stresses are renewed by updating the averaged stress of the RVE at the new strain states.   Next increment Figure 15. The nonlinear interaction between two scales in the Newton-Raphson framework.

Stress Distributions at the Macroscopic and Microscopic Levels
Let us consider the case with volume fraction ξ = 40.7%. The stress distribution of the bending beam with the prescribed displacement in the Y direction reaching 0.5 mm is illustrated in Figure 16. The deformations of both the macroscopic and microscopic structures are magnified by three times. The high compressive stresses in the Z direction are observed in the elements above the middle plane Y = 2.5 mm, while in contrast, the compressive stresses are observed in the elements below the middle plane. In the meantime, for different macroscopic points, the RVEs have different stress states in correspondence with their associated macroscopic strain states. To specify the microscopic structure more clearly, we introduce RVE kl to denote an RVE corresponding to the integration point l of the macroscopic element k. The RVEs above the middle plane, such as RVE A4 and RVE D4 , are mainly subjected to compression, while the RVEs below the middle plane, such as RVE E3 and RVE H3 , are mainly subjected to tension. The stresses in the RVEs far from face Z = 10 mm, such as RVE A4 and RVE E3 , are at a relatively low level compared to those RVEs near face Z = 10 mm, such as RVE D4 and RVE H3 .  Figure 17 shows the nonlinear response of the macroscopic structure with three different void fractions. The linear response is observed at the very beginning, and the microscopic structures deform without any phase transformation. As the loading increases, the forward phase transformation begins over the high stress area of the beam. For example, Figure 18 gives the stress strain relations of integration points in element D and element H with ξ = 40.7%, which also represent the effective behavior of the associated RVEs. Note that Integration Points 5 to 8 in each element are not illustrated considering the symmetry of the structure in the middle plane X = 2.5 mm. It is observed that the hysteresis effects in RVE D3 and RVE H3 are relatively much more obvious than those in RVE D4 and RVE H4 ; because RVE D3 and RVE H3 are located in the high stress area. The martensite transformation states of the RVEs with the prescribed displacement in the Y direction reaching 0.5 mm, shown in Figure 19, also give a reasonable confirmation from the RVE level. This figure depicts the RVEs corresponding to these integration points closest to the face Z = 10 mm. Only the RVEs far from the middle plane Y = 2.5 mm have an obvious phase transformation. Simultaneously, the phase transformation at the RVE level is accompanied by the softening of macroscopic stiffness. This softening goes gradually, because the thermomechanical phase transformations at the microscopic level are not synchronous, considering that the stress states of the macroscopic structure are not uniform over the whole beam and that the microscopic structure is architected. Once the unloading begins, the behaviors of the structure and RVEs immediately turn into the linear case. Later, the reverse transformation follows, and finally, all the phases transform back into austenite.   Figure 19. Distribution of martensitic volume fractionf t (named SDV7 in colorbar) on the RVEs with the prescribed displacement in the Y direction reaching 0.5 mm and ξ = 40.7%.

Structural Response with Different Microscopic Structures
The differences between structural responses with three void fractions are observed in Figure 17. With a higher void fraction, the macroscopic structure shows lower stiffness and a weaker hysteresis effect, because the macroscopic behavior is represented by the mean behavior of the microscopic structures. The stress-strain relations of Integration Points H3 and D4, which are also the averaged stress-strain relations of RVE H3 and RVE D4 , are depicted in Figure 20 to investigate the microscopic response of the structure. Figure 21 gives the stress distributions in the Z direction for associated RVEs at maximum loads corresponding to the displacement of the beam ends U 2 = 0.5 mm. It is observed that the averaged stiffness and the hysteresis effect of the RVE with a higher void fraction are lower than those of the RVE with a lower void fraction, which is consistent with the response at the macroscopic level. As we have discussed in Section 2, the high void fraction results in less SMA, which can provide low stiffness for the structure. For this architected structure, the stresses are concentrated over the pillars along the main loading direction when the void fraction is high. This also results in the martensite phase being mainly concentrated over the pillars.

Comments on the Computational Efficiency
The proposed numerical tool contributes to the design of innovative SMA applications, which can optimize the quantity of material with the required functionalities (actuation, recovery force, damping, energy absorption, etc.), lighten the structure and consider alternative manufacturing processes such as additive manufacturing. The latter can consider more complex geometries of the unit cell, which has a direct influence on the overall behavior of the structure. Note that this numerical multiscale procedure is very time consuming. To reduce the CPU time, the authors propose to associate this method with model reduction techniques. Among these techniques, we can introduce the proper orthogonal decomposition (POD) to build up a small number of basis functions where the solution can be computed (Yvonnet et al. [46]). When this technique is used in an optimization procedure requiring several repetitive calculations, it can significantly reduce the computation time. We can also consider the proper generalized decomposition (PGD) technique (Ammar et al. [47]; Kpogan et al. [48]), which can reduce the dimensionality of the problem to be solved and which is well adapted for unsteady problems, as well as parametric studies. Another technique that has proven its efficiency in solving non-linear problems and that was developed by our group is the asymptotic numerical method (ANM). This technique is based on the development of variables in the form of Taylor series truncated at large order, which allows one to reduce the computation time significantly by reducing the number of computation steps. This technique can be considered as a high order predictor algorithm without the need for any iteration phase. It is particularly suited for strong nonlinearities and for thin-structure problems involving instability phenomena (Nezamabadi et al. [49]; Assidi et al. [50]; Aggoune et al. [51]). In addition, bridging techniques (Hu et al. [52]; Yu et al. [53]) can also be considered to combine the advantages of different reduced finite element models. These tools will be associated with the present algorithm to reduce the computation time and to deal with complex structural geometry and complex response curves.

Conclusions
In this paper, the superelasticity behavior of the architected SMA structure is studied with a 3D multiscale finite element model. Firstly, RVEs with architected SMA are built to simulate the cellular responses of the structure. The behavior of SMA is described by the constitutive model proposed by Chemisky et al. [40]. Both tension and compression loading cycles are applied on the RVEs with different void fractions. The superelasticity responses and the hysteresis effects are observed in the RVEs. The effect of changing the void fraction on the stiffness, the maximum martensitic volume fraction and hysteresis effect are discussed in detail. Moreover, a multiscale approach is carried out to model the structural response, as well as the cellular response. The relations between the structural and cellular responses are studied in a three-point bending test. It is observed that the macroscopic response is related to the phase transformation in the RVE, which changes the effective constitutive behavior of the structure. The stress-strain state of the RVE directly depends on the stress strain state of the associated macroscopic point. Thus, the multiscale model is necessary and successful for simulating the nonlinear behavior of the architected SMA structure. Furthermore, structural responses with different void fractions are studied, which gives a good reference of the void fractions' influences on structural stiffness and hysteresis.

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

Abbreviations
The following abbreviations are used in this manuscript:

SMA
Shape memory alloy RVE Representative volume element FE 2 Multilevel finite element method FEM Finite element method UMAT User-defined materials MPCs Multi-point constraints RP Reference point C3D8 Continuum 3D solid element with full integration C3D8R Continuum 3D solid element with reduced integration C3D8I Continuum 3D solid element with incompatible modes