-
开放科学（资源服务）标识码（OSID）：
-
功能梯度材料广泛用于航空航天、生物医学工程、核工业、土木与机械工程等领域。功能梯度板作为一类特殊的线弹性板，其特殊之处在于材料参数(如杨氏模量)在厚度方向呈梯度变化。适用于线弹性板的求解理论在满足一定条件时也可应用于功能梯度板的求解。许多学者对功能梯度板的力学性能展开了研究。解决这类问题的方法主要基于特定的剪切变形理论、有限元法及其他特定理论。
文献[1]基于von Karman非线性运动学，采用一阶剪切变形理论，通过牛顿-拉夫森迭代技术提出了一种非线性有限元方法，研究了多向多孔材料功能梯度矩形板的弹塑性特性。文献[2]提出了一种广义单变量剪切变形板理论，用于分析多孔功能梯度矩形板的自由振动。文献[3]求解了一种具有通厚孔隙率的智能层复合材料和功能梯度材料矩形板的静态弯曲问题。文献[4]用高阶剪切法和法向变形板理论分析了不可压缩矩形功能梯度板的弯曲、屈曲和自由振动。文献[5-10]基于各种剪切变形理论提出了适用于分析不同状况的功能梯度矩形板的弯曲、自由振动和屈曲的无网格方法。文献[11]采用无单元kp-Ritz方法对具有任意边界条件的功能梯度矩形板在机械和热载荷作用下进行了静力分析。文献[12]提出了一种基于边缘的光滑有限元法，用于功能梯度材料正方形板的静态、自由振动和屈曲分析。文献[13]提出了基于节点的应变平滑的一种改进的有限元方法，应用于分析功能梯度正方形板的静态、自由振动和机械热屈曲问题。文献[14]提出了一种基于等几何方法和高阶变形板理论的分析方法，研究了功能梯度材料矩形与圆形板的力学行为。文献[15]提出了一种改进的无网格径向点插值方法，分析了功能梯度矩形板的非线性弯曲问题。文献[16]基于Levinson板块理论和一阶剪切变形板块理论对功能梯度厚圆扇形板进行了弯曲分析。文献[17]用正弦剪切变形理论分析了功能梯度陶瓷-金属矩形板的弯曲问题。文献[18]基于修正的耦合应力理论、Kirchhoff/Mindlin板理论以及von Karman几何非线性理论研究了弹性基础上简支功能梯度微孔矩形板的非线性弯曲和自由振动响应问题。文献[19]提出了一种简单的一阶剪切变形理论，用于分析功能梯度矩形板的弯曲和自由振动。文献[20]基于无网格自然单元法建立了功能梯度板的一阶剪切变形理论的自由振动分析格式，并通过矩形功能梯度板的基频优化了实例验证算法的可行性和有效性。文献[21]将能量辐射传递法推广到矩形功能梯度板模型中，预测结构的高频振动响应。文献[22]提出了基于S-R理论和分解定理的无网格法，用于求解矩形功能梯度板的非线性弯曲问题。文献[23]提出了一种简单精化板理论，利用此理论分析矩形功能梯度板的弯曲问题时只需3个未知量，而分析矩形功能梯度板的自由振动问题时只需1个未知量。文献[24]基于一阶剪切变形板理论提出了求解功能梯度板自由振动问题的Cel-Based光滑有限元方法。
总之，学者们对功能梯度板的研究成果主要集中在对功能梯度矩形板、圆扇形板、圆形板的研究，很少有对凹凸不规则复杂形状功能梯度板的研究。本研究基于一阶剪切变形理论，针对复杂形状功能梯度板弯曲问题展开探索。引入R-函数理论[25]，借助其隐函数形式构建方程，可精准构造出复杂形状板的边界条件，能够便捷地描绘含缺角、凹凸及非对称等特征的一般性复杂形状，在形状描述的灵活性方面更优，为凹凸不规则复杂形状功能梯度板的分析提供新途径。通过引入R-函数理论，可以方便地构造具有不规则凹凸形状的功能梯度板的挠度函数ω，确保这些函数满足边界条件。R-函数理论可以用隐式函数形式描述复杂区域，通过归一化方程ω=0表示问题边界，并通过不等式ω＞0定义区域本身。随后，使用变分法求解这些不规则形状的功能梯度板的挠度函数。
本研究将R-函数理论运用到功能梯度板的弯曲问题中，以隐函数形式构造了满足复杂形状板边界条件的方程，通过对功能梯度板弯曲总能量方程的变分，求解出了复杂形状功能梯度板弯曲问题的解，给出了不规则形状功能梯度板的详细数值分析案例。使用MATLAB软件求解变分公式，并将所得结果与有限元分析数据进行比较，以验证所提出方法的可靠性。
HTML
-
研究对象是板厚恒为h的功能梯度板。功能梯度板由陶瓷和钢复合而成，上表面为陶瓷，下表面为钢，其弹性模量符合幂律函数，则功能梯度板的弹性模量E(z)表示为该板的弹性参数随厚度变化的律分布，其值为：
式中：
$k(k \geqslant 0)$ 表示梯度指数；下标中的$c$ 和$s$ 分别代表陶瓷和钢，$E_{c}$ 是陶瓷的弹性模量，$E_{s}$ 是钢的弹性模量；$z$ 是厚度方向坐标$\left(-\frac{h}{2} \leqslant z \leqslant \frac{h}{2}\right)$ 。
-
在修正的一阶剪切变形理论中，横向位移w分为弯曲项和剪切项
式中：wb为弯曲项；ws为剪切项。
假设：
式中：下标x和y分别代表对弯曲项x和y的偏导数。得到以下位移场：
式中：u、v、w分别为x、y、z方向位移；u0、v0分别是中面(z=0)处x、y方向位移；wb为弯曲项位移，ws为剪切项位移。
由(4)-(6)式得：
式中：εx、εy分别为板面内沿x、y轴正应变；γxy为板面内剪应变，γxz、γyz为横向剪应变。
当考虑一阶剪切变形理论时，功能梯度板的线性本构方程为：
式中：
$\sigma_{x} 、\sigma_{y}$ 分别为板面内沿$x 、y$ 轴正应力；$\tau_{x y}$ 为板面内剪应力，$\tau_{y z} 、\tau_{x z}$ 为横向剪应力；$\boldsymbol{\sigma}=\left(\sigma_{x}, \sigma_{y}\right.$ ，$\left.\tau_{y z}, \tau_{x z}, \tau_{x y}\right)^{\mathrm{T}}$ 和$\boldsymbol{\varepsilon}=\left(\varepsilon_{x}, \varepsilon_{y}, \gamma_{y z}, \gamma_{x z}, \gamma_{x y}\right)^{\mathrm{T}}$ 分别是应力矢量和应变矢量；刚度系数$Q_{i j}$ 表示为：式中：μ为泊松比。
在直角坐标系oxyz下，整个弹性体势能为：
式中：Vε为弹性体势能；积分域V为功能梯度板所占空间区域。
由于剪切一阶变形理论是二维剪切变形理论，当考虑二维剪切变形理论时，εz=0，即：
对于(14)式，若为厚板或中厚板，则需变为：
式中：K是与板厚度相关的修正系数。
对于功能梯度板的应变能计算，需将厚度方向离散化。将板的厚度范围
$\left[-\frac{h}{2}, \frac{h}{2}\right]$ 离散化为$n$ 个小的厚度区间，每个区间的厚度增量为$\Delta z=\frac{h}{n}$ 。定义厚度离散点$z_{i}=-\frac{h}{2}+i \times \Delta z, i=0 、1 、2 、\cdots 、n$ 。计算每个厚度区间的应变能贡献：对于每个厚度区间$\left[z_{i}, z_{i+1}\right]$ ，根据$E(z)$ 随厚度变化的公式计算该区间中点
$z_{\text {mid }}=\frac{z_{i}+z_{i+1}}{2}$ 处的$E(z)$ 值。计算该厚度区间对应板平面内的应变能贡献：
式中：dz在该区间近似于Δz。得到应变能：
对所有厚度区间的应变能贡献进行求和，得到近似的总应变能：
将(7-8)式代入(19)式并利用广义胡克定律得：
式中：
$E\left(z_{\text {mid }}\right)$ 为厚度区间中点处的弹性模量。现在把
$u 、v 、w$ 的表达式设定为：式中：Am、Bm、Cm为待定系数；um、vm为设定函数；Cm为互不依赖的m个待定系数；wm为满足板位移边界条件(即约束条件)的设定函数，wbm、wsm分别为弯曲项设定函数和剪切项设定函数，wm=wbm+wsm，当研究对象为功能梯度薄板时，不考虑厚度拉伸效应，且横向位移中只考虑弯曲项，则可假设wm=wbm。这样，不论Cm如何取值，(21)式所示的挠度w总能满足位移边界条件。注意：挠度w的变分只是由系数Cm的变分来实现的；设定函数wm仅随坐标而变，与上述变分完全无关。
为了确定系数Cm，需应用下列式子：
式中：fz为方向z的体积力集度；fz为边界面力集度；S是表面积分区域。
在薄板的弯曲问题中，体力及面力都归入荷载q。由此可以得出Cm的m个线性方程用来确定Cm，从而由(20)式得出挠度w，求得板的内力。
根据最小势能原理，总势能取极值的条件为：
(23) 式又可以写成：
式中：
$F_{x} 、F_{y} 、F_{z}$ 分别为$x 、y 、z$ 方向集力；$\bar{f}_{x} 、\bar{f}_{y}$ 为$x 、y$ 方向边界面力集度。
-
根据
$R$ -函数理论，一个复杂区域$\omega_{0}$ 的各个子域可通过布尔操作$∨_{\alpha} 、\wedge_{\alpha}$ 来实现$\left(∨_{\alpha} 、\wedge_{\alpha}\right.$ 分别为求并集与求交集)。若各子域的函数$\omega_{l}$ 满足一阶规范化方程，则$\omega_{0}$ 也满足一阶规范化方程。一阶规范化方程定义为：
式中：
$\varOmega$ 为复杂形状区域；$\partial \varOmega$ 为区域边界线。设
$X$ 与$Y$ 满足一阶规范化方程，则布尔操作$X \vee_{\alpha} Y 、X \wedge_{\alpha} Y$ 也满足一阶规范化方程：式中：
$\alpha$ 是布尔操作参数$(-1 \leqslant \alpha \leqslant 1)$ 。
-
挠度的表达式为：
在弯曲问题中，将R-函数加入公式(31)，讨论固支边的情况，固支边挠度为0，转角为0。参考R-函数的性质，将得到的ω0进行平方使其满足固支边的边界条件。
将(31)式代入(20)式，得：
进一步对Ci求导(i=1、2、3、…、m-1、m)，可得：
将(32)-(40)式代入(22)式，化简之后得：
式中：
$\sigma_{1}=\frac{\partial}{\partial x} w_{b j}(x, y) ; \sigma_{2}=\frac{\partial}{\partial y} w_{b j}(x, y) ; \sigma_{3}=2\left(\mu^{2}-1\right) ; \sigma_{4}=\frac{\partial^{2}}{\partial y^{2}} w_{b j}(x, y) ; \sigma_{5}=\frac{\partial^{2}}{\partial x^{2}} w_{b j}(x, y)$ ；$\sigma_{6}=w_{b i}(x, y) \sigma_{8}+\sigma_{11} \frac{\partial^{2}}{\partial y^{2}} w_{b i}(x, y)+4 \sigma_{13} w_{0}(x, y) \sigma_{9} ; \sigma_{7}=w_{b i}(x, y) \sigma_{10}+\sigma_{11} \frac{\partial^{2}}{\partial x^{2}} w_{b i}(x, y)+ 4 \sigma_{14} w_{0}(x, y) \sigma_{12} $ ；$ \sigma_{8}=2 w_{0}(x, y) \frac{\partial^{2}}{\partial y^{2}} w_{0}(x, y)+2 \sigma_{13}{ }^{2} ; \sigma_{9}=\frac{\partial}{\partial y} w_{b i}(x, y) ; \sigma_{10}=2 w_{0}(x, y) \frac{\partial^{2}}{\partial x^{2}} w_{0}(x$ ，$y)+2 \sigma_{14}^{2} ; \sigma_{11}=w_{0}(x, y)^{2} ; \sigma_{12}=\frac{\partial}{\partial x} w_{b i}(x, y) ; \sigma_{13}=\frac{\partial}{\partial y} w_{0}(x, y)$ ；$\sigma_{14}=\frac{\partial}{\partial x} w_{0}(x, y)$ 。经推导可得系数方程组：
式中：
$\boldsymbol{A}=\left(a_{i j}\right)_{m \times m}$ 为刚度矩阵；$\boldsymbol{B}=\left(b_{i j}\right)_{m \times 1}$ 为荷载向量。式中：
$i=1 、2 、\cdots 、m, j=1 、2 、\cdots 、m ; \varOmega$ 为笛卡尔坐标系中横截面($x \circ y$ )所包含的区域。
-
例1(功能梯度薄板) 设U形功能梯度薄板泊松比取0.3，板厚度h取0.01 m，均布荷载大小q为5×103 Pa，如图 1所示，其中a=b=0.5 m，c=d=0.25 m。上表面为陶瓷，下表面为钢，其弹性模量符合幂律函数，则功能梯度薄板的弹性模量E(z)表示为该板的弹性参数随厚度变化的律分布，如公式(1)所示。陶瓷的弹性模量Ec取210 GPa，钢的弹性模量Es取390 GPa。功能梯度材料的组成材料确定后，板的弹性模量取决于梯度指数。随着梯度指数的增大，板的材料中陶瓷成分的比例逐渐减少而钢成分的比例逐渐增多。当梯度指数增大到无穷时，板的材料近似于纯钢板，从而最大挠度也近似于钢板的最大挠度。在此算例中分别取k=0、k=1、k=10和k=+∞。
令
$n=20$ ，将板的厚度范围$\left[-\frac{h}{2}, \frac{h}{2}\right]$ 离散化为20个小的厚度区间，每个区间的厚度增量为$\Delta z=\frac{h}{n}$ ，定义厚度离散点$z_{i}=-\frac{h}{2}+i \times \Delta z(i=0 、1 、2 、\cdots 、n)$ 。计算每个厚度区间的应变能贡献：对于每个厚度区间$\left[z_{i}, z_{i+1}\right]$ ，计算该区间中点$z_{\text {mid }}=\frac{z_{i}+z_{i+1}}{2}$ 处的$E(z)$ 值，代入(42)式即可得系数方程组。根据R-函数理论，U型功能梯度薄板的ω0表达式为：
式中：
$\omega_1=\frac{a^2-x^2}{2 a} \geqslant 0 ; \omega_2=\frac{b^2-y^2}{2 b} \geqslant 0 ; \omega_3=(c-x) \geqslant 0 ; \omega_4=\frac{y^2-d^2}{2 d} \geqslant 0$ 。当研究对象为功能梯度薄板时不考虑厚度拉伸效应，且横向位移只考虑弯曲项。假设wm=wbm，ws=0，则设定函数wm取1、x、y2、x3、y4、x5、y6、y8、xy4、x3y2。
数值计算中将U形功能梯度薄板划分为3个矩形，第一个矩形板划分为4N×N网格，第二个矩形板划分为3N×2N网格，第三个矩形板划分为4N×N网格。
当k=0和k=+∞，m=10，N取不同值时，中间点挠度结果如表 1。当k=0，N=12，m=10时，本文方法与有限元法计算y=0上的挠度曲线结果对比如图 2。当k=+∞，m取不同值时，计算y=0上的挠度曲线如图 3。说明了本文方法是收敛的、正确的。当取k=0、k=1、k=10和k=+∞，试函数数量m=10，网格数取N=12时，y=0上的挠度曲线如图 4。变化趋势显示，在其他条件相同的情况下，k值越大，挠度越小，验证了本文方法的可行性与正确性。
例2(功能梯度中厚板) 设矩形功能梯度中厚板的边长分别为2a和2b，a=b=0.5 m，泊松比取0.3，板厚度h取0.12 m，均布荷载大小q为2×107 Pa，如图 5所示。鉴于当前板的宽度为1 m，厚度达0.12 m，故在相关计算中需要引入剪切修正因子K，取值为
$\frac{5}{6}$ 。其他参数取值如例1。矩形功能梯度中厚板的上表面为陶瓷，下表面为钢，其弹性模量符合幂律函数，则功能梯度中厚板的弹性模量E(z)表示为该板的弹性参数随厚度变化的律分布，如公式(1)。
根据R-函数理论，取：
式中：
$\omega_{1}=\frac{a^{2}-x^{2}}{2 a} \geqslant 0, \omega_{2}=\frac{b^{2}-y^{2}}{2 b} \geqslant 0$ 。将厚度离散化为30个小区间。假设wm=wbm+wsm。此板的形状为对称图形，设定挠度w的函数为x及y的偶函数，所以wbm取7、7x2、7y2、7x4、7y4、7x2y2、7x6、7y6、7x2y4、7x4y2。wsm取1、x4、y4、x8、y8、x4y4、x12、y12、x4y8、x8y4。
数值计算中将矩形区域划分成N×N网格。当k=0和k=+∞，试函数数量m=10，网格数分别取10×10、15×15、20×20、35×35时，板中心点挠度结果如表 2。由表 2可知本文方法计算的功能梯度中厚板中心点挠度是有效的。当k=+∞，m=10时，本文方法与有限元法计算y=0上的挠度曲线吻合得很好(图 6)。由此证明了本文方法的可行性与正确性。当k=0，m取不同值时，功能梯度板y=0上的挠度曲线如图 7。当取k=0、k=1、k=10和k=+∞，试函数数量m=10，网格数取35×35时y=0上的挠度曲线如图 8。
例3(功能梯度厚板) 设功能梯度厚板的泊松比取0.3，板厚度h取0.25 m，均布荷载大小q为2×108 Pa，如图 9所示。a=b=0.5 m。鉴于当前板的宽度为1 m，厚度0.25 m，故在相关计算中需要引入剪切修正因子K，取值为
$\frac{5}{6}$ 。其他参数取值如例1。矩形功能梯度厚板的上表面为陶瓷，下表面为钢，其弹性模量符合幂律函数，则功能梯度厚板的弹性模量E(z)表示为该板的弹性参数随厚度变化的律分布，如公式(1)。
根据R-函数理论，取ω0为：
式中：
$\omega_{1}=\frac{a^{2}-x^{2}}{2 a} \geqslant 0 ; \omega_{2}=\frac{b^{2}-y^{2}}{2 b} \geqslant 0 ; \omega_{3}=\frac{(x-a)^{2}+y^{2}-r^{2}}{2 r} \geqslant 0$ 。将厚度离散化为50个小区间。此处板的形状为不对称图形，设定挠度
$w$ 的函数为$x$ 及$y$ 的偶函数和奇函数，所以$w_{b m}$ 取$1 、x^{2} 、y^{2} 、x 、y^{4} 、x^{3} 、y^{6} 、x y^{2} 、x^{2} y^{2} 、x^{2} y^{4}, ~ w_{s m}$ 取$1 、0.3 x^{12} 、0.3 y^{12}$ 、$0.3 x^{11}, 0.3 y^{14}, 0.3 x^{13}, 0.3 y^{16}, 0.3 x^{11} y^{2}, 0.3 x^{12} y^{2}, 0.3 x^{2} y^{14}$ 。将矩形区域划分成N×N网格，获取所有坐标后，去除半圆形所包含的坐标。当k=0和k=+∞，试函数数量m=10，网格数分别取10×10、15×15、35×35时，板中心点挠度结果如表 3，由表 3可知本文方法计算的功能梯度厚板中心点挠度是有效的。当k=0，m=10时，有限元法与本文方法计算y=0上的挠度曲线对比如图 10。由此证明了本文方法的可行性与正确性。当k=0，网格数为35×35，m取不同值时，功能梯度板y=0上的挠度曲线如图 11。当取k=0、k=1、k=10和k=+∞，试函数数量m=10，网格数取35×35时，y=0上的挠度曲线如图 12。
R-函数理论可精确描述复杂几何形状的边界条件，将不规则形状转化为数学表达式，极大地简化了数值计算过程并保障了几何信息的准确性。而变分法基于能量最小化原理，把挠度问题转化为求解能量泛函最小值的问题。R-函数理论与变分法的结合不仅能够妥善处理复杂的几何边界条件，还能确保计算结果具备高精度，在处理不同形状以及不同厚度的功能梯度板的挠度计算时展现出良好的有效性。在试函数数量和网格划分精细程度方面，本文方法对计算结果有着重要影响。对于试函数，本文方法用于描述功能梯度板的可能变形模式。当试函数数量增加时，能更精细地拟合实际的变形状况。例如在考虑功能梯度板的剪切变形时，添加包含剪切变形相关项的试函数，便可更精准地描述板在受载下的实际变形，进而提高了计算结果的精度。同样，网格划分越精细，计算结果越接近有限元解。与有限元法的结果的一致性进一步证实了本文方法在工程应用中的价值。
在计算实例中，本文方法与有限元法计算的结果存在一定偏差，可能由以下因素引起：
1) 本文选定的试函数可能在非边界区域意外地符合一些本不存在的边界条件，这种情形可能会干扰研究的计算结果。以x2为例，该函数在x=0的直线上恒为0，这显然与考虑的情况相悖。
2) 本文的探讨中可能存在一个局限性，即在选取试验函数时可能未能纳入充足的数量。在分析中，可以观察到当m值偏低，意味着试验函数的数量较少时，最终得出的计算结果与有限元法的结果相比，存在偏差。这一发现说明，为了通过变分法获得更为精确的解，必须尽可能地增加试验函数的数量，以确保结果的可靠性和精确性。
3) 在对R-函数进行积分操作时，可能会遇到分母在边界处为0的特殊情况。因此，必须采用类似于有限元分析的方法，通过划分网格并逐个求解各个网格单元的积分，然后累加以获得整体解。这一过程强调了网格划分的精细度对于最终结果准确性的重要性。在应用R-函数理论与变分原理相结合的数值求解策略时，为了提高解的精度，应当在计算过程中采用尽可能密集的网格划分。
4) 在本文的讨论过程中，可能存在关于剪切修正系数的局限性。由于剪切修正系数的局限性，在某些特定情况，如特定的板厚、荷载以及材料分布组合下，可能会导致预测结果的精度不足，这在一定程度上限制了该剪切修正方法在厚板精确分析中的应用。尤其是当工程设计与分析对精度有较高要求时，可能需要采用更精确的剪切修正方法，或考虑采用可变剪切修正系数的理论来满足需求。为有效解决这一问题，可针对不同厚度的厚板开展实验与模拟分析，以确定相应的剪切修正系数，从而提高挠度计算的精度。鉴于此，未来的一个重要研究方向是基于厚板的厚度差异构建剪切修正系数模型，以便更精准地预测不同厚度厚板的挠度行为，为工程实践提供更可靠的理论支撑。
-
变分法被广泛应用于解决工程和物理问题，特别是在计算功能梯度材料的变形和应力分布方面。然而，当遇到形状复杂的功能梯度板时，传统的变分法可能会遇到一些问题。这时，引入R-函数的概念，可以作为一种补充方法来处理边界条件的复杂性。
R-函数理论能够在不直接定义边界条件的情况下，通过隐函数来描述复杂的几何形状。将R-函数理论与变分法结合使用，可以有效地简化复杂边界条件下的功能梯度板弯曲问题。本文通过理论推导，将R-函数理论和变分原理结合，为求解挠度变化提供了一个新方法。这种方法不仅提高了求解过程的效率，而且也增加了模型的适用性和灵活性。本文通过结合具体的工程实例，展示了这种结合方法的有效性。通过对比分析，可以观察到本文方法的计算结果与有限元法计算结果的一致性，验证了本文方法的准确性和可靠性。
此外，本文还探讨了计算结果与变分法中使用的试函数选择、积分网格的划分等因素之间的关系。这种分析有助于优化计算过程，提高结果的精度，同时也为R-函数在功能梯度板弯曲问题中的应用提供了理论支持和实践指导。
总之，R-函数理论与变分法的结合为解决具有复杂边界条件的功能梯度板弯曲问题提供了一种新的视角和方法，这种结合能够提高计算效率。
DownLoad: