Screening significant properties of macro-fiber composite laminate substrate anisotropy with viscoelastic and thermal considerations at the quasi-static level

A macro-fiber composite (MFC) is a class of smart material actuator that employs piezoceramic fibers to extend or contract under an applied electric field. When embedded on a thin substrate, actuated MFCs induce surface pressure, resulting in out-of-plane motion. Design characteristics of the substrate, such as anisotropy, give rise to unique bend/twist behavior. Previous studies on MFC applications have focused on optimizing patch location in relation to the local design to achieve optimal bend/twist motion. However, the resolution to these approaches is restricted to the design space of the application, presenting an absence in piezoelectric laminate design intuition. This research aims to streamline the initial optimization process by analyzing the significant material properties associated with substrate design. To accomplish this, a viscoelastic piezoelectric plate theory model is developed to accurately capture the displacement behavior of the MFC laminate for any given substrate design. Experimental data is used to validate and refine the model, ensuring it closely represents real-world physics. A 2k fractional factorial design of experiments approach is used to identify significant substrate properties per MFC laminate configuration, assigning each material property a two-level coding system. Three configurations are observed, with each exhibiting distinct significant properties and effects compared to others. Furthermore, the study explores the impact of thermal factors on substrate behavior for these MFC laminates, highlighting how significant properties can be temperature dependent. The contributions to the kinematic response due to substrate property perturbation varies uniquely per material property and MFC laminate configuration. Statistical evidence supports the notion that the hysteretic behavior of a material is unaffected by thermal influences and is solely determined by the elastic constants of the substrate. The maximum remanent hysteretic response is shown to be proportional to the maximum actuation response at near room temperature.


Introduction
Smart materials possess the unique capability to serve as both sensors and induce strain.When integrated into loadbearing structures, these materials offer valuable insights into the stress state of the material or enable precise geometric control, provided they are securely bonded to the surface and do not significantly increase weight and stiffness [1].Over the past few decades, they have found increasing relevance in various modern applications, including active vibration control, noise attenuation, aeroelastic stability, damping, morphing, and stress distribution [2][3][4][5][6][7][8][9][10][11][12][13].
Among these applications, morphing structures using piezoelectric actuators is still a relevant topic in research [1,14].However, traditional piezoelectric materials are inherently brittle, making them susceptible to fractures and presenting challenges in handling [15].Furthermore, due to their bulkiness, monolithic piezoelectric materials are difficult to actuate, restricting significant motion when embedded onto a host surface [16].Piezoceramic fiber lamina composites play a significant role in improving the effectiveness of piezoelectric applications by addressing these weaknesses.An exemplary piezoceramic device is the macro-fiber composite (MFC).Originally developed by NASA Langley and later trademarked by Smart Material Corp., MFCs are designed to enhance the flexibility and robustness of piezoelectric actuators [15,17,18].They consist of modular material components, each improving the flexibility and actuation capability of the device (see figure 1).
Several examples can be drawn from the literature that elucidated the capability of MFCs as morphing actuators [6,7,[19][20][21][22].Prior works have observed the bending and twisting behavior of actuated MFCs, or like actuators, on thin structures, with the extent and direction of deformation influenced by the surface geometry, substrate material, and orientation of the actuator [23][24][25][26][27][28][29].However, these studies have neglected the crucial design aspects of substrates.Specifically, the selection of key substrate material properties for general entry-to-design purposes remains elusive, as existing literature tends to optimize substrates based on specific scenarios that may not readily apply to other situations.With the vast array of potential properties and combinations to evaluate for bend/twist kinematic responses, determining the most significant properties becomes somewhat challenging.Furthermore, the thermal and hygroscopic expansion characteristics of MFC laminates, along with their ferroelectric hysteretic behavior [30][31][32][33], have not been systematically classified in relation to substrate design.
The purpose of this work is to streamline the substrate design process by identifying substrate variables that have a significant impact on the bend-twist coupling effect of piezoelectric laminates, taking into account the thermal environment and positional error caused by hysteretic ferroelectricity; hysteresis is due to domain wall motion and nucleation at the boundaries of domains, resulting in dissipation of energy that cannot be recovered [30][31][32][33].The goal is to provide valuable insights for the design process, allowing readers to minimize blind guessing and optimization when working with this smart material.To do so, a finite element (FE) model of a cantilever MFC laminate plate simulating the electrical equivalent of viscoelasticity will be utilized to assess the impact of substrate properties on the resulting kinematic responses.To focus on the extreme bend and twist responses at the laminate's tip, a cantilever boundary condition is chosen; however, the selected boundary condition is arbitrary and is only meant to reflect the optimum design space for other types of boundary conditions.

Multi-surface switching method
To screen the viscoelastic response of an MFC laminate with respect to significant substrate properties, a proposed nonlinear method is established that models the path-dependent ferroelectric behavior using a surface switching method.The basis of this method is consistent with thermodynamic and conservation laws [34][35][36] and has been adapted to mimic kinematic hardening from plasticity models motivated by domain wall motion and nucleation at the mesoscopic scale [37][38][39][40][41][42].The presented method is chosen because of its numerical efficiency and moderate accuracy to real piezoelectric physics, though the author does not suggest it superior to prior hysteretic type models [43,44].This work is based on recent development of a single surface switching method by Laxman et al [42] but is adapted to better fit the convexity of the hysteretic curvature, motivated by Subramanian et al [41].
The rationale behind the developed algorithm should be motivated by deliberation using figure 2. Here, a biaxial plot is used to illustrate the electric field space with respect to inplane orthogonal axes.In this example, two surfaces are used to represent their own percentage group of domains inside a piezoelectric system.Each surface has its own set of properties.The center of each surface is translated by the intermediate value G e , otherwise known as the electric field equivalent of the irreversible spontaneous polarization, or the intensity of biased polarity, within a ferroelectric system.The radius is identified by the coercive field E c .
Point A is the center point of the envelope, and B is the field state of the material with respect to an applied electric field loading E. The magnitude of irreversible spontaneous polarization is represented by the distance between point A and  The two surfaces in figure 2 are concentric for simplicity of demonstration, i.e. they need not be concentric.When only a single surface bound is violated, the portion of domains associated with this surface begins to switch.Since the field vector Ê lies within the outer switching surface, those associated domains remain stationary, and the outer surface does not translate.With MFCs that are spontaneously polarized in a single direction due to the interdigitated electrodes, such as MFC type M8528-P1, the model in figure 2 simplifies from a 2D space to a 1D space.It is then simple to incorporate the nonlinear hysteresis of MFCs into the multi-surface algorithm.
The first portion of the algorithm identifies the initial intermediate value for each surface that needs to be updated in the next iteration.This is accomplished using the derivative of the dependency between the irreversible spontaneous polarization of the ferroelectric material to the internal electric field energy of the system (for some remarks see e.g.[45,46]).The resultant equation to identify the initial intermediate values is then given by, j G e k = a j arctanh where P o and P s are the initial polarization at time t = 0 and the irreversible saturation polarization of the system, respectively.The hardening parameter a j is an empirical parameter that must be determined.The subscript k presents the load step of the iteration while the denotation j selects the switching surface to update among the total number of surfaces N.
Temperature effects influence the actuation capabilities and polarization switching behavior of piezoelectric materials.Primarily, these effects are observed in the resultant coercive fields of the major loop.At elevated temperatures, while not approaching the Curie temperature, the system absorbs energy, reducing the required electric field energy for domain switching compared to lower temperatures.Consequently, the coercive field decreases as the temperature increases [47,48].The multi-surface algorithm can emulate this behavior by [41], where θ T is the temperature of the material system, T c is the Curie temperature, and Υ is a temperature saturation constant.
The factorization constant C is used such that, The effect of C is straightforward.As the temperature rises, C becomes larger and diminishes the coercive fields in the given set of domains.To best fit the actuation behavior of a piezoelectric system, the value for Υ must be appropriately selected.
From there, a yield criterion is developed to indicate if the bounds of the surfaces have been violated.A trial function is developed for the switching criteria of each surface at the (k + 1)th load step, given by, If j Φ e trial > 0, then the surface switching criteria of a given surface is violated, and the intermediate value must be updated accordingly.
To update the intermediate value j G e k at the (k + 1)th load step, the trial function in (3b) is expanded using a Taylor series expansion as a function of the intermediate value.Then by imposing the switching condition Φ e = 0 and truncating the higher-order terms, the form of the expanded equation becomes, where ⟨•⟩ is the Macauley's bracket defined as By algebraically adapting equation (1) to the (k + 1)th load step, G e k+1 is used to estimate the spontaneous irreversible polarization P i k+1 .The following expression is written as, j P i k+1 = α j P s tanh j G e k+1 a j (5) where α j is the percentage of domains in a system associated with the current surface j.The value of α j is less than unity unless only a single surface is used to represent the system, i.e.
N j α j = 1.The hyperbolic tangent term in (5) never exceeds unity regardless of sign.Therefore, the irreversible spontaneous polarization at the (k + 1)th load step can never exceed the saturation polarization P s .
Next, the irreversible spontaneous polarization per surface is summed to represent the total irreversible spontaneous polarization in an MFC, which is done by, From here, the total irreversible polarization is then applied to the proposed adaptive constitutive actuator equation for a ferroelectric system [37] with superimposed thermal expansion considerations [36], and is given by, where S is the compliance matrix of the ferroelectric material, σ is the applied stress vector, ε is the strain vector, d is the piezoelectric strain matrix, E is the electric field vector, α is the coefficient of thermal expansion (CTE) vector, and ∆T is the scalar change in temperature.The polarization ratio in the second term on the right-hand side was proposed by Kamlah [37] who postulated that the piezoelectric effect is non-trivial if there exists some irreversible polarization in the system.The change in irreversible ferroelectric strain ∆ε i e depicts the hysteretic strain contribution to the system.The ∆ signifies that the term is with respect to the poled reference frame in contrast to the unpoled frame; this is relevant to MFCs as they are, by design, poled during production to actuate at low electric fields [49].
∆ε i e has been extended from the works of Klinkel [40] to be implemented to MFCs as, where ε s is the irreversible saturation strain with respect to the unpoled system, and the power constant b is an empirical factor that must be determined.The column vector term in (8) relates the empirical irreversible polarization influence on the irreversible saturation strain vector, where the first, second, and third row represent the in-plane longitudinal normal, transverse normal, and shear strains of the system, given plane stress assumptions.It has been empirically observed that the transverse normal strain is proportional to the longitudinal normal strain by −1/2 [37], thus the present ratio.Actuated piezoelectric systems, like MFCs, typically do not present any form of shear coupling, hence the placement of zero contributions to the irreversible saturation strain.The term in (8) can be substituted into (7) to produce the updated strain vector in the (k + 1)th iteration.However, it is necessary to mend a disconnect in (7) from the real behavior of MFCs.When domains approach the maximum polarization, they align in a way that the entire crystal population is oriented towards the direction of the applied field.As the domain polarizations saturate, the material's strain rate starts to decrease in relation to the increasing electric field.However, the second and fourth terms on the right-hand side of equation (7) do not properly account for this saturation behavior.
To rectify this deviation from the actual physics of the system, it is proposed to introduce saturation constant terms, denoted as, Here, E sat is the electric field associated with the maximum polarization.ξ i , m i , and c i are constants that must be determined empirically.The saturation parameter H i is bounded between zero and unity and is meant to diminish the strain of the system as the electric field approaches E sat .A unique value of the parameter H i is multiplied with ( 6) and ( 8) then substituted into (7) to yield the updated constitutive equation set, where the numerical indexes (1, 2) notify the unique value of H for the particular term.The logic behind using different values of H is best exemplified in (10c).The second term on the righthand side would increase strain linearly as the electric field increases near saturation; however, the fourth term increases strain quadratically thus requiring two different evolving irreversible polarizations, and by extension two different saturation parameters, to better fit the hysteretic curve.

Classical laminate plate theory with hysteresis
One can adapt the viscoelastic constitutive model in (10) to the general plate theory approach in order to mimic the hysteretic behavior of thin MFC laminates.By employing a plane stress assumption, i.e. neglecting the through-the-thickness stress terms, and rearranging the stress vector to the left-hand side, the stress equivalent of the viscoelastic model is acquired.Indeed, the terms in equation (10) can be rewritten to the constitutive piezoelectric stress equation: where Q is known as the reduced stiffness matrix, and e is the piezoelectric stress coefficient matrix.The • indicator defines the plane stress assumption applied to the term.The indices (1, 2) are introduced onto the vector nomenclature to denote the linear margination from the constitutive equation.However, this restricts (11) to the directionality of individual plies at the lamina level.
With the purpose of accounting for the assortment of orientations from each orthotropic layer, ( 12) is transformed from the lamina to global level: The relations between the lamina and global vectors are presented as, where T ε and T are the second-order strain transformation matrix and first-order transformation matrix about the transverse normal axis.Transformation can be expressed as a function of ply angle θ, The remaining course of the derivation is to substitute the mandatory matrix elements into (11), transform into (13) per lamina, then integrate the stresses with respect to the laminate thickness.This process yields the finalized constitutive equation for laminates using CLPT: or for sake of brevity, Here, N and M are the mechanical force and moment resultants, N T and M T are their imaginary thermal cousins, and N P and M P are their imaginary piezoelectric cousins.[A], [B], and [D] are the extensional, coupling, and bending stiffness matrices that summarize the individual mechanical ply properties.As their titles imply, the [A] and [D] stiffnesses are respective to their application while [B] couples the bending and twisting stiffnesses of a laminate.Additionally, the kinematics of this coupling is congregated by M xy and κ xy , which describe the twisting moment and curvature.The full expression of each stiffness matrix and piezoelectric force resultants is seen as, where the denotation k identifies the top/bottom ply layer of the current iteration, and z indicates the through-the-thickness axial position with respect to the midplane.
In order to accurately simulate the displacements of an MFC laminate, an FE approach is used.The FE method is a popular choice for analytical models as the complexities of geometry and accounted boundary conditions become trivial by discretizing the model into several elements.By extending the CLPT model to an FE model, any configuration of MFC laminates at any given boundary condition can be available at the reader's interest.In this study, nonconforming rectangular bending elements will be utilized to estimate the displacements of various cantilevered MFC laminate configurations [50].The foundational derivation process follows a similar procedure as outlined by Reddy [51].
The MFC laminate model must be benchmarked against empirical data to ensure accuracy, which is one of many options for validation assurance.This comparison focuses on matching the displacement hysteresis profiles and the magnitudes of bend/twist.Both the experimental and simulated MFC laminates have dimensions of 103 mm × 28 mm, which replicate the dimensions of the MFC M8528-P1 samples [52] used in the experiments.The MFC laminate is fixed using a cantilevered setup.
The MFC properties employed for linear and nonlinear properties at temperatures above and below 40 • C are displayed in table 1.Three switching surfaces are used to emulate the hysteretic behavior of the MFCs.The thermal and piezoelectric constants were acquired through experimental methods, as presented in the appendix.Linear elastic constants were achieved using micromechanical and equivalent single layer plate models [53].Nonlinear parameters, aside from the saturation strain, were obtained through a Sequential Quadratic Programming optimization approach, where the root mean square error (RMSE) was minimized between the simulated and observed hysteresis profiles.Note that the MFCs will be modeled as a single lamina, and its properties effectively consider the pyro-electro-mechanical contributions of its copper interdigitated electrode, Kapton, acrylic, epoxy, and PZT-5A fiber components.Detailed description of this method can be found in [54].
It is worth noting that the piezoelectric strain constants (specifically, d 11 and d 12 ) were determined by averaging the strain-electric field slopes ranging from 0 MV m −1 to the maximum electric field 2.8 MV m −1 .As a consequence of this approach, the resultant values exhibit an elevation compared to those reported in previous literature [54,55] due to contributions from both low (linear) and high (nonlinear) electric field ferroelectric behavior.Importantly, this enhancement aligns more accurately with the viscoelastic behavior of the MFC, as presented in the appendix, particularly during the initial transition from its neutral to actuated electric field position.
Furthermore, beyond the aforementioned statement, the quantity of samples employed for approximating material properties proved insufficient in achieving a normal (Gaussian) distribution of data.It becomes apparent that the furnished material properties do not necessarily represent the expected mean and variance of the overall population.Consequently, the values presented in table 1 remain subject to modification pending the accumulation of additional data points to expand the dataset.
A converged FE model of an MFC laminate is generated using 120 elements (11 × 12 elements) and 143 nodes.This model is exclusively developed in Matlab utilizing custom function files and is available in a publicly available repository [56].The MFC, together with its substrate, is considered as an equivalent single layer based on CLPT, and it incorporates the viscoelastic behavior using the multi-surface switching method.The model employs nonconforming four-node rectangular elements with Lagrange and Hermite interpolation functions to calculate in-plane (u, v) and out-of-plane displacements and slopes (w, dw/dx, dw/dy).A vectorized electric field is applied to the elements, generating an electric body force that simulates the actuated MFC laminate; thermal expansion is modeled in a similar manner except that the thermal load, i.e. change in temperature, is applied as a scalar field.For a more comprehensive discussion of the derivation details, the authors will refer to [51], where it becomes evident how to integrate the viscoelastic model into FE methods.To ascertain the fidelity of the FE model to CLPT, midplane strains and curvatures were assessed using both approaches for an arbitrary cantilevered MFC laminate.The outcome produced a consistent set of strains and curvatures between the two methodologies.
Two MFC laminate configurations are observed in the benchmarked study to validate the accuracy of the viscoelastic FE model: one using a [0 MFC /0 HR40/321 /0 MFC ] T bending bimorph and the other with a [0 MFC /45 HR40/321 ] T unimorph; bimorphs refer to a configuration with two MFCs sandwiching a substrate whereas a unimorph consists of only one MFC and its substrate.The hysteretic displacement profiles of both the FE model and experimental data for a bending bimorph configuration are presented in figure 3. A stereovision digital image correlation (DIC) method was used to measure tip displacement and twist in both actuated MFC laminate configurations.DIC is a non-contact strain-displacement measurement technique that correlates speckle pattern movement from an initial reference image to the next image of interest, making it suitable for measuring both configurations with minimal adjustments to the test bed [57].The MFC laminates were manufactured by placing MFC M8528-P1 samples on HR40/321 prepreg and then subjecting them to vacuum bagging and oven curing.The asymmetric unimorph was created by attaching the MFC to HR40/321 carbon fiber after it had cured and reached room temperature.
The profiles show points of relatively good agreement, particularly in terms of the maximum and remanent displacement values, which are comparable between the two.However, when comparing the travel profiles, there are noticeable differences in the hysteresis path.The overall error between the two loops exhibits comparable magnitudes, with a RMSE of 2.11 mm.While the estimated loop predicts the general behavior moderately well, there is room for improvement to better reflect the real behavior.For this study, the estimate of the hysteretic loop is acceptable for the substrate design studies later in this work.
Next, the comparison between the experimental data and FE model are made for a [0 MFC /45 HR40/321 ] T unimorph.The transverse bending displacement is measured along the span of the laminate, while the twist is compared at the tip of the laminates.Figure 4 presents the results of the comparison between the two studies.Good agreement is observed, with a relative error of less than 2% for the transverse bending displacement.Thus, the MFC laminate model adequately predicts the bend/twist real behavior.

Results and discussion: screening significant factors
A factorial design is employed to streamline the design process and identify the key factors that significantly affect the kinematics of a specific MFC laminate configuration.This entails perturbation of the material properties of the substrate and the variance in response with property design interactions.The responses of interest encompass four distinct outputs: tip displacement at the midpoint of the plate attributed to bending, twist at the tip, and their corresponding hysteretic counterparts.Henceforth, these will be denoted as the bending, twist, remanent bending, and remanent twist responses.The bending and twist responses are assessed subsequent to actuating the MFC at its maximum actuation input, i.e. 1500 V. Twist is computed by determining the arctangent of the difference between the out-of-plane displacement at the vertices, divided by the width of the plate.The remanent responses are collected in a similar manner, with the sole difference being the release of the electric load, allowing for the retention of remanent displacement, attributed to persistent irreversible polarization.Each lamina is assumed to deform within their elastic regime.
The influence of factors is analyzed through an Analysis of Variance (ANOVA) table.Significant factors are determined by achieving p-values, or probability values, below the predetermined significance level and rejecting the null hypothesis in (18), which states that the mean response of all substrate design perturbation treatments are indistinguishable, and is given as, The ranking of factor significance is estimated through post hoc analysis and presented graphically to facilitate consideration in the design entry process.For further details, the author will defer to primary literature [58,59].
We apply the factorial design approach to investigate different substrate material designs for an arbitrary MFC laminate configuration and identify the significant factors and interactions.To construct the factorial design of experiments, we need to identify several material properties as independent factors and assign predetermined low and high levels to measure their response.Table 2 showcases each independent factor, their typical minimum and maximum, and the corresponding high and low levels selected for the factorial design of experiments in relation to the arbitrary substrate per MFC laminate configuration.These ranges were extracted from conventional elastic materials, such as ultra-high modulus fibers, standard fibers, aluminum, steel, and wood, from literary sources [51,53] and online databases [60].
Seven independent factors are considered to vary the material properties of the substrate.These properties were selected based on the relevant independent elastic and thermal constants in (11).Substrate ply angle and thickness are also considered as their design falls under the purview of free choice; ply angle, specifically, modifies the effective laminate properties at off-axis angles.It is worth noting that the low and high levels of the ply angle are chosen as 0 • and 30 • , respectively.If angles such as 45 • or 90 • were chosen, the symmetry of a textile composite substrate would not yield a coupling kinematic response, which defeats the purpose of the need for a high-level factor.
Concurrently with the examination of significant responses resulting from variations in material properties, we also consider the importance of interactions among these material properties, specifically, how changes in one property correlate with changes in another.Notably, this study omits the consideration of higher-order interactions, such as the response to perturbations of three properties.The rationale behind this omission lies in the challenge of ascribing causality to higher-order interactions, rendering them less conclusive in the context of the factorial design ANOVA.A notable advantage of excluding these higher-order interactions is the infusion of additional degrees of freedom for error in the analysis.This, in turn, enhances the replicability of the analysis and contributes to the overall accuracy of the results [58].
The current study necessitates 2 7 = 128 runs to complete the design resolution, which proves infeasible within the limitations of the typical operating system.To overcome this challenge, a fractional factorial design is implemented to reduce the design space.In this study, a 1 /4 fractional design is selected, resulting in a reduced number of runs to 32.Note that the resolution of the reduced design space is ensured to produce identical results as the full design space, with improvements only to the computation time being made from the reduction.The suggested perturbations for each run, after partitioning the factorial design space, were generated using R-statistics, while Matlab was employed to automate and input each perturbed case into the finite element program.
To isolate the extreme bend/twist responses at the tip of the laminate, the models for each configuration closely adhere to the cantilevered MFC laminate FE model.In thermal studies conducted at varying temperatures, the influence of substrate CTEs on the maximum potential kinematic response becomes apparent.However, the effects of substrate thermal expansion overshadow the remanent displacement response at the relief stage.Hence, each iteration of the FE model factors out the thermal warping effects to observe the pure remanent behavior of the MFC laminate.
A significance level of 0.05 is chosen for this study, though any other level is viable so long as the ANOVA is updated accordingly.A total of 25 main and interaction effects are observed in each study.The following sections present an observation of the factorial design analysis concerning three common configurations of MFCs: the pure bending bimorph, the pure twist bimorph, and the unimorph.

Bending bimorph substrate significant properties
A set of simulated experiments at a temperature of 25 • C is carried out using a bending bimorph configuration to examine how variations in substrate materials affect the kinematic response.In this setup, two parallel-oriented MFCs sandwich the design substrate, and the configuration is initially maintained at its maximum actuation state, with a voltage of 1500 V applied to the top MFC and −500 V to the other [52], to measure the maximum bend/twist displacements.Additionally, the remanent displacements are assessed after the electric loads are released, and the variations in these responses are recorded in the ANOVA table; for a pure bending response, the results are presented in table 3.Each parameter and interaction are assigned a single degree of freedom because the factorial design analysis encompasses only a minimum and maximum for each.Degrees of freedom related to higher-order interactions are attributed to the residual error degrees of freedom.The sum squares/mean squares value quantifies the deviation of each parameter from the mean, while the F-statistic, along with the corresponding probability value (Pr), assesses the significance of each property or interaction.The term η 2 gauges the fractional contribution to the output displacement responses and is subsequently employed to visually represent the significance of parameters for each MFC laminate configuration and its corresponding response.Note that only the ANOVA in table 3 is presented in this document for the sake of brevity.Readers should familiarize themselves with standard ANOVA analysis [58] and apply them to this study to attain the remaining results for full resolution of the data.The remaining statistical data can be found in a publicly available data repository [56].For this study, as well as other simultaneous studies, the residual error degrees of freedom amount to six.Given that the higher-order interactions are deemed insignificant, they are omitted to introduce more repetition to the ANOVA analysis, which is an advantage of factorial design approaches.To enhance clarity and minimize repetition, it is possible to disregard a larger number of insignificant factors in the postanalysis.However, the author's intention in excluding only the higher-order interaction terms is to provide a comprehensive representation of the relevant data resolution and allow the reader to determine if further reduction of relevant design variables, i.e. increasing the repetition of the run space, is appropriate.
To ascertain the impact of each significant effect on the response, a post hoc analysis utilizing η 2 is conducted to rank their influence [58].It should be noted that the conclusions presented in this work may vary for other applications depending on the chosen significance level.As an aside, the common nomenclature to present factorial design results include the degrees of freedoms used to determine the F-statistic, followed by the p-value, and is used ordinarily in this work.
Figure 5 presents the graphical representation of the η 2 statistics, indicating the influence of various effects on their corresponding observed output kinematics.The semicolon operator ':' indicates the interaction between two properties.Notably, the significant material properties and interactions differ depending on the desired type of response.In the case of a bending bimorph, the thickness of the substrate demonstrates a significant impact on the bending deflection response (F(1,6) = 1.09 × 10 5 , p = 5.04 × 10 −14 ).The η 2 analysis reveals a 0.992 valuation for thickness; that translates to approximately 99.2% of the bending displacement variation can be attributed to perturbations in thickness.Similarly, the longitudinal modulus (F(1,6) = 819, p = 1.10 × 10 −7 ) and ply angle (F(1,6) = 22.0, p = 0.003) exhibit significant effects on bending variation, although they contribute to less than 1% to the overall effect.Considering the significance of ply angle, it is logical that anisotropy (F(1,6) = 22.7, p = 0.004) and its interaction with ply angle (F(1,6) = 21.0,p = 0.004) also influence the bending displacement of a bimorph configuration.However, their descriptive effects on bending displacement are also below 1%.Consequently, from a design perspective, allocating resources to optimizing the thickness becomes crucial as it has the most substantial influence on the bending displacement of a bending bimorph.While the thickness alone does not have a significant effect, its interaction with anisotropy (F(1,6) = 6.00, p = 0.050) is found to be reasonably significant.The interaction between anisotropy and thickness perturbations in the substrate design accounts for 8.2% of the observed twisting variation.
The remanent bend and twist responses exhibit similar pvalues and η 2 values compared to their corresponding actuated counterparts, showing minimal differences in the observed numbers.Furthermore, all significant material properties and interactions are found to be shared between these two types of responses.This indicates that the hysteresis behavior of a bending bimorph is influenced by the same factors that affect the actuated kinematic response.In other words, it suggests that the hysteresis is not solely determined by the substrate properties, but rather follows the actuated response, making it dependent on the functionality of the MFCs.This trend is consistently observed across other configurations as well.

Twisting bimorph substrate significant properties
An ANOVA analysis is performed on a twisting bimorph configuration where the top and bottom MFCs are oriented at ±45 • , respectively.Each is actuated at 1500 V at 25 • C, and the remanent displacements are measured upon the release of the electric load.Figure 6 presents the η 2 values of each factor in the form of pie charts, illustrating the influence of each factor on the corresponding kinematic response.Notably, the results obtained for significant substrate properties in a bending bimorph configuration differ significantly from the trend observed here.
Regarding the twist response, most of the elastic material properties play a role.Specifically, the longitudinal modulus (F(1,6) = 2955, p = 3.35 × 10 −5 ), anisotropy (F(1,6) = 780, p = 4.35 × 10 −5 ), and the transverse stiffness ratio (F(1,6) = 109, p = 1.35 × 10 −4 ) influence the twisting response, along with the ply angle (F(1,6) = 14.2, p = 0.009).However, their combined contribution is less than 3% of the overall twisting response.Interestingly, the shear modulus, despite its intuitive relevance due to the transverse ratio, only influences 0.08% of the twisting response.Its impact is secondary to the longitudinal modulus (2%) and anisotropy (0.52%), but it takes precedence over the ply angle (0.01%).Similar to the bending bimorph configuration, the remanent displacements of the twisting configuration exhibit the same significant properties and corresponding p-/η 2 values as the actuation displacements.Once again, this suggests that the hysteresis response is not primarily influenced by the substrate design, but rather by the internal variables of the MFC itself.

Unimorph substrate significant properties
Analogous to the ANOVA analyses conducted on the bending and twisting bimorph configurations, the unimorph configurations, being a single MFC laminate, are subjected to an ANOVA at 25 • C to ascertain the significant substrate material properties influencing the response kinematics.The unimorph actuator is oriented at 0 • and activated with a voltage of 1500 V.The η 2 statistical responses to the corresponding kinematics for each substrate material property are depicted in figure 7.
In the study of unimorph actuation response, it is observed that only the bending response is significantly influenced by a particular property, namely the thickness (F(1,6) = 762, p = 3 × 10 −5 ), along with its interaction with the longitudinal modulus (F(1,6) = 12.375, p = 0.013).On the other hand, the twisting response does exhibit properties and interactions that contribute to approximately 10% of the response, but they do not attain statistical significance under the hypothesis test; it is important to emphasize that the conclusions drawn from these tests can be subjectively adjusted at the human level, but the final determination is firmly established through rigorous coding when a property surpasses the significance level.Once again, the remanent displacements of the unimorph configurations share the same significant properties as their actuation counterparts.Consequently, similar conclusions can be drawn from previous studies on different configurations.

Thermal impact on significant substrate properties
To investigate the thermal influence on the significant substrate properties of MFC laminate configurations, the aforementioned configurations were subjected to the same factorial design analysis, with the additional variation that the temperature was increased from 25 • C to 60 • C. The selection of 60 • C as the measurement temperature was arbitrary.It was observed in preliminary analysis that the significant factors, specific to each configuration and kinematic response, vary with temperature.Therefore, in this subsection, the discussion of numerical values will be omitted as these values are specific to the given temperature.Consequently, the results will be discussed in a more general manner, focusing on trends rather than specific values.Based on the methodology described above for factorial design, the reader should be able to conduct similar tests of their own to determine the significant properties at various temperatures within the operational range of MFCs.
In the case of the bending bimorph configuration at 60 • C, it was observed that there were no variations to significant properties for any kinematic response compared to those at 25 • C (refer to figure 8).Therefore, it can be concluded that the design of the substrate should be based on the results obtained  at a controlled temperature, as the results remain consistent for any environment within the temperature range of the MFCs.
However, different conclusions can be drawn for the twisting bimorph and unimorph configurations, as depicted in figures 9 and 10.The p-and η 2 values exhibit significant variations between the two temperature ranges.Furthermore, properties that were deemed significant at the control temperature of 25 • C no longer exhibit significance at elevated or  lower temperatures.Therefore, it can be concluded that the significant substrate properties for a twisting bimorph and unimorph configuration are influenced by temperature.
An intriguing discovery is worth noting concerning the remanent displacements in both twisting bimorph and unimorph configurations with respect to the substrate proper-ties.As depicted in figures 9 and 10, the remanent bend and twist responses exhibit the same influential properties as their counterparts at the control temperature.In light of this observation, a plausible conclusion can be drawn, suggesting that the hysteretic internal variable of the MFCs is resistant to the factors that influence both the substrate and kinematics.Consequently, novel evidence is presented, indicating that the hysteretic response of an MFC laminate, excluding the bending bimorph configuration, solely responds to substrate designs in relation to the pure elastic response, without being affected by thermal effects.

Graphical design space example: bending bimorph
The previous sections have addressed the objectives of this work, and for readers interested in applying this information, this section serves as a practical guide for conducting substrate design optimization analysis.For the sake of brevity, only an example using a bending bimorph MFC laminate configuration at 25 • C is discussed.To facilitate the design process, ANOVA and η 2 analysis are utilized to guide decisionmaking.A graphical approach will be adopted to visualize the optimal points of desired response kinematics based on several significant substrate variables; it should be noted that it is always good practice to pair statistical results, i.e.ANOVA, with visual trends before presenting a definitive design choice.The desired outcomes of these examples are subjective and dependent on the reader's preferences.This includes the selection of initial constants to hold all other substrate properties as the interested properties are varied in the designated design space, as presented in table 4.
It is important to prioritize allocating resources and efforts based on the highest η 2 values.If the desired response is bending, then the substrate's thickness and longitudinal modulus should be the primary focus of design, as they contribute to 99.2% and 0.75% of the variation in bending displacement response, respectively.These contributions outweigh those of other material properties or interactions.When two significant properties demonstrate significant interaction effects on the kinematic response, it is advantageous to consider both properties simultaneously to identify the optimal design within the given design space, if feasible within the available resources.If there are no interaction effects, evaluating the design space between the two properties is optional.Figure 11 provides a visualization of the actuation and remanent bending displacements in relation to varying thickness and longitudinal modulus values.This graph can aid in identifying the optimal design points within the design space for the desired bending response.
Figure 11 clearly illustrates that reducing the thickness of the substrate leads to higher bending displacements if that is the desired outcome.Since thickness has the higher η 2 values, its effect on the bending displacement is more pronounced than the longitudinal modulus.However, it is important to consider the tradeoff with hysteresis.The significant substrate properties that influence the bending response also affect the remanent displacement.Therefore, maximizing bending displacement will correspond to an increase in remanent displacement, and vice versa.Thus, the trend observed in figure 11 highlights a tradeoff scenario.
Once the initial set of significant properties has been determined, the subsequent focus shifts to studying the next set of significant material properties.In the case of bending displacement for a bending bimorph, these properties include substrate anisotropy and ply angle, whose interaction contributes to 0.018% of the bending displacement variation.It is advantageous to generate a contour plot showcasing the influence of these two properties on the bending response.However, the interaction between anisotropy and ply angle not only impacts the bending kinematic response but also exhibits statistically significant influence on the twisting response variation, as confirmed by the ANOVA analysis.When confronted with scenarios where both bending and twisting responses share similar significant properties or interactions, it is advisable to compare the respective response curves or contours.This approach facilitates the identification of regions that optimize both responses or prioritize one while mitigating the other.
Figure 12 presents the bend/twist contour plots illustrating the relationship between anisotropy and ply angle for a bending bimorph.Optimality can be achieved by comparing design points between plots.For instance, if the goal is to achieve maximum bending displacement with minimal twist, selecting a substrate with nearly unidirectional anisotropy and a ply angle of 90 • would be the optimal choice.However, for a well-balanced tradeoff between the two responses, it is advisable to opt for a similarly anisotropic substrate with a ply angle approximately ±30 • .It is observed that a similar trend is applicable to remanent displacement behavior.The ultimate selection depends on the reader's specific requirements and preferences.Nonetheless, the graphical methods provided offer valuable insights for making an informed decision.
Sometimes the optimal design point, as perceived by the reader, may fall within a region that requires fictitious material  properties.For instance, it could be found that an anisotropy value of 0.6 best aligns with the optimal performance of an MFC laminate, even though such a material does not exist naturally or commercially.Despite the absence of readily available materials with these specific properties, it does not preclude the possibility of manufacturing them.In this scenario, achieving an anisotropy of 0.6 could be accomplished by creating a fibered composite with variations in the fiber-volume matrix to achieve said anisotropy; alternatively, a laminate that is constructed to attain an anisotropy of 0.6 is also feasible.

Conclusion
This concludes the analysis and compilation of significant material properties that impact the displacement response of morphing MFC laminates.The identification of these factors employed a randomized block design strategy, utilizing a 2 k fractional factorial design of experiments method to linearize the process within an infinite design space.
This work offers valuable insights into the identification of significant properties that influence the actuation kinematics of MFC laminates, presenting a comprehensive guide for thin actuator technology beyond piezoceramic fibered actuators.It provides a step-by-step process for building ANOVA data and offers end conclusions that can serve as a design reference.
In conclusion, the following are accomplishments or findings in this work that can contribute to or motivate prospective research: • A set of ANOVA response data is published for the reader to use for entry into design of thin piezoelectric, or other applicable smart material, laminates.• The hysteretic behavior of an MFC, or any applicable smart actuation, laminate is proportional to the maximum actuation kinematics of the laminate.• The hysteretic behavior is only a function of elastic properties and is not influenced by thermal perturbations in the substrate.

Prospective work
For quick and efficient reference in the design of piezoelectric or similar smart material actuators, the present work provides valuable insights.However, there is room for further exploration within the present research space.For instance, the findings of this study have the potential to extend into selfconsistency studies, such as artificial intelligence learning or machine learning.An intriguing research area to explore is the plastic response of the substrate.While it has been established that MFCs do not undergo plastic deformation, particularly in tension, when fully actuated, the surface traction exerted on the substrate can induce stresses exceeding 100 MPa for a carbon fiber-epoxy substrate (which would not yield).This raises the possibility of substrate plasticity interacting with the piezoelectric actuation response and its inherent hysteresis.Consequently, there is a need to investigate how plastically deformed substrates can influence the resulting kinematics.Moreover, by incorporating mechanical loads alongside thermal and electrical considerations, we can gain a more comprehensive understanding of the substrate's influence on MFC polarization behavior, enhancing the overall analysis.The author also encourages the implementation of nonlinear ferroelastic actuator behavior, where applicable, specifically with the stress resultants that accumulate in each ply of an MFC-substrate laminate.
The experimental investigation and modeling of MFCs as actuators was to specifically explore their converse effect as piezoelectric actuators operating at the quasi-static state level.The dynamic behavior of piezoelectric materials exhibits distinct characteristics [1], which falls outside the scope of this study but can be expanded upon in prospective research.Additionally, the perturbation of substrate material design factors did not encompass geometric complexities or variations in boundary conditions.The author's rationale for this approach is that such conditions are problem-specific and do not fall under the purview of free design choice.This reasoning extends to mechanical loading scenarios as well.However, the basis of these studies has been established in this work and can be readily adapted to a more dynamic environment if the reader so chooses.The collected linear thermal and piezoelectric coefficients are presented in figures A2 and A3.A minimum of eight MFC type M8528-P1 samples were used for data acquisition.Nonetheless, the current sample size did not yield a normal distribution.Consequently, box-and-whisker plots have been chosen as the preferred method for visualizing the distribution of data points.It is important to note that these data have undergone post-processing to eliminate the influence of outliers and significant noise, ensuring the selection of appropriate properties.A trendline is presented between the means of each data subset but should not be interpreted as the actual trend of MFC temperature-dependent behavior until further data is collected.
The behavior of MFC minor loops is examined at various temperatures in figure A4.Understanding the minor loop behavior under different input parameters in tandem with isothermal environments is crucial for improving the accuracy of the subsequent model.The sinusoidal input journey of the minor loop starts from the neutral point, reaches the highest allowable voltage of 1500 V, then goes down to the lowest allowable voltage of −500 V, and finally returns back up to 1500 V. Figure A4 illustrates the average minor loop of MFCs at a 5 Hz actuation rate, showcasing the longitudinal normal strains at different temperatures.

Figure 1 .
Figure 1.Image of macro-fiber composite and its modular material components.

Figure 2 .
Figure 2. Illustration of multi-surface switching method.The domains of a localized piezoelectric crystal structure are color coded with their surface counterpart.

Figure 5 .
Figure 5. Pie chart of eta squared values for 25 main and interaction effects for a given bending bimorph [0 MFC /θ Substrate /0 MFC ] T configuration at 25 • C. The chart presents the kinematic response for tip (a) bending deflection, (b) twist, (c) remanent bending, and (d) remanent twist.

Figure 6 .
Figure 6.Pie chart of eta squared values for 25 main and interaction effects for a given twisting bimorph [45 MFC /θ Substrate /-45 MFC ] T configuration at 25 • C. The chart presents the kinematic response for tip (a) bending deflection, (b) twist, (c) remanent bending, and (d) remanent twist.

Figure 7 .
Figure 7. Pie chart of eta squared values for 25 main and interaction effects for a given unimorph [0 MFC /θ Substrate ] T configuration at 25 • C. The chart presents the kinematic response for tip (a) bending deflection, (b) twist, (c) remanent bending, and (d) remanent twist.

Figure 8 .
Figure 8. Pie chart of eta squared values for 25 main and interaction effects for a given bending bimorph [0 MFC /θ Substrate /0 MFC ] T configuration at 60 • C. The chart presents the kinematic response for tip (a) bending deflection, (b) twist, (c) remanent bending, and (d) remanent twist.

Figure 9 .
Figure 9. Pie chart of eta squared values for 25 main and interaction effects for a given twisting bimorph [45 MFC /θ Substrate /-45 MFC ] T configuration at 60 • C. The chart presents the kinematic response for tip (a) bending deflection, (b) twist, (c) remanent bending, and (d) remanent twist.

Figure 10 .
Figure 10.Pie chart of eta squared values for 25 main and interaction effects for a given unimorph [0 MFC /θ Substrate ] T configuration at 60 • C. The chart presents the kinematic response for tip (a) bending deflection, (b) twist, (c) remanent bending, and (d) remanent twist.

Figure 11 .
Figure 11.(a) Actuation and (b) remanent bending displacements of bending bimorph with respect to substrate thickness and longitudinal modulus at room temperature.

Figure 12 .
Figure 12.Actuation and remanent bend/twist overlayed contour plots with respect to anisotropy and ply angle for a bending bimorph configuration at room temperature.The plots are organized in a chronological row pattern: (a) actuation bending displacement, (b) remanent bending displacement, (c) actuation twist, and (d) remanent twist.

Figure A1 .
Figure A1.Fabricated test bed for thermal dependent MFC testing.An additional DIC camera is placed equidistant on the opposite side of the Rohacell chamber.

Figure A2 .
Figure A2.Linear piezoelectric strain constants for longitudinal (d 11 ), transverse (d 22 ), and shear (d 16 ) directions in box-and-whisker graphical format.The intercepting line represents the mean line from each temperature range.Outliers have been removed.

Figure A3 .
Figure A3.Thermal expansion box-and-whiskers plot for MFC specimens in the longitudinal (α 11 ), transverse (α 22 ), and shear (α 22 ) directions.The intercepting line represents the mean line from each temperature range.

Figure A4 .
Figure A4.Minor loop plots for MFC actuation strain in the longitudinal normal (top), transverse normal (bottom left), and shear direction (bottom right).Sample actuation rates at 5 Hz.Color coding is used to assign the profile to the correlated temperature.
12 A 16 A 12 A 22 A 26 A 16 A 26 A 66 B 11 B 12 B 16 B 12 B 22 B 26 B 16 B 26 B 66 B 11 B 12 B 16 B 12 B 22 B 26 B 16 B 26 B 66 D 11 D 12 D 16 D 12 D 22 D 26 D 16 D 26 D 66

Table 1 .
Recommended nonlinear variables to drive multi-surface method.

Table 2 .
Substrate property ranges for available linear elastic materials.

Table 4 .
Sample initial substrate material properties used for substrate design optimization example.