Active layered composites: A variable kinematics approach with stimulus expansion model for piezo-actuators and hydrogels

Abstract

Active materials are often applied in the form of plate-like smart composites. Therein, they are combined with other (active or passive) materials that serve different purposes, such as providing mechanical counterpoints, protection against multi-field influences, or electrical functionality. Furthermore, sensoric layers can be included. For an enhanced description of the mechanical response of Soft-Hard Active-Passive Embedded Structures (SHAPES) in plate configuration, we combine two approaches: (i) A variable kinematics approach called sublaminate Generalized Unified Formulation (sGUF), which allows the choice of adequate kinematics for each physical layer, and (ii) the Stimulus-Expansion-Model (SEM) for description of active behavior. We present results for SHAPES with polymer gel layers and for piezo-ceramic layers. The current sGUF-SEM modeling approach has been protypically implemented along with a Navier-type solution. This allows the assessment of combinations of different passive and active materials inside plate-like composites. Excellent agreement with reference solutions from literature and/or from three-dimensional Finite Element simulations in commercial software tools has been found. The promising preliminary results presented in this work suggest further steps to be considered: the implementation of sGUF-SEM into more general numerical solution methods, such as Finite Element or Ritz methods, and the integration of the physics of additional active materials.

1. Introduction: Active laminates made from soft and hard layers

Smart structures consisting of active and passive layers can be applied for various purposes, such as function-integration, vibration control, etc. (Bar-Cohen, 2000; Chopra and Sirohi, 2013; Ehrenhofer, 2023a; Gaudenzi, 2009; Leo, 2007). The wide range of the stimulus-responsive (active) materials in these structures goes from hard ceramic materials that react to electrical fields to soft polymeric materials that respond to a variety of stimuli (often at the same time, i.e. multisensitivity (Ehrenhofer et al., 2020a)). Integrated into structures, they provide a complex composite behavior. Efficient structural models can be applied for these structures to allow assessing and optimizing performance, incorporating the various design parameters that typically enter the definition of such complex structures. Beyond the conventional parameters related to the macroscale response (dictated by geometry and boundary conditions), also mesoscale and microscale parameters can be of paramount importance. Parameters can be, for example, the stacking sequence of the active and passive layers, stiffness mismatches, the joining between the layers, spatially distributed activation stimuli, or the chemo-physical constitution of the smart material. Currently, classical macro-scale models available in commercial simulation software are inadequate in accounting for such lower scale effects. Furthermore, the solution of the 3D fully resolved continuum mechanics equations—including the stimulus-to-mechanics coupling that is typical of smart materials—makes an efficient exploration of the design space prohibitively expensive.

These challenges get magnified when layers with different activities are included, see Figure 1. Furthermore, when passive layers are integrated that power the active layers, such as accumulators (Winkler et al., 2020), the local stresses and strains in these layers must be well understood. To meet these challenges, we propose a variable kinematics approach for the definition of structural models (beams, plates, shells) of composites with arbitrary compositions of active and passive layers. The active behavior is included with a physics-based regression method called Stimulus-Expansion-Model (SEM) (Ehrenhofer and Wallmersperger, 2018; Ehrenhofer et al., 2018), allowing the definition of homogenized active properties of the smart material in every layer. The variable kinematics approach provides a virtually infinite number of structural models, ranging from the well known and simple Classical Laminate Theory (CLT) and First-order Shear Deformation Theory (FSDT), up to quasi-3D models, encompassing Equivalent Single Layer (ESL) as well as Layer-Wise (LW) models (Carrera, 2003; Carrera et al., 2011). This way, it is possible to adapt the computational cost of the structural model depending on the problem at hand and on the output quantity of interest, for example, a global displacement or a local stress.

Figure 1.

Example smart laminate with layers of different activity under different stimulus fields. From left to right: Electric field $E_k$ , chemical concentration field $c_{\alpha }$ of the species $\alpha$ , temperature $\vartheta$ , humidity $\varphi$ .

General Literature Review The class of composites comprising layers with different stiffness and activity can be called Soft-Hard Active-Passive Embedded Structures (SHAPES) (Ehrenhofer, 2023a, 2023c). In addition to the reaction to outer mechanical loads, their behavior heavily depends on the environment in which they are embedded, because this environment provides the stimulus for their activation/actuation (Ehrenhofer et al., 2021). Different classes of active materials can be used in SHAPES. Comprehensive overviews over active material classes and their application in structures are given in important literature works, for example, of Bar-Cohen (2000), Leo (2007), Gaudenzi (2009), or Chopra and Sirohi (2013). Since in the current work our attention will be devoted to hydrogels and piezo-ceramics, the following literature review will be restricted to these two active material classes.

Piezo-ceramics are the most classical electro-active materials with a wide variety of industrial applications (Erturk and Inman, 2011; Gao et al., 2020; Huang and Li, 2021). They are polycrystalline materials that reversibly react to electrical fields, due to their Perovskite-type lattice structure.

On the modeling side, there have been extensive efforts to include piezoelectric material behavior, in particular piezo-ceramic actuators and sensors, in structural analysis tools such as the Finite Element Method (FEM) (Benjeddou, 2000). It has been recognized early that dedicated structural theories of higher order (compared to the classical laminate models) are recommended for grasping the complex electro-mechanical response induced by the heterogeneous stack of passive substructures and active piezo-layers (Ballhause et al., 2005; Gopinathan et al., 2000; Saravanos and Heyliger, 1999; Wang and Yang, 2000). Over the last 20 years, a large number of plate/shell FE for composite structures with embedded piezo-layers have been proposed (Balamurugan and Narayanan, 2008; Kögl and Bucalem, 2005; Kulikov and Plotnikova, 2008; Kulikov et al., 2020; Lage et al., 2004; Lammering and Mesecke-Rischmann, 2003; Marinkovich et al., 2006; Vidal et al., 2011; Zemčík et al., 2007). Extensions to geometric nonlinearities have been taken into account as well (Klinkel and Wagner, 2006; Pasquali and Gaudenzi, 2016; Varelis and Saravanos, 2004). Variable kinematics-based approaches have been developed for smart piezoelectric composites in order to reduce the computational effort upon limiting the required high-order models to those local regions in which steep gradients of the response are to be expected (D’Ottavio et al., 2018; Robbins and Reddy, 1996; Zappino and Carrera, 2020).

Hydrogels are polymeric materials with the capability to reversibly swell and deswell in water, which makes them suitable for small-scale actuator and sensor purposes (Gerlach and Arndt, 2009; Gerlach et al., 2021).

Polyelectrolyte gels and other hydrogels are often modeled based on the multi-physics processes, for example, using multi-field theory (Li, 2009; Lucantonio et al., 2013; Wallmersperger et al., 2004, 2008), or poromechanical models (de Boer, 2005) like the Theory of Porous Media (Bluhm et al., 2016; Ehlers and Bluhm, 2002; Ehlers, 2002; Leichsenring and Wallmersperger, 2017). In our previous works, we focused on describing the active response of hydrogels based on a physics-based regression using the analogy concept (or: isomorphy in the describing equations) between active expansion and thermal expansion, called the Stimulus-Expansion-Model (SEM) (Ehrenhofer et al., 2018). Besides a single stimulus, hydrogels can respond to multiple stimuli with reversible swelling (Gerlach and Arndt, 2009). Concepts of including this multisensitivity were shown as well (Ehrenhofer et al., 2020a).

The present work proposes a novel approach to SHAPES, based on the combination of the variable kinematics approach, which allows an efficient representation of the heterogeneous constituent material plies, and the Stimulus Expansion Model, which allows to represent the active behavior of a large variety of smart materials. A prototypical implementation is carried out that is—in the current work—limited to geometrical linear response of simply-supported rectangular plates. This simplification nevertheless allows to illustrate the capabilities of the proposed approach and leads to various opportunities for extension, such as other active materials, geometrical and material nonlinearity, sensory behavior, functionally grading of material properties, and locally changing field properties, which will be part of future research.

Overview over the current work The current work is an extended version of our previous conference paper for the SMART 2023 conference in Patras, Greece (Ehrenhofer et al., 2023). It is structured as follows: In the following Section 2, the Stimulus-Expansion-Model (SEM) and the sublaminate Generalized Unified Formulation (sGUF) are presented. Explicit SEM formulations for piezo-ceramics and hydrogels are provided there, as well. In the results Section 3, some active-passive composites are calculated with our combined sGUF-SEM model. The results are compared to literature results and/or 3D finite element simulations. Finally, the main conclusions are summarized and an outlook towards future work is suggested in Section 4.

2. Methods: Modeling of active layers

In the current section, the Stimulus-Expansion-Model in its general form is explained, see Section 2.1. The model is then applied to two different active material classes: piezo-ceramics and hydrogels. This is followed by a description of the sublaminate Generalized Unified Formulation (sGUF) theory and its assumptions, see Section 2.2. The solution with Navier’s method is shown in Section 2.3, and results for selected testcases are discussed in the next Section 3.

2.1. Stimulus-Expansion-Model for the inclusion of active behavior

The Stimulus-Expansion-Model (SEM) is a “greybox approach” that allows the inclusion of active behavior by using the analogy between active expansion and thermal expansion (Ehrenhofer et al., 2018, 2020b). Due to the isomorphy of the continuum mechanics equations describing stimulus-active behavior and thermal expansion, the SEM can be implemented in the same way as thermal expansion, for example, in the context of Finite Element Analysis software, or in structural mechanics approaches. The temperature $\vartheta$ (as well as its change $\mathrm{\Delta }\vartheta$ ) and the expansion coefficient ${\alpha }^{\mathrm{therm}}$ are then reinterpreted (physical parameters of active behavior → thermal expansion coefficient(s) and temperature change). This transfer/reinterpretation of parameters is the key step of applying the SEM; in our previous works, we provided a software toolbox to facilitate the step of deriving the expansion curves based on the physical behavior (Ehrenhofer, 2023b). Besides monosenstive behavior (reaction to a single stimulus), the SEM can also be applied for multisensitive reactions with any number of (cross-)sensitivities (Ehrenhofer et al., 2020a).

Active behavior can be modeled linearly or nonlinearly. There are different types of nonlinearity that must be considered, depending on the scope of the regarded physical system. In the current work, we use linearity in the kinematics (small deformation gradients), and in both the active and passive material behavior, as already done for hydrogels in permeation control membranes (Ehrenhofer et al., 2016) and Dielectric Elastomer Actuators (DEAs) in antagonistic setup (Franke et al., 2020). The constituents of the model are the following:

$$\mathrm{Mechanical}\mathrm{balance}\mathrm{laws}{\sigma }_{kl,k}+f_l=0,{\sigma }_{kl}={\sigma }_{lk}$$

(1)

$$\mathrm{Stimulus}\mathrm{field}\mathrm{balance}\mathrm{law}(\mathrm{s}){\Phi }_{,ii}=\varphi$$

(2)

$$\mathrm{Kinematics}{\epsilon }_{kl}=\frac{1}{2}\left(u_{k,l}+u_{l,k}\right)$$

(3)

$$\mathrm{Material}\mathrm{law}{\sigma }_{kl}=C_{klmn}\left({\epsilon }_{mn}-\underset{\mathrm{Active}\mathrm{expansion}\mathrm{term}{\epsilon }_{\mathrm{mn}}^{\mathrm{active}}}{\underbrace{{\alpha }_{mn}\mathrm{\Delta }S}}\right)$$

(4)

Here, $\mathbf{\sigma }$ denotes the symmetric stress tensor, $\mathit{f}$ is the body force vector, $\mathbf{\epsilon }$ is the linear technical strain tensor (in reference to an initial state, $\mathit{u}$ is the displacement vector, $\mathit{C}$ is the elasticity tensor, and $\mathbf{\alpha }$ is the stimulus (strain) expansion tensor. The change of the stimulus $\mathrm{\Delta }S$ is defined according to a reference state and can take various forms. It is generally derived from the stimulus field variable. The governing field variable $\Phi$ (here a scalar) of the stimulus is field resolved, that is, its distribution is governed by a balance law equation (2) (here for an equilibrium state, leading to a Poisson’s equation) with a possible source term $\varphi$ . The concrete form of the field equation and therefore the derivation of $S=\mathcal{D}\left(\Phi \right)$ depends on the physics of the described active materials. However, this equation can also be simplified, for example, for plate electrodes on piezo-ceramics, in which the electric field can be assumed homogeneous, see Section 2.1.2.

Please note that Einstein’s summation convention holds. Isotropic expansion is represented by setting ${\alpha }_{\mathrm{mn}}={\alpha }^{\mathrm{S},\mathrm{i}}{\delta }_{\mathrm{mn}}$ , where the expansion coefficient ${\alpha }^{\mathrm{S},\mathrm{i}}$ is a scalar and ${\delta }_{\mathrm{mn}}$ denotes the Kronecker-delta. Please note that equation (4) can also be written in terms of the stress expansion coefficients ${\lambda }_{\mathrm{kl}}=C_{\mathrm{klmn}}{\alpha }_{\mathrm{mn}}$ .

Time-dependency considerations In the current work, the behavior is presented in its static form according to the following considerations: (i) no external mechanical loads are considered that introduce inertial effects (hence the right hand side of the mechanical balance law (1) is zero); (ii) the transient, time-dependent processes in the stimulus field are assumed to have already reached an equilibrium (hence, there is no time derivative in the balance law for the stimulus field (2)); (iii) any time-dependent material behavior (such as viscoelasticity) is neglected as well, as can be seen in the constitutive law (4).

Different physics for different materials The specialty of the Stimulus-Expansion-Model is the definition of the active expansion term ${\epsilon }_{\mathrm{mn}}^{\mathrm{active}}$ in the material law (4), according to the physics of the active material. This term denotes a one-way coupling from the non-mechanical field to the mechanical field through a stress-free expansion term. In the current work, two active material classes are considered. The physical phenomena lead to active expansion terms as described in the respective sections:

• Hydrogels are considered in Section 2.1.1. Various previous works of the authors uncover different effects of these active materials with nonlinear isotropic actuation behavior. The role of the reference state on the sensitivity of the active material can be well illustrated from this material class.

• Piezo-ceramics are considered in Section 2.1.2 as linear actuator materials. Owing to their crystalline micro-structure, however, this linear behavior is orthotropic.

In combination with the active materials introduced above, the present work considers composite structures comprising some passive material layers as well. For the hydrogel testcase, we chose a soft material with low elastic modulus. For the piezoelectrically actuated composite structure we consider a passive base laminate consisting of Graphite-Epoxy orthotropic plies with 0° and $90°$ orientation.

2.1.1. Stimulus Expansion Model for hydrogels

Active hydrogels are synthethic polymers that can reversibly swell and deswell in water (Gerlach and Arndt, 2009). A commonly applied hydrogel is the temperature-sensitive poly(N-isopropyl-acrylamide) (Ehrenhofer et al., 2016) (pNiPAAm), which is considered in the current work. However, many other active polymers with different sensitivities can be found in literature and be represented within the framework of the SEM (Ehrenhofer et al., 2018). The basic principle and the definition of reference and actuated state are shown in Figure 2.

Figure 2.

(a) Principle and states for hydrogels under homogeneous swelling conditions, (b) linearization of the passive material behavior at the stress-free reference point, (c) linearization of the active swelling behavior at the inflexion point of the sigmoidal swelling curve. An example behavior for a Lower Critical Solution Temperature (LCST)-type swelling behavior, as shown here, can be found, for example, in our previous work (Ehrenhofer et al., 2018).

The isotropic expansion coefficient for hydrogel (HG) swelling ${\alpha }^{\mathrm{HG}}$ is derived from free swelling experimental data, that is, the swelling curve of the hydrogel diameter $d^{\exp }\left(\vartheta \right)$ , with the (kinematically linear) technical strain

$${\epsilon }^{{\vartheta }_{\mathrm{ref}},\exp }\left(\vartheta \right)=\frac{d^{\exp }\left(\vartheta \right)-d_{\mathrm{ref}}}{d_{\mathrm{ref}}}.$$

(5)

The reference state is with reference diameter $d_{\mathrm{ref}}$ and reference temperature ${\vartheta }_{\mathrm{ref}}$ . By applying the normalization concept from our previous works (Ehrenhofer and Wallmersperger, 2018; Ehrenhofer et al., 2018), other hydrogel sensitivities can be included in the Stimulus-Expansion-Model. Then, the stimulus difference can be gained, for example, by using a difference in concentrations (or pH) (Binder et al., 2022), or light intensity (Mählich et al., 2023). Please note that the expansion coefficients can also be derived directly from physics, for example, from statistical theories, such as Flory-Rehner-Theory.

The stimulus distribution can be derived by solving the corresponding stimulus field equation (2). For thermo-responsive hydrogels with lower critical solution temperature (LCST) behavior, such as pNiPAAm, $\Phi$ is the temperature $\vartheta$ and $\varphi$ equals zero. This is the stationary form of Fourier’s heat transfer equation without heat sources/sinks ( ${\vartheta }_{,\mathrm{kk}}=0$ ). For chemical stimuli, Fick’s law of diffusion is used instead; in this case, sources and sinks can represent chemical reactions.

The final equations to be inserted in equation (4) are:

$${\alpha }_{\mathrm{mn}}^{\mathrm{HG}}={\delta }_{\mathrm{mn}}\frac{{\epsilon }^{{\vartheta }_{\mathrm{ref}},\exp }}{\vartheta -{\vartheta }_{\mathrm{ref}}};\mathrm{\Delta }{S}^{\mathrm{HG}}=\vartheta -{\vartheta }_{\mathrm{ref}}.$$

(6)

2.1.2. Stimulus Expansion Model for piezo-ceramics

Piezo-ceramics are crystalline materials that can deform due to electrical potential differences. A typical piezoelectric-material is lead zirconate titanate (PZT). Depending on the mode in which it is actuated, different effects of the anisotropic behavior are relevant. The basic principle is shown in Figure 3.

Figure 3.

(a) Principle and states for a piezo-ceramic layer with polarization in 3-direction, (b) passive orthotropic material behavior, (c) linearization of the active behavior for small electric fields (compared to the electric field, when the saturation effect leads to a significant deviation from the linear approximation)

Strictly speaking, piezoelectricity is inherently a linearized model valid for small electric fields. The equations are written according to the IEEE standard on piezoelectricity (IEEE, 1988) with ${\tilde{E}}_i$ denoting the electric field and ${\tilde{D}}_i$ the dielectric displacement vector. For further reference, we introduce the vectors of generalized stresses ${\mathbf{\sigma }}^{\mathit{g}}=[\mathbf{\sigma },\tilde{\mathit{D}}{]}^{\top }$ and generalized strains ${\mathbf{\epsilon }}^{\mathit{g}}=[\mathbf{\epsilon },\tilde{\mathit{E}}{]}^{\top }$ comprising both the mechanical and the electrical field variables. The electro-mechanically strongly coupled material behavior is expressed in Voigt notation with indices $P,Q\in \{1,6\}$ in terms of the enthalpy function $H({\sigma }_P,E_i)$ , which leads to the following expression with a symmetric constitutive matrix:

$$\begin{array}{c}\multicolumn{1}{c}{\left(\begin{array}{c}{\epsilon }_1\hfill \\ \multicolumn{1}{c}{{\epsilon }_2}\\ \multicolumn{1}{c}{{\epsilon }_3}\\ \multicolumn{1}{c}{{\epsilon }_4}\\ \multicolumn{1}{c}{{\epsilon }_5}\\ \multicolumn{1}{c}{{\epsilon }_6}\\ \multicolumn{1}{c}{{\tilde{D}}_1}\\ \multicolumn{1}{c}{{\tilde{D}}_2}\\ \multicolumn{1}{c}{{\tilde{D}}_3}\end{array}\right)=\left(\begin{array}{ccccccccc}{S}_{11}^{\mathrm{E}}\hfill & S_{12}^{\mathrm{E}}\hfill & S_{13}^{\mathrm{E}}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & {\tilde{d}}_{31}\hfill \\ \multicolumn{1}{c}{S_{12}^{\mathrm{E}}}& S_{22}^{\mathrm{E}}\hfill & S_{23}^{\mathrm{E}}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & {\tilde{d}}_{32}\hfill \\ \multicolumn{1}{c}{S_{13}^{\mathrm{E}}}& S_{23}^{\mathrm{E}}\hfill & S_{33}^{\mathrm{E}}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & {\tilde{d}}_{33}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & S_{44}^{\mathrm{E}}\hfill & 0\hfill & 0\hfill & 0\hfill & {\tilde{d}}_{24}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & S_{55}^{\mathrm{E}}\hfill & 0\hfill & {\tilde{d}}_{15}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & S_{66}^{\mathrm{E}}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & {\tilde{d}}_{15}\hfill & 0\hfill & {ϵ}_{11}^T\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & {\tilde{d}}_{24}\hfill & 0\hfill & 0\hfill & 0\hfill & {ϵ}_{22}^T\hfill & 0\hfill \\ \multicolumn{1}{c}{{\tilde{d}}_{31}}& {\tilde{d}}_{32}\hfill & {\tilde{d}}_{33}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & {ϵ}_{33}^T\hfill \end{array}\right)\left(\begin{array}{c}{\sigma }_1\hfill \\ \multicolumn{1}{c}{{\sigma }_2}\\ \multicolumn{1}{c}{{\sigma }_3}\\ \multicolumn{1}{c}{{\sigma }_4}\\ \multicolumn{1}{c}{{\sigma }_5}\\ \multicolumn{1}{c}{{\sigma }_6}\\ \multicolumn{1}{c}{{\tilde{E}}_1}\\ \multicolumn{1}{c}{{\tilde{E}}_2}\\ \multicolumn{1}{c}{{\tilde{E}}_3}\end{array}\right)}\end{array}$$

(7)

The generalized constitutive matrix is here expressed in terms of the elastic compliance coefficients $S_{\mathrm{PQ}}^E$ at constant electric field, the dielectric permittivity coefficients at constant stress ${ϵ}_{\mathrm{ij}}^{\mathrm{T}}$ and the piezoelectric coupling coefficients ${\tilde{d}}_{\mathrm{iP}}$ . The corresponding mechanical stiffness coefficients at constant electric field are $C_{\mathrm{PQ}}^{\mathrm{E}}={\left[S_{\mathrm{PQ}}^{\mathrm{E}}\right]}^{-1}$ . Note that the above piezoelectric constitutive law has been expressed in the principal material frame, of which the $x_3-$ axis is the polarization direction.

The equivalent expansion coefficients of the SEM can be derived by considering a stress-free actuation (i.e. all generalized stresses ${\mathbf{\sigma }}^{\mathit{g}}=0$ ). Assuming a driving electric field that is aligned with the polarization, that is, ${\tilde{E}}_3\ne 0$ , ${\tilde{E}}_1={\tilde{E}}_2=0$ , equation (7) simplifies to

$${\epsilon }_1={\tilde{d}}_{31}{\tilde{E}}_3{\epsilon }_2={\tilde{d}}_{32}{\tilde{E}}_3{\epsilon }_3={\tilde{d}}_{33}{\tilde{E}}_3,$$

(8)

where the terms associated to ${\tilde{d}}_{24}$ and ${\tilde{d}}_{15}$ drop out due to ${\tilde{E}}_1={\tilde{E}}_2=0$ . The electric field strength can be derived by solving the stimulus field equation (2) with $\Phi$ being the electrostatic scalar potential $\tilde{\Psi }$ , that is, following the definition ${\tilde{E}}_i=-{\tilde{\Psi }}_{,\mathrm{i}}$ ). Please note that $\varphi =-q/ϵ=0$ because piezo-ceramics do not have any net volume charge $q$ .

The equivalent expansion coefficients of the SEM are eventually obtained as

$$\begin{array}{cc}\multicolumn{1}{c}{{\alpha }_{\mathrm{xx}}^{\mathrm{PZT}}}& ={\alpha }_{11}^{\mathrm{PZT}}=\frac{{\epsilon }_{11}\left({\tilde{E}}_3\right)}{{\tilde{E}}_3}={\tilde{d}}_{31};\hfill \\ \multicolumn{1}{c}{{\alpha }_{\mathrm{yy}}^{\mathrm{PZT}}}& ={\alpha }_{22}^{\mathrm{PZT}}=\frac{{\epsilon }_{22}\left({\tilde{E}}_3\right)}{{\tilde{E}}_3}={\tilde{d}}_{32};\hfill \\ \multicolumn{1}{c}{{\alpha }_{\mathrm{zz}}^{\mathrm{PZT}}}& ={\alpha }_{33}^{\mathrm{PZT}}=\frac{{\epsilon }_{33}\left({\tilde{E}}_3\right)}{{\tilde{E}}_3}={\tilde{d}}_{33}\hfill \end{array}$$

(9)

along with

$${\alpha }_{\mathrm{mn}}=0\mathrm{for}m\ne n\mathrm{and}\mathrm{\Delta }S={\tilde{E}}_3.$$

(10)

Please note that the proposed model does not include a hysteretic behavior as reported, for example, by Goldfarb and Celanovic (1997).

2.2. Variable kinematics model in sublaminate generalized unified formulation

Reduced-order structural models (1D beam, 2D plates/shells) are formulated in axiomatic sense upon introducing a priori assumptions for the kinematic field across the composite stack. A whole family of models can be encompassed by expressing the kinematic constraints in the compact index notation referred to as Carrera Unified Formulation (CUF) (Carrera et al., 2014). The present work adopts the formal extension of CUF for plate models towards the possibility of regrouping plies into sublaminates, in order to allow a physically meaningful definition of the smart composite, in which passive sub-structures are functionally separated from the active components, which in turn may be composed of driving plies (e.g. electrodes), energy storage plies (e.g. accumulators), and active material plies. Furthermore, active-active and active-passive combinations can be grouped into sublaminates. In Figure 4, an example for the subdivision into sublaminates is provided.

Figure 4.

Definition of layers and coordinate systems of the example from Figure 1. The full details of coordinate definition can be found in our previous work (D’Ottavio et al., 2016).

A brief description of this sublaminate Generalized Unified Formulation (sGUF) is given next, for more details reference is made to the original paper (D’Ottavio, 2016) and its extension to piezoelectric composites (D’Ottavio et al., 2018).

We shall consider a plate occupying the volume $V=\mathrm{\Omega }\times [-\frac{h}{2}\le z\le \frac{h}{2}]$ , where $h$ is the uniform plate thickness and $\mathrm{\Omega }=[0,a]\times [0,b]$ is the rectangular reference surface lying on the $(x_1,x_2)-$ plane. The plate may be subjected to external traction loads ${\overline{\sigma }}_{i3}{\nu }_3={\overline{t}}_i$ acting at the plate top and/or bottom surfaces ${\mathrm{\Omega }}_{\pm }$ of outward normal ${\nu }_3=\pm z$ . Unless otherwise stated, minuscule Latin indices ( $i,j,k,l,m$ ) range in $\{1,2,3\}$ while Greek indices ( $\alpha ,\beta$ ) take values in $\{1,2\}$ , and a bar/overline $\overline{(\mathrm{\square })}$ denotes prescribed quantities. Adopting index notation and the repeated index summation convention (Einstein’s summation convention), the weak form of the equilibrium equations of the considered problem is expressed by the following virtual work principle:

$${\int }_{\mathrm{\Omega }}\left\{{\int }_h\delta {\epsilon }_{\mathrm{ij}}{\sigma }_{\mathrm{ij}}\mathrm{d}z\right\}\mathrm{d}{x}_1\mathrm{d}{x}_2={\int }_{{\mathrm{\Omega }}_{\pm }}\delta u_i{\overline{t}}_i\mathrm{d}{x}_1\mathrm{d}{x}_2,$$

(11)

where $\delta {\epsilon }_{\mathrm{ij}}$ represents the kinematically compatible virtual strain field ( $\delta$ is the variational operator); it satisfies the linear geometric strain-displacement relations within the small perturbations assumption (see equation (3)), as well as the kinematic boundary conditions defined at the lateral boundary ${\Gamma }_u=\partial \mathrm{\Omega }\times h$ . The stress tensor ${\sigma }_{\mathrm{ij}}$ is symmetric and defined by the material constitutive law, which may include an active deformation expressed by the adopted Stimulus-Expansion-Model (see equation (4)) and the prescribed stimulus difference $\overline{\mathrm{\Delta }S}$ :

$${\sigma }_{\mathrm{ij}}=C_{\mathrm{ijkl}}({\epsilon }_{\mathrm{kl}}-{\alpha }_{\mathrm{kl}}\overline{\mathrm{\Delta }S})=C_{\mathrm{ijkl}}{\epsilon }_{\mathrm{kl}}-{\lambda }_{\mathrm{ij}}\overline{\mathrm{\Delta }S}$$

(12)

Referring to the conventional Voigt notation instead of the tensorial notation, capital indices $P,Q\in \{1,6\}$ are used to identify the six independent strain and stress components. The principle of virtual work can be thus re-written as

$$\begin{array}{cc}\multicolumn{1}{c}{{\int }_{\mathrm{\Omega }}\left\{{\int }_h\delta {\epsilon }_P{C}_{\mathrm{PQ}}{\epsilon }_Q\mathrm{d}z\right\}\mathrm{d}{x}_1\mathrm{d}{x}_2}& ={\int }_{\mathrm{\Omega }}\left\{{\int }_h\delta {\epsilon }_P{\lambda }_P\overline{\mathrm{\Delta }S}\mathrm{d}z\right\}\mathrm{d}{x}_1\mathrm{d}{x}_2\hfill \\ +{\int }_{{\mathrm{\Omega }}_{\pm }}\delta u_i{\overline{t}}_i\mathrm{d}{x}_1\mathrm{d}{x}_2\hfill \end{array}$$

(13)

Let the plate be composed of $p=1,2,\dots N_p$ homogeneous physical plies, each of thickness $h_p$ , stacked along the $x_3\equiv z$ direction. In SGUF, the stack is represented as an assembly of $k=1,2,\dots N_k$ numerical layers or sublaminates, each having a constant thickness $h_k$ composed of $N_p^k$ adjacent physical plies. In each ply of the sublaminate, through-the-thickness assumptions for the displacements $u_i\left(z\right)$ are introduced a priori according to the GUF-type compact notation (Demasi, 2008):

$$u_i^{p,k}(x_1,x_2,z)=\sum _{{\alpha }_{u_i}=0}^{N_{u_i}^k}{F}_{{\alpha }_{u_i}}\left(z\right){\hat{u}}_{{\alpha }_{u_i}i}^{p,k}(x_1,x_2)$$

(14)

For the p^th ply within the k^th sublaminate, the distribution across the thickness of each component of the displacement field is thus expressed as a sum of $N_{u_i}^k$ “thickness functions” $F_{{\alpha }_{u_i}}\left(z\right)$ . The thickness functions are actually expressed in terms of a non-dimensional coordinate ${\zeta }_{\square }\in [-1,+1]$ , which is related to the physical coordinate $z$ depending on whether an Equivalent Single Layer (ESL) or a Layer-Wise (LW) description is adopted for describing the kinematics of the k^th sublaminate:

• In the case of an ESL description, the approximation equation (14) is defined over the whole sublaminate’s thickness $h_k=z_t^k-z_b^k$ , that is, it encompasses all the plies and the non-dimensional coordinate is

$${\zeta }_k=\frac{2}{h_k}z-\frac{z_t^k+z_b^k}{z_t^k-z_b^k}$$

(15)

• In the case of an LW description, the approximation equation (14) is introduced over the thickness $h_p=z_t^p-z_b^p$ of each individual ply of the k^th sublaminate and the non-dimensional coordinate is

$${\zeta }_p=\frac{2}{h_p}z-\frac{z_t^p+z_b^p}{z_t^p-z_b^p}$$

(16)

In either case, the “thickness functions” $F_{\alpha }\left(\zeta \right)$ are built from orthogonal Legendre polynomials that allow to use the classical assembly procedure of, for example, FEM, for constructing the resultant assumption across the whole composite stack, upon enforcing the continuity of the displacement field at adjacent plies and layers. Explicit expressions have been already reported elsewhere (D’Ottavio, 2016; D’Ottavio et al., 2016) and are omitted for the sake of brevity.

The resulting structural model is then defined by indicating the kinematic assumption introduced for each sublaminate the composite stack has been subdivided into. Therefore, each displacement variable within a sublaminate can be either described in an ESL or LW manner, and its polynomial approximation is given by the integer $N_{u_i}^k$ . The stimulus is described always LW and is constant across each active material ply:

$${\overline{\mathrm{\Delta }S}}^{p,k}(x_1,x_2,z)={\overline{\mathrm{\Delta }\hat{S}}}^{p,k}(x_1,x_2)$$

(17)

Note that, in the current framework, gradients of the stimulus may be introduced by further subdividing each active ply into a number of fictitious numerical sublayers, so to provide a step-wise approximation of the stimulus’ distribution across the thickness.

The two-dimensional (2D) model of the smart composite plate is constructed upon introducing the constrained kinematics expressed by equation (14) into the weak formulation equation (11) and carrying out the derivatives and integrals along the thickness direction $z$ . Hence, a generic expression for the virtual internal work contribution of the p^th ply in the k^th sublaminate, defined by the work-conjugate displacement components $\delta u_i$ and $u_j$ reads

$$\begin{array}{c}\multicolumn{1}{c}{\delta W_{\mathrm{int}}^{p,k}={\int }_{\mathrm{\Omega }}\delta \left({\partial }_{x_{\alpha }}\right){\hat{u}}_{{\alpha }_{u_i}i}(x_1,x_2)\left[C_{\mathrm{PQ}}^{p,k}{Z}_{\left({\partial }_z\right)u_i\left({\partial }_z\right)u_j}^{p{\alpha }_{u_i}{\beta }_{u_j}}\right]}\\ \left({\partial }_{x_{\alpha }}\right){\hat{u}}_{{\beta }_{u_j}j}(x_1,x_2)\mathrm{d}{x}_1\mathrm{d}{x}_2\hfill \end{array}$$

(18)

in which the notations $\left({\partial }_{x_{\alpha }}\right)$ and $\left({\partial }_z\right)$ indicate that the displacement components $u_i$ and $u_j$ may require differentiation (along the in-plane directions $x_{\alpha }$ and/or the thickness direction $z$ ) for defining the strain components ${\epsilon }_P$ and ${\epsilon }_Q$ . The thickness integrals have been introduced according to the notation

$$Z_{\left({\partial }_z\right)u_i\left({\partial }_z\right)u_j}^{p{\alpha }_{u_i}{\beta }_{u_j}}={\int }_{h_p}\left({\partial }_z\right)F_{{\alpha }_{u_i}}\left({\partial }_z\right)F_{{\beta }_{u_j}}\mathrm{d}z$$

(19)

The virtual internal work for the whole plate is then obtained by cycling over all indices ${\alpha }_{u_i}$ and ${\beta }_{u_j}$ pertaining to the p^th ply of the k^th sublaminate, and upon assembling the contributions for the whole stack. The virtual external work due to the prescribed traction loads acting at ${\mathrm{\Omega }}_{\pm }$ does not require the computation of thickness integrals and is obtained by only one term by virtue of the orthogonality of the Legendre polynomials. Details for these contributions and for the assembly procedures are provided as indicated in (D’Ottavio, 2016) and are here omitted for the sake of brevity.

The virtual external work given by actuation through the prescribed stimulus ${\overline{\mathrm{\Delta }\hat{S}}}^{p,k}$ of the smart material ply $p$ in the k^th sublaminate is expressed in analogous manner as

$$\delta W_{\mathrm{act}}^{p,k}={\int }_{\mathrm{\Omega }}\delta \left({\partial }_{x_{\alpha }}\right){\hat{u}}_{{\alpha }_{u_i}i}(x_1,x_2)\left[{\lambda }_P^{p,k}{Z}_{\left({\partial }_z\right)u_i}^{p{\alpha }_{u_i}}\right]{\overline{\mathrm{\Delta }\hat{S}}}^{p,k}(x_1,x_2)\mathrm{d}{x}_1\mathrm{d}{x}_2$$

(20)

with

$$Z_{\left({\partial }_z\right)u_i}^{p{\alpha }_{u_i}}={\int }_{h_p}\left({\partial }_z\right)F_{{\alpha }_{u_i}}\mathrm{dz}$$

(21)

The virtual external work of the actuation for the whole stack is obtained by cycling over all indices ${\alpha }_{u_i}$ of all active plies, and upon assembling the contributions for the whole stack. Note that the stimulus needs not to be continuous across adjacent active plies or layers.

It is worth noticing that models based on a thickness-wise constant transverse displacement $u_3$ , that is, models with $N_{u_3}^k=0$ , are known to suffer Poisson locking if the 3D constitutive law is retained (Carrera and Brischetto, 2008). Therefore, for those models the modified “reduced” constitutive coefficients are used in order to introduce the plane stress condition ${\sigma }_{33}=0$ :

$$C_{\mathrm{PQ}}^{*}=C_{\mathrm{PQ}}-\frac{C_{P3}{C}_{3Q}}{C_{33}};{\lambda }_P^{*}={\lambda }_P-\frac{{\lambda }_3{C}_{P3}}{C_{33}}.$$

(22)

2.3. Navier-type solution

Weak-form solutions of the resulting two-dimensional Boundary Value Problem (BVP) defined over the reference surface $\mathrm{\Omega }$ can be obtained by introducing some shape functions

$${\hat{u}}_{{\alpha }_{i}i}(x_1,x_2)=N_r(x_1,x_2)U_{r{\alpha }_{i}i}$$

(23)

and carrying out the integrals over $\mathrm{\Omega }$ of equation (18). Ritz and FE methods have been presented in D’Ottavio et al. (2016) and Di Cara et al. (2023), respectively. In this work, however, we adopt the strong-form solution of Navier type, which solves in strong form of the BVP of the plate, as presented in D’Ottavio (2016). The differential equations over the 2D domain are obtained from the virtual work principle built from equation (13) by integrating by parts those terms whose virtual variations appear differentiated with respect to the in-plane axes $x_1$ or $x_2$ . The Navier-type solution adopts the following trigonometric shape functions (Reddy, 2003):

$${\hat{u}}_{{\alpha }_{1}1}(x_1,x_2)=U_{{\alpha }_1}\cos \frac{m\pi x_1}{a}\sin \frac{n\pi x_2}{b}$$

(24a)

$${\hat{u}}_{{\alpha }_{2}2}(x_1,x_2)=V_{{\alpha }_2}\sin \frac{m\pi x_1}{a}\cos \frac{n\pi x_2}{b}$$

(24b)

$${\hat{u}}_{{\alpha }_{3}3}(x_1,x_2)=W_{{\alpha }_3}\sin \frac{m\pi x_1}{a}\sin \frac{n\pi x_2}{b}$$

(24c)

with $m,n$ integers defining the number of half-waves. The prescribed tractions and the stimulus are required to have analogous distributions:

$${\overline{t}}_1(x_1,x_2,z_{\pm })={\overline{T}}_{1\pm }\cos \frac{m\pi x_1}{a}\sin \frac{n\pi x_2}{b}$$

(25a)

$${\overline{t}}_2(x_1,x_2,z_{\pm })={\overline{T}}_{2\pm }\sin \frac{m\pi x_1}{a}\cos \frac{n\pi x_2}{b}$$

(25b)

$$\{{\overline{t}}_3\left(z_{\pm }\right),\overline{\mathrm{\Delta }\hat{S}}\}(x_1,x_2)=\{{\overline{T}}_{3\pm },\overline{\mathrm{\Delta }S}\}\sin \frac{m\pi x_1}{a}\sin \frac{n\pi x_2}{b}$$

(25c)

This solution fulfils in strong form the 2D differential equations over the domain $\mathrm{\Omega }$ and the boundary conditions at $\partial \mathrm{\Omega }$ only under some stringent limitations: the material plies are required to have orthotropy axes aligned with the plate axes $x_1,x_2$ , and the kinematic boundary conditions must correspond to the “simply-support” conditions:

$$u_1\left(x_1,x_2,z\right)=0\mathrm{at}{x}_1=0,a$$

(26)

$$u_2(x_1,x_2,z)=0\mathrm{at}{x}_2=0,b$$

(27)

$$u_3\left(x_1,x_2,z\right)=0\mathrm{at}{x}_1=0,a\mathrm{and}{x}_2=0,b$$

(28)

This way, the proposed model introduces approximations in weak form only for the through-the-thickness distributions. If the external load is distributed over $\mathrm{\Omega }$ as a single harmonic function, foe example, as given by equation (25c) with $m=n=1$ , the Navier solution furnishes a single-term closed-form solution. For other distributions of the external loading, the Fourier series development is invoked for representing the resulting response as a sum of a finite number of harmonics. A single harmonic, as well as a Fourier series with a varying number of terms, is used for an active-passive hydrogel composite in Section 3.1. Afterwards, the convergence of the solution towards the Finite Element solution (with a constant stimulus) for an increasing number of harmonics is shown.

The discrete system The virtual work principle yields the following linear algebraic system for each harmonic $m,n$ :

$$\left[\begin{array}{ccc}{\mathit{K}}_{\mathrm{UU}}\hfill & {\mathit{K}}_{\mathrm{UV}}\hfill & {\mathit{K}}_{\mathrm{UW}}\hfill \\ \multicolumn{1}{c}{{\mathit{K}}_{\mathrm{UV}}^{\top }}& {\mathit{K}}_{\mathrm{VV}}\hfill & {\mathit{K}}_{\mathrm{VW}}\hfill \\ \multicolumn{1}{c}{{\mathit{K}}_{\mathrm{UW}}^{\top }}& {\mathit{K}}_{\mathrm{VW}}^{\top }\hfill & {\mathit{K}}_{\mathrm{WW}}\hfill \end{array}\right]\left[\begin{array}{c}\mathit{U}\hfill \\ \multicolumn{1}{c}{\mathit{V}}\\ \multicolumn{1}{c}{\mathit{W}}\end{array}\right]=\left[\begin{array}{c}{\mathit{F}}_{\mathrm{US}}\hfill \\ \multicolumn{1}{c}{{\mathit{F}}_{\mathrm{VS}}}\\ \multicolumn{1}{c}{{\mathit{F}}_{\mathrm{WS}}}\end{array}\right]+\left[\begin{array}{c}{\mathit{T}}_U\hfill \\ \multicolumn{1}{c}{{\mathit{T}}_V}\\ \multicolumn{1}{c}{{\mathit{T}}_W}\end{array}\right]$$

(29)

whose unknowns are the vectors of the generalized displacement variables that define the through-thickness approximations for the whole composite stack. The external work is defined by the actuation stimuli present in the whole stack, which is represented by the vectors ${\mathit{F}}_{\mathrm{US}}$ , ${\mathit{F}}_{\mathrm{VS}}$ , and ${\mathit{F}}_{\mathrm{WS}}$ , and by the external tractions represented by the vectors $\mathit{T}$ . The size of the system only depends on the approximation orders used to define the model. For instance, for a plate composed of $N_p$ plies, an ESL model adopting the classical FSDT kinematics will yield a system with only 5 unknown variables, whereas if FSDT is used for each ply independently ( $L{D}_{1,0}$ model) there will be $2(N_p+1)+1=2N_p+3$ variables.

Please note that the Navier solution requires the expansion coefficients to be constant over the reference surface $\mathrm{\Omega }$ , see equation (20). However, since the coefficients for the active behavior are stimulus-dependent, they should vary locally due to the sine distribution as defined in equation (25c). This inconsistency is also accounted for in the Finite Element simulations that will be used as reference solution, in which the coefficients are all considered as constant inside each physical ply. For the sake of completeness and further reference, we mention that an alternative strong form solution method that allows accounting for in-plane variation of coefficients is the Differential Quadrature Method (DQM) (Groh and Weaver, 2014; Tornabene and Fantuzzi, 2021).

3. Results: Simulation results for selected active-passive composites

The physics of the two different active material classes (piezo-ceramics and hydrogels) were implemented into the sublaminate Generalized Unified Formulation presented in Section 2.2. The three testcases illustrated in Figure 5 are investigated. While a multitude of other realizations is of course possible, the actual testcases are chosen for the following reasons:

• Testcase (a) in Figure 5(a) (Section 3.1) is chosen due to the simple setup and the direct link to our previous works on hydrogel composites (Ehrenhofer, 2023a) (including the availability of experimental material data in (Ehrenhofer et al., 2016)).

• Testcase (b) in Figure 5(b) (Section 3.2) concerns a standard benchmark for piezoelectric actuators (Fernandes and Pouget, 2002, 2003) and allows to investigate the effect of the SEM with respect to the full piezoelectric coupling.

• Testcase (c) in Figure 5(c) (Section 3.3) comprises a more complex passive structure (an orthotropic heterogeneous laminate), for which reference results are available in literature (Cheng and Batra, 2000), and it allows to showcase the capability of the sGUF approach to adapt the modeling and computational effort depending on the output of interest.

The various models that can be obtained by the sGUF approach are named according to the following system of acronyms. For each sublaminate, the expansion orders of the variables are $N_{u_1}=N_{u_2}:=N_{u_{\alpha }}$ for the in-plane displacements and $N_{u_3}$ for the out-of-plane deflection. A sublaminate with an LW description is named $L{D}_{N_{u_{\alpha }}{N}_{u_3}}$ . If $N_{u_{\alpha }}=N_{u_3}:=N$ the CUF acronym $\mathrm{LDN}$ is used. If an ESL description is used for the sublaminate, the capital “L” is replaced by “E”. FSDT is the First order Shear Deformation Theory and is therefore denoted in GUF notation as $E{D}_{10}$ (ESL with through-thickness constant $u_3$ and through-thickness linear $u_1$ and $u_2$ ). If several sublaminates are used for a composite stack, the acronyms for each sublaminate are given separated by a slash.

Figure 5.

Chosen testcases (not on scale) for active-passive composites: (a) testcase 1: Hydrogel-PET-composite, (b) testcase 2: PZT bimorph, (c) testcase 3: PZT bimorph with laminated basis structure.

3.1. Validation with a hydrogel composite

As shown in our previous work (Ehrenhofer et al., 2023), a two-layer composite with an active top layer made of the hydrogel pNiPAAm and a passive bottom layer can be used for validation of the approach. This kind of setup was also considered to derive interaction cases in Soft-Hard Active-Passive Embedded Structures (Ehrenhofer, 2023a). The hydrogel is actuated with a very small distributed temperature change of ΔS = −0.1 K (cooling leads to positive expansion) to preserve small deformations. Along the axial direction, two stimulus distributions will be considered, namely a sinusoidal variation as defined by equation (25c) and a uniform distribution, which is obtained by summing up the individual sinusoidal harmonics according to Fourier series decomposition. The reference temperature is ${\vartheta }_{\mathrm{ref}}=32°\mathrm{C}$ . At $31.9°\mathrm{C}$ , the elastic modulus is E^HG = 0.0014 MPa and the isotropic expansion coefficient α^HG = −7.5316°C⁻¹. This leads to an isotropic active strain of ${\epsilon }^{\mathrm{active}}=0.0235$ for the top layer. For the passive soft bottom layer, the elastic modulus is constant with E^passive = 2.1 MPa.

The active hydrogel layer has a thickness of h₂ = 2 mm, the passive layer of h₁ = 1 mm. The plate dimensions are ℓ₁ = ℓ₂ = 100 mm. A plane strain model in the $(x_1,x_3)-$ plane is considered without loss of generality by setting $n=0$ in equations (25) and (25). The simply-supported boundary conditions of Navier’s solution are reproduced in the FE model by setting

$$\begin{array}{c}{u}_z(x=0,z)=0;u_z(x={\ell }_1,z)=0;\hfill \\ u_x(x=\frac{{\ell }_1}{2},z)=0;z\in [0,h_1+h_2]\hfill \end{array}$$

(30)

A total of 3000 biquadratic quadrilateral elements of type CPE8H (hybrid with linear pressure to account for quasi-incompressibility, quadratic for avoiding locking issues) were used in a linear perturbation step in Abaqus. A refined sGUF model has been used which consists of one sublaminate per physical ply; the highly refined, quasi-3D model $\mathrm{LD}6$ has been used for the active hydrogel ply, whereas the simple $E{D}_{10}$ (FSDT) is used for the passive layer.

Figure 6 displays the displacement and stress distributions obtained by the Abaqus FEM with Stimulus-Expansion-Model and by our new sGUF-SEM model. The case of a sine-distributed stimulus (thus of a single-term closed-form Navier solution) is displayed in Figure 6(a). An excellent agreement for the axial displacement $u_x$ and the stress ${\sigma }_{\mathrm{xx}}$ is obtained. A small discrepancy of less than $1\%$ is found for the out-of-plane deflection $u_z$ , which can be attributed the handling of quasi-incompressibility with hybrid elements. These results show that our sGUF-SEM model can adequately predict these distributions for the considered active-passive-layer case.

Figure 6.

Distributions of $u_x(x=0,z)$ , $u_z(x={\ell }_1/2,z)$ , and ${\sigma }_{\mathrm{xx}}(x={\ell }_1/2,z)$ obtained from our sGUF model (solid line –) and from Abaqus FEM with SEM (dots °). Results for (a) sine-distributed temperature load and (b) constant temperature load. The unscaled deformed geometry at the top is shown to highlight the positions for the results $u_x,u_z$ and ${\sigma }_{\mathrm{xx}}$ .

Figure 6(b) reports the results for the case of a stimulus that is constant along the $x_1-$ direction. In this case, the sGUF-SEM results depend on the number of harmonics that have been considered to build the Fourier-series solution. For a number of $M=\{51,101,201\}$ Fourier coefficients, a convergence towards the Finite Element solution is obtained. Furthermore, a constant stimulus distribution leads to a higher bending deflection compared to the sine-distributed stimulus, as expected. Please note that the local stress distribution in the central region $x=0$ is very similar between the two cases, because it is the region in which the maximum stimulus occurs in both cases.

3.2 PZT bimorph actuator

The configuration of a PZT parallel bimorph actuator is studied as depicted in Figure 5(b) and investigated by Fernandes and Pouget (2002, 2003). It consists of two equal PZT-4 plies of length ℓ₁ = 25 mm and thickness $h_{\mathrm{PZT}}=h/2$ , whose polarization axis is aligned with $x_3$ . The material properties of the PZT-4 piezo-ceramics are listed in Table 1. The PZT material is transversely isotropic in the $(x_1,x_2)-$ plane perpendicular to the polarization axis. For the sake of completeness, we recall the relation allowing to obtain the piezoelectric strain coefficients $\tilde{\mathit{d}}={\tilde{d}}_{\mathrm{iP}}$ from the piezoelectric stress coefficients $\mathit{e}=e_{\mathrm{iP}}$ listed in Table 1: $\tilde{\mathit{d}}=\mathit{e}{S}^E={\mathit{e}{C}^E}^{-1}$ . The shape of the bi-layer is defined by the length-to-thickness ratio ${\ell }_1/h$ . The structure is investigated in cylindrical bending by referring to the plane strain assumption in the $(x_1,x_3)-$ plane. Both PZT layers are actuated with a voltage difference $\mathrm{\Delta }\tilde{\Psi }=\pm 50\mathrm{V}$ (sign according to the direction in Figure 4(b)) acting uniformly across the thickness. This is accomplished by an electrical circuit that imposes the voltage difference between the electrode located at the bi-layer interface and the two electrodes located at the outer (top and bottom) surfaces. The electrodes being equipotential surfaces, the electrical stimulus is uniform along the Navier solution and for each harmonic $m$ , the stimulus is given as

Table 1.

Material properties of PZT-4 and graphite/epoxy.

Property	PZT-4	Gr/Ep
$E_1$ [GPa]	81.239	132.39
$E_2$ [GPa]	81.239	10.756
$E_3$ [GPa]	64.073	10.756
${\nu }_{12}$ [-]	0.3274	0.24
${\nu }_{13}$ [-]	0.4345	0.24
${\nu }_{23}$ [-]	0.4345	0.49
$G_{44}$ [GPa]	25.60	3.606
$G_{55}$ [GPa]	25.60	5.654
$G_{66}$ [GPa]	30.60	5.654
$e_{31}$ [ $C/m^2$ ]	-5.2	0
$e_{33}$ [ $C/m^2$ ]	15.1	0
$e_{15}$ [ $C/m^2$ ]	12.7	0

$$\mathrm{\Delta }{S}_m=\frac{\mathrm{\Delta }{\tilde{\Psi }}_m}{h_{\mathrm{PZT}}}=\frac{1}{h_{\mathrm{PZT}}}\frac{4{\tilde{\Psi }}_0}{m\pi }$$

(31)

where ${\tilde{\Psi }}_0$ =50 V is the electric potential prescribed at the outer electrodes (the central electroded is considered grounded, ( $\tilde{\Psi }$ =0 V)) (Fernandes and Pouget, 2002, 2003).

The results for ${\ell }_1/h=10$ are illustrated in Figure 7, where present results obtained with $M=101$ terms of the Fourier series development and the LD6 model are compared against a 3D FEM solution obtained with the biquadratic PLANE 223 elements of ANSYS with full piezoelectric coupling (Polit et al., 2016).

Figure 7.

Validation of the present modeling approach for a moderately thick ( ${\ell }_1/h=10$ ) PZT bimorph structure: (a) displacement $u_x(x=0,z)$ , (b) displacement $u_z(x=\frac{{\ell }_1}{2},z)$ , (c) stress ${\sigma }_{\mathrm{xx}}(x=\frac{{\ell }_1}{2},z)$ .

Similarly to the previous hydrogel-based actuator, an excellent agreement between the FEM solution and the present sGUF-SEM results is found for the in-plane displacement $u_x$ and the bending stress ${\sigma }_{\mathrm{xx}}$ . However, a certain discrepancy can again be found with regard to the out-of-plane deflection $u_z$ . In particular, present results are shown to overestimate the deflection by 4.5%. This is due to the fundamental assumption of the SEM approach, which disregards the back-coupling effect from the mechanical deformation to the electric field. In fact, the electrical stimulus, assumed to be constant in SEM, is overestimated because the piezoelectric back-coupling actually reduces the effective electric field. To illustrate this, two additional curves are plotted in Figure 6(b). They are obtained by two LD4 models retaining the full piezoelectric coupling (i.e. direct and converse effect) as presented by Ballhause et al. (2005). The first model adopts a quadratic electric potential $\tilde{\Psi }\left(z\right)$ across the thickness: the quadratic approximation allows to represent the electric field induced by the mechanical bending. The second model relies on a linear distribution of $\tilde{\Psi }\left(z\right)$ , which thus neglects the bending-induced field. Comparing these curves it is apparent that the induced field reduces the effective driving field, thus yielding a lower bending deflection. Therefore, since the present SEM approach disregards the back-coupling effects, it is expected to provide upper bounds estimates for the actuation strain.

Influence of plate thickness The structural modeling of the presented sGUF-SEM approach is applied here for investigating PZT bimorph configurations of different length-to-thickness ratios. The previously studied moderately thick bimorph ( ${\ell }_1/h=10$ ) is compared against bimorphs characterized by ${\ell }_1/h=50$ (thin) and ${\ell }_1/h=100$ (very thin). The driving voltage difference is kept fixed to $\mathrm{\Delta }\tilde{\Psi }={\tilde{\Psi }}^{\mathrm{outer}}-{\tilde{\Psi }}^{\mathrm{mid}}=50\mathrm{V}$ . Due to the equations of the electric field in a capacitance, this corresponds for the very thin actuator to an electric field intensity of ${\tilde{E}}_3=400\mathrm{V}/\mathrm{mm}$ . The results for $u_x,u_z$ and ${\sigma }_{\mathrm{xx}}$ are reported in Figure 8. As expected, the active bending becomes larger as the bimorph structure becomes thinner. This is on the one hand due to the higher field intensity obtained at iso-voltage, on the other hand due to the reduced bending stiffness. As a closure to this testcase, it is worth emphasizing that in the present modeling approach, it is sufficient to input the updated geometric dimensions and to re-run the simulation for obtaining the updated results, that is, there are no issues with slenderness ratios. This is not the case for commercial FE packages, in which only solid elements are equipped with piezoelectric coupling capabilities and which are known to require adequate meshing strategies for avoiding excessive element distortions.

Figure 8.

PZT bimorph actuated by $\mathrm{\Delta }\tilde{\Psi }=50\mathrm{V}$ : influence of the length-to-thickness ratio on the actuation: (a) displacement $u_x(x=0,z)$ , (b) displacement $u_z(x=\frac{{\ell }_1}{2},z)$ , (c) stress ${\sigma }_{\mathrm{xx}}(x=\frac{{\ell }_1}{2},z)$ .

3.3 PZT bimorph with laminated base structure

After having validated the proposed approach for the case of piezoelectric actuation, a more complex active-passive structure is considered in order to showcase the model reduction capability of sGUF-SEM. To this aim, an extension of the previously considered PZT-4 bimorph is considered, in which the parallel polarized piezo-ceramic layers are bonded to a composite basis structure consisting of two Graphite-Epoxy (Gr/Ep) plies. The resulting smart laminate is thus a $[p/90/0/p]$ composite stack, with $p$ denoting the piezo-layers. As for the bimorph actuator, each PZT-4 layer is actuated by an imposed electric field, which is driven by the voltage difference between the outer electrodes set at ${\tilde{\Psi }}^{\mathrm{out}}=V_0$ and the internal electrodes separating the passive laminate and the active layers that are grounded ( $\tilde{\Psi }$ ^int = 0 V). As a result, a bending movement of the smart laminate is attained. This problem has been solved in strong-form by Cheng and Batra (2000) in the framework of an asymptotic approach. The smart laminate is a square with a side length of ${\ell }_1={\ell }_2$ and a total thickness of $h$ . The thickness of each ply of the basis laminate and of the piezoelectric plies is $0.4h$ and $0.1h$ , respectively. The thickness of the electrodes as well of the bonding adhesive is neglected. Material properties in the respective principal directions are listed in 1; both materials are considered as transversely isotropic: the plane of isotropy is perpendicular to the polarization direction $\left(x_3\right)$ for PZT and to the fiber direction $\left(x_1\right)$ for Gr/Ep. For allowing a direct comparison with the reference results of Cheng and Batra (2000), only the first harmonic with $m=n=1$ of the Fourier series development is considered (i.e. the stimulus is taken as distributed bi-sinusoidally over the the surfaces), thus a one-term closed-form Navier solution is obtained. The results are presented in non-dimensional form as

$${\overline{u}}_i=\frac{c^{*}}{V_0{e}^{*}}{u}_i;{\overline{\sigma }}_{\mathrm{ij}}=\frac{L}{V_0{e}^{*}}{\sigma }_{\mathrm{ij}}$$

(32)

with $c^{*}=1\mathrm{N}/{\mathrm{m}}^2$ and $e^{*}=1\mathrm{C}/{\mathrm{m}}^2$ .

The out-of-plane deflection obtained for different length-to-thickness ratios ${\ell }_1/h$ and with various SGUF models are reported in Section 2 and compared to the reference solution by Cheng and Batra (2000).

The most refined model employed in this case adopts 6^th order polynomials in each layer (LD6). A very close agreement with the reference results is obvious. The other two sGUF models considered in Table 2 adopt FSDT for the piezo-layers while the basis cross-ply laminate is modeled with the quasi-3D model LD6, the high-order ESL model $E{D}_{32}$ and FSDT. It should be noted that FSDT calls for the use of reduced stiffnesses and reduced piezoelectric coupling coefficients in order to account for the plane stress conditions ${\sigma }_{33}=0$ . Finally, the most simple model is considered in which the whole smart laminate is modeled according to FSDT. No relevant discrepancy is appreciable between all these models: even the most simple model (FSDT) grasps the actuator’s global deflection with good accuracy.

Table 2.

Piezoelectric composite bimorph: non-dimensional out-of-plane deflection ${\overline{u}}_3(\frac{{\ell }_1}{2},\frac{{\ell }_2}{2},0)\times {10}^8$ (the factor ${10}^8$ was chosen for easier presentation).

${\ell }_1/h$	Ref (Cheng and Batra, 2000)	LD6	FSDT/LD6	FSDT/ED32	FSDT/FSDT	FSDT
4	0.2349	0.2379	0.2375	0.2374	0.2332	0.2261
10	1.4166	1.4257	1.4253	1.4251	1.4204	1.4132
20	5.6341	5.6655	5.6652	5.6647	5.6601	5.6530
50	35.1561	35.3435	35.3432	35.3410	35.3381	35.3310

The distributions across the thickness of the in-plane stresses obtained by the most refined model (LD6) and the most simple model (FSDT) are reported in Figure 9 for ${\ell }_1/h=4$ and ${\ell }_1/h=50$ . The result confirms that FSDT can accurately predict the bending motion, thus with accurate bending strains and associated stresses.

Figure 9.

Non-dimensional in-plane stress (a, b) and transverse stress (c, d) distributions across the thickness of the piezoelectric composite bimorph: (a) in-plane normal displacement ${\overline{u}}_x$ , (b) out-of-plane displacement ${\overline{u}}_z$ , (c) in-plane normal stress ${\overline{\sigma }}_{\mathrm{xx}}$ , (d) in-plane shear stress ${\overline{\sigma }}_{\mathrm{xy}}$ , (e) transverse shear stress ${\overline{\sigma }}_{\mathrm{xz}}$ , (f) transverse normal stress ${\overline{\sigma }}_{\mathrm{zz}}$ .

The situation is more complicated if the transverse stress field is considered. Figure 9(e) and (f) allow to assess the accuracy of the transverse stress distributions across the thickness obtained by several models. It is obvious that only the refined model LD6 accurately grasps the steep gradients characterizing the response across the composite stack. All equilibrium conditions are seen to be correctly met, with a perfect match of the strong form solution by Cheng and Batra (2000). A high transverse shear stress peak is seen to occur at the top interface. This appears to be the most critical interface where delamination is likely to occur. Simpler models do not allow to perform such an accurate stress analysis, unless dedicated post-processing techniques are adopted (Rohwer et al., 2005).

The presented cases show (i) the adaptability of the model to active materials with very different properties, and (ii) the influences of the model’s kinematics order. A smart structure consisting of a combination of these two chosen materials—hydrogels and piezo-ceramics—could be simulated by the present prototypical implementation. However, it should be mentioned that such a setup would require unrealistically thin piezo-layer (Case II SHAPES of Ehrenhofer (2023a): With the ratio in elastic modulus of $E_{\mathrm{PZT}}/E_{\mathrm{Hydrogel}}\approx 81\mathrm{GPa}/1.4\mathrm{kPa}\approx 6·{10}^7$ , a practically infeasible high ratio in layer thickness would in fact be required. Therefore, combining piezo-electric layers with hydrogels would rather be meaningful by adopting piezo-polymers instead of piezo-ceramics: the model for piezo-polymers will be derived in future works.

4. Conclusions

In the current paper, we presented a combination of approaches to model active layered composite structures. The aim was to provide an efficient simulation tool dedicated to the design of composite actuators, in which active material layers are combined with passive substructures. The active layers are modeled according to the Stimulus Expansion Model (SEM): based on the analogy principle, the active deformation is represented as an equivalent thermal expansion. Equivalent thermal expansion coefficients are exemplarily derived for two classes of active materials, namely hydrogels and piezo-ceramics, and included within a variable kinematics structural modeling approach referred to as sublaminate Generalized Unified Formulation (sGUF). This provides a unique software framework that allows to formulate reduced order models of hierarchic complexity, which can efficiently resolve global and local gradients depending on the designer’s intended purpose.

The combined sGUF-SEM model has been protypically implemented along with the analytical Navier solution to demonstrate the efficiency of the approach. In conjunction with the Fourier series development, arbitrary stimulus distributions along the in-plane directions can be taken into consideration. Three example applications have been presented for showcasing the main features of the sGUF-SEM approach:

• Numerical validation against reference solutions obtained from 3D FEM or the open literature has been successful; the proposed reduced order structural modeling allows to reduce the computational effort required by full three-dimensional FEM models.

• The variable kinematics approach allows the designer to adapt the modeling and computational effort depending on the requested output quantity: low-fidelity models can be used to identify gross response parameters (e.g. deflection, bending stress) while refined high-order models arrive at accurately grasping local stress gradients required for a strength analysis.

• The approach is easily extensible for the inclusion of additional active materials, as well as additional and distributed field influences.

As a result, the sGUF-SEM approach can be effectively applied for high-throughput screening of different active-passive and active-active material pairings, eventually leading to the identification of novel smart composites. The first promising results presented in this paper call for further developments of the approach. Future works shall be directed towards the integration of other classes of active layers (e.g. Dielectric Elastomer Actuators, Ionic Polymer Metal Composites, Conductive Polymers, Shape Memory Alloys). Furthermore, a more refined representation of stimuli distributions inside the active layers shall be considered. In particular, while still refraining from resolving the two-way multifield coupling, the stimulus-dependent mechanical parameters and expansion coefficients of the active material shall be taken into consideration. To this aim, the approach developed for Functionally Graded Materials (FGMs) will be carried over to represent continuously varying constitutive properties. Weak-form solutions over the in-plane domain shall be considered as well in order to broaden the range of applicability of the sGUF-SEM approach. The extension towards geometric nonlinearities is eventually required for obtaining meaningful results for practice-oriented actuator systems implying large displacements. With these additions, the proposed combination of sGUF and SEM will be a valuable tool in the design of new actuator concepts and real-world application of Soft-Hard Active Passive Embedded Structures (SHAPES) as functionally integrated lightweight structures.