A Numerical Framework for Flexible–Electrical Coupled Analysis of Piezoelectric Structures with Large Deformations

1  Introduction

Piezoelectric smart structures exploit their inherent electromechanical coupling capability to enable both motion actuation and state sensing, making them highly attractive for active control and monitoring applications [1–3]. As a result, piezoelectric materials have been widely employed in flexible electronic devices and soft robotic actuators [4,5], sensing and vibration control of large-scale flexible aerospace appendages [6], as well as energy conversion and energy harvesting systems [7].

Conventional piezoelectric ceramics, such as lead zirconate titanate (PZT), exhibit high electromechanical efficiency and strong coupling coefficients. However, their intrinsic brittleness and limited strain capacity significantly restrict their applicability in highly deformable smart structures [8]. In contrast, flexible piezoelectric materials, including polyvinylidene fluoride (PVDF) and macro fiber composites (MFC), exhibit recoverable deformation, high flexibility, and superior fatigue resistance, making them more promising for applications involving large deformation, compliant actuation, and long-term durability [9,10]. In practical applications, however, these smart materials and structures are often subjected to pronounced geometric nonlinearities [11]. Accurate modeling and reliable prediction of the large-deformation behavior of smart structures, along with precise characterization of their sensing and actuation performance, are crucial for improving control accuracy and ensuring functional reliability. Consequently, the development of robust modeling approaches capable of capturing large-deformation kinematics, while simultaneously accounting for strong bidirectional electromechanical coupling, has become a key challenge in the analysis and design of flexible piezoelectric smart structures.

Existing electromechanical modeling approaches can generally be classified into two categories: equivalent force models [12] and numerical formulations based on the finite element method [13]. Regarding the equivalent force method, the piezoelectric effect is represented by distributed surface or body forces acting on the host structure. This approach has been widely adopted in vibration control and sensing analyses of smart beams, plates, and shells due to its simplicity and high computational efficiency. For instance, Liu et al. [14] modeled the control effect of piezoelectric actuators as equivalent bending moments and applied this process to the active control of nonlinear vibrations of the membrane. Lu et al. [15] adopted a comparable simplification strategy for PVDF actuators and investigated the influence of control parameters on vibration attenuation efficiency. However, because the equivalent force method does not explicitly account for the actual energy conversion process between the structural field and electric fields, it may lead to inaccurate predictions when dealing with large deformations and strong electromechanical coupling effects.

Another mainstream modeling approach incorporates electric potential degrees of freedom directly into conventional finite element formulations, thereby establishing electromechanical coupled finite element models [16]. Such methods enable a rigorous expression of the piezoelectric constitutive relations. For example, Ref. [17] employed quadrilateral shell elements to establish the model for piezoelectric laminated plates and curvilinear shell structures. Electromechanical coupling behavior and distributed active vibration control were subsequently investigated. Nevertheless, most conventional electromechanical coupled finite element models are developed under small deformation assumptions. As a result, their applicability remains limited for highly flexible structures undergoing large rotations and large strains, particularly in soft robotic structures, and lightweight aerospace appendages, where geometric nonlinearity and strong electromechanical coupling coexist.

The absolute nodal coordinate formulation (ANCF) has emerged as an effective modeling method for flexible structures with large deformation and large rotation, owing to its exact representation of nodal positions and slopes in the global coordinate system [18]. Recently, several studies have extended the ANCF framework to multiphysics modeling in order to investigate electromechanical behaviors by incorporating piezoelectric effects [19,20]. Typical approaches introduce scalar electric potential variables into the ANCF formulation to account for the influence of the electric field on the structural response. While these electromechanical models based on the ANCF formulation represent a significant step toward multiphysics coupling, they often rely on simplified electric potential descriptions [21]. The reciprocal effect of structural deformation on the electric field is not fully captured, limiting their ability to describe the complete bidirectional electromechanical coupling mechanism. Moreover, the displacement and electric fields are discretized using different shape functions or interpolation orders in many existing studies [22]. Such inconsistent discretization strategies inevitably introduce interpolation errors and data-mapping inaccuracies during the coupling process. These inconsistencies may degrade numerical accuracy and potentially lead to stability issues in strongly coupled nonlinear analyses, particularly when large deformations and strong electromechanical interactions coexist.

To address these challenges, unified mesh descriptions based on the ANCF framework have been proposed for multiphysics problems, demonstrating superior numerical accuracy and robustness by enforcing consistent high-order discretization across multiple physical fields [23,24]. Cui et al. [25,26] proposed a thermo-flexible coupled element based on the unified description, in which the displacement and temperature fields are discretized using the same mesh. This modeling approach was subsequently applied to the rigid–flexible–thermal coupled analysis of a large-aperture paraboloid antenna. Tian et al. [27–29] further advanced the unified-mesh modeling concept for multiphysics analysis and developed corresponding model order reduction techniques. These developments clearly demonstrate the effectiveness of unified interpolation strategies in multiphysics modeling. However, these developments have not yet been extended to establish fully bidirectional coupling between the structural and electric fields within the ANCF framework. In particular, existing studies rarely adopt an interpolation strategy in which the displacement and electric fields are discretized using identical shape functions. The absence of a unified discretization scheme may introduce field-mapping inconsistencies and thereby restrict the accurate description of voltage-induced deformation and nonlinear dynamic behavior under large-deformation conditions.

To address this limitation, the present study develops an integrated flexible–electrical coupling computational framework within the ANCF method. By employing a unified interpolation scheme for both displacement and electric fields, the proposed formulation effectively eliminates interpolation inconsistencies and enhances computational accuracy. Furthermore, the governing equations are derived from a unified energy functional and solved synchronously within a fully coupled nonlinear system. This strategy enables an accurate and comprehensive representation of bidirectional electromechanical interactions under large deformations. The remainder of this paper is organized as follows. Section 2 presents the integrated flexible–electrical coupled formulation. In Section 3, the governing equations of the coupled system are established and the corresponding solution procedures are developed. Section 4 provides several numerical examples to validate the proposed electromechanical coupling framework. Finally, the main conclusions are summarized in Section 5.

2  Integrated Flexible–Electrical Coupling Formulation Based on ANCF

This section develops an integrated flexible–electrical coupling formulation within the ANCF framework. The displacement field and the electric field are discretized using the same interpolation mesh and shape functions. Based on continuum mechanics and linear piezoelectric constitutive relations, the generalized elastic force vector and the electromechanical coupling matrices of the system are derived.

2.1 Unified Description of the Displacement and Electric Fields

Although fully parameterized solid elements provide general modeling capability, thin plates and membranes represent the dominant structural forms for thin piezoelectric materials such as PVDF films and MFC laminates. For such slender structures, thin plate elements offer a more efficient modeling strategy. Fig. 1 illustrates the flexible–electrical coupled ANCF thin plate element formulated in the global coordinate system. Specifically, Fig. 1a depicts the displacement field of the ANCF thin plate element, while Fig. 1b shows the electromechanical coupled field with a unified mesh description. All nodal coordinates of the element are defined in the global Cartesian coordinate system O-XYZ. In addition, a local coordinate system o-xyz is attached to the element to describe the position of an arbitrary point within the element. Each ANCF thin plate element consists of four nodal coordinate vectors, which are defined as follows:

eP=[ePT1,ePT2,ePT3,ePT4]T,(1)

in which

ePj=[rTj,r,xTj,r,yTj]T(j=1,2,3,4),(2)

where jeP denotes the nodal coordinate vector of the j-th node, r is the global position vector, r,x and r,y represent the gradient vectors with respect to the x and y directions, respectively.

Figure 1: Description of flexible–electrical coupled model: (a) ANCF thin plate element displacement field; (b) Electromechanical coupled field with a unified mesh description.

The displacement field of the ANCF thin element can be expressed in the global coordinate system as:

rP(x,y,t)=[r1r2r3]T=SP(ξ,η)eP(t),(3)

where ri, i = 1, 2, 3 represent the components of the position coordinates in the X, Y, and Z directions, and SP(ξ, η) is the shape function matrix of the displacement field, given by:

SP=[s1I3×3,s2I3×3,…,s12I3×3]∈R3×36,(4)

in which

{s1=−(ξ−1)(η−1)(2ξ2−ξ+2η2−η−1)s2=−lξ(ξ−1)2(η−1)s3=−wη(ξ−1)(η−1)2s4=ξ(η−1)(2ξ2−3ξ+2η2−η)s5=−lξ2(ξ−1)(η−1)s6=wξη(η−1)2s7=−ξη(2ξ2−3ξ+2η2−3η+1)s8=lξ2η(ξ−1)s9=wξη2(η−1)s10=η(ξ−1)(2ξ2−ξ+2η2−3η)s11=lξη(ξ−1)2s12=−wη2(ξ−1)(η−1),(5)

where l and w denote the length and width of the undeformed element, and ξ=x/l,η=y/w are dimensionless parameters.

The element mesh of the thin plate is defined on the neutral surface. Accordingly, the electric potential φ for each node represents the electric potential difference between the upper and lower electrodes at that location. The electric field is discretized using the same shape functions as those employed for the displacement field in the flexible–electrical coupled ANCF thin plate element. The electrical nodal coordinate vector eE is constructed from the nodal electric potential difference and its in-plane gradient components. Consequently, the generalized electrical nodal coordinates at each node are defined as follows:

eEj=[φj,φ,xjφ,yj]T(j=1,2,3,4).(6)

Thus, the electric field of the flexible–electrical coupled ANCF thin plate element is described by the following interpolation process:

φ(x,y,t)=SE(ξ,η)eE(t),(7)

where φ represents the electric potential difference at an arbitrary point within the coupled element, SE = [s1, s2, …, s12] ∈ ℝ1×12 is the shape function of the electric field, and eE is the electrical nodal coordinate vector.

By integrating the displacement and electric field formulations, the flexible–electrical coupled thin plate model with a unified mesh description is developed within the ANCF framework through Eqs. (3) and (7). Consequently, the global position vector and the electric potential difference at an arbitrary point of the element can be written as:

[rPφ]=[SP03×1201×36SE][ePeE]=Su(ξ,η)eu(t),(8)

where Su(ξ, η) and eu(t) denote the augmented shape function matrix and the unified nodal coordinate vector, respectively.

2.2 Flexible–Electrical Coupled Constitutive Model

In this work, isotropic linear piezoelectric constitutive relations are adopted and expressed as:

{σ=cε−eTED=eε+χE,(9)

where σ and ε are the stress and strain vectors of the coupled element, respectively, c is the elastic coefficient matrix, e represents the piezoelectric stress constant matrix, and D, E, and χ denote the electric displacement vector, electric field intensity vector, and dielectric permittivity matrix, respectively. It should be noted that, although most piezoelectric materials are anisotropic and non-centrosymmetric, an isotropic linear elastic stiffness matrix is adopted here for methodological clarity and numerical verification. While this simplification may introduce quantitative deviations for specific materials, it does not limit the generality of the present formulation.

The ANCF thin plate element is formulated based on Kirchhoff–Love plate theory, which assumes that the strain at an arbitrary point of the plate is related to the middle surface strain and changes in curvature and torsion [30], as follows:

ε=ε0+zκ,(10)

where ε0 and κ denote the Lagrangian strain tensor and curvature vector of the element mid-surface, respectively, which are defined by the following equations:

ε0=[εxxεyyεxy]T=[12(rP,xTrP,x−1)12(rP,yTrP,y−1)12(rP,xTrP,y)]T,(11)

κ=[κxxκyyκxy]T=[rP,xxTn‖n‖3rP,yyTn‖n‖3rP,xyTn‖n‖3]T,(12)

where n = rP,x × rP,y represents the normal vector of the neutral surface, and ‖•‖ denotes the norm of a vector.

According to continuum mechanics theory, the elastic coefficient matrix c can be expressed in terms of Lamé constants as:

c=[λ+2μλ0λλ+2μ0002μ],(13)

in which

λ=Eν(1+ν)(1−2ν),μ=E2(1+ν),(14)

where E is the elastic modulus of the material, and ν is Poisson’s ratio.

The relationship between the piezoelectric stress constant matrix e and the piezoelectric strain constant matrix d is given by:

e=dc.(15)

In the present study, the polarization direction of the piezoelectric material is assumed to be aligned with the thickness direction of the model. Besides, the electric field within the electromechanical coupled element is assumed to be uniform along the thickness direction. Accordingly, the electric field intensity E is obtained through the gradient operator as follows:

E=−∇φ=−[∂∂x∂∂y∂∂z]Tφ.(16)

Substituting Eq. (7) into the above equation yields:

E=−BEeE,(17)

in which

BE=[00SE/h]T∈R3×12.(18)

Based on the constitutive relations defined in Eq. (9), the strain energy of the flexible–electrical coupled ANCF thin plate element can be stated as:

U=12∫VεTσdV=12∫V[(ε0+zκ)Tc(ε0+zκ)−(ε0+zκ)TeTE]dV=12∫V(ε0Tcε0+z2κTcκ−ε0TeTE)dV.(19)

Furthermore, the work done by the electric force is given by:

WE=12∫VETDdV=12∫V(ETeε0+ETχE)dV.(20)

3  Integrated Modeling Framework for Flexible–Electrical Coupling and Its Solution Scheme

This section derives the governing equations of motion for flexible–electrical coupled multibody systems based on the Lagrange equations of the second kind. The presented coupled formulation systematically incorporates the structural displacement variables and the electrical variables into a unified set of nonlinear equations. In addition, numerical solution procedures for both static and dynamic flexible–electrical coupled problems are developed.

3.1 Governing Equations of Flexible–Electrical Coupled Dynamics

Based on the Lagrange equations of the second kind, we have the following dynamic equations for the displacement field and the electric field:

{ddt(∂L∂e˙P)−∂L∂eP=0ddt(∂L∂e˙E)−∂L∂eE=0,(21)

where L = T − U + WE + WF is the Lagrangian function.

The kinetic energy theorem is used to derive the element mass matrix. The velocity of an arbitrary infinitesimal volume element in the global coordinate system is r˙P=SPe˙P. Thus, the element kinetic energy is expressed as:

T=12∫Vr˙PTρr˙PdV=12∫V(SPe˙P)Tρ(SPe˙P)dV=12e˙PT(∫VSPTρSPdV)e˙P=12e˙PTMe˙P,(22)

where ρ is the material density, and M is the mass matrix of the element.

The work done by external forces can be written as:

Wext=∫VrPTFbdV+∫ArPTFsdA+rPTFc−∫AφqdA,(23)

where Fb, Fs, and Fc are the body force, surface force, and concentrated force vectors, respectively, q is the surface charge density applied to the coupled element, defined as:

q=CpeφcA=χ33φch,(24)

where Cpe = χ33A/h represents the capacitance of the piezoelectric material, and φc denotes the externally applied control voltage.

Next, taking partial derivatives of the Lagrangian L with respect to the displacement coordinates eP and velocities e˙P yields:

ddt(∂L∂e˙P)=Me¨P,(25a)

∂L∂eP=−∂U∂eP+∂WE∂eP+∂Wext∂eP=−∫V[(∂ε0∂eP)Tcε0+z2(∂κ∂eP)Tcκ−12(∂ε0∂eP)TeTE]dV+12∫VETe∂ε0∂ePdV+∫VSPTFbdV+∫ASPTFsdA+SPTFc,(25b)

in which

∂ε0∂eP=[∂εxx∂eP∂εyy∂eP∂εxy∂eP]=[SP,xTrP,xSP,yTrP,y12(SP,xTrP,y+rP,xTSP,y)],(26a)

∂κ∂eP=[∂κxx∂eP∂κyy∂eP∂κxy∂eP]=[SP,xxTn+rP,xxT∂n∂eP‖n‖3−3(rP,xxTn){(∂n∂eP)Tn}‖n‖5SP,yyTn+rP,yyT∂n∂eP‖n‖3−3(rP,yyTn){(∂n∂eP)Tn}‖n‖5SP,xyTn+rP,xyT∂n∂eP‖n‖3−3(rP,xyTn){(∂n∂eP)Tn}‖n‖5].(26b)

Similarly, taking derivatives of the Lagrangian L with respect to the electrical coordinates eE and their time derivatives leads to:

ddt(∂L∂e˙E)=0,(27a)

∂L∂eE=−∂U∂eE+∂WE∂eE+∂Wext∂eE=−∫V(−12ε0TeT∂E∂eE)dV+∫V[12(∂E∂eE)Teε0+(∂E∂eE)TχE]dV−∫ASETqdA,(27b)

where the partial derivative of the electric field intensity E with respect to the electrical coordinates eE is given by Eq. (17).

Substituting the Eqs. (25a)–(27b) into Eq. (21), one can derive the unconstrained equations of motion for the flexible–electrical coupled system as:

{Me¨P+Qs+QPE−QF=0QEP+QEE+Qφ=0,(28)

where M, Qs, QF, QEE, and Qφ represent the mass matrix, elastic force vector, external force vector, dielectric force vector, and electrical charge force associated with the applied electric field, respectively. The terms QPE and QEP denote the generalized electromechanical coupling force vectors. The specific expressions of these matrices and vectors are provided as follows.

M=∫VSPTρSPdV,(29a)

Qs=∫V[(∂ε0∂eP)Tcε0+z2(∂κ∂eP)Tcκ]dV,(29b)

QPE=−∫V(∂ε0∂eP)TeTEdV,(29c)

QF=∫VSPTFbdV+∫ASPTFsdA+SPTFc,(29d)

QEP=−∫V(∂E∂eE)Teε0dV,(29e)

QEE=−∫V(∂E∂eE)TχEdV,(29f)

Qφ=∫ASETqdA.(29g)

3.2 Multiphysics Solution Strategy Based on the Generalized-α Method

The Lagrange multiplier method is employed to establish the constrained governing equations of the flexible–electrical coupled system, which can be written as:

{Me¨P+Qs+QPE−QF+ΦP,qPTλP=0ΦP(eP,t)=0QEP+QEE+Qφ+ΦE,qETλE=0ΦE(eE,t)=0,(30)

where Φp denotes the dynamic constraint equation and ΦE represents the electrical constraint equation, λP and λE are the corresponding dynamic and electrical Lagrange multipliers, respectively.

In this work, the coupled governing equations of the system are solved by the generalized-α method. The computational flowchart of the solving algorithm is illustrated in Fig. 2. The temporal discretization of the global generalized coordinates can be expressed as:

{qn+1=qn+hq˙n+h2(0.5−β)anq˙n+1=q˙n+h(1−γ)anan+1=(αfq¨n−αman)/(1−αm)qn+1=qn+1+h2βan+1q˙n+1=q˙n+1+hγan+1q¨n=0,(31)

where q is global generalized coordinate vector, h is the time step, the subscripts n and n + 1 indicate the current and next time steps, respectively. αm and αf are parameters related to spectral radius ρ ∈ [0, 1], which are set as follows:

αm=2ρ−1ρ+1,αf=ρρ+1,β=14(1+αf−αm)2,γ=12+αf−αm.(32)

Figure 2: Computational flowchart of the solving scheme for the flexible–electrical coupled governing equations based on the generalized-α method.

Accordingly, the residual vector of the coupled system is written as:

G(qP,qE)=[GP(qP,qE)GE(qP,qE)ΦP(qP)ΦE(qE)]=[Mq¨P+Qs(qP)+QPE(qP,qE)−QF+ΦP,qPTλPQEP(qP,qE)+QEE(qE)+Qφ+ΦE,qETλEΦPΦE],(33)

in which

ΦP,qP=∂ΦP∂qP,ΦE,qE=∂ΦE∂qE.(34)

The Newton-Raphson method is then applied at each time step to solve the resulting system of nonlinear algebraic equations:

−JG[ΔqPΔqEΔλPΔλE]=−[∂GP∂qP∂GP∂qE∂GP∂λP∂GP∂λE∂GE∂qP∂GE∂qE∂GE∂λP∂GE∂λE∂ΦP∂qP∂ΦP∂qE∂ΦP∂λP∂ΦP∂λE∂ΦE∂qP∂ΦE∂qE∂ΦE∂λP∂ΦE∂λE][ΔqPΔqEΔλPΔλE]=G,(35)

where Jacobian matrix JG can be written as:

JG=[β^M+JQs+JQPE,PJQPE,EΦP,qPT0JQEP,PJQEP,E+JQEE0ΦE,qETΦP,qP0000ΦE,qE00],(36)

in which

β^=1−αmh2β(1−αf).(37)

The explicit expressions of the Jacobian matrices JQs, JQPE,P, JQPE,E, JQEP,P, JQEP,E, and JQEE are provided in Appendix A.

Upon convergence of the Newton–Raphson iteration, the generalized coordinates of the displacement field and the electrical field, together with the corresponding Lagrange multipliers, are updated according to the following relations:

{qP,n+1=qP,n+1+ΔqP,n+1q˙P,n+1=q˙P,n+1+γ^ΔqP,n+1q¨P,n+1=q¨P,n+1+β^ΔqP,n+1λP,n+1=λP,n+1+ΔλP,n+1qE,n+1=qE,n+1+ΔqE,n+1λE,n+1=λE,n+1+ΔλE,n+1,(38)

in which

γ^=γhβ.(39)

4  Numerical Examples

In this section, several numerical examples are presented to verify the accuracy and effectiveness of the proposed flexible–electrical coupled modeling and analysis framework. By comparing the numerical results obtained by the present formulation with those computed from the finite element software, the capability of the proposed approach in addressing multiphysics coupled problems is systematically validated.

4.1 In-Plane Deformation under External Electric Field

To verify the accuracy and convergence of the proposed flexible–electrical coupled ANCF thin plate element in electromechanical response analysis, the in-plane deformation of a piezoelectric patch subjected to a uniform surface charge is investigated, as shown in Fig. 3. The geometric dimensions and material properties of the rectangular piezoelectric plate are provided in Table 1. The material constants correspond to representative values of PVDF films, rather than those of a specific commercial product. The boundary conditions are specified such that one corner of the plate is clamped, while all remaining edges are free. Such configuration is intentionally selected to induce pronounced in-plane deformations along both the x- and y-directions under an electric field applied through the thickness. Points A and B on the patch are selected as the observation points to quantitatively evaluate the in-plane deformation response. A comparison test is conducted in ABAQUS using the C3D8E piezoelectric solid element, where each edge of the plate is discretized into 20 elements. The applied surface charge density is set to q = 0.01 C/m2.

Figure 3: In-plane deformation of the piezoelectric patch under the external electric field, where the electric field is applied along the Z direction.

Table 2 presents the variation of the in-plane displacement magnitudes with the mesh density at observation point A. The results clearly indicate that, with mesh refinement, the displacements predicted by the proposed method converge steadily and monotonically toward the reference solution. When a 6 × 6 mesh is adopted, the relative error with respect to the ABAQUS solution is less than 10%, which confirms the convergence and accuracy of the proposed flexible–electrical coupled formulation. Fig. 4 illustrates the spatial distribution of the in-plane displacement magnitude for different mesh resolutions. The displacement exhibits a smooth increase from the clamped corner toward the free edges. Different colors represent different displacement magnitudes. As expected, the displacement is zero at the clamped corner and gradually increases along both the x and y directions. The displacement contours are continuous and free of oscillations, indicating a stable and well-resolved numerical solution. Moreover, the displacement distributions predicted by the proposed method are in excellent agreement with the ABAQUS results, which further demonstrates that the proposed flexible–electrical coupled framework can accurately capture charge-driven in-plane electromechanical responses.

Figure 4: Contour plots of in-plane displacement distributions: (a) 2 × 2 elements; (b) 3 × 3 elements; (c) 4 × 4 elements; (d) 5 × 5 elements; (e) 6 × 6 elements; (f) ABAQUS results.

Furthermore, a parametric study is conducted to investigate the influence of surface charge density on the in-plane deformation response. The applied surface charge density is varied from 0.01 to 0.05 C/m2, and the corresponding in-plane displacements in the x- and y-directions predicted by the proposed method and ABAQUS at the observation points A and B are summarized in Table 3. Fig. 5 compares the results obtained from the two methods under different surface charge densities. It can be observed that the in-plane displacement increases approximately proportionally with increasing surface charge density, and the results obtained from the present formulation agree well with the ABAUQS reference solutions over the entire range of electrical loading. The displacement–charge curves exhibit an approximately linear relationship, which is consistent with the linear piezoelectric constitutive model adopted in this work. Besides, a small discrepancy between the two methods can be observed, and the deviation slightly increases as the surface charge density increases. This trend is attributed to the accumulation of discretization and numerical approximation errors under stronger electromechanical coupling conditions. Fig. 6 illustrates the in-plane deformation of the piezoelectric thin plate under a surface charge density of 1 C/m2, clearly indicating the large deformation characteristics induced by high electrical excitation. Overall, these comparative and parametric studies demonstrate the robustness and accuracy of the proposed formulation for charge-driven in-plane deformation analysis, further validating the predictive capability of the proposed flexible–electrical coupled element.

Figure 5: Comparison of in-plane displacements induced by surface charge density obtained by the proposed method and ABAQUS.

Figure 6: In-plane deformation under the surface charge density of 1 C/m2.

4.2 Flexible–Electrical Energy Conversion under Uniaxial Tension

To further verify the flexible–electrical conversion capability of the proposed piezoelectric coupled element, a piezoelectric plate subjected to uniaxial tensile loading is investigated to demonstrate that mechanical deformation induces an electrical response, as shown in Fig. 7. A rectangular piezoelectric plate with the same geometric and material parameters as those listed in Table 1 is considered. One side of the plate is fully clamped, while a uniformly distributed tensile force is applied along the opposite side in the in-plane direction. The applied tensile load σ is varied from 10 to 50 kPa to examine its influence on the induced electric potential distribution.

Figure 7: Schematic diagram of the piezoelectric patch under uniaxial tension, where Φ is the sensing electric potential difference, and σ is the applied tensile stress.

The Fig. 8 presents the spatial distributions of the induced electric potential under different tensile load levels. The electric potential contours exhibit smooth and continuous distributions over the plate surface. A clear potential gradient is observed along the loading direction, with the maximum potential magnitude occurring near the clamped boundary, where the strain is more pronounced. As the tensile load increases, the overall potential level increases significantly, while the spatial distribution pattern remains essentially unchanged, indicating a physically consistent electromechanical conversion behavior. To assess the accuracy of the proposed approach in capturing the inverse piezoelectric effect, the numerical results are compared with those obtained using ABAQUS with C3D8E elements. It can be observed that both methods exhibit similar electric potential magnitudes and spatial distributions. The two examples confirm that the proposed piezoelectric coupled element can accurately capture the bidirectional electromechanical coupling behavior, including both electrically induced mechanical deformations and mechanically induced electrical responses.

Figure 8: Sensing electric potential distributions under different external tensile levels: (a) σ = 10 kPa, present model; (b) σ = 10 kPa, ABAQUS result; (c) σ = 20 kPa, present model; (d) σ = 20 kPa, ABAQUS result; (e) σ = 30 kPa, present model; (f) σ = 30 kPa, ABAQUS result; (g) σ = 40 kPa, present model; (h) σ = 40 kPa, ABAQUS result; (i) σ = 50 kPa, present model; (j) σ = 50 kPa, ABAQUS result.

4.3 Dynamic Characteristics under Harmonic Voltage Excitation

The dynamic response of a cantilevered piezoelectric plate subjected to harmonic voltage excitation is investigated to assess the effectiveness of the proposed formulation in dynamic electromechanical analyses. The geometric and material parameters of the structure are identical to those listed in Table 1, with one side of the plate fully clamped. A harmonic voltage is applied across the electrodes, with a voltage amplitude of 10 kV and the excitation frequency of 10 Hz. Dynamic in-plane displacements at observation point A are shown in Fig. 9. An excellent agreement between the two solutions is observed in terms of both amplitude and phase. This comparison demonstrates the accuracy and effectiveness of the proposed modeling framework in capturing the dynamic electromechanical response of piezoelectric structures.

Figure 9: Time histories of dynamic in-plane displacements under harmonic voltage excitation. Displacement component in the (a) x-direction and (b) y-direction.

To further illustrate the spatial characteristics of the dynamic deformation, Fig. 10 depicts the displacement distributions of the piezoelectric patch at several instants within one period. The displacement contours clearly show the evolution of the in-plane deformation pattern over time, providing an intuitive visualization of the electrically induced dynamic responses. These results confirm that the proposed piezoelectric coupled element is capable of accurately predicting both the temporal and spatial characteristics of the dynamic response of piezoelectric structures.

Figure 10: In-plane displacement distributions and evolution of the dynamic deformation pattern: (a) x-direction; (b) y-direction.

4.4 Energy Conservation under Gravity

Next, the proposed flexible–electrical coupled modeling framework is further validated from the perspective of energy conversion by investigating the motion of a piezoelectric flexible pendulum under gravity. Table 4 provides the geometric and material parameters. One corner of the plate is constrained by a spherical joint, and the flexible pendulum is initially held stationary in the horizontal configuration. Under the action of gravity, the structure undergoes large-amplitude flexible motion, accompanied by the evolution of the induced electric potential due to the electromechanical coupling effect. Fig. 11 illustrates the deformed configurations of the flexible pendulum and the corresponding electric potential distributions during the motion. To quantitatively examine the energy transfer process, Fig. 12 presents the variation curves of kinetic energy EK, gravitational potential energy EG, elastic potential energy EE, and the total mechanical energy ET. It can be observed that the gravitational potential energy is gradually converted into kinetic energy and elastic potential energy during the falling process. Meanwhile, the total mechanical energy remains nearly constant, which indicates that the presented work accurately preserves the energy balance of the system.

Figure 11: Configurations and induced electric potential distributions of the flexible pendulum at different times.

Figure 12: Energy variation curves of the kinetic energy EK, gravitational potential energy EG, elastic potential energy EE, and total mechanical energy ET.

5  Conclusions

This work presents an integrated computational framework for modeling and nonlinear dynamic analysis of strongly coupled multiphysics systems with large deformation. Based on the ANCF formulation, a flexible–electrical coupled thin plate element is developed, in which the displacement and the electric fields are discretized with a unified mesh and incorporated into a coupled governing equation. To efficiently solve the nonlinear multiphysics system, a synchronized solution strategy based on the generalized-α method is further presented. Numerical examples systematically validate the effectiveness and robustness of the proposed framework for nonlinear analysis of piezoelectric structures. The results demonstrate that the present approach is capable of accurately capturing both the direct and inverse piezoelectric effects, as well as the bidirectional coupling mechanisms between the structural and the electric fields. In particular, the strong electromechanical coupling behavior under large deformation, together with energy conservation and conversion characteristics, is clearly revealed.

This study establishes a strongly coupled electromechanical modeling methodology that provides a unified theoretical foundation for future investigations of piezoelectric laminated structures, electrically induced bending, and active vibration control. These research directions facilitate systematic exploration of smart structures undergoing large deformations and their nonlinear dynamic behavior under coupled electromechanical excitation. It should be noted that the present formulation is developed under the classical thin-plate assumption with a linear isotropic elastic constitutive model, which defines the current scope of applicability. Nevertheless, the coupled multiphysics modeling framework is extensible and can be generalized to other plate and shell theories, anisotropic material models, and additional classes of finite elements.

Acknowledgement: The authors acknowledge the support of the China Postdoctoral Science Foundation (2025M774301).

Funding Statement: This work was supported by the China Postdoctoral Science Foundation (2025M774301), the Natural Science Foundation of Chongqing, China (CSTB2025NSCQ-GPX0430), the National Postdoctoral Researcher Program of China (GZC20252730), the Postdoctoral Research Startup Fund of Harbin Institute of Technology (AUGA5710027225), the Development and Construction Fund of Harbin Institute of Technology (CBQQ8880101025), and the Foreign Experts Fund of Harbin Institute of Technology (AUIQ8840004125).

Author Contributions: Xuan Sun proposed the methodology and drafted the main manuscript text. Yueying Zhu and Jiaxi Jin conducted data analysis. Zhitong Li reviewed the manuscript and provided supervision. Leizhi Wang performed visualization and validation. Zhaobo Chen was responsible for project administration and supervision. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data that support the findings of this study are available from the corresponding authors upon reasonable request.

Ethics Approval: Not applicable.

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

Nomenclature

Appendix A

The Jacobian matrices JQs, JQPE,P, JQPE,E, JQEP,P, JQEP,E, and JQEE are the partial derivatives of the corresponding generalized force vectors with respect to the element displacement nodal coordinates or electrical nodal coordinates. The expressions are given as follows:

JQs=∂Qs∂ep=∫V(∂2ε0T∂eP,i∂eP,jcε0+∂ε0T∂eP,ic∂ε0∂eP,j+z2∂2κT∂eP,i∂eP,jcκ+z2∂κT∂eP,ic∂κ∂eP,j)dV,(A1)

JQPE,P=∂QPE∂ep=−∫V∂2ε0T∂eP,i∂eP,jeTEdV,(A2)

JQPE,E=∂QPE∂eE=−∫V∂ε0T∂eP,ieT∂E∂eE,jdV,(A3)

JQEP,P=∂QEP∂ep=−∫V∂ET∂eE,ie∂ε0∂eP,jdV,(A4)

JQEP,E=∂QEP∂eE=−∫V∂2ET∂eE,i∂eE,jeε0dV,(A5)

JQEE=∂QEE∂eE=−∫V(∂2ET∂eE,i∂eE,jχE+∂ET∂eE,iχ∂E∂eE,j)dV,(A6)

where the detailed derivations of the second-order derivatives of the mid-surface Lagrangian strain tensor ε0, the curvature term κ and the electric field intensity E with respect to the generalized coordinates are provided below:

∂2ε0∂eP,i∂eP,j=[∂2εxx∂eP,i∂eP,j∂2εyy∂eP,i∂eP,j∂2εxy∂eP,i∂eP,j]=[SP,xiTSP,xjSP,yiTSP,yj12(SP,xiTSP,yj+SP,xjTSP,yi)],(A7)

∂2κ∂eP,i∂eP,j=[∂2κxx∂eP,i∂eP,j∂2κyy∂eP,i∂eP,j∂2κxy∂eP,i∂eP,j]T,(A8)

in which

∂2κxx∂eP,i∂eP,j=1‖n‖3(SP,xxiT∂n∂eP,j+SP,xxjT∂n∂eP,i+rP,xxT∂2n∂eP,i∂eP,j)−31‖n‖5[(SP,xxiTn+rP,xxT∂n∂eP,i)(∂nT∂eP,jn)+(SP,xxjTn+rP,xxT∂n∂eP,j)(∂nT∂eP,in)+(rP,xxTn)(∂nT∂eP,i∂eP,jn+∂nT∂eP,i∂n∂eP,j)]+151‖n‖7[(rP,xxTn)(∂nT∂eP,in)(∂nT∂eP,jn)],(A9)

∂2κyy∂eP,i∂eP,j=1‖n‖3(SP,yyiT∂n∂eP,j+SP,yyjT∂n∂eP,i+rP,yyT∂2n∂eP,i∂eP,j)−31‖n‖5[(SP,yyiTn+rP,yyT∂n∂eP,i)(∂nT∂eP,jn)+(SP,yyjTn+rP,yyT∂n∂eP,j)(∂nT∂eP,in)+(rP,yyTn)(∂nT∂eP,i∂eP,jn+∂nT∂eP,i∂n∂eP,j)]+151‖n‖7[(rP,yyTn)(∂nT∂eP,in)(∂nT∂eP,jn)],(A10)

∂2κxy∂eP,i∂eP,j=1‖n‖3(SP,xyiT∂n∂eP,j+SP,xyjT∂n∂eP,i+rP,xyT∂2n∂eP,i∂eP,j)−31‖n‖5[(SP,xyiTn+rP,xyT∂n∂eP,i)(∂nT∂eP,jn)+(SP,xyjTn+rP,xyT∂n∂eP,j)(∂nT∂eP,in)+(rP,xyTn)(∂nT∂eP,i∂eP,j+∂nT∂eP,i∂n∂eP,j)]+151‖n‖7[(rP,xyTn)(∂nT∂eP,in)(∂nT∂eP,jn)],(A11)

∂2E∂eE,i∂eE,j=0∈R12×12.(A12)

References

1. Hasim KA, Kefal A. A novel isogeometric layerwise element for piezoelectric analysis of laminated plates with straight/curvilinear fibers. Comput Meth Appl Mech Eng. 2022;399:115440. doi:10.1016/j.cma.2022.115440. [Google Scholar] [CrossRef]

2. Li J, Deng J, Zhang S, Liu Y. Development of a miniature quadrupedal piezoelectric robot combining fast speed and nano-resolution. Int J Mech Sci. 2023;250:108276. doi:10.1016/j.ijmecsci.2023.108276. [Google Scholar] [CrossRef]

3. Zhang Y, Huang Y, Zhao S, Jiao Z, Lü C, Yang J. Nonlinear dynamic response and stability of piezoelectric shells with piezoelectric nonlinearities. Int J Mech Sci. 2025;304:110731. doi:10.1016/j.ijmecsci.2025.110731. [Google Scholar] [CrossRef]

4. Zhang J, Liu D, Han Q, Jiang L, Shao H, Tang B, et al. Mechanically stretchable piezoelectric polyvinylidene fluoride (PVDF)/Boron nitride nanosheets (BNNSs) polymer nanocomposites. Compos Part B Eng. 2019;175:107157. doi:10.1016/j.compositesb.2019.107157. [Google Scholar] [CrossRef]

5. Tao B, Luo K, Tian Q, Hu H. Inverse dynamic design for motion control of soft machines driven by dielectric elastomer actuators. Int J Mech Sci. 2024;273:109213. doi:10.1016/j.ijmecsci.2024.109213. [Google Scholar] [CrossRef]

6. Zhang J, Fan M, Tzou H. Rigid–flexible coupled modeling and flexoelectric control for flexible rotating structures. AIAA J. 2024;62(12):4605–14. doi:10.2514/1.j064188. [Google Scholar] [CrossRef]

7. Yu D, Yang Y, Hu G, Zhou Y, Hong J. Energy harvesting from thermally induced vibrations of antenna panels. Int J Mech Sci. 2022;231:107565. doi:10.1016/j.ijmecsci.2022.107565. [Google Scholar] [CrossRef]

8. Wang G, Wang W, Liu C, Chen W, Mei Y. Extended locally exact homogenization theory for effective coefficients and localized responses of piezoelectric composites. Adv Eng Mater. 2022;24(5):2101194. doi:10.1002/adem.202101194. [Google Scholar] [CrossRef]

9. Lv S, Meng L, Zhang Q, Shi Y, Gao C. Numerical framework for anisotropic flexible piezoelectrics with large deformation. Int J Mech Sci. 2023;258:108564. doi:10.1016/j.ijmecsci.2023.108564. [Google Scholar] [CrossRef]

10. Gao C, Hu C, Gao S, Yu H. A high-order composite plate element for mechanical behavior of multi-drive coupled dielectric elastomer actuator via ANCF. Int J Appl Mech. 2026;18:2550122. doi:10.1142/s1758825125501224. [Google Scholar] [CrossRef]

11. Guo Y, Li L, Zhang D, Liao WH. Dynamic modeling and analysis for dielectric elastomer tube actuators. Int J Mech Sci. 2025;288:109994. doi:10.1016/j.ijmecsci.2025.109994. [Google Scholar] [CrossRef]

12. Liu X, Cai G, Peng F, Zhang H. Piezoelectric actuator placement optimization and active vibration control of a membrane structure. Acta Mech Solida Sin. 2018;31(1):66–79. doi:10.1007/s10338-018-0005-y. [Google Scholar] [CrossRef]

13. Aktaş KG, Esen İ. Dynamic response analysis and active vibration control of the smart sandwich composite plate with FGM core layers and MIMO FGPM actuators and sensors. Int J Mech Syst Dyn. 2025;5(1):3–19. doi:10.1002/msd2.70001. [Google Scholar] [CrossRef]

14. Liu X, Cai G, Peng F, Zhang H. Active control of large-amplitude vibration of a membrane structure. Nonlinear Dyn. 2018;93(2):629–42. doi:10.1007/s11071-018-4214-1. [Google Scholar] [CrossRef]

15. Lu Y, Shao Q, Amabili M, Yue H, Guo H. Nonlinear vibration control effects of membrane structures with in-plane PVDF actuators: a parametric study. Int J Non Linear Mech. 2020;122:103466. doi:10.1016/j.ijnonlinmec.2020.103466. [Google Scholar] [CrossRef]

16. Sharma A, Kumar R, Vaish R, Chauhan VS. Active vibration control of space antenna reflector over wide temperature range. Compos Struct. 2015;128:291–304. doi:10.1016/j.compstruct.2015.03.062. [Google Scholar] [CrossRef]

17. Balamurugan V, Narayanan S. Shell finite element for smart piezoelectric composite plate/shell structures and its application to the study of active vibration control. Finite Elem Anal Des. 2001;37(9):713–38. doi:10.1016/S0168-874X(00)00070-6. [Google Scholar] [CrossRef]

18. Shabana AA. Flexible multibody dynamics: review of past and recent developments. Multibody Syst Dyn. 1997;1(2):189–222. doi:10.1023/A:1009773505418. [Google Scholar] [CrossRef]

19. Yu X, You B, Wen J, Wei C, Zhao Y. Active vibration control and experimental investigations of the composite parabolic shell with variable cross-section under thermal load. Compos Struct. 2023;315:117011. doi:10.1016/j.compstruct.2023.117011. [Google Scholar] [CrossRef]

20. Guo Y, Yu Y, Li L, Zhang D, Liao WH, Du C, et al. Modeling and vibration control of a rotating flexible plate actuated by MFC. Compos Struct. 2024;331:117907. doi:10.1016/j.compstruct.2024.117907. [Google Scholar] [CrossRef]

21. Nada AA, El-Assal AM. Absolute nodal coordinate formulation of large-deformation piezoelectric laminated plates. Nonlinear Dyn. 2012;67(4):2441–54. doi:10.1007/s11071-011-0158-4. [Google Scholar] [CrossRef]

22. Nemov AS, Matikainen MK, Wang T, Mikkola A. Analysis of electromechanical systems based on the absolute nodal coordinate formulation. Acta Mech. 2022;233(3):1019–30. doi:10.1007/s00707-022-03153-2. [Google Scholar] [CrossRef]

23. Yu Z, Shen Y. Dynamics analysis of multibody systems based on flexible thermal coupling solid elements. Actuators. 2025;14(12):606. doi:10.3390/act14120606. [Google Scholar] [CrossRef]

24. Sun X, Jin J, Chen Z. Thermo-flexible coupled modeling and active control of thermally induced vibrations for a flexible plate. Thin Walled Struct. 2025;211:113092. doi:10.1016/j.tws.2025.113092. [Google Scholar] [CrossRef]

25. Cui YQ, Yu ZQ, Lan P. A novel method of thermo-mechanical coupled analysis based on the unified description. Mech Mach Theory. 2019;134:376–92. doi:10.1016/j.mechmachtheory.2019.01.001. [Google Scholar] [CrossRef]

26. Cui Y, Lan P, Zhou H, Yu Z. The rigid-flexible-thermal coupled analysis for spacecraft carrying large-aperture paraboloid antenna. J Comput Nonlinear Dyn. 2020;15(3):031003. doi:10.1115/1.4045890. [Google Scholar] [CrossRef]

27. Tian Q, Yu Z, Lan P, Cui Y, Lu N. Model order reduction of thermo-mechanical coupling flexible multibody dynamics via free-interface component mode synthesis method. Mech Mach Theory. 2022;172:104786. doi:10.1016/j.mechmachtheory.2022.104786. [Google Scholar] [CrossRef]

28. Tian Q, Sun H, Yu Z. Model order reduction of thermal-dynamic coupled flexible multibody system with multiple varying parameters. Appl Math Model. 2024;136:115634. doi:10.1016/j.apm.2024.08.006. [Google Scholar] [CrossRef]

29. Tian Q, Pan C, Liu Z, Chen X. Model order reduction for rigid–flexible–thermal coupled viscoelastic multibody system via the modal truncation with complex global modes. Actuators. 2025;14(10):479. doi:10.3390/act14100479. [Google Scholar] [CrossRef]

30. Amabili M. Nonlinear vibrations and stability of shells and plates. Cambridge, UK: Cambridge University Press; 2008. [Google Scholar]