Computational Modeling of Flexoelectricity—A Review

Zhuang, Xiaoying; Nguyen, Binh Huy; Nanthakumar, Subbiah Srivilliputtur; Tran, Thai Quoc; Alajlan, Naif; Rabczuk, Timon

doi:10.3390/en13061326

开放式访问审查

挠曲电的计算模型——综述

¹

越南胡志明市Ton Duc Thang大学计算力学系，邮编758307

²

越南胡志明市东德成大学土木工程学院758307

^三

阿佩尔斯特汉诺威莱布尼茨大学连续介质力学研究所。11，30167德国汉诺威

⁴

沙特阿拉伯利雅得国王沙特大学计算机与信息科学学院计算机工程系，邮编：11543

^*

信件应寄给的作者。

能源 2020,13(6), 1326;https://doi.org/10.3390/en13061326

收到的提交文件：2020年1月21日/修订日期：2020年2月23日/接受时间：2020年2月26日/发布日期：2020年3月12日

（本文属于特刊多物理问题的计算方法Ⅱ)

下载

浏览地物

版本注释

摘要

:

机电耦合器件在现代工程中发挥着不可或缺的作用。尤其是柔性电，一种涉及应变梯度的机电耦合效应，在未来的小型化机电耦合器件中显示出巨大的潜力。因此，对柔性电进行仿真是必要的和不可避免的。在这篇文章中，我们概述了模拟柔性电的数值程序。具体来说，我们总结了一个包括静电应力张量的广义公式，可以对其进行简化以检索文献中的其他公式。我们进一步说明了不同数值方法（包括等几何分析和混合有限元法）的边值问题的弱形式和离散形式。给出了几个基准问题来演示数值实现。实现的源代码可用于分析和开发更复杂的柔性电纳米器件。

关键词：

挠曲电;数值方法;建模

1.简介

近年来，柔性电学在材料科学、机械工程和化学等各个研究领域受到了极大的关注[1,2]，所有电介质材料中的机电耦合现象（无论材料对称性如何）。它描述了应变梯度引起的激发极化。尽管与压电效应相比，其系数较小，但由于微型化趋势，柔性电效应得到了更多的关注，并且非常显著，许多机电系统在亚微米或纳米尺度上运行，导致应变梯度大幅增加。由于应变梯度的尺寸效应，柔性电在传感器等微/纳米机电系统中具有广阔的应用前景[三,4,5]，执行器[6,7]和纳米能源收割机[8,9,10,11,12,13]. 它还可以解释诸如应变梯度驱动的极性控制等物理现象[14,15,16]，柔性光伏效应[17]和弯曲热量效应[18]. 感兴趣的读者可以参考其他有关柔性电的优秀评论文章，了解更多详细信息[19,20,21,22].

据我们所知，Maranganti等人首次提出了介质固体中挠曲电的连续介质理论[23]他以应变、应变梯度、极化和极化梯度为自变量，考虑了内能中的直接和反向柔性电效应。正如他们在工作中指出的那样，该公式可以被视为对Mindlin提出的极化粒度理论的补充[24]随后由Sahin和Host主持[25]. 基于他们的公式，Sharma的小组进一步研究了纳米结构中挠曲电的增强[26,27,28,29,30]. 应该指出，在上述公式中，没有考虑麦克斯韦应力发散的电力。Hu和Shen指出，在纳米电介质固体中，考虑麦克斯韦应力可能很重要[31]，他利用了Kuang提出的物理变化原理[32,33,34,35]使静电应力自然地从电焓的最小化中出现。在Liu的工作中，在连续磁电弹性的统一能量公式中也讨论了麦克斯韦应力的重要性[36,37].

然而，由于非定域性，即应变梯度或极化梯度，固体电介质中的挠曲电由四阶偏微分方程控制。因此，挠曲电的建模大多简化为一维模型或具有轴对称（例如圆柱体或圆盘）或板的简单几何问题。与应变梯度弹性力学一样，柔性电的全多维建模的困难在于

{C类}^{1}

位移场近似的连续性。这一问题首先由Abdollahi等人解决[38]使用局部最大熵（LME）无网格方法。Mao等人还采用了类似于应变梯度弹性的混合有限元法来解决柔性电问题[39,40]. 等几何分析（IGA）也是柔性电学的一种稳健方法，如[41,42,43,44]这得益于它对高阶连续性的直接处理。或者，Argyris三角形元素也可以用于补救

{C类}^{1}

要求，如Yvonnet和Liu的工作[45]. 最近，本着与IGA方法类似的精神，Codony等人[46]针对柔性电问题，提出了基于分层B样条的浸没边界法，具有较高的效率。此外，Roy和Roy还用周动力学方法探讨了挠曲电效应及其相关的损伤行为[47].

随着柔性电学研究领域的迅速发展，对数值模拟的需求显然是不可避免的。然而，目前还没有对挠曲电的计算方面进行综述。因此，在本文中，我们旨在提供一个关于柔性电建模的概述和详细过程。采用等几何分析、混合有限元和无网格等不同的数值方法求解边值问题。本文的其余部分结构如下。在第2节基于介电固体中麦克斯韦应力张量的存在，我们总结了物理变分原理，在该原理中，静电应力张量自然地表现为胡和沈工作中电焓的最小化[31]. 此外，我们提出了另一种公式，该公式由电吉布斯自由能导出（但具有不同的自变量），它也会产生静电应力张量。值得注意的是，这两个公式可以被视为一个广义公式，对于该公式，当省略静电应力时，可以恢复以前文献中的强形式。因此，对于不同的数值方法，我们从Galerkin方法得到了弱形式。有关离散化形式的数值实现的详细信息，请参阅第3节尽可能与源代码一起使用。下面给出了几个数值例子来验证和说明数值实现第4节.

2.地面工程制定

2.1. 物理变分原理

胡和沈[31]从Toupin推广变分原理[48]根据况提出的物理变分原理[32,33,34,35]推导了考虑麦克斯韦应力的挠曲电控制方程。我们在本节中简要总结了理论和公式。假设内能是应变的函数

ϵ_{我 j个}

，应变梯度

ϵ_{我 j个, k个}

，极化

{P（P）}_{我}

和极化梯度

{P（P）}_{我, j个}

,

\begin{matrix} U型 (ϵ_{我 j个}, ϵ_{我 j个, k个}, {P（P）}_{我}, {P（P）}_{我, j个}) = \frac{1}{2} 一_{我 j个} {P（P）}_{我} {P（P）}_{j个} + \frac{1}{2} {b条}_{我 j个 k个 我} {P（P）}_{我, j个} {P（P）}_{k个, 我} + \frac{1}{2} {c（c）}_{我 j个 k个 我} ε_{我 j个} ε_{k个 我} + {d日}_{我 j个 k个} {P（P）}_{我} ε_{j个 k个} + {（f）}_{我 j个 k个 我} ε_{我 j个} {P（P）}_{k个, 我} + {小时}_{我 j个 k个 我} {P（P）}_{我} ε_{j个 k个, 我} + \frac{1}{2} 克_{我 j个 k个 我 米 n个} ε_{我 j个, k个} ε_{我 米, n个}, \end{matrix}

(1)

从中获得本构关系

\begin{matrix} σ_{我 j个} & = \frac{\partial U型}{{\partial ε}_{我 j个}} = {c（c）}_{我 j个 k个 我} ε_{k个 我} + {d日}_{我 j个 k个} {P（P）}_{k个} + {（f）}_{我 j个 k个 我} {P（P）}_{k个, 我} \end{matrix}

（2a）

\begin{matrix} τ_{我 j个 k个} & = \frac{\partial U型}{{\partial ε}_{我 j个, k个}} = {小时}_{我 j个 k个 我} {P（P）}_{我} + 克_{我 j个 k个 我 米 n个} ε_{我 米, n个} \end{matrix}

（2b）

\begin{matrix} {E类}_{我} & = \frac{\partial U型}{{\partial P（P）}_{我}} = 一_{我 j个} {P（P）}_{j个} + {d日}_{我 j个 k个} ε_{j个 k个} + {小时}_{我 j个 k个 我} ε_{j个 k个, 我} \end{matrix}

（2c）

\begin{matrix} {V（V）}_{我 j个} & = \frac{\partial U型}{{\partial P（P）}_{我, j个}} = {b条}_{我 j个 k个 我} {P（P）}_{k个, 我} + {（f）}_{我 j个 k个 我} ε_{k个 我}, \end{matrix}

（2天）

哪里

σ_{我 j个}

,

{E类}_{我}

,

τ_{我 j个 k个}

和

{V（V）}_{我 j个}

分别是（局部）应力、电场、高阶应力和高阶电场，而

一_{我 j个}

,

{b条}_{我 j个 k个 我}

,

{c（c）}_{我 j个 k个 我}

,

{d日}_{我 j个 k个}

,

{（f）}_{我 j个 k个 我}

、和

{小时}_{我 j个 k个 我}

是材料张量[32]. 换句话说，内能的变化可以写成

δ U型 = σ_{我 j个} δ ε_{我 j个} + τ_{我 j个 k个} δ ε_{我 j个, k个} + {E类}_{我} δ {P（P）}_{我} + {V（V）}_{我 j个} δ {P（P）}_{我, j个} .

（3）

能量密度是内能和麦克斯韦自场内部的总和[24]

\begin{matrix} W公司 = U型 + \frac{1}{2} ε_{0} ϕ_{, 我} ϕ_{, 我}, \end{matrix}

(4)

其中麦克斯韦电场可以定义为电势的梯度

ϕ

{E类}_{我}^{M（M） S公司} = - ϕ_{, 我} .

（5）

接下来，我们可以将电焓定义为

H（H） = W公司 - {E类}_{我} {D类}_{我} = U型 - \frac{1}{2} ε_{0} ϕ_{, 我} ϕ_{, 我} + ϕ_{, 我} {P（P）}_{我}

（6）

我们考虑介质固体占据一个体积v（v）以曲面为界一与周围环境分离

{v（v）}^{'}

（在本研究中被视为空气）。拿

{v（v）}^{*} = v（v） + {v（v）}^{'}

，该系统在等温过程中的平衡可以用变分原理定义

δ {¦Β}_{{v（v）}^{*}} H（H） {d日 v（v）}^{*} = δ W公司,

(7)

哪里W公司外功是由体力做的吗

{（f）}_{我}

，外部电场

{E类}_{我}^{提取}

和免费

ρ^{（f）}

论固体与牵引

{\bar{t吨}}_{我}

，双牵引

{\bar{τ}}_{我}

，表面电荷

σ^{*}

和高阶电场

{\bar{v（v）}}_{我}

在表面边界上。因此，它的变化

δ W公司

定义为

δ W公司 = {¦Β}_{v（v）} ({（f）}_{我} δ {u个}_{我} + {E类}_{我}^{提取} δ {P（P）}_{我}) d日 v（v） + {¦Β}_{一^{σ}} {t吨}_{我} δ {u个}_{我} d日 一 + {¦Β}_{一^{τ}} {\bar{τ}}_{我} {\tilde{D类}}^{n个} (δ {u个}_{我}) d日 一 + {¦Β}_{一^{V（V）}} {\bar{v（v）}}_{我} δ {P（P）}_{我} d日 一 - {¦Β}_{一^{ϕ}} σ^{*} δ ϕ d日 一 - {¦Β}_{v（v）} ρ^{（f）} δ ϕ d日 v（v）,

(8)

其中操作员

{\tilde{D类}}^{n个}

将在下一节中进行说明。

在物理变分原理方法中，电焓的变化可以进一步表示为[32,33,34,35].

δ {¦Β}_{{v（v）}^{*}} H（H） {d日 v（v）}^{*} = {¦Β}_{{v（v）}^{*}} δ H（H） {d日 v（v）}^{*} + {¦Β}_{{v（v）}^{*}} {H（H）}^{e（电子）} δ {u个}_{k个, k个} {d日 v（v）}^{*},

(9)

哪里

{H（H）}^{e（电子）} = \frac{1}{2} ({E类}_{我} {P（P）}_{我} + {V（V）}_{我 j个} {P（P）}_{我, j个})

是从总能量减去纯变形能量后减去的能量。此外，迁移变化是物理变分原理的基本特征，即虚位移不仅产生应变变化，还引起电位和极化变化，例如[31]

\begin{matrix} δ ϕ & = δ_{ϕ} ϕ + δ_{u个} ϕ = δ_{ϕ} ϕ + ϕ_{, 我} δ {u个}_{我}, \end{matrix}

（10年）

\begin{matrix} δ ϕ_{, 我} & = δ_{ϕ} ϕ_{, 我} + δ_{u个} ϕ_{, 我} = δ_{ϕ} ϕ_{, 我} + ϕ_{, 我 j个} δ {u个}_{j个}, \end{matrix}

（10亿）

\begin{matrix} δ {P（P）}_{我} & = δ_{P（P）} {P（P）}_{我} + δ_{u个} {P（P）}_{我} = δ_{P（P）} {P（P）}_{我} + {P（P）}_{我, j个} δ {u个}_{j个}, \end{matrix}

（10美分）

\begin{matrix} δ {P（P）}_{我, j个} & = δ_{P（P）} {P（P）}_{我, j个} + δ_{u个} {P（P）}_{我, j个} = δ_{P（P）} {P（P）}_{我, j个} + {P（P）}_{我, j个 k个} δ {u个}_{k个} = δ_{P（P）} {P（P）}_{我, j个} + {P（P）}_{我, k个 j个} δ {u个}_{k个}, \end{matrix}

（10天）

哪里

δ_{ϕ} (\circ)

和

δ_{P（P）} (\circ)

分别是电位和极化的（局部）变化产生的变化，而

δ_{u个} (\circ)

是由虚拟位移引起的迁移变化。随后，将方程式（10）代入方程式(9)使用方程式(6)和(三)利用高斯发散定理，通过一个冗长但直接的操作，电焓的变化可以重新表示为[31]

\begin{matrix} δ {¦Β}_{v（v）} H（H） d日 v（v） & = {¦Β}_{v（v）} [- σ_{我 j个, j个} + τ_{我 j个 米, 米 j个} - σ_{我 j个, j个}^{E类 S公司}] δ {u个}_{我} d日 v（v） \\ + {¦Β}_{一} {[σ_{我 j个} - τ_{我 j个 米, 米} + σ_{我 j个}^{E类 S公司}] {n个}_{j个} + [{\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{米} {n个}_{j个} τ_{我 j个 米} - {\tilde{D类}}_{j个}^{t吨} ({n个}_{米} τ_{我 j个 米})]} δ {u个}_{我} d日 一 \\ + {¦Β}_{一} τ_{我 j个 米} {n个}_{米} {n个}_{j个} {\tilde{D类}}^{n个} (δ {u个}_{我}) d日 一 \\ + {¦Β}_{v（v）} ({E类}_{我} - {V（V）}_{我 j个, j个} + ϕ_{, 我} + {E类}_{我}^{提取}) δ {P（P）}_{我} d日 v（v） + {¦Β}_{一} {V（V）}_{我 j个} {n个}_{j个} δ {P（P）}_{我} d日 一 \\ - {¦Β}_{v（v）} {(- ε_{0} ϕ_{, 我} + {P（P）}_{我})}_{, 我} δ ϕ d日 v（v） + {¦Β}_{一} (- ε_{0} ϕ_{, 我} + {P（P）}_{我}) {n个}_{我} δ ϕ d日 一, \end{matrix}

(11)

哪里

{n个}_{我}

是法向量。

{\tilde{D类}}_{k个}^{t吨} (\circ) = (δ_{k个 我} - {n个}_{k个} {n个}_{我}) \partial_{, 我} (\circ)

和

{\tilde{D类}}^{n个} = {n个}_{我} \partial_{, 我} (\circ)

分别是从位移梯度变化的处理中产生的切向和法向曲面微分算子

δ {u个}_{我, j个}

在边界上（另请参见[31,49]). 值得注意的是，在上述方程中，作为变分过程的直接结果，广义静电应力

σ_{我 j个}^{E类 S公司}

显示为

\begin{matrix} σ_{我 j个}^{E类 S公司} & = - {V（V）}_{k个 j个} {P（P）}_{k个, 我} + \frac{1}{2} ({E类}_{k个} {P（P）}_{k个} + {V（V）}_{k个 我} {P（P）}_{k个, 我}) δ_{我 j个} - {D类}_{j个} ϕ_{, 我} + (- \frac{1}{2} ε_{0} ϕ_{, k个} ϕ_{, k个} + ϕ_{, k个} {P（P）}_{k个}) δ_{我 j个} \\ = σ_{我 j个}^{M（M）} + τ_{我 j个 米, 米}^{M（M）}, \end{matrix}

(12)

具有

σ_{我 j个}^{M（M）}

和

τ_{我 j个 米, 米}^{M（M）}

对应于应变梯度和极化梯度的麦克斯韦应力和静电应力分别定义为

\begin{matrix} σ_{我 j个}^{M（M）} & = \frac{1}{2} {E类}_{k个} {P（P）}_{k个} δ_{我 j个} - {D类}_{j个} ϕ_{, j个} + ϕ_{, 我} + (- \frac{1}{2} ε_{0} ϕ_{, k个} ϕ_{, k个} + ϕ_{, k个} {P（P）}_{k个}) δ_{我 j个} \end{matrix}

（13年a）

\begin{matrix} τ_{我 j个 米, 米}^{M（M）} & = - {V（V）}_{k个 j个} {P（P）}_{k个, j个} + \frac{1}{2} {V（V）}_{k个 我} {P（P）}_{k个, 我} . \end{matrix}

（13亿）

代入方程式后(11)和(8)到方程式中(7)并利用

δ {u个}_{我}

,

δ ϕ

和

δ {P（P）}_{我}

注意到周围空气环境中不存在机械位移和极化，我们可以得到控制方程和边界条件

\begin{matrix} {(σ_{我 j个} - τ_{我 j个 米, 米} + σ_{我 j个}^{E类 S公司})}_{, j个} + {（f）}_{我} & = 0, 在里面 v（v） \end{matrix}

（14a）

\begin{matrix} σ_{我 j个, j个}^{E类 S公司} & = 0, 在里面 {v（v）}^{'} \end{matrix}

（14b）

\begin{matrix} {V（V）}_{我 j个, j个} + {E类}_{我}^{L（左）} - ϕ_{, 我} + {E类}_{我}^{e（电子） x个 t吨} & = 0, 在里面 v（v） \end{matrix}

（14美分）

\begin{matrix} - ε_{0} ϕ_{, 我 我} + {P（P）}_{我, 我} & = ρ^{（f）}, 在里面 v（v） \end{matrix}

（14天）

\begin{matrix} ϕ_{, 我 我} & = 0, 在里面 {v（v）}^{'} \end{matrix}

（14e）

Dirichlet边界条件

\begin{matrix} {u个}_{我} & = {\bar{u个}}_{我}, 在 一^{u个} \end{matrix}

（15年）

\begin{matrix} ϕ & = \bar{ϕ}, 在 一^{ϕ} \end{matrix}

（15亿）

\begin{matrix} {\tilde{D类}}^{n个} ({u个}_{我}) = {u个}_{我, 我} {n个}_{我} & = {\bar{u个}}_{我, 我} {n个}_{我}, 在 一^{u个, n个} \end{matrix}

（15美分）

\begin{matrix} {P（P）}_{我} & = {\bar{P（P）}}_{我}, 在 一^{P（P）} \end{matrix}

（15天）

Neumann边界条件

\begin{matrix} (σ_{我 j个} - τ_{我 j个 米, 米} + [[σ_{我 j个}^{E类 S公司}]]) {n个}_{j个} + {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{米} {n个}_{j个} τ_{我 j个 米} - {\tilde{D类}}_{j个}^{t吨} ({n个}_{米} τ_{我 j个 米}) & = {t吨}_{我}, 在 一^{σ} \end{matrix}

（16年a）

\begin{matrix} τ_{我 j个 米} {n个}_{米} {n个}_{j个} & = {\bar{τ}}_{我}, 在 一^{τ} \end{matrix}

（16亿）

\begin{matrix} {n个}_{我} (- ε_{0} [[ϕ_{, 我}]] + {P（P）}_{我}) & = - σ^{*}, 在 一^{σ^{*}} \end{matrix}

（16c）

\begin{matrix} {V（V）}_{我 j个} {n个}_{j个} & = {\bar{v（v）}}_{我}, 在 一^{V（V）} \end{matrix}

（16天）

2.2. 替代配方

上述边值问题由总能量密度表示U型谁的应变

ε_{我 j个}

，应变梯度

ε_{我 j个, k个}

，极化

{P（P）}_{我}

和极化梯度

{P（P）}_{我, j个}

是自变量。或者，我们也可以选择电吉布斯自由能（与等温系统中的电焓相同），其自变量为应变

ε_{我 j个}

，应变梯度

ε_{我 j个, k个}

，电场

{E类}_{我}

和电场梯度

{E类}_{我, j个}

[38,40,50]:

G公司 (ε_{我 j个}, ε_{我 j个, k个}, {E类}_{我}, {E类}_{我, j个}) = - \frac{1}{2} κ_{我 j个} {E类}_{我} {E类}_{j个} - \frac{1}{2} {b条}_{我 j个 k个 我} {E类}_{我, j个} {E类}_{k个, 我} + \frac{1}{2} {c（c）}_{我 j个 k个 我} ε_{我 j个} ε_{k个 我} - {e（电子）}_{我 j个 k个} {E类}_{我} ε_{j个 k个} - \frac{1}{2} μ_{我 j个 k个 我} ({E类}_{k个} ε_{我 j个, 我} - ε_{我 j个} {E类}_{k个, 我})

(17)

因此，本构关系如下所示

\begin{matrix} σ_{我 j个} & = \frac{\partial G公司}{{\partial ε}_{我 j个}} = {c（c）}_{我 j个 k个 我} ε_{k个 我} + \frac{1}{2} μ_{我 j个 k个 我} {E类}_{k个, 我} \end{matrix}

（18年）

\begin{matrix} τ_{我 j个 k个} & = \frac{\partial G公司}{{\partial ε}_{我 j个, k个}} = - \frac{1}{2} μ_{我 j个 k个 我} {E类}_{我} \end{matrix}

（18亿）

\begin{matrix} {D类}_{我} & = - \frac{\partial G公司}{{\partial E类}_{我}} = κ_{我 j个} {E类}_{j个} + \frac{1}{2} μ_{我 j个 k个 我} ε_{k个, 我} \end{matrix}

（18美分）

\begin{matrix} 问_{我 j个} & = \frac{\partial G公司}{{\partial E类}_{我, j个}} = {b条}_{我 j个 k个 我} {E类}_{k个, 我} - \frac{1}{2} μ_{我 j个 k个 我} ε_{k个 我}, \end{matrix}

（18天）

哪里

σ_{我 j个}

,

{D类}_{我}

,

τ_{我 j个 k个}

和

{V（V）}_{我 j个}

分别是（局部）应力、电位移、高阶应力和高阶电位移，而

κ_{我 j个}

,

{b条}_{我 j个 k个 我}

,

{c（c）}_{我 j个 k个 我}

,

{e（电子）}_{我 j个 k个}

和

μ_{我 j个 k个 我}

是材料张量[38]. 这些关系还意味着电吉布斯自由能的变化可以表示为

δ G公司 = σ_{我 j个} δ ε_{我 j个} + τ_{我 j个 k个} δ ε_{我 j个, k个} - {D类}_{我} δ {E类}_{我} + 问_{我 j个} δ {E类}_{我, j个},

(19)

然后可以用它从变分原理中确定平衡

δ {¦Β}_{{v（v）}^{*}} G公司 {d日 v（v）}^{*} = δ W公司,

(20)

哪里W公司同样是外功，但不像前面的公式那样有外电场和高阶电场的贡献，高阶表面电荷

\bar{q个}

在边界上规定，因此其变化

δ W公司

读取

δ W公司 = {¦Β}_{v（v）} ({（f）}_{我} δ {u个}_{我} + {E类}_{我}^{提取} δ {P（P）}_{我}) d日 v（v） + {¦Β}_{一^{σ}} {t吨}_{我} δ {u个}_{我} d日 一 + {¦Β}_{一^{τ}} {\bar{τ}}_{我} {\tilde{D类}}^{n个} (δ {u个}_{我}) d日 一 + {¦Β}_{一^{问}} \bar{q个} {\tilde{D类}}^{n个} (δ ϕ) d日 一 - {¦Β}_{一^{ϕ}} σ^{*} δ ϕ d日 一 - {¦Β}_{v（v）} ρ^{（f）} δ ϕ d日 v（v） .

(21)

类比方程（10），电场和电场梯度的迁移变化如下

\begin{matrix} δ {E类}_{我} & = δ_{ϕ} {E类}_{我} + δ_{u个} {E类}_{我} = - δ_{ϕ} ϕ_{, 我} + {E类}_{我, 对} δ {u个}_{对} = - δ_{ϕ} ϕ_{, 我} + {E类}_{对, 我} δ {u个}_{对}, \end{matrix}

（22年a）

\begin{matrix} δ {E类}_{我, j个} & = δ_{ϕ} {E类}_{我, j个} + δ_{u个} {E类}_{我, j个} = δ_{ϕ} {E类}_{我, j个} + {E类}_{我, j个 对} δ {u个}_{对} = δ_{ϕ} {E类}_{我, j个} + {E类}_{我, 对 j个} δ {u个}_{对} \end{matrix}

（22亿）

现在，将方程式（22）代入方程式(20)使用方程式(19)利用高斯发散定理，可以得到电吉布斯自由能的物理变分

\begin{matrix} δ {¦Β}_{v（v）} G公司 d日 v（v） & = {¦Β}_{v（v）} (- σ_{我 j个, j个} + τ_{我 j个 k个, k个 j个} - σ_{我 j个, j个}^{M（M）} - τ_{我 j个 米, 米 j个}^{M（M）}) δ {u个}_{我} d日 v（v） \\ + {¦Β}_{一} (σ_{我 j个} {n个}_{j个} - τ_{我 j个 k个, k个} {n个}_{j个} + σ_{我 j个}^{M（M）} {n个}_{j个} + τ_{我 j个 米, 米}^{M（M）} {n个}_{j个} + {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{米} {n个}_{j个} τ_{我 j个 米} - {\tilde{D类}}_{j个}^{t吨} ({n个}_{米} τ_{我 j个 米})) δ {u个}_{我} d日 一 \\ + {¦Β}_{一} τ_{我 j个 k个} {n个}_{k个} {n个}_{j个} {\tilde{D类}}^{n个} (δ {u个}_{我}) d日 一 + {¦Β}_{v（v）} (- {D类}_{我, 我} + {D类}_{我 j个, j个 我}) δ ϕ d日 v（v） \\ + {¦Β}_{一} ({D类}_{我} {n个}_{我} - {D类}_{我 j个, j个} {n个}_{我} + {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{我} {n个}_{j个} 问_{我 j个} - {\tilde{D类}}_{我}^{t吨} (问_{我 j个} {n个}_{j个})) δ ϕ d日 一 + {¦Β}_{一} 问_{我 j个} {n个}_{我} {n个}_{j个} {\tilde{D类}}^{n个} (δ ϕ) d日 一, \end{matrix}

(23)

哪里

σ_{我 j个}^{M（M）}

和

τ_{我 j个 米, 米}^{M（M）}

麦克斯韦应力和高阶静电应力是变分过程的直接结果，表示为

\begin{matrix} σ_{我 j个}^{M（M）} & = {E类}_{我} {D类}_{j个} - \frac{1}{2} {E类}_{k个} {D类}_{k个} δ_{我 j个} \end{matrix}

（24a）

\begin{matrix} τ_{我 j个 米, 米}^{M（M）} & = {E类}_{我, k个} 问_{k个 j个} - \frac{1}{2} {E类}_{k个, 我} 问_{k个 我} δ_{我 j个} . \end{matrix}

（24b）

注意，方程式(24个)和(第13页)是相同的。最后，通过替换方程式(23)和(21)到方程式中(20)，我们可以得到控制方程和边界条件

\begin{matrix} {(σ_{我 j个} - τ_{我 j个 米, 米} + σ_{我 j个}^{E类 S公司})}_{, j个} + {（f）}_{我} & = 0, 在里面 v（v） \end{matrix}

（25年）

\begin{matrix} σ_{我 j个, j个}^{E类 S公司} & = 0, 在里面 {v（v）}^{'} \end{matrix}

（25亿）

\begin{matrix} {({D类}_{我} - 问_{我 j个, j个})}_{, 我} & = ρ^{（f）}, 在里面 v（v） \end{matrix}

（25c）

\begin{matrix} {({D类}_{我} - 问_{我 j个, j个})}_{, 我} = ϕ_{, 我 我} & = 0, 在里面 {v（v）}^{'} \end{matrix}

（25天）

Dirichlet边界条件

\begin{matrix} {u个}_{我} & = {\bar{u个}}_{我}, 在 一^{u个} \end{matrix}

（26年a）

\begin{matrix} ϕ & = \bar{ϕ}, 在 一^{ϕ} \end{matrix}

（26亿）

\begin{matrix} {\tilde{D类}}^{n个} ({u个}_{我}) = {u个}_{我, 我} {n个}_{我} & = {\bar{u个}}_{我, 我} {n个}_{我}, 在 一^{u个, n个} \end{matrix}

（26美分）

\begin{matrix} {\tilde{D类}}^{n个} (ϕ) = ϕ_{, 我} {n个}_{我} & = {\bar{ϕ}}_{, 我} {n个}_{我}, 在 一^{ϕ, n个} \end{matrix}

（26天）

Neumann边界条件

\begin{matrix} (σ_{我 j个} - τ_{我 j个 米, 米} + [[σ_{我 j个}^{E类 S公司}]]) {n个}_{j个} + {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{米} {n个}_{j个} τ_{我 j个 米} - {\tilde{D类}}_{j个}^{t吨} ({n个}_{米} τ_{我 j个 米}) & = {t吨}_{我}, 在 一^{σ} \end{matrix}

（27年a）

\begin{matrix} τ_{我 j个 k个} {n个}_{k个} {n个}_{j个} & = {\bar{τ}}_{我}, 在 一^{τ} \end{matrix}

（27亿）

\begin{matrix} [[{D类}_{我} - 问_{我 j个, j个}]] {n个}_{我} + {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{我} {n个}_{j个} 问_{我 j个} - {\tilde{D类}}_{我}^{t吨} (问_{我 j个} {n个}_{j个}) & = - σ^{*}, 在 一^{σ^{*}} \end{matrix}

（27 c）

\begin{matrix} 问_{我 j个} {n个}_{我} {n个}_{j个} & = \bar{q个}, 在 一^{问} \end{matrix}

（27天）

在上面的方程中，虽然它不是平凡的，但我们需要记住，在真空中，位移和极化是不存在的；因此，数量

{D类}_{我} - 问_{我 j个, j个}

退化为

- ϕ_{, 我}

真空中，这也会导致跳跃

[[{D类}_{我} - 问_{我 j个, j个}]]

在边界上

一^{ϕ}

.

为了清楚起见，我们比较了表1.让我们表示由内能导出的强形式U型和电吉布斯能量G公司在里面第2.1节和第2.2节分别为P公式和E公式。与经典（局部）理论相比，在这两种公式中，应变梯度

\nabla ϵ

结果其共轭变量，双应力

τ_{我 j个 k个}

在线性动量、表面牵引力和表面高阶牵引力的平衡中进行铸造。除了经典的Dirichlet边界条件外，边界上还规定了法向位移梯度。同样，电场梯度依赖于E公式的电吉布斯能量，导致其共轭变量，即高阶电位移

问_{我 j个}

，出现在高斯定律和表面电荷的修正中。在边界上也规定了法向的电位梯度。另一方面，在P公式中，当考虑极化梯度时，其共轭变量

{V（V）}_{我 j个}

进入分子内力平衡方程和高阶表面极化，而高斯定律和表面电荷条件与经典理论相比没有变化。此外，如上所述，两种公式给出了相同的麦克斯韦应力张量

σ_{我 j个}^{M（M）}

，而高阶麦克斯韦应力是在其相应的电自变量中给出的。值得注意的是，静电应力张量

σ_{我 j个}^{E类 S公司}

在真空中减少到

σ_{我 j个}^{E类 S公司} = σ_{我 j个}^{M（M）} = ε_{0} ({E类}_{我} {E类}_{j个} - \frac{1}{2} {E类}_{k个} {E类}_{k个} δ_{我 j个}),

(28)

这是电磁学中常见的麦克斯韦应力张量。

此外，我们注意到，可以从文中给出的广义公式推导出不同的公式。例如，在没有静电应力的情况下

σ_{我 j个}^{E类 S公司}

和极化梯度

{P（P）}_{我, j个}

（在P公式中）或电场梯度

{E类}_{我, j个}

（在E公式中），我们可以检索Sharma给出的强形式[30]、毛[39]阿卜杜拉希[38]，邓[40]和阮[43]. 我们总结了以下两种简化配方表2还应该注意的是，在P-制剂中，当外部电场也被省略时，分子内的力平衡被降低到

{E类}_{我}^{L（左）} = ϕ_{, 我}

或电场和电势之间的通常关系

{E类}_{我} = - ϕ_{, 我}

可以检索，并且两个公式的强形式是相同的。为此，用于数值计算的公式只是一个选择问题，因为这两个公式的弱形式是相同的，可以写成

{¦Β}_{v（v）} [σ_{我 j个} δ {u个}_{我, j个} + τ_{我 j个 米} δ {u个}_{我, j个 米} + (ε_{0} {E类}_{我} + {P（P）}_{我}) δ ϕ_{, 我}] d日 v（v） = {¦Β}_{v（v）} {（f）}_{我} δ {u个}_{我} d日 v（v） + {¦Β}_{一^{σ}} {t吨}_{我} δ {u个}_{我} d日 一 - {¦Β}_{v（v）} ρ^{（f）} δ ϕ d日 v（v） - {¦Β}_{一^{σ^{*}}} σ^{*} δ ϕ d日 一

(29)

用于P-配方和

{¦Β}_{v（v）} (σ_{我 j个} δ {u个}_{我, j个} + τ_{我 j个 米} δ {u个}_{我, j个 米} + {D类}_{我} δ ϕ_{, 我}) d日 v（v） = {¦Β}_{v（v）} {（f）}_{我} δ {u个}_{我} d日 v（v） + {¦Β}_{一^{σ}} {t吨}_{我} δ {u个}_{我} d日 一 - {¦Β}_{v（v）} ρ^{（f）} δ ϕ d日 v（v） - {¦Β}_{一^{σ^{*}}} σ^{*} δ ϕ d日 一

(30)

用于E-配方。注意，方程式(29)将根据位移场和电场进行求解。因此，需要通过利用方程（2c）中的本构关系，即。，

{P（P）}_{我} = 一_{我 j个}^{- 1} ϕ_{, j个} - 一_{我 j个}^{- 1} {d日}_{j个 k个 我} ε_{k个 我} - 一_{我 j个}^{- 1} {小时}_{j个 k个 我 米} ε_{k个 我, 米}

另一方面，作为电场

{E类}_{我} = - ϕ_{, 我}

是E公式中的一个独立变量，从方程中直接求解位移场和电位场(30). 因此，在下面的章节中，我们使用不同的数值方法来求解方程中的弱形式(30)简化E-公式。完整配方时(表1)如果选择的话，电吉布斯自由能还是比较有利的。在这种情况下，用电势表示的极化在域积分和表面积分中都可以实现。换句话说，表面极化在P公式中的应用等效于E公式中的电位法向导数。

3.计算建模

挠曲电的建模是为了解决表2或最小化弱形式（方程式(29)或方程式(30)). 为了进行数值计算，应变梯度的存在要求所选的数值方法必须至少确保

{C类}^{0}

-位移场二阶导数的高阶连续性（至少

{C类}^{1}

-连续性）形状函数。在下文中，我们总结了柔性电学中使用的数值方法，包括无网格、混合有限元和IGA。

3.1. 无网格方法

对于柔性电结构的计算评估，一些最初的突出工作是采用局部最大熵（LME）近似的无网格方法，正如Abdollahi等人所做的那样[38]. 通过比较一维问题中的解析和数值能量转换因子，验证了无网格公式。受Ma和Cross实验工作的启发，采用基于LME的无网格方法分析了一个二维金字塔结构[38]和三维金字塔结构[50]. 由于原子力显微镜尖端压入SrTiO而导致应变的逆挠曲电

_{三}

使用基于LME的无网格方法对样品进行了实验研究和模拟[51]. 除了LME的工作外，Bo He等人[52]采用无单元伽辽金法（MLS）对二维柔性电结构进行了分析。然而，由于MLS和LME形状函数都缺乏Kroneckerδ特性，因此仍需要使用无网格方法对柔性电进行更多的研究，这导致对施加Dirichlet边界条件进行额外的处理。此外，在具有材料界面的多材料问题中，一些特殊的技术，如域划分[53]，拉格朗日乘数[54]和跳转功能[55]用于捕捉材料的不连续性。另一方面，高阶Dirichlet边界条件在柔性电学的无网格方法中没有提及。

3.2. 混合有限元法

混合有限元法具有捕捉

{C类}^{1}

利用域中的连续性

{C类}^{0}

有限元。Mao等人首次利用混合有限元[39]为了分析周期性柔性电结构中的柔性电板和裂纹，他们提出了一种混合有限元单元，表示为I9-87。中采用的柔性电配方[39]基于P公式，因此电自由度包括极化和电势。机械自由度包括位移和松弛位移梯度，运动约束由拉格朗日乘子施加。混合FE也由以下作者用于E-配方[40,56,57]. Nanthakumar等人[56]采用一个具有54个自由度的九节点四边形单元来分析和优化二维柔性电结构。邓等人[40]提出了两个三角形单元T37和T45，以及两个四边形单元Q47和Q59。Deng等人提出了一种六面体单元[57]他利用混合有限元分析了微金字塔结构。

混合有限元的主要优点是由于位移梯度自由度的缘故，它在施加高阶Dirichlet边界条件方面很简单；因此，它可以直接在FE商业软件中实现。此外，由于其

{C类}^{0}

-近似下，混合有限元法可用于简化界面和复合柔性电结构的建模。然而，该方法有一个主要缺点：每个元素的自由度多，计算成本昂贵，尤其是在3D中。

3.3. IGA方法

由于操纵基函数连续性的灵活性，IGA被用于求解柔性电问题以进行拓扑优化[41]，软电介质材料[42]和半导体[44]. 在这些工作中，选择NURBS基函数来近似几何和场变量（位移场和电位场）。通过使用节点插入技术[58,59]，可以获得所需的NURBS基函数连续顺序，如所示图1此外，通过控制连续性，IGA方法可以用于接口建模[42]，其中跨界面的几何近似为

{C类}^{0}

-连续性，同时

{C类}^{1}

-域中像往常一样需要连续性。此外，在拓扑优化工作中[41]，假定孔隙为密度比主体材料小得多的固体材料；然而，

{C类}^{1}

在空隙-固相之间施加连续性。对于跨内边界的建模不连续性，可以从[60,61,62]. 该方法已应用于多孔介质、复合材料、多场和多材料问题。在这些作品中，界面元素已经被应用于构建

{C类}^{0}

-沿重新定义界面两侧的连续插值。该单元通过在平衡方程中添加由内部不连续处的本构关系定义的机械牵引力，实现了强制边界条件的可能性。

完美的界面可产生连续的位移和牵引场沿着接口，即。，

\begin{matrix} [[{u个}_{我}]] & = 0 \end{matrix}

（31a）

\begin{matrix} [[{t吨}_{我}]] & = 0, \end{matrix}

（31亿）

哪里

[[]]

表示跳转。弱不连续位移的条件可以通过适当的富集函数实现，如XFEM[63,64]. 此外，在应变梯度弹性问题中，与位移法向导数和双重牵引等高阶项相关的附加连续性条件[65,66]

\begin{matrix} [[{u个}_{我, 我} {n个}_{我}]] & = 0 \end{matrix}

（32a）

\begin{matrix} [[τ_{我 j个 米} {n个}_{米} {n个}_{j个}]] & = 0, \end{matrix}

（32亿）

应该沿着接口满足。此外，微观结构问题中的界面通常被认为是不完美的（位移和/或牵引场沿界面不连续）[67,68]. 然而，这些不完善的界面模型仅限于局部弹性，没有考虑应变梯度。虽然界面元素被应用于力学问题，但同样的概念可以继承并应用于压电和柔性电复合材料等多场问题。

值得讨论IGA方法中边界条件的施加。而Dirichlet边界条件可以通过最小二乘法施加[69]由于NURBS基函数缺乏Kronecker-delta特性，可以通过拉格朗日乘子对位移自由度施加约束来应用位移梯度条件。这个想法经常用于求解应变-颗粒弹性[70]但在柔性电建模工作中尚未讨论。位移梯度条件的数值实现见下一节。

在不考虑电场或极化梯度的简化公式中，不需要施加电位梯度或表面极化。然而，为了建模柔性电传感器，有必要通过施加等电位条件来建模结构表面上的电极。这可以通过使用拉格朗日乘数来实现，如我们在下一节中所示。请注意，在以前使用IGA研究挠曲电的工作中未使用等电位。值得注意的是，最近Codony等人[46]针对柔性电问题，提出了基于分层B样条的浸没边界法，在处理复杂域形状的边界条件方面具有更高的效率。此外，还处理了非光滑边界上的非局部条件（我们忽略了）。

在表3，我们比较了不同数值方法在柔性电建模中的一些特殊特性，包括连续性处理、计算成本（通过每个单元的自由度）、执行高阶Dirichlet边界条件的能力和材料不连续性近似。每种数值方法都有优缺点，应根据所需的应用进行考虑。例如，混合有限元更适合于简化的二维挠曲电问题，而由空洞或夹杂物组成的挠曲电域可以用无网格方法建模。

3.4. 启动位置

在本节中，我们将介绍IGA和混合有限元方法的计算过程，以解决受挠曲电控制的问题，该方法可以用作模拟纳米能量收割机、传感器或致动器的数值工具。基于上述原因，我们考虑在表2，这也是作者的选择[38,40]，因此非常适合进行比较和验证。

混合有限元法

要求

{C类}^{1}

柔性电的四阶偏微分方程导致的连续基函数排除了经典拉格朗日多项式的使用（或至少使其复杂化）。过于确保

{C类}^{1}

连续位移梯度，

ψ

在公式和运动约束中引入，

ψ = \nabla u个

，是使用拉格朗日乘子强加的。机械和静电平衡的弱形式是通过求

u个 \in {u个 = \bar{u个} 在 一^{u个}, u个 \in {H（H）}^{1} (v（v）)}

,

ψ \in {ψ = \bar{ψ} 在 一^{u个, n个}, ψ \in {H（H）}^{1} (v（v）)}

,

ϕ \in {ϕ = \bar{ϕ} 在 一^{ϕ}, ϕ \in {H（H）}^{1} (v（v）)}

,

λ \in {L（左）}^{2} (v（v）)

，以及类似的

外国金融情报机构 \in {外国金融情报机构 = 0 在 一^{u个}, 外国金融机构联合会 \in {H（H）}^{1} (v（v）)}

,

δ ψ \in {δ ψ = 0 在 一^{u个, n个}, δ ψ \in {H（H）}^{1} (v（v）)}

,

δ ϕ \in {δ ϕ = 0 在 一^{ϕ}, δ ϕ \in {H（H）}^{1} (v（v）)}

和

δ λ \in {L（左）}^{2} (v（v）)

，这样，

\begin{matrix} δ Π & = {¦Β}_{Ω} (σ : ε (δ u个) + \hat{τ} ⋮ \hat{η} (δ ψ) - D类 \cdot E类 (δ ϕ)) d日 Ω + {¦Β}_{Ω} λ : (δ ψ - \nabla δ u个) d日 Ω \\ + {¦Β}_{Ω} (ψ - \nabla u个) : δ λ d日 Ω - {¦Β}_{Γ_{N个}} δ u个 \cdot \bar{t吨} d日 Γ - {¦Β}_{Ω} δ u个 \cdot b条 d日 Ω = 0 \end{matrix}

（33）

\begin{matrix} {¦Β}_{Ω} ε (δ u个) : C类 : ε (u个) d日 Ω - {¦Β}_{Ω} ε (δ u个) : e（电子） \cdot E类 (ϕ) d日 Ω - {¦Β}_{Ω} \hat{η} (δ ψ) ⋮ 小时 \cdot E类 (ϕ) d日 Ω \\ + {¦Β}_{Ω} \hat{η} (δ ψ) ⋮ 克 ⋮ \hat{η} (ψ) d日 Ω - {¦Β}_{Ω} E类 (δ ϕ) \cdot e（电子） : ε (u个) d日 Ω - {¦Β}_{Ω} E类 (δ ϕ) \cdot 小时 ⋮ \hat{η} (ψ) d日 Ω \\ - {¦Β}_{Ω} E类 (δ ϕ) \cdot κ \cdot E类 (ϕ) d日 Ω + {¦Β}_{Ω} δ ψ : λ d日 Ω - {¦Β}_{Ω} \nabla δ u个 : λ d日 Ω - {¦Β}_{Ω} δ λ : \nabla u个 d日 Ω \\ + {¦Β}_{Ω} δ λ : ψ d日 Ω = {¦Β}_{Γ_{N个}} δ u个 \cdot \bar{t吨} d日 Γ + {¦Β}_{Ω} δ u个 \cdot b条 d日 Ω \end{matrix}

(34)

哪里

\hat{τ}

=

\frac{\partial U型}{\partial \hat{η}}

和

\hat{η}

是与松弛位移梯度相关的三阶张量

ψ

作为

\hat{η}

=

\frac{1}{2} (ψ_{j个 k个, 我} + ψ_{我 k个, j个})

.约束条件

ψ = \nabla u个

以加权剩余方式施加，包括

δ {¦Β}_{Ω} (ψ - \nabla u个) : λ d日 Ω

在方程式中(33).

Nanthakumar等人提出了九节点等参元[56]. 自由度为

{u个}_{1}

和

{u个}_{2}

在所有节点，松弛位移梯度

ψ_{11}

,

ψ_{12}

,

ψ_{21}

和

ψ_{22}

四角节点和电势

ϕ

在四个角节点处。除了这些自由度之外，该元素还有四个拉格朗日乘数

λ_{11}

,

λ_{12}

,

λ_{21}

和

λ_{22}

在四个角节点处。柔性电气元件如所示图2.双二次形状函数，

{N个}^{q个}

，用于插值位移自由度

{u个}_{1}

和

{u个}_{2}

.双线性形状函数，

{N个}^{我}

，用于插值松弛位移梯度、拉格朗日乘子和电势，

ϕ

.有限元近似如下所示，

\begin{matrix} {u个}^{小时} (x个) = \sum_{我 = 1}^{9} {N个}_{我}^{q个} {u个}_{我}; ψ^{小时} (x个) = \sum_{我 = 1, 三, 5, 7} {N个}_{我}^{我} ψ_{我}; \\ ϕ^{小时} (x个) = \sum_{我 = 1, 三, 5, 7} {N个}_{我}^{我} ϕ_{我}; λ^{小时} (x个) = \sum_{我 = 1, 三, 5, 7} {N个}_{我}^{我} λ_{我} \end{matrix}

(35)

\begin{matrix} δ {u个}^{小时} (x个) = \sum_{我 = 1}^{9} {N个}_{我}^{q个} δ {u个}_{我}; δ ψ^{小时} (x个) = \sum_{我 = 1, 三, 5, 7} {N个}_{我}^{我} δ ψ_{我}; \\ δ ϕ^{小时} (x个) = \sum_{我 = 1, 三, 5, 7} {N个}_{我}^{我} δ ϕ_{我}; δ λ^{小时} (x个) = \sum_{我 = 1, 三, 5, 7} {N个}_{我}^{我} δ λ_{我} \end{matrix}

(36)

替代方程式(35)和(36)，到方程式中(34)根据节点自由度得出以下表达式：

\begin{matrix} δ {u个}^{T型} ({¦Β}_{Ω} {B类}_{u个}^{T型} C类 {B类}_{u个} d日 Ω) u个 + δ {u个}^{T型} (\underset{Ω}{¦Β} {B类}_{u个}^{T型} {e（电子）}^{T型} {B类}_{ϕ} d日 Ω) ϕ + {δ ψ}^{T型} ({¦Β}_{Ω} {H（H）}^{T型} {小时}^{T型} {B类}_{ϕ} d日 Ω) ϕ \\ + {δ ψ}^{T型} ({¦Β}_{Ω} {H（H）}^{T型} 克 H（H） d日 Ω) ψ + {δ ϕ}^{T型} (\underset{Ω}{¦Β} {B类}_{ϕ}^{T型} e（电子） {B类}_{u个} d日 Ω) u个 + {δ ϕ}^{T型} ({¦Β}_{Ω} {B类}_{ϕ}^{T型} 小时 H（H） d日 Ω) ψ \\ - {δ ϕ}^{T型} ({¦Β}_{Ω} {B类}_{ϕ}^{T型} κ {B类}_{ϕ} d日 Ω) ϕ - {δ u个}^{T型} ({¦Β}_{Ω} {B类}_{ψ u个}^{T型} {N个}^{我} d日 Ω) λ + {δ ψ}^{T型} ({¦Β}_{Ω} {N个}^{我}^{T型} {N个}^{我} d日 Ω) λ \\ - {δ λ}^{T型} ({¦Β}_{Ω} {N个}^{我}^{T型} {B类}_{ψ u个} d日 Ω) u个 + {δ λ}^{T型} ({¦Β}_{Ω} {N个}^{我}^{T型} {N个}^{我} d日 Ω) ψ = {δ u个}^{T型} ({¦Β}_{Γ_{N个}} {N个}^{q个} \bar{t吨} d日 Γ) + {δ u个}^{T型} ({¦Β}_{Ω} {N个}^{q个} b条 d日 Ω) \end{matrix}

（37）

哪里

{B类}_{u个} = [\begin{matrix} \frac{\partial {N个}^{q个}}{\partial x个} & 0 \\ 0 & \frac{\partial {N个}^{q个}}{\partial 年} \\ \frac{\partial {N个}^{q个}}{\partial 年} & \frac{\partial {N个}^{q个}}{\partial x个} \end{matrix}]

,

{B类}_{ϕ} = [\begin{matrix} \frac{\partial {N个}^{我}}{\partial x个} \\ \frac{\partial {N个}^{我}}{\partial 年} \end{matrix}]

,

H（H） = [\begin{matrix} \frac{\partial {N个}^{我}}{\partial x个} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{\partial {N个}^{我}}{\partial x个} \\ 0 & \frac{1}{2} \frac{\partial {N个}^{我}}{\partial x个} & \frac{1}{2} \frac{\partial {N个}^{我}}{\partial x个} & 0 \\ \frac{\partial {N个}^{我}}{\partial 年} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{\partial {N个}^{我}}{\partial 年} \\ 0 & \frac{1}{2} \frac{\partial {N个}^{我}}{\partial 年} & \frac{1}{2} \frac{\partial {N个}^{我}}{\partial 年} & 0 \end{matrix}]

,

{B类}_{ψ u个} = [\begin{matrix} \frac{\partial {N个}^{我}}{\partial x个} & 0 \\ 0 & \frac{\partial {N个}^{我}}{\partial x个} \\ \frac{\partial {N个}^{我}}{\partial 年} & 0 \\ 0 & \frac{\partial {N个}^{我}}{\partial 年} \end{matrix}]

为了获得柔性电结构中的位移和电势，需要求解的代数方程如下：，

[\begin{matrix} {K（K）}_{你} & 0 & {K（K）}_{u个 ϕ} & {K（K）}_{u个 λ} \\ 0 & {K（K）}_{ψ ψ} & {K（K）}_{ψ ϕ} & {K（K）}_{ψ λ} \\ {K（K）}_{ϕ u个} & {K（K）}_{ϕ ψ} & - {K（K）}_{ϕ ϕ} & 0 \\ {K（K）}_{λ u个} & {K（K）}_{λ ψ} & 0 & 0 \end{matrix}] [\begin{matrix} u个 \\ ψ \\ ϕ \\ λ \end{matrix}] = [\begin{matrix} {F类}_{u个} \\ 0 \\ 0 \\ 0 \end{matrix}]

(38)

\begin{matrix} {K（K）}_{你} = {¦Β}_{Ω} {B类}_{u个}^{T型} C类 {B类}_{u个} d日 Ω; \\ {K（K）}_{u个 ϕ} = \underset{Ω}{¦Β} {B类}_{u个}^{T型} {e（电子）}^{T型} {B类}_{ϕ} d日 Ω = {K（K）}_{ϕ u个}^{T型}; \\ {K（K）}_{u个 λ} = - {¦Β}_{Ω} {B类}_{ψ u个}^{T型} {N个}_{λ} d日 Ω = {K（K）}_{λ u个}^{T型} \\ {K（K）}_{ψ ϕ} = {¦Β}_{Ω} {H（H）}^{T型} 小时 {B类}_{ϕ} d日 Ω = {K（K）}_{ϕ ψ}^{T型} \\ {K（K）}_{ψ ψ} = {¦Β}_{Ω} {B类}_{ψ}^{T型} 克 {B类}_{ψ} d日 Ω \\ {K（K）}_{ψ λ} = {¦Β}_{Ω} {N个}^{我^{T型}} {N个}_{λ} d日 Ω = {K（K）}_{λ ψ}^{T型} \\ {K（K）}_{ϕ ϕ} = {¦Β}_{Ω} {B类}_{ϕ}^{T型} κ {B类}_{ϕ} d日 Ω; \\ {F类}_{u个} = {¦Β}_{Γ} {N个}^{q个}^{T型} \bar{t吨} d日 Γ + {¦Β}_{Ω} {N个}^{q个} b条 d日 Ω \end{matrix}

3.5. IGA方法

在Galerkin方法中，试验函数和测试函数均由NURBS基函数近似

\begin{matrix} {u个}^{(e（电子）)} = {N个}_{u个} \tilde{u个}, δ {u个}^{(e（电子）)} = {N个}_{u个} δ \tilde{u个}, \end{matrix}

（39a）

\begin{matrix} ϕ^{(e（电子）)} = {N个}_{ϕ} \tilde{ϕ}, δ ϕ^{(e（电子）)} = {N个}_{ϕ} δ \tilde{ϕ}, \end{matrix}

（39亿）

在哪儿

{N个}_{u个}

和

{N个}_{ϕ}

是包含元素的局部NURBS基函数的矩阵e（电子）定义为

\begin{matrix} {N个}_{u个}^{我} & = [\begin{matrix} {N个}_{我} & 0 & 0 \\ 0 & {N个}_{我} & 0 \\ 0 & 0 & {N个}_{我} \end{matrix}], {N个}_{u个} = [\begin{matrix} {N个}_{u个}^{1} & {N个}_{u个}^{2} & \dots & {N个}_{u个}^{n个} \end{matrix}], \end{matrix}

（40年）

\begin{matrix} {N个}_{ϕ} & = [\begin{matrix} {N个}_{1} & {N个}_{2} & \dots & {N个}_{n个} \end{matrix}], \end{matrix}

（40亿）

哪里

{N个}_{我}

是控制点处的NURBS基函数我th和n个元素的控制点数量。在上面的方程式中，

\tilde{u个}

,

\tilde{ϕ}

,

δ \tilde{u个}

和

δ \tilde{u个}

分别是与试验函数和测试函数的位移和电势相关的节点控制变量，可以用矩阵形式表示为

\begin{matrix} \tilde{u个} & = {[{\tilde{u个}}_{x个}^{1} {\tilde{u个}}_{年}^{1} {\tilde{u个}}_{z（z）}^{1} \dots {\tilde{u个}}_{x个}^{n个} {\tilde{u个}}_{年}^{n个} {\tilde{u个}}_{z（z）}^{n个}]}^{T型}, \end{matrix}

(41)

\begin{matrix} δ \tilde{u个} & = {[δ {\tilde{u个}}_{x个}^{1} δ {\tilde{u个}}_{年}^{1} δ {\tilde{u个}}_{z（z）}^{1} \dots δ {\tilde{u个}}_{x个}^{n个} δ {\tilde{u个}}_{年}^{n个} δ {\tilde{u个}}_{z（z）}^{n个}]}^{T型}, \end{matrix}

(42)

\begin{matrix} \tilde{ϕ} & = {[{\tilde{ϕ}}^{1} {\tilde{ϕ}}^{2} \dots {\tilde{ϕ}}^{n个}]}^{T型}, \end{matrix}

(43)

\begin{matrix} δ \tilde{ϕ} & = {[δ {\tilde{ϕ}}^{1} δ {\tilde{ϕ}}^{2} \dots δ {\tilde{ϕ}}^{n个}]}^{T型} \end{matrix}

(44)

应变场的近似

ε = {[ε_{11} ε_{22} ε_{33} 2 ε_{23} 2 ε_{13} 2 ε_{12}]}^{T型}

，电场

E类 = {[{E类}_{1} {E类}_{2} {E类}_{三}]}^{T型}

，应变梯度

\nabla ε = {[\frac{\partial ε}{\partial x个} \frac{\partial ε}{\partial 年} \frac{\partial ε}{\partial z（z）}]}^{T型}

和电场梯度

\nabla E类 = {[\frac{\partial E类}{\partial x个} \frac{\partial E类}{\partial 年} \frac{\partial E类}{\partial z（z）}]}^{T型}

可以通过使用微分算子获得

{B类}_{u个}

,

{B类}_{ϕ}

,

{H（H）}_{u个}

和

{H（H）}_{ϕ}

根据方程式（39）得出：

\begin{matrix} ε^{(e（电子）)} & = {B类}_{u个} \tilde{u个}, δ ε^{(e（电子）)} = {B类}_{u个} δ \tilde{u个}, \end{matrix}

（45年）

\begin{matrix} {E类}^{(e（电子）)} & = {B类}_{ϕ} \tilde{ϕ}, δ {E类}^{(e（电子）)} = {B类}_{ϕ} δ \tilde{ϕ}, \end{matrix}

（45亿）

\begin{matrix} \nabla ε^{(e（电子）)} & = {H（H）}_{u个} \tilde{u个}, δ \nabla ε^{(e（电子）)} = {H（H）}_{u个} δ \tilde{u个}, \end{matrix}

（45美分）

\begin{matrix} \nabla {E类}^{(e（电子）)} & = {H（H）}_{ϕ} \tilde{ϕ}, δ \nabla {E类}^{(e（电子）)} = {H（H）}_{ϕ} δ \tilde{ϕ}, \end{matrix}

（45天）

在哪儿

{B类}_{u个}

,

{B类}_{ϕ}

,

{H（H）}_{u个}

和

{H（H）}_{ϕ}

包含NURBS基函数的一阶和二阶偏导数，分别定义为

{B类}_{u个}^{我} = [\begin{matrix} \frac{\partial {N个}_{我}}{\partial x个} & 0 & 0 \\ 0 & \frac{\partial {N个}_{我}}{\partial 年} & 0 \\ 0 & 0 & \frac{\partial {N个}_{我}}{\partial z（z）} \\ 0 & \frac{\partial {N个}_{我}}{\partial z（z）} & \frac{\partial {N个}_{我}}{\partial 年} \\ \frac{\partial {N个}_{我}}{\partial z（z）} & 0 & \frac{\partial {N个}_{我}}{\partial x个} \\ \frac{\partial {N个}_{我}}{\partial 年} & \frac{\partial {N个}_{我}}{\partial x个} & 0 \end{matrix}], {B类}_{u个} = [\begin{matrix} {B类}_{u个}^{1} & {B类}_{u个}^{2} & \dots & {B类}_{u个}^{n个} \end{matrix}], {H（H）}_{u个} = [\begin{matrix} \frac{\partial {B类}_{u个}^{1}}{\partial x个} & \frac{\partial {B类}_{u个}^{2}}{\partial x个} & \dots & \frac{\partial {B类}_{u个}^{n个}}{\partial x个} \\ \frac{\partial {B类}_{u个}^{1}}{\partial 年} & \frac{\partial {B类}_{u个}^{2}}{\partial 年} & \dots & \frac{\partial {B类}_{u个}^{n个}}{\partial 年} \\ \frac{\partial {B类}_{u个}^{1}}{\partial z（z）} & \frac{\partial {B类}_{u个}^{2}}{\partial z（z）} & \dots & \frac{\partial {B类}_{u个}^{n个}}{\partial z（z）} \end{matrix}]

(46)

{B类}_{ϕ}^{我} = [\begin{matrix} - \frac{\partial {N个}_{我}}{\partial x个} \\ - \frac{\partial {N个}_{我}}{\partial 年} \\ - \frac{\partial {N个}_{我}}{\partial z（z）} \end{matrix}], {B类}_{ϕ} = [\begin{matrix} {B类}_{ϕ}^{1} & {B类}_{ϕ}^{2} & \dots & {B类}_{ϕ}^{n个} \end{matrix}], {H（H）}_{ϕ} = [\begin{matrix} \frac{\partial {B类}_{ϕ}^{1}}{\partial x个} & \frac{\partial {B类}_{ϕ}^{2}}{\partial x个} & \dots & \frac{\partial {B类}_{ϕ}^{n个}}{\partial x个} \\ \frac{\partial {B类}_{ϕ}^{1}}{\partial 年} & \frac{\partial {B类}_{ϕ}^{2}}{\partial 年} & \dots & \frac{\partial {B类}_{ϕ}^{n个}}{\partial 年} \\ \frac{\partial {B类}_{ϕ}^{1}}{\partial z（z）} & \frac{\partial {B类}_{ϕ}^{2}}{\partial z（z）} & \dots & \frac{\partial {B类}_{ϕ}^{n个}}{\partial z（z）} \end{matrix}]

（47）

注意矩阵

{H（H）}_{u个}

和

{H（H）}_{ϕ}

包含关于基函数物理坐标的二阶导数。详细信息见附录A现在，基于离散化有限元，方程中的弱形式(30)可以重新表示为

⋃_{e（电子） = 1}^{{n个}_{e（电子）}} [{¦Β}_{{v（v）}^{(e（电子）)}} (δ ϵ^{T型} σ + δ \nabla ϵ^{T型} τ - δ {E类}^{T型} D类) d日 v（v） = {¦Β}_{v（v）} δ {u个}^{T型} （f） d日 v（v） + {¦Β}_{一^{σ, (e（电子）)}} δ {\tilde{u个}}^{T型} \bar{t吨} d日 一 + {¦Β}_{v（v）} {\tilde{ϕ}}^{T型} ρ^{（f）} d日 v（v） - {¦Β}_{一^{ϕ, (e（电子）)}} δ {\tilde{ϕ}}^{T型} {\bar{σ}}^{*} d日 一] .

（48）

最后，将方程（39）和（45）代入方程(48)利用测试函数的任意性，可以得到线性方程组

[\begin{matrix} {K（K）}_{u个 u个} & {K（K）}_{u个 ϕ} \\ {K（K）}_{ϕ u个} & {K（K）}_{ϕ ϕ} \end{matrix}] [\begin{matrix} \tilde{u个} \\ \tilde{ϕ} \end{matrix}] = [\begin{matrix} {（f）}_{u个} \\ {（f）}_{ϕ} \end{matrix}],

(49)

其中，刚度矩阵和力向量计算为

\begin{matrix} {K（K）}_{u个 u个} & = ⋃_{e（电子） = 1}^{{n个}_{e（电子）}} {¦Β}_{{v（v）}^{(e（电子）)}} {B类}_{u个}^{T型} c（c） {B类}_{u个} d日 v（v）, \end{matrix}

（50年）

\begin{matrix} {K（K）}_{u个 ϕ} & = ⋃_{e（电子） = 1}^{{n个}_{e（电子）}} {¦Β}_{{v（v）}^{(e（电子）)}} (- {B类}_{u个}^{T型} e（电子） {B类}_{ϕ} - {H（H）}_{u个}^{T型} μ {B类}_{ϕ}) d日 v（v）, \end{matrix}

（50亿）

\begin{matrix} {K（K）}_{ϕ u个} & = ⋃_{e（电子） = 1}^{{n个}_{e（电子）}} {¦Β}_{{v（v）}^{(e（电子）)}} (- {B类}_{ϕ}^{T型} {e（电子）}^{T型} {B类}_{u个} - {B类}_{ϕ}^{T型} μ^{T型} {H（H）}_{u个}) d日 v（v）, \end{matrix}

（50美分）

\begin{matrix} {K（K）}_{ϕ ϕ} & = ⋃_{e（电子） = 1}^{{n个}_{e（电子）}} {¦Β}_{{v（v）}^{(e（电子）)}} - {B类}_{ϕ}^{T型} κ {B类}_{ϕ} d日 v（v）, \end{matrix}

（50天）

\begin{matrix} {（f）}_{u个} & = ⋃_{e（电子） = 1}^{{n个}_{e（电子）}} ({¦Β}_{v（v）} {N个}_{u个}^{T型} （f） d日 v（v） + {¦Β}_{一^{σ, (e（电子）)}} {N个}_{u个}^{T型} \bar{t吨} d日 一), \end{matrix}

（50秒）

\begin{matrix} {（f）}_{ϕ} & = ⋃_{e（电子） = 1}^{{n个}_{e（电子）}} ({¦Β}_{v（v）} {N个}_{u个}^{T型} ρ^{（f）} d日 v（v） + {¦Β}_{一^{ϕ, (e（电子）)}} - {N个}_{ϕ}^{T型} \bar{τ} d日 一) . \end{matrix}

（50英尺）

注意挠曲电张量

μ

对于立方对称晶体

κ

表示为

\begin{matrix} μ & = [\begin{matrix} μ_{11} & μ_{12} & μ_{12} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & μ_{44} & 0 & 0 & 0 & 0 & μ_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & μ_{44} & μ_{12} & μ_{11} & μ_{12} & 0 & 0 & 0 & 0 & 0 & 0 & μ_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & μ_{44} & 0 & 0 & 0 & 0 & μ_{44} & 0 & 0 & μ_{12} & μ_{12} & μ_{11} & 0 & 0 & 0 \end{matrix}] \end{matrix}

(51)

挠曲电张量的对称性可以在[71,72]. 对于立方晶体，挠曲电张量的非零分量为

μ_{1111}

,

μ_{1221}

和

μ_{1122}

而其他具有循环指数的系数相等（例如。，

μ_{1221} = μ_{1331}

)索引的奇数为零（例如。，

μ_{1112} = 0

) [73]. 在矩阵表示法中，系数

μ_{我 我 我 我}

,

μ_{j个 我 我 j个}

和

μ_{我 我 j个 j个}

表示为

μ_{11}

,

μ_{12}

和

μ_{44}

分别为[38].

4.数值示例

4.1. 空心圆柱体

在第一个数值例子中，我们考虑一个无限长的空心圆柱体，其内外表面承受机械位移和电势，如图所示图3a.利用轴对称性，对四分之一的空心圆柱体进行建模，如所示图3b、对称边界条件在

{x个}_{1} = 0

和

{x个}_{2} = 0

这个问题是梯度弹性中经常使用的基准[74,75]，并已扩展到柔性电[39,40]. 在我们的工作中，我们调整了两种内能U型来自[39]和电吉布斯能G公司来自[40]如下

U型 = \frac{1}{2} λ ε_{j个 j个} ε_{k个 k个} + μ ε_{j个 k个} ε_{j个 k个} + \frac{1}{2} 我^{2} (λ ε_{j个 j个, 我} ε_{k个 k个, 我} + 2 μ ε_{j个 k个, 我} ε_{j个 k个, 我}) + {\hat{（f）}}_{12} ε_{j个 j个, 我} {P（P）}_{我} + 2 {\hat{（f）}}_{44} ε_{我 j个, 我} {P（P）}_{j个} + \frac{1}{2} α {P（P）}_{我} {P（P）}_{我}

（52）

和

G公司 = \frac{1}{2} λ ε_{j个 j个} ε_{k个 k个} + μ ε_{j个 k个} ε_{j个 k个} + \frac{1}{2} 我^{2} (λ ε_{j个 j个, 我} ε_{k个 k个, 我} + 2 μ ε_{j个 k个, 我} ε_{j个 k个, 我}) - μ_{12} ε_{j个 j个, 我} {E类}_{我} - 2 μ_{44} ε_{我 j个, 我} {E类}_{j个} - \frac{1}{2} κ {E类}_{我} {E类}_{我},

(53)

哪里

λ

和

μ

是两个Lamé参数，我是长度刻度，

μ_{12}

和

μ_{44}

是两个非零挠曲电系数，以及

{\hat{（f）}}_{1}

和

{\hat{（f）}}_{2}

是两个挠曲电耦合系数，通过磁化率与挠曲电系数相关

χ = κ / ε_{0} - 1

使得

μ_{12} = χ {\hat{（f）}}_{12}

,

μ_{44} = χ {\hat{（f）}}_{44}

,

α

和

κ

是相互的磁化率和介电常数，也与磁化率有关。这些材料参数可在表4此外，由于几何形状和边界条件，位移和电势仅取决于径向第页空心圆柱体的

\begin{matrix} {u个}_{第页} (第页) & = {c（c）}_{1} 第页 + {c（c）}_{2} \frac{1}{第页} + {c（c）}_{三} 我_{1} (\frac{第页}{我_{0}}) + {c（c）}_{4} {K（K）}_{1} (\frac{第页}{我_{0}}), \end{matrix}

（54年a）

\begin{matrix} ϕ (第页) & = {c（c）}_{5} 自然对数 (第页) + {c（c）}_{6} + \frac{μ_{11}}{κ} (\frac{{d日 u个}_{第页}}{d日 第页} + \frac{{u个}_{第页} (第页)}{第页}), \end{matrix}

（54亿）

哪里

我_{0}^{2} = 我^{2} + \frac{μ_{11}^{2}}{(λ + μ) κ}

,

我_{1}

和

{K（K）}_{1}

分别是第一类和第二类修正贝塞尔函数。系数

{c（c）}_{1}, {c（c）}_{2}, {c（c）}_{三}, {c（c）}_{4}, {c（c）}_{5}, {c（c）}_{6}

由边界条件决定

\begin{matrix} {{u个}_{第页} │}_{第页 = {第页}_{我}} = {u个}_{{第页}_{我}}, {{u个}_{第页} │}_{第页 = {第页}_{o个}} & = {u个}_{{第页}_{o个}}, \end{matrix}

（55a）

\begin{matrix} {ϕ_{第页} │}_{第页 = {第页}_{我}} = ϕ_{{第页}_{我}}, {ϕ_{第页} │}_{第页 = {第页}_{o个}} & = ϕ_{{第页}_{o个}}, \end{matrix}

（55亿）

\begin{matrix} {τ_{第页 第页 第页} │}_{第页 = {第页}_{我}} = {τ_{第页 第页 第页} │}_{第页 = {第页}_{o个}} & = 0 \end{matrix}

（55美分）

对于数值建模，两种公式表2用IGA求解，此后表示为

我 G公司 A类 (P（P）)

和

我 G公司 A类 (E类)

分别从电焓和吉布斯能量得到的公式。此外，我们还用混合有限元法演示了第二种公式的结果。

对于IGA方法，我们采用

25 \times 25

具有以下特性的二阶NURBS基函数

{C类}^{1}

对四分之一空心圆柱体进行连续建模，如所示图3b、并在左边缘和下边缘施加额外的对称边界条件。而齐次和非齐次Dirichlet边界条件很容易施加[69]位移梯度的边界条件要求约束与边界控制点及其相邻点相关的自由度。例如，中底部边缘的对称边界条件图3b读取

\begin{matrix} {{u个}_{2} │}_{{x个}_{2} = 0} & = 0, \end{matrix}

（56a）

\begin{matrix} {{u个}_{1, 2} │}_{{x个}_{2} = 0} & = 0, \end{matrix}

（56亿）

\begin{matrix} {{u个}_{2, 1} │}_{{x个}_{2} = 0} & = 0, \end{matrix}

（56美分）

式中，方程式（56c）根据方程式自动完成(第56页). 为了应用方程（56b），节点参数（与

{x个}_{1}

方向

{u个}_{{x个}_{1}}

)两个底层控制点的相同，即。，

{u个}_{{x个}_{1}}^{底部} = {u个}_{{x个}_{1}}^{相邻的 - 底部}

这可以通过使用拉格朗日乘数来实现。另一方面，对于具有节点位移梯度的混合FEM，可以毫不费力地施加位移梯度边界条件。

图4给出了三种方法得到的位移和电势的数值结果，显示了它们之间的良好一致性。为了更仔细地观察，径向位移的数值结果

{u个}_{R（右）}

，电势

ϕ

，径向应变

ε_{第页 第页}

和周向应变

ε_{θ θ}

沿径向有角度

θ = 45^{°}

给出并与中的解析解进行比较图5可以看出，数值结果和分析结果之间存在良好的一致性。

4.2. 悬臂梁

在这个例子中，我们考虑了一个受到右端点荷载作用的二维悬臂梁，它具有两种不同的电极配置，作为纳米发电机的演示[38]. 如所示图6，光束在左端受到约束，而对于两个电边界条件设置了两种不同的电极配置。在类型1边界条件下，光束右端的电极接地，即电势

ϕ

规定为零。在类型2条件下，在顶部和底部电极之间施加电位差，例如，顶部的电位固定为零，而底部被诱导具有未知的等电位V（V）利用拉格朗日乘子求解。本研究中的材料参数为单一钛酸钡（BaTi

{O（运行）}_{三}

)晶体采用自[38]并在中呈现表5机械变形可以通过纳米发电机中的机电耦合效应转换为电能。这种能量转换的性能可以通过机电耦合系数进行评估

{k个}_{效率}

表示为静电能之间的比率

{E类}_{电子}

和弹性能

{E类}_{埃拉斯}

{k个}_{效率} = \frac{{E类}_{电子}}{{E类}_{埃拉斯}} = \frac{\frac{1}{2} {¦Β}_{v（v）} E类 \cdot κ \cdot E类 d日 v（v）}{\frac{1}{2} {¦Β}_{v（v）} ϵ \cdot C类 \cdot ϵ d日 v（v）},

（57）

因此，归一化压电常数定义为

\bar{e（电子）} = \frac{{k个}_{效率}}{{k个}_{压电的}},

(58)

哪里

{k个}_{压电的}

获得方式为

{k个}_{效率}

不考虑挠曲电效应。物理上，归一化压电常数表明挠曲电效应对能量转换的增强。图7证明了归一化压电常数

\bar{e（电子）}

作为标准化厚度

\bar{小时} = - {e（电子）}_{31} 小时 / μ_{12}

在类型1配置下会有所不同。简化情况下的数值结果，其中

κ_{33}

和

μ_{12}

非零以模拟一维问题，与解析解进行比较，其中

{k个}_{效率} = \sqrt{\frac{{e（电子）}_{31}^{2} + 12 {(μ_{12} / 小时)}^{2}}{κ_{33} Y（Y）}}

结果表明，两者吻合良好，这说明随着束层厚度的减小，柔性电性大大增强。此外，这种增强也可以在完整的2D模型中看到，这也与Abdollahi等人的工作非常一致[38].

我们进一步研究了归一化有效压电常数的差异

\bar{e（电子）}

使用两种不同的电气边界条件，并在图8在这两种边界条件设置中，随着梁厚度的减小，归一化有效压电常数增加，这是由于柔性电的尺寸效应性质。值得注意的是，1型边界条件似乎比2型更有效。此外，IGA和混合有限元方法对两种边界条件类型的性能给出了相同的预测。为了更仔细地观察，我们给出了两种电势边界条件的电势分布，这两种边界条件是由IGA和混合FEM在图9梁的高度和长度为

小时 = 1.2468 \times 10^{- 7} (米)

和

L（左） = 10 小时

分别是。两种方法的结果是一致的。对于1类边界条件，电势差在梁的左端最大，应变梯度最大。对于第2类边界条件，顶部电极的感应电势约为

- 9 μ V（V）

.

此外，我们还使用IGA对2型电边界条件下的悬臂梁进行了三维模拟。假设3D梁具有方形横截面，其中宽度等于高度小时。梁的长度以及施加在右侧表面上的力与二维情况相同。图10显示了感应电势。在2D情况下观察到了类似的响应，其中最大电势位于左端。然而，顶部电极的电势值为

3.21 V（V）

.

我们还从电边界条件研究了逆挠曲电效应。当省略机械点荷载时，均匀电场的大小

│ E类 │ = 8 中压 / 米

通过将电势设置为

V（V） = - 8 小时 (中压)

在底部电极上。在这种情况下，梁起到了致动器的作用，并产生了弯曲曲率，如所示图11b.由逆挠曲电效应引起的弯曲梁的计算曲率与[38]以及来自Bursian和OI的实验[76]. 对光束厚度方向电场分布的进一步研究如所示图11a、这也与[38]. 可以观察到电场的非均匀分布，特别是在光束顶部和底部边缘附近的区域，这会产生较高的电场梯度，从而导致变形。注意，所考虑的挠曲电束不考虑电阻负载[10]或整流电路[77,78]. 在[77]分析了压电结构及其电路连接。此外，在[78]对电荷型哈密顿量和电压型哈密尔顿量进行了比较，以识别压电能量采集器的输出功率和电压。开发一种完全组合的柔性电发生器和接口电路在未来可能具有重要意义。

4.3. 截断金字塔

本节介绍了在不同类型边界条件下压缩的二维截断金字塔的数值结果[38].图12描述了一个截顶金字塔的示意图，其顶面用零电势接地，而底面连接有一个电极，从而在其上产生等电位条件。考虑了两个机械边界条件，即底面约束在的刚性边界条件

{x个}_{2}

方向和顶面承受均匀载荷F类和柔性边界条件，其中底面和顶面都承受均匀的压缩载荷F类对于数值建模，问题通过二次型离散化

{C类}^{1}

-连续性NURBS基函数，如所示图12b以及混合有限元法图12c。

电势的数值结果

ϕ

和压缩应变

ε_{年 年}

IGA和混合有限元法的刚性边界条件见图13IGA的柔性边界条件如所示图14.两种方法之间以及与参考结果之间可以观察到良好的一致性[38]. 在刚性模型中，我们的IGA和混合FEM方法预测底部电极的等电位为

2.7

和

2.5

mV，而参考结果为

2.6

毫伏[38]. 同样，在柔性模型中，底部等电位为

5.7

mV，符合

5.6

mV来自[38].

此外，还研究了截断金字塔的三维版本。几何参数保持为2D情况，因此顶部和底部曲面的尺寸现在为

一_{1} \times 一_{1}

和

一_{2} \times 一_{2}

分别是。考虑柔性边界条件，其中施加的力F类与2D案例相同。IGA的数值结果如所示图15如图所示图15c、三维模型中的电势与二维模型中的非常相似。然而，顶部电极上的感应电势现在是

6.918 V（V）

，比2D模型预测的结果大三个数量级。

5.讨论

自发现以来，柔性电因其有助于开发新型机电耦合器件而具有重要意义。为了充分利用这一效应，必须采用实验和模拟方法。从仿真的角度出发，本文给出了弯曲电的广义公式，作为梯度弹性的延伸。值得注意的是，在柔性电学中涉及四阶偏微分方程的边值问题

{C类}^{1}

连续性是由电焓和电吉布斯自由能得到的。为了克服高阶连续性，我们采用了等几何分析和混合有限元方法，并给出了详细的实现方法。我们在二维和三维柔性电问题上使用这些方法解决了基准问题。提供了源代码，可以用于研究更复杂的问题。

作者贡献

概念化、X.Z.、B.H.N.、S.S.N.、T.Q.T.、N.A.和T.R。；方法、X.Z.、B.H.N.、S.S.N.、T.Q.T.、N.A.和T.R。；软件、B.H.N.、S.S.N.和T.Q.T。；验证、B.H.N.、S.S.N.和T.Q.T。；形式分析，B.H.N.、S.S.N.和T.Q.T。；调查、B.H.N.、S.S.N.和T.Q.T。；书面、X.Z.、B.H.N.、S.S.N.、T.Q.T.、N.A.和T.R。；可视化、B.H.N.、S.S.N.和T.Q.T。；和监督，X.Z.和T.R.所有作者都已阅读并同意手稿的出版版本。

基金

这项研究没有得到外部资助。

致谢

作者对沙特国王大学杰出科学家研究计划（DSFP）资助这项工作表示感谢。

利益冲突

作者声明没有利益冲突。

附录A.NURBS基函数的二阶导数

关于基函数的物理坐标的二阶导数的计算可以通过首先注意到基函数关于参数坐标的一阶导数是根据雅可比矩阵计算的来实现

[\begin{matrix} \frac{\partial {N个}_{我}}{\partial ξ} \\ \frac{\partial {N个}_{我}}{\partial η} \\ \frac{\partial {N个}_{我}}{\partial ζ} \end{matrix}] = [\begin{matrix} \frac{\partial x个}{\partial ξ} & \frac{\partial 年}{\partial ξ} & \frac{\partial z（z）}{\partial ξ} \\ \frac{\partial x个}{\partial η} & \frac{\partial 年}{\partial η} & \frac{\partial z（z）}{\partial η} \\ \frac{\partial x个}{\partial ζ} & \frac{\partial 年}{\partial ζ} & \frac{\partial z（z）}{\partial ζ} \end{matrix}] [\begin{matrix} \frac{\partial {N个}_{我}}{\partial x个} \\ \frac{\partial {N个}_{我}}{\partial 年} \\ \frac{\partial {N个}_{我}}{\partial z（z）} \end{matrix}]

（A1）

接下来，通过对参数坐标求二阶导数并重新排列，根据以下公式确定了对物理坐标的二阶导数

[\begin{matrix} {(x个,_{ξ})}^{2} & {(年,_{ξ})}^{2} & {(z（z）,_{ξ})}^{2} & 2 年_{, ξ} {z（z）}_{, ξ} & 2 {x个}_{, ξ} {z（z）}_{, ξ} & 2 {x个}_{, ξ} 年_{, ξ} \\ {(x个,_{η})}^{2} & {(年,_{η})}^{2} & {(z（z）,_{ξ})}^{2} & 2 年_{, η} {z（z）}_{, η} & 2 {x个}_{, η} {z（z）}_{, η} & 2 {x个}_{, η} 年_{, η} \\ {(x个,_{ζ})}^{2} & {(年,_{ζ})}^{2} & {(z（z）,_{ζ})}^{2} & 2 年_{, ζ} {z（z）}_{, ζ} & 2 {x个}_{, ζ} {z（z）}_{, ζ} & 2 {x个}_{, ζ} 年_{, ζ} \\ {x个}_{, η} {x个}_{, ζ} & 年_{, η} 年_{, ζ} & {z（z）}_{, η} {z（z）}_{, ζ} & 年_{, η} {z（z）}_{, ζ} + 年_{, ζ} {z（z）}_{, η} & {x个}_{, η} {z（z）}_{, ζ} + {x个}_{, ζ} {z（z）}_{, η} & {x个}_{, η} 年_{, ζ} + {x个}_{, ζ} 年_{, η} \\ {x个}_{, ζ} {x个}_{, ξ} & 年_{, ζ} 年_{, ξ} & {z（z）}_{, ζ} {z（z）}_{, ξ} & 年_{, ζ} {z（z）}_{, ξ} + 年_{, ξ} {z（z）}_{, ζ} & {x个}_{, ζ} {z（z）}_{, ξ} + {x个}_{, ξ} {z（z）}_{, ζ} & {x个}_{, ζ} 年_{, ξ} + {x个}_{, ξ} 年_{, ζ} \\ {x个}_{, ξ} {x个}_{, η} & 年_{, ξ} 年_{, η} & {z（z）}_{, ξ} {z（z）}_{, η} & 年_{, ξ} {z（z）}_{, η} + 年_{, η} {z（z）}_{, ξ} & {x个}_{, ξ} {z（z）}_{, η} + {x个}_{, η} {z（z）}_{, ξ} & {x个}_{, ξ} 年_{, η} + {x个}_{, η} 年_{, ξ} \end{matrix}] [\begin{matrix} {N个}_{我, x个 x个} \\ {N个}_{我, 年 年} \\ {N个}_{我, z（z） z（z）} \\ {N个}_{我, 年 z（z）} \\ {N个}_{我, x个 z（z）} \\ {N个}_{我, x个 年} \end{matrix}] = [\begin{matrix} {N个}_{我, ξ ξ} \\ {N个}_{我, η η} \\ {N个}_{我, ζ ζ} \\ {N个}_{我, η ζ} \\ {N个}_{我, ξ ζ} \\ {N个}_{我, ξ η} \end{matrix}] - [\begin{matrix} {x个}_{, ξ ξ} & 年_{, ξ ξ} & {z（z）}_{, ξ ξ} \\ {x个}_{, η η} & 年_{, η η} & {z（z）}_{, η η} \\ {x个}_{, ζ ζ} & 年_{, ζ ζ} & {z（z）}_{, ζ ζ} \\ {x个}_{, η ζ} & 年_{, η ζ} & {z（z）}_{, η ζ} \\ {x个}_{, ζ ξ} & 年_{, ζ ξ} & {z（z）}_{, ζ ξ} \\ {x个}_{, ξ η} & 年_{, ξ η} & {z（z）}_{, ξ η} \end{matrix}] [\begin{matrix} {N个}_{我, x个} \\ {N个}_{我, 年} \\ {N个}_{我, z（z）} \end{matrix}] .

（A2）

参考文献

Tolpygo，K。金刚石晶体的长波振荡，包括远程力。苏联。物理学-固态 1963,4, 1297–1305. [谷歌学者]
Kogan，S.M.晶体中载流子非均匀变形和声散射期间的压电效应。苏联。物理学-固态 1964,5, 2069–2070. [谷歌学者]
黄，W。；严，X。；Kwon，S.R.公司。；张，S。；袁，F.G。；Jiang，X.使用Ba0。64Sr0.36TiO3用于传感。申请。物理学。莱特。 2012,101, 252903. [谷歌学者] [交叉参考]
Kwon，S.R.公司。；黄，W。；张，S。；袁，F.G。；Jiang，X.使用钛酸锶钡研究柔性电传声器。J.微型机械。马略能。 2016,26, 045001. [谷歌学者] [交叉参考] [绿色版本]
梅鲁波，V.I。；吉法德，B。；Seveno，R。；Tabellout，M。；Kassiba，A.软聚氨酯薄膜中的挠曲电响应及其在大曲率传感中的应用。J.应用。物理学。 2017,122, 144101. [谷歌学者] [交叉参考]
英国巴斯卡。；北班纳吉。；阿卜杜拉希，A。；王，Z。；施洛姆，D.G。；Rijnders，G。；Catalan，G.硅上的柔性电微电子机械系统。自然纳米技术。 2016,11, 263. [谷歌学者] [交叉参考]
张，S。；刘凯。；徐，M。；Shen，S.一种弯曲谐振柔性电动执行器。申请。物理学。莱特。 2017,111, 082904. [谷歌学者] [交叉参考]
公元雷伊。；塞维奥，P。；Herrera-Valencia，E.基于柔性电膜的机械能量收集的仿生模型。物理学。版本E 2013,87, 022505. [谷歌学者] [交叉参考] [绿色版本]
蒋，X。；黄，W。；Zhang，S.柔性电纳米发电机：材料、结构和器件。纳米能源 2013,2, 1079–1092. [谷歌学者] [交叉参考]
邓，Q。；M.Kammoun。；埃尔图尔克，A。；Sharma，P.纳米柔性电能采集。国际固体结构杂志。 2014,51, 3218–3225. [谷歌学者] [交叉参考] [绿色版本]
Choi，S.B。；Kim，G.W.压电振动能量采集器用聚偏氟乙烯薄膜中挠曲电响应的测量。《物理学杂志》。D申请。物理学。 2017,50, 075502. [谷歌学者] [交叉参考]
X·梁。；胡，S。；Shen，S.利用压电和挠曲电收集纳米级机械能。聪明的马特。结构。 2017,26，035050。[谷歌学者] [交叉参考] [绿色版本]
朱，R。；王，Z。；马，H。；袁，G。；Wang，F。；郑，Z。；Kimura，H.基于坚固的自极化铁电纤维的无极化能量采集器。纳米能源 2018,50, 97–105. [谷歌学者] [交叉参考]
卢，H。；Bark，C.W。；De Los Ojos，D.E。；Alcala，J。；Eom、C.B。；加泰罗尼亚语，G。；Gruverman，A.铁电极化的机械书写。科学类 2012,336, 59–61. [谷歌学者] [交叉参考] [公共医学] [绿色版本]
卢，H。；刘，S。；Ye，Z。；Yasui，S。；Funakubo，H。；Rappe，A.M。；Gruverman，A.机械极化开关中的不对称。申请。物理学。莱特。 2017,110, 222903. [谷歌学者] [交叉参考] [绿色版本]
郭，R。；沈，L。；Wang，H。；Lim，Z。；卢·W。；杨，P。；Gruverman，A。；Venkatesan，T。；冯，Y.P。；Chen，J.通过界面工程和柔性电效应定制BaTiO3薄膜的自极化。高级主管。接口 2016,三, 1600737. [谷歌学者] [交叉参考]
杨，M.M。；Kim，D.J。；Alexe，M.Flexo光伏效应。科学类 2018,360, 904–907. [谷歌学者] [交叉参考] [绿色版本]
刘，Y。；陈，J。；邓，H。；胡，G。；朱，D。；Dai，N.应变Bi 2 Te 3薄膜中的异常热电性。科学。代表。 2016,6，32661页。[谷歌学者] [交叉参考] [绿色版本]
Nguyen，T.D。；毛，S。；Yeh，Y.W。；Purohit，P.K。；McAlpine，M.C.纳米柔性电。高级主管。 2013,25, 946–974. [谷歌学者] [交叉参考]
Zubko，P。；加泰罗尼亚语，G。；Tagantsev，A.K.固体中的挠曲电效应。每年。修订版材料。物件。 2013,43. [谷歌学者] [交叉参考] [绿色版本]
尤丁，P。；Tagantsev，A.固体挠曲电基础。纳米技术 2013,24, 432001. [谷歌学者] [交叉参考] [公共医学]
Krichen，S。；Sharma，P.Flexoelectricity:一种不同寻常的机电耦合透视图。J.应用。机械。 2016,83, 030801. [谷歌学者] [交叉参考] [绿色版本]
Maranganti，R。；北卡罗来纳州夏尔马。；Sharma，P.纳米尺度非局部尺寸效应导致的非压电材料中的机电耦合：格林函数解和嵌入夹杂物。物理学。版本B 2006,74, 014110. [谷歌学者] [交叉参考] [绿色版本]
Mindlin，R.D.弹性电介质中的极化梯度。国际固体结构杂志。 1968,4, 637–642. [谷歌学者] [交叉参考]
Sahin，E。；Dost，S.A空间色散弹性电介质的应变理论。国际工程科学杂志。 1988,26, 1231–1245. [谷歌学者] [交叉参考]
Majdoub，M。；Sharma，P。；在压电纳米结构的窄尺寸范围内，能量收集显著增强。物理学。版本B 2008,78, 121407. [谷歌学者] [交叉参考] [绿色版本]
Majdoub，M。；Sharma，P。；乔安·伊恩，T.Erratum：压电纳米结构窄尺寸范围内能量收集的显著增强【物理评论B 78，121407（R）（2008）】。物理学。版本B 2009,79, 159901. [谷歌学者] [交叉参考]
Majdoub，M。；Sharma，P。；Cagin，T.由于柔性电效应，纳米结构中的尺寸依赖性压电和弹性增强。物理学。版本B 2008,77, 125424. [谷歌学者] [交叉参考] [绿色版本]
Majdoub，医学硕士。；Sharma，P。；乔安，T.Erratum：由于柔性电效应，纳米结构中的尺寸依赖性压电和弹性增强【物理评论B77，125424（2008）】。物理学。版本B 2009,79. [谷歌学者] [交叉参考]
北卡罗来纳州夏尔马。；兰迪斯，C。；Sharma，P.不使用压电材料的压电薄膜超晶格。J.应用。物理学。 2010,108, 024304. [谷歌学者] [交叉参考] [绿色版本]
胡，S。；Shen，S.《柔性电效应纳米电介质中的变分原理和控制方程》。科学。中国物理、机械。阿童木。 2010,53, 1497–1504. [谷歌学者] [交叉参考]
Kuang，Z.B.电致伸缩和磁致伸缩材料中的一些问题。机械学报。索里达·辛。 2007,20，219–227页。[谷歌学者] [交叉参考]
Kuang，Z。有限变形下电弹性介质中的一些变分原理。科学。中国Ser。物理、机械。阿童木。 2008,51, 1390–1402. [谷歌学者] [交叉参考]
Kuang，Z.b。电弹性分析中的内能变分原理和控制方程。国际固体结构杂志。 2009,46, 902–911. [谷歌学者] [交叉参考] [绿色版本]
Kuang，Z.B.先生。理论电弹性; 施普林格：德国柏林/海德堡，2014年。[谷歌学者]
Liu，L.关于连续介质静电的能量公式。J.机械。物理学。固体 2013,61, 968–990. [谷歌学者] [交叉参考]
Liu，L.连续磁电弹性的能量公式及其应用。J.机械。物理学。固体 2014,63, 451–480. [谷歌学者] [交叉参考]
阿卜杜拉希，A。；佩科，C。；Millan，D。；阿罗约，M。；Arias，I.电介质固体中挠曲电效应的计算评估。J.应用。物理学。 2014,116, 093502. [谷歌学者] [交叉参考]
毛，S。；Purohit，P.K。；Aravas，N.压电性和挠曲电性的混合有限元公式。程序。R.Soc.A公司 2016,472, 20150879. [谷歌学者] [交叉参考] [绿色版本]
邓，F。；邓，Q。；于伟（Yu，W.）。；Shen，S.柔性电固体的混合有限元。J.应用。机械。 2017,84, 081004. [谷歌学者] [交叉参考]
加西米，H。；帕克，H.S。；Rabczuk，T.基于水平集的IGA公式，用于柔性电材料的拓扑优化。计算。方法应用。机械。工程师。 2017,313，239–258。[谷歌学者] [交叉参考]
泰国，T.Q。；Rabczuk，T。；庄，X。柔性电和软材料的大变形等几何方法。计算。方法应用。机械。工程师。 2018,341，718–739页。[谷歌学者] [交叉参考]
Nguyen，B。；庄，X。；Rabczuk，T.压电和柔性电复合材料中Maxwell-Wagner松弛特性的数值模型。计算。结构。 2018,208, 75–91. [谷歌学者] [交叉参考]
Nguyen，B。；庄，X。；Rabczuk，T.基于NURBS的半导体非线性电梯度弹性公式。计算。方法应用。机械。工程师。 2019,346, 1074–1095. [谷歌学者] [交叉参考]
Yvonnet，J。；Liu，L.有限应变下软电介质中挠曲电和麦克斯韦应力建模的数值框架。计算。方法应用。机械。工程师。 2017,313，450–482。[谷歌学者] [交叉参考]
Codony，D。；奥马尔科。；Fernández-Méndez，S。；Arias，I.柔性电学的浸没边界层次B样条方法。arXiv公司 2019，arXiv:1902.02567。[谷歌学者] [交叉参考] [绿色版本]
罗伊·P。；Roy，D.挠曲电和损伤的周动力学模型。申请。数学。模型。 2019,68, 82–112. [谷歌学者] [交叉参考]
Toupin，R.A.弹性电介质。J.配给。机械。分析。 1956,5, 849–915. [谷歌学者] [交叉参考]
R.D.明德林。；Eshel，N.关于线性弹性中的第一个应变-颗粒理论。国际固体结构杂志。 1968,4, 109–124. [谷歌学者] [交叉参考]
阿卜杜拉希，A。；Millán，D。；佩科，C。；阿罗约，M。；Arias，I.重新审视金字塔压缩以量化挠曲电：三维模拟研究。物理学。版本B 2015,91, 104103. [谷歌学者] [交叉参考] [绿色版本]
阿卜杜拉希，A。；北多明戈。；阿里亚斯一世。；Catalan，G.Converse柔性电在非压电材料中产生较大的压电响应力显微镜信号。国家公社。 2019,10, 1266. [谷歌学者] [交叉参考] [绿色版本]
他，B。；贾夫瓦吉，B。；Zhuang，X.使用无单元伽辽金方法表征复合材料中的柔性。能源 2019,12, 271. [谷歌学者] [交叉参考] [绿色版本]
刘总。无网格方法：超越有限元方法; CRC出版社：美国佛罗里达州博卡拉顿，2009年。[谷歌学者]
科尔德斯，L。；Moran，B.无元素Galerkin方法中材料不连续性的处理。计算。方法应用。机械。工程师。 1996,139, 75–89. [谷歌学者] [交叉参考]
Krongauz，Y。；Belytschko，T.不连续导数的EFG近似。国际期刊数字。方法工程。 1998,41, 1215–1233. [谷歌学者] [交叉参考]
南塔库马，S。；庄，X。；帕克，H.S。；Rabczuk，T.柔性电结构的拓扑优化。J.机械。物理学。固体 2017,105, 217–234. [谷歌学者] [交叉参考]
邓，F。；邓，Q。；Shen，S.柔性电三维混合有限元。J.应用。机械。 2018,85, 031009. [谷歌学者] [交叉参考]
莱斯·皮格尔（Les Piegl，W.T.）。NURBS手册; 施普林格：德国柏林/海德堡，1997年。[谷歌学者]
科特雷尔，J。；休斯·T。；巴兹列夫斯，Y。等几何分析：走向CAD和FEA的集成; John Wiley&Sons Inc.：美国新泽西州霍博肯，2009年。[谷歌学者]
G.T.卡马乔。；Ortiz，M.脆性材料冲击损伤的计算模型。国际固体结构杂志。 1996,33, 2899–2938. [谷歌学者] [交叉参考]
Schrefler，文学学士。；塞奇，S。；Simoni，L.关于多场问题（包括内聚断裂）中的自适应细化技术。计算。方法应用。机械。工程师。 2006,195, 444–461. [谷歌学者] [交叉参考]
de Borst，R.《裂缝和压裂多孔介质中的流体流动：统一观点》。机械。Res.Commun公司。 2017,80, 47–57. [谷歌学者] [交叉参考]
北苏库马尔。；肖普，D.L。；莫尔斯，N。；Belytschko，T.在扩展有限元方法中通过水平集建模孔和夹杂物。计算。方法应用。机械。工程师。 2001,190, 6183–6200. [谷歌学者] [交叉参考] [绿色版本]
北苏库马尔。；黄，Z。；Prévost，J.H。；Suo，Z.双材料界面裂纹的单位富集划分。国际期刊数字。方法工程。 2004,59, 1075–1102. [谷歌学者] [交叉参考]
Shu，J.Y。；金，W.E。；Fleck，N.A.具有应变梯度效应的材料的有限元。国际期刊数字。方法工程。 1999,44, 373–391. [谷歌学者] [交叉参考] [绿色版本]
Soh，A.K。；Wanji，C.微观结构应变梯度理论的有限元公式和C0-1斑片试验。国际期刊数字。方法工程。 2004,61, 433–454. [谷歌学者] [交叉参考]
顾S。；He，Q.C.耦合多场现象的界面不连续关系及其在薄界面作为非完美界面建模中的应用。J.机械。物理学。固体 2011,59, 1413–1426. [谷歌学者] [交叉参考]
刘，J.T。；他，B。；顾，S.T。；He，Q.C.在XFEM和LSM上下文中实现基于物理的一般弹性非完美界面模型。国际期刊数字。方法工程。 2018,115, 1499–1525. [谷歌学者] [交叉参考]
Nguyen，副总裁。；Anitescu，C。；博尔达斯公司。；Rabczuk，T.等几何分析：概述和计算机实现方面。数学。计算。模拟。 2015,117, 89–116. [谷歌学者] [交叉参考] [绿色版本]
菲舍尔，P。；Klassen，M。；Mergheim，J。；斯坦曼，P。；Müller，R.二维梯度弹性的等几何分析。计算。机械。 2011,47, 325–334. [谷歌学者] [交叉参考]
舒，L。；魏，X。；Pang，T。；姚，X。；Wang，C.晶体介质中挠曲电系数的对称性。J.应用。物理学。 2011,110, 104106. [谷歌学者] [交叉参考]
Le Quang，H。；He，Q.C.挠曲电张量所有可能的旋转对称的数量和类型。程序。R.Soc.伦敦。数学。物理学。工程科学。 2011,467, 2369–2386. [谷歌学者] [交叉参考]
Hong，J。；Vanderbilt，D.第一原理理论和挠曲电的计算。物理学。版本B 2013,88, 174107. [谷歌学者] [交叉参考] [绿色版本]
Gao，X.L。；Park，S.简化应变梯度弹性理论的变分公式及其在受压厚壁圆筒问题中的应用。国际固体结构杂志。 2007,44, 7486–7499. [谷歌学者] [交叉参考]
Aravas，N.一类梯度弹性模型的平面应变问题——应力函数方法。J.弹性。 2011,104, 45–70. [谷歌学者] [交叉参考]
E.布赖恩。；由于极化引起的铁电膜的曲率变化。苏联。物理学。固态，苏联 1968,10, 1121–1124. [谷歌学者]
Lumentut，M。；Shu，Y.多分段智能平板能量采集器的统一机电有限元动力学分析：电路连接模式。机械表演。 2018,229, 4575–4604. [谷歌学者] [交叉参考] [绿色版本]
Lumentut，M。；Howard，I.压电收割机的机电动力学方程。机械表演。 2017,228，631–650。[谷歌学者] [交叉参考] [绿色版本]

样品可用性：作者可根据要求提供详细数据和代码。

图1。(一)

{C类}^{0}

基于节点向量的连续二次NURBS基函数

Ξ = \{0, 0, 0, 1 / 三, 1 / 三, 2 / 三, 2 / 三, 1, 1, 1\}

; 和(b条)

{C类}^{1}

基于节点向量的连续二次NURBS基函数

Ξ = \{0, 0, 0, 1 / 三, 2 / 三, 1, 1, 1\}

图1。(一)

{C类}^{0}

基于节点向量的连续二次NURBS基函数

Ξ = \{0, 0, 0, 1 / 三, 1 / 三, 2 / 三, 2 / 三, 1, 1, 1\}

; 和(b条)

{C类}^{1}

基于节点向量的连续二次NURBS基函数

Ξ = \{0, 0, 0, 1 / 三, 2 / 三, 1, 1, 1\}

图2。九节点柔性电元件。圆圈的自由度为

{u个}_{1}

和

{u个}_{2}

; 在方格处

ψ_{11}

,

ψ_{12}

,

ψ_{21}

,

ψ_{22}

和

ϕ

; 在三角形处

λ_{11}

,

λ_{12}

,

λ_{21}

、和

λ_{22}

.

图2。九节点柔性电元件。圆圈的自由度为

{u个}_{1}

和

{u个}_{2}

; 在方格处

ψ_{11}

,

ψ_{12}

,

ψ_{21}

,

ψ_{22}

和

ϕ

; 在三角形处

λ_{11}

,

λ_{12}

,

λ_{21}

、和

λ_{22}

.

图3。无限长空心圆柱的示意图，其中内半径和外半径为

{第页}_{我} = 10 μ 米

和

{第页}_{0} = 20 μ 米

分别是。内表面受到径向位移的影响

{u个}_{{第页}_{我}} = 0.045 μ 米

接地，即。，

ϕ_{{第页}_{我}} = 0 V（V）

，而径向位移

{u个}_{{第页}_{o个}} = 0.05 μ 米

和电势

ϕ_{{第页}_{o个}} = 1 V（V）

应用于外表面。实施对称边界条件，以便

{{u个}_{1} │}_{{x个}_{1} = 0} = {{u个}_{1, 2} │}_{{x个}_{1} = 0} = {{u个}_{2, 1} │}_{{x个}_{1} = 0} = 0

和

{{u个}_{2} │}_{{x个}_{2} = 0} = {{u个}_{1, 2} │}_{{x个}_{2} = 0} = {{u个}_{2, 1} │}_{{x个}_{2} = 0} = 0

.

图3。无限长空心圆柱体的示意图，其中内半径和外半径为

{第页}_{我} = 10 μ 米

和

{第页}_{0} = 20 μ 米

分别是。内表面受到径向位移的影响

{u个}_{{第页}_{我}} = 0.045 μ 米

接地，即。，

ϕ_{{第页}_{我}} = 0 V（V）

，而径向位移

{u个}_{{第页}_{o个}} = 0.05 μ 米

和电势

ϕ_{{第页}_{o个}} = 1 V（V）

应用于外表面。实施对称边界条件，以便

{{u个}_{1} │}_{{x个}_{1} = 0} = {{u个}_{1, 2} │}_{{x个}_{1} = 0} = {{u个}_{2, 1} │}_{{x个}_{1} = 0} = 0

和

{{u个}_{2} │}_{{x个}_{2} = 0} = {{u个}_{1, 2} │}_{{x个}_{2} = 0} = {{u个}_{2, 1} │}_{{x个}_{2} = 0} = 0

.

图4。径向位移和电位分布。

图5。径向位移的精确结果与数值结果的比较

{u个}_{第页}

，电势

ϕ

，径向应变

ε_{第页 第页}

和周向应变

ε_{ϕ ϕ}

作为半径的函数第页.

图5。径向位移的精确结果与数值结果的比较

{u个}_{第页}

，电势

ϕ

，径向应变

ε_{第页 第页}

和周向应变

ε_{ϕ ϕ}

作为半径的函数第页.

图6。(一)空心圆柱体示意图。施加了两种不同类型的电气边界条件。上图描述了在右边施加零电势的1类边界。2类边界如下图所示，其中上边缘施加零电位，下边缘施加未知等电位（由于电极）。(b条)IGA离散化（顶部）和混合FEM离散化（底部）。

图7。利用类型1边界条件获得的有效归一化压电常数。

图7。利用类型-1边界条件获得的有效归一化压电常数。

图8。不同电边界条件下的有效压电常数。

图9。1类和2类电气边界条件下的电位分布。

图9。1型和2型电边界条件下的电势分布。

图10。三维弯曲柔性电子束的电势。

图11。电气负载下的梁。

图12。截断金字塔的示意图和IGA离散化。

图13。刚性边界条件。

图14。柔性边界条件。

图15。3D金字塔。

表1。强形式的比较。

{L（左）}_{我} = {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{米} {n个}_{j个} τ_{我 j个 米} - {\tilde{D类}}_{j个}^{t吨} ({n个}_{米} τ_{我 j个 米})

,

M（M） = {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{我} {n个}_{j个} 问_{我 j个} - {\tilde{D类}}_{我}^{t吨} (问_{我 j个} {n个}_{j个})

.

表1。强形式的比较。

{L（左）}_{我} = {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{米} {n个}_{j个} τ_{我 j个 米} - {\tilde{D类}}_{j个}^{t吨} ({n个}_{米} τ_{我 j个 米})

,

M（M） = {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{我} {n个}_{j个} 问_{我 j个} - {\tilde{D类}}_{我}^{t吨} (问_{我 j个} {n个}_{j个})

.

		$U型 (ϵ, P（P）, \nabla ϵ, \nabla P（P）)$	$G公司 (ϵ, E类, \nabla ϵ, \nabla E类)$
余额	动量inv（v）	${(σ_{我 j个} - τ_{我 j个米, 米} + σ_{我 j个}^{E类 S公司})}_{, j个} + {（f）}_{我} = 0$	${(σ_{我 j个} - τ_{我 j个 k个, k个} + σ_{我 j个}^{E类 S公司})}_{, j个} + {（f）}_{我} = 0$
	麦克斯韦公司 ${v（v）}^{'}$	$σ_{我 j个, j个}^{E类 S公司} = 0$	$σ_{我 j个, j个}^{E类 S公司} = 0$
	高斯定律v（v）	$- ε_{0} ϕ_{, 我我} + {P（P）}_{我, 我} = ρ^{（f）}$	${({D类}_{我} - 问_{我 j个, j个})}_{, 我} = ρ^{（f）}$
	分子内作用力[24]	${V（V）}_{我 j个, j个} + {E类}_{我}^{L（左）} - ϕ_{, 我} + {E类}_{我}^{e（电子） x个 t吨} = 0$
	高斯定律 ${v（v）}^{'}$	$ϕ_{, 我我} = 0$	$ϕ_{, 我我} = 0$
迪里克莱BC	位移开 $一^{u个}$	${u个}_{我} = {\bar{u个}}_{我}$	${u个}_{我} = {\bar{u个}}_{我}$
	电势开 $一^{ϕ}$	$ϕ = \bar{ϕ}$	$ϕ = \bar{ϕ}$
	位移法向导数开 $一^{u个, n个}$	${u个}_{我, 我} {n个}_{我} = {\bar{u个}}_{我, 我} {n个}_{我}$	${u个}_{我, 我} {n个}_{我} = {\bar{u个}}_{我, 我} {n个}_{我}$
	上的潜在正态导数 $一^{ϕ, n个}$		$ϕ_{, 我} {n个}_{我} = {\bar{ϕ}}_{, 我} {n个}_{我}$
	表面极化开启 $一^{P（P）}$	${P（P）}_{我} = {\bar{P（P）}}_{我}$
纽曼BCs	地面牵引力 $一^{σ}$	$(σ_{我 j个} - τ_{我 j个米, 米} + [[σ_{我 j个}^{E类 S公司}]]) {n个}_{j个} + {L（左）}_{我} = {t吨}_{我}$	$(σ_{我 j个} - τ_{我 j个 k个, k个} + [[σ_{我 j个}^{E类 S公司}]]) {n个}_{j个} + {L（左）}_{我} = {t吨}_{我}$
	表面电荷 $一^{σ^{*}}$	${n个}_{我} (- ε_{0} [[ϕ_{, 我}]] + {P（P）}_{我}) = - σ^{*}$	$[[{D类}_{我} - {D类}_{我 j个, j个}]] {n个}_{我} + M（M） = - σ^{*}$
	高阶抓地力 $一^{τ}$	$τ_{我 j个米} {n个}_{米} {n个}_{j个} = {\bar{τ}}_{我}$	$τ_{我 j个米} {n个}_{米} {n个}_{j个} = {\bar{τ}}_{我}$
	高阶表面电荷 $一^{问}$		$问_{我 j个} {n个}_{我} {n个}_{j个} = \bar{q个}$
	高阶电场开启 $一^{V（V）}$	${V（V）}_{我 j个} {n个}_{j个} = {\bar{v（v）}}_{我}$

表2。简化强形式的比较。

{L（左）}_{我} = {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{米} {n个}_{j个} τ_{我 j个 米} - {\tilde{D类}}_{j个}^{t吨} ({n个}_{米} τ_{我 j个 米})

,

M（M） = {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{我} {n个}_{j个} 问_{我 j个} - {\tilde{D类}}_{我}^{t吨} (问_{我 j个} {n个}_{j个})

.

表2。简化强形式的比较。

{L（左）}_{我} = {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{米} {n个}_{j个} τ_{我 j个 米} - {\tilde{D类}}_{j个}^{t吨} ({n个}_{米} τ_{我 j个 米})

,

M（M） = {\tilde{D类}}_{我}^{t吨} ({n个}_{我}) {n个}_{我} {n个}_{j个} 问_{我 j个} - {\tilde{D类}}_{我}^{t吨} (问_{我 j个} {n个}_{j个})

.

		$U型 (ϵ, P（P）, \nabla ϵ)$ 参考文献[30,39]	$G公司 (ϵ, E类, \nabla ϵ)$ 参考文献[38,40,43]
余额	动量inv（v）	${(σ_{我 j个} - τ_{我 j个米, 米} + σ_{我 j个}^{E类 S公司})}_{, j个} + {（f）}_{我} = 0$	${(σ_{我 j个} - τ_{我 j个 k个, k个} + σ_{我 j个}^{E类 S公司})}_{, j个} + {（f）}_{我} = 0$
	高斯定律v（v）	$- ε_{0} ϕ_{, 我我} + {P（P）}_{我, 我} = ρ^{（f）}$	${({D类}_{我} - 问_{我 j个, j个})}_{, 我} = ρ^{（f）}$
	分子内作用力[24]	${V（V）}_{我 j个, j个} + {E类}_{我}^{L（左）} - ϕ_{, 我} + {E类}_{我}^{e（电子） x个 t吨} = 0$
	高斯定律 ${v（v）}^{'}$	$ϕ_{, 我我} = 0$	$ϕ_{, 我我} = 0$
迪里克莱BC	位移开 $一^{u个}$	${u个}_{我} = {\bar{u个}}_{我}$	${u个}_{我} = {\bar{u个}}_{我}$
	电势开 $一^{ϕ}$	$ϕ = \bar{ϕ}$	$ϕ = \bar{ϕ}$
	位移法向导数开 $一^{u个, n个}$	${u个}_{我, 我} {n个}_{我} = {\bar{u个}}_{我, 我} {n个}_{我}$	${u个}_{我, 我} {n个}_{我} = {\bar{u个}}_{我, 我} {n个}_{我}$
纽曼BCs	地面牵引力 $一^{σ}$	$(σ_{我 j个} - τ_{我 j个米, 米} + [[σ_{我 j个}^{E类 S公司}]]) {n个}_{j个} + {L（左）}_{我} = {\bar{t吨}}_{我}$	$(σ_{我 j个} - τ_{我 j个 k个, k个} + [[σ_{我 j个}^{E类 S公司}]]) {n个}_{j个} + {L（左）}_{我} = {\bar{t吨}}_{我}$
	表面电荷 $一^{σ^{*}}$	${n个}_{我} (- ε_{0} [[ϕ_{, 我}]] + {P（P）}_{我}) = - σ^{*}$	$[[{D类}_{我} - {D类}_{我 j个, j个}]] {n个}_{我} + M（M） = - σ^{*}$
	开启更高阶牵引力 $一^{τ}$	$τ_{我 j个米} {n个}_{米} {n个}_{j个} = {\bar{τ}}_{我}$	$τ_{我 j个米} {n个}_{米} {n个}_{j个} = {\bar{τ}}_{我}$

表3。不同数值方法的比较。

(一)

和

(b条)

分别表示2D和3D中的二次元素。

表3。不同数值方法的比较。

(一)

和

(b条)

分别表示2D和3D中的二次元素。

	无网格	政府间协议	混合有限元法
连续性	取决于内核函数	结操作	节点位移梯度
国防部	取决于影响范围	27（a）或108（b）	54（问题54[56]); 37（T37）、47（Q47）、45（T45）、59（Q59）[40]；233（3D）[57]
高阶Dirichlet BC	-	拉格朗日倍增管	直接征收
${C类}^{0}$ 接口建模	拉格朗日乘法器，丰富的函数	多批次	-

表4。中描述的问题的材料参数图3a。

Y（Y）	$ν$	我	$μ_{12}$	$μ_{44}$	$κ$
$139 \times 10^{9} 帕$	0.3	$2 \times 10^{- 6} 米$	$1 \times 10^{- 6} C类 / 米$	$1 \times 10^{- 6} C类 / 米$	$1 \times 10^{- 9} F类 / 米$

表5。中描述的问题的材料参数图6.

Y（Y）	$ν$	${e（电子）}_{31}$	$μ_{12}$	$κ_{11}$	$κ_{33}$	F类
$100 \times 10^{9} 帕$	0.37	−4.4 $C类 / 米^{2}$	$1 \times 10^{- 6} C类 / 米$	11 $\times 10^{- 9} F类 / 米$	12.48 $\times 10^{- 9} F类 / 米$	100 $\times 10^{- 6} N个$

分享和引用

MDPI和ACS样式

庄，X。；Nguyen，B.H。；南塔库马，S.S。；Tran，T.Q。；阿拉伊兰，N。；拉布祖克，T。挠曲电的计算模型——综述。能源 2020,13, 1326.https://doi.org/10.3390/en13061326

AMA风格

Zhung X、Nguyen BH、Nanthakumar SS、Tran TQ、Alajlan N、Rabczuk T。挠曲电的计算模型——综述。能源. 2020; 13(6):1326.https://doi.org/10.3390/en13061326

芝加哥/图拉宾风格

庄、小英、平惠阮、Subbiah Srivilliputtur Nanthakumar、Thai Quoc Tran、Naif Alajlan和Timon Rabczuk。2020年，《挠曲电计算模型——综述》能源13，编号6:1326。https://doi.org/10.3390/en13061326

请注意，从2016年第一期开始，该杂志使用文章编号而不是页码。请参阅更多详细信息在这里.

文章菜单

挠曲电的计算模型——综述

摘要

1.简介

2.地面工程制定

2.1. 物理变分原理

2.2. 替代配方

3.计算建模

3.1. 无网格方法

3.2. 混合有限元法

3.3. IGA方法

3.4. 启动位置

混合有限元法

3.5. IGA方法

4.数值示例

4.1. 空心圆柱体

4.2. 悬臂梁

4.3. 截断金字塔

5.讨论

作者贡献

基金

致谢

利益冲突

附录A.NURBS基函数的二阶导数

参考文献

分享和引用

文章指标

文章访问统计

更多信息

指导方针

MDPI计划

遵循MDPI