基于

摘要

采用C¹连续性Bell三角形单元，进行了挠曲电效应的机理研究，求解了厚壁筒的平面应变边值问题，并利用挠曲电结构的非对称性产生压电效应。通过算例分析，验证了Bell单元的准确性和收敛性，研究了梯度弹性理论中内禀尺度对挠曲电结构变形的影响，分析了挠曲电的尺寸效应。该工作为挠曲电效应分析和结构设计提供了基础和依据。

关键词

C¹型有限元; 挠曲电; 梯度弹性理论; 尺寸效应

力电耦合效应广泛存在与各种材料之中，包括人工和天然材料。自然界的很多生物功能都是力电耦合的，如听觉感知、与动作有关的神经元突起等。这些现象激发了对此相关机理的研究，并开发可以模仿它们的新材料。在所有力电耦合效应中，研究最多的就是压电效应，它是应变和极化之间的力电耦合效应。

$$p_i=e_{ijk}{\epsilon }_{jk}$$

（1）

式中：$p_i$是电极化；$e_{ijk}$为三阶压电系数；${\epsilon }_{jk}$为应变。压电效应是由晶体内原子的相对位移引起的，只存在于非中心对称晶体中。

不同于压电效应的应变与极化的耦合，挠曲电效应是应变梯度与极化或者应变与极化梯度的耦合。

$$p_i=f_{ijkl}{\eta }_{jkl}$$

（2）

式中：$f_{ijkl}$为四阶挠曲电系数；${\eta }_{jkl}$为高阶应变。挠曲电效应产生的原因是晶体材料局部对称性的变化。传统宏观尺度上应变梯度较小，挠曲电效应可以忽略不计，但在微纳尺度上结构可以产生很大的应变梯度，使得挠曲电效应在纳米尺度上和压电效应相当。挠曲电效应存在于所有晶体类型中，并在超过居里温度时仍然出现。相比压电效应，挠曲电效应更为普遍，而且它还能产生压电不具备的其他力电功能，关于挠曲电效应理论及相关应用可见以下文献综述^［1-7］。

固体材料结构随着特征尺度的减小表现出明显的尺度效应，但经典连续介质力学本构关系中不包含任何尺度参量，无法描述微构件力学性能的尺寸效应，因此需要引入长度标度，如梯度弹性理论。Mindlin在1965年提出了一个梯度弹性理论^［8］，该理论中，应变能密度被认为是应变和一阶、二阶应变梯度的函数，从而应力$\mathit{\sigma }$、应变$\mathit{\epsilon }$的本构关系可以写成^［9］：

$$\mathit{\sigma }=\lambda (\mathit{I}-c_{\mathrm{1}}\mathit{I}{\nabla }^{\mathrm{2}}-c_{\mathrm{2}}\nabla \otimes \nabla )\mathrm{t}\mathrm{r}\mathit{\epsilon }+\mathrm{2}G(\mathrm{1}-c_{\mathrm{3}}{\nabla }^{\mathrm{2}})\mathit{\epsilon }$$

（3）

式中：$c_{\mathrm{1}}$、$c_{\mathrm{2}}$、$c_{\mathrm{3}}$是三个具有长度平方量纲的独立梯度参数；$\lambda$和$G$为Lamé常数；$\mathit{I}$为单位张量；$\nabla$为微分算子；${\nabla }^{\mathrm{2}}$为拉普拉斯算子；$\mathrm{t}\mathrm{r}$为迹运算；$\otimes$为张量积运算。当$c_{\mathrm{2}}=\mathrm{0}$，$c_{\mathrm{1}}=c_{\mathrm{3}}=l^{\mathrm{2}}$时，式（3）就是在挠曲电效应分析中常见的拉普拉斯型梯度弹性本构方程：

$$\mathit{\sigma }=\lambda (\mathrm{t}\mathrm{r}\mathit{\epsilon })\mathit{I}+\mathrm{2}G\mathit{\epsilon }-l^{\mathrm{2}}{\nabla }^{\mathrm{2}}\left[\mathrm{2}G\mathit{\epsilon }+\lambda (\mathrm{t}\mathrm{r}\mathit{\epsilon })\mathit{I}\right]$$

（4）

式中：$l$即为材料单参数的内禀尺度。由于挠曲电效应中涉及应变梯度，要求插值函数至少C¹连续，传统的有限元方法采用拉格朗日插值往往仅具备C⁰连续。常见的应用于挠曲电效应分析的方法有混合有限元方法^［10-11］、无网格方法^［12-13］、等几何方法^［14-15］。混合有限元方法中，位移和位移导数是相互独立的变量，它们的运动学关系需要额外的自由度施加约束，从而使得格式较为复杂。对于无网格方法和等几何方法，它们不是插值类型的单元，因此不容易处理复杂的边界条件，在商用软件中也不容易扩展。C¹型有限单元开始用于平板的受弯分析，相关学者开发了一系列的C¹型有限单元^［16］，其中Argyris三角形单元也用于挠曲电效应的分析^［17］，但对于C¹型有限单元在挠曲电效应分析的研究还很少。

本文主要研究Bell三角形单元，它是由Argyris单元简化而来，其节点自由度包含位移及全部的一阶和二阶导数，Bell单元已经被用于许多梯度弹性问题的求解^［18-20］。基于梯度弹性理论，给出了Bell单元分析挠曲电效应的一般格式，并通过数值算例验证了Bell单元的准确性和收敛性，分析了内禀尺度对结构变形的影响。此外，利用挠曲电结构的非对称性设计压电结构，阐述了相关原理，并说明了挠曲电的尺寸效应。

1 挠曲电效应

不考虑压电效应，对于各向同性挠曲电材料，电焓密度$H$^［10］可以表示为

$$H=\frac{\mathrm{1}}{\mathrm{2}}\lambda {\epsilon }_{jj}{\epsilon }_{kk}+G{\epsilon }_{jk}{\epsilon }_{jk}+\frac{\mathrm{1}}{\mathrm{2}}{l}^{\mathrm{2}}(\lambda {\epsilon }_{jj,i}{\epsilon }_{kk,i}+\mathrm{2}G{\epsilon }_{jk,i}{\epsilon }_{jk,i})-f_{\mathrm{1}}{\epsilon }_{jj,i}{E}_i-\mathrm{2}{f}_{\mathrm{2}}{\epsilon }_{ij,i}{E}_j-\frac{\mathrm{1}}{\mathrm{2}}\kappa E_i{E}_i$$

（5）

式中：${\epsilon }_{ij}$为应变，定义为${\epsilon }_{ij}=\frac{\mathrm{1}}{\mathrm{2}}(u_{i,j}+u_{j,i})$，其中u为位移；$E_i$为电场强度，定义为$E_i=-{\varphi }_{,i}$，其中$\varphi$为电势；$f_{\mathrm{1}}$和$f_{\mathrm{2}}$为挠曲电材料两个相互独立的系数；$\kappa$为介电系数；$l$为材料内禀长度，它在梯度弹性理论（式（4））中引入。

本构方程可表示为

$${\sigma }_{jk}=\frac{\partial H}{\partial {\epsilon }_{jk}}=\lambda {\epsilon }_{ii}{\delta }_{jk}+\mathrm{2}G{\epsilon }_{jk}$$

（6）

$$\begin{array}{l}{\tau }_{ijk}=\frac{\partial H}{\partial {\eta }_{ijk}}=\frac{\partial H}{\partial {\epsilon }_{jk,i}}=l^{\mathrm{2}}(\lambda {\epsilon }_{mm,i}{\delta }_{jk}+\mathrm{2}G{\epsilon }_{jk,i})-\\ f_{\mathrm{1}}{E}_i{\delta }_{jk}-\mathrm{2}{f}_{\mathrm{2}}{E}_j{\delta }_{ik}\end{array}$$

（7）

$$D_i=\frac{\partial H}{\partial (-E_i)}=\kappa E_i+f_{\mathrm{1}}{\epsilon }_{mm,i}+\mathrm{2}{f}_{\mathrm{2}}{\epsilon }_{ji,j}$$

（8）

式中：${\sigma }_{jk}$为应力；${\tau }_{ijk}$为高阶应力；$D_i$为电位移；${\eta }_{ijk}$为高阶应变，定义为${\eta }_{ijk}={\epsilon }_{jk,i}$。

平衡方程可推得：

$${\sigma }_{jk,j}-{\tau }_{ijk,ij}+b_k=\mathrm{0}$$

（9）

$$D_{i,i}-\rho =\mathrm{0}$$

（10）

式中：$b_k$为体积力；$\rho$为体自由电荷。

2 C¹型三角形单元

Argyris三角形单元（图1）和Bell 三角形单元（图2）是两种重要的C¹型三角形单元。对于Argyris单元，单元自由度包括节点位移w及其一阶、二阶导数，以及三边中点法向导数。位移插值为五次函数，边界上的法向导数插值为四次函数。然而，边中节点法向导数自由度会显著增大单元刚度矩阵的带宽，这三个自由度可通过角节点的位移导数消去，即为Bell单元。对于Bell单元，位移的插值依然是五次函数，而边界上的法向导数插值变为三次函数。挠曲电效应是四阶偏微分方程问题，故Argyris和Bell单元均满足完备性和协调性要求，当划分单元足够小，有限元解趋近于精确解。

图1  Argyris三角形单元

Fig. 1  Argyris triangle element

图2  Bell三角形单元

Fig. 2  Bell triangle element

Bell三角形单元是基于插值函数的协调单元，便于建模和施加复杂的边界条件，也很容易在商用软件中扩展。 对于挠曲电效应分析，Bell单元每个节点有18个自由度，分别是位移、电势以及它们的一阶、二阶导数。

挠曲电效应弱形式可表达为

$${\int }_V(\mathit{\sigma }:\delta \mathit{\epsilon }+\mathit{\tau }⋮\delta \mathit{\eta }-\mathit{D}\cdot \delta \mathit{E})\mathrm{d}V={\int }_{V}b\cdot \mathrm{\delta }u\mathrm{d}V+{\int }_{{\mathrm{\Gamma }}_t}t\cdot \mathrm{\delta }u\mathrm{d}\mathrm{\Gamma }-{\int }_V\rho \mathrm{\delta }\varphi \mathrm{d}V-{\int }_{{\Gamma }_D}\omega \delta \varphi d\Gamma$$

（11）

式中：$\mathit{b}$、$\mathit{t}$、$\rho$和$\omega$分别为体积力、面积力、体自由电荷和表面电荷密度。

离散后，单元内部的位移场和电势场可表示为

$$\mathit{u}={\mathit{N}}_u{\mathit{v}}_{\mathrm{e}}$$

（12a）

$$\varphi ={\mathit{N}}_{\varphi }{\mathit{\phi }}_{\mathrm{e}}$$

（12b）

式（12）中：${\mathit{v}}_{\mathrm{e}}$和${\mathit{\phi }}_{\mathrm{e}}$分别为单元节点位移和电势；${\mathit{N}}_u$和${\mathit{N}}_{\varphi }$分别为单元位移和电势形函数，定义为

$${\mathit{N}}_u=\left[\begin{array}{cccccc}{N}_{\mathrm{1}}\cdots N_{\mathrm{6}}& \mathrm{0}& N_{\mathrm{7}}\cdots N_{\mathrm{12}}& \mathrm{0}& N_{\mathrm{13}}\cdots N_{\mathrm{18}}& \mathrm{0}\\ \mathrm{0}& N_{\mathrm{1}}\cdots N_{\mathrm{6}}& \mathrm{0}& N_{\mathrm{7}}\cdots N_{\mathrm{12}}& \mathrm{0}& N_{\mathrm{13}}\cdots N_{\mathrm{18}}\end{array}\right]$$

（13）

$${\mathit{N}}_{\varphi }=\left[\begin{array}{ccc}{N}_{\mathrm{1}}\cdots N_{\mathrm{6}}& N_{\mathrm{7}}\cdots N_{\mathrm{12}}& N_{\mathrm{13}}\cdots N_{\mathrm{18}}\end{array}\right]$$

（14）

式中：$N_{\mathrm{1}}~N_{\mathrm{18}}$取值为

$$N_{\mathrm{1}}=L_{\mathrm{1}}^{\mathrm{5}}+\mathrm{5}{L}_{\mathrm{1}}^{\mathrm{4}}{L}_{\mathrm{2}}+\mathrm{5}{L}_{\mathrm{1}}^{\mathrm{4}}{L}_{\mathrm{3}}+\mathrm{10}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{10}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{3}}^{\mathrm{2}}+\mathrm{20}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}+\mathrm{30}{r}_{\mathrm{21}}{L}_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}^{\mathrm{2}}+\mathrm{30}{r}_{\mathrm{31}}{L}_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}^{\mathrm{2}}{L}_{\mathrm{3}}$$

（15）

$$\begin{array}{l}{N}_{\mathrm{2}}=c_{\mathrm{3}}{L}_{\mathrm{1}}^{\mathrm{4}}{L}_{\mathrm{2}}-c_{\mathrm{2}}{L}_{\mathrm{1}}^{\mathrm{4}}{L}_{\mathrm{3}}+\mathrm{4}{c}_{\mathrm{3}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}^{\mathrm{2}}-\mathrm{4}{c}_{\mathrm{2}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{3}}^{\mathrm{2}}+\\ \mathrm{4}(c_{\mathrm{3}}-c_{\mathrm{2}})L_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}-(\mathrm{3}{c}_{\mathrm{1}}+\mathrm{15}{r}_{\mathrm{21}}{c}_{\mathrm{2}})L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}^{\mathrm{2}}+\\ (\mathrm{3}{c}_{\mathrm{1}}+\mathrm{15}{r}_{\mathrm{31}}{c}_{\mathrm{3}})L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}^{\mathrm{2}}{L}_{\mathrm{3}}\end{array}$$

（16）

$$\begin{array}{l}{N}_{\mathrm{3}}=-b_{\mathrm{3}}{L}_{\mathrm{1}}^{\mathrm{4}}{L}_{\mathrm{2}}+b_{\mathrm{2}}{L}_{\mathrm{1}}^{\mathrm{4}}{L}_{\mathrm{3}}-\mathrm{4}{b}_{\mathrm{3}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{4}{b}_{\mathrm{2}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{3}}^{\mathrm{2}}+\\ \mathrm{4}(b_{\mathrm{2}}-b_{\mathrm{3}})L_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}+(\mathrm{3}{b}_{\mathrm{1}}+\mathrm{15}{r}_{\mathrm{21}}{b}_{\mathrm{2}})L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}^{\mathrm{2}}-\\ (\mathrm{3}{b}_{\mathrm{1}}+\mathrm{15}{r}_{\mathrm{31}}{b}_{\mathrm{3}})L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}^{\mathrm{2}}{L}_{\mathrm{3}}\end{array}$$

（17）

$$N_{\mathrm{4}}=\frac{\mathrm{1}}{\mathrm{2}}{c}_{\mathrm{3}}^{\mathrm{2}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}{c}_{\mathrm{2}}^{\mathrm{2}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{3}}^{\mathrm{2}}-c_{\mathrm{2}}{c}_{\mathrm{3}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}+(c_{\mathrm{1}}{c}_{\mathrm{2}}+\frac{\mathrm{5}}{\mathrm{2}}{r}_{\mathrm{21}}{c}_{\mathrm{2}}^{\mathrm{2}})L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}^{\mathrm{2}}+(c_{\mathrm{1}}{c}_{\mathrm{3}}+\frac{\mathrm{5}}{\mathrm{2}}{r}_{\mathrm{31}}{c}_{\mathrm{3}}^{\mathrm{2}})L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}^{\mathrm{2}}{L}_{\mathrm{3}}$$

（18）

$$\begin{array}{l}{N}_{\mathrm{5}}=-b_{\mathrm{3}}{c}_{\mathrm{3}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}^{\mathrm{2}}-b_{\mathrm{2}}{c}_{\mathrm{2}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{3}}^{\mathrm{2}}+(b_{\mathrm{2}}{c}_{\mathrm{3}}+b_{\mathrm{3}}{c}_{\mathrm{2}})L_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}-\\ (b_{\mathrm{1}}{c}_{\mathrm{2}}+b_{\mathrm{2}}{c}_{\mathrm{1}}+\mathrm{5}{r}_{\mathrm{21}}{b}_{\mathrm{2}}{c}_{\mathrm{2}})L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}^{\mathrm{2}}-\\ (b_{\mathrm{1}}{c}_{\mathrm{3}}+b_{\mathrm{3}}{c}_{\mathrm{1}}+\mathrm{5}{r}_{\mathrm{31}}{b}_{\mathrm{3}}{c}_{\mathrm{3}})L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}^{\mathrm{2}}{L}_{\mathrm{3}}\end{array}$$

（19）

$$N_{\mathrm{6}}=\frac{\mathrm{1}}{\mathrm{2}}{b}_{\mathrm{3}}^{\mathrm{2}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}{b}_{\mathrm{2}}^{\mathrm{2}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{3}}^{\mathrm{2}}-b_{\mathrm{2}}{b}_{\mathrm{3}}{L}_{\mathrm{1}}^{\mathrm{3}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}+(b_{\mathrm{1}}{b}_{\mathrm{2}}+\frac{\mathrm{5}}{\mathrm{2}}{r}_{\mathrm{21}}{b}_{\mathrm{2}}^{\mathrm{2}})L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}{L}_{\mathrm{3}}^{\mathrm{2}}+(b_{\mathrm{1}}{b}_{\mathrm{3}}+\frac{\mathrm{5}}{\mathrm{2}}{r}_{\mathrm{31}}{b}_{\mathrm{3}}^{\mathrm{2}})L_{\mathrm{1}}^{\mathrm{2}}{L}_{\mathrm{2}}^{\mathrm{2}}{L}_{\mathrm{3}}$$

（20）

式中：$L_{\mathrm{1}}$、$L_{\mathrm{2}}$、$L_{\mathrm{3}}$为面积坐标；$b_i$、$c_i$、$r_{ij}$定义为

$$b_i=y_j-y_k$$

（21）

$$c_i=x_k-x_j$$

（22）

$$r_{ij}=\frac{b_i{b}_j+c_i{c}_j}{b_i^{\mathrm{2}}+c_i^{\mathrm{2}}}$$

（23）

式中：$x_j$和$y_j$为坐标值；i、j、k的取值采用指标轮换，即1→2，2→3，3→1。$N_{\mathrm{7}}~N_{\mathrm{18}}$的取值也采用同样的指标轮换。

应变$\mathit{\epsilon }$、高阶应变$\mathit{\eta }$及电场$\mathit{E}$的取值为$\mathit{\epsilon }=({\epsilon }_{\mathrm{11}}{\epsilon }_{\mathrm{22}}\mathrm{2}{\epsilon }_{\mathrm{12}}{)}^{\mathrm{T}}$，$\mathit{\eta }=({\eta }_{\mathrm{111}}{\eta }_{\mathrm{122}}\mathrm{2}{\eta }_{\mathrm{112}}{\eta }_{\mathrm{211}}{\eta }_{\mathrm{222}}\mathrm{2}{\eta }_{\mathrm{212}}{)}^{\mathrm{T}}$，$\mathit{E}=(E_{\mathrm{1}}{E}_{\mathrm{2}}{)}^{\mathrm{T}}$，可计算为

$$\begin{array}{l}\mathit{\epsilon }={\mathit{B}}_u{\mathit{v}}_{\mathrm{e}}(\mathrm{24}\mathrm{a})\\ \mathit{\eta }={\mathit{H}}_u{\mathit{v}}_{\mathrm{e}}(\mathrm{24}\mathrm{b})\\ -\mathit{E}={\mathit{B}}_{\varphi }{\mathit{\phi }}_{\mathrm{e}}\end{array}$$

（24c）

式中：${\mathit{B}}_u$、${\mathit{H}}_u$和${\mathit{B}}_{\varphi }$分别为位移形函数的梯度矩阵、海森矩阵和电势形函数的梯度矩阵。

将式（12）和（24）代入弱形式方程（11），可得如下单元节点平衡方程：

$$\left[\begin{array}{cc}{\mathit{K}}_{uu}& {\mathit{K}}_{u\varphi }\\ {\mathit{K}}_{\varphi u}& {\mathit{K}}_{\varphi \varphi }\end{array}\right]\left[\begin{array}{c}\mathit{v}\\ \mathit{\phi }\end{array}\right]=\left[\begin{array}{c}{\mathit{F}}_v\\ {\mathit{F}}_{\varphi }\end{array}\right]$$

（25）

式中：$\mathit{v}$和$\mathit{\phi }$分别为结构整体位移和电势；结构整体耦合刚度矩阵和节点载荷矩阵表达为

$${\mathit{K}}_{uu}=\sum _e{\int }_v({\mathit{B}}_u^{\mathrm{T}}\mathit{C}{\mathit{B}}_u+l^{\mathrm{2}}{\mathit{H}}_u^{\mathrm{T}}\mathit{Q}{\mathit{H}}_u)\mathrm{d}v$$

（26）

$${\mathit{K}}_{u\varphi }=\sum _e{\int }_v({\mathit{H}}_u^{\mathrm{T}}{\mathit{f}}^{\mathrm{T}}{\mathit{B}}_{\varphi })\mathrm{d}v$$

（27）

$${\mathit{K}}_{\varphi u}=\sum _e{\int }_v({\mathit{B}}_{\varphi }^{\mathrm{T}}\mathit{f}{\mathit{H}}_u)\mathrm{d}v$$

（28）

$${\mathit{K}}_{\varphi \varphi }=\sum _e{\int }_v(-{\mathit{B}}_{\varphi }^{\mathrm{T}}\mathit{\kappa }{\mathit{B}}_{\varphi })\mathrm{d}v$$

（29）

$${\mathit{F}}_v=\sum _e{\int }_v({\mathit{N}}_u^{\mathrm{T}}\mathit{b})\mathrm{d}v+\sum _e{\int }_{{\mathrm{\Gamma }}_t}({\mathit{N}}_u^{\mathrm{T}}\mathit{t})\mathrm{d}\mathrm{\Gamma }$$

（30）

$${\mathit{F}}_{\varphi }=-\sum _e{\int }_v({\mathit{N}}_{\varphi }^{\mathrm{T}}\rho )\mathrm{d}v-\sum _e{\int }_{{\mathrm{\Gamma }}_D}({\mathit{N}}_{\varphi }^{\mathrm{T}}\omega )\mathrm{d}\mathrm{\Gamma }$$

（31）

式中：材料参数$\mathit{C}$、$\mathit{Q}$、$\mathit{f}$、$\mathit{\kappa }$所对应的矩阵为

$$\mathit{C}=\left[\begin{array}{ccc}\lambda +\mathrm{2}G& \lambda & \mathrm{0}\\ \lambda & \lambda +\mathrm{2}G& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& G\end{array}\right]$$

（32）

$$\mathit{Q}=l^{\mathrm{2}}\left[\begin{array}{cc}\mathit{C}& \mathrm{0}\\ \mathrm{0}& \mathit{C}\end{array}\right]$$

（33）

$$\mathit{f}=\left[\begin{array}{cccccc}{f}_{\mathrm{1}}+\mathrm{2}{f}_{\mathrm{2}}& f_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& f_{\mathrm{2}}\\ \mathrm{0}& \mathrm{0}& f_{\mathrm{2}}& f_{\mathrm{1}}& f_{\mathrm{1}}+\mathrm{2}{f}_{\mathrm{2}}& \mathrm{0}\end{array}\right]$$

（34）

$$\mathit{\kappa }=\left[\begin{array}{cc}\kappa & \mathrm{0}\\ \mathrm{0}& \kappa \end{array}\right]$$

（35）

3 算例

Bell三角形单元通过Abaqus用户单元子程序UEL实现，采用TECPLOT进行结果的后处理。本文首先通过厚壁筒内外表面受压下变形的数值解与理论值的对比，验证Bell单元的准确性，讨论内禀尺度对结果的影响。然后，利用挠曲电结构的非对称性设计产生压电效应，并研究挠曲电的尺寸效应。

3.1　厚壁筒平面应变问题

图3为一厚壁筒，它的内径r₀为10 μm，外径r₁为20 μm，其受到内外表面电势差1.0 V作用。此外，内外径施加指定的位移载荷，内径位移u_r0为0.045 μm，外径位移u_r1为0.05 μm，挠曲电材料参数取值参照文献［10］，见表1。

图3  厚壁筒模型

Fig. 3  Model of cylinder

表1  挠曲电材料参数

Tab. 1  Material properties of flexoelectricity

$\lambda$ / GPa	G / GPa	l / μm	f1 / （μC·m-1）	f2 / （μC·m-1）	κ / （nC·V-1·m-1）
53.461 5	80.192 3	2	1.0	1.0	1.0

不考虑体自由电荷，将本构方程（6）—（8）代入平衡方程（9）、（10），可得：

$${\varphi }_{,ii}-\frac{f}{\kappa }{\nabla }^{\mathrm{2}}{u}_{i,i}=\mathrm{0}$$

（36）

$$(\lambda +G)(\mathrm{1}-l_{\mathrm{0}}^{\mathrm{2}}{\nabla }^{\mathrm{2}})u_{i,ik}+G(\mathrm{1}-l^{\mathrm{2}}{\nabla }^{\mathrm{2}})u_{k,ii}=\mathrm{0}$$

（37）

将方程（37）用极坐标系表示为

$$(\mathrm{1}-l_{\mathrm{0}}^{\mathrm{2}}{\nabla }^{\mathrm{2}}+\frac{l_{\mathrm{0}}^{\mathrm{2}}}{r^{\mathrm{2}}})({\nabla }^{\mathrm{2}}u(r)-\frac{u(r)}{r^{\mathrm{2}}})=\mathrm{0}$$

（38）

其中：

$$f=f_{\mathrm{1}}+\mathrm{2}{f}_{\mathrm{2}}$$

（39）

$$l_{\mathrm{0}}^{\mathrm{2}}=l^{\mathrm{2}}-\frac{f^{\mathrm{2}}}{(\lambda +G)\kappa }$$

（40）

由于厚壁筒结构的对称性，位移和电势只与半径有关，方程（38）的解可表示为^［10-11］

$$u(r)=C_{\mathrm{1}}r+C_{\mathrm{2}}\frac{\mathrm{1}}{r}+C_{\mathrm{3}}{I}_{\mathrm{1}}(\frac{r}{l_{\mathrm{0}}})+C_{\mathrm{4}}{K}_{\mathrm{1}}(\frac{r}{l_{\mathrm{0}}})$$

（41）

从而，电势的解可表示为

$$\varphi (r)=C_{\mathrm{5}}\mathrm{l}\mathrm{n}(r)+C_{\mathrm{6}}+\frac{f}{\kappa }(\frac{\mathrm{d}u(r)}{\mathrm{d}r}+\frac{u(r)}{r})$$

（42）

其中：$I_{\mathrm{1}}$和$K_{\mathrm{1}}$为一阶修正的第一类和第二类贝塞尔（Bessel）函数，方程（41）和（42）有$C_{\mathrm{1}}~C_{\mathrm{6}}$共6个未知数，可以通过如下6个边界条件求得：

$${\left.u_r\right|}_{r=r_{\mathrm{0}}}=u_{r_{\mathrm{0}}},{\left.u_r\right|}_{r=r_{\mathrm{1}}}=u_{r_{\mathrm{1}}}$$

（43）

$${\left.{\varphi }_r\right|}_{r=r_{\mathrm{0}}}=\mathrm{0},{\left.{\varphi }_r\right|}_{r=r_{\mathrm{1}}}=\mathrm{1}\mathrm{V}$$

（44）

$${\left.{\tau }_{rrr}\right|}_{r=r_{\mathrm{0}}}={\left.{\tau }_{rrr}\right|}_{r=r_{\mathrm{1}}}=\mathrm{0}$$

（45）

厚壁筒的切向位移和电势均为零，沿筒壁厚度方向径向位移的理论解以及不同网格密度的数值解见图4。对于比较稀疏的网格（环向48个单元，径向5个单元），数值解也与理论解吻合得很好。呈环状的位移云图也反映了位移仅与半径有关。图5为不同网格密度下，径向位移数值解与理论解误差绝对值的对数坐标图 （两端点由于施加位移约束，误差过小，因此略去）。总体来讲，随着网格密度的增加，数值解与理论解的误差减少。当网格比较稀疏，误差在对数坐标下的曲线呈现比较明显的振荡。径向网格数量为20和40的误差曲线大致重合。图6为筒壁厚度方向电势的理论解和数值解对比，稀疏网格会造成数值解与理论值略有偏差，当径向网格数量为20时，数值解与理论值位移曲线几乎重合。图7为不同网格密度下，径向电势数值解与理论解误差绝对值的对数坐标图（端点值略去）。随着网格的加密，误差同样减少，当径向网格数量加密到20时，误差曲线已较为光滑，数值解收敛。从图5和图7中可以看出，此时径向位移和电势数值解与理论解的误差绝对值分别小于10^-5 μm和10^-2 V。位移和电势的数值解与理论值对比说明了Bell三角形单元的准确性与解的收敛性。

图4  位移$u_r$的理论解和不同网格密度的数值解

Fig. 4  Theory solution and numerical solutions of different mesh grids of $u_r$

图5  位移$u_r$的理论解和不同网格密度的数值解误差绝对值

Fig. 5  Absolute value of $u_r$$$error between theory solution and numerical solutions of different mesh grids

图6  电势${\varphi }_r$的理论解和不同网格密度的数值解

Fig. 6  Theory solution and numerical solutions of different mesh grids of ${\varphi }_r$

图7  电势${\varphi }_r$的理论解和不同网格密度的数值解误差绝对值

Fig. 7  Absolute value of ${\varphi }_r$ error between theory solution and numerical solutions of different mesh grids

材料的内禀尺度$l$对径向位移的影响见图8，网格密度设为188×20，即环向188个单元，径向20个单元。从图中可以看出，数值解同样与理论解完全重合。此外，内禀尺度$l$越小，径向位移曲线弯曲程度越大，即径向应变（$\partial u_r/\partial r$）的梯度越大，这也可以从图9中更直观地看出。材料系数Q反映了高阶应力与高阶应变的关系，其表达式中含有内禀尺度的平方项，因此较大的内禀尺度会减小结构的应变梯度。在挠曲电材料计算中，内禀尺度大都直接设定，其取值方法需要进一步研究。

图8  不同内禀尺度下的位移$u_r$的解

Fig. 8  Solutions of $u_r$ for different intrinsic lengths

图9  不同内禀尺度下的应变梯度$u_r\mathrm{\text{'}}\mathrm{\text{'}}$

Fig. 9  Strain gradient $u_r\mathrm{\text{'}}\mathrm{\text{'}}$ for different intrinsic lengths

3.2　利用挠曲电材料结构的非对称性产生压电效应

相比压电材料，挠曲电效应存在于所有电介质中，而且不受居里温度的限制。挠曲电效应的一个重要应用就是利用挠曲电材料（非压电材料）设计结构产生压电效应^［21-22］。在正方形中心开孔，当受到均匀拉力时，结构内部由于应力分布不均匀就会产生应变梯度。但如果开孔是上下左右均对称，如圆形等，结构的平均电极化仍为零。因此，产生非零电极化需要利用开孔的不对称性，如图10所示的三角形结构。结构边长S为1 μm，开孔三角形的内接圆半径r为0.4 μm，左边位移和电势均为零，右边施加1.02×10^-4 N·μm^-1的线荷载F，此问题设为平面应变问题。

图10  三角形开孔的正方形结构

Fig. 10  Square structure with a triangular hole

不同内禀尺度下，结构右边电位移D₁（水平方向）如图11所示。由于结构的对称性，电位移曲线沿y = 0.5 μm轴对称。当内禀尺度较大时，电位移很小，而当内禀尺度较小时，电位移则比较显著，这是由于较大的内禀尺度减弱了结构的应变梯度。右侧电位移D₂（竖直方向）如图12所示，此时电位移曲线沿中轴反对称，电位移D₂沿边线的积分为零。由于结构是上下对称、左右非对称的，电位移才会出现上述的对称和反对称性，这也说明了开孔需是非对称的才可以产生压电效应。与图11类似，当内禀尺度较小时，结构中会产生更大的应变梯度，电位移D₂会更为显著。

图11  不同内禀尺度下的电位移D₁

Fig. 11  Electric displacement D₁ for different intrinsic lengths

图12  不同内禀尺度下的电位移D₂

Fig. 12  Electric displacement D₂ for different intrinsic lengths

挠曲电材料一个重要特点是尺寸效应，即挠曲电效应在小的尺度下尤为显著。对于许多材料，挠曲电系数很小，但在纳米尺度下，挠曲电效应非常强，而不能忽略。当逐渐缩小结构边长S，由微米尺度到纳米尺度（内禀尺度和荷载以相同比例缩小），结构的电位移D₁如图13所示。随着结构尺寸的减小，结构电位移随之增大。结构尺寸每减小10倍，电位移相应增大10倍左右。这一特点有利于微纳尺度下挠曲电材料结构设计。

图13  不同结构大小的电位移D₁

Fig. 13  Electric displacement D₁ for different structure sizes

4 结语

本文基于Bell三角形单元进行了挠曲电效应的分析。Bell三角形单元是C¹型协调单元，便于处理复杂边界条件，也很容易在商用软件中扩展。给出了在挠曲电分析中Bell单元的一般格式，并在数值算例中验证了它的准确性和收敛性。梯度弹性理论中引入的内禀尺度，可以用来预测微小尺度材料和结构中的尺度效应。较大的内禀尺度会减小结构的应变梯度，减弱正向挠曲电效应产生的电位移。此外，本文还利用结构的非对称性产生压电效应，通过电位移的分布阐述了相应的原理。挠曲电效应会随着结构尺度的减小显著增大，这有利于其在微纳尺度中的应用。

Bell三角形单元的一个劣势是其节点自由度较多，另外它扩展到三维单元还存在一定难度。本文采用了单参数的内禀尺度，其取值大都直接给出，此外还有学者采用多个内禀尺度参数进行挠曲电效应的分析^［23］。Bell单元的计算效率、扩展问题，以及内禀尺度的选取，都需要进一步研究。

作者贡献声明

庄晓莹： 学术指导，研究方法，撰写论文，论文修改。

李彬： 理论推导，数值计算，撰写论文。