变厚度各向异性功能梯度转动圆盘的弹性分析*

彭旭龙，黄海平，李进宝，陈 央，陈子光

（1.长沙理工大学 土木工程学院 力学系，长沙 410114；2.华中科技大学 航空航天学院 工程力学系，武汉 430074）

引言

轴对称结构作为一类常见的结构，一直是工程科学中的重点研究对象[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 模型建立

h(r)

如图1 所示，考虑一固结于刚性轴上的变厚度功能梯度圆盘，以匀角速度ω 转动.圆盘内径为a，外径为b.圆盘由正交各向异性材料制成，其厚度、径向和环向弹性模量、密度沿径向任意变化，分别表示为，Er(r)，Eθ(r)，ρ(r).

图1 模型剖切示意图（左）与俯视示意图（右）Fig.1 Schematic diagram of the model:profile(left);top view(right)

假设圆盘厚度h(r)＜b/10，则该问题为平面应力问题.建立如图1 所示的极坐标系，则径向位移ur≠0，环向位移uθ=0.其物理方程和几何方程为

其中，σr，σθ和 εr，εθ分别代表径向、环向应力和径向、环向应变，νrθ，νθr为Poisson 比，与径向和环向弹性模量之间的关系为

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

且边界条件为

2 问题求解

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

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

式中h0，Er0，Eθ0，ρ0分别表示圆盘内侧对应的物理量；δ，β1，β2分别为厚度、弹性模量、密度的梯度参数.因常见材料的Poisson 比大都在0.3 ～ 0.5 之间，对弹性场的影响很小，故假设为常数.

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

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

式中C1，C2为待定常数，可由边界条件确定；C3，n1和n2的具体形式如下：

因此，应力分量 σr和 σθ继而可表示为

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

其中C,Y,Z分别为

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

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

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

上式可看成是关于径向位移ur的微分方程.因此，ur可由径向应力 σr表示为

其中，D1为可由边界条件确定的待定常数，

同时，环向应力 σθ也用径向应力 σr和径向位移ur来表示.联立式(1)和式(2)，可得

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

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

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

其中函数f1，f2，f3的具体表达式为

考虑式(14)、(18)和边界条件(5)，可得D1和D2：

将得到的D1和D2回代入式(18)，整理后即可得到关于径向应力 σr的第二类Fredholm 积分方程：

其中K(r,s)是核函数,

有关第二类Fre dholm 积分方程的求解方法有很多种，接下来的求解中将主要采用文献[20]中的方法.一旦径向应力 σr求出，即可由式(14)和(16)确定出ur和 σθ.

3 数值算例和讨论

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

3.1 积分方程方法的验证

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

图2 幂函数时数值解和解析解的比较（α=δ=β1=β2）Fig.2 A comparison of the exact and numerical results of the power law function （α=δ=β1=β2）

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

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

考虑转动圆盘的内表面（r=a）材料为Al2O3（Er= 90.43 GPa，Eθ= 116.36 GPa，ρ = 3 980 kg/m3，νθr= 0.28），外表面（r=b）材料为ZrO2（Er= 151 GPa，Eθ= 151 GPa，ρ = 5 700 kg/m3，νθr= 0.3）.厚度参数ha= 0.03，hb=0.01.同时假设式(23)中的各梯度参数相同，δ=β1=β2=α.采用ANSYS 有限元软件，创建变厚度圆盘的板壳结构实体模型，采用四节点壳单元Shell 63 对结构进行划分网格，环向网格密度为400 份，径向网格密度为100 份.得到由ANSYS 软件计算得到的结果和本文积分方程方法求得数值解的对比，如图3 所示.对于不同厚度变化情况和梯度参数α，本文获得的数值解和ANSYS 有限元分析结果吻合得非常好，进一步说明了本文方法可行、有效.

图3 Voigt 函数时本文解和有限元解的比较（α=δ=β1=β2）Fig.3 A comparison of the present method results and the FEM solutions for the Voigt model （α=δ=β1=β2）

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

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

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

图4 厚度梯度参数δ 对应力场和位移场的影响Fig.4 The effects of δ on the stress and radial displacement components

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

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

β1是控制结构内部弹性模量函数的梯度指标.当 β2=0，δ=1（厚度按式（23）呈线性分布）时，考虑不同 β1，结构应力场和位移场的变化如图5 所示.由于使用了Voigt 函数来进行定义，当 β1＞0 且有限时，整个圆盘结构内外表面的弹性模量值不变.当β1=0时，整个结构的弹性模量Er=Erb，Eθ=Eθb；当β1=1时，圆盘结构密度变化呈线性变化；而当β1=∞时，整个结构的弹性模量Er=Era，Eθ=Eθa.由图5 可知，结构内部的环向应力和径向位移整体随着 β1的增大而单调增大，而径向应力几乎没有改变.所以在进行结构设计时，可以尽量减小结构的弹性模量值以减小模型内部的环向应力和位移值.

图5 弹性模量梯度参数β1 对应力场和位移场的影响Fig.5 The effects of β1 on the stress and radial displacement components

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

图6 密度梯度参数β2 对应力场和位移场的影响Fig.6 The effects of β2 on the stress and radial displacement components

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

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

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

图7 各向异性度参数λ 对应力场和位移场的影响Fig.7 The effects of λ on the stress and radial displacement components

4 结论

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

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

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

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

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

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