Large deformation modeling of flexible piezoelectric materials

Abstract

The applications of piezoelectric materials in the field of smart structures have received significant attention from both the communities of science and engineering. Numerous experimental studies have been carried out to endow smart structures with good flexibility. The flexible/stretchable piezoelectric materials are developed to fit this emerging trend. Generally, these materials undergo significant deformation before reaching fracture failure, and they often exhibit a stress-softening phenomenon during the deformation process. However, the traditional linear constitutive model, typically used for rigid piezoelectric ceramics, continues to dominate theoretical and modeling processes in many scenarios. Existing nonlinear constitutive models usually introduce additional coefficients besides elastic, piezoelectric, and dielectric coefficients. Determining these coefficients requires a substantial number of experiments. In this work, based on a Neo-Hookean material model and electromechanical theory, a novel model for flexible piezoelectric material considering large deformation has been established. In contrast with existing models, the present model describes the nonlinear behavior of flexible piezoelectric material without the need for introducing additional parameters. Furthermore, this model exhibits a quadratic dependence of stress on the electric field. To facilitate practical applications, the constitutive model has been implemented using the commercial simulation software ABAQUS through a user subroutine. The accuracy of the subroutine is validated by comparing simulations with analytical solutions for uniaxial stretching of a flexible piezoelectric ribbon. Several numerical examples are followed to demonstrate the robustness of the elements. The proposed model offers a valuable tool for the analysis and design of flexible piezoelectric material.

Introduction

Recently, the developments of flexible electronics have greatly facilitated the applications of soft robots [1], flexible wearable devices [2], as well as soft skins [3, 4–5]. These flexible devices are usually fabricated by integrating multifunctional materials with flexible organic substrates [6]. Among them, flexible devices incorporating piezoelectric materials have attracted widespread attention due to their unique mechanical–electrical energy conversion capabilities. Compared with other flexible materials, piezoelectric materials exhibit superior mechanical robustness [7], noise resistance [8], and rapid mechanical response [9, 10].

Numerous experimental studies have been done to design piezoelectric structures [6, 11, 12]. With the increasing applications of piezoelectric materials in smart sensors [7, 12] and actuators [9], there is a particular focus on predicting and modeling their electromechanical performances. In the previous research, piezoelectric ceramics dominated the field of device applications, and a linear constitutive model considering small deformation was efficient for analyses. However, when it comes to piezoelectric materials used in lightweight space structures [13], actuators [14], and interdigital transducers [13], they are often subjected to large external loads that induce significant deformations. For instance, Yi et al. [15] derived the dynamic responses of piezoelectric structures with consideration of large geometrical deformation based on the updated Lagrangian formulation. Mukherjee and Chaudhuri [16] studied the large deformation of a piezoelectric bimorph beam subjected to transverse and axial loadings. However, these studies still relied on traditional piezoelectric linear constitutive relations.

To meet the growing demands for designing dexterous mechanisms, smart structures usually incorporate flexible piezoelectric materials with high stretchability [17, 18]. These materials exhibit significant deformation before experiencing fracture failure. In the pursuit of enhanced performance, Yoon and Kelarakis [19] utilized organically modified Lucentite nanoclay to optimize electrospun PVDF membranes. The inclusion of nanoclay resulted in the achievement of high-performance membranes. The neat PVDF membranes exhibited an elongation at a break of 87%. However, when incorporating 1 wt% and 10 wt% of Lucentite, the elongation at break increased to approximately 108% and 145%, respectively. Notably, these materials exhibit stress-softening phenomena during large deformation [19, 20–21]. This phenomenon cannot be adequately described by the widely used linear piezoelectric constitutive model. There is a need to develop new constitutive models that can accurately capture the nonlinear behavior exhibited by these materials. In response to this challenge, several attempts have been made to establish nonlinear constitutive models. For example, Joshi [22] extended the linear constitutive formulations by considering higher-order effects between physical quantities. Wagner and Hagedorn [23] proposed nonlinear constitutive equations based on the beam vibration experiments, taking into account the dependence of Young’s modulus and piezoelectric factor on strain. Tajeddini and Muliana [24] developed a second-order electric field constitutive model to investigate the large deformation of a flexible piezoelectric cantilever induced by a large electric field stimulation. However, these studies introduced additional coefficients beyond the elastic, piezoelectric, and dielectric coefficients. Determining these coefficients often presents challenges as it requires a large number of experiments, thus posing difficulties for practical applications.

The finite element method is the main design tool for mechanical analyses and structure design [25]. In order to simulate the use of piezoelectric devices in engineering, several studies have been carried out to develop the user subroutine element based on the piezoelectric elements available in commercial software such as ABAQUS and ANSYS. Nestorović et al. [26] proposed a biquadratic piezoelectric composite shell element, while Js and Ahaa [27] developed the thermo-piezoelectric element. Serrao and Kozinov [28] presented a three-dimensional mixed finite element method specifically designed for high-order electromechanical applications, where flexoelectricity in piezoelectric solids was taken into account. With increasing attention to flexible piezoelectric materials, numerical capabilities need to be enhanced to accurately simulate the large deformations exhibited by these materials [29]. Smart thin-walled structures have significant potential in practical applications, leading to extensive investigations on shell elements capable of handling large deformations [30, 31]. Notably, a persistent challenge in shell element analysis is the handling of rotation degrees of freedom. In contrast, solid elements demonstrate superior behavior in finite element analyses. The solid element proposed by Yi et al. [15], which was developed using the combination of Green strain and linear constitutive relations, may fail to accurately simulate the large deformations observed by flexible piezoelectric devices. Consequently, there is a clear necessity to develop piezoelectric solid elements that can effectively account for and accurately represent these large deformations.

Motivated by large deformations of flexible piezoelectrics, the objective of this work is to establish a theoretical model and a user subroutine for piezoelectric materials capable of accommodating such deformations. As previously reported in the literature [32, 33, 34–35], soft/flexible materials are widely recognized to exhibit anisotropic and viscoelastic behavior. For simplification, we have assumed neglecting anisotropic and viscoelastic effects in the subsequent derivations. To implement the established constitutive equations in simulations, a user element subroutine (UEL) in commercial software ABAQUS is utilized. The remainder of this paper is organized as follows. Section 2 provides an overview of the continuum framework employed to model coupled electromechanical behavior. In Sect. 3, the constitutive equations for a flexible piezoelectric material are derived. To validate the accuracy of the constitutive model, an analytic solution of a uniaxial stretching case is presented in Sect. 4. In addition, the tensile experiment is also performed on polyvinylidene fluoride (PVDF) ribbon. In Sect. 5, the building process of the new element is described in detail. Section 6 presents several numerical examples to demonstrate the capabilities of the developed user element. Finally, the work is concluded in Sect. 7.

Continuum framework

Consider a homogeneous body with its reference state B_r, which corresponds to the state without mechanical and electric loadings, defined in Cartesian coordinates. If X represents an arbitrary material point in B_r, then x represents its current position in the deformed configuration B_t. The deformation gradient F is defined as the derivative of the current position vector x with respect to the reference position vector X, given by $\mathbf{F}=\partial \mathbf{x}\left./\phantom{\partial \mathbf{x}\partial \mathbf{X}}\right)\partial \mathbf{X}$. The right Cauchy–Green deformation tensor is denoted as C = F^TF, whose first principal invariant is I₁ = tr C.

Assuming an electric potential $\varphi (\mathbf{X},t)$ is applied to the material particle. Define the referential electric field intensity as the gradient of electric potential, i.e., ${\tilde{E}}_i=-\frac{\partial \varphi (\mathbf{X},t)}{\partial X_i}$, where the negative sign conforms to the conventional direction for the electric field, going from high potential to low. Likewise, the spatial electric field intensity is defined as $E_i=-\frac{\partial \varphi (\mathbf{x},t)}{\partial x_i}$. Herein X_i and x_i represent the components of X and x in their Cartesian coordinates, respectively.

The equilibrium of forces in the deformed body B_t is expressed as follows:

1

$$\text{div}\mathbf{\sigma }+\mathbf{\text{f}}=0,$$where σ is the Cauchy stress. f represents an external body force per unit deformed volume. Herein, Cauchy stress includes the mechanical stress as well as the Maxwell stress. On an element of the deformed surface $\partial B_t$, the external surface traction is given by

2

$$\mathbf{\text{t}}(\mathbf{n})=[[\mathbf{\sigma }]]\mathbf{n},$$where the notation [[·]] denotes the jump operator, defined as the difference between the physical quantity inside and outside domain, i.e., $[[\mathbf{\sigma }]]={\mathbf{\sigma }}^{\text{in}}-{\mathbf{\sigma }}^{\text{out}}$ on boundary $\partial B_t$. n represents the unit outer normal vector.

By introducing the spatial electric displacement D, Gauss’s law then may be expressed spatially as follows [36]:

3

$$\text{div}\mathbf{D}=q,$$where q is the free charge density per unit deformed volume. Analogous to Eq. (2), the free charge density on the deformed surface is given by

4

$$w(\mathbf{n})=-[[\mathbf{D}]]·\mathbf{n},$$where $[[\mathbf{D}]]={\mathbf{D}}^{\text{in}}-{\mathbf{D}}^{\text{out}}$ on boundary $\partial B_t$.

Nominal quantities in reference configuration (i.e., referential electric displacement $\tilde{\mathbf{D}}$ and referential electric field $\tilde{\mathbf{E}}$) are related to true quantities in deformed configuration (i.e., spatial electric displacement D and spatial electric field E), by [37]

5

$${\tilde{D}}_i=J{F}_{\mathit{ij}}^{-1}{D}_j,\begin{array}{ccc}& & \end{array}{\tilde{E}}_i=F_{\mathit{ji}}{E}_j.$$

Constitutive theory

With the consideration of the nonlinearity of the material, the potential energy W per unit of reference volume is a function of the deformation gradient F and referential electric field $\tilde{\mathbf{E}}$. Further, we assume that $W\left(\mathbf{F},,,\tilde{\mathbf{E}}\right)$ can be written in a separable form of

6

$$W\left(\mathbf{F},,,\tilde{\mathbf{E}}\right)=W_m\left(\mathbf{F}\right)+W_c\left(\mathbf{F},,,\tilde{\mathbf{E}}\right)+W_e\left(\tilde{\mathbf{E}}\right).$$where W_m(F) is the change induced by mechanical deformation, $W_c\left(\mathbf{F},,,\tilde{\mathbf{E}}\right)$ corresponds to the coupled energy density, and $W_e\left(\tilde{\mathbf{E}}\right)$ represents the electric energy density, respectively.

Taking into account thermodynamics restrictions, the energy density determines the first Piola–Kirchhoff stress S and referential electric displacement $\tilde{\mathbf{D}}$ through the state relations [37]

7

$$\mathbf{S}=\frac{\partial W\left(\mathbf{F},\tilde{\mathbf{E}}\right)}{\partial \mathbf{F}},\begin{array}{ccc}& & \end{array}\tilde{\mathbf{D}}=-\frac{\partial W\left(\mathbf{F},\tilde{\mathbf{E}}\right)}{\partial \tilde{\mathbf{E}}}.$$

For simplicity, the anisotropy of flexible material is neglected. A Neo-Hookean material model is adopted in the form of

8

$$W_m(\mathbf{F})=\frac{\mu }{2}\left({\overline{I}}_1-3\right)+\frac{K_0}{2}{\left(J-1\right)}^2,$$where μ and K₀ are the initial shear and bulk moduli, respectively. ${\overline{I}}_1$ corresponds to the first invariant of the distortional right Cauchy–Green deformation tensor, defined by $J^{-\frac{2}{3}}{I}_1$ [25]. J represents the determinant of the deformation gradient.

By applying the chain rule, the first Piola–Kirchhoff stress $S_{\mathit{ij}}^m$ related to material deformation is given by

9

$$S_{\mathit{ij}}^m=\frac{\partial W_m}{\partial F_{\mathit{ij}}}=\mu J^{-\frac{2}{3}}\left(F_{\mathit{ij}}-\frac{1}{3}{I}_1{F}_{\mathit{ji}}^{-1}\right)+K_{0}J\left(J-1\right)F_{\mathit{ji}}^{-1}.$$

The corresponding Cauchy stress ${\sigma }_{\mathit{ij}}^m$ has the form of

10

$${\sigma }_{\mathit{ij}}^m=J^{-1}{S}_{\mathit{ik}}^m{F}_{\mathit{jk}}=\mu J^{-\frac{5}{3}}\left(F_{\mathit{ik}}{F}_{\mathit{jk}}-\frac{1}{3}{I}_1{\delta }_{\mathit{ij}}\right)+K_0\left(J-1\right){\delta }_{\mathit{ij}}.$$

In the linear theory of piezoelectricity [38], the relationship between stress and the electric field is linear, implying that reversing the electric field also reverses the stress. For this to hold true, the condition of [39]

11

$$\frac{\partial W_c\left(\mathbf{F},-\tilde{\mathbf{E}}\right)}{\partial \mathbf{F}}=-\frac{\partial W_c\left(\mathbf{F},\tilde{\mathbf{E}}\right)}{\partial \mathbf{F}}$$must be satisfied. Supposing the electric energy density can be represented as

12

$$W_c\left(\mathbf{F},\mathbf{E}\right)=\mathbf{A}\left(\mathbf{C}\right)\tilde{\mathbf{E}},$$where A(C) is a vector function of right Cauchy–Green deformation tensor C, defined in the reference configuration, and therefore automatically objective [39].

The referential electric displacement corresponding to the coupled energy density can be calculated by

13

$${\tilde{D}}_i^c=-\frac{\partial W_c}{\partial {\tilde{E}}_i}=-A_i.$$

When the material is in its undeformed state (i.e., C = I with I representing the identity tensor), the above expression should be evaluated to zero. One obtains

14

$$\mathbf{A}\left(\mathbf{I}\right)=0$$

The corresponding nominal stress is

15

$$S_{\mathit{ij}}^c=\frac{\partial W_c}{\partial F_{\mathit{ij}}}=\frac{\partial A_k}{\partial F_{\mathit{ij}}}{\tilde{E}}_k=\frac{\partial A_k}{\partial C_{\mathit{lm}}}\frac{\partial C_{\mathit{lm}}}{\partial F_{\mathit{ij}}}{\tilde{E}}_k=\frac{\partial A_k}{\partial C_{\mathit{lm}}}\left({\delta }_{\mathit{jl}}{F}_{\mathit{im}}+{\delta }_{\mathit{jm}}{F}_{\mathit{il}}\right){\tilde{E}}_k=\left(\frac{\partial A_k}{\partial C_{\mathit{jm}}}{F}_{\mathit{im}}+\frac{\partial A_k}{\partial C_{\mathit{lj}}}{F}_{\mathit{il}}\right){\tilde{E}}_k.$$

Assuming the piezoelectric tensor takes the form given by

16

$$e_{\mathit{ijk}}=-2\frac{\partial A_i}{\partial C_{\mathit{jk}}}=e_{\mathit{ikj}}$$

With the help of Eqs. (14) and (16), the vector function A(C) is constructed as a function of Green–Lagrange strain $\mathbf{\gamma }$, which is defined by $\mathbf{A}\left(\mathbf{C}\right)=-\mathbf{\gamma }e$. The strain tensor has the form of $\mathbf{\gamma }=\frac{1}{2}\left(\mathbf{C}-\mathbf{I}\right)$.

Substituting Eqs. (16)–(15) yields

17

$$S_{\mathit{ij}}^c=-e_{\mathit{kjm}}{F}_{\mathit{im}}{\tilde{E}}_k$$

Regarding the electric energy of $W_e\left(\tilde{\mathbf{E}}\right)$, assuming it depends on the invariant of $I_5=\tilde{\mathbf{E}}·\left({\mathbf{C}}^{-1},\tilde{\mathbf{E}}\right)$. This energy function takes the form of

18

$$W_e\left(\tilde{\mathbf{E}}\right)=-\frac{1}{2}\epsilon J{I}_5=-\frac{1}{2}\epsilon J{\tilde{E}}_i{\tilde{E}}_j{C}_{\mathit{ij}}^{-1}.$$where $\epsilon$ is the dielectric tensor.

Recall Eqs. (6), (8), (12), and (18), one obtains the total energy density

19

$$W(\mathbf{F},\widetilde{\mathbf{\text{E}}})=\frac{\mu }{2}\left({\overline{I}}_1-3\right)+\frac{K_0}{2}{\left(J-1\right)}^2-\frac{1}{2}\epsilon J{\tilde{E}}_i{\tilde{E}}_j{C}_{\mathit{ij}}^{-1}-e_{\mathit{ijk}}{\gamma }_{\mathit{jk}}{\tilde{E}}_i.$$

Straightforward calculations yield the first Piola–Kirchhoff stress as

20

$$S_{\mathit{ij}}=\mu J^{-\frac{2}{3}}\left(F_{\mathit{ij}}-\frac{1}{3}{I}_1{F}_{\mathit{ji}}^{-1}\right)+K_{0}J\left(J-1\right)F_{\mathit{ji}}^{-1}-\frac{1}{2}\epsilon J{F}_{\mathit{ji}}^{-1}{I}_5+\epsilon J{\tilde{E}}_p{\tilde{E}}_q{C}_{\mathit{pk}}^{-1}{C}_{\mathit{qj}}^{-1}{F}_{\mathit{ik}}-e_{\mathit{kjm}}{F}_{\mathit{im}}{\tilde{E}}_k,$$where the quadratic dependence of stress on the referential electric field has been introduced. Further, the Cauchy stress takes the form of

21

$$\begin{array}{cc}\hfill \sigma _{\mathit{ij}}^{}& =J^{-1}S_{\mathit{ik}}^{}{F}_{\mathit{jk}}=\mu J^{-\frac{5}{3}}\left(F_{\mathit{ik}}{F}_{\mathit{jk}}-\frac{1}{3}{I}_1{\delta }_{\mathit{ij}}\right)+K_0\left(J-1\right){\delta }_{\mathit{ij}}-\frac{1}{2}\epsilon E_m{E}_m{\delta }_{\mathit{ij}}+\frac{1}{2}\epsilon E_i{E}_j-J^{-1}{e}_{\mathit{lkm}}{F}_{\mathit{im}}{F}_{\mathit{jk}}{\tilde{E}}_l\hfill \\ \hfill \end{array}$$

The referential electric displacement can be expressed as follows:

22

$${\tilde{D}}_i=-\frac{\partial W}{\partial {\tilde{E}}_i}=-A_i+\epsilon J{\tilde{E}}_j{C}_{\mathit{ij}}^{-1}=e_{\mathit{ijk}}{\gamma }_{\mathit{jk}}+\epsilon J{\tilde{E}}_j{C}_{\mathit{ij}}^{-1}.$$

The spatial electric displacement D can be computed by

23

$$D_i=J^{-1}F_{\mathit{ij}}^{}{\tilde{D}}_j=J^{-1}F_{\mathit{ij}}^{}{e}_{\mathit{jkl}}{\gamma }_{\mathit{kl}}+\epsilon E_i,$$

Notably, the derivation of electromechanical energy density follows the references [25, 40, 41–42]. The electromechanical energy density of a dielectric elastomer can be obtained by setting the piezoelectric parameter to zero in Eq. (19), as described in the literature [25, 41]. In addition, Horák et al. [43] proposed a polyconvex phenomenological invariant-based formulation for simulating transversely isotropic electro-active polymers (EAPs) under large strains. This formulation provides a robust computational framework for accurately modeling the behavior of transversely isotropic EAPs, taking into account their anisotropic orientation and response to three-dimensional stretching, bending, and torsion. Furthermore, the twisting behavior of piezoelectric-type materials is also explored.

Uniaxial stretching theory and experiment

A schematic diagram of uniaxial stretching in the X₁-direction for a flexible piezoelectric membrane is shown in Fig. 1. Considering the scenario, the X₃-direction is the polarization direction, and a constant voltage ϕ is applied in the same direction. After loadings, the length, width, and thickness of the piezoelectric membrane change from L₁, L₂, and L₃ to l₁, l₂, and l₃, respectively. It is assumed that the piezoelectric membrane is completely incompressible, i.e., J = 1. The nonzero deformation gradients in principal directions are expressed as F₁₁ = l₁/L₁ = λ, F₂₂ = l₂/L₂ = λ^−1/2, and F₃₃ = l₃/L₃ = λ^−1/2.

[See PDF for image]

Fig. 1

Schematic illustration of the uniaxial stretching along the X₁-direction for a flexible piezoelectric membrane coated with electrodes on both sides, where an electric field is applied in the X₃-direction

The energy density for the flexible piezoelectric undergoing the uniaxial stretching is given by

24

$$W\left(\lambda ,{\tilde{E}}_3\right)=\frac{\mu }{2}\left({\lambda }^2+2{\lambda }^{-1}-3\right)-\frac{1}{2}{\epsilon }_{33}\lambda {\tilde{E}}_3^2-\frac{1}{2}{\tilde{E}}_3\left[e_{311}\left({\lambda }^2-1\right)+\left(e_{322}+e_{333}\right)\left({\lambda }^{-1}-1\right)\right].$$

The corresponding constitutive relation becomes

25

$$S=\mu \left(\lambda -{\lambda }^{-2}\right)-\frac{1}{2}{\epsilon }_{33}{\tilde{E}}_3^2-\frac{1}{2}{\tilde{E}}_3\left[2e_{311}\lambda -\left(e_{322}+e_{333}\right){\lambda }^{-2}\right],$$

26

$${\tilde{D}}_3={\epsilon }_{33}\lambda {\tilde{E}}_3+\frac{1}{2}\left[e_{311}\left({\lambda }^2-1\right)+\left(e_{322}+e_{333}\right)\left({\lambda }^{-1}-1\right)\right].$$

The uniaxial stretching experiments are performed using the tensile instrument (Instron 5943) as shown in Fig. 2. The flexible piezoelectric film made of polyvinylidene fluoride (PVDF, from MEAS, USA) has been polarized along the thickness direction and coated with 6-μm-thick ink electrodes on both sides. The dimensions of the PVDF sample film are 52 mm (length) × 15 mm (width) × 52 μm (thickness, excluding electrode thickness). Three PVDF samples are prepared for the uniaxial stretching tests. Notably, the contribution of the compliant silver electrodes to the mechanical response is neglected due to their softness. In addition, some environmental factors, such as temperature, humidity, and chemical exposure, have the potential to affect the material properties of PVDF. However, in the present work, our main focus is on investigating the influence of the electric field on the performance of flexible piezoelectric materials, with PVDF being considered mainly in its piezoelectric effects. Herein, we neglected the potential effects of other environmental factors on PVDF's performance.

[See PDF for image]

Fig. 2

Uniaxial stretching experiment of PVDF performed by a tensile instrument

In order to determine the shear modulus of PVDF, Eq. (25) is employed to fit the result of the uniaxial tensile experiment. As illustrated in Fig. 3, the fitting curve obtained from the equation is in good agreement with the experiments. The fitted shear modulus μ is found to be 0.78 GPa, which is close to the reference value reported by [44]. Particularly, for convenience, it is assumed that the dielectric constant ${\epsilon }_{11}$ is equal to ${\epsilon }_{22}$, the piezoelectric constant e₃₁₁ is equal to e₃₂₂, and e₁₁₃ is equal to e₂₂₃ [45].

[See PDF for image]

Fig. 3

Fitting result for the uniaxial stretching experiment of PVDF. Solid and dashed lines represent the experimental data and fitting results, respectively. Experiments 1, 2, and 3 correspond to the curves measured from different samples

All material parameters of PVDF used in this work are listed in Table 1, where piezoelectric constants and dielectric constants are obtained from the literature [46].

Table 1. Material properties of PVDF used in simulations

Shear modulus (GPa)	Piezoelectric constants (C/m2)			Dielectric constants (10–11 F/m)
μ	e311	e333	e113	${\epsilon }_{11}$	${\epsilon }_{33}$
0.78	0.11	-0.165	0.07	4.87	5.51

Finite element implementation

Neglecting body forces and volumetric free charge density, the governing partial differential equation for Cauchy stress in the deformed body B_t is expressed as follows [36]:

27

$$\text{div}\mathbf{\sigma }=0\begin{array}{cc}& \text{in}\begin{array}{c}\end{array}\begin{array}{c}{B}_t\end{array}\end{array}.$$

Gauss’s law for spatial electric displacement is

28

$$\text{div}\mathbf{D}=0\begin{array}{cc}& \text{in}\begin{array}{c}\end{array}\begin{array}{c}{B}_t\end{array}\end{array},$$where the expressions for $\mathbf{\sigma }$ and D are given by Eqs. (21) and (23), respectively.

The mechanical boundary conditions on the surface $\partial B_t$ can be defined in the following forms:

29

$$\mathbf{u}=\stackrel{⌣}{\mathbf{u}}\begin{array}{cc}& \text{on}\begin{array}{c}\end{array}\end{array}{S}_{\mathbf{u}},\text{and}\mathit{\sigma }\mathbf{n}=\stackrel{⌣}{\mathbf{t}}\begin{array}{cc}& \text{on}\begin{array}{c}\end{array}\end{array}{S}_{\mathbf{\text{t}}},$$where n is the outward unit normal on the boundary $\partial B_t$. $\stackrel{⌣}{\mathbf{u}}$ and $\stackrel{⌣}{\mathbf{\text{t}}}$ are the prescribed displacements and surface tractions, respectively. S_u and S_t denote the complementary subsurfaces of $\partial B_t$, satisfying $S_{\mathbf{u}}\cup S_{\mathbf{\text{t}}}=\partial B_t$ and $S_{\text{u}}\cap S_{\text{t}}=\mathrm{\varnothing }$. The symbols ∪ and ∩ stand for union and intersection operations, respectively.

To specify electric boundary conditions, we introduce another set of subsurfaces $S_{\varphi }$ and S_w. The electric potential is specified on $S_{\varphi }$ while spatial surface charge density on S_w, given by

30

$$\varphi =\stackrel{⌣}{\varphi }\begin{array}{cc}& \text{on}\begin{array}{c}\end{array}\end{array}{S}_{\varphi },\text{and}\mathbf{D}·\mathbf{n}=-\stackrel{⌣}{w}\begin{array}{cc}& \text{on}\begin{array}{c}\end{array}\end{array}{S}_w,$$with the prescribed functions $\stackrel{⌣}{\varphi }$ and $\stackrel{⌣}{w}$.

The deformed body is approximated using finite elements, $B_t=\cup B_t^e$. The displacement and electric potential serve as the basic variables. Within each element, these variables are interpolated using $u_i=\sum u_i^A{N}^A(i=1,2,3)$ and $\varphi =\sum {\varphi }^A{N}^A$. The superscript A = 1, 2, …8 denotes the number of nodes within an element, and N^A represents the shape functions.

A combination of Eqs. (27)–(30) represents the strong forms of the spatial boundary value problem for the displacement and electric potential. Notably, an essential step for finite element implementation is to derive the weak forms for the boundary value problem. By employing the Galerkin approach, one obtains

31

$$\begin{array}{cc}& {\int }_{B_t^e}\mathbf{\sigma }\frac{\partial N^A}{\partial \mathbf{x}}\text{d}V={\int }_{S_{\mathbf{\text{t}}}^e}{N}^A\stackrel{⌣}{\text{t}}\text{d}A,\hfill \\ \hfill & {\int }_{B_t^e}\mathbf{D}·\frac{\partial N^A}{\partial \mathbf{x}}\text{d}V+{\int }_{S_w^e}{N}^A\stackrel{⌣}{w}\text{d}A=0.\hfill \\ \hfill \end{array}$$

Define the element residuals for the displacements and electric potential as

32

$$\begin{array}{cc}\hfill R_{u_i}^A& =-{\int }_{B_t^e}{\sigma }_{\mathit{ij}}\frac{\partial N^A}{\partial x_j}\text{d}V+{\int }_{S_{\text{t}}^e}{N}^A{\stackrel{⌣}{t}}_i\text{d}A,\hfill \\ \hfill R_{\varphi }^A& ={\int }_{B_t^e}{D}_i\frac{\partial N^A}{\partial x_i}\text{d}V+{\int }_{S_w^e}{N}^A\stackrel{⌣}{w}\text{d}A,\hfill \\ \hfill \end{array}$$which are solved by Newton–Raphson’s procedure.

The tangents are defined in index notation as

33

$$\begin{array}{cc}\hfill K_{u_i^A{u}_k^B}& =-\frac{\partial R_{u_i}^A}{\partial u_k^B},K_{u_i^A{\varphi }^B}=-\frac{\partial R_{u_i}^A}{\partial {\varphi }^B},\hfill \\ \hfill K_{{\varphi }^A{u}_i^B}& =-\frac{\partial R_{\varphi }^A}{\partial u_i^B},K_{{\varphi }^A{\varphi }^B}=-\frac{\partial R_{\varphi }^A}{\partial {\varphi }^B}.\hfill \\ \hfill \end{array}$$

Evaluating the corresponding tangents in the absence of surface tractions and surface charge density, we have

34

$$\begin{array}{cc}\hfill K_{u_i^A{u}_k^B}& ={\int }_{B_t^e}\frac{\partial N^A}{\partial x_a}\left(J^{-1}{F}_{\mathit{aj}}{F}_{\mathit{mn}}\frac{\partial S_{\mathit{ij}}}{\partial F_{\mathit{kn}}}\right)\frac{\partial N^B}{\partial x_m}\text{d}V,\begin{array}{cc}& \end{array}{K}_{u_i^A{\varphi }^B}=-{\int }_{B_t^e}\frac{\partial N^A}{\partial x_a}\left(J^{-1}{F}_{\mathit{aj}}{F}_{\mathit{ml}}\frac{\partial S_{\mathit{ij}}}{\partial {\tilde{E}}_l}\right)\frac{\partial N^B}{\partial x_m}\text{d}V,\hfill \\ \hfill K_{{\varphi }^A{u}_i^B}& =-{\int }_{B_t^e}\frac{\partial N^A}{\partial x_a}\left(J^{-1}{F}_{\mathit{al}}{F}_{\mathit{jn}}\frac{\partial {\tilde{D}}_l}{\partial F_{\mathit{in}}}\right)\frac{\partial N^B}{\partial x_j}\text{d}V,\begin{array}{cc}& \end{array}{K}_{{\varphi }^A{\varphi }^B}={\int }_{B_t^e}\frac{\partial N^A}{\partial x_a}\left(J^{-1}{F}_{\mathit{al}}{F}_{\mathit{ji}}\frac{\partial {\tilde{D}}_l}{\partial {\tilde{E}}_i}\right)\frac{\partial N^B}{\partial x_j}\text{d}V.\hfill \\ \hfill \end{array}$$

Detailed derivations of these expressions can be found in the provided Appendix.

For the finite element analyses, the implementation has been carried out in ABAQUS using the user element subroutine. This newly developed three-dimensional hexahedral element consists of eight nodes, with each node having four degrees of freedom. In the UEL subroutine, RHS and AMATRX are necessary. The RHS matrix contains the residual vectors for each degree of freedom at each node and can be expressed as follows:

35

$$\text{R}={\left[R_{u_1}^1\begin{array}{c}\end{array}{R}_{u_2}^1\begin{array}{c}\end{array}{R}_{u_3}^1\begin{array}{c}\end{array}{R}_{\varphi }^1\begin{array}{c}\end{array}{R}_{u_1}^2\begin{array}{c}\end{array}{R}_{u_2}^2\begin{array}{c}\end{array}{R}_{u_3}^2\begin{array}{c}\end{array}{R}_{\varphi }^2\cdots R_{u_1}^8\begin{array}{c}\end{array}{R}_{u_2}^8\begin{array}{c}\end{array}{R}_{u_3}^8\begin{array}{c}\end{array}{R}_{\varphi }^8\right]}^{\text{T}}.$$

The stiffness matrix AMATRX is given by

36

$$\begin{array}{c}\hfill \mathbf{K}=\left[\begin{array}{ccccccccccccc}{K}_{u_1^1{u}_1^1}& K_{u_1^1{u}_2^1}& K_{u_1^1{u}_3^1}& K_{u_1^1\varphi _{}^1}^{}& K_{u_1^1{u}_1^2}& K_{u_1^1{u}_2^2}& K_{u_1^1{u}_3^2}& K_{u_1^1\varphi _{}^2}^{}& & K_{u_1^1{u}_1^8}& K_{u_1^1{u}_2^8}& K_{u_1^1{u}_3^8}& K_{u_1^1\varphi _{}^8}^{}\\ K_{u_2^1{u}_1^1}& K_{u_2^1{u}_2^1}& K_{u_2^1{u}_3^1}& K_{u_2^1\varphi _{}^1}^{}& K_{u_2^1{u}_1^2}& K_{u_2^1{u}_2^2}& K_{u_2^1{u}_3^2}& K_{u_2^1\varphi _{}^2}^{}& & K_{u_2^1{u}_1^8}& K_{u_2^1{u}_2^8}& K_{u_2^1{u}_3^8}& K_{u_2^1\varphi _{}^8}^{}\\ K_{u_3^1{u}_1^1}& K_{u_3^1{u}_2^1}& K_{u_3^1{u}_3^1}& K_{u_3^1\varphi _{}^1}^{}& K_{u_3^1{u}_1^2}& K_{u_3^1{u}_2^2}& K_{u_3^1{u}_3^2}& K_{u_3^1\varphi _{}^2}^{}& & K_{u_3^1{u}_1^8}& K_{u_3^1{u}_2^8}& K_{u_3^1{u}_3^8}& K_{u_3^1\varphi _{}^8}^{}\\ K_{{\varphi }^1{u}_1^1}& K_{{\varphi }^1{u}_2^1}& K_{{\varphi }^1{u}_3^1}& K_{{\varphi }^1\varphi _{}^1}^{}& K_{{\varphi }^1{u}_1^2}& K_{{\varphi }^1{u}_2^2}& K_{{\varphi }^1{u}_3^2}& K_{{\varphi }^1\varphi _{}^2}^{}& & K_{{\varphi }^1{u}_1^8}& K_{{\varphi }^1{u}_2^8}& K_{{\varphi }^1{u}_3^8}& K_{{\varphi }^1\varphi _{}^8}^{}\\ K_{u_1^2{u}_1^1}& K_{u_1^2{u}_2^1}& K_{u_1^2{u}_3^1}& K_{u_1^2\varphi _{}^1}^{}& K_{u_1^2{u}_1^2}& K_{u_1^2{u}_2^2}& K_{u_1^2{u}_3^2}& K_{u_1^2\varphi _{}^2}^{}& \cdots & K_{u_1^2{u}_1^8}& K_{u_1^2{u}_2^8}& K_{u_1^2{u}_3^8}& K_{u_1^2\varphi _{}^8}^{}\\ K_{u_2^2{u}_1^1}& K_{u_2^2{u}_2^1}& K_{u_2^2{u}_3^1}& K_{u_2^2\varphi _{}^1}^{}& K_{u_2^2{u}_1^2}& K_{u_2^2{u}_2^2}& K_{u_2^2{u}_3^2}& K_{u_2^2\varphi _{}^2}^{}& & K_{u_2^2{u}_1^8}& K_{u_2^2{u}_2^8}& K_{u_2^2{u}_3^8}& K_{u_2^2\varphi _{}^8}^{}\\ K_{u_3^2{u}_1^1}& K_{u_3^2{u}_2^1}& K_{u_3^2{u}_3^1}& K_{u_3^2\varphi _{}^1}^{}& K_{u_3^2{u}_1^2}& K_{u_3^2{u}_2^2}& K_{u_3^2{u}_3^2}& K_{u_3^2\varphi _{}^2}^{}& & K_{u_3^2{u}_1^8}& K_{u_3^2{u}_2^8}& K_{u_3^2{u}_3^8}& K_{u_3^2\varphi _{}^8}^{}\\ K_{{\varphi }^2{u}_1^1}& K_{{\varphi }^2{u}_2^1}& K_{{\varphi }^2{u}_3^1}& K_{{\varphi }^2\varphi _{}^1}^{}& K_{{\varphi }^2{u}_1^2}& K_{{\varphi }^2{u}_2^2}& K_{{\varphi }^2{u}_3^2}& K_{{\varphi }^2\varphi _{}^2}^{}& & K_{{\varphi }^2{u}_1^8}& K_{{\varphi }^2{u}_2^8}& K_{{\varphi }^2{u}_3^8}& K_{{\varphi }^2\varphi _{}^8}^{}\\ & & & ⋮& & & & & \ddots & & ⋮& & \\ K_{u_1^8{u}_1^1}& K_{u_1^8{u}_2^1}& K_{u_1^8{u}_3^1}& K_{u_1^8\varphi _{}^1}^{}& K_{u_1^8{u}_1^2}& K_{u_1^8{u}_2^2}& K_{u_1^8{u}_3^2}& K_{u_1^8\varphi _{}^2}^{}& & K_{u_1^8{u}_1^8}& K_{u_1^8{u}_2^8}& K_{u_1^8{u}_3^8}& K_{u_1^8\varphi _{}^8}^{}\\ K_{u_2^8{u}_1^1}& K_{u_2^8{u}_2^1}& K_{u_2^8{u}_3^1}& K_{u_2^8\varphi _{}^1}^{}& K_{u_2^8{u}_1^2}& K_{u_2^8{u}_2^2}& K_{u_2^8{u}_3^2}& K_{u_2^8\varphi _{}^2}^{}& \cdots & K_{u_2^8{u}_1^8}& K_{u_2^8{u}_2^8}& K_{u_2^8{u}_3^8}& K_{u_2^8\varphi _{}^8}^{}\\ K_{u_3^8{u}_1^1}& K_{u_3^8{u}_2^1}& K_{u_3^8{u}_3^1}& K_{u_3^8\varphi _{}^1}^{}& K_{u_3^8{u}_1^2}& K_{u_3^8{u}_2^2}& K_{u_3^8{u}_3^2}& K_{u_3^8\varphi _{}^2}^{}& & K_{u_3^8{u}_1^8}& K_{u_3^8{u}_2^8}& K_{u_3^8{u}_3^8}& K_{u_3^8\varphi _{}^8}^{}\\ K_{{\varphi }^8{u}_1^1}& K_{{\varphi }^8{u}_2^1}& K_{{\varphi }^8{u}_3^1}& K_{{\varphi }^8\varphi _{}^1}^{}& K_{{\varphi }^8{u}_1^2}& K_{{\varphi }^8{u}_2^2}& K_{{\varphi }^8{u}_3^2}& K_{{\varphi }^8\varphi _{}^2}^{}& & K_{{\varphi }^8{u}_1^8}& K_{{\varphi }^8{u}_2^8}& K_{{\varphi }^8{u}_3^8}& K_{{\varphi }^8\varphi _{}^8}^{}\end{array}\right]\end{array}$$

During the analysis, UEL is called for each iteration within a given increment. The subroutine receives the initial nodal coordinates of each element node, as well as the current guesses of the nodal displacements and electric potential. In addition, the nodal residuals and tangents (Eqs. (35) and (36)) are calculated for each element. To accommodate the nearly incompressible material behavior and mitigate volumetric locking behavior, the F-bar method proposed by Neto et al. [47] for fully integrated elements is implemented. In this method, a modified deformation gradient is constructed as $\overline{\mathbf{F}}={\left(det{\mathbf{F}}_0/det\mathbf{F}\right)}^{\frac{1}{3}}\mathbf{F}$, where F₀ is the deformation gradient at the centroid of the element. To account for the material behavior, Cauchy stress and spatial electric displacement at each Gauss point calculated by $\overline{\mathbf{F}}$ instead of deformation gradient F are taken into account, namely, $\mathbf{\sigma }=\widehat{\mathit{\sigma }}\left(\overline{\mathbf{F}}\right)$ and $\mathbf{D}=\widehat{\mathbf{D}}\left(\overline{\mathbf{F}}\right)$. The corresponding residuals in Eq. (32) are also calculated by $\overline{\mathbf{F}}$. However, tangents $K_{u_i^A{u}_k^B}$ and $K_{{\varphi }^A{u}_i^B}$ in Eq. (34) must be replaced by

37

$$\begin{array}{cc}\hfill K_{u_i^A{u}_k^B}& ={\int }_{B_t^e}\frac{\partial N^A}{\partial x_a}\left(\frac{1}{3}{\overline{J}}^{-1}{\overline{F}}_{\mathit{aj}}{\overline{F}}_{\mathit{rn}}\frac{\partial {\widehat{S}}_{\mathit{ij}}}{\partial {\overline{F}}_{\mathit{mn}}}{\delta }_{\mathit{mr}}{\delta }_{\mathit{bp}}-\frac{2}{3}{\widehat{\sigma }}_{\mathit{ia}}{\delta }_{\mathit{bp}}\right)\left(\frac{\partial N_0^B}{\partial x_k}-\frac{\partial N^B}{\partial x_k}\right)\text{d}V\hfill \\ \hfill & +{\int }_{B_t^e}\frac{\partial N^A}{\partial x_a}\left({\overline{J}}^{-1}{\overline{F}}_{\mathit{aj}}{\overline{F}}_{\mathit{ln}}\frac{\partial {\widehat{S}}_{\mathit{ij}}}{\partial {\overline{F}}_{\mathit{kn}}}\right)\frac{\partial N^B}{\partial x_l}\text{d}V,\hfill \\ \hfill K_{{\varphi }^A{u}_i^B}& ={\int }_{B_t^e}\frac{\partial N^A}{\partial x_a}\left(\frac{2}{3}{\widehat{D}}_a{\delta }_{\mathit{jk}}-\frac{1}{3}{\overline{J}}^{-1}{\overline{F}}_{\mathit{al}}{\overline{F}}_{\mathit{rn}}\frac{\partial {\underset{\sim }{\widehat{D}}}_l}{\partial {\overline{F}}_{\mathit{mn}}}{\delta }_{\mathit{mr}}{\delta }_{\mathit{jk}}\right)\left(\frac{\partial N_0^B}{\partial x_i}-\frac{\partial N_{}^B}{\partial x_i}\right)\text{d}V\hfill \\ \hfill & -{\int }_{B_t^e}\frac{\partial N^A}{\partial x_a}\left({\overline{J}}^{-1}{\overline{F}}_{\mathit{al}}{\overline{F}}_{\mathit{pn}}\frac{\partial {\underset{\sim }{\widehat{D}}}_l}{\partial {\overline{F}}_{\mathit{in}}}\right)\frac{\partial N^B}{\partial x_p}\text{d}V.\hfill \\ \hfill \end{array}$$

The details of derivations are shown in the Appendix.

Simulations

In this section, several representative numerical examples are presented to demonstrate the robustness of the developed user element subroutine. The following simulations employ the material properties listed in Table 1. If not specifically stated, it is assumed that the bulk modulus is 100 times the shear modulus, namely, K₀ = 10²μ, and a quasi-incompressible response is then effectively obtained.

Single-element verification

Firstly, single-element simulations of uniaxial stretching with an applied electric potential ϕ along the thickness direction are performed to validate the UEL. The analytical solution (Eq. (25)) under the complete incompressible condition is employed for comparison. Simulation results have also been compared with the corresponding widely accepted results from the built-in eight-noded linear piezoelectric brick element (C3D8E) in ABAQUS.

As illustrated in Fig. 4a, a hexahedron with length L₁, width L₂, and thickness L₃ is considered as the initial configuration. The direction of the applied electric field is consistent with the direction of the polarization. The stress is non-dimensionalized by the shear modulus μ, and the electric field is normalized by $\left(\varphi /L_3\right)\sqrt{{\epsilon }_{33}/\mu }$ as used in the literature [25].

[See PDF for image]

Fig. 4

a Schematic diagram of a flexible piezoelectric subjected to mechanical loading in the X₁-direction and electric loading in the X₃-direction. b Relationship between normalized stress S/μ and stretching ratio λ for uniaxial stretching with a zero constraint of electric potential. c Relation of normalized stress S/μ versus stretching ratio λ with a normalized applied electric field $\left(\varphi /L_3\right)\sqrt{{\epsilon }_{33}/\mu }=-0.1$. Dots and solid lines indicate the results from finite element analyses (red dots from developed user elements and black from C3D8E) and theoretical models, respectively. d Relationship between normalized stress S/μ and stretching ratio λ at varying normalized electric potentials ranging from − 0.1 to 0.1. (Color figure online)

In the case of L₁ = L₂ and L₂/L₃ = 10, the simulated response using the developed user element for uniaxial stretching with a zero constraint on electric potential is compared with the analytical solution (Eq. (25)) in Fig. 4b. The comparisons verify that the element accurately captures the nonlinear response of a Neo-Hookean material. In the small deformation scenario, the above results are in good agreement with those calculated by C3D8E. Notably, with the increase in the deformation, both the stress obtained from the developed user element and analytical solution are significantly smaller than that from C3D8E. This discrepancy arises due to the softening effects exhibited by flexible piezoelectric materials, which can be captured by the proposed constitutive model but not by the traditional linear constitutive model. Figure 4c demonstrates the response of stress and stretch ratio at a normalized electric potential of -0.1. When subjected to this electric field stimulation, the magnitude of stress is higher compared with the case without an electric field (Fig. 4b). In addition, the simulation result of the user element closely matches the theoretical solution, although they are slightly smaller than those from C3D8E. In Fig. 4d, the relationship between stress and stretch ratio is shown, under various normalized electric potentials of − 0.1, − 0.05, 0, 0.05, and 0.1. Initially, only electric loading is applied without displacement loading, resulting in stress due to the piezoelectric effect. It can be observed that the negative electric potential leads to higher stress, compared with the scenario with zero electric potential. This observation highlights the influence of the applied electric field on the material's electromechanical response.

Notably, the direct piezoelectric effect described in Eq. (26) can also be validated using the uniaxial stretching case. When the piezoelectric block is subjected to uniaxial tensile loads, an electric potential difference will be generated between its top and bottom surfaces. The electric potential difference can be obtained by $\varphi ={\tilde{E}}_3{L}_3$, where ${\tilde{E}}_3$ has the form of

38

$${\tilde{E}}_3=\frac{1}{{\epsilon }_{33}\lambda }\left({\tilde{D}}_3-\frac{1}{2}\left[e_{311}\left({\lambda }^2-1\right)+\left(e_{322}+e_{333}\right)\left({\lambda }^{-1}-1\right)\right]\right).$$

As shown in Fig. 5, the relationship between the normalized electric potential difference and stretch ratio is nonlinear. This nonlinearity presents an exciting opportunity for designing complex flexible piezoelectric sensors. The ability to manipulate and control the stretch ratio through the electric potential difference allows for the creation of sophisticated sensor designs. In addition, the analytical solution fits well with the numerical simulation, which validates the accuracy of the developed numerical framework.

[See PDF for image]

Fig. 5

Relation of generated normalized electric potential difference $\left(\varphi /L_3\right)\sqrt{{\epsilon }_{33}/\mu }$ and stretch ratio $\lambda$. Solid line and dots represent the results from finite element analysis and theoretical model, respectively

Piezoelectric sensor

Piezoelectric materials are commonly used in the design of sensors because they can convert mechanical deformation into an electric potential difference through the direct piezoelectric effect. As shown in Fig. 6a, a cubic block with a side length of L is subjected to compressive loadings along the X₃-direction. An electric potential difference will be generated between the upper and lower surfaces. This model has meshed with 1000 elements. Denote the normalized electric field by $\left(\varphi /L\right)\sqrt{{\epsilon }_{33}/\mu }$.

[See PDF for image]

Fig. 6

Schematic diagram and sensing response of flexible piezoelectric cube. a Schematic diagram of simple compression along the X₃-direction for a cube with a side length of L. b Relation of normalized electric field $\left(\varphi /L\right)\sqrt{{\epsilon }_{33}/\mu }$ versus normalized compressive displacement u₃/L. Solid and dashed lines represent the simulation results from the user element and C3D8E, respectively

Figure 6b illustrates the relationship between normalized electric field $\left(\varphi /L\right)\sqrt{{\epsilon }_{33}/\mu }$ and normalized compressive displacement u₃/L, where u₃ represents the total compressive displacement. When subjected to compressive load, the electric potential on the bottom surface of the material is higher than that on the top surface. For small compressive displacements, the normalized electric field exhibits a linear dependence on the normalized displacement. Furthermore, this finding is consistent with the results obtained from C3D8E, further validating the accuracy and reliability of our proposed model. However, with the increase in compression, the normalized electric field increases in an approximately logarithmical model. This nonlinear behavior is accurately captured by the developed user element subroutine, which accurately models the electric signal stimulated by mechanical effects. In contrast, the simulation result obtained from the C3D8E element is linear, failing to capture the nonlinear response. This capability to accurately model the nonlinear response of the electric field stimulated by mechanical effects is crucial in the design of smart sensors, where accurately characterizing the relationship between mechanical loading and the resulting electric field is essential.

Notably, the assumption of complete incompressibility (J = 1) may not hold true for all flexible piezoelectric materials. This assumption could introduce errors in simulation results. Figure 7 presents the sensor effects for flexible piezoelectric materials with varying bulk moduli, which correspond to different Passion ratios. As depicted in Fig. 7, a decrease in the bulk modulus is associated with a decrease in the generated electric potential difference. This observation suggests that the mechanical properties of the material, specifically the bulk modulus, have a significant influence on the resulting electric potential. Therefore, accurate characterization and consideration of the material's mechanical properties, including the bulk modulus, are crucial for obtaining reliable simulation results and designing effective flexible piezoelectric sensors.

[See PDF for image]

Fig. 7

Comparison of the generated normalized electric potential differences for a flexible piezoelectric material with different bulk moduli K₀

Piezoelectric actuator

Consider a simple piezoelectric actuator capable of converting electric work into mechanical motion. As shown in Fig. 8, a bilayer structure is composed of a piezoelectric material and an elastomeric substrate. It has a length of L₁, a width of L_2, and a thickness of 2L₃. A tie constraint is imposed between the contact surfaces of these two layers. A fixed displacement constraint is imposed on the left side of the structure, and an electric loading is applied in the X₁-direction of the piezoelectric layer. The substrate is taken to be an isotropic elastomer with a shear modulus μ_s. For the case of L₁/L₂ = 10 and L₂/L₃ = 10, the flexible piezoelectric and elastomeric substrates have been meshed with 750 developed user elements and 750 eight-noded linear brick elements (C3D8), respectively.

[See PDF for image]

Fig. 8

Schematic of the actuator structure consisting of a flexible piezoelectric and an elastomeric substrate

With modulus ratio μ_s/μ = 1 and normalized electric field $\left(\varphi /L_3\right)\sqrt{{\epsilon }_{33}/\mu }=0.3$, Fig. 9 illustrates the displacement and strain contours, as well as stress contours, for the flexible piezoelectric material and the elastomeric substrate.

[See PDF for image]

Fig. 9

With the normalized actuated electric field $\left(\varphi /L_3\right)\sqrt{{\epsilon }_{33}/\mu }=0.3$ and modulus ratio μ_s/μ = 1, displacement contours of (a) u₁ and (b) u₃, nominal strain contours of (c) NE₁₁ and (d) NE₃₃, and stress contours S₁₁ and S₃₃, along the X₁ and X₃-directions, respectively. Where (i) and (ii) represent the stress contours of flexible piezoelectric material and elastomeric substrate, respectively

Figure 9a and b shows the actuated displacement u₁ and u₃ along the X₁ and X₃-directions. In the deformed configuration, the top piezoelectric layer undergoes deformation at the applied electric field, causing a corresponding deformation in the elastomeric substrate that is tied to it. The maximum actuated displacement, u₃, occurs at the right end of the beam. This indicates that the right end experiences the highest displacement in the X₃-direction. Conversely, the maximum displacement in the X₁-direction is observed near the left end. Subsequently, the nominal strain contours are presented in Fig. 9c and d. The maximum nominal strain NE₁₁ is observed on the left side, which aligns with the maximum displacement u₁. Figure 9e and f illustrates the stress distributions along the X₁ and X₃-directions, represented by S₁₁ and S₃₃, respectively. In these figures, (i) and (ii) differentiate the results between the flexible piezoelectric material and the elastomeric substrate, respectively. It can be observed that the maximum stresses occur within flexible piezoelectric material and are higher than those within the substrate. This difference in stress levels can be attributed to the additional parameters associated with the piezoelectric material, which contribute to the overall stress within the structure. Despite assuming the same elastic parameters for the elastomeric substrate and the flexible piezoelectric material, the presence of piezoelectric properties leads to increased stress levels in the piezoelectric material.

In Fig. 10a, actuation performance for flexible piezoelectric materials with varying bulk moduli is investigated. A higher bulk modulus generally provides better actuation capabilities. On the other hand, a lower bulk modulus may result in slightly reduced actuation performance. However, this reduction is relatively small. This implies that even with a lower bulk modulus, the flexible piezoelectric material still retains a significant portion of its actuation capabilities. Figure 10b illustrates the relationship between normalized maximum actuated displacement and modulus ratio for the normalized electric field of 0.1, 0.2, and 0.3, respectively. A higher normalized electric field results in a larger actuated displacement. Compared with the harder elastomeric substrate, the presence of a softer one leads to a greater actuation response. The softer elastomeric substrate causes a greater actuation response. The relationship between the maximum actuated displacement and the modulus ratio follows an approximate negative exponential trend. As expected, when the substrate is significantly stiffer than the piezoelectric material, the actuation effect is not obvious. Notably, after a certain modulus ratio, further increases in the substrate modulus do not significantly reduce the actuation performance. The model described can be regarded as a simple piezoelectric device that converts electric stimulation into beam deflection for actuation.

[See PDF for image]

Fig. 10

Actuation performance under the applied electric field. a With the normalized actuated electric field $\left(\varphi /L_3\right)\sqrt{{\epsilon }_{33}/\mu }=0.3$ and modulus ratio μ_s/μ = 1, normalized actuated displacement u₃/L₃ across the structure length, where flexible piezoelectric material with different bulk moduli K₀ is considered. b Relations between maximum normalized actuated displacement and the shear modulus ratio μ_s/μ for the normalized electric field of 0.1, 0.2, and 0.3, respectively

Concluding remarks

In this work, based on the Neo-Hookean material model and electromechanical theory, we develop a large deformation model for flexible piezoelectric material. The model involves the quadratic dependence of stress on the electric field. To facilitate the structural design, the corresponding eight-noded solid elements have been developed in ABAQUS. To validate the accuracy of the developed elements, numerical examples of uniaxial stretching with various applied electric fields are performed. The results obtained from the developed user element under small deformations show good agreement with those obtained from the C3D8E. Moreover, the results for large deformations match well with the analytical solution. The current model reveals the stress-softening phenomenon encountered by flexible piezoelectric materials. Other numerical simulations are performed to demonstrate the robustness of the developed elements.

Taking into account that flexible materials often exhibit anisotropic and viscoelastic behavior, we acknowledge the importance of including these effects in our analysis. Moving forward, our future research will build upon the foundation established in this work and focus on incorporating the anisotropic and viscoelastic effects. Through these endeavors, our aim is to provide a more comprehensive understanding of the mechanical and electromechanical behavior of flexible piezoelectric materials. By incorporating anisotropic and viscoelastic effects, we will enhance the accuracy and applicability of our models, enabling more realistic predictions and analyses.

Acknowledgements

The authors would like to thank the supports from National Natural Science Foundation of China (NSFC) (12072150). This work is also supported by Joint Fund of Advanced Aerospace Manufacturing Technology Research (U1937601), Beijing Science and Technology Rising Star Program, National Natural Science Foundation of China for Creative Research Groups (No. 51921003).

Author Contribution

S. L. and Y. S. wrote the main manuscript. B. L. carried out the simulations and prepared Figs. 4–6. C. G. derived the constitutive model. All authors reviewed the manuscript.

Declarations

Conflict of interest

The authors declare no conflict of interest.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Li, T; Li, G; Liang, Y et al. Fast-moving soft electronic fish. Sci. Adv.; 2017; 3, 4 e1602045. [DOI: https://dx.doi.org/10.1126/sciadv.1602045]

2. Lyu, Q; Gong, S; Yin, J et al. Soft wearable healthcare materials and devices. Adv. Healthc. Mater.; 2021; 10, 17 2100577. [DOI: https://dx.doi.org/10.1002/adhm.202100577]

3. Chen, N; Tian, C. Design, modeling and testing of a 3-DOF flexible piezoelectric thin sheet nanopositioner. Sens. Actuator A Phys.; 2021; 323, 112660. [DOI: https://dx.doi.org/10.1016/j.sna.2021.112660]

4. Grosu, S; De Rijcke, L; Grosu, V et al. Driving robotic exoskeletons using cable-based transmissions: a qualitative analysis and overview. Appl. Mech. Rev.; 2018; [DOI: https://dx.doi.org/10.1115/1.4042399]

5. Yu, Y; Nassar, J; Xu, C et al. Biofuel-powered soft electronic skin with multiplexed and wireless sensing for human–machine interfaces. Sci. Robot.; 2020; 5, 41 eaaz7946. [DOI: https://dx.doi.org/10.1126/scirobotics.aaz7946]

6. Zhao, Z; Dai, Y; Dou, SX et al. Flexible nanogenerators for wearable electronic applications based on piezoelectric materials. Mater. Today Energy; 2021; 20, 100690. [DOI: https://dx.doi.org/10.1016/j.mtener.2021.100690]

7. Qiu, Y; Sun, S; Xu, C et al. The frequency-response behaviour of flexible piezoelectric devices for detecting the magnitude and loading rate of stimuli. J. Mater. Chem. C; 2021; 9, 2 pp. 584-594. [DOI: https://dx.doi.org/10.1039/D0TC02949A]

8. Kim, K; Kim, J; Jiang, X et al. Static force measurement using piezoelectric sensors. J Sensors; 2021; [DOI: https://dx.doi.org/10.1155/2021/6664200]

9. Gohari, S; Mozafari, F; Moslemi, N et al. Analytical solution of the electro-mechanical flexural coupling between piezoelectric actuators and flexible-spring boundary structure in smart composite plates. Arch. Civ. Mech. Eng.; 2021; 21, 1 pp. 1-25. [DOI: https://dx.doi.org/10.1007/s43452-021-00180-z]

10. Fang, D; Li, F; Liu, B et al. Advances in developing electromechanically coupled computational methods for piezoelectrics/ferroelectrics at multiscale. Appl. Mech. Rev.; 2013; [DOI: https://dx.doi.org/10.1115/1.4025633]

11. Choi, W; Kim, J; Lee, E et al. Asymmetric 2D MoS2 for scalable and high-performance piezoelectric sensors. ACS Appl. Mater. Interfaces; 2021; 13, 11 pp. 13596-13603. [DOI: https://dx.doi.org/10.1021/acsami.1c00650]

12. Luo, J; Zhang, L; Wu, T et al. Flexible piezoelectric pressure sensor with high sensitivity for electronic skin using near-field electrohydrodynamic direct-writing method. Extreme Mech. Lett.; 2021; 48, 101279. [DOI: https://dx.doi.org/10.1016/j.eml.2021.101279]

13. Dong, K; Wang, X. Influences of large deformation and rotary inertia on wave propagation in piezoelectric cylindrically laminated shells in thermal environment. Int. J. Solids Struct.; 2006; 43, 6 pp. 1710-1726. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2005.05.044]

14. Hu, KM; Li, H. Large deformation mechanical modeling with bilinear stiffness for Macro-Fiber Composite bimorph based on extending mixing rules. J. Intell. Mater. Syst. Struct.; 2020; 32, pp. 127-139. [DOI: https://dx.doi.org/10.1177/1045389X20951257]

15. Yi, S; Ling, SF; Ying, M. Large deformation finite element analyses of composite structures integrated with piezoelectric sensors and actuators. Finite Elem. Anal. Des.; 2000; 35, pp. 1-15. [DOI: https://dx.doi.org/10.1016/S0168-874X(99)00045-1]

16. Mukherjee, A; Chaudhuri, AS. Piezolaminated beams with large deformations. Int. J. Solids Struct.; 2002; 39, 17 pp. 4567-4582. [DOI: https://dx.doi.org/10.1016/S0020-7683(02)00341-4]

17. Guo, X; Sun, J; Li, L et al. Large deformations of piezoelectric laminated beams based on the absolute nodal coordinate formulation. Compos. Struct.; 2021; 275, 114426. [DOI: https://dx.doi.org/10.1016/j.compstruct.2021.114426]

18. Tani, J; Speziale, CG; Qiu, J. Intelligent material systems: application of functional materials. Appl. Mech. Rev.; 1998; 51, 8 pp. 505-521. [DOI: https://dx.doi.org/10.1115/1.3099019]

19. Yoon, K; Kelarakis, A. Nanoclay-directed structure and morphology in PVDF electrospun membranes. J. Nanomater.; 2014; [DOI: https://dx.doi.org/10.1155/2014/367671]

20. Zhang, J; Liu, D; Han, Q et al. Mechanically stretchable piezoelectric polyvinylidene fluoride (PVDF)/Boron nitride nanosheets (BNNSs) polymer nanocomposites. Compos. Part B Eng.; 2019; 175, 107157. [DOI: https://dx.doi.org/10.1016/j.compositesb.2019.107157]

21. Wang, Y; Hong, M; Venezuela, J et al. Expedient secondary functions of flexible piezoelectrics for biomedical energy harvesting. Bioact. Mater.; 2023; 22, pp. 291-311.

22. Joshi, SP. Non-linear constitutive relations for piezoceramic materials. J. Smart Mater. Struct.; 1992; 1, 1 80. [DOI: https://dx.doi.org/10.1088/0964-1726/1/1/012]

23. Wagner, UV; Hagedorn, P. Piezo-beam systems subjected to weak electric field: experiments and modeling of non-linearities. J. Sound Vib.; 2002; 256, 5 pp. 861-872. [DOI: https://dx.doi.org/10.1006/jsvi.2002.5024]

24. Tajeddini, V; Muliana, A. Nonlinear deformations of piezoelectric composite beams. Compos. Struct.; 2015; 132, pp. 1085-1093. [DOI: https://dx.doi.org/10.1016/j.compstruct.2015.06.041]

25. Henann, DL; Chester, SA; Bertoldi, K. Modeling of dielectric elastomers: Design of actuators and energy harvesting devices. J. Mech. Phys. Solids; 2013; 61, 10 pp. 2047-2066.3092108 [DOI: https://dx.doi.org/10.1016/j.jmps.2013.05.003]

26. Nestorović, T; Marinković, D; Shabadi, S; Trajkov, M. User defined finite element for modeling and analysis of active piezoelectric shell structures. Meccanica; 2014; 49, 8 pp. 1763-1774.3237966 [DOI: https://dx.doi.org/10.1007/s11012-014-9925-x]

27. Js, A; Ahaa, B. Architected cellular piezoelectric metamaterials: Thermo-electro-mechanical properties. Acta Mater.; 2019; 163, pp. 91-121. [DOI: https://dx.doi.org/10.1016/j.actamat.2018.10.001]

28. Serrao, PH; Kozinov, S. A novel 3D mixed finite element for flexoelectricity in piezoelectric materials. Int. J. Numer. Meth. Eng.; 2024; 4781405 [DOI: https://dx.doi.org/10.1002/nme.7500]

29. Gao, J; Cao, SQ. Second-order approximation of primary resonance of a disk-type piezoelectric stator for traveling wave vibration. Nonlinear Dyn.; 2010; 61, 4 pp. 591-603. [DOI: https://dx.doi.org/10.1007/s11071-010-9673-y]

30. Tan, XG; Vu-Quoc, L. Efficient and accurate multilayer solid-shell element: non-linear materials at finite strain. Int. J. Numer. Methods Eng.; 2005; 63, 15 pp. 2124-2170. [DOI: https://dx.doi.org/10.1002/nme.1360]

31. Tan, XG; Vu-Quoc, L. Optimal solid shell element for large deformable composite structures with piezoelectric layers and active vibration control. Int. J. Numer. Methods Eng.; 2010; 64, 15 pp. 1981-2013. [DOI: https://dx.doi.org/10.1002/nme.1433]

32. Schlögl, T; Leyendecker, S. Electrostatic-viscoelastic finite element model of dielectric actuators. Comput. Methods Appl. Mech. Eng.; 2016; 299, pp. 421-439.3434922 [DOI: https://dx.doi.org/10.1016/j.cma.2015.10.017]

33. Nandan, S., Sharma, D., Sharma, A.K.: Viscoelastic effects on the nonlinear oscillations of hard-magnetic soft actuators. J. Appl. Mech. (2023).

34. Sharma, AK; Joglekar, MM. Effect of anisotropy on the dynamic electromechanical instability of a dielectric elastomer actuator. Smart Mater. Struct.; 2019; 28, 1 015006. [DOI: https://dx.doi.org/10.1088/1361-665X/aaeaa5]

35. Kumar, A; Khurana, A; Sharma, AK et al. Dynamic analysis of anisotropic dielectric viscoelastomers incorporating humidity effect. J. Braz. Soc. Mech. Sci. Eng.; 2022; [DOI: https://dx.doi.org/10.1007/s40430-022-03646-0]

36. Wang, S; Decker, M; Henann, DL et al. Modeling of dielectric viscoelastomers with application to electromechanical instabilities. J. Mech. Phys. Solids; 2016; 95, pp. 213-229.3541870 [DOI: https://dx.doi.org/10.1016/j.jmps.2016.05.033]

37. Vu, DK; Steinmann, P; Possart, G. Numerical modelling of non-linear electroelasticity. Int. J. Numer. Methods Eng.; 2007; 70, 6 pp. 685-704.2316183 [DOI: https://dx.doi.org/10.1002/nme.1902]

38. Gao, CF; Noda, N. Faber series method for two-dimensional problems of an arbitrarily shaped inclusion in piezoelectric materials. Acta Mech.; 2004; 171, 1–2 pp. 1-13.

39. Dorfmann, L; Ogden, RW. Nonlinear Theory of Electroelastic and Magnetoelastic Interactions; 2014; Boston, Springer: [DOI: https://dx.doi.org/10.1007/978-1-4614-9596-3]

40. Sharma, AK; Joglekar, MM. A numerical framework for modeling anisotropic dielectric elastomers. Comput. Methods Appl. Mech. Eng.; 2019; 344, pp. 402-420.3873864 [DOI: https://dx.doi.org/10.1016/j.cma.2018.10.005]

41. Zhao, X; Hong, W; Suo, Z. Stretching and polarizing a dielectric gel immersed in a solvent. Int. J. Solids Struct.; 2008; 45, 14–15 pp. 4021-4031. [DOI: https://dx.doi.org/10.1016/j.ijsolstr.2008.02.023]

42. Wilson, ZA; Borden, MJ; Landis, CM. A phase-field model for fracture in piezoelectric ceramics. Int. J. Fract.; 2013; 183, 2 pp. 135-153. [DOI: https://dx.doi.org/10.1007/s10704-013-9881-9]

43. Horák, M; Gil, AJ; Ortigosa, R et al. A polyconvex transversely-isotropic invariant-based formulation for electro-mechanics: stability, minimisers and computational implementation. Comput. Methods Appl. Mech. Eng.; 2023; 403, 4502987 [DOI: https://dx.doi.org/10.1016/j.cma.2022.115695] 115695.

44. Lang, SB; Sollish, BD; Moshitzky, M et al. Model of a pvdf piezoelectric transducer for use in biomedical studies. Ferroelectr; 1980; 24, 1 pp. 289-292. [DOI: https://dx.doi.org/10.1080/00150198008238657]

45. Nan, CW; Li, M; Huang, JH. Calculations of giant magnetoelectric effects in ferroic composites of rare-earth–iron alloys and ferroelectric polymers. Phys. Rev. B; 2001; 63, 14 144415. [DOI: https://dx.doi.org/10.1103/PhysRevB.63.144415]

46. Omote, K; Ohigashi, H. Shear piezoelectric properties of vinylidene fluoride trifluoroethylene copolymer, and its application to transverse ultrasonic transducers. Appl. Phys. Lett.; 1995; 66, 17 pp. 2215-2217. [DOI: https://dx.doi.org/10.1063/1.113171]

47. Neto, E; Peri, D; Dutko, M et al. Design of simple low order finite elements for large strain analysis of nearly incompressible solids. Int. J. Solids Struct.; 1996; 33, 20–22 pp. 3277-3296.1400785 [DOI: https://dx.doi.org/10.1016/0020-7683(95)00259-6]