An analytical model for nanoscale flexoelectric doubly curved shells

Abstract

The phenomenon of flexoelectricity, wherein the generation of electrical polarization results from strain gradients, has garnered significant attention in electromechanics. This phenomenon holds immense potential for diverse applications in nanoelectromechanical systems. Doubly curved shells present a particularly compelling structure for implementing flexoelectric devices, representing the most general scenario. To comprehensively understand the mechanical behavior of nanoscale flexoelectric doubly curved shells, we present a novel analytical model in this study. Our model incorporates a reformulated strain gradient elasticity theory, a modified expression of electric enthalpy density considering Maxwell’s self-field gradient, and a moderately thick shell configuration. To further elucidate this model, we establish complete formulations and analytical solutions of static bending and free vibration problems. This study reveals the crucial roles that various parameters such as principal curvature radii, length scale parameters, transverse structural sizes, and the Winkler–Pasternak foundation play in controlling the mechanical behavior of nanoscale flexoelectric doubly curved shells. These theoretical findings offer valuable insights into the design of nanodevices based on the direct flexoelectric effect.

Keywords

Flexoelectricity doubly curved shells strain gradient nanoscale size effect

1. Introduction

Flexoelectricity is a size-dependent phenomenon that has garnered increasing attention in recent years [1 –6]. It encompasses two distinct phenomena, namely, the direct and the converse flexoelectric effect [7 –10]. The direct flexoelectric effect bridges strain gradient and electrical polarization, which is applied in energy harvesters [11 –14] and sensors [15 –18]. On the contrary, the converse flexoelectric effect links mechanical stresses and electric field gradients, which has a high practical value in actuators [19 –22]. The flexoelectric effect breaks the limitation that piezoelectricity solely occurs in crystals which are not centrosymmetric [2], thus significantly broadening the range of materials suitable for electromechanical coupling. In addition, the flexoelectric effect outweighs piezoelectricity in nanoscale structures due to the considerable strain and polarization gradient [23]. As a result, flexoelectricity is more widespread and dominant than piezoelectricity at the nanoscale and merits further exploration.

As a classic electromechanical coupling effect, piezoelectricity has been studied by many scholars systematically [24 –30]. The intriguing electromechanical coupling phenomenon of flexoelectricity has spurred many theoretical studies in the past few years. Previous research by Mao and Purohit [31] employed strain gradient elasticity theory to gain insights into flexoelectric solids. They analytically solved simple problems such as a disk under pressure and the in-plane shear of a disk. Deng et al. [13] formulated a model for flexoelectric slender beams, which enabled a thorough investigation of the size effect in nanoscale flexoelectric energy harvesters. Li et al. [32] put forward a modified flexoelectric theory that applied to isotropic dielectrics. Their work differentiated the independent contributions of individual strain gradient components. Using the Kirchhoff–Love assumption and simplified strain gradient elasticity theory, Tzou and Zhang [33] examined the direct flexoelectric effect of energy harvesters appearing in the forms of thin doubly curved shells. Furthermore, Qi et al. [34 –36] investigated the mechanical behaviors of flexoelectric curved microbeams, spherical microshells, and annular microplates based on a modified expression of the electric enthalpy density and strain gradient elasticity theory. Finally, Qu et al. [37] developed a circular cylindrical Kirchhoff–Love shell model utilizing couple stress theory and provided analytical solutions to both static and dynamic problems.

To date, there is a scarcity of research on nanoscale flexoelectric shells. This study proposes a novel approach to model such shells by incorporating a reformulated strain gradient elasticity theory, a modified electric enthalpy density expression, all within the context of Reissner–Mindlin shell theory. The formulation possesses several advantages over other flexoelectric shell models. First, the reformulated theory of strain gradient elasticity is highly unified, where the microstructure effect is only characterized by three strain gradient parameters. Compared with the strain gradient elasticity theory of Mindlin where there are 16 intrinsic length constants [38], this theory is much simpler while without loss of any generality. This is because it is based on pure mathematical reformulations without additional conditions and allows for the treatment of the couple stress effect simultaneously. Second, this modified electric enthalpy density expression considers the effect of electric field gradients, providing a more comprehensive and reasonable approach to studying flexoelectric effects. Finally, the proposed flexoelectric shell model on the basis of the Reissner–Mindlin shell theory is suitable for analyzing moderately thick shells in addition to Kirchhoff–Love thin shells. The double-curvature property allows for a comparison of the electromechanical behaviors of flexoelectric shells belonging to various geometries.

This paper’s principal objective is to study the electromechanical behavior of nanoscale flexoelectric doubly curved shells in a systematic and general approach. Driven by this motivation, we formulate the nanoscale flexoelectric doubly curved shell model using Hamilton’s principle in section 2. In section 3, we perform a static bending analysis by applying the newly developed model. In section 4, we conduct a free vibration analysis. Finally, we conclude this study in section 5.

Concerning notation, we use the common Cartesian index notation and employ summation convention where applicable. Latin indices take on the values ${1, 2, 3}$ , and partial derivatives will be denoted by a comma.

2. Formulation of flexoelectric doubly curved shells

2.1. Kinematic relations

The geometric configuration of a Reissner–Mindlin double-curvature shell with elastic foundation and its related coordinate system are depicted in Figure 1. We assume the shell has a uniform thickness $h$ and possesses two constant principal curvature radii, $R_{1}$ and $R_{2}$ , respectively. The boundary planes are normal to the middle surface. This coordinate system is defined by ${θ_{1}, θ_{2}, z}$ positioning coordinates. The two longitudinal axes are represented by $θ_{1}$ and $θ_{2}$ , with the dimensions measured in angles. The shell thickness is oriented along the $z$ direction, with the unit of measurement being length. The parameters $α_{0}$ and $β_{0}$ denote the dimensions of the doubly curved shells in the two longitudinal directions.

In the context of doubly curved shells with orthogonal curvilinear coordinates, the length of a microelement is characterized as:

\begin{matrix} {(ds)}^{2} = g_{11} {(d θ_{1})}^{2} + g_{22} {(d θ_{2})}^{2} + g_{33} {(dz)}^{2} . \end{matrix}

(1)

The covariant metric tensor components are represented by $g_{11}$ , $g_{22}$ , and $g_{33}$ , and their definitions are:

\begin{matrix} g_{11} = R_{1}^{2} {(1 + \frac{z}{R_{1}})}^{2}, g_{22} = R_{2}^{2} {(1 + \frac{z}{R_{2}})}^{2}, g_{33} = 1 . \end{matrix}

(2)

The second type of Christoffel symbol in this coordinate system is expressed as:

\begin{matrix} Γ_{13}^{1} = Γ_{31}^{1} = \frac{1}{R_{1} (1 + z / R_{1})}, Γ_{23}^{2} = Γ_{32}^{2} = \frac{1}{R_{2} (1 + z / R_{2})}, \\ Γ_{11}^{3} = - R_{1} (1 + \frac{z}{R_{1}}), Γ_{22}^{3} = - R_{2} (1 + \frac{z}{R_{2}}), \end{matrix}

(3)

while all other occurrences of $Γ_{ij}^{k}$ are nullified.

The displacement field is formulated in accordance with the shell theory of Reissner–Mindlin, which assumes two independent rotation components. This displacement field can be represented as:

U_{1} (θ_{1}, θ_{2}, z, t) = u_{1} (θ_{1}, θ_{2}, t) + z ϕ_{1} (θ_{1}, θ_{2}, t),

(4)

U_{2} (θ_{1}, θ_{2}, z, t) = u_{2} (θ_{1}, θ_{2}, t) + z ϕ_{2} (θ_{1}, θ_{2}, t),

(5)

U_{3} (θ_{1}, θ_{2}, z, t) = u_{3} (θ_{1}, θ_{2}, t),

(6)

where $U_{1}, U_{2},$ and $U_{3}$ represent the displacements along the $θ_{1}$ -, $θ_{2}$ -, and $z$ -axes, respectively. $u_{1}, u_{2},$ and $u_{3}$ correspond to the displacements of the middle surface, while the remaining two components, $ϕ_{1}, ϕ_{2}$ , are the rotation components around the $θ_{2}$ - and $θ_{1}$ -axes, respectively. These five basic displacement components, $u_{1}, u_{2}, u_{3}, ϕ_{1}, ϕ_{2}$ are the functions of only $θ_{1}, θ_{2}, t$ .

The strain gradient tensor contains three independent portions [39]: the deviatoric stretch and rotation parts as well as the gradient of dilatation, which are calculated in the following procedure. First, the classical strain tensor and the strain gradient tensor are:

ε_{ij} = \frac{1}{2} (u_{i, j} + u_{j, i}),

(7)

η_{ijk} = ε_{jk, i} = \frac{1}{2} (u_{j, ki} + u_{k, ji}) .

(8)

Furthermore, the quantities $η_{ijk}^{h}$ , $η_{ijk}^{'}$ , and $η_{ijk}^{s}$ represent the tensors that describe the hydrostatic part, deviatoric part, and symmetric part, respectively. These tensors can be expressed as:

η_{ijk}^{h} = \frac{1}{3} δ_{jk} η_{inn},

(9)

η_{ijk}^{'} = η_{ijk} - \frac{1}{3} η_{ijk}^{h},

(10)

η_{ijk}^{s} = \frac{1}{3} (η_{ijk} + η_{jki} + η_{kij}) .

(11)

The deviatoric stretch tensor is denoted by $η_{ijk}^{(1)}$ , the dilatation gradient vector is represented by $γ_{i}$ , and the deviatoric rotation tensor is indicated by $χ_{ij}^{'}$ . These quantities are defined as:

η_{ijk}^{(1)} = η_{ijk}^{s} - \frac{1}{5} (δ_{ij} η_{mmk}^{s} + δ_{jk} η_{mmi}^{s} + δ_{ki} η_{mmj}^{s}),

(12)

γ_{i} = ε_{nn, i},

(13)

χ_{ij}^{'} = e_{ipq} η_{pqj}^{'} .

(14)

Within the context of doubly curved shells, expressed in terms of orthogonal curvilinear coordinates, the first-order and second-order displacement gradients are characterized by their physical components [40]:

\begin{matrix} u_{i} |_{j} = \frac{1}{\sqrt{g_{ii} g_{jj}}} [{(\sqrt{g_{ii}} u_{i})}_{, j} - \sum_{k = 1}^{3} \sqrt{g_{kk}} u_{k} Γ_{ij}^{k}] (No sum on i, j), \end{matrix}

(15)

\begin{matrix} u_{i} |_{jk} = \frac{1}{\sqrt{g_{ii} g_{jj} g_{kk}}} [{(\sqrt{g_{ii} g_{jj}} u_{i} |_{j})}_{, k} - \sum_{r = 1}^{3} \sqrt{g_{ii} g_{rr}} u_{i} |_{r} Γ_{jk}^{r} - \sum_{r = 1}^{3} \sqrt{g_{rr} g_{jj}} u_{r} |_{j} Γ_{ik}^{r}] \\ (No sum on i, j, k) . \end{matrix}

(16)

The covariant derivative is represented by the vertical line. It is noteworthy to mention that in the above equations, the notation $u_{i}$ corresponds to the set of variables ${U_{1}, U_{2}, U_{3}}$ . Neglecting higher-order components and acknowledging that $1 / (1 + z / R_{1}) \approx 1$ and $1 / (1 + z / R_{2}) \approx 1$ , the strain components within the context of doubly curved shells, expressed in terms of orthogonal curvilinear coordinates, are:

ε_{11} = \frac{1}{R_{1}} u_{1, 1} + \frac{1}{R_{1}} u_{3} + \frac{z}{R_{1}} ϕ_{1, 1},

(17)

ε_{22} = \frac{1}{R_{2}} u_{2, 2} + \frac{1}{R_{2}} u_{3} + \frac{z}{R_{2}} ϕ_{2, 2},

(18)

ε_{33} = 0,

(19)

ε_{23} = ε_{32} = \frac{1}{2} (ϕ_{2} + \frac{1}{R_{2}} u_{3, 2} - \frac{1}{R_{2}} u_{2}),

(20)

ε_{13} = ε_{31} = \frac{1}{2} (ϕ_{1} + \frac{1}{R_{1}} u_{3, 1} - \frac{1}{R_{1}} u_{1}),

(21)

ε_{12} = ε_{21} = \frac{1}{2} (\frac{1}{R_{1}} u_{2, 1} + \frac{1}{R_{2}} u_{1, 2}) + \frac{z}{2} (\frac{1}{R_{1}} ϕ_{2, 1} + \frac{1}{R_{2}} ϕ_{1, 2}) .

(22)

The subscripts $(),_{1}$ , $(),_{2}$ , and $(),_{3}$ denote the partial derivatives with respect to $θ_{1}, θ_{2},$ and $z$ , respectively. The independent strain gradient components $η_{ijk}^{(1)}$ , $γ_{i}$ , and $χ_{ij}^{'}$ can be computed accordingly, as presented in Appendix 1. In addition, an uneven distribution of $ε_{13}$ and $ε_{23}$ throughout the shell thickness necessitates the introduction of a shear correction factor, denoted as $k_{s}$ . Consequently, the thickness shear strains $ε_{13}$ and $ε_{23}$ are replaced by $k_{s} ε_{13}$ and $k_{s} ε_{23}$ , respectively [41]. Therefore, we have obtained the complete formulation of the kinematic relations of doubly curved shell structures. In the next section, the constitutive relations of isotropic flexoelectric solids will be discussed.

2.2. Constitutive relations

In accordance with the linear flexoelectric theory, the internal energy density is related to several variables, including strain, polarization, and their respective gradients [42]:

\begin{matrix} U = U (ε_{ij}, η_{ijk}, P_{i}, Q_{ij}), \end{matrix}

(23)

where the electric polarization is indicated by $P_{i}$ , and $Q_{ij}$ represents the polarization gradient.

To separate independent strain gradient components, a reformulated theory of flexoelectricity for isotropic dielectrics is introduced [32]. Accordingly, the function of internal energy density is:

\begin{matrix} U = \frac{1}{2} σ_{ij} ε_{ij} + \frac{1}{2} τ_{ijk}^{(1)} η_{ijk}^{(1)} + \frac{1}{2} p_{i} γ_{i} + \frac{1}{2} m_{ij}^{'} χ_{ij}^{'} + \frac{1}{2} E_{i} P_{i} + \frac{1}{2} V_{ij} Q_{ij} \\ = \frac{1}{2} k ε_{ii} ε_{jj} + μ ε_{ij}^{'} ε_{ij}^{'} + μ l_{0}^{2} γ_{i} γ_{i} + μ l_{1}^{2} η_{ijk}^{(1)} η_{ijk}^{(1)} + μ (l_{2}^{2} + \frac{9}{5} l_{0}^{2}) χ_{ij}^{'} χ_{ij}^{'} + μ (l_{2}^{2} - \frac{9}{5} l_{0}^{2}) χ_{ij}^{'} χ_{ji}^{'} \\ + \frac{1}{2} a (P_{i} P_{i} + β_{1}^{2} Q_{ii} Q_{jj} + β_{2}^{2} Q_{ij} Q_{ij} + β_{3}^{2} Q_{ij} Q_{ji}) \\ + (f_{1} + \frac{2}{3} f_{2}) P_{i} γ_{i} - f_{1} Q_{ii} ε_{jj} + 2 f_{2} (P_{i} e_{ijk} χ_{jk}^{'} - Q_{ij} ε_{ij}) . \end{matrix}

(24)

In the aforementioned equation, the classical stress tensor is denoted by $σ_{ij}$ , while $E_{i}$ and $V_{ij}$ are the local and higher-order local electric fields. In addition, $τ_{ijk}^{(1)}$ , $p_{i}$ , and $m_{ij}^{'}$ are defined as conjugated non-classical stresses. The material parameters $k$ , $μ$ , and $f_{i}$ correspond to the bulk volume modulus, shear modulus, and flexocoupling coefficients, respectively. The electric parameter $a$ represents the reciprocal dielectric susceptibility. Furthermore, the equation involves two sets of material length parameters: $l_{i}$ indicates the intrinsic strain gradient length parameters, while $β_{i}$ represents the material length parameters relevant to electric polarization.

The constitutive relations for flexoelectric isotropic dielectrics can be obtained using the internal energy density function presented above:

σ_{ij} = \frac{\partial U}{\partial ε_{ij}} = k ε_{nn} δ_{ij} + 2 μ ε_{ij}^{'} - f_{1} Q_{kk} δ_{ij} - f_{2} (Q_{ij} + Q_{ji}),

(25)

τ_{ijk}^{(1)} = \frac{\partial U}{\partial η_{ijk}^{(1)}} = 2 μ l_{1}^{2} η_{ijk}^{(1)},

(26)

p_{i} = \frac{\partial U}{\partial γ_{i}} = 2 μ l_{0}^{2} γ_{i} + (f_{1} + \frac{2}{3} f_{2}) P_{i},

(27)

m_{ij}^{'} = \frac{\partial U}{\partial χ_{ij}^{'}} = 2 μ (l_{2}^{2} + \frac{9}{5} l_{0}^{2}) χ_{ij}^{'} + 2 μ (l_{2}^{2} - \frac{9}{5} l_{0}^{2}) χ_{ji}^{'} + 2 f_{2} e_{ijk} P_{k},

(28)

E_{i} = \frac{\partial U}{\partial P_{i}} = a P_{i} + (f_{1} + \frac{2}{3} f_{2}) γ_{i} + 2 f_{2} e_{ijk} χ_{jk}^{'},

(29)

V_{ij} = \frac{\partial U}{\partial Q_{ij}} = a (β_{1}^{2} Q_{nn} δ_{ij} + β_{2}^{2} Q_{ij} + β_{3}^{2} Q_{ji}) - f_{1} ε_{nn} δ_{ij} - 2 f_{2} ε_{ij} .

(30)

The Reissner–Mindlin theory suggests that the constitutive relations for stress components induced by mechanical loads, such as $σ_{11}$ , $σ_{22}$ , and $σ_{12}$ , should have the same form as that of the plane stress problem. This is because the influence of $σ_{33}$ on the strain components is neglected. Therefore, the constitutive relationships of $σ_{11}$ and $σ_{22}$ should follow the reduced material stiffness for plane stress and can be expressed as:

σ_{11} = \frac{E}{1 - ν^{2}} (ε_{11} + ν ε_{22}) - f_{1} Q_{kk} - 2 f_{2} Q_{11} (Sum on k),

(31)

σ_{22} = \frac{E}{1 - ν^{2}} (ν ε_{11} + ε_{22}) - f_{1} Q_{kk} - 2 f_{2} Q_{22} (Sum on k),

(32)

where $E$ is Young’s modulus and $ν$ is Poisson’s ratio. So far, the constitutive equations of flexoelectric solids have been fully obtained. In the next section, we will utilize the above kinematic and constitutive relations to derive the complete governing equations of flexoelectric doubly curved shells.

2.3. Governing equations

In order to derive the complete mechanical and electrical governing equations of flexoelectric doubly curved shells, Hamilton’s variational principle is utilized:

\begin{matrix} δ \int_{t_{1}}^{t_{2}} (T - H + W) dt = 0, \end{matrix}

(33)

where $T, H,$ and $W$ refer to total kinematic energy, electric enthalpy, and external work, respectively. The equilibrium equations of flexoelectric doubly curved shells can be formulated with the kinematic and constitutive relationships mentioned before. In this regard, the total kinematic energy is:

\begin{matrix} T = \int_{V} \frac{1}{2} ρ (z) [({\overset{\cdot}{u}}_{1} + z {\overset{\cdot}{ϕ}}_{1})^{2} + {({\overset{\cdot}{u}}_{2} + z {\overset{\cdot}{ϕ}}_{2})}^{2} + {\overset{\cdot}{u}}_{3}^{2}] dV . \end{matrix}

(34)

In the presence of polarization gradients, the electric enthalpy density must account for Maxwell’s self-field gradients, since the total electrical Gibbs free energy of flexoelectric dielectrics incorporates the contribution of macroscopic electrical field gradients [43]. Consequently, a revised expression for the electric enthalpy density is proposed as [36]:

\begin{matrix} H = \int_{V} {U - \frac{1}{2} ε_{0} [φ_{, i} φ_{, i} + β_{1}^{2} φ_{, ii} φ_{, jj} + (β_{2}^{2} + β_{3}^{2}) φ_{, ij} φ_{, ij}] \\ + [φ_{, i} P_{i} + β_{1}^{2} φ_{, ii} P_{j, j} + \frac{1}{2} (β_{2}^{2} + β_{3}^{2}) φ_{, ij} (P_{i, j} + P_{j, i})]} dV, \end{matrix}

(35)

in which $ε_{0}$ denotes the permittivity of vacuum, and $φ$ represents the Maxwell self-field, also known as the macroscopic electric field. We have to note that the gradient of these electrical variables is also covariant derivatives. Here, it is assumed that the polarization vector $P_{i}$ depends on $θ_{1}, θ_{2}, z,$ and $t$ , and the electric variables in the orthogonal curvilinear coordinates are expressed as:

\begin{matrix} Q_{ij} = P_{i} |_{j} = \frac{1}{\sqrt{g_{ii} g_{jj}}} [{(\sqrt{g_{ii}} P_{i})}_{, j} - \sum_{k = 1}^{3} \sqrt{g_{kk}} P_{k} Γ_{ij}^{k}] (No sum on i, j), \end{matrix}

(36)

\begin{matrix} φ |_{ij} = \frac{1}{\sqrt{g_{ii} g_{jj}}} [φ_{, ij} - \sum_{k = 1}^{3} φ_{, k} Γ_{ij}^{k}] (No sum on i, j) . \end{matrix}

(37)

We have to note that nanodevices are usually highly sensitive to external factors, such as vibrations and mechanical stresses. However, the Winkler–Pasternak foundation can serve as an effective approach to provide a stable base for nanodevices and mitigate the external disturbances. Therefore, considering the influence of the Winkler–Pasternak foundation, the comprehensive external work is delineated in the existing literature [44, 45]:

\begin{matrix} W = W_{d} + W_{b} + W_{s} + W_{f}, \end{matrix}

(38)

where $W_{d}, W_{b}, W_{s},$ and $W_{f}$ signify the work attributed to body forces, classical and higher-order boundary forces, higher-order sharp-edge line forces, and the impact of the Winkler–Pasternak foundation, respectively. The explicit expressions of these external work constituents are:

W_{d} = \int_{0}^{α_{0}} \int_{0}^{β_{0}} [b_{1} u_{1} + b_{2} u_{2} + b_{3} u_{3}] R_{1} R_{2} d θ_{2} d θ_{1},

(39)

\begin{matrix} W_{b} = \int_{0}^{β_{0}} [{\bar{N}}_{11} u_{1} + {\hat{T}}_{11} u_{1, 1} + {\bar{N}}_{12} u_{2} + {\hat{T}}_{12} u_{2, 1} + {\bar{N}}_{13} u_{3} + {\hat{T}}_{13} u_{3, 1} \\ + {\bar{M}}_{11} ϕ_{1} + {\hat{S}}_{11} ϕ_{1, 1} + {\bar{M}}_{12} ϕ_{2} + {\hat{S}}_{12} ϕ_{2, 1}] {R_{2} d θ_{2} |}_{0}^{α_{0}} \\ + \int_{0}^{α_{0}} [{\bar{N}}_{21} u_{1} + {\hat{T}}_{21} u_{1, 2} + {\bar{N}}_{22} u_{2} + {\hat{T}}_{22} u_{2, 2} + {\bar{N}}_{23} u_{3} + {\hat{T}}_{23} u_{3, 2} \\ + {\bar{M}}_{21} ϕ_{1} + {\hat{S}}_{21} ϕ_{1, 2} + {\bar{M}}_{22} ϕ_{2} + {\hat{S}}_{22} ϕ_{2, 2}] {R_{1} d θ_{1} |}_{0}^{β_{0}}, \end{matrix}

(40)

W_{s} = \sum [{\hat{N}}_{s 1} u_{1} + {\hat{N}}_{s 2} u_{2} + {\hat{N}}_{s 3} u_{3} + {\hat{M}}_{s 1} ϕ_{1} + {\hat{M}}_{s 2} ϕ_{2}],

(41)

W_{f} = \int_{0}^{α_{0}} \int_{0}^{β_{0}} [- (k_{1} u_{3} - k_{2} \nabla^{2} u_{3}) u_{3}] R_{1} R_{2} d θ_{2} d θ_{1},

(42)

where $k_{1}$ and $k_{2}$ are the two parameters of the Winkler–Pasternak foundation, ${b_{i}}$ refer to classical body forces, ${{\bar{N}}_{ij}, {\bar{M}}_{ij}}$ represent the classical boundary forces, and ${{\hat{T}}_{ij}, {\hat{S}}_{ij}}$ are the higher-order boundary forces. In addition, ${{\hat{N}}_{si}, {\hat{M}}_{si}}$ denote the higher-order sharp-edge line forces. By utilizing the variation principle for flexoelectricity, with the arbitrariness of $δ u_{i}$ , $δ ϕ_{i}$ , $δ P_{i}$ , and $δ φ$ , the mechanical governing equations are:

δ u_{1} : \frac{1}{R_{1}} N_{11, 1} + \frac{1}{R_{1}} N_{13} + \frac{1}{R_{2}} N_{12, 2} + {\tilde{T}}_{1} + b_{1} = m_{0} {\overset{\cdot\cdot}{u}}_{1} + m_{1} {\overset{\cdot\cdot}{ϕ}}_{1},

(43)

δ u_{2} : \frac{1}{R_{2}} N_{22, 2} + \frac{1}{R_{2}} N_{23} + \frac{1}{R_{1}} N_{12, 1} + {\tilde{T}}_{2} + b_{2} = m_{0} {\overset{\cdot\cdot}{u}}_{2} + m_{1} {\overset{\cdot\cdot}{ϕ}}_{2},

(44)

δ u_{3} : \frac{1}{R_{1}} N_{13, 1} + \frac{1}{R_{2}} N_{23, 2} - \frac{1}{R_{1}} N_{11} - \frac{1}{R_{2}} N_{22} + {\tilde{T}}_{3} + b_{3} - k_{1} u_{3} + k_{2} \nabla^{2} u_{3} = m_{0} {\overset{\cdot\cdot}{u}}_{3},

(45)

δ ϕ_{1} : \frac{1}{R_{1}} M_{11, 1} + \frac{1}{R_{2}} M_{12, 2} - N_{13} + {\tilde{S}}_{1} = m_{1} {\overset{\cdot\cdot}{u}}_{1} + m_{2} {\overset{\cdot\cdot}{ϕ}}_{1},

(46)

δ ϕ_{2} : \frac{1}{R_{2}} M_{22, 2} + \frac{1}{R_{1}} M_{12, 1} - N_{23} + {\tilde{S}}_{2} = m_{1} {\overset{\cdot\cdot}{u}}_{2} + m_{2} {\overset{\cdot\cdot}{ϕ}}_{2},

(47)

where ${{\tilde{T}}_{i}, {\tilde{S}}_{i}}$ are the higher-order terms related to strain gradient elasticity, which are presented in Appendix 2 in detail. Here, the equivalent masses ${m_{0}, m_{1}, m_{2}}$ are defined as:

\begin{matrix} {m_{0}, m_{1}, m_{2}} = \int_{- \frac{h}{2}}^{\frac{h}{2}} {ρ (z), z ρ (z), z^{2} ρ (z)} dz . \end{matrix}

(48)

In addition, classical and non-classical mechanical resultants are written as:

{N_{ij}, N_{i}^{p}, N_{ijk}^{τ}, N_{ij}^{m}} = \int_{- \frac{h}{2}}^{\frac{h}{2}} {σ_{ij}, p_{i}, τ_{ijk}^{(1)}, m_{ij}^{'}} dz,

(49)

{M_{ij}, M_{i}^{p}, M_{ijk}^{τ}, M_{ij}^{m}} = \int_{- \frac{h}{2}}^{\frac{h}{2}} {z σ_{ij}, z p_{i}, z τ_{ijk}^{(1)}, {zm}_{ij}^{'}} dz .

(50)

The higher-order mechanical resultants are incorporated in ${{\tilde{T}}_{i}, {\tilde{S}}_{i}}$ . Furthermore, the electrical governing equations are:

δ φ : ε_{0} (\frac{1}{R_{1}^{2}} φ_{, 11} + \frac{1}{R_{2}^{2}} φ_{, 22} + φ_{, 33}) - \frac{1}{R_{1}} P_{1, 1} - \frac{1}{R_{2}} P_{2, 2} - P_{3, 3} + \tilde{D} = 0,

(51)

δ P_{1} : E_{1} - \frac{1}{R_{1}} V_{11, 1} - \frac{1}{R_{2}} V_{12, 2} - V_{13, 3} - \frac{1}{R_{1}} V_{31} + \frac{1}{R_{1}} φ_{, 1} + {\tilde{E}}_{1} = 0,

(52)

δ P_{2} : E_{2} - \frac{1}{R_{1}} V_{21, 1} - \frac{1}{R_{2}} V_{22, 2} - V_{23, 3} - \frac{1}{R_{2}} V_{32} + \frac{1}{R_{2}} φ_{, 2} + {\tilde{E}}_{2} = 0,

(53)

δ P_{3} : E_{3} - \frac{1}{R_{1}} V_{31, 1} - \frac{1}{R_{2}} V_{32, 2} - V_{33, 3} + \frac{1}{R_{1}} V_{11} + \frac{1}{R_{2}} V_{22} + φ_{, 3} + {\tilde{E}}_{3} = 0,

(54)

where ${\tilde{D}, {\tilde{E}}_{i}}$ are the higher-order terms related to the gradient of electric field.

This model represents a comprehensive framework that encompasses numerous specific cases. It is grounded on several fundamental theories such as Reissner–Mindlin shell theory, reformulated strain gradient elasticity theory, modified electric enthalpy expression, and general doubly curved shell geometry. The broad applicability of this model is a significant advantage as it allows for the inclusion of various physical phenomena. Consequently, the model’s ability to account for multiple scenarios and configurations makes it a valuable contribution to the field of flexoelectricity. This model’s relationship with other models currently found in the literature is summarized as follows:

When the internal energy does not account for polarization and polarization gradients, the model for isotropic doubly curved shells simplifies to a Reissner–Mindlin shell theory combined with a strain gradient elasticity theory, assuming isotropy and a doubly curved geometry.

When the length scale parameters $l_{0} = l_{1} = l_{2} = 0$ are assumed in conjunction with the simplified model described above, the resulting model corresponds to an isotropic doubly curved shell based on the Reissner–Mindlin shell theory. Interestingly, if $l_{0} = l_{1} = 0$ , it represents the modified couple stress theory scenario, resulting in another non-classical model for doubly curved shells.

By imposing the constraint that the rotations $ϕ_{1}$ and $ϕ_{2}$ are related to the displacement gradients in a specific way, i.e., $ϕ_{1} = - (u_{3, 1} - u_{1}) / R_{1}$ and $ϕ_{2} = - (u_{3, 2} - u_{2}) / R_{2}$ , the simplified model reduces to the classical isotropic thin Kirchhoff–Love doubly curved shells.

By varying the values of $R_{1}$ and $R_{2}$ , the simplified model can describe different thin shell geometries, including cylindrical and spherical shells, as well as the Kirchhoff plates.

3. Static bending analysis

3.1. Boundary conditions

This research is primarily focused on investigating the direct flexoelectric response of doubly curved shells that are simply supported. Prior to examining the mechanical behaviors of such structures, it is crucial to establish the pertinent boundary conditions that correspond to the direct flexoelectric effect.

Here, a simply supported homogeneous flexoelectric doubly curved shell is considered as this type is most commonly encountered in practice. To simplify the analysis, only $P_{3}, Q_{33}, φ_{, 3}$ , and $φ_{, 33}$ are taken into consideration for the electrical enthalpy density, as suggested by previous studies [34]. This implies that electrical polarization along the two other directions is neglected. In addition, only transverse loading $b_{3}$ is taken into account, thus ignoring the effects of $b_{1}$ and $b_{2}$ . Moreover, the homogeneity of the shell means a constant mass density.

Regarding the classical mechanical boundary conditions of simply supported doubly curved shells, it is typically assumed that each boundary of the shell is fully constrained in the boundary plane, with the displacement components in every boundary plane being zero. However, in the out-of-plane direction of each boundary, there are no classical displacement constraints imposed [46]. While there are theoretically $2^{5} = 32$ possible selections of non-classical boundary conditions on each boundary, only one of these can yield an analytical solution to the problem at hand. Therefore, the assumed mechanical boundary conditions for simply supported flexoelectric doubly curved shells are:

{\begin{matrix} u_{2} = 0, u_{3} = 0, ϕ_{2} = 0, u_{1, 1} = 0, ϕ_{1, 1} = 0, \\ {\bar{N}}_{11} = 0, {\bar{M}}_{11} = 0, {\hat{T}}_{12} = 0, {\hat{T}}_{13} = 0, {\hat{S}}_{12} = 0, on θ_{1} = const, \end{matrix}

(55)

{\begin{matrix} u_{1} = 0, u_{3} = 0, ϕ_{1} = 0, u_{2, 2} = 0, ϕ_{2, 2} = 0, \\ {\bar{N}}_{22} = 0, {\bar{M}}_{22} = 0, {\hat{T}}_{21} = 0, {\hat{T}}_{23} = 0, {\hat{S}}_{21} = 0, on θ_{2} = const . \end{matrix}

(56)

Notably, in the case of mechanical sharp-edge conditions, each sharp edge corresponds to the intersection of two boundaries. Therefore, any mechanical sharp-edge conditions must satisfy both displacement boundary conditions simultaneously:

\begin{matrix} u_{1} = u_{2} = u_{3} = ϕ_{1} = ϕ_{2} = 0 . \end{matrix}

(57)

Furthermore, according to the variational principle for flexoelectricity, the complete electrical boundary conditions are:

{\begin{array}{l} φ (- h / 2) = 0, \\ {[P_{3} - l_{φ}^{2} P_{3, 33} - ε_{0} φ_{, 3} + ε_{0} l_{φ}^{2} φ_{, 333}]}_{z = h / 2} = 0, \\ {[V_{33} + l_{φ}^{2} φ_{, 33}]}_{z = \pm h / 2} = 0, \\ {[P_{3, 3} - ε_{0} φ_{, 33}]}_{z = \pm h / 2} = 0, \end{array}

(58)

where $l_{φ}^{2} = β_{1}^{2} + β_{2}^{2} + β_{3}^{2}$ . Consequently, the solutions for $φ$ and $P_{3}$ can be obtained as:

\begin{matrix} φ = \frac{f_{1}}{1 + a ε_{0}} {[\frac{1}{R_{1}} u_{1, 1} + \frac{1}{R_{2}} u_{2, 2} + (\frac{1}{R_{1}} + \frac{1}{R_{2}}) u_{3}] (\frac{\cosh (z / l_{φ})}{\cosh (h / 2 l_{φ})} - 1) \\ + \frac{h}{2} (\frac{1}{R_{1}} ϕ_{1, 1} + \frac{1}{R_{2}} ϕ_{2, 2}) (\frac{\sinh (z / l_{φ})}{\sinh (h / 2 l_{φ})} + 1)} \\ - \frac{z + h / 2}{1 + a ε_{0}} [- (f_{1} + 3 f_{2}) (\frac{1}{R_{1}^{2}} u_{1, 1} + \frac{1}{R_{2}^{2}} u_{2, 2}) - (f_{1} + 2 f_{2}) (\frac{1}{R_{1}^{2}} + \frac{1}{R_{2}^{2}}) u_{3} \\ + f_{2} (\frac{1}{R_{1}^{2}} u_{3, 11} + \frac{1}{R_{2}^{2}} u_{3, 22}) + (2 f_{1} + f_{2}) (\frac{1}{R_{1}} ϕ_{1, 1} + \frac{1}{R_{2}} ϕ_{2, 2})], \end{matrix}

(59)

\begin{matrix} P_{3} = \frac{f_{1} ε_{0}}{l_{φ} (1 + a ε_{0})} {[\frac{1}{R_{1}} u_{1, 1} + \frac{1}{R_{2}} u_{2, 2} + (\frac{1}{R_{1}} + \frac{1}{R_{2}}) u_{3}] \frac{\sinh (z / l_{φ})}{\cosh (h / 2 l_{φ})} \\ + \frac{h}{2} (\frac{1}{R_{1}} ϕ_{1, 1} + \frac{1}{R_{2}} ϕ_{2, 2}) \frac{\cosh (z / l_{φ})}{\sinh (h / 2 l_{φ})}} \\ - \frac{ε_{0}}{1 + a ε_{0}} [- (f_{1} + 3 f_{2}) (\frac{1}{R_{1}^{2}} u_{1, 1} + \frac{1}{R_{2}^{2}} u_{2, 2}) - (f_{1} + 2 f_{2}) (\frac{1}{R_{1}^{2}} + \frac{1}{R_{2}^{2}}) u_{3} \\ + f_{2} (\frac{1}{R_{1}^{2}} u_{3, 11} + \frac{1}{R_{2}^{2}} u_{3, 22}) + (2 f_{1} + f_{2}) (\frac{1}{R_{1}} ϕ_{1, 1} + \frac{1}{R_{2}} ϕ_{2, 2})] . \end{matrix}

(60)

By substituting the expression for $P_{3}$ and constitutive relations, the mechanical governing equations with respect to basic variables can be obtained.

3.2. Closed-form solutions

For the static bending problem, terms related to time derivatives are disregarded, resulting in the automatic vanishing of all acceleration terms in the mechanical governing equations. Moreover, the transverse loading $b_{3}$ is considered to be uniform and equal to $- q$ . The exact solutions for $u_{1}, u_{2}, u_{3}, ϕ_{1},$ and $ϕ_{2}$ can be obtained by taking into account governing equations and boundary conditions in their respective forms:

u_{1} (θ_{1}, θ_{2}) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} U_{mn} \cos (\frac{m π θ_{1}}{α_{0}}) \sin (\frac{n π θ_{2}}{β_{0}}),

(61)

u_{2} (θ_{1}, θ_{2}) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} V_{mn} \sin (\frac{m π θ_{1}}{α_{0}}) \cos (\frac{n π θ_{2}}{β_{0}}),

(62)

u_{3} (θ_{1}, θ_{2}) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} W_{mn} \sin (\frac{m π θ_{1}}{α_{0}}) \sin (\frac{n π θ_{2}}{β_{0}}),

(63)

ϕ_{1} (θ_{1}, θ_{2}) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} X_{mn} \cos (\frac{m π θ_{1}}{α_{0}}) \sin (\frac{n π θ_{2}}{β_{0}}),

(64)

ϕ_{2} (θ_{1}, θ_{2}) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} Y_{mn} \sin (\frac{m π θ_{1}}{α_{0}}) \cos (\frac{n π θ_{2}}{β_{0}}) .

(65)

It is evident that the aforementioned solutions satisfy all of the mechanical boundary conditions. Furthermore, the transverse loading $q$ can also be expanded as:

q = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} Q_{mn} \sin (\frac{m π θ_{1}}{α_{0}}) \sin (\frac{n π θ_{2}}{β_{0}}),

(66)

where coefficients $Q_{mn}$ are calculated as:

Q_{mn} = \frac{4 q}{π^{2} mn} [1 - \cos (m π)] [1 - \cos (n π)] .

(67)

Upon substituting the aforementioned Fourier expansions into the mechanical governing equations, five linear algebraic equations with respect to $U_{mn}, V_{mn}, W_{mn}, X_{mn}$ , and $Y_{mn}$ are obtained. This can be organized in matrix notation as:

[K] {Δ} + {F} = 0,

(68)

where $[K]$ is a $5 \times 5$ symmetric matrix, and ${Δ},$ and ${F}$ are expressed as:

{Δ}^{T} = {U_{mn}, V_{mn}, W_{mn}, X_{mn}, Y_{mn}},

(69)

{F}^{T} = {0, 0, - Q_{mn}, 0, 0} .

(70)

The linear algebraic equations obtained by substituting the Fourier expansions into the mechanical governing equations can be solved to obtain $U_{mn}, V_{mn}, W_{mn}, X_{mn},$ and $Y_{mn}$ . These solutions can then be utilized to obtain exact solutions to $u_{1}, u_{2}, u_{3}, ϕ_{1},$ and $ϕ_{2}$ . Substitution of the analytical solutions for these displacement components into the expressions for $φ$ and $P_{3}$ yields the exact solutions for the electrical variables. Therefore, closed-form solutions to the static bending problem are completely obtained, based on which a systematic parameter study is made in the following section.

3.3. Case studies and discussions

This study can be validated through two different scenarios: classical homogeneous isotropic doubly curved shells and microscale doubly curved shells incorporating couple stress effects. This model is highly versatile and can be adapted to achieve two simplified circumstances by setting various parameters to zero. Specifically, when all length scale parameters and electrical variables are zero and electrical variables are zero with $l_{0} = l_{1} = 0, l_{2} = l / \sqrt{2}$ , the two simplified scenarios can be achieved. All cases considered in the study adopt simply supported boundary conditions. The reliability of this model is validated by comparing with other literature in Appendix 3. It is shown that this study has a very satisfactory agreement with previous studies.

The flexoelectric nanoshell considered in this study is assumed to be composed of barium strontium titanate ceramic $(B a_{0.5} S r_{0.5} Ti O_{3})$ . This material exhibits a paraelectric phase above its Curie temperature, which eliminates the piezoelectric effect induced by stress polarization and allows for the observation of a pure flexoelectric effect [2]. Furthermore, this ceramic material possesses $\infty \infty m$ Curie group symmetry within its volume [47], which renders it an isotropic dielectric with two independent flexoelectric coefficients. The mechanical material parameters considered in this study are obtained from Assali et al. [48], where the density is $ρ = 5640 kg / m^{3}$ , the bulk modulus is $k = 182.9 GPa$ , and the shear modulus is $μ = 118.6 GPa$ . As suggested by Shu et al. [49], the flexocoupling coefficients $f_{1}$ and $f_{2}$ of BST ceramics have approximately equal absolute values but opposite signs. The electrical material parameters adopted in this study are taken from Hou et al. [50], where the permittivity of vacuum is $ε_{0} = 8.854 \times 10^{- 12} C / (V \cdot m)$ , the reciprocal dielectric susceptibility is $a = 6.6 \times 10^{7} (V \cdot m) / C$ , and the flexocoupling coefficients are $f_{1} = - f_{2} = 66 V$ . In addition, the material length parameters of strain gradient are the same, such that $l_{0} = l_{1} = l_{2} = l$ .

Regarding the structural geometry parameters, as previously mentioned, the values of $R_{1}$ and $R_{2}$ can be adjusted to simulate various structural forms. Specifically, by setting $1 / R_{1} = 1 / R_{2} = 0$ , $1 / R_{1} = 0$ , $R_{2} = R$ , and $R_{1} = R_{2} = R$ , the doubly curved shell model can be specialized to plates and cylindrical and spherical shells, respectively. Unless otherwise stated, the structural geometry parameters are: $h = 20 nm$ , $L_{α} = R_{1} α_{0} = L_{β} = R_{2} β_{0} = 25 h$ , $h / l = 2$ , and $h / l_{φ} = 10$ . The Winkler–Pasternak foundation parameters, $k_{1}$ and $k_{2}$ , are not considered at this time. The shear correction factor is $k_{s} = 5 / 6$ , and the transverse uniform loading is $q = 1 MPa$ .

The shell’s stiffness and the upper-surface electrical potential are two primary parameters of interest for the bending problem. Analysis of Figure 2 reveals that the ratio of the two principal curvature radii, $R_{2} / R_{1}$ , is crucial in determining the direct flexoelectric effect. Specifically, for a fixed $R_{2}$ , an increase in $R_{1}$ leads to more significant transverse deformation, indicating lower structural stiffness. Furthermore, the electrical potential on the upper surface is zero at simply supported edges, but it is higher in the middle region of the shell surface. However, increased transverse deformation results in a greater electrical polarization extent. Figure 3 illustrates the effect of the absolute value of the radius of curvature $R$ on the direct flexoelectric effect, with two principal curvature radii representing spherical flexoelectric shells. As the radius of curvature of spherical shells increases, the structural deflection and the magnitude of upper-surface electrical potential also increase. Consequently, we can conclude that, with all other parameters being equal, flexoelectric plates exhibit the least stiffness, followed by flexoelectric cylindrical shells. In contrast, flexoelectric spherical shells possess significant stiffness.

As depicted in Figure 4, it is observed that a general trend exists where the stiffness of the shells decreases as the material length parameter of the strain gradient decreases. However, when $l$ is sufficiently small, the maximum transverse deformation will decline slightly, and the deflection in the central area of the shell will become uniform. Notably, as $l$ decreases, the maximum absolute electrical potential will increase, leading to more significant fluctuations of electrical potential along the shell length direction. When $h / l$ reaches 10, the middle upper-surface potential approaches zero. However, Figure 5 shows that the structural stiffness and electrical polarization extent of flexoelectric shells are not highly sensitive to the length scale parameter of electrical polarization $l_{φ}$ .

Furthermore, the influence of the dimensions of the shells’ two transverse directions on the direct flexoelectric effect cannot be disregarded. As shown in Figure 6, it is observed that for spherical flexoelectric shells, larger side lengths $L$ lead to more significant transverse deformation. However, when $L$ is sufficiently large, the maximum transverse deflection will decrease slightly, with the central deflection being slightly smaller than that of the surrounding area. It is noteworthy that when the transverse size of the spherical shell greatly exceeds its thickness, the electrical potential distribution on the upper surface becomes quite non-uniform, with a smaller maximum absolute value of the electrical potential. However, Figure 7 demonstrates that more similar sizes of the two transverse dimensions of flexoelectric doubly curved shells lead to greater structural stiffness. However, the electrical polarization extent of such shells is not significantly affected by the ratio $L_{α} / L_{β}$ .

Moreover, it is necessary to consider the Winkler–Pasternak foundation’s impact on the direct flexoelectric effect. Because energy harvesters and sensors are typically supported by an elastic medium. As demonstrated in Figure 8, the first parameter of the Winkler–Pasternak foundation, $k_{1}$ , results in less transverse deflection and electrical polarization extent of flexoelectric shells. However, larger values of $k_{1}$ simultaneously generate more significant fluctuations of upper-surface electrical potential. Similarly, Figure 9 illustrates that greater values of $k_{2}$ , representing the second parameter of the Winkler-Pasternak foundation, contribute to smaller transverse deflection and upper-surface electrical potential, similar to the effect of $k_{1}$ . Notably, a key difference between $k_{1}$ and $k_{2}$ is that higher values of $k_{2}$ result in a flatter distribution of upper-surface electrical potential, which contrasts with the role of $k_{1}$ in this regard.

The theoretical findings discussed above have significant implications in designing flexoelectric energy harvesting devices and sensors. It is advisable to maintain equal dimensions in the two transverse directions $(L_{α} = L_{β})$ , as this would result in a more substantial structural rigidity without significantly compromising the upper-surface electrical potential. Moreover, it is advisable to avoid much larger lateral dimensions of flexoelectric shells than their thickness, as this would lead to significant transverse deflection and greater fluctuations of the upper-surface electrical potential, resulting in decreased device efficiency. In addition, the thickness of the flexoelectric shells should not be excessively large, and it is optimal to take the shell thickness twice the length scale parameter. Finally, the Winkler–Pasternak foundation’s parameters must be suitably adjusted to ensure moderate transverse deflection and a relatively large upper-surface electrical potential.

4. Free vibration analysis

4.1. Closed-form solutions

As for the dynamic problem, the transverse loading $b_{3}$ is absent. In this scenario, the functions $u_{1}, u_{2}, u_{3}, ϕ_{1},$ and $ϕ_{2}$ can also be expanded into Fourier series, which satisfy all the mechanical boundary conditions:

u_{1} (θ_{1}, θ_{2}, t) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} U_{mn}^{D} \cos (\frac{m π θ_{1}}{α_{0}}) \sin (\frac{n π θ_{2}}{β_{0}}) e^{i ω t},

(71)

u_{2} (θ_{1}, θ_{2}, t) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} V_{mn}^{D} \sin (\frac{m π θ_{1}}{α_{0}}) \cos (\frac{n π θ_{2}}{β_{0}}) e^{i ω t},

(72)

u_{3} (θ_{1}, θ_{2}, t) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} W_{mn}^{D} \sin (\frac{m π θ_{1}}{α_{0}}) \sin (\frac{n π θ_{2}}{β_{0}}) e^{i ω t},

(73)

ϕ_{1} (θ_{1}, θ_{2}, t) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} X_{mn}^{D} \cos (\frac{m π θ_{1}}{α_{0}}) \sin (\frac{n π θ_{2}}{β_{0}}) e^{i ω t},

(74)

ϕ_{2} (θ_{1}, θ_{2}, t) = \sum_{m = 1}^{\infty} \sum_{n = 1}^{\infty} Y_{mn}^{D} \sin (\frac{m π θ_{1}}{α_{0}}) \cos (\frac{n π θ_{2}}{β_{0}}) e^{i ω t},

(75)

where $ω$ denotes the frequency of natural vibration. Substituting the expressions above into the mechanical governing equations will yield the following linear algebraic equations:

[K] {Δ} + ω^{2} [M] {Δ} = 0,

(76)

in which the coefficients $K_{ij}$ are identical to those of the static bending problem. The nonzero components of $[M]$ are:

M_{11} = M_{22} = M_{33} = ρ h, M_{44} = M_{55} = \frac{ρ h^{3}}{12} .

(77)

An untrivial solution ${Δ}$ requires the following relation:

det [[K] + ω^{2} [M]] = 0,

(78)

from which the frequency of natural vibration $ω$ can be obtained. Therefore, we have obtained the exact solutions to the free vibration problem, which is the basis of the following discussions.

4.2. Case studies and discussions

Based on the results presented in Figure 10, it is seen that the first natural frequency of flexoelectric doubly curved shells is dependent on the curvature radii. Specifically, as the curvature radii gradually increase, the first natural frequency experiences an initial rapid decline, followed by a slower decline until it reaches a constant value. Furthermore, the ratio of the two curvature radii also influences the first natural frequency, with a smaller difference between the two radii resulting in a higher first natural frequency, indicating a greater degree of structural stiffness. Therefore, holding other parameters constant, the spherical flexoelectric shell exhibits the greatest stiffness, consistent with the findings from the bending analysis conducted earlier.

Figure 11 demonstrates the impact of material length parameters on the first natural frequency. A larger ratio of thickness to material length parameter leads to a decrease in the first natural frequency. This behavior is primarily evident in the range of $h / l < 10$ , beyond which the decline becomes less pronounced. As a consequence, the first-order natural frequency of flexoelectric nanoshells is found to be subject to a significant size-dependent effect. Conversely, the variation of the length scale parameter of electrical polarization $l_{φ}$ has only a negligible impact on the first natural frequency. Thus, it can be inferred that the first-order natural frequency is insensitive to $l_{φ}$ .

From Figure 12, it is observed that the transverse sizes of spherical flexoelectric shells also impact their first natural frequency. Larger lateral sizes result in lower natural vibration frequencies, albeit the decline is gradual without a significant drop. In addition, when the two lateral structural sizes are close to each other, the first natural frequency of flexoelectric spherical shells also increases. However, when $L_{α} / L_{β} \geq 2$ , the increase in $L_{α}$ no longer triggers a significant difference in the first natural frequency. Furthermore, Figure 13 demonstrates that increasing both the Winkler–Pasternak foundation parameters results in a higher first natural frequency, thus increasing the flexoelectric shells’ stiffness. Therefore, in the actual design of flexoelectric energy harvesters and sensors, it is crucial to adjust the flexoelectric doubly curved shell parameters to ensure that its first natural frequency is higher than the loading frequency to avoid resonance.

5. Conclusion

This study introduces a novel analytical model for nanoscale flexoelectric doubly curved shells, which for the first time incorporates the reformulated strain gradient elasticity theory, the modified expression for electric enthalpy density, and the Reissner–Mindlin shell theory. The complete governing equations, boundary conditions, and sharp-edge conditions are obtained using Hamilton’s principle. This highly general model can be reduced to numerous specified models and exhibits size-dependent effects.

To further illustrate this new model, we establish a complete formulation of the static bending and free vibration problems of flexoelectric doubly curved shells and provide exact solutions. The impact of principal curvature radii, length scale parameters, transverse structural sizes, and the Winkler–Pasternak foundation parameters on the direct flexoelectric effect is systematically investigated. Our results indicate that the structural stiffness of flexoelectric cylindrical shells is greater than that of flexoelectric plates and less than that of flexoelectric spherical shells. In addition, larger length scale parameters of strain gradient contribute to greater structural stiffness and a higher first natural frequency. However, these properties are insensitive to the length scale parameter of electrical polarization. Smaller transverse structural sizes of flexoelectric doubly curved shells will lead to greater stiffness and higher first natural frequency while simultaneously enhancing the magnitude of electrical polarization. The Winkler–Pasternak foundation will enhance the bending stiffness and first natural frequency of flexoelectric doubly curved shells.

This highly general flexoelectric doubly curved shell model has important implications for the design of flexoelectric energy harvesters and sensors based on the direct flexoelectric effect, as well as serving as a benchmark for numerical simulations and experimental validations.

Footnotes

Appendix 1

Appendix 2

Appendix 3 Acknowledgements

The authors thank Professor David Steigmann and two anonymous reviewers for their encouragement and helpful comments.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Jinchen Xie

References

Krichen

Sharma

Flexoelectricity: a perspective on an unusual electromechanical coupling. J Appl Mech 2016; 83(3): 030801.

Yudin

Tagantsev

AK.

Fundamentals of flexoelectricity in solids. Nanotechnology 2013; 24(43): 432001.

Nguyen

Mao

Yeh

, et al. Nanoscale flexoelectricity. Adv Mater 2013; 25(7): 946–974.

Wang

Zhang

, et al. Flexoelectricity in solids: progress, challenges, and perspectives. Prog Mater Sci 2019; 106: 100570.

Enakoutsa

Corte

Giorgio

A model for elastic flexoelectric materials including strain gradient effects. Math Mech Solids 2016; 21(2): 242–254.

Enakoutsa

Vescovo

Scerrato

Combined polarization field gradient and strain field gradient effects in elastic flexoelectric materials. Math Mech Solids 2017; 22(5): 938–951.

Shu

Huang

, et al. Relationship between direct and converse flexoelectric coefficients. J Appl Phys 2014; 116(14): 144105.

Liu

Deng

Xie

, et al. Electromechanical analysis of direct and converse flexoelectric effects under a scanning probe tip. J Mech Phys Solids 2020; 142: 104020.

Springolo

Royo

Stengel

Direct and converse flexoelectricity in two-dimensional materials. Phys Rev Lett 2021; 127: 216801.

10.

Mao

Xiang

, et al. Theory for dielectrics considering the direct and converse flexoelectric effects and its finite element implementation. Appl Math Model 2016; 40(15): 7115–7137.

11.

Yan

Modeling of a nanoscale flexoelectric energy harvester with surface effects. Phys E Low Dimens Syst Nanostruct 2017; 88: 125–132.

12.

Deng

Shen

Flexoelectric energy harvesting using circular thin membranes. J Appl Mech 2020; 87(9): 091004.

13.

Deng

Kammoun

Erturk

, et al. Nanoscale flexoelectric energy harvesting. Int J Solids Struct 2014; 51(18): 3218–3225.

14.

Moura

Erturk

Electroelastodynamics of flexoelectric energy conversion and harvesting in elastic dielectrics. J Appl Phys 2017; 121(6): 064110.

15.

Tzou

Distributed flexoelectric structural sensing: theory and experiment. J Sound Vib 2015; 348: 126–136.

16.

Kwon

Huang

Zhang

, et al. Flexoelectric sensing using a multilayered barium strontium titanate structure. Smart Mater Struct 2013; 22(11): 115017.

17.

Huang

Kwon

Zhang

, et al. A trapezoidal flexoelectric accelerometer. J Intell Mater Syst Struct 2014; 25(3): 271–277.

18.

Huang

Yan

Kwon

, et al. Flexoelectric strain gradient detection using Ba0.64Sr0.36TiO3 for sensing. Appl Phys Lett 2012; 101(25): 252903.

19.

Zhang

Liu

, et al. A curved resonant flexoelectric actuator. Appl Phys Lett 2017; 111(8): 082904.

20.

Liu

Zhang

, et al. An actuation method by a biconcave beam structure with converse flexoelectric effect. Smart Mater Struct 2019; 28(11): 115025.

21.

Ray

MC.

Analysis of smart nanobeams integrated with a flexoelectric nano actuator layer. Smart Mater Struct 2016; 25(5): 055011.

22.

Chen

Yan

Investigation of pull-in behaviors of a nanoswitch tuned by piezoelectric and flexoelectric effects. Int J Mech Sci 2019; 161–162: 105032.

23.

Huang

Kim

Hou

, et al. Flexoelectric characterization of BaTiO3-0.08Bi(Zn1/2Ti1/2)O3. Appl Phys Lett 2017; 110(22): 222904.

24.

Linder

An analysis of the exponential electric displacement saturation model in fracturing piezoelectric ceramics. Tech Mech 2012; 32: 53–69.

25.

Linder

Miehe

Effect of electric displacement saturation on the hysteretic behavior of ferroelectric ceramics and the initiation and propagation of cracks in piezoelectric ceramics. J Mech Phys Solids 2012; 60(5): 882–903.

26.

Linder

A complex variable solution based analysis of electric displacement saturation for a cracked piezoelectric material. J Appl Mech 2014; 81(9): 091006.

27.

Kargarnovin

Hashemi

, et al. Axisymmetric thermo-electro-elastic analysis of a piezoelectric half space. Math Mech Solids 2012; 17(5): 500–515.

28.

Aouadi

Aspects of uniqueness in micropolar piezoelectric bodies. Math Mech Solids 2008; 13(6): 499–512.

29.

Kozinov

Sheveleva

Loboda

Fracture behavior of periodically bonded interface of piezoelectric bi-material under compressive–shear loading. Math Mech Solids 2019; 24(10): 3216–3230.

30.

Bhargava

Jangid

A mathematical strip-saturation model for a piezoelectric plane weakened by two collinear equal cracks. Math Mech Solids 2014; 19(6): 713–724.

31.

Mao

Purohit

PK.

Insights into flexoelectric solids from strain-gradient elasticity. J Appl Mech 2014; 81(8).

32.

Zhou

, et al. A reformulated flexoelectric theory for isotropic dielectrics. J Phys D Appl Phys 2015; 48(46): 465502.

33.

Tzou

Zhang

XF.

A Flexoelectric double-curvature nonlinear shell energy harvester. J Vib Acoust 2016; 138(3): 031006.

34.

Huang

, et al. On the mechanics of curved flexoelectric microbeams. Int J Eng Sci 2018; 124: 1–15.

35.

Zhou

A flexoelectric spherical microshell model incorporating the strain gradient effect. Appl Math Model 2019; 75: 692–708.

36.

Huang

, et al. Modeling of the flexoelectric annular microplate based on strain gradient elasticity theory. Mech Adv Mater Struct 2019; 26(23): 1958–1968.

37.

Guo

Zhang

, et al. A new model for circular cylindrical Kirchhoff–Love shells incorporating microstructure and flexoelectric effects. J Appl Mech 2022; 89(12): 121010.

38.

Mindlin

RD.

Microstructure in linear elasticity. Technical report, Department of Civil Engineering and Engineering Mechanics, Columbia University, New York, 1 October 1963.

39.

Zhou

Wang

A reformulation of constitutive relations in the strain gradient elasticity theory for isotropic materials. Int J Solids Struct 2016; 80: 28–37.

40.

Zhou

SJ.

The strain gradient elasticity theory in orthogonal curvilinear coordinates and its applications. J Mech 2018; 34(3): 311–323.

41.

Mindlin

High frequency vibrations of piezoelectric crystal plates. Int J Solids Struct 1972; 8(7): 895–906.

42.

Maranganti

Sharma

Electromechanical coupling in nonpiezoelectric materials due to nanoscale nonlocal size effects: Green’s function solutions and embedded inclusions. Phys Rev B 2006; 74: 014110.

43.

Shen

Electric field gradient theory with surface effect for nano-dielectrics. Comput Mater Contin 2009; 13(1): 63–88.

44.

Lam

Yang

Chong

, et al. Experiments and theory in strain gradient elasticity. J Mech Phys Solids 2003; 51(8): 1477–1508.

45.

Zhang

Gao

, et al. A transversely isotropic magneto-electro-elastic Timoshenko beam model incorporating microstructure and foundation effects. Mech Mater 2020; 149: 103412.

46.

Reddy

JN.

Exact solutions of moderately thick laminated shells. J Eng Mech 1984; 110(5): 794–809.

47.

Cross

LE.

Flexoelectric polarization of barium strontium titanate in the paraelectric state. Appl Phys Lett 2002; 81(18): 3440–3442.

48.

Assali

Kanouni

Laidoudi

, et al. Structural and electromechanical properties of Sr-substituted barium titanate (BST) as potential material for high performance electroacoustic devices. Mater Today Commun 2020; 25: 101643.

49.

Shu

Huang

Kwon

, et al. Converse flexoelectric coefficient f1212 in bulk Ba0.67Sr0.33TiO3. Appl Phys Lett 2014; 104(23): 232902.

50.

Hou

Tian

Chu

Flexoelectric response of (1-x)BaTiO3-xSrTiO3 ceramics. Ceram Int 2020; 46(9): 12928–12932.

51.

Leissa

The free vibration of rectangular plates. J Sound Vib 1973; 31(3): 257–293.

52.

Zhang

Liu

, et al. A non-classical Mindlin plate finite element based on a modified couple stress theory. Eur J Mech A Solids 2013; 42: 63–80.