随着现代科学技术的发展，对材料使用性能的要求越来越高，非均匀材料也越来越多地引起了国内外业界和学界的重视。在宏观上，非均匀材料可以通过连续地控制各组分含量的分布，使材料特性在空间位置上呈现梯度变化, 从而满足结构元件不同部位对材料特性的不同要求^[1]。对于非均匀电弹性材料，通过调整材料的组分可以使其力电特性按照优化设计进行分配，相应的元器件可以产生预期的变形和运动, 相比于普通的压电元器件有更显著的效果和更广泛的应用。

文献[2]中研究了本构非均匀的电弹性材料，这种材料的各物理参数由原来均匀情况下的常数变成为空间坐标的函数，相应其整体的物理特性也发生了变化。刘殿锋等^[3]通过改变压电相的体积百分比，研制表面振动为某些特殊函数分布的压电复合材料，进而制作出具有应用前景的非均匀振动换能器。文献[4]解得了功能梯度压电材料板条中移动裂纹的力电场和能量释放率。文献[5]研究了指数型功能梯度压电板条的裂纹问题。张亚芳等^[6]对在无穷远处受非均匀反平面剪应力场及非均匀平面电位移场作用下的内含螺旋位错和椭圆夹杂的压电材料作了理论分析，并给出了映射平面上压电材料基质和夹杂中的解析解。Zenkour等^[7-8]分别对非均匀压电空心圆柱体的湿-热-弹性响应、以及在热梯度作用下的压电行为进行了研究。

由于材料参数的非均匀性，其相应结构的控制方程通常为变系数的偏微分方程，解析求解难度大, 获得精确解非常困难^{[1, 9]}。板壳是工程结构的主要元件，对其自由振动精确解的研究不论在学术还是工程上都具有重要的意义。邢誉峰等^[10]利用空间坐标分离变量方法(非逆法)给出了板壳在简支和固支边界任意组合情况下自由振动的封闭(非级数)形式精确解，其中有一些精确解过去普遍认为不存在或难以求解的。仲政等^[11]基于状态空间列式给出了四边简支压电矩形板的精确解，伍晓红等^[12]利用幂级数展开法得到了压电功能梯度板自由振动的三维精确解公式，对四边简支压电功能梯度矩形板进行了分析，上述这两项工作都是假设材料参数沿厚度方向以同一指数形式变化。本文也采取相同的假设，给出了非均匀矩形电弹性板的基本方程。应用分离变量方法^[10]分别给出了四边简支和四边固支非均匀电弹性板的横向自由振动的精确解，这里所说的自由振动精确解指的是由本征微分方程得到的精确本征解。最后，通过算例讨论了相关问题。

1 非均匀电弹性矩形板的基本方程

在笛卡尔坐标系(x, y, z)中，非均匀电弹性矩形板的基本方程如下：

运动方程

$ \begin{array}{*{20}{c}} {\eta \left( k \right)\left( {D_1^E\frac{{{\partial ^2}{M_x}}}{{\partial {x^2}}} + 2D_3^E\frac{{{\partial ^2}{M_{xy}}}}{{\partial x\partial y}} + D_2^E\frac{{{\partial ^2}{M_y}}}{{\partial {y^2}}}} \right) + }\\ {\rho h\frac{{{\partial ^2}w}}{{\partial {t^2}}} = 0} \end{array} $

(1)

式中: ρ为密度, h为板厚，η(k)为非均匀性函数，k为不均匀系数，D^E为恒定电场下板的弯曲刚度，D^E=Y₀^Eh³/[12(1-υ²)], M_x、M_y为弯距, M_xy为扭矩, w为板的挠度，Y₀^E为恒定电场下的杨氏模量, υ为泊松比

$ \begin{array}{*{20}{c}} {\eta \left( k \right) = \frac{{12}}{{{h^3}}}\int_{ - h/2}^{h/2} {{z^2}{{\rm{e}}^{kz/h}}{\rm{d}}z} = }\\ {\frac{6}{{{k^3}}}\left[ { - 4k{\rm{ch}}\frac{k}{2} + \left( {8 + {k^2}} \right){\rm{sh}}\frac{k}{2}} \right],} \end{array} $

本构方程：

$ \begin{array}{*{20}{c}} {{M_x} = \int_{ - h/2}^{h/2} {{\sigma _{xx}}z{\rm{d}}z} = }\\ { - \eta \left( k \right)\left( {D_{11}^E\frac{{{\partial ^2}w}}{{\partial {x^2}}} + D_{12}^E\frac{{{\partial ^2}w}}{{\partial {y^2}}}} \right) - \int_{ - h/2}^{h/2} {z{d_{31}}{E_z}{{\rm{e}}^{kz/h}}{\rm{d}}z} }\\ {{M_y} = \int_{ - h/2}^{h/2} {{\sigma _{yy}}z{\rm{d}}z} = }\\ { - \eta \left( k \right)\left( {D_{12}^E\frac{{{\partial ^2}w}}{{\partial {x^2}}} + D_{22}^E\frac{{{\partial ^2}w}}{{\partial {y^2}}}} \right) - \int_{ - h/2}^{h/2} {z{d_{31}}{E_z}{{\rm{e}}^{kz/h}}{\rm{d}}z} }\\ {{M_{xy}} = \int_{ - h/2}^{h/2} {{\sigma _{xy}}z{\rm{d}}z} = - 2\eta \left( k \right)D_{66}^E\frac{{{\partial ^2}w}}{{\partial x\partial y}}} \end{array} $

(2)

式中：σ_xx和σ_yy为正应力, σ_xy为剪应力, E_z为电场强度, d₃₁为压电应变常数。

简支边界条件：

$ \begin{array}{*{20}{c}} {w\left( {0,y} \right) = 0,w\left( {a,y} \right) = 0,}\\ {w\left( {x,0} \right) = 0,w\left( {x,b} \right) = 0,}\\ {\frac{{{\partial ^2}w\left( {0,y} \right)}}{{\partial {x^2}}} = 0,\frac{{{\partial ^2}w\left( {a,y} \right)}}{{\partial {x^2}}} = 0,}\\ {\frac{{{\partial ^2}w\left( {x,0} \right)}}{{\partial {y^2}}} = 0,\frac{{{\partial ^2}w\left( {x,b} \right)}}{{\partial {y^2}}} = 0} \end{array} $

(3)

式中a和b为矩形板的边长。

固支边界条件：

$ \begin{array}{*{20}{c}} {w\left( {0,y} \right) = 0,w\left( {a,y} \right) = 0,}\\ {w\left( {x,0} \right) = 0,w\left( {x,b} \right) = 0,}\\ {\frac{{\partial w\left( {0,y} \right)}}{{\partial x}} = 0,\frac{{\partial w\left( {a,y} \right)}}{{\partial x}} = 0,}\\ {\frac{{\partial w\left( {x,0} \right)}}{{\partial y}} = 0,\frac{{\partial w\left( {x,b} \right)}}{{\partial y}} = 0} \end{array} $

(4)

考虑式(2)，则式(1)可以表示为

$ \eta \left( k \right)\left( D_{1}^{E}\frac{{{\partial }^{4}}w}{\partial {{x}^{4}}}+2D_{3}^{E}\frac{{{\partial }^{4}}w}{\partial {{x}^{2}}\partial {{y}^{2}}}+D_{2}^{E}\frac{{{\partial }^{4}}w}{\partial {{y}^{4}}} \right)+\rho h\frac{{{\partial }^{2}}w}{\partial {{t}^{2}}}=0 $

(5)

非均匀电弹性薄板自由振动或固有振动的解可写为

$ w = W\left( {x,y} \right)\left[ {A\cos \left( {\omega t} \right) + B\sin \left( {\omega t} \right)} \right] $

(6)

式中：W(x, y)为振型函数, ω为固有振动频率。

把式(6)代入式(5)，可得仅包含振型函数W(x, y)的本征微分方程为

$ D_1^E\frac{{{\partial ^4}w}}{{\partial {x^4}}} + 2D_3^E\frac{{{\partial ^4}w}}{{\partial {x^2}\partial {y^2}}} + D_2^E\frac{{{\partial ^4}w}}{{\partial {y^4}}} - W{\omega ^2}\rho h/\eta \left( k \right) = 0 $

(7)

2 非均匀电弹性板横向自由振动的精确解 2.1 四边简支

设满足边界条件式(3)的振型函数为

$ {W_{mn}} = \sin \frac{{m{\rm{ \mathsf{ π} }}x}}{a}\sin \frac{{n{\rm{ \mathsf{ π} }}y}}{b} $

(8)

式中m及n为整数。

将式(8)代入本征微分方程(7)，可得

$ \begin{array}{*{20}{c}} {\left[ {{{\rm{ \mathsf{ π} }}^4}\left( {D_1^E\frac{{{m^4}}}{{{a^4}}} + 2D_3^E\frac{{{m^2}{n^2}}}{{{a^2}{b^2}}} + D_2^E\frac{{{n^4}}}{{{b^4}}}} \right) - {\omega ^2}\rho h/\eta \left( k \right)} \right].}\\ {\sin \frac{{m{\rm{ \mathsf{ π} }}x}}{a}\sin \frac{{n{\rm{ \mathsf{ π} }}y}}{b} = 0} \end{array} $

(9)

式(9)对任意x和y都成立，因此有

$ {{\rm{ \mathsf{ π} }}^4}\left( {D_1^E\frac{{{m^4}}}{{{a^4}}} + 2D_3^E\frac{{{m^2}{n^2}}}{{{a^2}{b^2}}} + D_2^E\frac{{{n^4}}}{{{b^4}}}} \right) - {\omega ^2}\rho h/\eta \left( k \right) = 0 $

(10)

由此得固有振动频率

$ {\omega _{mn}} = {{\rm{ \mathsf{ π} }}^2}\sqrt {\frac{{\eta \left( k \right)\left( {D_1^E\frac{{{m^4}}}{{{a^4}}} + 2D_3^E\frac{{{m^2}{n^2}}}{{{a^2}{b^2}}} + D_2^E\frac{{{n^4}}}{{{b^4}}}} \right)}}{{\rho h}}} $

(11)

当非均匀系数k取不同值时，薄板的振动频率也不相同。

2.2 四边固支

设满足边界条件式(4)的振型函数为

$ W\left( {x,y} \right) = \phi \left( x \right)\psi \left( y \right) $

(12)

令ϕ(x)=e^μx，ψ(y)=e^λy。其中，变量μ和λ分别是与本征函数ϕ(x)和ψ(y)对应的本征值。将式(12)代入本征微分方程(7)，可得

$ D_1^E{\mu ^4} + 2D_3^E{\mu ^2}{\lambda ^2} + D_2^E{\lambda ^4} - {\omega ^2}\rho h/\eta \left( k \right) = 0 $

(13)

式(13)为式(7)的本征值代数方程。式(13)的解可表示为

$ \begin{array}{*{20}{c}} {{\mu _{1,2}} = \pm {\rm{i}}\sqrt {{\vartheta _1} + {\delta _1}} = \pm {\rm{i}}{\alpha _1},}\\ {{\mu _{3,4}} = \pm \sqrt {{\vartheta _1} - {\delta _1}} = \pm {\beta _1}} \end{array} $

(14)

式中：$ {\delta _1} = {\lambda ^2}\frac{{D_3^E}}{{D_1^E}},{\vartheta _1} = \sqrt {{\lambda ^4}\left[ {{{\left( {\frac{{D_3^E}}{{D_1^E}}} \right)}^2} - \frac{{D_2^E}}{{D_1^E}}} \right] + \frac{{{\omega ^2}\rho h}}{{D_1^E\eta (k)}}} $。

式(13)的解也可以表示为

$ \left\{ \begin{array}{l} {\lambda _{1,2}} = \pm {\rm{i}}\sqrt {{\vartheta _2} + {\delta _2}} = \pm {\rm{i}}{\alpha _2},\\ {\lambda _{3,4}} = \pm \sqrt {{\vartheta _2} - {\delta _2}} = \pm {\beta _2} \end{array} \right. $

(15)

式中$ {{\delta }_{2}}={{\mu }^{2}}\frac{D_{3}^{E}}{D_{2}^{E}} $,

$ {\vartheta _2} = \sqrt {{\mu ^4}\left[ {{{\left( {\frac{{D_3^E}}{{D_2^E}}} \right)}^2} - \frac{{D_1^E}}{{D_2^E}}} \right] + \frac{{{\omega ^2}\rho h}}{{D_2^E\eta \left( k \right)}}} $

把μ=iα₁代入式(15)，可得

$ \begin{array}{l} {\alpha _2} = \sqrt {\sqrt {\alpha _1^4\left[ {{{\left( {\frac{{D_3^E}}{{D_2^E}}} \right)}^2} - \frac{{D_1^E}}{{D_2^E}}} \right] + \frac{{{\omega ^2}\rho h}}{{D_2^E\eta \left( k \right)}}} - \alpha _1^2\frac{{D_3^E}}{{D_2^E}}} \\ {\beta _2} = \sqrt {\sqrt {\alpha _1^4\left[ {{{\left( {\frac{{D_3^E}}{{D_2^E}}} \right)}^2} - \frac{{D_1^E}}{{D_2^E}}} \right] + \frac{{{\omega ^2}\rho h}}{{D_2^E\eta \left( k \right)}}} - \alpha _1^2\frac{{D_3^E}}{{D_2^E}}} \end{array} $

(16)

从式(14)中消去λ可得

$ {\left( {\alpha _1^2 + \beta _1^2} \right)^2} + \frac{{D_1^ED_2^E - {{\left( {D_3^E} \right)}^2}}}{{{{\left( {D_3^E} \right)}^2}}}{\left( {\alpha _1^2 - \beta _1^2} \right)^2} = \frac{{4{\omega ^2}\rho h}}{{D_1^E\eta \left( k \right)}} $

(17)

根据式(14)和(15)的本征根，式(12)中的两个本征函数的解可写为

$ \begin{array}{*{20}{c}} {\varphi \left( x \right) = {A_1}\cos \left( {{\alpha _1}x} \right) + {B_1}\sin \left( {{\alpha _1}x} \right) + }\\ {{C_1}\cosh \left( {{\beta _1}x} \right) + {H_1}\sinh \left( {{\beta _1}x} \right)}\\ {\psi \left( y \right) = {A_2}\cos \left( {{\alpha _2}y} \right) + {B_2}\sin \left( {{\alpha _2}y} \right) + }\\ {{C_2}\cosh \left( {{\beta _2}y} \right) + {H_2}\sinh \left( {{\beta _2}y} \right)} \end{array} $

(18)

由边界条件式(4)可得

$ \begin{array}{*{20}{c}} {\phi \left( 0 \right) = 0,\phi \left( a \right) = 0,\phi '\left( 0 \right) = 0,\phi '\left( a \right) = 0}\\ {\psi \left( 0 \right) = 0,\psi \left( b \right) = 0,\psi '\left( 0 \right) = 0,\psi '\left( b \right) = 0} \end{array} $

(19)

将式(19)代入式(18)，分别令两个系数行列式为零，经简化可得两个本征方程为

$ \begin{array}{*{20}{c}} {\frac{{1 - \cos \left( {{\alpha _1}a} \right)\cosh \left( {{\beta _1}a} \right)}}{{\sin \left( {{\alpha _1}a} \right)\sinh \left( {{\beta _1}a} \right)}} = \frac{{\alpha _1^2 - \beta _1^2}}{{2{\alpha _1}{\beta _1}}}}\\ {\frac{{1 - \cos \left( {{\alpha _2}b} \right)\cosh \left( {{\beta _2}b} \right)}}{{\sin \left( {{\alpha _2}b} \right)\sinh \left( {{\beta _2}b} \right)}} = \frac{{\alpha _2^2 - \beta _2^2}}{{2{\alpha _2}{\beta _2}}}} \end{array} $

(20)

3 非均匀电弹性板横向自由振动数值分析

非均匀电弹性矩形薄板厚度h=0.02 m，维度a×b=1 m×1.2 m。BaTiO₃材料参数参照文献[13]。根据式(16)、(17)和(20)，采用Newton-Raphson方法可以求得相应的频率和空间本征值。

对于四边简支和四边固支的情况，不均匀系数与频率的关系分别如图 1和图 2所示。其中，不均匀系数的绝对值大小反映了材料本身不均匀程度的强弱，正负号只代表材料性能变化的方向。

由图 1和图 2可看出，在简支和固支的情况下，当不均匀系数取正数时，随着数值的增大，频率也随之增大，频率增速由频率的阶数决定，阶数越高，频率增长的越快；当不均匀系数取负数时，随着数值的减小，频率随之增大，频率增速由频率阶数决定，阶数越高，频率增长的越快。在相同的频率阶数、相同的不均匀系数的情况下，固支边界的固有频率远高于简支边界的固有频率。

当k=0时，η(0)=1，非均匀电弹性矩形薄板退化为均匀电弹性矩形薄板，频率可以表示为

$ {\omega _{mn}} = {{\rm{ \mathsf{ π} }}^2}\sqrt {\frac{{D_1^E\frac{{{m^4}}}{{{a^4}}} + 2D_3^E\frac{{{m^2}{n^2}}}{{{a^2}{b^2}}} + D_2^E\frac{{{n^4}}}{{{b^4}}}}}{{\rho h}}} $

(21)

根据文献[12]中材料1的参数和板的尺寸，采用上述方法，可得均匀电弹性板横向自由振动的精确解，如表 1所示。

由表 1可以看出，四边简支均匀电弹性板自由振动前6阶频率与文献[12]的结果相差很小, 可以证明本文所给精确解的正确性。

4 结论

1) 分析了非均匀电弹性矩形板的不均匀系数与振动频率之间的变化规律。在简支和固支的情况下，当不均匀系数取正数时，随着数值的增大，频率也随之增大，频率增速由频率阶数决定，阶数越高，频率增长的越快；当不均匀系数取负数时，随着数值的减小，频率随之增大，频率增速由频率阶数决定，阶数越高，频率增长的越快。在相同的阶数、相同的不均匀系数的情况下，固支边界的固有频率远大于简支边界的固有频率。

2) 将非均匀电弹性板的横向自由振动的精确解退化为均匀电弹性板横向自由振动的精确解，并验证了其精确解的正确性。这也说明了均匀电弹性板横向自由振动的精确解是非均匀电弹性板的横向自由振动的精确解的特例。同样，也可以应用推广与退化原理，将均匀电弹性板的横向自由振动的精确解推广为非均匀电弹性板横向自由振动的精确解。