變厚度各向異性功能梯度轉(zhuǎn)動(dòng)圓盤的彈性分析*

彭旭龍，黃海平，李進(jìn)寶，陳 央，陳子光

（1.長沙理工大學(xué) 土木工程學(xué)院 力學(xué)系，長沙 410114；2.華中科技大學(xué) 航空航天學(xué)院 工程力學(xué)系，武漢 430074）

引言

軸對稱結(jié)構(gòu)作為一類常見的結(jié)構(gòu)，一直是工程科學(xué)中的重點(diǎn)研究對象[1-3].軸對稱結(jié)構(gòu)常作為結(jié)構(gòu)中的核心部件，廣泛應(yīng)用于通信、機(jī)械制造、航空航天以及土木工程等諸多領(lǐng)域[4-6].對于單一材料的均勻軸對稱結(jié)構(gòu)，它在各種荷載下的彈性行為已經(jīng)有了很多研究[7].但在實(shí)際工程情況中，軸對稱結(jié)構(gòu)的不同位置需要滿足強(qiáng)度、耐磨性、耐熱性等不同需求，而單一材料難以滿足要求.

為了使軸對稱結(jié)構(gòu)更好地滿足實(shí)際使用需求，有學(xué)者引入了功能梯度材料.功能梯度材料（functionally graded materials，簡稱為FGMs）是指一種特殊的復(fù)合材料，其材料的微觀組成和性能在空間上呈連續(xù)變化[8].部分學(xué)者對功能梯度材料的軸對稱結(jié)構(gòu)的力學(xué)問題進(jìn)行了研究[9-21].Huang[9]研究了功能梯度圓柱梁的彎曲和自由振動(dòng)問題.Zenkour[10]研究了一個(gè)帶有功能梯度過渡帶的夾層旋轉(zhuǎn)圓盤，在假設(shè)功能梯度過渡帶的材料性能呈指數(shù)函數(shù)梯度形式變化的情況下，得到了固支和自由邊界條件下轉(zhuǎn)動(dòng)圓盤的解析解.Dai 等[11]計(jì)算了功能梯度材料空心圓環(huán)在溫度場變化下的位移場和應(yīng)力場，但只求得了半解析解.Abdalla 等[12]假定該材料在徑向上梯度變化，分析了功能梯度轉(zhuǎn)動(dòng)軸對稱變厚度空心盤的熱機(jī)械應(yīng)力行為.Kalali 等[13]提出了一種功能梯度材料軸對稱問題（轉(zhuǎn)動(dòng)盤、圓柱和球形容器）的彈塑性應(yīng)力解.Bose 和Rattan[14]計(jì)算了線性變化功能梯度材料熱梯度轉(zhuǎn)動(dòng)圓盤的穩(wěn)態(tài)蠕變模型應(yīng)力和應(yīng)變率的分布.沈景鳳等[15]基于熱彈性耦合理論建立了統(tǒng)一溫度場的熱耦合本構(gòu)方程，求得了不同溫度場下圓環(huán)碟片的位移控制方程.Madan 等[16]計(jì)算了基于S-FGM（材料參數(shù)呈S 型變化形式）的旋轉(zhuǎn)圓環(huán)在不同材料參數(shù)以及不同高寬比下的極限速度.張瑩等[17]利用England-Spencer 板理論研究了材料梯度因子、板的厚度以及無量綱正應(yīng)力對FGM 圓環(huán)板的影響.Li 等[18-19]對功能梯度圓板在橫向荷載作用下的彈性、磁電彈等問題做了系統(tǒng)的研究.Peng 和Li[20-21]使用Fredholm 積分方程對等厚度各向同性功能梯度轉(zhuǎn)動(dòng)圓盤進(jìn)行了彈性分析和熱彈性分析.

上述研究工作中模型材料設(shè)置均為各向同性材料.由于功能梯度材料的制備工藝特點(diǎn)，在實(shí)際情況下，它們很少為各向同性材料[22].Yildirim[23]使用互補(bǔ)函數(shù)法（CFM）來求得正交各向異性圓盤的數(shù)值解.Mahdavi等[24]利用變材料性能理論研究了變厚度、等角速度功能梯度材料在熱-機(jī)械荷載作用下的應(yīng)力和應(yīng)變.Yildirim 等[25]對變厚度極坐標(biāo)正交各向異性轉(zhuǎn)動(dòng)圓環(huán)/圓盤的臨界轉(zhuǎn)速和穩(wěn)定性進(jìn)行了研究.目前有關(guān)各向異性功能梯度圓盤結(jié)構(gòu)的相關(guān)力學(xué)問題研究已經(jīng)取得了一定成果，但大多是針對材料性能呈固定梯度變化的情況，有關(guān)任意梯度變化的變厚度各向異性圓盤的研究還較少，而材料性能任意梯度變化時(shí)的通用求解方法無疑將為結(jié)構(gòu)和材料的優(yōu)化設(shè)計(jì)提供重要的理論指導(dǎo)意義.

鑒于此，我們考慮厚度和材料性能沿徑向任意變化的轉(zhuǎn)動(dòng)圓盤，建立其繞剛性軸勻速轉(zhuǎn)動(dòng)的力學(xué)理論模型，給出一種有效的積分方程方法，將問題轉(zhuǎn)化為關(guān)于徑向應(yīng)力的Fredholm 積分方程，從而通過對Fredholm 積分方程的數(shù)值求解來求出圓盤內(nèi)應(yīng)力場和位移場的分布情況，并研究了不同的材料梯度參數(shù)變化對應(yīng)力場和位移場的影響.

1 模型建立

h(r)

如圖1 所示，考慮一固結(jié)于剛性軸上的變厚度功能梯度圓盤，以勻角速度ω 轉(zhuǎn)動(dòng).圓盤內(nèi)徑為a，外徑為b.圓盤由正交各向異性材料制成，其厚度、徑向和環(huán)向彈性模量、密度沿徑向任意變化，分別表示為，Er(r)，Eθ(r)，ρ(r).

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

假設(shè)圓盤厚度h(r)＜b/10，則該問題為平面應(yīng)力問題.建立如圖1 所示的極坐標(biāo)系，則徑向位移ur≠0，環(huán)向位移uθ=0.其物理方程和幾何方程為

其中，σr，σθ和 εr，εθ分別代表徑向、環(huán)向應(yīng)力和徑向、環(huán)向應(yīng)變，νrθ，νθr為Poisson 比，與徑向和環(huán)向彈性模量之間的關(guān)系為

對于勻速轉(zhuǎn)動(dòng)的變厚度圓盤，其應(yīng)力滿足下面的平衡關(guān)系：

且邊界條件為

2 問題求解

2.1 冪函數(shù)變化的變厚度功能梯度圓盤

首先考慮一種特殊情況，假設(shè)厚度、彈性模量和密度沿半徑呈冪函數(shù)形式變化：

式中h0，Er0，Eθ0，ρ0分別表示圓盤內(nèi)側(cè)對應(yīng)的物理量；δ，β1，β2分別為厚度、彈性模量、密度的梯度參數(shù).因常見材料的Poisson 比大都在0.3 ～ 0.5 之間，對彈性場的影響很小，故假設(shè)為常數(shù).

對于冪函數(shù)梯度分布模型，我們可得到圓盤內(nèi)應(yīng)力和位移場的精確解.將幾何關(guān)系(2)和材料性能(6)代入物理方程(1)，繼而代入平衡方程(4)，可得關(guān)于徑向位移的微分方程：

上式為關(guān)于徑向位移的Euler 方程，其解為

式中C1，C2為待定常數(shù)，可由邊界條件確定；C3，n1和n2的具體形式如下：

因此，應(yīng)力分量 σr和 σθ繼而可表示為

將式(10)代入邊界條件(5)，可進(jìn)一步由下列關(guān)系確定出待定常數(shù)C1和C2：

其中C,Y,Z分別為

2.2 任意梯度變化的變厚度功能梯度圓盤

值得說明的是，上述精確解(8)和(10)僅僅只適用于材料性能和厚度沿半徑呈冪函數(shù)變化的形式，而實(shí)際工程和實(shí)驗(yàn)中很少存在這樣的情況.為更符合實(shí)際情況，且基于對圓盤的結(jié)構(gòu)和材料優(yōu)化設(shè)計(jì)提供理論指導(dǎo)的角度，接下來我們考慮結(jié)構(gòu)的厚度、材料常數(shù)和密度沿半徑任意梯度分布的情況，并給出一種通用的求解方法——積分方程方法.

不同于前面的推導(dǎo)，首先將幾何關(guān)系式(2)代入本構(gòu)方程(1)，可得

上式可看成是關(guān)于徑向位移ur的微分方程.因此，ur可由徑向應(yīng)力 σr表示為

其中，D1為可由邊界條件確定的待定常數(shù)，

同時(shí)，環(huán)向應(yīng)力 σθ也用徑向應(yīng)力 σr和徑向位移ur來表示.聯(lián)立式(1)和式(2)，可得

只要確定了徑向應(yīng)力 σr，則徑向位移ur和環(huán)向應(yīng)力 σθ即可由式(14)和式(15)確定.尤其值得說明的是，本文采用的方法不僅適用于材料性能連續(xù)分布的情況，對不連續(xù)或分段連續(xù)的情況也是可用的.同時(shí)由式(14)和式(15)可以看出，對于理想黏結(jié)的層合圓環(huán)或圓筒結(jié)構(gòu)，不連續(xù)或分段連續(xù)的材料常數(shù)情況，結(jié)構(gòu)的徑向位移總是連續(xù)的，而環(huán)向應(yīng)力則會(huì)因在每個(gè)界面上有明顯的跳躍現(xiàn)象，這是導(dǎo)致層合結(jié)構(gòu)出現(xiàn)脫層和開裂的主要原因.而具有連續(xù)變化材料性能的功能梯度圓盤可有效防止這一現(xiàn)象的出現(xiàn).

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

上式為關(guān)于徑向應(yīng)力 σr的變系數(shù)微分-積分方程，其求解過程可參考文獻(xiàn)[20].式(17)兩邊同時(shí)對半徑積分，可得關(guān)于 σr的積分方程：

其中函數(shù)f1，f2，f3的具體表達(dá)式為

考慮式(14)、(18)和邊界條件(5)，可得D1和D2：

將得到的D1和D2回代入式(18)，整理后即可得到關(guān)于徑向應(yīng)力 σr的第二類Fredholm 積分方程：

其中K(r,s)是核函數(shù),

有關(guān)第二類Fre dholm 積分方程的求解方法有很多種，接下來的求解中將主要采用文獻(xiàn)[20]中的方法.一旦徑向應(yīng)力 σr求出，即可由式(14)和(16)確定出ur和 σθ.

3 數(shù)值算例和討論

上節(jié)提出的積分方程方法適用于變厚度圓盤的材料常數(shù)沿半徑任意梯度變化的情況.本節(jié)首先考慮厚度、彈性模量等參數(shù)為特殊的冪函數(shù)形式的情況，對此將由Fredholm 積分方程求出的數(shù)值解與對應(yīng)精確解進(jìn)行對比，以及針對常見的Voigt 模型，將由本方法算得的數(shù)值解和ANSYS 有限元計(jì)算結(jié)果進(jìn)行對比，來驗(yàn)證該方法的準(zhǔn)確性和精度.然后，對于Voigt 模型，分析厚度變化、材料性能梯度參數(shù)、各向異性度等對應(yīng)力場和位移場的影響.在下面的算例中，均假設(shè)圓盤的外徑為內(nèi)徑的兩倍.

3.1 積分方程方法的驗(yàn)證

為驗(yàn)證方法的有效性以及精度，本小節(jié)考慮特殊的冪函數(shù)梯度變化情況.將式(6)代入徑向應(yīng)力 σr的第二類Fredholm 積分方程，對由積分方程方法算得的數(shù)值解和解析解（式（8）、（10））進(jìn)行對比（令α=δ=β1=β2），對比結(jié)果（歸一化的徑向應(yīng)力、環(huán)向應(yīng)力和徑向位移）如圖2 所示.由圖2 可以看出，首先積分方程方法數(shù)值解已滿足邊界條件，同時(shí)當(dāng) α取不同值時(shí)，由積分方程方法算得的位移、應(yīng)力數(shù)值解（+ 表示）與解析解（實(shí)線表示）吻合得非常好.值得說明的是當(dāng)圓盤僅轉(zhuǎn)動(dòng)的時(shí)候，歸一化的徑向位移與冪函數(shù)的梯度參數(shù)無關(guān)（可由式（8）和式（9）確定）.

圖2 冪函數(shù)時(shí)數(shù)值解和解析解的比較（α=δ=β1=β2）Fig.2 A comparison of the exact and numerical results of the power law function （α=δ=β1=β2）

上述的驗(yàn)證通過在特殊情況與精確解的比較來進(jìn)行.接下來，我們考慮一般的形式，假設(shè)圓盤的厚度和材料參數(shù)沿半徑呈Voigt 函數(shù)變化：

這種梯度變化形式下，無法獲得精確解.我們將通過和有限元計(jì)算結(jié)果進(jìn)行對比來繼續(xù)驗(yàn)證本文方法的精度和有效性.

考慮轉(zhuǎn)動(dòng)圓盤的內(nèi)表面（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）.厚度參數(shù)ha= 0.03，hb=0.01.同時(shí)假設(shè)式(23)中的各梯度參數(shù)相同，δ=β1=β2=α.采用ANSYS 有限元軟件，創(chuàng)建變厚度圓盤的板殼結(jié)構(gòu)實(shí)體模型，采用四節(jié)點(diǎn)殼單元Shell 63 對結(jié)構(gòu)進(jìn)行劃分網(wǎng)格，環(huán)向網(wǎng)格密度為400 份，徑向網(wǎng)格密度為100 份.得到由ANSYS 軟件計(jì)算得到的結(jié)果和本文積分方程方法求得數(shù)值解的對比，如圖3 所示.對于不同厚度變化情況和梯度參數(shù)α，本文獲得的數(shù)值解和ANSYS 有限元分析結(jié)果吻合得非常好，進(jìn)一步說明了本文方法可行、有效.

圖3 Voigt 函數(shù)時(shí)本文解和有限元解的比較（α=δ=β1=β2）Fig.3 A comparison of the present method results and the FEM solutions for the Voigt model （α=δ=β1=β2）

3.2 厚度變化對結(jié)構(gòu)變形和應(yīng)力分布的影響

本小節(jié)研究厚度變化對轉(zhuǎn)動(dòng)圓盤應(yīng)力場和位移場的影響.本小節(jié)以及后面的分析中，均假設(shè)厚度參數(shù)ha=0.03，hb= 0.01.考慮厚度沿徑向呈Voigt 函數(shù)變化（如式（23）所示）.顯然當(dāng)δ ＜1時(shí)，厚度變化沿徑向?yàn)榘己瘮?shù)線型；當(dāng)δ ＞1時(shí)，厚度變化沿徑向?yàn)橥购瘮?shù)線型.無論δ 如何變化，整個(gè)圓盤內(nèi)部厚度的最大值最小值都不會(huì)超過內(nèi)外徑處的值.當(dāng)δ＞0 且有限時(shí)，h(a)=ha，h(b)=hb；當(dāng)δ=0時(shí)，圓盤厚度一定，為hb；當(dāng)δ=∞時(shí)，圓盤厚度也為常數(shù)，但值變?yōu)閔a；當(dāng)δ=1時(shí)，圓盤厚度則沿徑向呈線性變化.圓盤材料采用均質(zhì)ZrO2.

由圖4 可知，無論δ 如何變化，徑向應(yīng)力的最大值都存在于內(nèi)徑表面，這一點(diǎn)和均勻材料是相同的.但是在δ＜1 時(shí)，由于在靠近內(nèi)徑處的厚度變化較為顯著，所以會(huì)有一個(gè)應(yīng)力集中的現(xiàn)象.環(huán)向應(yīng)力的最大值出現(xiàn)位置隨著厚度參數(shù)δ的增大，有向內(nèi)徑處移動(dòng)的趨勢，環(huán)向應(yīng)力的最大值出現(xiàn)在結(jié)構(gòu)內(nèi)部，而不是在內(nèi)外表面處，這是與等厚均勻材料存在區(qū)別的地方.當(dāng)δ=1 時(shí)，即當(dāng)厚度沿徑向呈線性變化時(shí)，應(yīng)力場和位移場分布更加均勻.

圖4 厚度梯度參數(shù)δ 對應(yīng)力場和位移場的影響Fig.4 The effects of δ on the stress and radial displacement components

3.3 材料性能梯度參數(shù)對結(jié)構(gòu)變形和應(yīng)力分布的影響

本小節(jié)研究材料性能梯度參數(shù)對轉(zhuǎn)動(dòng)圓盤變形和應(yīng)力的影響.材料參數(shù)設(shè)置參照3.2 小節(jié).與3.2 小節(jié)中考慮厚度分布一樣，將材料性能變化函數(shù)設(shè)置為Voigt 函數(shù)（見式（23））.

β1是控制結(jié)構(gòu)內(nèi)部彈性模量函數(shù)的梯度指標(biāo).當(dāng) β2=0，δ=1（厚度按式（23）呈線性分布）時(shí)，考慮不同 β1，結(jié)構(gòu)應(yīng)力場和位移場的變化如圖5 所示.由于使用了Voigt 函數(shù)來進(jìn)行定義，當(dāng) β1＞0 且有限時(shí)，整個(gè)圓盤結(jié)構(gòu)內(nèi)外表面的彈性模量值不變.當(dāng)β1=0時(shí)，整個(gè)結(jié)構(gòu)的彈性模量Er=Erb，Eθ=Eθb；當(dāng)β1=1時(shí)，圓盤結(jié)構(gòu)密度變化呈線性變化；而當(dāng)β1=∞時(shí)，整個(gè)結(jié)構(gòu)的彈性模量Er=Era，Eθ=Eθa.由圖5 可知，結(jié)構(gòu)內(nèi)部的環(huán)向應(yīng)力和徑向位移整體隨著 β1的增大而單調(diào)增大，而徑向應(yīng)力幾乎沒有改變.所以在進(jìn)行結(jié)構(gòu)設(shè)計(jì)時(shí)，可以盡量減小結(jié)構(gòu)的彈性模量值以減小模型內(nèi)部的環(huán)向應(yīng)力和位移值.

圖5 彈性模量梯度參數(shù)β1 對應(yīng)力場和位移場的影響Fig.5 The effects of β1 on the stress and radial displacement components

β2是控制結(jié)構(gòu)內(nèi)部密度分布的梯度指標(biāo).當(dāng) β1=0，δ=1 時(shí)，圖6 中考慮不同 β2對結(jié)構(gòu)應(yīng)力場和位移場的影響.同樣的，當(dāng)β2=0時(shí)，整個(gè)結(jié)構(gòu)的密度為 ρb；當(dāng)β2=1時(shí)，圓盤結(jié)構(gòu)密度變化沿徑向呈線性變化；而當(dāng)β2=∞時(shí)，整個(gè)結(jié)構(gòu)的密度為 ρa(bǔ).由圖6 可知，隨著密度參數(shù) β2的減小，整個(gè)結(jié)構(gòu)的徑向應(yīng)力、環(huán)向應(yīng)力以及徑向位移均單調(diào)增大，這種現(xiàn)象主要是由于 β2減小后造成結(jié)構(gòu)自重的增加，使得結(jié)構(gòu)旋轉(zhuǎn)時(shí)的慣性力增加造成的，所以我們在滿足結(jié)構(gòu)承載要求的情況下要盡量增大 β2的值.

圖6 密度梯度參數(shù)β2 對應(yīng)力場和位移場的影響Fig.6 The effects of β2 on the stress and radial displacement components

3.4 各向異性度對結(jié)構(gòu)變形和應(yīng)力分布的影響

本小節(jié)討論材料各向異性度對轉(zhuǎn)動(dòng)圓盤變形和應(yīng)力的影響.各向異性度為兩個(gè)方向彈性模量之比：

算例通過固定Er，調(diào)整Eθ的值來改變各向異性度的大小.圓盤的結(jié)構(gòu)和其他材料參數(shù)均與3.3 小節(jié)中給定的值一致.由圖7 可知，當(dāng)固定Er的值，增大各向異性度時(shí)，結(jié)構(gòu)內(nèi)部的徑向應(yīng)力分布均單調(diào)遞減，而環(huán)向應(yīng)力在靠近內(nèi)壁處單調(diào)遞減，在遠(yuǎn)離內(nèi)壁處單調(diào)遞增，并且整個(gè)結(jié)構(gòu)內(nèi)環(huán)向應(yīng)力的最大值有向外徑移動(dòng)的趨勢.在環(huán)向應(yīng)力分布中，存在著一個(gè)截面，在該截面處，無論各向異性度如何變化，環(huán)向應(yīng)力值不變.徑向位移方面，在λ取值到1 附近時(shí)會(huì)出現(xiàn)徑向位移的最大值.由此可知我們在進(jìn)行模型設(shè)計(jì)時(shí)，不能盲目增大或者減小各向異性度，要考慮側(cè)重點(diǎn)是在徑向應(yīng)力方面還是環(huán)向應(yīng)力等其他方面，再根據(jù)優(yōu)化目標(biāo)進(jìn)行調(diào)整.

圖7 各向異性度參數(shù)λ 對應(yīng)力場和位移場的影響Fig.7 The effects of λ on the stress and radial displacement components

4 結(jié)論

本文采用積分方程方法，推導(dǎo)了變厚度各向異性功能梯度轉(zhuǎn)動(dòng)圓盤的軸對稱平面應(yīng)力問題中關(guān)于徑向應(yīng)力的Fredholm 積分方程，并采用數(shù)值方法對該積分方程進(jìn)行了求解.此方法不限定結(jié)構(gòu)幾何參數(shù)和材料性能參數(shù)的具體變化形式.通過對比冪函數(shù)材料參數(shù)梯度分布條件下的精確解、Voigt 模型時(shí)的ANSYS 有限元計(jì)算結(jié)果驗(yàn)證了本方法的合理性.在此基礎(chǔ)上，分別分析討論了呈Voigt 函數(shù)變化的厚度分布和彈性模量密度梯度分布、各向異性度等對結(jié)構(gòu)應(yīng)力場和位移場的影響.主要結(jié)論有：

1）對于厚度呈Voigt 函數(shù)變化的情況，當(dāng)Voigt 函數(shù)為線性函數(shù)時(shí)，結(jié)構(gòu)應(yīng)力場和位移場整體最小.

2）結(jié)構(gòu)的變形和應(yīng)力分布隨材料的材料性能（彈性模量、密度）單調(diào)變化.具體的依賴關(guān)系，需要考慮結(jié)構(gòu)內(nèi)外壁的材料參數(shù).

3）各向異性度（環(huán)向和徑向模量的比值）的增大會(huì)降低圓盤最大徑向應(yīng)力，但是會(huì)顯著增加圓盤的最大環(huán)向應(yīng)力.

4）本文提出的求解方法適用于結(jié)構(gòu)幾何參數(shù)和材料性能沿徑向任意梯度變化的情況.

致謝本文作者衷心感謝長沙理工大學(xué)研究生科研創(chuàng)新項(xiàng)目（CX2021SS133）對本文的資助.