引　　言

轴对称结构作为一类常见的结构，一直是工程科学中的重点研究对象^[1-3]。轴对称结构常作为结构中的核心部件，广泛应用于通信、机械制造、航空航天以及土木工程等诸多领域^[4-6]。对于单一材料的均匀轴对称结构，它在各种荷载下的弹性行为已经有了很多研究^[7]。但在实际工程情况中，轴对称结构的不同位置需要满足强度、耐磨性、耐热性等不同需求，而单一材料难以满足要求。

为了使轴对称结构更好地满足实际使用需求，有学者引入了功能梯度材料。功能梯度材料（functionally graded materials，简称为FGMs）是指一种特殊的复合材料，其材料的微观组成和性能在空间上呈连续变化^[8]。部分学者对功能梯度材料的轴对称结构的力学问题进行了研究^[9-21]。Huang^[9]研究了功能梯度圆柱梁的弯曲和自由振动问题。Zenkour ^[10]研究了一个带有功能梯度过渡带的夹层旋转圆盘，在假设功能梯度过渡带的材料性能呈指数函数梯度形式变化的情况下，得到了固支和自由边界条件下转动圆盘的解析解。Dai等^[11]计算了功能梯度材料空心圆环在温度场变化下的位移场和应力场，但只求得了半解析解。Abdalla等^[12]假定该材料在径向上梯度变化，分析了功能梯度转动轴对称变厚度空心盘的热机械应力行为。Kalali等^[13]提出了一种功能梯度材料轴对称问题（转动盘、圆柱和球形容器）的弹塑性应力解。Bose和Rattan^[14]计算了线性变化功能梯度材料热梯度转动圆盘的稳态蠕变模型应力和应变率的分布。沈景凤等^[15]基于热弹性耦合理论建立了统一温度场的热耦合本构方程，求得了不同温度场下圆环碟片的位移控制方程。Madan等^[16]计算了基于S-FGM（材料参数呈S型变化形式）的旋转圆环在不同材料参数以及不同高宽比下的极限速度。张莹等^[17]利用England-Spencer板理论研究了材料梯度因子、板的厚度以及无量纲正应力对FGM圆环板的影响。Li等^[18-19]对功能梯度圆板在横向荷载作用下的弹性、磁电弹等问题做了系统的研究。Peng和Li^[20-21]使用Fredholm积分方程对等厚度各向同性功能梯度转动圆盘进行了弹性分析和热弹性分析。

上述研究工作中模型材料设置均为各向同性材料。由于功能梯度材料的制备工艺特点，在实际情况下，它们很少为各向同性材料^[22]。Yildirim^[23]使用互补函数法（CFM）来求得正交各向异性圆盘的数值解。Mahdavi等^[24]利用变材料性能理论研究了变厚度、等角速度功能梯度材料在热-机械荷载作用下的应力和应变。Yildirim等^[25]对变厚度极坐标正交各向异性转动圆环/圆盘的临界转速和稳定性进行了研究。目前有关各向异性功能梯度圆盘结构的相关力学问题研究已经取得了一定成果，但大多是针对材料性能呈固定梯度变化的情况，有关任意梯度变化的变厚度各向异性圆盘的研究还较少，而材料性能任意梯度变化时的通用求解方法无疑将为结构和材料的优化设计提供重要的理论指导意义。

鉴于此，我们考虑厚度和材料性能沿径向任意变化的转动圆盘，建立其绕刚性轴匀速转动的力学理论模型，给出一种有效的积分方程方法，将问题转化为关于径向应力的Fredholm积分方程，从而通过对Fredholm积分方程的数值求解来求出圆盘内应力场和位移场的分布情况，并研究了不同的材料梯度参数变化对应力场和位移场的影响。

1.   模 型 建 立

如图1所示，考虑一固结于刚性轴上的变厚度功能梯度圆盘，以匀角速度ω转动。圆盘内径为a，外径为b。圆盘由正交各向异性材料制成，其厚度、径向和环向弹性模量、密度沿径向任意变化，分别表示为$h\left( r \right)$，${E_r}\left( r \right)$，${E_\theta }\left( r \right)$，$\rho \left( r \right)$。

图  1  模型剖切示意图（左）与俯视示意图（右）

Figure  1.  Schematic diagram of the model: profile (left); top view (right)

下载: 全尺寸图片 幻灯片

假设圆盘厚度$h\left( r \right) < b/10$，则该问题为平面应力问题。建立如图1所示的极坐标系，则径向位移${u_r} \ne 0$，环向位移${u_\theta } = 0$。其物理方程和几何方程为

其中，${\sigma _r}$，${\sigma _\theta }$和${\varepsilon _r}$，${\varepsilon _\theta }$分别代表径向、环向应力和径向、环向应变，${\nu _{r\theta }}$，${\nu _{\theta r}}$为Poisson比，与径向和环向弹性模量之间的关系为

对于匀速转动的变厚度圆盘，其应力满足下面的平衡关系：

且边界条件为

2.   问 题 求 解

2.1   幂函数变化的变厚度功能梯度圆盘

首先考虑一种特殊情况，假设厚度、弹性模量和密度沿半径呈幂函数形式变化：

式中${h_0}$，${E_{r0}}$，${E_{\theta 0}}$，${\rho _0}$分别表示圆盘内侧对应的物理量；$\delta $，${\beta _1}$，${\beta _2}$分别为厚度、弹性模量、密度的梯度参数。因常见材料的Poisson比大都在0.3 ~ 0.5之间，对弹性场的影响很小，故假设为常数。

对于幂函数梯度分布模型，我们可得到圆盘内应力和位移场的精确解。将几何关系(2)和材料性能(6)代入物理方程(1)，继而代入平衡方程(4)，可得关于径向位移的微分方程：

上式为关于径向位移的Euler方程，其解为

式中C₁，C₂为待定常数，可由边界条件确定；C₃，n₁和n₂的具体形式如下：

因此，应力分量${\sigma _r}$和${\sigma _\theta }$继而可表示为

将式(10)代入边界条件(5)，可进一步由下列关系确定出待定常数C₁和C₂：

其中$ {\boldsymbol{C}},\;{\boldsymbol{Y}},\;{\boldsymbol{Z}} $分别为

2.2   任意梯度变化的变厚度功能梯度圆盘

值得说明的是，上述精确解(8)和(10)仅仅只适用于材料性能和厚度沿半径呈幂函数变化的形式，而实际工程和实验中很少存在这样的情况。为更符合实际情况，且基于对圆盘的结构和材料优化设计提供理论指导的角度，接下来我们考虑结构的厚度、材料常数和密度沿半径任意梯度分布的情况，并给出一种通用的求解方法——积分方程方法。

不同于前面的推导，首先将几何关系式(2)代入本构方程(1)，可得

上式可看成是关于径向位移${u_r}$的微分方程。因此，${u_r}$可由径向应力${\sigma _r}$表示为

其中，${D_1}$为可由边界条件确定的待定常数，

同时，环向应力${\sigma _\theta }$也用径向应力${\sigma _r}$和径向位移${u_r}$来表示。联立式(1)和式(2)，可得

只要确定了径向应力${\sigma _r}$，则径向位移${u_r}$和环向应力${\sigma _\theta }$即可由式(14)和式(15)确定。尤其值得说明的是，本文采用的方法不仅适用于材料性能连续分布的情况，对不连续或分段连续的情况也是可用的。同时由式(14)和式(15)可以看出，对于理想黏结的层合圆环或圆筒结构，不连续或分段连续的材料常数情况，结构的径向位移总是连续的，而环向应力则会因在每个界面上有明显的跳跃现象，这是导致层合结构出现脱层和开裂的主要原因。而具有连续变化材料性能的功能梯度圆盘可有效防止这一现象的出现。

将式(14)代入式(16)，再代入平衡方程(4)，有

上式为关于径向应力${\sigma _r}$的变系数微分-积分方程，其求解过程可参考文献[20]。式(17)两边同时对半径积分，可得关于${\sigma _r}$的积分方程：

其中函数${f_1}$，${f_2}$，${f_3}$的具体表达式为

考虑式(14)、(18)和边界条件(5)，可得${D_1}$和${D_2}$：

将得到的${D_1}$和${D_2}$回代入式(18)，整理后即可得到关于径向应力${\sigma _r}$的第二类Fredholm积分方程：

其中$ K\left( {r,s} \right) $是核函数,

有关第二类Fredholm积分方程的求解方法有很多种，接下来的求解中将主要采用文献[20]中的方法。一旦径向应力${\sigma _r}$求出，即可由式(14)和(16)确定出${u_r}$和${\sigma _\theta }$。

3.   数值算例和讨论

上节提出的积分方程方法适用于变厚度圆盘的材料常数沿半径任意梯度变化的情况。本节首先考虑厚度、弹性模量等参数为特殊的幂函数形式的情况，对此将由Fredholm积分方程求出的数值解与对应精确解进行对比，以及针对常见的Voigt模型，将由本方法算得的数值解和ANSYS有限元计算结果进行对比，来验证该方法的准确性和精度。然后，对于Voigt模型，分析厚度变化、材料性能梯度参数、各向异性度等对应力场和位移场的影响。在下面的算例中，均假设圆盘的外径为内径的两倍。

3.1   积分方程方法的验证

为验证方法的有效性以及精度，本小节考虑特殊的幂函数梯度变化情况。将式(6)代入径向应力${\sigma _r}$的第二类Fredholm积分方程，对由积分方程方法算得的数值解和解析解（式（8）、（10））进行对比（令$\alpha = \delta = {\beta _1} = {\beta _2}$），对比结果（归一化的径向应力、环向应力和径向位移）如图2所示。由图2可以看出，首先积分方程方法数值解已满足边界条件，同时当$\alpha $取不同值时，由积分方程方法算得的位移、应力数值解（ + 表示）与解析解（实线表示）吻合得非常好。值得说明的是当圆盘仅转动的时候，归一化的径向位移与幂函数的梯度参数无关（可由式（8）和式（9）确定）。

图  2  幂函数时数值解和解析解的比较（$\alpha = \delta = {\beta _1} = {\beta _2}$）

Figure  2.  A comparison of the exact and numerical results of the power law function （$\alpha = \delta = {\beta _1} = {\beta _2}$）

下载: 全尺寸图片 幻灯片

上述的验证通过在特殊情况与精确解的比较来进行。接下来，我们考虑一般的形式，假设圆盘的厚度和材料参数沿半径呈Voigt函数变化：

这种梯度变化形式下，无法获得精确解。我们将通过和有限元计算结果进行对比来继续验证本文方法的精度和有效性。

考虑转动圆盘的内表面（r=a）材料为Al₂O₃（E_r = 90.43 GPa，E_θ = 116.36 GPa，ρ = 3980 kg/m³，ν_θr = 0.28），外表面（r=b）材料为ZrO₂（E_r = 151 GPa，E_θ = 151 GPa，ρ = 5700 kg/m³，ν_θr = 0.3）。厚度参数h_a = 0.03，h_b = 0.01。同时假设式(23)中的各梯度参数相同，$\delta = {\beta _1} = {\beta _2} = \alpha $。采用ANSYS有限元软件，创建变厚度圆盘的板壳结构实体模型，采用四节点壳单元Shell 63对结构进行划分网格，环向网格密度为400份，径向网格密度为 100份。得到由ANSYS软件计算得到的结果和本文积分方程方法求得数值解的对比，如图3所示。对于不同厚度变化情况和梯度参数$\alpha $，本文获得的数值解和ANSYS有限元分析结果吻合得非常好，进一步说明了本文方法可行、有效。

图  3  Voigt函数时本文解和有限元解的比较（$\alpha = \delta = {\beta _1} = {\beta _2}$）

Figure  3.  A comparison of the present method results and the FEM solutions for the Voigt model （$\alpha = \delta = {\beta _1} = {\beta _2}$）

下载: 全尺寸图片 幻灯片

3.2   厚度变化对结构变形和应力分布的影响

本小节研究厚度变化对转动圆盘应力场和位移场的影响。本小节以及后面的分析中，均假设厚度参数h_a = 0.03，h_b = 0.01。考虑厚度沿径向呈Voigt函数变化（如式（23）所示）。显然当$\delta < 1$时，厚度变化沿径向为凹函数线型；当$\delta > 1$时，厚度变化沿径向为凸函数线型。无论δ如何变化，整个圆盘内部厚度的最大值最小值都不会超过内外径处的值。当$\delta $>0且有限时，h(a) = h_a，h(b) = h_b；当$\delta = 0$时，圆盘厚度一定，为h_b；当$\delta = \infty $时，圆盘厚度也为常数，但值变为h_a；当$\delta = 1$时，圆盘厚度则沿径向呈线性变化。圆盘材料采用均质ZrO₂。

由图4可知，无论δ如何变化，径向应力的最大值都存在于内径表面，这一点和均匀材料是相同的。但是在δ<1时，由于在靠近内径处的厚度变化较为显著，所以会有一个应力集中的现象。环向应力的最大值出现位置随着厚度参数$\delta $的增大，有向内径处移动的趋势，环向应力的最大值出现在结构内部，而不是在内外表面处，这是与等厚均匀材料存在区别的地方。当δ=1时，即当厚度沿径向呈线性变化时，应力场和位移场分布更加均匀。

图  4  厚度梯度参数δ对应力场和位移场的影响

Figure  4.  The effects of δ on the stress and radial displacement components

下载: 全尺寸图片 幻灯片

3.3   材料性能梯度参数对结构变形和应力分布的影响

本小节研究材料性能梯度参数对转动圆盘变形和应力的影响。材料参数设置参照3.2小节。与3.2小节中考虑厚度分布一样，将材料性能变化函数设置为Voigt函数（见式（23））。

${\,\beta _1}$是控制结构内部弹性模量函数的梯度指标。当${\,\beta _2}$=0，δ=1（厚度按式（23）呈线性分布）时，考虑不同${\,\beta _1}$，结构应力场和位移场的变化如图5所示。由于使用了Voigt函数来进行定义，当${\,\beta _1}$>0且有限时，整个圆盘结构内外表面的弹性模量值不变。当${\,\beta _1} = 0$时，整个结构的弹性模量${E_r} = {E_{rb}}$，${E_\theta } = {E_{\theta b}}$；当${\,\beta _1} = 1$时，圆盘结构密度变化呈线性变化；而当${\,\beta _1} = \infty $时，整个结构的弹性模量${E_r} = {E_{ra}}$，${E_\theta } = {E_{\theta a}}$。由图5可知，结构内部的环向应力和径向位移整体随着${\,\beta _1}$的增大而单调增大，而径向应力几乎没有改变。所以在进行结构设计时，可以尽量减小结构的弹性模量值以减小模型内部的环向应力和位移值。

图  5  弹性模量梯度参数β₁对应力场和位移场的影响

Figure  5.  The effects of β₁ on the stress and radial displacement components

下载: 全尺寸图片 幻灯片

${\,\beta _2}$是控制结构内部密度分布的梯度指标。当${\,\beta _1}$=0，δ=1时，图6中考虑不同${\,\beta _2}$对结构应力场和位移场的影响。同样的，当${\,\beta _2} = 0$时，整个结构的密度为${\rho _b}$；当${\,\beta _2} = 1$时，圆盘结构密度变化沿径向呈线性变化；而当${\,\beta _2} = \infty $时，整个结构的密度为${\rho _a}$。由图6可知，随着密度参数${\,\beta _2}$的减小，整个结构的径向应力、环向应力以及径向位移均单调增大，这种现象主要是由于${\,\beta _2}$减小后造成结构自重的增加，使得结构旋转时的惯性力增加造成的，所以我们在满足结构承载要求的情况下要尽量增大${\,\beta _2}$的值。

图  6  密度梯度参数β₂对应力场和位移场的影响

Figure  6.  The effects of β₂ on the stress and radial displacement components

下载: 全尺寸图片 幻灯片

3.4   各向异性度对结构变形和应力分布的影响

本小节讨论材料各向异性度对转动圆盘变形和应力的影响。各向异性度为两个方向弹性模量之比：

算例通过固定${E_r}$，调整${E_\theta }$的值来改变各向异性度的大小。圆盘的结构和其他材料参数均与3.3小节中给定的值一致。由图7可知，当固定${E_r}$的值，增大各向异性度时，结构内部的径向应力分布均单调递减，而环向应力在靠近内壁处单调递减，在远离内壁处单调递增，并且整个结构内环向应力的最大值有向外径移动的趋势。在环向应力分布中，存在着一个截面，在该截面处，无论各向异性度如何变化，环向应力值不变。径向位移方面，在$\lambda $取值到1附近时会出现径向位移的最大值。由此可知我们在进行模型设计时，不能盲目增大或者减小各向异性度，要考虑侧重点是在径向应力方面还是环向应力等其他方面，再根据优化目标进行调整。

图  7  各向异性度参数λ对应力场和位移场的影响

Figure  7.  The effects of λ on the stress and radial displacement components

下载: 全尺寸图片 幻灯片

4.   结　　论

本文采用积分方程方法，推导了变厚度各向异性功能梯度转动圆盘的轴对称平面应力问题中关于径向应力的Fredholm积分方程，并采用数值方法对该积分方程进行了求解。此方法不限定结构几何参数和材料性能参数的具体变化形式。通过对比幂函数材料参数梯度分布条件下的精确解、Voigt模型时的ANSYS有限元计算结果验证了本方法的合理性。在此基础上，分别分析讨论了呈Voigt函数变化的厚度分布和弹性模量密度梯度分布、各向异性度等对结构应力场和位移场的影响。主要结论有：

1） 对于厚度呈Voigt函数变化的情况，当Voigt函数为线性函数时，结构应力场和位移场整体最小。

2） 结构的变形和应力分布随材料的材料性能（弹性模量、密度）单调变化。具体的依赖关系，需要考虑结构内外壁的材料参数。

3） 各向异性度（环向和径向模量的比值）的增大会降低圆盘最大径向应力，但是会显著增加圆盘的最大环向应力。

4） 本文提出的求解方法适用于结构几何参数和材料性能沿径向任意梯度变化的情况。

致谢　本文作者衷心感谢长沙理工大学研究生科研创新项目（CX2021SS133）对本文的资助。