1 Introduction↩︎

Dielectric elastomers (DEs) represent a class of electroactive polymers that exhibit significant deformation under electric fields, making them promising for actuators, sensors and energy harvesters. The theoretical foundations of electroelasticity were established by [1], who formulated a general theory for elastic dielectrics, followed by contributions from [2]–[4]. These works laid the groundwork for nonlinear continuum electromechanics, incorporating Maxwell stresses and polarisation effects. [5], [6] developed a comprehensive nonlinear theory for electroelasticity, introducing constitutive models based on free energy functions. This was extended to incompressible materials by [7], who derived variational principles for coupled problems. The modeling of DEs often involves hyperelastic potentials, such as the neo-Hookean and Gent models, to capture large deformations [8], [9].

Constitutive modeling for DEs has evolved to account for material nonlinearity and electromechanical coupling. [8] proposed a finite element formulation for EAPs using a free energy-based approach. [10] introduced a multiplicative formulation for nonlinear electro-elasticity, while [9] incorporated micromechanically motivated network models. For nearly incompressible materials, [11] addressed the challenges of volumeric locking by using mixed formulations. [12] extended these ideas to large deformations, employing augmented free energy functions. Experimental characterisation by [13] and [14] provided data for model validation, highlighting the importance of accurate permittivity and hyperelastic parameters.

The development of robust finite element methods for DEs has been subject to focus due to locking issues in thin structures. Early work by [15] adopted special stress elements for piezoelectric materials. [16] proposed a mixed formulation with six independent fields, and subsequently extended it to nonlinear dielectrics [17]. The Tangential Displacement Normal Normal Stress (TDNNS) method, introduced by [18], avoided shear and volume locking by using mixed elements with tangential displacement continuity. This was applied to piezoelectric solids [19] and later to large-deformation electro-elasticity [12]. For shells, [20] developed a Hellan-Herrmann-Johnson-type formulation, which [21] adapted to dielectric elastomer shells with independent thickness deformation. [22] and [23] provided foundations for nonlinear shell theory. For DEs, [24] presented a convex multi-variable potential for large strains, while  [17] developed a solid shell element with through-thickness electric field approximation. [12] introduced relaxed Kirchhoff-Love kinematics with independent thickness stretch, validated against 3D benchmarks. Applications include buckling actuators [12], peristaltic pumps [25], and spherical grippers [26], demonstrating the ability to capture complex instabilities.

Computational efficiency is a paramount concern in the numerical simulation of dielectric elastomers, as these materials exhibit complex behaviours like near-incompressibility and geometric nonlinearities that challenge conventional finite element methods. A significant contribution comes from [26], who introduced a novel framework employing Bézier elements within a mixed displacement-pressure formulation. This approach leverages monolithic solving strategies to overcome key limitations of traditional elements, such as Q1/P0 and F-bar elements [27], which often suffer from volumetric locking and poor convergence for incompressible materials. The framework’s effectiveness lies in its ability to maintain stability under large deformations.

Conventional finite element methods, based on Lagrange polynomials, often struggle with the Kirchhoff-Love shell formulation due to the requirement for \(C^1\)-continuous discretisations. This continuity condition ensures proper representation of bending effects without rotational degrees of freedom. To achieve \(C^1\)-continuity in Kirchhoff-Love shell discretisations, specialised approaches include exotic finite elements like Argyris spaces on triangles [28] and quadrilaterals [29], TUBA plate elements for adaptive \(C^1\) discretisation [30], as well as discontinuous Galerkin [31] and meshless methods [32], [33] for weak continuity imposition. These methods provide robust alternatives to conventional formulations while maintaining computational efficiency. Isogeometric analysis (IGA) [34] has emerged as a powerful alternative, offering smooth basis functions that naturally satisfy these continuity requirements. Recent advances in IGA Kirchhoff–Love shell formulations have addressed membrane locking through computationally efficient discretisations [35] and have extended the framework to trimmed multi-patch geometries using reduced-order methods [36]. Among IGA approaches, subdivision surfaces [37], [38] provide particularly attractive features, including the ability to handle complex geometries with arbitrary topology while maintaining the desired smoothness [39].

However, accurately capturing symmetry-breaking and post-buckling paths induced by electromechanical coupling instabilities in the numerical simulation of dielectric elastomer shells remains a challenge [40]–[42]. Recent studies have explored tunable morphing of DE balloons [43] and exploited instabilities for large shape transformations in DEs [44], yet a unified numerical framework capable of robustly tracing post-bifurcation paths in thin-shell geometries under combined electromechanical loading is still lacking. This requires high-fidelity computational models capable of handling large deformations, and geometric and material nonlinearities. The IGA analysis framework based on subdivision surfaces proposed in this paper, owing to its high-order continuity and accurate geometric representation, is particularly well-suited for simulating such nonlinear phenomena involving smooth deformation modes and complex instability patterns. This paper presents a comprehensive framework for the isogeometric analysis of electroelastic thin shells based on subdivision surfaces. The main contributions of this work are:

1. A nonlinear Kirchhoff–Love shell formulation specifically developed for dielectric elastomers, incorporating finite deformation kinematics and electromechanical coupling effects, with a consistent treatment of Maxwell stress and material incompressibility.

2. A systematic numerical framework for inducing bifurcation and tracing post-buckling equilibrium paths in electroelastic thin shells. This is achieved through a staged arc-length procedure combined with eigenmode perturbation, enabling the robust detection of symmetry-breaking instabilities and the stable traversal of unstable equilibrium branches.

3. The use of Catmull–Clark subdivision surfaces provides the \(C^1\)-continuity required by Kirchhoff–Love shell theory, ensuring a smooth representation of deformation fields even in the presence of severe localisation and self-contact during post-buckling.

4. Through comprehensive numerical examples, including spherical membranes, prestretched circular plates, and toroidal membranes, the proposed method is validated against analytical solutions, demonstrating its unique capability to capture bifurcation onset, mode switching, and post-buckling responses under coupled electromechanical loading.

The remainder of this paper is organised as follows. Section 2 establishes the theoretical foundation for electromechanical coupling in dielectric elastomers. Section 3 provides a comprehensive electroelastic shell formulation, addressing geometric description, kinematics, constitutive modelling, stress decomposition, incompressibility and plane stress constraints, stress resultants, consistent tangent moduli via static condensation, and the treatment of initial prestretch. Section 4 details the numerical implementation using subdivision surfaces. Section 5 describes the specialised techniques and algorithms developed for analysing bifurcation and post-buckling behaviour in electroelastic shells. Section 6 provides numerical examples that validate and showcase the proposed approach, and Section 7 concludes with a summary of key findings and future research directions.

Notations↩︎

Brackets↩︎

Square brackets \([ \, ]\) are used to group algebraic expressions. Round brackets \(( \,)\) are used to denote the dependencies of a function. If brackets are used to denote an interval, then \((\,)\) stands for an open interval and \([\,]\) is a closed interval. Curly brackets \(\{\,\}\) are used to define sets.

Symbols↩︎

A variable typeset in a normal weight font represents a scalar. A bold weight font denotes a first- or second-order tensor. An overline indicates that the variable is defined with respect to the reference configuration. If absent, the variable is defined with respect to the deformed configuration. A scalar variable with superscript or subscript indices normally represents the components of a vector or second-order tensor. Upright font is used to denote matrices and vectors.

Indices \(i,j,k,\dots\) vary from \(1\) to \(3\), while \(\alpha, \beta, \gamma,\dots\), used to indicate surface variable components, vary from \(1\) to \(2\). Einstein summation convention is used throughout.

The comma symbol in a subscript represents a partial derivative, for example, \(A_{,\beta}\) is the partial derivative of \(A\) with respect to the \(\beta^{\text{th}}\) coordinate.

To ensure clarity and avoid ambiguity in the nonlinear formulation and numerical implementation, the following notation is adopted throughout this work:

1. Layer Difference: \(\Delta\) denotes the physical difference between the upper and lower layers of the shell (e.g., \(\Delta \Phi = \Phi_{\text{top}} - \Phi_{\text{bottom}}\)).

2. Variation: \(\delta\) denotes the first variation of a variable (e.g., \(\delta \mathbf{u}\) is the virtual displacement).

3. Linearisation: \(\varDelta\) denotes the linearisation (total increment) over a solution step (e.g., \(\varDelta \mathbf{u}\) is the displacement increment for the current arc-length step).

4. Iterative Improvements: Symbols with a tilde (\(\tilde{\cdot}\)) are reserved for iterative updates within the nonlinear solver. Within a Newton–Raphson iteration, the iterative correction \(\tilde{\delta}\mathbf{u}\) is added to the accumulated increment: \[\tilde{\varDelta}\mathbf{u}^{(k+1)} = \tilde{\varDelta}\mathbf{u}^{(k)} + \tilde{\delta}\mathbf{u},\] so that the total increment \(\tilde{\varDelta}\mathbf{u}\) is the sum of all iterative improvements.

Coordinates↩︎

\(x\), \(y\), and \(z\) denote the Cartesian coordinates of a three-dimensional Euclidean space. \({\theta}^i\) denotes coordinates in the local element space. The three covariant basis vectors for a surface point are denoted as \(\mathbf{a}_i\), where \(\mathbf{a}_1, \mathbf{a}_2\) are tangential vectors and \(\mathbf{a}_3\) is the normal vector.

2 Electromechanical Coupling in Dielectric Elastomers↩︎

This section presents the three-dimensional electroelastic framework that serves as the foundation of the proposed shell formulation. The electrostatic governing equations and constitutive relations are first introduced in both the reference and current configurations. Subsequently, the transformation of electric field quantities under finite deformation is described, followed by the formulation of a general electroelastic strain-energy density function to characterise the coupling between mechanical deformation and electric fields.

2.1 Configurations↩︎

To analyse the electromechanical coupling problem involving a dielectric elastomer, it is necessary to carefully consider the governing equations and the relationship between the electric field and electric displacement in reference and deformed configurations. Consider the dielectric elastomer occupying a region of space \(\bar{\Omega}\) and \({\Omega}\) in \({\mathbb{R}}^3\) in its reference and deformed configuration, respectively (shown in Fig. 1). One introduces the deformation map \(\mathbf{r} = \chi(\mathbf{\bar{r}})\) to describe the motion from the reference to the deformed configuration, where \(\mathbf{r}\) and \(\mathbf{\bar{r}}\) are the position vectors in deformed and reference configurations, respectively. Also, the three-dimensional region \(\Omega\) lies inside region \(\mathcal{V}\), so that the surrounding free space is \({\Omega}^{\mathrm{'}} = \mathcal{V} \setminus \Omega \cup \partial \Omega.\) For \(\bar{\mathcal{V}}\) is the referential region corresponding to \(\mathcal{V}\) in \({\mathbb{R}}^3\), then \(\bar{\Omega}^{'}= \bar{\mathcal{V}} \setminus \bar{\Omega} \cup \partial \bar{\Omega}.\)

2.2 Electrostatic Governing Equations↩︎

The electrostatic problem for a dielectric elastomer is governed by Maxwell’s equations in the deformed domain \(\mathcal{V}\). The spatial electric field \(\mathbf{E}\) satisfying Faraday’s Law for electrostatics as \[\mathrm{\nabla} \times\,\mathbf{E}=\mathbf{0}.\] This implies that the electric field is irrotational and it can be expressed as the gradient of a scalar potential \(\Phi\) expressed as \[\mathbf{E}=-\mathrm{\nabla}\,\Phi. \label{eq:H-Phi95relation}\tag{1}\] Also, in the absence of volume charges, the spatial electric displacement \(\mathbf{D}\) is governed by the Gauss’s Law for Electricity as \[\mathrm{\nabla} \cdot\,\mathbf{D} = 0. \label{eq:Mag95field}\tag{2}\]

2.3 Constitutive Relationship↩︎

In a dielectric material, the electric displacement is related to the electric field by a constitutive equation as \[\mathbf{D}=\epsilon \mathbf{E} + \mathbf{P}, \label{eq:magneticc95constitutive}\tag{3}\] where \(\epsilon\) is the constant electric permittivity of free space and \(\mathbf{P}\) is the spatial polarisation, which vanishes in \(\Omega^{'}\).

2.4 Transformation Between Configurations↩︎

To relate the electric displacement and electric field in the reference and deformed configurations, the deformation gradient tensor \(\mathbf{F} := \partial\mathbf{r}/\partial\bar{\mathbf{r}}\), is adopted to perform a pull-back transformation, where \(\bar{\mathbf{r}}\) and \(\mathbf{r}\) are the position vectors of a material point in the reference and deformed configurations, respectively. Thus, the referential electric displacement \(\bar{\mathbf{D}}\) and electric field \(\bar{\mathbf{E}}\) (defined over the reference domain \(\bar{\mathcal{V}}\)) are given by \[\begin{align} \bar{\mathbf{D}} &=& \mathcal{J} \mathbf{F}^{-1} \mathbf{D}, \nonumber \\ \bar{\mathbf{E}} &=& \mathbf{F}^{\mathrm{T}}\mathbf{E}, \label{eq:transformations} \end{align}\tag{4}\] where \(\mathcal{J} = \mathrm{det}(\mathbf{F})\) represents the local volume change due to deformation. The electric field does not carry the Jacobian factor \(\mathcal{J}\) under the pull-back because it is a vector quantity that transforms covariantly with the deformation gradient. In contrast, the electric displacement is a flux density and must be scaled by \(\mathcal{J}\) to account for the change in cross-sectional area during deformation. These referential quantities satisfy the following Maxwell’s equations in the reference configuration: \[\begin{align} \bar{\nabla}\cdot\,\bar{\mathbf{D}} &= 0, \nonumber \\ \bar{\nabla}\times\,\bar{\mathbf{E}} &= \mathbf{0}, \end{align}\] where \(\bar{\nabla}\) denotes the gradient operator with respect to the reference coordinates \(\bar{\mathbf{r}}\). The electric field in the reference configuration can be expressed as the gradient of a scalar potential \(\bar\Phi\): \[\bar{\mathbf{E}} = -\bar{\nabla} \bar\Phi.\] The scalar potential \(\bar{\Phi} = \bar{\Phi}(\bar{\mathbf{r}})\) is the referential counterpart of the spatial potential \(\Phi = \Phi(\mathbf{r})\), obtained via the pull-back operation through the deformation map \(\mathbf{r} = \boldsymbol{\chi}(\bar{\mathbf{r}})\): \[\bar{\Phi}(\bar{\mathbf{r}}) = \Phi\big(\boldsymbol{\chi}(\bar{\mathbf{r}})\big). \label{eq:potential95map}\tag{5}\] This equality ensures that the potentials coincide at each material point, despite being expressed as functions of different coordinates.

Using the transformations 4 , the constitutive relationship 3 in the reference configuration is derived from its spatial form as \[\mathcal{J}^{-1} \mathbf{C}\,\bar{\mathbf{D}} = \epsilon\,\bar{\mathbf{E}} + \bar{\mathbf{P}} \quad \text{in } \bar{\mathcal{V}}, \label{eq:magneticc95constitutive95ref}\tag{6}\] where \(\bar{\mathbf{P}} = \mathbf{F}^{\mathrm{T}} \mathbf{P}\) and \(\mathbf{C} = \mathbf{F}^{\mathrm{T}}\mathbf{F}\) is the right Cauchy-Green deformation tensor. Since the polarization \(\mathbf{P}\) vanishes identically in the free space \(\Omega^{\prime}\) surrounding the material, it follows that \(\bar{\mathbf{P}} = \mathbf{F}^{\mathrm{T}}\mathbf{P} = \mathbf{0}\) in the corresponding referential region \(\bar{\Omega}^{\prime}\). Consequently, the constitutive relationship simplifies to \[\bar{\mathbf{D}} = \epsilon \mathcal{J} \mathbf{C}^{-1} \bar{\mathbf{E}} \quad \text{in } \bar{\Omega}^{\prime}. \label{free95sp95constitutive}\tag{7}\]

2.5 Strain Energy Density Function↩︎

For a general incompressible electroelastic solid one may regard the strain energy density function in the reference configuration as a function of the right Cauchy–Green tensor \(\mathbf{C}\) and the reference electric field vector \(\bar{\mathbf{E}}\) as:

\[W(\mathbf{C},\bar{\mathbf{E}}) = \widetilde{W}(I_1,I_2,I_3,I_4,I_5,I_6), \label{eq:energy95density}\tag{8}\] where \(\widetilde{W}\) depends on the six scalar invariants of \(\mathbf{C}\) and \(\bar{\mathbf{E}}\). The first three invariants are purely mechanical: \[\begin{align} I_1 &= \operatorname{tr}\mathbf{C}, & I_2 &= \frac{1}{2}\big[(\operatorname{tr}\mathbf{C})^2 - \operatorname{tr}(\mathbf{C}^2)\big], & I_3 &= \det\mathbf{C} = \mathcal{J}^2, \label{eq:mechanical95invariants} \end{align}\tag{9}\] The remaining invariants involve the electric field: \[\begin{align} I_4 &= \bar{\mathbf{E}} \cdot \bar{\mathbf{E}}, & I_5 &= \bar{\mathbf{E}} \cdot [\mathbf{C}\bar{\mathbf{E}}], & I_6 &= \bar{\mathbf{E}} \cdot [\mathbf{C}^2\bar{\mathbf{E}}], \label{eq:electrical95invariants} \end{align}\tag{10}\] where \(I_4\) captures purely electric effects, while \(I_5\) and \(I_6\) represent electromechanical coupling. For incompressible materials, \(\mathcal{J} \equiv 1\), and consequently \(I_3 \equiv 1\). A common simplifying assumption [5], [9], adopted in the present work, is that the mechanical and electrical contributions are additively separable: \[\widetilde{W}= \widetilde{W}_{\mathrm{mech}} + \widetilde{W}_{\mathrm{elec}}. \label{eq:energy95decomposed}\tag{11}\]

2.5.0.1 Electric energy for voltage‑controlled elastomers

In the present work, the dielectric material is assumed to be isotropic and linearly polarizable. Under voltage control, the electric energy density [9] is most naturally written in the current configuration as \[\widetilde{W}_{\mathrm{elec}}(\mathbf{C},\bar{\mathbf{E}}) = -\frac{1}{2}\,\epsilon\,\mathbf{E}\cdot\mathbf{E} = -\frac{1}{2}\,\epsilon\, \bigl[\bar{\mathbf{E}}\otimes\bar{\mathbf{E}}\bigr]:\mathbf{C}^{-1}. \label{eq:elec95energy951}\tag{12}\] Using the Cayley–Hamilton theorem for incompressible materials, this expression can be rewritten entirely in terms of the invariants: \[\widetilde{W}_{\mathrm{elec}} = -\frac{1}{2}\,\epsilon\, \bigl[I_6 - I_1 I_5 + I_2 I_4\bigr]. \label{eq:elec95energy95invariants}\tag{13}\] For comprehensive derivations, refer to 8.

3 Electroelastic Shells↩︎

The analysis of electroelastic thin shells presents several challenges that fundamentally distinguish it from purely mechanical shell problems. First, the presence of an electric field introduces additional stress contributions, namely Maxwell stresses, which are strongly coupled with mechanical deformation and give rise to pronounced nonlinear electromechanical interactions [16], [17]. Second, the thin-shell geometry necessitates a careful treatment of through-thickness kinematics and the plane stress condition, particularly under large deformations where thickness stretching becomes non-negligible. Third, the incompressibility constraint characteristic of dielectric elastomers must be incorporated consistently within the shell framework to ensure physically admissible deformation states.

To address these challenges, this section develops a comprehensive electroelastic shell formulation based on the Kirchhoff–Love hypothesis. The formulation systematically covers the geometric description and kinematics of the shell, energetic principles, stress decomposition, plane stress reduction, stress resultants, and the derivation of consistent tangent moduli. Finally, the framework is extended to incorporate the effects of initial prestretch, enabling the analysis of prestrained electroelastic thin-shell structures.

3.1 Geometric Description↩︎

Consider a shell made of dielectric elastomer in its reference configuration as the physical domain \(\bar\Omega \subset \mathbb{R}^3\), as shown in Fig. 2. The Kirchhoff–Love shell formulations is adopted to describe the mechanical behaviour of thin shell structures, which assumes that the lines that are perpendicular to the mid-surface before deformation remain straight after deformation. Each shell point \(\bar{\mathbf{r}} \in \bar\Omega\) is mapped from the parametric domain defined by the coordinates \(\{ \theta^1, \theta^2, \theta^3\}\). Assuming the shell has a uniform thickness \(\bar h\) in the reference configuration, the point \(\bar{\mathbf{r}}\) can be defined using a point on the mid-surface \(\bar\Gamma\), denoted \(\bar{\mathbf{x}} \in \bar\Gamma\), and the associated unit normal vector \(\bar{\mathbf{n}}\) as \[\bar {\mathbf{r}} (\theta^1,\theta^2,\theta^3)= \bar{\mathbf{x}}(\theta^1,\theta^2) + \theta^3 \bar{\mathbf{n}}(\theta^1,\theta^2), \label{eq:32r32x32n32relation}\tag{14}\] where \(\theta^3 \in [-{\bar h}/{2}, {\bar h}/{2}]\).

Both the reference and deformed configurations of the shell mid-surface are mapped from the mid-surface of the parametric domain. The corresponding mid-surface points in the reference and deformed configurations are denoted by \(\bar{\mathbf{x}}\) and \(\mathbf{x}\), respectively. The position vector of a mid-surface point in the deformed configuration, \(\mathbf{x}\), is related to its counterpart in the reference configuration, \(\bar{\mathbf{x}}\), through \[\mathbf{x} = \bar{\mathbf{x}} + \mathbf{u}, \label{eq:reference95to95deformed}\tag{15}\] where \(\mathbf{u}\) denotes the displacement vector of the mid-surface. Moreover, the covariant basis vectors on the mid-surface in the reference and deformed configurations are defined by \[\bar{\mathbf{a}}_\alpha = \frac{\partial \bar{\mathbf{x}}}{\partial \theta^\alpha} \quad \text{and} \quad {\mathbf{a}}_\alpha = \frac{\partial {\mathbf{x}}}{ \partial \theta^\alpha}.\] The corresponding unit normal vectors are then given by \[\bar{\mathbf{n}} = \bar{\mathbf{a}}^3 = \frac{\bar{\mathbf{a}}_1 \times \bar{\mathbf{a}}_2}{\bar{J}} \quad \text{and} \quad {\mathbf{n}} = {\mathbf{a}}^3 = \frac{{\mathbf{a}}_1 \times {\mathbf{a}}_2}{J},\] where \(\bar{J}\) and \(J\) denote the mid-surface Jacobians in the reference and deformed configurations, respectively, defined by \[\bar{J} = |\bar{\mathbf{a}}_1 \times \bar{\mathbf{a}}_2| \quad \text{and} \quad J = |{\mathbf{a}}_1 \times {\mathbf{a}}_2|.\]

The covariant components of the metric tensors for the mid-surface points \(\bar{\mathbf{x}}\) and \(\mathbf{x}\) are respectively given by \[\bar{a}_{ij} = \bar{\mathbf{a}}_i \cdot \bar{\mathbf{a}}_j \quad \text{and} \quad {a}_{ij} = {\mathbf{a}}_i \cdot {\mathbf{a}}_j.\] The corresponding contravariant metric tensors \(\bar{a}^{ik}\) and \(a^{ik}\) are defined by \[\bar{a}^{ik}\bar{a}_{kj} = \delta^{i}_{j} \quad \text{and} \quad a^{ik}a_{kj} = \delta^{i}_{j}, \label{eq:co95and95contra95metric}\tag{16}\] where \(\delta^{i}_{j}\) denotes the Kronecker delta.

The thickness stretch \(\lambda_3\) for a finitely deformed shell is defined by \[\lambda_3 = \frac{h}{\bar h},\] where \(h(\theta^1,\theta^2)\) is the shell thickness in the deformed configuration. We introduce a vector \(\mathbf{d}\) combining the thickness stretch and normal vector as \[\mathbf{d} = \lambda_3 \mathbf{a}_3,\] to write the position vector \(\mathbf{r}\) of a point in the deformed configuration of the shell-space as \[{\mathbf{r}}(\theta^1,\theta^2,\theta^3) = {\mathbf{x}}(\theta^1,\theta^2) + \theta^3 {\mathbf{d}}(\theta^1,\theta^2).\] Thus, the three-dimensional covariant basis vectors in the shell-space of the reference and the deformed configurations, respectively, follow as \[\bar{\mathbf{g}}_\alpha = \frac{\partial \bar{\mathbf{r}}}{\partial \theta^\alpha} = \bar{\mathbf{a}}_\alpha + \theta^3 \bar{\mathbf{a}}_{3,\alpha}, \quad \bar{\mathbf{g}}_3 = \frac{\partial \bar{\mathbf{r}}}{\partial \theta^3} = \bar{\mathbf{a}}_3, \label{eq:cov95tensors95r}\tag{17}\] and \[{\mathbf{g}}_\alpha = \frac{\partial {\mathbf{r}}}{\partial \theta^\alpha} = {\mathbf{a}}_\alpha + \theta^3 \mathbf{d}_{,\alpha}, \quad {\mathbf{g}}_3 = \frac{\partial {\mathbf{r}}}{\partial \theta^3} = {\mathbf{d}}. \label{eq:cov95tensors95d}\tag{18}\]

The components of the covariant metric tensors in the shell-space are given by \[\bar{g}_{ij} = \bar{\mathbf{g}}_i \cdot \bar{\mathbf{g}}_j\quad \text{and} \quad g_{ij} = \mathbf{g}_i \cdot \mathbf{g}_j,\] and the contravariant components of the metric tensor at point \(\mathbf{r}\) follow as \[\bar{g}^{ij} = \bar{\mathbf{g}}^{i} \cdot \bar{\mathbf{g}}^{j} \quad \text{and} \quad {g}^{ij} = {\mathbf{g}}^{i} \cdot {\mathbf{g}}^{j}, \label{eq:contra95metric}\tag{19}\] where \(\bar{\mathbf{g}}^i\) and \({\mathbf{g}}^i\) denotes the contravariant basis vectors in reference and deformed configuration of the shell-space defined by \[\bar{\mathbf{g}}^i \cdot \bar{\mathbf{g}}_j = \delta^i_j \quad \text{and} \quad {\mathbf{g}}^i \cdot {\mathbf{g}}_j = \delta^i_j. \label{eq:deform95tensor}\tag{20}\]

3.2 Kinematics↩︎

The deformation gradient tensor \(\mathbf{F}\) is defined by \[\mathbf{F} = \frac{\partial \mathbf{r}}{\partial \bar{\mathbf{r}}}= \frac{\partial\mathbf{r}}{\partial\theta^{i}}\otimes\frac{\partial\theta^{i}}{\partial\bar{\mathbf{r}}}=\mathbf{g}_{i}\otimes\mathbf{\bar{g}}^{i},\] thus the right Cauchy-Green deformation tensor is computed as \[\mathbf{C}=\mathbf{F}^\mathrm{T}\mathbf{F}=g_{ij}\mathbf{\bar{g}}^i\otimes\mathbf{\bar{g}}^j.\] For thin shells undergoing moderate deformations, the out-of-plane shear terms (\(C_{\alpha 3}\)) are negligible. This simplifies the tensor to a membrane-dominated form: \[\mathbf{C} = g_{\alpha\beta} \bar{\mathbf{g}}^\alpha \otimes \bar{\mathbf{g}}^\beta + [\lambda_3]^2 \bar{\mathbf{g}}^3 \otimes \bar{\mathbf{g}}^3. \label{eq:right95cauchy95green95simplified}\tag{21}\] Consequently, the inverse right Cauchy-Green tensor adopts a simplified block-diagonal form: \[[\mathbf{C}^{-1}] = \begin{bmatrix} C^{-1}_{11} & C^{-1}_{12} & 0 \\ C^{-1}_{12} & C^{-1}_{22} & 0 \\ 0 & 0 & C^{-1}_{33} \end{bmatrix}, \label{eq:block95diagonal}\tag{22}\] where \(C^{-1}_{33} = [\lambda_3]^{-2}\). The Green-Lagrange strain tensor \(\boldsymbol{\mathcal{E}}\) can be expressed as \[\boldsymbol{\mathcal{E}}=\frac{1}{2}[\mathbf{C}-\mathbf{I}]=\underbrace{\frac{1}{2}[g_{\alpha\beta}-\bar{g}_{\alpha\beta}]\bar{\mathbf{g}}^\alpha\otimes\bar{\mathbf{g}}^\beta}_{\tilde{\boldsymbol{\mathcal{E}}}}+\frac{1}{2}\left[[\lambda_3]^2-1\right]\bar{\mathbf{g}}^3\otimes\bar{\mathbf{g}}^3, \label{eq:green95lagrange95strain}\tag{23}\] where \(\tilde{\boldsymbol{\mathcal{E}}}\) is the in-plane strain tensor. It is further decomposed into two parts: \[\tilde{\boldsymbol{\mathcal{E}}} = \boldsymbol{\epsilon} + \theta^3 \boldsymbol{\kappa}, \label{eq:strain95decomp}\tag{24}\] where \(\boldsymbol{\epsilon}\) represents the membrane strain tensor, capturing in-plane stretching and shearing deformations, while \(\boldsymbol{\kappa}\) represents the bending strain tensor, describing curvature changes due to bending or twisting. Their components are computed as \[\epsilon_{\alpha\beta} = \frac{1}{2} [a_{\alpha \beta} - \bar{a}_{\alpha \beta}]\quad \text{and} \quad \kappa_{\alpha\beta} =[- b_{\alpha \beta} + \bar {b}_{\alpha \beta}], \label{eq:GL95strian952}\tag{25}\] with \[\bar{a}_{\alpha\beta} = \bar{\mathbf{a}}_\alpha \cdot \bar{\mathbf{a}}_\beta,\, {a}_{\alpha\beta} = {\mathbf{a}}_\alpha \cdot {\mathbf{a}}_\beta\quad \text{and} \quad \bar{b}_{\alpha\beta} = \bar{\mathbf{a}}_{\alpha,\beta} \cdot \bar{\mathbf{a}}_3,\,{b}_{\alpha\beta} = {\mathbf{a}}_{\alpha,\beta} \cdot {\mathbf{a}}_3.\]

3.3 Energetic Formulation and Weak Form↩︎

The governing equations for the electromechanical equilibrium of the shell are derived from the principle of stationary potential energy. The total potential energy \(\Pi_{\text{tot}}\) of the system consists of the internal energy \(\Pi_{\text{int}}\) and the external work \(\Pi_{\text{ext}}\):

\[\Pi_{\text{tot}} = \Pi_{\text{int}} + \Pi_{\text{ext}}.\]

3.3.1 Internal Energy Functional↩︎

The internal energy accounts for the stored mechanical and electrical energy, and enforces the incompressibility constraint via a Lagrange multiplier \(\tilde{p}_0\): \[\Pi_{\text{int}} = \int_{\bar{\Omega}} \widetilde{W}_{\text{mech}} \, dV + \int_{\bar{\Omega}} \widetilde{W}_{\text{elec}} \, dV - \int_{\bar{\Omega}} \tilde{p}_0 [\mathcal{J} - 1] \, dV, \label{eq:internal95energy}\tag{26}\] This functional is the direct application of the strain energy density formulation from Section 2.5 to the principle of stationary potential energy. Equation 11 is adopted to separate the mechanical and electrical contributions. \(\widetilde{W}_{\mathrm{mech}}\) and \(\widetilde{W}_{\mathrm{elec}}\) are the mechanical and electrical energy densities further defined in Equations 30 and 32 , respectively, and \(\mathcal{J} = \det(\mathbf{F}) = 1\) is the incompressibility constraint.

3.3.2 First Variation and Weak Form↩︎

The equilibrium state corresponds to a stationary point of the total energy. Taking the first variation \(\delta\Pi_{\text{tot}} = 0\) yields the weak form of the balance laws: \[\delta\Pi_{\text{tot}} = \delta\Pi_{\text{int}} + \delta\Pi_{\text{ext}} = 0. \label{eq:weak95form}\tag{27}\] The internal virtual work is obtained as the variation of the internal energy: \[\delta\Pi_{\text{int}} = \int_{\bar{\Omega}} \delta W_{\text{int}} \, \mathrm{d}V = \int_{\bar{\Omega}} {\mathbf{S}} : \delta{\boldsymbol{\mathcal{E}}} \, \mathrm{d}V, \label{eq:variation95int95energy}\tag{28}\] where \({\mathbf{S}}\) is the total second Piola–Kirchhoff stress tensor and \(\delta{\boldsymbol{\mathcal{E}}}\) is the variation of the Green-Lagrange strain tensor. Remark: Consistent with the Kirchhoff–Love hypothesis and the plane stress assumption (\({S}^{33} = 0\)), the thickness strain \({\mathcal{E}}_{33}\) and transverse shear strains \({\mathcal{E}}_{\alpha 3}\) do not contribute to the virtual work. While the thickness stretch \(\lambda_3\) is kinematically determined by the in-plane deformation and incompressibility, the vanishing stress \({S}^{33}\) ensures that \({\mathcal{E}}_{33}\) performs no work. Consequently, the internal energy depends only on the in-plane strains, which are the sole contributors to the weak form. Thus, Equation 28 reduces to \[\delta\Pi_{\text{int}} = \int_{\bar{\Omega}} \tilde{\mathbf{S}} : \delta\tilde{\boldsymbol{\mathcal{E}}} \, \mathrm{d}V,\] where \(\tilde{\mathbf{S}}\) and \(\delta\tilde{\boldsymbol{\mathcal{E}}}\) are only in-plane tensors.

The external virtual work \(\delta\Pi_{\text{ext}}\) incorporates contributions from applied mechanical tractions and electrical boundary conditions, the specifics of which depend on the problem setup.

Equation 27 constitutes the nonlinear variational equation to be solved. In the present isogeometric discretisation (Section 4.1), it leads to a system of nonlinear algebraic equations for the control point displacements and electric potential.

3.3.3 Linearisation and Material Tangent↩︎

To solve Equation 27 using the Newton-Raphson method, consistent linearisation is required. The directional derivative (linearisation) of the weak form yields the tangent stiffness operator. This involves the linearisation of the stress, which introduces the material tangent moduli (detailed in Section 3.6).

The final discrete tangent stiffness matrix assembled from the finite element discretisation therefore comprises both geometric stiffness contributions (from the linearisation of the strain variation \(\delta\boldsymbol{\mathcal{E}}\)) and material stiffness contributions (from the tangent moduli \(\mathbb{C}^{ijkl}\)). The specific expressions for the condensed plane-stress tangent moduli \(\hat{\mathbb{C}}^{\alpha\beta\gamma\delta}\), which are used in the shell resultant formulation, will be provided in Section 3.6.3.

3.4 Total Stress in Electroelastic Thin Shells↩︎

The total stress within the dielectric elastomer shell originates from three contributions, expressed through the constitutive relationship: \[{\mathbf{S}} ={2\frac{\partial \widetilde{W}_{\mathrm{mech}}}{\partial \mathbf{C}}}+ {2\frac{\partial \widetilde{W}_{\mathrm{elec}}}{\partial \mathbf{C}}} - {\tilde{p}_0 \mathbf{C}^{-1}}. \label{eq:total95stress95decomposition}\tag{29}\] where the first term is the hyperelastic mechanical stress arising from the finite deformation of the dielectric elastomer. The second term corresponds to the Maxwell stress generated by the interaction between the applied electric field and dielectric material, and the third term accounts for the hydrostatic pressure that maintains the incompressibility constraint.

3.4.1 Mechanical Stress Contribution↩︎

The present framework is general and can be applied to any hyperelastic constitutive equation. Here, the Mooney–Rivlin model is employed as a representative example. Assuming a mechanical energy density \(\widetilde{W}_{\mathrm{mech}}\) for finite deformations, the Mooney–Rivlin form is: \[\widetilde{W}_{\mathrm{mech}}(\mathbf{C}) = c_1[I_1 - 3] + c_2[I_2 - 3], \label{eq:mech95energy}\tag{30}\] where the invariants \(I_1\) and \(I_2\) are defined in Equation 9 . The corresponding stress contribution derives from the derivative of \(\widetilde{W}_{\mathrm{mech}}\) with respect to \(\mathbf{C}\): \[{\mathbf{S}}_{\text{mech}} \equiv 2\frac{\partial \widetilde{W}_{\mathrm{mech}}}{\partial \mathbf{C}} = 2c_1\frac{\partial I_1}{\partial \mathbf{C}} + 2c_2\frac{\partial I_2}{\partial \mathbf{C}}. \label{eq:mech95stress95deriv}\tag{31}\]

3.4.2 Electrically Induced Stress↩︎

The electric energy density \(\widetilde{W}_{\mathrm{elec}}\) captures the energy stored in the dielectric material due to polarisation under an applied electric field. For an isotropic voltage‑controlled elastomers, we recall the expression 12 here: \[\widetilde{W}_{\mathrm{elec}}(\mathbf{C},\bar{\mathbf{E}}) = -\frac{1}{2}\epsilon [\bar{\mathbf{E}}\otimes \bar{\mathbf{E}}]: \mathbf{C}^{-1}. \label{eq:elec95energy}\tag{32}\]

The electrically induced stress contribution to the total stress, the Maxwell stress, is derived as the work conjugate to the material strain measure. Applying the chain rule, its derivative the Maxwell stress yields: \[{\mathbf{S}}_{\text{elec}} \equiv 2\frac{\partial \widetilde{W}_{\mathrm{elec}}}{\partial \mathbf{C}} = \epsilon \mathbf{C}^{-1}[\bar{\mathbf{E}} \otimes \bar{\mathbf{E}}]\mathbf{C}^{-1}. \label{eq:elec95stress95deriv}\tag{33}\] This expression represents the Maxwell stress in material coordinates, which arises from electrostatic interactions within the dielectric medium. For thin shell applications, a key simplification occurs when an electric potential difference \(\Delta\Phi\) is applied across the thickness. The spatial electric field then simplifies to: \[\mathbf{E} = -\frac{\Delta\Phi}{h}\mathbf{a}_3 = -\frac{\Delta\Phi}{\lambda_3 \bar h}\mathbf{a}_3,\] The corresponding material electric field, computed via the inverse deformation gradient, becomes: \[\bar{\mathbf{E}} = \mathbf{F}^{-\mathrm{T}}\mathbf{E} = -\frac{\Delta\Phi}{\bar{h}}\bar{\mathbf{a}}_3 . \label{eq:electric95field95shell}\tag{34}\] This formulation confirms that the electric field remains aligned with the out-of-plane direction throughout deformation. Substituting the material electric field expression into Equation 33 yields the simplified Maxwell stress tensor: \[{\mathbf{S}}_{\text{elec}} = \epsilon\left[\frac{\Delta\Phi}{\bar{h}}\right]^2 \mathbf{C}^{-1}[\bar{\mathbf{a}}_3 \otimes \bar{\mathbf{a}}_3]\mathbf{C}^{-1}. \label{eq:se95tensor}\tag{35}\] This expression is valid for the general anisotropic case and captures the full electromechanical coupling.

Applying the transverse isotropy simplification to the Maxwell stress expression, the quadratic product \(\mathbf{C}^{-1}[\bar{\mathbf{a}}_3 \otimes \bar{\mathbf{a}}_3]\mathbf{C}^{-1}\) reduces to: \[\mathbf{C}^{-1}[\bar{\mathbf{a}}_3 \otimes \bar{\mathbf{a}}_3]\mathbf{C}^{-1} = [\lambda_3]^{-4} [\bar{\mathbf{a}}_3 \otimes \bar{\mathbf{a}}_3]. \label{eq:out95of95plane95term}\tag{36}\] This dramatic simplification holds because \(\bar{\mathbf{a}}_3\) is an eigenvector of \(\mathbf{C}^{-1}\), a property guaranteed by the assumption of transverse isotropy, satisfying \(\mathbf{C}^{-1}\bar{\mathbf{a}}_3 = [\lambda_3]^{-2}\bar{\mathbf{a}}_3\). Exploiting the orthogonality of the shell director \(\bar{\mathbf{a}}_3\) to the mid-surface, the primary stress component normal to the mid-surface is obtained by taking the tensor contraction: \[S_{\text{elec}}^{33} = [\bar{\mathbf{a}}_3 \otimes \bar{\mathbf{a}}_3] : \mathbf{S}_{\text{elec}} = \epsilon [\lambda_3]^{-4}\left[\frac{\Delta\Phi}{\bar{h}}\right]^2. \label{eq:s3395result}\tag{37}\] This is the principal electromechanical stress component driving thickness changes in actuation. The transverse isotropy assumption and normal electric field orientation cause all off-diagonal and in-plane components of \(\mathbf{S}_{\text{elec}}\) to vanish. Mathematically, this occurs because: \[{S}_{\text{elec}}^{\alpha\beta} = [\bar{\mathbf{a}}_\alpha \otimes \bar{\mathbf{a}}_\beta] : \mathbf{S}_{\text{elec}} = 0 \quad \text{for} \quad \alpha,\beta = 1,2,\] due to orthogonality between \(\bar{\mathbf{a}}_3\) and \(\bar{\mathbf{a}}_\alpha\). Physically, this suppression of in-plane electromechanical coupling arises because the applied electric field is oriented exclusively normal to the mid-surface.

3.4.3 Explicit Enforcement of Plane Stress and Incompressibility↩︎

For thin-shell structures, three physical considerations justify the plane stress assumption:

1. Dimensional disparity: The thickness dimension is orders of magnitude smaller than in-plane dimensions

2. Boundary conditions: Both top and bottom surfaces are traction-free (\(\mathbf{t} = \mathbf{0}\))

3. Stress magnitude: Through-thickness stresses are negligible compared to in-plane stresses

Remark: The Maxwell stress is an internal electromechanical coupling effect and does not constitute an external mechanical traction on the boundaries. The electrodes impose only the electric potential, leaving the mechanical traction \(\mathbf{t}=\mathbf{0}\). The thickness component of the Maxwell stress is internally balanced by the Lagrange multiplier \(\tilde{p}_0\) to satisfy \({S}^{33}=0\).

These conditions remain valid for electroelastic shells under actuation, leading to the equilibrium condition: \[{S}^{33} = S^{33}_{\text{elec}} + S^{33}_{\text{mech}} - \tilde{p}_0 C^{33} = 0, \label{eq:plane95stress95condition}\tag{38}\] where \(C^{33} = [\lambda_3]^{-2}\). This equation expresses the vanishing normal stress in the thickness direction. The Lagrange multiplier \(\tilde{p}_0\) is explicitly determined by solving Equation 38 . The resulting expression is: \[\begin{align} \tilde{p}_0 &= [\lambda_3]^2\left[S^{33}_{\text{mech}} + S^{33}_{\text{elec}}\right] \nonumber \\ &= 2[\lambda_3]^2\frac{\partial \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}} + [\lambda_3]^{-2} \epsilon \left[\frac{\Delta\Phi}{\bar{h}}\right]^2. \label{eq:pressure95expression} \end{align}\tag{39}\] This solution strategy eliminates \(\tilde{p}_0\) as an additional unknown variable while ensuring the exact enforcement of both constraints. It clearly separates into mechanical and electrical contributions here. The electrical term \([\lambda_3]^{-2} \epsilon [\Delta\Phi/\bar{h}]^2\) represents ‘electrostatic pressure’, which decreases with thickness stretch.

Substituting \(\tilde{p}_0\) into the general stress expression yields the working in-plane stress components: \[\begin{align} \tilde{S}^{\alpha\beta} &= \underbrace{2\left[\frac{\partial \widetilde{W}_{\mathrm{mech}}}{\partial C_{\alpha\beta}} - [\lambda_3]^2\frac{\partial \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}}C^{\alpha\beta}\right]}_{\text{Mechanical stress}} - \underbrace{[\lambda_3]^{-2} \epsilon \left[\frac{\Delta\Phi}{\bar{h}}\right]^2C^{\alpha\beta}}_{\text{Electrically induced stress}} \label{eq:inplane95stress} \end{align}\tag{40}\] This formulation reveals that the electrical term scales with \((\Delta\Phi)^2\) and \(C^{\alpha\beta}\), creating voltage-dependent in-plane stresses and they increase with thickness reduction. The decoupled structure enables efficient implementation while capturing the essential electromechanical coupling mechanisms.

3.4.4 Stress Resultants for Electroelastic Thin Shells↩︎

The thin shell formulation reduces the three-dimensional continuum to a two-dimensional surface with a thickness. The internal forces are expressed as the resultants of the integrated stress through the thickness, defined by the in-plane second Piola–Kirchhoff stress \(\tilde{\mathbf{S}}\). The resultant of the in-plane stress \(\boldsymbol{\mathcal{N}}\) and the bending moment \(\boldsymbol{\mathcal{M}}\) have their components calculated as: \[\begin{align} n^{\alpha\beta} &= \int_{-\frac{\bar{h}}{2}}^{\frac{\bar{h}}{2}} \tilde{S}^{\alpha\beta} J_c \, \mathrm{d}\theta^3, \\ m^{\alpha\beta} &= \int_{-\frac{\bar{h}}{2}}^{\frac{\bar{h}}{2}} \tilde{S}^{\alpha\beta} \theta^3 J_c \, \mathrm{d}\theta^3, \end{align} \label{eq:stress95resultants}\tag{41}\] where \(\theta^3 \in [-\bar{h}/2, \bar{h}/2]\) is the coordinate along the reference thickness and \(J_c\) is the thickness-direction Jacobian correction factor: \[J_c = \frac{\left| [\bar{\mathbf{g}}_1 \times \bar{\mathbf{g}}_2] \cdot \bar{\mathbf{g}}_3 \right|}{\left| [\bar{\mathbf{a}}_1 \times \bar{\mathbf{a}}_2] \cdot \bar{\mathbf{a}}_3 \right|}, \label{eq:jacobian95correction}\tag{42}\] accounting for volume changes between the reference base vectors \(\bar{\mathbf{g}}_i\) and mid-surface basis \(\bar{\mathbf{a}}_\alpha\).

The stress increments are linearised via the constitutive relation between \({\mathbf{S}}\) and the Green-Lagrange strain \(\boldsymbol{\mathcal{E}}\) (Equation 29 ) as \[\mathrm{d} {S}^{ij}=\frac{\partial S^{ij}}{\partial \mathcal{E}_{kl}}\mathrm{d}\mathcal{E}_{kl}=2\frac{\partial S^{ij}}{\partial C_{kl}}\mathrm{d}\mathcal{E}_{kl}=\hat{\mathbb{C}}^{ijkl}\mathrm{d}\mathcal{E}_{kl} \label{eq:stress95increment}\tag{43}\] where \(\hat{\mathbb{C}}^{ijkl}\) are the components of the fourth-order elasticity tensor.

3.5 Variation of Total Energy↩︎

The total potential energy \(\Pi_{\text{tot}}\) of the electroelastic thin shell system comprises internal energy from finite deformation and polarisation (\(\Pi_{\text{int}}\)) and external work contributions (\(\Pi_{\text{ext}}\)). Its first variation is \[\delta\Pi_{\text{tot}} = \delta\Pi_{\text{int}} + \delta\Pi_{\text{ext}} = \int_{\bar{\Omega}} \delta W_{\text{int}} \mathrm{~d} V + \int_{\Omega}\delta W_{\text{ext}} \mathrm{~d} V.\] The tangent stiffness required for Newton-Raphson iterations derives from consistent linearisation of this variation, are detailed in Section 4.3.

3.6 Material Tangent Moduli and Static Condensation↩︎

3.6.1 Pressure Derivatives for Constraint Enforcement↩︎

The Lagrange multiplier \(\tilde{p}_0\), which enforces incompressibility (\(\mathcal{J} = 1\)), depends implicitly on deformation through the plane stress condition 39 . Its derivatives with respect to the deformation components are essential for consistent tangent moduli: \[\begin{align} \frac{\partial \tilde{p}_0}{\partial C_{\alpha\beta}} &= 2[\lambda_3]^2\frac{\partial^2 \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}\partial C_{\alpha\beta}}, \\ \frac{\partial \tilde{p}_0}{\partial C_{33}} &= 2[\lambda_3]^2\frac{\partial^2 \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}^2} + 2\frac{\partial \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}} - [\lambda_3]^{-4} \epsilon \left[\frac{\Delta\Phi}{\bar{h}}\right]^2. \end{align}\] These capture how pressure responds to: (1) in-plane stretching (\(\alpha,\beta=1,2\)) through mechanical-kinematic coupling, and (2) thickness changes (\(C_{33}\)) with explicit electromechanical contributions.

3.6.2 Material Tangent Moduli Derivation↩︎

The fourth-order elasticity tensor \(\mathbb{C}^{ijkl} \equiv 2\frac{\partial {S}^{ij}}{\partial C_{kl}}\) is derived by consistent differentiation of the total stress (Equation 29 ) with respect to the right Cauchy-Green tensor. Its partitions exhibit distinct symmetries:

In-plane Moduli (\(\alpha\beta\gamma\delta\))

\[\begin{align} \mathbb{C}^{\alpha\beta\gamma\delta} &= 4\frac{\partial^2 \widetilde{W}_{\mathrm{mech}}}{\partial C_{\alpha\beta} \partial C_{\gamma\delta}} - 4[\lambda_3]^2\left[\frac{\partial^2 \widetilde{W}_{\mathrm{mech}}}{\partial C_{33} \partial C_{\gamma\delta}}C^{\alpha\beta} + \frac{\partial^2 \widetilde{W}_{\mathrm{mech}}}{\partial C_{33} \partial C_{\alpha\beta}}C^{\gamma\delta}\right] \nonumber \\ &\quad - {\left[2[\lambda_3]^2\frac{\partial \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}} + [\lambda_3]^{-2} \epsilon \left[\frac{\Delta\Phi}{\bar{h}}\right]^2\right]}\left[C^{\alpha\beta}C^{\gamma\delta} - C^{\alpha\gamma}C^{\beta\delta} - C^{\alpha\delta}C^{\beta\gamma}\right], \label{eq:C95alpha95beta95gamma95delta} \end{align}\tag{44}\] where the last term arises from the product rule applied to \(-\tilde{p}_0\mathbf{C}^{-1}\). This expression has the major symmetry \(\mathbb{C}^{\alpha\beta\gamma\delta} = \mathbb{C}^{\gamma\delta\alpha\beta}\).

Thickness Coupling Moduli (\(\alpha\beta33\))

\[\begin{align} \mathbb{C}^{\alpha\beta33} &= -C^{\alpha\beta}\left[6\frac{\partial \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}} + 4[\lambda_3]^2\frac{\partial^2 \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}^2} - [\lambda_3]^{-4} \epsilon \left[\frac{\Delta\Phi}{\bar{h}}\right]^2\right], \label{eq:C95alpha95beta9533} \end{align}\tag{45}\] quantifying how in-plane stresses change with thickness stretch.

Thickness Moduli (\(3333\))

\[\begin{align} \mathbb{C}^{3333} &= -[\lambda_3]^{-2}\left[6\frac{\partial \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}} + 4[\lambda_3]^2\frac{\partial^2 \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}^2} - [\lambda_3]^{-4} \epsilon \left[\frac{\Delta\Phi}{\bar{h}}\right]^2\right], \label{eq:C953333} \end{align}\tag{46}\] governing thickness-direction stiffness. Its negative definiteness reflects the kinematic constraint.

3.6.3 Static Condensation for Plane Stress↩︎

To enforce the plane stress condition (\(\tilde{S}^{33}=0\)), the moduli undergo static condensation: \[\begin{align} \hat{\mathbb{C}}^{\alpha\beta\gamma\delta} &= \mathbb{C}^{\alpha\beta\gamma\delta} - \frac{\mathbb{C}^{\alpha\beta33}\mathbb{C}^{33\gamma\delta}}{\mathbb{C}^{3333}}. \label{eq:condensation} \end{align}\tag{47}\] This simplifies the 3D constitutive relation to 2D, eliminating explicit \(\tilde{p}_0\) dependence. Substituting Eqs. 44 46 yields: \[\begin{align} \hat{\mathbb{C}}^{\alpha\beta\gamma\delta} &= 4\frac{\partial^2 \widetilde{W}_{\mathrm{mech}}}{\partial C_{\alpha\beta} \partial C_{\gamma\delta}} - 4[\lambda_3]^2\left[\frac{\partial^2 \widetilde{W}_{\mathrm{mech}}}{\partial C_{33} \partial C_{\gamma\delta}}C^{\alpha\beta} + \frac{\partial^2 \widetilde{W}_{\mathrm{mech}}}{\partial C_{33} \partial C_{\alpha\beta}}C^{\gamma\delta}\right] \nonumber \\ &\quad + C^{\alpha\beta}C^{\gamma\delta}\left[6[\lambda_3]^2\frac{\partial \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}} + 4[\lambda_3]^4\frac{\partial^2 \widetilde{W}_{\mathrm{mech}}}{\partial C_{33}^2} - [\lambda_3]^{-2} \epsilon \left[\frac{\Delta\Phi}{\bar{h}}\right]^2\right], \label{eq:reduced95tangent} \end{align}\tag{48}\] where \(\hat{\mathbb{C}}^{\alpha\beta\gamma\delta}\) is the plane-stress-reduced tangent modulus, which incorporates both mechanical and electrical effects while satisfying plane stress constraints intrinsically.

3.7 Prestretch in Electroelastic Thin Shells↩︎

In manufacturing and application, dielectric elastomers are often subjected to prestretch to enhance actuation performance or achieve specific configurations. Figure 3 illustrates three states for prestretched electroelastic shells. The prestretched state, denoted as \(\bar\Omega\), serves as the reference configuration for subsequent electromechanical analysis. The deformation from the stress-free natural state \(\Omega_{\text{nat}}\) to the prestretched reference configuration \(\bar\Omega\) is described by the prestretch gradient \(\mathbf{F}_\mathrm{p}\). The total deformation gradient from \(\Omega_{\text{nat}}\) to the current configuration \(\Omega\) decomposes multiplicatively as: \[\mathbf{F}_{\text{total}} = \mathbf{F} \cdot \mathbf{F}_\mathrm{p}, \label{eq:F95total95decomp}\tag{49}\] where \(\mathbf{F}\) is the deformation gradient from the prestretched reference configuration to the current configuration, consistent with Section 3.2.

The total right Cauchy-Green deformation tensor relative to the natural state is:

\[\mathbf{C}_{\text{total}} = \mathbf{F}_\mathrm{p}^\top \mathbf{C} \mathbf{F}_\mathrm{p}. \label{eq:C95total95relation}\tag{50}\]

The mechanical strain energy density \(\widetilde{W}_{\mathrm{mech}}\) remains defined with respect to the natural state, making it a function of \(\mathbf{C}_{\text{total}}\).

3.7.1 Isotropic Prestretch↩︎

A common prestretch pattern is isotropic in-plane stretching accompanied by thickness reduction to satisfy incompressibility. This is represented by \[\mathbf{F}_\mathrm{p} = \operatorname{diag}(\lambda_p, \lambda_p, \lambda_p^{-2}), \label{eq:Fp95isotropic}\tag{51}\] where \(\lambda_p\) is the in-plane prestretch ratio. The thickness component \(\lambda_p^{-2}\) ensures \(\det(\mathbf{F}_\mathrm{p}) = 1\).

Under this assumption, the components of \(\mathbf{C}_{\text{total}}\) relate to those of \(\mathbf{C}\) (defined in the prestretched configuration) as:

\[\mathbf{C}_{\text{total}} = \begin{bmatrix} \lambda_p^2 C_{11} & \lambda_p^2 C_{12} & 0 \\ \lambda_p^2 C_{12} & \lambda_p^2 C_{22} & 0 \\ 0 & 0 & \lambda_p^{-4} C_{33} \end{bmatrix}, \label{eq:C95total95isotropic}\tag{52}\] where \(C_{33} = [\lambda_3]^2\) and \(\lambda_3\) is the through-thickness stretch from the prestretched reference state.

3.7.2 Stress Resultants with Prestretch↩︎

The in-plane stress components in the prestretched reference configuration follow the same general form as Equation 40 , but with the derivatives of \(\widetilde{W}_{\mathrm{mech}}\) now taken with respect to \(\mathbf{C}_{\text{total}}\). For the Mooney–Rivlin model, the derivatives become:

\[\begin{align} \frac{\partial\widetilde{W}_{\mathrm{mech}}}{\partial C_{\alpha\beta}} &= c_1\lambda_p^2\bar{g}^{\alpha\beta} + c_2\Bigl[\lambda_p^4\bigl[C_{\gamma\delta}\bar{g}^{\gamma\delta}\bar{g}^{\alpha\beta} - \bar{g}^{\alpha\gamma}C_{\gamma\delta}\bar{g}^{\delta\beta}\bigr] + \lambda_p^{-2}\bar{g}^{\alpha\beta}C_{33}\Bigr], \tag{53} \\ \frac{\partial\widetilde{W}_{\mathrm{mech}}}{\partial C_{33}} &= c_1\lambda_p^{-4} + c_2\lambda_p^{-2}C_{\gamma\delta}\bar{g}^{\gamma\delta}. \tag{54} \end{align}\]

Substituting these into Equation 40 yields the in-plane stress components explicitly accounting for prestretch:

\[\begin{align} \tilde{S}^{\alpha\beta} = 2\Bigg[ &c_1\lambda_p^2\bar{g}^{\alpha\beta} + c_2\Bigl[\lambda_p^4\bigl[C_{\gamma\delta}\bar{g}^{\gamma\delta}\bar{g}^{\alpha\beta} - \bar{g}^{\alpha\gamma}C_{\gamma\delta}\bar{g}^{\delta\beta}\bigr] + \lambda_p^{-2}\bar{g}^{\alpha\beta}C_{33}\Bigr] \nonumber \\ &- [\lambda_3]^2 \bigl[c_1\lambda_p^{-4} + c_2\lambda_p^{-2}C_{\gamma\delta}\bar{g}^{\gamma\delta}\bigr] C^{\alpha\beta} \Bigg] - [\lambda_3]^{-2} \epsilon \left[ \frac{\Delta\Phi}{\bar{h}} \right]^2 C^{\alpha\beta}. \label{eq:inplane95stress95with95prestretch} \end{align}\tag{55}\]

3.7.3 Tangent Moduli with Prestretch↩︎

The consistent tangent moduli for the prestressed configuration are derived in the same manner as in Equation 48 , but with the derivatives of \(\widetilde{W}_{\mathrm{mech}}\) now evaluated with respect to \(\mathbf{C}_{\text{total}}\) and then transformed to the prestressed reference configuration. The reduced in-plane tangent moduli \(\hat{\mathbb{C}}^{\alpha\beta\gamma\delta}\) are given by the same expression as Equation 48 , but with the following derivatives for the Mooney–Rivlin model under isotropic prestretch: \[\begin{align} \frac{\partial^2\widetilde{W}_{\mathrm{mech}}}{\partial C_{\alpha\beta}\partial C_{\gamma\delta}} &= c_2\lambda_p^4\left[\bar{g}^{\alpha\beta}\bar{g}^{\gamma\delta} \frac{1}{2}\left[\bar{g}^{\alpha\gamma}\bar{g}^{\beta\delta} + \bar{g}^{\alpha\delta}\bar{g}^{\beta\gamma}\right]\right], \tag{56} \\ \frac{\partial^2\widetilde{W}_{\mathrm{mech}}}{\partial C_{\alpha\beta}\partial C_{33}} &= c_2\lambda_p^{-2}\bar{g}^{\alpha\beta}, \tag{57} \\ \frac{\partial^2\widetilde{W}_{\mathrm{mech}}}{\partial C_{33}^2} &= 0. \tag{58} \end{align}\]

These expressions, along with the first derivatives given in Equations 53 and 54 , are used to compute the tangent moduli in the prestressed configuration. The electrical contribution remains unchanged, as it is independent of the prestretch.

3.7.4 Remarks on Implementation↩︎

Incorporating prestretch within the isogeometric shell framework requires the following modifications:

1. The reference geometry is defined in the prestretched configuration \(\bar\Omega\).

2. The prestretch tensor \(\mathbf{F}_\mathrm{p}\) is stored as a field (constant or varying spatially).

3. The strain energy derivatives in the weak form and tangent stiffness are evaluated using Equations 53 –58 .

4. The thickness stretch \(\lambda_3\) is determined from the incompressibility constraint \(\mathcal{J}=1\), which now includes the prestretch effect.

This formulation enables the analysis of prestrained dielectric elastomer shells undergoing large electromechanical deformations while maintaining the \(C^1\)-continuity requirements of the Kirchhoff–Love shell theory.

4 Numerical Implementation↩︎

The nonlinear finite element implementation leverages the deal.II library [45] to realize the electroelastic thin shell formulation. The primary numerical challenge lies in the coupling between the large-deformation mechanics of a Kirchhoff–Love shell and the electrostatic forces arising from an applied voltage. This section outlines the key components of the numerical discretisation, the linearisation of the weak form, and the resulting solution procedure.

4.1 Subdivision Surface Discretisation↩︎

As established in our previous work on hyperelastic thin shells [39], the mid-surface geometry and the displacement field are discretised using Catmull–Clark subdivision surfaces. This choice is motivated by the \(C^1\)-continuity requirement of the Kirchhoff–Love theory, which necessitates basis functions for the displacement in the Sobolev space \(H^2(\Omega)\). Subdivision surfaces provide smooth, \(C^1\)-continuous limit surfaces everywhere, even on unstructured control meshes containing extraordinary vertices. The mid-surface in the reference configuration \(\bar{\mathbf{x}}\) and the displacement field \(\mathbf{u}\) are approximated using the same set of subdivision basis functions \[\begin{align} \bar{\mathbf{x}}(\theta^1, \theta^2) &\approx \sum_{A=1}^{n_{\text{node}}} N^A(\theta^1, \theta^2) \bar{\mathbf{X}}^A, \\ \mathbf{u}(\theta^1, \theta^2) &\approx \sum_{A=1}^{n_{\text{node}}} N^A(\theta^1, \theta^2) \mathbf{u}^A, \label{eq:displacement95interpolation} \end{align}\tag{59}\] where \(\bar{\mathbf{X}}^A\) and \(\mathbf{u}^A\) are the reference position and displacement vector of control point \(A\), respectively, and \(n_{\text{node}}\) is the number of control points in the support of a given parametric location. The deformed mid-surface position is then \({\mathbf{x}} = \bar{\mathbf{x}} + {\mathbf{u}}\). The functions \(N^A(\theta^1, \theta^2)\) are the cubic Catmull-Clark subdivision bases. For regular patches (interior elements with valence 4), they are equivalent to bi-cubic B-splines. For patches containing extraordinary vertices, Stam’s algorithm [46] is employed for fast and exact evaluation.

4.2 Discretised Kinematic Fields↩︎

All strain and curvature measures are derived from the discretised mid-surface geometry. The covariant basis vectors on the deformed mid-surface are computed as the partial derivatives of the interpolated position: \[\mathbf{a}_\alpha = \mathbf{x}_{,\alpha} = \bar{\mathbf{a}}_\alpha + \sum_{A} N^A_{,\alpha} \mathbf{u}^A.\] The metric components \(\epsilon_{\alpha\beta}\) and curvature components \(\kappa_{\alpha\beta}\) (and their referential counterparts) are then calculated according to their definitions in Section 2.2. The Green–Lagrange membrane and bending strain components, as per Equation 25 , become functions of the nodal displacements: \[\begin{align} \epsilon_{\alpha\beta} &= \frac{1}{2}[\mathbf{a}_\alpha \cdot \mathbf{a}_\beta - \bar{\mathbf{a}}_\alpha \cdot \bar{\mathbf{a}}_\beta], \\ \kappa_{\alpha\beta} &= -[\mathbf{a}_{\alpha,\beta} \cdot \mathbf{a}_3]+ [\bar{\mathbf{a}}_{\alpha,\beta} \cdot \bar{\mathbf{a}}_3]. \end{align}\] The thickness stretch \(\lambda_3\), required for the stress evaluation, is computed from the area change of the in-plane basis vectors: \[\lambda_3 = \frac{|\bar{\mathbf{a}}_1 \times \bar{\mathbf{a}}_2|}{|{\mathbf{a}}_1 \times {\mathbf{a}}_2|} = \frac{\bar{J}}{J}.\] The variations of the strains, \(\delta\epsilon_{\alpha\beta}\) and \(\delta\kappa_{\alpha\beta}\), are obtained by taking the directional derivative (or Gateaux derivative) of the above expressions with respect to the displacement field \(\mathbf{u}\), leading to expressions linear in \(\delta \mathbf{u}\).

4.3 Weak Form and Consistent Linearisation↩︎

The principle of virtual work for the shell, including the internal electromechanical stresses and external loads, is \[\delta \Pi_{\text{tot}} = \delta\Pi_{\text{int}} + \delta\Pi_{\text{ext}} = 0,\] where the internal virtual work is the integral of the stress resultants (Equation 41 ) working through the virtual strains: \[\delta\Pi_{\text{int}} = \int_{\bar{\Gamma}} \left[ n^{\alpha\beta} \, \delta\epsilon_{\alpha\beta} + m^{\alpha\beta} \, \delta\kappa_{\alpha\beta} \right] d\bar{\Gamma}. \label{eq:internal95virtual95work}\tag{60}\] Here, \(\bar{\Gamma}\) is the shell mid-surface in the reference configuration, and \(d\bar{\Gamma} = \bar{J} \, d\theta^1 d\theta^2\). The stress resultants \(n^{\alpha\beta}\) and \(m^{\alpha\beta}\) encapsulate the full electromechanical coupling. They are computed by integrating the Piola–Kirchhoff stress \(\tilde{S}^{\alpha\beta}\) (Equation 40 ) through the thickness, as defined in Equation 41 . The integration is performed numerically using a Gauss quadrature rule along \(\theta^3\). The Jacobian correction factor \(J_c\) in Equation 42 accounts for the curvature of the shell space. The external virtual work \(\delta\Pi_{\text{ext}}\) accounts for mechanical surface tractions, pressure loads, and follower forces, which are standard in shell formulations.

The solution is obtained via the Newton–Raphson method, which requires the linearisation of Equation 60 . The linearised increment of the internal virtual work is: \[\varDelta(\delta\Pi_{\text{int}}) = \int_{\bar{\Gamma}} \left[ \varDelta n^{\alpha\beta} \, \delta\epsilon_{\alpha\beta} + n^{\alpha\beta} \, \varDelta(\delta\epsilon_{\alpha\beta}) + \varDelta m^{\alpha\beta} \, \delta\kappa_{\alpha\beta} + m^{\alpha\beta} \, \varDelta(\delta\kappa_{\alpha\beta}) \right] d\bar{\Gamma}.\] The stress resultant increments \(\varDelta n^{\alpha\beta}\) and \(\varDelta m^{\alpha\beta}\) are related to the increments of the strain measures via the constitutive tangent moduli: \[\begin{align} \varDelta n^{\alpha\beta} &= \mathbb{A}^{\alpha\beta\gamma\delta} \, \varDelta\epsilon_{\gamma\delta} + \mathbb{B}^{\alpha\beta\gamma\delta} \, \varDelta\kappa_{\gamma\delta}, \\ \varDelta m^{\alpha\beta} &= \mathbb{B}^{\gamma\delta\alpha\beta} \, \varDelta\epsilon_{\gamma\delta} + \mathbb{D}^{\alpha\beta\gamma\delta} \, \varDelta\kappa_{\gamma\delta}. \end{align}\] The tangent stiffness tensors \(\mathbb{A}\), \(\mathbb{B}\), and \(\mathbb{D}\) are the integrated counterparts of the reduced material tangent modulus \(\hat{\mathbb{C}}^{\alpha\beta\gamma\delta}\) (Equation 48 ): \[\begin{align} \mathbb{A}^{\alpha\beta\gamma\delta} &= \int_{-\frac{\bar{h}}{2}}^{\frac{\bar{h}}{2}} \hat{\mathbb{C}}^{\alpha\beta\gamma\delta} \, J_c \, d\theta^3, \\ \mathbb{B}^{\alpha\beta\gamma\delta} &= \int_{-\frac{\bar{h}}{2}}^{\frac{\bar{h}}{2}} \hat{\mathbb{C}}^{\alpha\beta\gamma\delta} \, \theta^3 \, J_c \, d\theta^3, \\ \mathbb{D}^{\alpha\beta\gamma\delta} &= \int_{-\frac{\bar{h}}{2}}^{\frac{\bar{h}}{2}} \hat{\mathbb{C}}^{\alpha\beta\gamma\delta} \, [\theta^3]^2 \, J_c \, d\theta^3. \label{eq:tangent95tensors} \end{align}\tag{61}\] Note that \(\hat{\mathbb{C}}^{\alpha\beta\gamma\delta}\) depends on the current deformation (via \(\mathbf{C}\) and \(\lambda_3\)) and the applied voltage \(\Delta\Phi\). Its explicit form, provided in Equation 48 , is evaluated at each integration point during the assembly of the stiffness matrix. The terms \(\varDelta(\delta\epsilon_{\alpha\beta})\) and \(\varDelta(\delta\kappa_{\alpha\beta})\) arise from the geometric non-linearity and contribute to the geometric stiffness matrix. Their detailed derivation involves the second variation of the kinematic quantities, which is standard in non-linear shell formulations (see, e.g., [47] for the Kirchhoff-Love case and [48] for the general geometrically nonlinear shell framework).

4.4 Finite Element Formulation↩︎

Substituting the discretised kinematic fields into the weak form yields the discrete nonlinear equilibrium equations. For a control point \(A\), the residual force is defined by \[\boldsymbol{\mathsf{R}}^{A} = \boldsymbol{\mathsf{F}}^{A,\mathrm{int}} - \boldsymbol{\mathsf{F}}^{A,\mathrm{ext}},\] where the internal force vector is given by \[\boldsymbol{\mathsf{F}}^{A,\mathrm{int}} = \int_{\bar{\Gamma}} \left[ n^{\alpha\beta} \frac{\partial \epsilon_{\alpha\beta}}{\partial \boldsymbol{\mathsf{u}}^{A}} + m^{\alpha\beta} \frac{\partial \kappa_{\alpha\beta}}{\partial \boldsymbol{\mathsf{u}}^{A}} \right] d\bar{\Gamma}. \label{eq:internal95force}\tag{62}\]

The tangent stiffness matrix is obtained by consistent linearisation of the residual vector with respect to the displacement degrees of freedom, \[\boldsymbol{\mathsf{K}}^{AB} = \frac{\partial \boldsymbol{\mathsf{R}}^{A}}{\partial \boldsymbol{\mathsf{u}}^{B}}.\]

Substituting the constitutive relations for the stress resultant increments yields \[\begin{align} \boldsymbol{\mathsf{K}}^{AB} = \int_{\bar{\Gamma}} \Bigg[ & \frac{\partial \epsilon_{\alpha\beta}}{\partial \boldsymbol{\mathsf{u}}^{A}} \, \mathbb{A}^{\alpha\beta\gamma\delta} \, \frac{\partial \epsilon_{\gamma\delta}}{\partial \boldsymbol{\mathsf{u}}^{B}} + \frac{\partial \epsilon_{\alpha\beta}}{\partial \boldsymbol{\mathsf{u}}^{A}} \, \mathbb{B}^{\alpha\beta\gamma\delta} \, \frac{\partial \kappa_{\gamma\delta}}{\partial \boldsymbol{\mathsf{u}}^{B}} \nonumber\\ & + \frac{\partial \kappa_{\alpha\beta}}{\partial \boldsymbol{\mathsf{u}}^{A}} \, \mathbb{B}^{\gamma\delta\alpha\beta} \, \frac{\partial \epsilon_{\gamma\delta}}{\partial \boldsymbol{\mathsf{u}}^{B}} + \frac{\partial \kappa_{\alpha\beta}}{\partial \boldsymbol{\mathsf{u}}^{A}} \, \mathbb{D}^{\alpha\beta\gamma\delta} \, \frac{\partial \kappa_{\gamma\delta}}{\partial \boldsymbol{\mathsf{u}}^{B}} \nonumber\\ & + n^{\alpha\beta} \frac{\partial^{2}\epsilon_{\alpha\beta}}{\partial \boldsymbol{\mathsf{u}}^{A}\partial \boldsymbol{\mathsf{u}}^{B}} + m^{\alpha\beta} \frac{\partial^{2}\kappa_{\alpha\beta}}{\partial \boldsymbol{\mathsf{u}}^{A}\partial \boldsymbol{\mathsf{u}}^{B}} \Bigg] d\bar{\Gamma}. \label{eq:total95stiffness95matrix} \end{align}\tag{63}\] The first four terms in Eq. 63 constitute the material stiffness arising from the constitutive tangent operators \(\mathbb{A}\), \(\mathbb{B}\) and \(\mathbb{D}\), whereas the remaining terms form the geometric stiffness associated with the current stress state. These geometric contributions are essential for accurately capturing large-deformation effects, bifurcation phenomena and snap-through instabilities.

After assembly of all element contributions, the global nonlinear equilibrium equations take the form \[\boldsymbol{\mathsf{R}}(\boldsymbol{\mathsf{u}}) = \boldsymbol{\mathsf{F}}^{\mathrm{int}}(\boldsymbol{\mathsf{u}}) - \boldsymbol{\mathsf{F}}^{\mathrm{ext}} = \boldsymbol{\mathsf{0}}, \label{eq:global95equilibrium}\tag{64}\] where \(\boldsymbol{\mathsf{u}}\) denotes the vector collecting all displacement degrees of freedom. Linearisation of Eq. 64 yields the global tangent system \[\boldsymbol{\mathsf{K}}\,\varDelta\boldsymbol{\mathsf{u}} = -\boldsymbol{\mathsf{R}}, \label{eq:newton95system}\tag{65}\] with \[\boldsymbol{\mathsf{K}} = \boldsymbol{\mathsf{K}}_{\mathrm{mat}} + \boldsymbol{\mathsf{K}}_{\mathrm{geo}},\] where \(\boldsymbol{\mathsf{K}}_{\mathrm{mat}}\) and \(\boldsymbol{\mathsf{K}}_{\mathrm{geo}}\) denote the assembled material and geometric stiffness matrices, respectively.

The resulting nonlinear system is solved iteratively using a Newton–Raphson procedure. For problems involving limit points and snap-through instabilities, the Newton iterations are combined with the arc-length continuation strategy described in the following section.

5 Bifurcation Tracking and Post-Buckling Analysis↩︎

Algorithm 4 outlines a staged arc-length procedure combined with eigenmode perturbation for bifurcation tracking. A detailed description is given in the subsequent sections.

5.1 Staged Arc-Length Procedure↩︎

To capture symmetry-breaking bifurcations and post-buckling equilibrium paths, the nonlinear equilibrium equations are solved using a Newton–Raphson scheme combined with a Crisfield spherical arc-length method. The arc-length strategy enables the solution path to be traced through limit points and snap-through instabilities that cannot be followed by conventional load-control procedures.

The external loading is parameterized by a scalar load multiplier \[\boldsymbol{\mathsf{F}}^{\mathrm{ext}}(\mathcal{L}) = \mathcal{L}\,\overline{\boldsymbol{\mathsf{F}}},\] where \(\mathcal{L}\) denotes the active loading parameter. Depending on the loading stage, \(\mathcal{L}\) represents either the pressure scaling factor \(p=\mathcal{L}p_0\) or the voltage scaling factor \(\Delta\Phi=\mathcal{L}\Delta\Phi_0.\)

The nonlinear equilibrium condition is expressed as \[\boldsymbol{\mathsf{R}}(\boldsymbol{\mathsf{u}},\mathcal{L}) = \boldsymbol{\mathsf{F}}^{\mathrm{int}}(\boldsymbol{\mathsf{u}}) - \boldsymbol{\mathsf{F}}^{\mathrm{ext}}(\mathcal{L}) = \boldsymbol{\mathsf{0}}.\]

Since both the displacement field \(\boldsymbol{\mathsf{u}}\) and the load multiplier \(\mathcal{L}\) are unknown, an additional arc-length constraint is introduced, \[g(\tilde{\varDelta}\boldsymbol{\mathsf{u}},\tilde{\varDelta}\mathcal{L}) = \tilde{\varDelta}\boldsymbol{\mathsf{u}}^{\mathsf T}\tilde{\varDelta}\boldsymbol{\mathsf{u}} + \psi^2(\tilde{\varDelta}\mathcal{L})^2 - \tilde{\varDelta} s^2 = 0,\] where \(\tilde{\varDelta} s\) is the prescribed arc length and \(\psi\) is a scaling parameter controlling the relative contribution of displacement and load increments.

Here, \(\tilde{\delta}(\cdot)\) denotes the Newton correction within the current iteration, whereas \(\tilde{\varDelta}(\cdot)\) denotes the cumulative increment over the current arc-length step.

At Newton iteration \(k\), the equilibrium equations and the arc-length constraint are linearized simultaneously, resulting in the extended system \[\begin{bmatrix} \boldsymbol{\mathsf{K}}^{(k)} & -\overline{\boldsymbol{\mathsf{F}}} \\ 2\tilde{\varDelta}\boldsymbol{\mathsf{u}}^{(k)\mathsf T} & 2\psi^2\tilde{\varDelta}\mathcal{L}^{(k)} \end{bmatrix} \begin{bmatrix} \tilde{\delta}\boldsymbol{\mathsf{u}} \\ \tilde{\delta}\mathcal{L} \end{bmatrix} = - \begin{bmatrix} \boldsymbol{\mathsf{R}}^{(k)} \\ g^{(k)} \end{bmatrix},\] where \(\boldsymbol{\mathsf{K}}^{(k)}\) is the consistent tangent stiffness matrix and \[g^{(k)} = \tilde{\varDelta}\boldsymbol{\mathsf{u}}^{(k)\mathsf T} \tilde{\varDelta}\boldsymbol{\mathsf{u}}^{(k)} + \psi^2 \left( \tilde{\varDelta}\mathcal{L}^{(k)} \right)^2 - \tilde{\varDelta} s^2\] is the residual of the arc-length constraint.

After solving for the Newton corrections, the cumulative increments are updated as \[\tilde{\varDelta}\boldsymbol{\mathsf{u}}^{(k+1)} = \tilde{\varDelta}\boldsymbol{\mathsf{u}}^{(k)} + \tilde{\delta}\boldsymbol{\mathsf{u}}, \qquad \tilde{\varDelta}\mathcal{L}^{(k+1)} = \tilde{\varDelta}\mathcal{L}^{(k)} + \tilde{\delta}\mathcal{L}.\]

The total solution is then updated according to \[\boldsymbol{\mathsf{u}}^{(k+1)} = \boldsymbol{\mathsf{u}}_{n} + \tilde{\varDelta}\boldsymbol{\mathsf{u}}^{(k+1)}, \qquad \mathcal{L}^{(k+1)} = \mathcal{L}_{n} + \tilde{\varDelta}\mathcal{L}^{(k+1)},\] where \((\boldsymbol{\mathsf{u}}_{n},\mathcal{L}_{n})\) denotes the converged solution at the previous arc-length step.

The Newton iterations continue until both the residual norm and the displacement increment satisfy the prescribed convergence tolerances.

To reproduce the experimental loading protocol and improve numerical robustness, the equilibrium path is traced through three consecutive continuation stages:

1. Pressure loading from the undeformed configuration to an intermediate pressure level \(p_{\mathrm{mid}}\) with \(\Delta\Phi=0\);

2. Voltage loading from \(0\) to \(\Delta\Phi_{\mathrm{target}}\) while maintaining the pressure fixed at \(p_{\mathrm{mid}}\);

3. Further pressure loading from \(p_{\mathrm{mid}}\) to the target pressure while keeping \(\Delta\Phi=\Delta\Phi_{\mathrm{target}}\) constant.

The converged solution of each stage serves as the initial state for the subsequent stage. This staged arc-length procedure allows the complete electro-mechanical equilibrium manifold to be followed, including symmetry-breaking bifurcations, limit points, and snap-through transitions.

5.2 Mechanism of Inducing Bifurcation↩︎

The electromechanical response of a perfectly symmetric shell under internal pressure and voltage loading is inherently susceptible to symmetry-breaking bifurcations. Due to the electro-softening effect, the membrane’s stiffness degrades as the electric potential increases. Upon reaching a critical load, the axisymmetric equilibrium path loses stability, and the structure seeks a lower-energy configuration by transitioning to a non-axisymmetric deformed state. Standard incremental solvers fail at the bifurcation point because the perfectly symmetric discretization lacks the necessary perturbation to trigger the descent onto the asymmetric branch. Several techniques exist to overcome this, including eigenmode perturbation, geometric imperfections, and material imperfections. In this work, we employ a specific numerical strategy to navigate this singularity and capture the post-buckling response.

Remarks: The perturbation is scaled to a small fraction of the shell’s characteristic dimension. However, determining the optimal scale factor is non-trivial: if too small, the solver may remain trapped on the principal branch; if too large, it may overshoot the true post-buckling path. In practice, the scale requires careful calibration based on mesh resolution and problem-specific sensitivity to ensure robust convergence onto the asymmetric branch without artificially altering the physical response.

6 Numerical Examples↩︎

This section presents three numerical examples designed to systematically validate and demonstrate the capabilities of the proposed electroelastic thin shell formulation. The examples are sequenced to progress from fundamental verification to the exploration of complex, coupled nonlinear phenomena. First, the inflation of a spherical membrane serves to benchmark the model’s accuracy against a known analytical solution. Second, the electromechanical response of a prestretched circular plate under combined pressure and voltage loading is investigated to elucidate the competitive interplay between mechanical prestress and electric-field-induced softening. Finally, the primary focus is on the nonlinear buckling and symmetry-breaking behaviour of toroidal membranes, showcasing the formulation’s ability to trace complex post-bifurcation equilibrium paths under strong electromechanical coupling. All simulations are performed using the numerical implementation based on subdivision surfaces described in Section 4.

6.1 Validation: Electroelastic Spherical Shells↩︎

To validate the proposed method, the inflation behaviour of an electroelastic spherical shell is examined against available analytical solutions. The pure-mechanical case serves as a widely adopted benchmark, as demonstrated in our previous work [39] and by others [38], [47], [49]. The analytical electroelastic solution employed here is derived by extending this benchmark to incorporate electromechanical coupling effects. The geometric setup and the corresponding control mesh are presented inFig. 5. One considers an electroelastic spherical shell of radius \(\bar{R} = 10\) and thickness \(\bar{h} = 0.1\), subjected to both internal pressure and an electrical potential. An inflating pressure \(p\) is applied to the inner surface of the shell, while different electric potentials \(\Delta\Phi \in \{0, 10, 15, 20\}\) are individually imposed across its thickness. The material parameters are taken as \(c_1 = 0.4375\mu\) and \(c_2 = 0.0625\mu\), while \(\mu = 4.225\times10^5\). For the spherical shell, the stretching of the mid-surface is the same in all directions, thus \(\lambda_1 = \lambda_2 = \lambda\). The analytical solution for the inner pressure \({p}\) is given by: \[p = \frac{4\bar{h}}{\bar{R}} \left[ c_1 [\lambda^{-1} - \lambda^{-7}] - c_2 [\lambda^{-5} - \lambda] - \frac{\varepsilon [\Delta\Phi]^2 \lambda}{2 \bar{h}^2} \right].\] Fig. 6 compares the numerical results obtained from the proposed method with the analytical solution. The plot shows the relationship between the internal pressure \(p\) and the maximum radial displacement \(u\) under different applied electric potentials. The solid lines represent the analytical prediction derived from the variational principle, while the markers denote the numerical results from the present finite element implementation. As shown in the figure, the numerical results are in excellent agreement with the analytical curves across all voltage levels. The electrostatic softening effect, manifested as a leftward and downward shift of the pressure–displacement curve with increasing \(\Delta\Phi\), is accurately captured by the numerical model. This validation confirms the correctness and robustness of the proposed electroelastic thin shell formulation for large-deformation problems involving strong electromechanical coupling.

6.2 Electro-Softening: Inflation of Prestreched Circular Plates↩︎

This example investigates the coupled response of a prestretched dielectric elastomer membrane subjected to simultaneous inflation pressure and electric potential (shown in Fig. 7). The primary objective is to analyse the competing mechanisms between mechanical prestretch, which increases the effective stiffness, and the applied electric field, which induces a pronounced electro-softening effect. A circular dielectric elastomer membrane with an initial undeformed radius \(\bar{R} = 7.5 \times 10^{-2}\) and thickness \(\bar{h} = 5\times10^{-4}\) is considered. The material parameters for the Mooney–Rivlin model are \(c_1 = 8\times10^4\), \(c_2 = 2\times10^4\).

Fig. 8 illustrates the load–displacement curves and corresponding deformation profiles of the electro-elastic circular membrane under various levels of prestrech. The electric potential \(\Delta\Phi\) increases from 0 to 15000 with a step size of 5000, while the stretch ratio \(\lambda_p\) induced by the prestress varies within the set \(\{1.1, 1.2, 1.5\}\). The pressure–displacement curve flattens as the electric potential rises, indicating a loss of structural stiffness. Accordingly, the maximum displacement at a constant internal pressure increases substantially. Introducing a prestrech elevates the membrane’s load-bearing capacity by enhancing its effective stiffness. This increased structural rigidity results in a substantial suppression of the maximum displacement. As illustrated in Fig. 8 (b), upon applying an electric potential of 15000, the maximum displacement under a given internal pressure exhibits a marked surge compared to the zero-potential scenario. This behaviour stems from the degradation of material stiffness induced by the electro-softening effect. This is attributed to the reduction in material stiffness caused by the electro-softening effect. In contrast, in Fig. 8 (d), under the same applied potential of 15000, the increase in the maximum displacement for the same internal pressure is relatively small. This indicates that, for the present numerical example, an increase in prestress suppresses the electro-softening effect.

6.3 Symmetric Breaking of Electroelastic Toroidal Membranes↩︎

The symmetry-breaking behaviour of electroelastic toroidal membranes was investigated in our previous work [50] using a semi-analytical framework. Subsequently, a simplified discrete model for axisymmetric dielectric elastomer membranes was developed [51] and successfully applied to the analysis of toroidal structures. Although the axisymmetric formulation greatly reduces the computational cost by exploiting geometric symmetry, it inherently excludes non-axisymmetric deformation modes. Consequently, neither symmetry-breaking bifurcations nor the associated post-buckling branches can be resolved. The present three-dimensional electroelastic shell formulation overcomes this limitation, providing a general framework for analysing symmetry-breaking instabilities, branch switching, and post-buckling behaviour in toroidal dielectric elastomer structures.

6.3.1 Symmetric Inflation and Limit Point Instability↩︎

We first investigate the deformation behaviour and principal equilibrium pathways under varying electric fields. A numerical example is first established to simulate the inflation process, employing material parameters \(c_1 = 0.4375\mu\) and \(c_2 = 0.0625\mu\) for a torus with major radius \(\bar{R}_b = 10\), minor radius \(\bar{R}_s = 2\), and thickness \(\bar{h} = 0.01\). The control mesh and limiting surface are shown in Fig. 9.

Figure 10 (a) illustrates the variation of the pressure \(p\) with respect to the maximum displacement \(u_m\) for different electrical loads. The results indicate that the application of an electric load effectively reduces the load-bearing capacity of the structure.

From a stability perspective, the limit point corresponds to the loss of tangent stiffness and the onset of snap-through instability. Under pure mechanical loading, the limit pressure \(p_s\) marks the transition from a stable equilibrium branch to an unstable post-buckling path. When an electric field is applied, the Maxwell stress introduces an additional tensile component in the membrane, which reduces the effective structural stiffness and shifts the equilibrium path towards lower pressures. Consequently, the limit point moves to smaller displacements and lower critical loads, indicating that the electric field acts as a destabilising agent.

Furthermore, Figure 10 (b) depicts the variation of the limit-point pressure \(p_s\) with the aspect ratio \(\bar{R}_s/\bar{R}_b\). It is observed that \(p_s\) decreases monotonically with increasing aspect ratio and that higher electrical loads further suppress the critical buckling pressure. These numerical observations are in excellent agreement with existing analytical benchmarks [51].

6.3.2 Electromechanical Softening and the Role of Material Nonlinearity↩︎

The influence of material nonlinearity on the electromechanical response is investigated by varying the Mooney–Rivlin parameter ratio \(c_2/c_1 \in \{0, 0.1, 0.2, 0.3\}\) under electric potentials \(\Delta\Phi \in \{0, 10000, 15000, 20000\}\). Figure 11 illustrates the resulting pressure–displacement curves. In the absence of an electric field (Fig. 11 (a)), an increase in \(c_2/c_1\) enhances the post-limit-point stiffening, indicative of a stronger strain-hardening response at large deformations. However, the application of an electric load exerts a destabilising effect, progressively suppressing this strain-hardening behaviour. For instance, at \(c_2/c_1 = 0.1\), increasing \(\Delta\Phi\) from \(0\) to \(20000\) gradually erodes the material’s hardening capacity until it vanishes entirely.

6.3.3 Electrically Induced Localisation and Post‑Buckling Path Switching↩︎

Tracing the post-buckling equilibrium path of a geometrically perfect toroidal membrane beyond the bifurcation point presents a significant computational challenge due to the loss of uniqueness at the critical state. To circumvent this numerical instability, a symmetry-breaking perturbation, either in the form of an eigenmode or a geometric imperfection, is strategically imposed. This artificial imperfection serves as a trigger, stabilising the nonlinear solver and enabling controlled access to the post-bifurcation branches.

Figures 12 and 13 (a) depict the pressure–volume response alongside the corresponding deformation modes. The bifurcation points are observed to emerge in proximity to the limit point, with the post-buckling trajectories intersecting the principal equilibrium path. For the slender torus (\(\bar{R}_s/\bar{R}_b = 0.2\)) under high electric potential (\(\Delta\Phi = 20000\)), intersection is suppressed due to severe localisation: excessive deformation on one side induces self-contact, preventing further volume expansion along the symmetric branch. Conversely, for the thicker torus (\(\bar{R}_s/\bar{R}_b = 0.4\)) without an electric field, the asymmetric response is marginal, resulting in negligible deviation from the principal path.

Critically, under identical volumetric constraints, increased electrical loading reduces the internal pressure, underscoring the dominance of electrostatic softening over geometric stiffening in the post-buckling regime. Furthermore, the electric field modulates the spatial extent of the localised mode: higher voltages confine the deformation to a smaller region, intensifying the strain localisation.

7 Conclusions↩︎

This manuscript has presented a comprehensive isogeometric framework for the analysis of electroelastic thin shells, with a particular focus on capturing symmetry-breaking instabilities and post-buckling behaviour under strong electromechanical coupling. The key contributions are summarised as follows:

1. A nonlinear Kirchhoff–Love shell formulation for dielectric elastomers was developed, incorporating finite deformation kinematics, electromechanical coupling, and the plane stress condition. The total stress was decomposed into hyperelastic, Maxwell, and hydrostatic contributions, and the consistent tangent moduli were derived in closed form with static condensation to eliminate the thickness stress component.

2. The numerical implementation employs Catmull–Clark subdivision surfaces, providing the \(C^1\)-continuous basis functions required by Kirchhoff–Love shell theory, which enables accurate geometric representation and smooth deformation fields essential for capturing complex instability patterns.

3. By employing eigenmode perturbation to break symmetry, the evolution from axisymmetric to non-axisymmetric deformation modes was successfully captured, and the staged arc-length procedure allowed robust tracing of post-buckling equilibrium paths.

4. Numerical examples on spherical, prestretched circular plate, and toroidal membranes validated the framework’s ability to capture large deformations, bifurcation onset, mode switching, and post-buckling responses under coupled electromechanical loading.

Three numerical examples validated the proposed framework. The spherical membrane benchmark confirmed excellent agreement with the analytical solution across multiple potential differences, demonstrating the correct implementation of the electromechanical coupling. The prestretched circular plate example revealed the competitive interplay between mechanical prestress and electro-softening: higher prestretch suppresses the voltage-induced softening effect, while larger electric potentials significantly reduce the membrane’s load-bearing capacity. The toroidal membrane example showcased the framework’s unique capability to trace post-bifurcation equilibrium paths under electromechanical loading. The results demonstrated that electrical loading not only reduces the limit point pressure but also influences the post-buckling deformation pattern, with larger electric potentials producing more localised asymmetric deformations.

The proposed isogeometric framework provides a robust and accurate tool for analysing the complex nonlinear behaviour of dielectric elastomer thin shells. Its ability to handle large deformations, electromechanical coupling, and symmetry-breaking instabilities makes it well-suited for the design and optimisation of soft actuators, sensors, and energy harvesters. Future work will focus on extending the framework to incorporate dynamic effects, viscoelastic material behaviour, and multi-layer dielectric elastomer laminates.

References↩︎