Optimal Design of Shape Memory Alloy Composite under Deflection Constraint

Shape-adaptive or morphing capability in both aerospace structures and wind turbine blade design is regarded as significant to increase aerodynamic performance and simplify mechanisms by reducing the number of moving parts. The underlying bistable behavior of asymmetric cross-ply composites makes them a suitable candidate for morphing applications. To date, various theoretical and experiential studies have been carried out to understand and predict the bistable behavior of asymmetric laminates and especially the curvature obtained in their stable configurations. However, when the bi-stable composite plate is integrated with shape memory alloy wires to control the curvature and to snap from a stable configuration to the other (shape memory alloy composite, SMAC), the identification of the design parameters, namely laminate edge length, ply thickness and ply orientation, is not straightforward. The aim of this article is to present the formulation of an optimization problem for the parameters of an asymmetric composite laminate integrated with pre-stressed shape memory alloys (SMA) wires under bi-stability and a minimum deflection requirement. Wires are modeled as an additional ply placed at the mid-plane of the composite host plate. The optimization problem is solved numerically in MATLAB and optimal design variables are then used to model the SMAC in ABAQUS™. Finite element results are compared against numerical results for validation.


Introduction
Morphing structures applications are of interest in both aerospace structures [1] and wind turbines [2], which take benefit of significant shape changes to increase performance and efficiency. Asymmetric composite laminates (ACL) may be formulated such that they possess two stable shapes, each of which is a natural equilibrium position. Therefore, ACL can settle at either of them without the need of an external aid. The bi-stable behavior of ACL is related to the residual stress field that develops after curing due to a mismatch in coefficient of thermal expansion in an unsymmetric stacking sequence. The application of bending and twisting moments with respect to mid-plane, can result in a new internal stress equilibrium at a second stable configuration. The state change can be obtained by inducing in-plane strains, typically by smart actuators such as piezoelectric and shape memory alloys (SMA). The intrinsic anisotropy of composite structures was first exploited by Hyer [3,4] using unsymmetric layups to produce bi-stable morphing structure at room temperature. Jun and Hong [5] modified Hyer's theory by taking into account in-plane shear strain. Schlecht and Schulte [6] firstly provided the bi-stable behavior of asymmetric laminates using MARC Finite Element Analysis (FEA) software. Tawfik et al. [7] presented a finite element approach using ABAQUS™ to predict the unsymmetric laminate shapes under thermal curing stresses.
Shape memory alloys have functional properties associated with the pseudo-elastic behavior and shape memory effect (SME). The pseudoelasticity effect allows the alloy to deform elastically up to 7-8% strain, that means about forty times the elastic regime of most steels. Shape-memory alloys can also remember their original shape (SME), after being pseudoplastically deformed below a given temperature (transformation temperature), just deformed by increasing their temperature. The SME is due to the thermoelastic austenite-martensite transformation, which can be described as a diffusion-less process. One of the major advantages of SMAs over other actuator materials is their large recovery force per unit area generated by the phase transformation following the initial deformation. In addition, due to its high tensile strength, SMAs are appropriate for use in small-sized, high-output actuators. Several SMA constitutive models have been proposed and developed over the past few decades to describe the SMA phase transformation phenomenon [8][9][10][11][12][13][14]. However, SMA based actuators [15] have a slow actuation speed which restricts their use in most of the nanotechnology applications such as RF filters, fluid sensors, or devices for the quantum state measurement. Stachiv et al. [16] design tunable high frequency SMA microcantilevers resonators and demonstrated that the phase changeable SMA film can be used by enabling tunability of several consecutive resonant frequencies. Here, in contrast to the usual thermally actuated SMA, the high frequency actuation is generated by the elastic substrate material. Furthermore, Stachiv et al. [17] also proposed nanocantilevers with a tunable spectrum of the resonant frequencies and changeable static deflection utilizing the phase transformation of NiTi film sputtered on the elastic substrate material.
An SMA-embedded composite could be used as an integrated smart actuator-structure without any additional or separated actuation device. As a result, it offers a significant contribution to the design of lightweight structures. Ryu et al. [18] verified actuation of asymmetric laminate using SMA spring actuator through the comparison between experiment and numerical simulation. Dano and Hyer [19] used a mechanism wherein, after SMA wires were stretched between a system of support placed above the laminate and upon electrical heating of the SMA wires, it was concluded that is possible to use SMA wires to change the shape of unsymmetric laminates and to predict reasonably well the overall response of the laminate as a function of the SMA wire temperature. Hassanli and Samali [20] investigated the buckling of curved laminated composite panels reinforced with SMA fibers. A different finite element method was used by Tawfik et al. [21] to examine the stability behavior of SMA composite panels. In addition, Von Karman nonlinear strains were considered in the formulations. Turner [22] examined the thermoelastic response of SMA hybrid structures by finite element analysis. Niknami et al. [23] investigated the effect of induced heat generations on impact responses and phase transformations of hyrbid SMA composite plate through proposing a refined Helmholtz free energy expression and refined constitutive and contact laws, in addition to employing a return-map Newton-Raphson method for enhancement of the numerical solution algorithm. Birman et al. [24] illustrate that SMA fibres embedded within the layers of a composite plate can significantly enhance its global resistance to low-velocity impact and the effectivness of SMA fibre can be further improved by optimizing their distribution throughout the plate. Gandhi et al. [25] presented the possibility to trigger the snap-through from one stable configuration to another by means of SMA wires embedded into the laminate.
Moving on from the work done in [25], the aim of the present study is to predict the optimal design variables for an asymmetric, square laminate with embedded SMA wires such that the laminate can attain a desired stable configuration after the cure cycle and a second stable configuration is admissible, that can be attained by snap from the first one actuated with the help of SMA wires. However, the need of a detailed investigation to identify the optimal design variables necessary to obtain the actuation, requires the definition of an optimization technique in order to avoid a time-consuming trial-and-error procedure, that is based on intuition with no prior knowledge about the obtained configuration and out-of-plane displacement. While a square, asymmetric laminate was selected based on the fact it allows a bistable behavior and a theoretical framework for optimization without SMA is already present in the literature [26], it is acknowledged that this will not cover all laminate configurations which may be considered for practical applications as, for example, rectangular laminates. The optimization study considers the design of bi-stable laminates through variation in ply orientations, ply thickness (same for all plies), and laminate edge length. An objective function has been defined to optimize the ratio of bending stiffness in two chosen direction of the laminate in order to control the shape of the laminate after the curing process. The search space comprises of multiple local and global optima that vary with the bounds posed on the design variables and the deflection constraint. MATLAB's sequential quadratic programming method [27,28] is used to solve the optimization problem. Given the complex and multimodal nature of the laminate design problem, MultiStart Algorithm is used to generate random points within the bounds to capture all local optima and search thoroughly for the global minimum. The optimization problem is solved, as an example, for [45 • , 0] and [75 • , 0 • ] low and high stiffness bending direction to find optimal geometric design variable when the laminate is subjected to deflections and bi-stability constraints. The accuracy of the optimal design is verified against SMAC model implemented in ABAQUS to predict the shape after curing and the response of the model.

Constitutive Model of SMA Wires
The constitutive model for SMA response to mechanical and thermal load is available as a user-defined material model (UMAT) for the ABAQUS™ 6.13 software [29], that is a FORTRAN (Formula Translation) code numerical implementation compiled against the Intel Fortran compiler v13.1 and linked into an ABAQUS executable via Microsoft Visual Studio 2012 Professional (Microsoft, Redmond, WA, USA). This particular constitutive model is rigorously developed to ensure that the key constitutive relations and evolution equations are thermodynamically sound. A detailed discussion of the model derivation process can be found in the works of Lagoudas et al. [13], and so will not be covered in this work. However, a brief summary of the resulting constitutive relations and transformation conditions is outlined here.
The total strain tensor ε is additively decomposed into elastic strain ε el , thermoelastic strain ε th , and transformation strain t : For simplicity, this work assumes that the only inelastic strain component is the transformation strain t and does not take in account for the evolution of transformation-induced plastic strain in SMA or consider the reorientation of variants of martensite.
The constitutive equation is stated as: where α is a second-order thermal expansion coefficient tensor and S is a fourth-order compliance tensor. The compliance tensor is a function of the evolving martensite volume fraction ξ and the compliance tensors S A and S M for the austenite and martensite phases, respectively, and is formulated as a rule of mixtures given by: The evolution equation, which is responsible for relating the time rate of change of transformation strain with the time rate of change of the internal state variable ξ, assumes the form of a flow rule: . Here, the transformation tensor Λ indicates the direction of transformation with the branching function: where H max is a material parameter associated with the maximum transformation strain, and t−r and t−r are the transformation strain and the effective transformation strain at the point of transformation reversal (i.e., the point at which the material stops forward transforming and begins reverse transforming). For reference, the deviatoric stress σ is given by: while the associated effective von Mises scalar measure of the stress tensor σ is defined by: The effective transformation strain at the reversal of the phase transformation is given by: The phase transformation of the SMA is described by a transformation function, Φ = Φ(σ, T, ξ), such that: where Y is the critical thermodynamic driving force necessary to initiate transformation, the current value (π) of which is given by: where ∆(·) refers to difference in the material properties between the martensite and austenite phases, i.e.; ∆(·) = (·) M − (·) A , ρ is the material density, c is the specific heat capacity, s 0 and u 0 are the respective reference entropy and internal energies, and f (ξ) is a chosen transformation hardening function. As seen in Equation (10), transformation occurs when π value equals Y. Thus, during transformation Φ is necessarily zero. From the above relations, the evolution of the martensitic volume fraction is said to be governed by the set of constraints called the Kuhn-Tucker conditions. These conditions are concisely stated as: With the equations described here, the constitutive model is implemented in the finite element framework.

Mathematical Modelling for Composite Laminate
The model to predict the cured shape is based on a nonlinear extension of the classical laminated plate theory with approximated midplane strain functions and nonzero in-plane shear strain. The coordinate system used is such that origin is placed at the geometric center of the laminate and plies are defined in order starting from the bottom (negative thickness coordinate) surface. The midplane strains including geometrical nonlinearity according to the von Karman hypothesis are defined as: where u 0 , v 0 and w 0 are the in-plane displacements in the x-, y-and z-directions respectively. Therefore, the total strains distribution is given by the Kirchhoff hypothesis as: where ε 0 xx , ε 0 yy , γ 0 xy and κ 0 x , κ 0 y , κ 0 xy are the total midplane strains and curvature respectively. Since the composite laminate is subjected to the thermal load, the constitutive law of off-axis stress-strain relation in the macro-mechanical behavior of a lamina can be presented as stress-strain relations of the kth layer in a multilayered laminate.
where α's are the thermal expansion coefficients in the x-y plane and are defined as: cos 2 θ sin 2 θ −2 cos θ sin θ sin 2 θ cos 2 θ 2 cos θ sin θ cos θ sin θ − cos θ sin θ cos 2 θ − sin 2 θ where α 1 and α 2 are the material thermal expansion coefficients in the fiber and transverse directions respectively and ∆T denotes the change in temperature from cure to operating temperature. The transformation stiffness Q ij and the transformation matrix presented in the following matrix form: [Q] k is the reduced stiffness matrix of the laminate defined as follows: where E 11 is the longitudinal Young's modulus, E 22 is the transverse Young's modulus, G 12 is the shear modulus, υ 12 is the major Poisson's ratio, and the minor Poisson's ratio υ 21 . The respective values used in this work are shown in Table 1. Table 1. Material properties of M21/T800 prepreg sheet data from [30].

Material Properties Value
The resultant forces and moments acting on a laminate are obtained by integration of the stresses in each lamina through the laminate thickness. The entire collection of force and moment resultants for an N layered laminate for mechanical and thermal strain are defined as, Since the middle surface strains and curvatures are not a function of z (because these values are always at the middle surface z = 0), they need not be included in the integration. Also, the laminate stiffness matrix is constant for a given ply so it will be a constant over the integration of a lamina thickness, too. Substituting the stress-strain relation of Equation (14) and pulling these constants to the front of the integral, we obtain the following constitutive relation for laminate in matrix notation: where, The A ij are extensional stiffnesses, the B ij are bending-extension coupling stiffnesses, and the D ij are bending stiffnesses. The presence of the B ij implies coupling between bending and extension of a laminate, because both forces and curvatures as well as moments and strains simultaneously exist.

Integration of SMA Wire in Composite Laminate
The additional SMA layer is accommodated at the mid-plane of the composite host plate, assuming that epoxy resin fills in the space within wires such that even and an equal number of composite plies are placed on the top and bottom surface of epoxy/SMA layer as shown in Figure 1a. This approach requires a volume fraction method to combine the properties of the SMA and epoxy taken from [24]. Accordingly, the thermo-elastic properties of an epoxy/SMA layer are written as where the "m" and "s" subscripts stand for the composite matrix and SMA fibers, respectively, and the value of V s is where n, d s and t s stands for number of SMA wires, the diameter of SMA wire and epoxy/SMA ply thickness respectively. Material parameters E s 11 , E s 22 , G s 12 , υ s 12 , α s 11 , α s 22 , Q s ij , denote therefore the values of the corresponding SMA/epoxy ply. The homogenization of SMA and epoxy properties across a ply according to Equations (23) and (24) has been chosen as a simple way to introduce SMA wires contribution in the plate. Of course, it is recognized that, when a small number of wires is considered, this is a simplifying assumption since the layer containing SMA wires is strongly heterogenous. The material parameters of the epoxy/SMA ply are evaluated using the values reported in Tables 2  and 3. The stress influence coefficients values as shown in Table 2 has been modified to keep the consistency of sign in the article. However, for FEA simulation these coefficients must be modified back to its original form of ρ∆s A and ρ∆s M which can be achieved as follows Since there is a large change in Young's modulus of SMA with phase transformation, which is reflected in Equation (23), leads to established a relation between the SMAC stiffness matrix and ξ. However, for the number of wires that can be embedded in practice, as considered in this paper, the value of V s makes the contribution of SMA wires to the stiffness of the laminate negligible. Based on this, it is possible to assume that the change in the E s due to phase transformation can be neglected and E s can be set to a constant value equal to that at the end of the transformation, simplifying the optimization problem. The appearance of the Martensite Volume Fraction (MVF), , and the Transformation Strain (TRNS), , in the stress-strain relation for the epoxy/SMA ply in Equation (26) manifest the presence of the SMA wires on the resultant force developed inside the epoxy/SMA ply. Hence, the total forces per unit length of the SMAC can be expressed as: where, stands for resultant force in the epoxy/SMA ply which can be calculated by substituting Equation (26) into Equation (19). The linear superposition principle can be used to calculate the total internal force based on the taken assumption that stiffness computed for epoxy/SMAC is constant. Hence, the net force caused by laminate and epoxy/SMA ply is the sum of the responses that would have been caused by each stimulus individually.

Optimization of SMAC
Betts [26] presented an optimization technique for the design of bistable laminates based on an analytical solution for an unsymmetric laminate design. A laminate of arbitrary orthogonal layup with two stable shapes with equal and opposite curvatures was considered by applying the following design rules: 1. Even number of groups of plies, where those below the laminate midplane are rotated 90° with respect to the corresponding ones above the midplane (see Figure 1a); in this work, the number of groups of plies above the midplane (n in Figure 1a) is set to 2; 2. Square edge length ; 3. Each ply thickness about the laminate midplane, , , …, ; 4. Each ply is made of the same material.    [29].

Material Parameter Value
Elastic stiffness of the austenite E A 70 GPa Elastic stiffness of the martensite E M 30 GPa Poisson's ratio (equal for both phases) υ 0.33 Coefficient of thermal expansion for the austenite Table 3. Material properties of Epoxy Resin data from [31].

Mechanical Properties Value
Tensile Modulus The constitutive equations of off-axis stress-strain relation for epoxy/SMA ply with arbitrary wire orientation can be presented as follows [23]: The appearance of the Martensite Volume Fraction (MVF), ξ, and the Transformation Strain (TRNS), t , in the stress-strain relation for the epoxy/SMA ply in Equation (26) manifest the presence of the SMA wires on the resultant force developed inside the epoxy/SMA ply. Hence, the total forces per unit length of the SMAC can be expressed as: where, N s x stands for resultant force in the epoxy/SMA ply which can be calculated by substituting Equation (26) into Equation (19). The linear superposition principle can be used to calculate the total internal force based on the taken assumption that stiffness computed for epoxy/SMAC is constant. Hence, the net force caused by laminate and epoxy/SMA ply is the sum of the responses that would have been caused by each stimulus individually.

Optimization of SMAC
Betts [26] presented an optimization technique for the design of bistable laminates based on an analytical solution for an unsymmetric laminate design. A laminate of arbitrary orthogonal layup with two stable shapes with equal and opposite curvatures was considered by applying the following design rules: 1.
Even number of groups of plies, where those below the laminate midplane are rotated 90 • with respect to the corresponding ones above the midplane (see Figure 1a); in this work, the number of groups of plies above the midplane (n in Figure 1a) is set to 2; 2.
Square edge length L;
Each ply is made of the same material.
The stacking sequence, illustrated in Figure 1a, was selected as it provided scope for tailoring the directional stiffness properties, while the orthogonal alignment of plies and the square laminate profile enabled maximum useful deflection between states.

Constrained Optimization Problem
The constrained optimization formulation is introduced as follows: The optimization problem is solved using MATLAB's sequential quadratic programming (SQP) method, fmincon. The essential idea of SQP is to model (27) at the current iterate by a quadratic programming subproblem and to use the minimizer of this subproblem to define a new iterate. The keynotes about SQP algorithm is discussed in Appendix B. Furthermore, due to the presence of multiple global and local minima in the optimization domain, MATLAB's MultiStart algorithm is used to generate potential initial points that are uniformly distributed within the domain in order to identify all minima.
The SMAC structure considered to be fixed at its center when the temperature change from cured ∆T = −160 K (cure 453 K, ambient 293 K) is applied to the structure. In the case of FE simulation, also the out-of-plane displacement was constrained in order to simulate the presence of the mold.

Objective Function
The formulation of the optimization problem maximizes the bending stiffness in a given direction of loading while at the same time the bending stiffness in the direction of snap to the second stable configuration is minimized (see Figure 1c). The objective function that optimizes the ratio of bending stiffness in two chosen direction of the laminate is represented by Equation (28), where ϕ 1 is the direction of low bending stiffness, ϕ 2 is the direction of high bending stiffness and δ is the differential operator, The bending stiffness is characterized as a function of ϕ by the change of a shape coefficient a ϕ with respect to a moment applied in that direction M xϕ . The definition of the shape coefficients is given in Appendix C. The bending stiffness is evaluated by considering the plate constitutive equations which relate applied forces and moments, N's and M's, to the plate curvatures (or shape coefficient), a ϕ , b ϕ and c ϕ , and midplane strains, ε 0 's, using the A, B and D matrices of SMAC and substituting into the plate equation lead to Equation (29). The presence of epoxy/SMA ply is taken in consideration by evaluating the overall ABD matrices of SMAC, i.e., by replacing reduced stiffness terms of composite (Q ij ) in Equation (22) with Q s ij terms of reduced stiffness of epoxy/SMA ply at the mid-plane of the composite host laminate. With that, it is acknowledged that the incorporation of SMA wires as an additional layer epoxy/SMA layer at the mid-plane of the composite host laminate violates the above design rules 1 and 4 but, since the contribution of this layer to the overall bending stiffness of the SMAC laminate is negligible (see Section 2.2.2), the approximation related to the violation of the design rules is negligible, too.
Equation (29) is then transformed to align with the direction either ϕ 1 or ϕ 2 . Note that the subscript ϕ has been removed from plate curvatures a, b and c in order to avoid confusion in Equation (31). The forces and moments in the transformed direction except M x are set to zero and the remaining system is solved to give an expression in the terms of A, B, D and M x . This is used to determine the objective function of Equation (28) which is to be minimized.

Deflection and Bistability Constraint
The out-of-plane deflection, expressed by Equation (30), is applied as a constraint for the optimization, For a deflection constraint, the out-of-plane displacement at a corner between stable configurations I and II was used (see Figure 1b). The corner deflection between states w de f at x = y = L/2 can be expressed by: The curvatures a ϕ , b ϕ and c ϕ can be expressed in terms of the first state (I) alone, as the two states have equal and opposite curvatures. Therefore, the equation had been rewritten in terms of a ϕI and b ϕI only, A minimum deflection requirement w de f > w de f is used to set a restriction on the allowable design space. The inequality was solved as a part of the numerical solution of the optimization problem. It is noted that, as a matter of fact, Equation (32) incorporates the bistability constraint, a + b > 0.

Simulation of Uniaxial Behavior of an SMA Wire
The inelastic strain in the SMA wire is the result of the stress-induced martensite phase transformation while it is loaded and unloaded at a temperature below A 0s . By subsequently heating it above A 0f , it is possible to recover inelastic strain and return to its parent phase (austenite) by virtue of the reverse transformation. The above phenomenon described as one-way shape memory effect or, simply SME. In this section, SME is exploited to control the shape of a composite plate and to snap its curvature between two stable configurations, therefore a wire of length 300 mm and diameter of 0.1 mm, with material properties listed in Table 2, is first modeled using truss element that can carry only axial loads alike thin, flexible NiTi wires. This in the attempt to emulate the cure cycle of the SMAC and to gain insight about stress distribution, the evolution of martensite volume fraction and transformation strain. The analysis consists of three steps: (i) Pre-strain of the wire at environmental temperature up to 6% elongation; (ii) heating of the wire up to the cure temperature of the composite plate T cure = 453 K, at constant strain; (iii) cooling down to an environmental temperature at constant strain. The results obtained at the end of the first step are also used as an initial condition for SMA wires in the SMAC laminate model as described in Section 3.3. The finite element mesh is composed of 30 linear truss elements of type T3D2 in ABAQUS. The ambient temperature is set to 293 K, which is in between martensite start M 0s , and austenite starts A 0s temperatures (see Table 2).
Notice that after the cool-down step the strain is kept as if strain recovery was hampered by the presence of the host composite plate.
As shown in Figure 2a, upon heating thermal expansion yields first a drop in axial stress from 241 MPa to 228 MPa, then the axial stress monotonically increases to a value of 884 MPa as the temperature brought up to T cure = 453 K due to the transformation of de-twinned martensite into austenite. This can be observed in Figure 2b,c, as Martensite Volume Fraction (MVF) decreases from ξ = 1 to ξ = 0.631 and transformation strain (TRNS) from t = 0.05 to t = 0.0314. As shown in Figure 2a, upon heating thermal expansion yields first a drop in axial stress from 241 MPa to 228 MPa, then the axial stress monotonically increases to a value of 884 MPa as the temperature brought up to Tcure = 453 K due to the transformation of de-twinned martensite into austenite. This can be observed in Figure 2b, c, as Martensite Volume Fraction (MVF) decreases from ξ = 1 to ξ = 0.631 and transformation strain (TRNS) from = 0.05 to = 0.0314.
In the cool-down step, the stress initially increases from 884 MPa to 899 MPa due to the thermal contraction of SMA ( Figure 2a). As the temperature gradually decreases, high recovery stress which is generated in SMA drops down the stress level which to about 241 MPa at the end of the cooling process. At the same time, MVF gradually increases to unity by accomplishing full phase transformation to de-twinned martensite. It has to be said that, even if the framework chosen to model the constitutive behavior of SMA does not address the transformation of de-twinned into twinned martensite upon cooling, in this case, the martensite remains de-twinned since the final stress level lies above the stress transformation range as shown in Figure 2a, i.e., the martensite is in its fully de-twinned state.

SMAC Numerical Optimization
The material properties used for SMA and M21/T800 are those shown in Tables 2 and 3, respectively. The in-plane dimension of the additional SMA ply is equal to the square composite laminate of edge length L, which is a geometric design variable among with the ply orientation Ө and Ө defining the stacking sequence of a four-ply laminate of constant ply thickness . The SMA ply contains virtually four SMA wires of 0.1mm diameter and length equal to L. The thickness of the SMA ply is set to = 0.2 mm.  In the cool-down step, the stress initially increases from 884 MPa to 899 MPa due to the thermal contraction of SMA ( Figure 2a). As the temperature gradually decreases, high recovery stress which is generated in SMA drops down the stress level which to about 241 MPa at the end of the cooling process. At the same time, MVF gradually increases to unity by accomplishing full phase transformation to de-twinned martensite. It has to be said that, even if the framework chosen to model the constitutive behavior of SMA does not address the transformation of de-twinned into twinned martensite upon cooling, in this case, the martensite remains de-twinned since the final stress level lies above the stress transformation range as shown in Figure 2a, i.e., the martensite is in its fully de-twinned state.

SMAC Numerical Optimization
The material properties used for SMA and M21/T800 are those shown in Tables 2 and 3, respectively. The in-plane dimension of the additional SMA ply is equal to the square composite laminate of edge length L, which is a geometric design variable among with the ply orientation θ 1 and θ 2 defining the stacking sequence of a four-ply laminate of constant ply thickness t. The SMA ply contains virtually four SMA wires of 0.1 mm diameter and length equal to L. The thickness of the SMA ply is set to t s = 0.2 mm.
The direction of low bending stiffness ϕ 1 is set alternatively to 45 • or 75 • and the direction of high bending stiffness ϕ 2 is set based on the orientation of SMA wires, that is 0 • . The pattern of optimum solutions is detailed in Tables 4 and 5 for, where the range of design variables and minimum deflection requirement is also shown. The motive behind the chosen set for w de f is to draw a clear distinction between the curvature acquired by various SMAC configuration simulated in the Abaqus as shown in Section 3.3. for all optimization sets except B2, the last optimal solution in the list shows a small angle between low bending stiffness, ϕ 1 , and ply fiber directions. This results obviously a value of the objective function that is significantly higher, as expected. Therefore, these solutions can be no longer considered global optima. The minimum deflection requirement and the range of L are gradually increased from a first guess value in order to explore a range of design configurations where the optimal solutions tend to set progressively to the highest value of L and the lowest of t. In this way it is possible to determine a larger set of optimal parameters including local and global optima, that is represented by all the results included in Tables 4 and 5. Figure 3 represents the graphical illustration of the objective function for sets D2 and E2 where the global optimum is located at the point θ 1 = −85 • , θ 2 = 26 • is the correct ply orientation in order to obtain low bending stiffness in ϕ 1 in comparison with the rotated cross-ply solution [−80 On the other hand, the rotated cross-ply solution in set D2 and E2 settled down for higher ply thickness as in comparison with the global optima and at these ply orientations, there are many solutions achieving almost equal objective function values. The deflection between states for these solutions vary from 18 to 93 mm. The laminate edge length, which does not affect the objective function, but only the deflection constraint. Therefore, a rotated cross-ply is chosen over the optimal solution with an allowable loss of 12.65% in laminate bending stiffness. Notice that that solution θ 1 = 5 • , θ 2 = −85 • while achieving an objective function value of 0.5272, is discarded since the bi-stability constraint as not satisfied. The minimum deflection requirement and the range of are gradually increased from a first guess value in order to explore a range of design configurations where the optimal solutions tend to set progressively to the highest value of L and the lowest of t. In this way it is possible to determine a larger set of optimal parameters including local and global optima, that is represented by all the results included in Tables 4 and 5. Figure 3 represents the graphical illustration of the objective function for sets D2 and E2 where the global optimum is located at the point θ1 = −85°, θ2 = 26° is the correct ply orientation in order to obtain low bending stiffness in in comparison with the rotated cross-ply solution [−80° /10° ] . On the other hand, the rotated cross-ply solution in set D2 and E2 settled down for higher ply thickness as in comparison with the global optima and at these ply orientations, there are many solutions achieving almost equal objective function values. The deflection between states for these solutions vary from 18 to 93 mm. The laminate edge length, which does not affect the objective function, but only the deflection constraint. Therefore, a rotated cross-ply is chosen over the optimal solution with an allowable loss of 12.65% in laminate bending stiffness. Notice that that solution θ1 = 5°, θ2 = −85° while achieving an objective function value of 0.5272, is discarded since the bi-stability constraint as not satisfied.

FE Simulation of SMAC
In this section, the commercial finite element code ABAQUS is employed to verify the shape Figure 3. Optimization function as a function of θ 1 , θ 2 for D2 and E2 sets with ply thickness of t 1 = t 2 = 0.4 mm.

FE Simulation of SMAC
In this section, the commercial finite element code ABAQUS is employed to verify the shape after the manufacturing cycle of the SMAC with the optimal design variables determined in the previous section. In particular, the value of w def of the optimization problem is compared with the one coming from FEA in order to assess the accuracy of the numerical solution. An L × L mm 2 laminate is modeled in ABAQUS using 4-node, reduced integration shell elements S4R with 4-plies of uniform ply thickness t and a θ 1 /θ 2 /90 • + θ 1 /90 • + θ 2 stacking sequence. The embedded SMA wires are modeled using T3D2 truss elements such that the nodes of SMA wires coincide with the nodes of the laminate and a tie constraint is established between the node region of SMA wires and laminate. The tie constraint between the nodes of SMA wires and the laminate physically implies that pre-strained wires are consolidated in the composite matrix and based on this assumption, initial conditions for SMA wire at the beginning of cool-down step are set to the values listed down in Table 6. The Finite Element mesh consists of 100 linear quadrilateral elements of type S4R and 40 linear line elements of type T3D2. The out-of-mold shape of the laminate after curing is obtained by cooling down from an initial temperature T cure applied to all the nodes of the model to room temperature. All the nodes of the laminate are constrained in the z-direction while the center node of SMAC is fixed to suppress rigid body motions. At the end of the cooling phase, the constraint in the z-direction is deactivated except at the center node to simulate the opening of the mold. Geometric nonlinearity is accounted for in this step. Figure 4 shows the out-of-plane displacement of the SMAC after the cool-down stage for the optimal sets A1-1, B1-1, D1-1, E1-1 in Table 4, while in Figure 5 it is shown for the sets A2-1, B2-1, D2-2 and E2-2. Even though the latter two are not global optima, they are considered since the values of θ 1 and θ 2 are the same of A2-1 and B2-1 cases, therefore facilitating the comparison. It must be underlined that neither L, nor t, nor L/t shows a unique relationship with the deflection, therefore it would have been hard to find the optimal solution for a given deflection simply by trial-and-error. For the results shown in Figure 5, a similar trend can be observed as variations in ply orientations between the laminates are small. Finally, a significant variation in z-displacement can be registered along the y-direction (I stable configuration) in comparison with the x-direction (II stable configuration). Furthermore, a direct correspondence between the attained stable configuration I by SMAC and the objective function criterion to have higher bending stiffness in ϕ 2 = 0 • direction is hard to validate since stable configurations are sensitive to modelled imperfections and uncertainties such as the effect of an additional resin layer which is not modeled in FEA simulations.
The results of simulations have been particularly useful in gaining detailed insight into the distribution of longitudinal and transversal stresses in the SMAC after curing. The case of set B1-1-is reported in Figure 6. The 0 • (180 • )-oriented bottom layer developed compressive stress of −49.31 MPa in fiber direction due to the curvature acquired by the laminate after cool-down step. Note that the peripheral region of the layer is under tension and the inner region under compression, and this state of stress facilitates the snap of the laminate to the second stable configuration [25]. The stresses developed in 90 • -oriented top layer is subject to stress concentrations near to the longitudinal edges: the peak tensile stress is +146.4 MPa right along the edges and +68.2 MPa in narrow regions that run longitudinally a small distance into the shell as shown in Figure 6b.   Table 4.  Table 4.  Table 5. Figure 4 shows the out-of-plane displacement of the SMAC after the cool-down stage for the optimal sets A1-1, B1-1, D1-1, E1-1 in Table 4, while in Figure 5 it is shown for the sets A2-1, B2-1, D2-2 and E2-2. Even though the latter two are not global optima, they are considered since the values of θ1 and θ2 are the same of A2-1 and B2-1 cases, therefore facilitating the comparison. It must be underlined that neither L, nor t, nor L/t shows a unique relationship with the deflection, therefore it would have been hard to find the optimal solution for a given deflection simply by trial-and-error. For the results shown in Figure 5, a similar trend can be observed as variations in ply orientations between the laminates are small. Finally, a significant variation in z-displacement can be registered along the y-direction (I stable configuration) in comparison with the x-direction (II stable configuration). Furthermore, a direct correspondence between the attained stable configuration I by SMAC and the objective function criterion to have higher bending stiffness in = 0° direction is hard to validate since stable configurations are sensitive to modelled imperfections and uncertainties such as the effect of an additional resin layer which is not modeled in FEA simulations.  Table 5.  Table 5. Figure 4 shows the out-of-plane displacement of the SMAC after the cool-down stage for the optimal sets A1-1, B1-1, D1-1, E1-1 in Table 4, while in Figure 5 it is shown for the sets A2-1, B2-1, D2-2 and E2-2. Even though the latter two are not global optima, they are considered since the values of θ1 and θ2 are the same of A2-1 and B2-1 cases, therefore facilitating the comparison. It must be underlined that neither L, nor t, nor L/t shows a unique relationship with the deflection, therefore it would have been hard to find the optimal solution for a given deflection simply by trial-and-error. For the results shown in Figure 5, a similar trend can be observed as variations in ply orientations between the laminates are small. Finally, a significant variation in z-displacement can be registered along the y-direction (I stable configuration) in comparison with the x-direction (II stable configuration). Furthermore, a direct correspondence between the attained stable configuration I by SMAC and the objective function criterion to have higher bending stiffness in = 0° direction is hard to validate since stable configurations are sensitive to modelled imperfections and uncertainties such as the effect of an additional resin layer which is not modeled in FEA simulations.    Table 3. In Figure 7a, the stress falls back from 884 MPa to 228 MPa, which is above the stress value requires to transform de-twinned to twinned martensite hence, the initial value of MVF can be preserved for actuation applications. It can be also deduced that the composition of austenite and de-twinned martensite at the beginning of cool-down transforms into full detwinned martensite phase as the final stress level lies above the stress transformation range. Due to cooling, thermal contraction of SMA occurs, and the downward peak, in Figure 7a,b, indicates the release of the constraint in out-of-plane direction (demolding). The results of simulations have been particularly useful in gaining detailed insight into the distribution of longitudinal and transversal stresses in the SMAC after curing. The case of set B1-1-is reported in Figure 6. The 0° (180°)-oriented bottom layer developed compressive stress of −49.31 MPa in fiber direction due to the curvature acquired by the laminate after cool-down step. Note that the peripheral region of the layer is under tension and the inner region under compression, and this state of stress facilitates the snap of the laminate to the second stable configuration [25]. The stresses developed in 90°-oriented top layer is subject to stress concentrations near to the longitudinal edges: the peak tensile stress is +146.4 MPa right along the edges and +68.2 MPa in narrow regions that run longitudinally a small distance into the shell as shown in Figure 6b. Figure 7 depicts the behavior of embedded SMA wires when the initial conditions are set to the values shown in Table 3. In Figure 7a, the stress falls back from 884 MPa to 228 MPa, which is above the stress value requires to transform de-twinned to twinned martensite hence, the initial value of MVF can be preserved for actuation applications. It can be also deduced that the composition of austenite and de-twinned martensite at the beginning of cool-down transforms into full detwinned martensite phase as the final stress level lies above the stress transformation range. Due to cooling, thermal contraction of SMA occurs, and the downward peak, in Figure 7a,b, indicates the release of the constraint in out-of-plane direction (demolding).  The comparison between the results obtained in ABAQUS and MATLAB is carried out in Tables 7 and 8 that represents the search space given by ϕ 1 = 45 • , ϕ 2 = 0 • and ϕ 1 = 75 • , ϕ 2 = 0 • , respectively. The resulting difference is contained with 10% and it decreases more sensibly for increasing edge length in comparison to the laminate ply thickness. The incorporation of an SMA ply at the mid-plane in the theoretical model of [26], which is exact in the case of pairs of orthogonal plies, is probably at the basis of the difference between FE and MATLAB results. Since the volume fraction of SMA with respect to whole laminate decreases with the square of edge length and linearly with the total thickness the deviation from the original model due to the presence of SMA ply vanishes for increasing laminate size.

Conclusions
The efficient modelling technique for the asymmetric, bi-stable composite laminates becomes even more significant when the possibility of actuation by SMA wires is introduced. Therefore, in this work, a design optimization technique for bi-stable composite plates embedded with SMA wires has been developed to find the suitable design variables for a given initial out-of-plane deflection, without having prior knowledge about a "good initial guess". The optimization problem formulation is based on the minimization of the bending stiffness ratio in two directions, namely ϕ 1 (low stiffness) and ϕ 2 (high stiffness) in order to have a plate that is stiff under loading (membrane) and/or bending in ϕ 2 direction but easy to snap to the other stable configuration by membrane/bending loading in in ϕ 1 direction. The procedure was implemented in MATLAB the analysis was carried out for ϕ 1 = 45 • , ϕ 2 = 0 • and ϕ 1 = 75 • , ϕ 2 = 0 • low and high stiffness bending direction, respectively, demonstrating the efficiency in finding both global and local optima. Furthermore, the initial deflection w de f estimated with MATLAB was compared with a FE analysis done with the software ABAQUS, that yielded a difference within 10%, decreasing down to 1% for increasing laminate size, that confirms the procedure developed in this work as a fast and reliable tool to estimate the optimal design of SMAC bi-stable laminates subjected to an initial deflection constraint. Future work is dedicated to the include in the objective function the SMA wires actuation force necessary to obtain the snap of the laminate to its second stable shape.

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

Appendix A
The stress and temperature conditions for phase transformation involves the determination of transformation "surfaces," or boundaries of the transformation regions in a stress-temperature space. These indicate where a given transformation begins or ends. These four surfaces, and their equations (given below) point out the following, from first to last: The beginning of forward phase transformation; the end of forward phase transformation; the beginning of reverse phase transformation and end of reverse phase transformation, respectively, at zero stress. For a stress-space (σ 11 , σ 22 ), the three-dimensional plot of the transformation surfaces for forward or reverse phase transformation can be presented in the σ-T space. In Figure A1a, a slice of the transformation surface boundaries of a 3D plot at a constant temperature is shown.
Similarly, considering σ 22 = 0 the boundaries of the transformation surfaces in the uniaxial σ-T space is shown in Figure A1b. The form of these boundaries is assumed to be linear (it might be quadratic which depend on the behavior of material) and irrespective of their form, these boundaries are partially described by their intersections with the zero-stress axis and by the slopes at some stress level (e.g., zero stress). Such slopes are known as "stress influence coefficients" [32].

⎩
, , = −Y Similarly, considering = 0 the boundaries of the transformation surfaces in the uniaxialspace is shown in Figure A1b. The form of these boundaries is assumed to be linear (it might be quadratic which depend on the behavior of material) and irrespective of their form, these boundaries are partially described by their intersections with the zero-stress axis and by the slopes at some stress level (e.g., zero stress). Such slopes are known as "stress influence coefficients" [32]. The value of entropy changes in the forward or reverse transformation (exothermic and endothermic processes) can be calculated through a Clausius-Clapeyron-like relation [33]: where, ∆ is a change in the entropy, is transformation strain, and represent stress influence coefficient for austenite and martensite. The entropy relation holds for forward transformation ( → ) by accounting for a decrease in entropy as the process is exothermic in nature however, it violates Le Chatelier-Braun principle for reverse transformation ( → ). This is because does not provide any information about the reversal of the phase transformation hence, it failed to evaluate the increase in entropy during the endothermic process ( → ). This additional information about the phase reversal can be understood by observing the rate of transformation at each time increment therefore, the above relation is modified as follows: Hence, a Clausius-Clapeyron-like relation justifies Le Chatelier-Braun principle for both forward and reverse transformation.
The value of entropy changes in the forward or reverse transformation (exothermic and endothermic processes) can be calculated through a Clausius-Clapeyron-like relation [33]: where, ∆S is a change in the entropy, t is transformation strain, C A and C M represent stress influence coefficient for austenite and martensite. The entropy relation holds for forward transformation ( A → M ) by accounting for a decrease in entropy as the process is exothermic in nature however, it violates Le Chatelier-Braun principle for reverse transformation ( M → A ). This is because t does not provide any information about the reversal of the phase transformation hence, it failed to evaluate the increase in entropy during the endothermic process ( M → A ). This additional information about the phase reversal can be understood by observing the rate of transformation at each time increment therefore, the above relation is modified as follows: . .
Hence, a Clausius-Clapeyron-like relation justifies Le Chatelier-Braun principle for both forward and reverse transformation.

Appendix B.1. Karush-Kuhn-Tucker (KKT) Conditions and the Lagrangian Function
The general Nonlinear Programming (NLP) problem and g i (x) ≥ 0 i ∈ I (active).
The Lagrangian function formulates the problem into one function using Lagrangian multipliers λ for equality constraints and µ for inequality constraints: L(x, λ, µ) = F(x) + i λ i h i (x) + i µ i g i (x). A single function can be optimized by finding critical points where the gradient is zero. This procedure now includes λ and µ as variables. The system formed from this gradient yields the KKT conditions:

Appendix B.2 The SQP Algorithm
Critical points of the objective function will also be critical points of the Lagrangian function and vice versa because the Lagrangian function is equal to the objective function at a KKT point; all constraints are either equal to zero or inactive. The algorithm is thus simply iterating Newton's method to find critical points of the Lagrangian function. Since the Lagrangian multipliers are additional variables, the iteration forms a system: The improvement direction p for Newton's method iterations is found with a quadratic minimization sub-problem that is solved using quadratic algorithms. min p F k (x) + ∇ f T k p + 1 2 p T ∇ 2 xx L k p s.t ∇h k p + h k = 0 and ∇g k p + g k = 0. This problem is quadratic and thus must be solved with non-linear methods, which once again introduces the need to solve a non-linear problem into the algorithm, but this predictable sub-problem with one variable is much easier to tackle than the parent problem. This subproblem can be solved using any Quadratic Programming algorithm.
The midplane strain in Equation (12) are approximated by third order polynomial. Based on, Dano and Hyer [34] findings the form of the midplane strains can be reduced as follows: ε 0 x = d 1 + d 1 x 2 + d 3 xy + d 4 y 2 ε 0 y = d 5 + d 6 x 2 + d 7 xy + d 8 y 2 ε 0 xy = 2d 9 + (ab − c 2 4 + 2d 4 + 2d 6 )x 2 + ( 1 2 ( ac 2 + d 3 ) + d 10 )xy + ( 1 2 ( bc 2 + d 7 ) + d 11 )y 2 , The spatial integrations results in an expression for the total strain energy of the lamintae as follows: This results in 14 nonlinear equilibrium equations to be solved to find the stable shapes defined by the 14 shape coefficients p i . The solutions of the equilibrium equations give the configurations of the laminate at the room temperature. Based on the design rules (Section 2.3), it possible to reduce the above system to just 3 equations and 3 unknows a, b and c and further details about the analytical solution is presented in [26].