朱竑禎，王緯波，高存法，殷學(xué)文

（1.南京航空航天大學(xué) 機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點(diǎn)實(shí)驗室，南京210016；2.中國船舶科學(xué)研究中心 船舶振動噪聲重點(diǎn)實(shí)驗室，江蘇 無錫214082）

0 引 言

船舶的振動噪聲影響船舶結(jié)構(gòu)的安全性、設(shè)備的靈敏性、船員的舒適性以及艦艇的隱聲性能等。開展船舶板類結(jié)構(gòu)的振動與噪聲控制是一項十分重要的工作。隨著艦艇低頻噪聲問題的日益突出，智能結(jié)構(gòu)在船舶結(jié)構(gòu)物中的應(yīng)用與智能結(jié)構(gòu)的主動控制研究，受到了各國的重視。美國最新型的海浪級核潛艇中裝備了振動噪聲主動控制系統(tǒng)與其全艦噪聲監(jiān)控系統(tǒng)配合使用，取得了良好的聲隱身性能。波蘭JERZY WICIAK對承受流體載荷力的水下結(jié)構(gòu)的聲輻射進(jìn)行研究，對粘貼有壓電作動器和傳感器的圓形板的主動噪聲控制進(jìn)行了試驗驗證。美國麻省理工學(xué)院正在研制一種增強(qiáng)型主動降噪潛艇殼體模型，殼體周圍粘貼有壓電陶瓷。

功能梯度材料（Functionally graded material,FGM）作為一種材料特性隨空間位置連續(xù)梯度變化的新型多相材料，可使水下航行器殼體兩側(cè)呈現(xiàn)截然不同的性能與功能，使結(jié)構(gòu)能夠更好地適應(yīng)超高溫、大溫度梯度等特殊極限環(huán)境。雖然國內(nèi)外針對采用壓電傳感器、驅(qū)動器的智能結(jié)構(gòu)主動控制進(jìn)行了廣泛研究，但是將壓電材料用于FGM結(jié)構(gòu)的主動控制研究相對較少，針對水下功能梯度材料板殼的主動控制文獻(xiàn)則更少。代峰等[1]建立了一種用于溫度梯度作用下功能梯度材料板殼力學(xué)性能分析的固體殼單元，針對壓電傳感器、驅(qū)動器的功能梯度材料板殼的熱變形及振動的主動控制進(jìn)行了數(shù)值仿真，探討了組份材料體積分?jǐn)?shù)冪指數(shù)分布、控制增益對功能梯度材料力學(xué)性能和控制效果的影響。Liew等[2-4]首先推導(dǎo)了含有壓電傳感器與作動器的FGM板的振動控制方程，針對溫度變化和隨厚度變化的功能梯度材料，采用位移—速度反饋算法從動態(tài)與靜態(tài)兩方面進(jìn)行壓電功能梯度材料的主動控制，并進(jìn)行了試驗驗證。為了深入研究FGM材料應(yīng)用于水下結(jié)構(gòu)的主動控制，其力學(xué)行為是目前重要研究方向。

功能梯度梁、板和殼結(jié)構(gòu)均有相應(yīng)的解析解、半解析解和數(shù)值解的研究[5]。例如，于天崇等[6]研究了功能梯度材料面內(nèi)變剛度矩形板的彎曲問題，該研究假定彎曲剛度沿板寬方向呈冪指數(shù)變化，設(shè)撓度為Levy解的形式求解基本方程。Zhang等[7]研究了含圓孔的功能梯度材料板的熱應(yīng)力問題，以圓孔中心為原點(diǎn)，假設(shè)楊氏模量、泊松比以及熱膨脹系數(shù)關(guān)于原點(diǎn)的距離呈級數(shù)展開形式，以彈性力學(xué)中極坐標(biāo)下的平衡方程、物理方程和本構(gòu)關(guān)系為基本方程，求解出應(yīng)力值。Yang和Gao[8]利用復(fù)變函數(shù)法和分層均勻化的思想，將平板在面內(nèi)分為若干個均勻?qū)?，相鄰層間位移與應(yīng)力連續(xù)，求得了無限大含圓孔功能梯度板的應(yīng)力解。楊正光[9]由壓電動力學(xué)平衡方程導(dǎo)出以位移、電勢及其一階導(dǎo)數(shù)為狀態(tài)變量的狀態(tài)方程，由邊界條件求得功能梯度壓電圓板軸對稱自由振動問題的精確解。Horgan與Chan[10]由平衡方程出發(fā)，假設(shè)彈性模量沿徑向冪律變化，研究了功能梯度旋轉(zhuǎn)圓盤及內(nèi)外壁受壓圓筒的軸對稱變形問題。

可是，對于面內(nèi)變剛度和小撓度功能梯度板彎曲問題的研究工作較少，這個問題的有限元解法也很少見。一般功能梯度材料的有限元解法多是將材料屬性也構(gòu)造為單元內(nèi)的插值函數(shù)形式[11-12]的梯度有限元法，或?qū)⑵浞謱泳鶆蚧嬎鉡13]，或者直接用商業(yè)有限元軟件進(jìn)行計算[14-15]。分層方法與有限元軟件方法類似，都是劃分為很多個均勻單元，但是各單元內(nèi)的材料參數(shù)不同，有一定的漸變關(guān)系，只有單元數(shù)量大才能不斷接近準(zhǔn)確結(jié)果。

本文嘗試從經(jīng)典薄板小撓度理論直接求解面內(nèi)變剛度的功能梯度圓板的薄板彎曲問題，即利用經(jīng)典薄板小撓度理論，首先推導(dǎo)出極坐標(biāo)下的控制方程，然后針對無孔圓板和含圓孔圓板兩種情況分別求解。而在有限元求解上，不使用梯度有限元法，而將材料參數(shù)的變化函數(shù)代替原來的常數(shù)，直接代入剛度矩陣計算，并提出了軸對稱圓板問題的簡易解法。通過將理論求解與有限元方法進(jìn)行對比，發(fā)現(xiàn)兩者吻合得非常好，表明本文的理論解具有一定的理論意義與應(yīng)用價值。

1 基于薄板小撓度理論的FGM圓板理論求解

1.1 基本方程

如圖1所示，以圓形薄板的中面為xy平面，圓板中心為原點(diǎn)，建立直角坐標(biāo)系。

假設(shè)板厚一定，半徑為a，上表面受均布力，不計體力，其楊氏模量沿圓板的半徑方向變化，泊松比的數(shù)值量級小，與其他參數(shù)相比波動較小，因此其變化可以忽略不計，視為常數(shù)，則薄板的彎曲剛度可以寫為如下形式：

由于考慮圓板問題，采用極坐標(biāo)系更便于計算。極坐標(biāo)系下坐標(biāo)與直角坐標(biāo)系下坐標(biāo)（x，y）的關(guān)系為：

在本文研究的問題中，材料的楊氏模量僅沿半徑方向變化，泊松比視為常數(shù) μ，則應(yīng)變 ερ、εφ、τρφ由應(yīng)力 σρ、σφ、τρφ表示為

由經(jīng)典板殼理論中的薄板小撓度彎曲理論假設(shè)[16]，可得ρ和φ方向的位移uρ、uφ：

其中：z為沿厚度方向的坐標(biāo)，w（ ρ， ）φ 為z方向的位移，即撓度。

本文考慮軸對稱圓板受軸對稱載荷的情況，顯然撓度函數(shù)也應(yīng)為軸對稱，僅與半徑有關(guān)，與角度無關(guān)：

因此，uφ=0。由幾何方程可得應(yīng)變與位移的關(guān)系為

由本構(gòu)方程（3）和（6）式得到應(yīng)力的表達(dá)式：

在薄板的每單位厚度上，將應(yīng)力分量σρ、σφ沿厚度方向積分合成為彎矩Mρ、Mφ，應(yīng)力分量τρφ合成為扭矩Mρφ，注意到積分過程中E（）ρ是面內(nèi)的函數(shù)，與z方向的積分無關(guān)，于是得到：

接下來由內(nèi)力推導(dǎo)平衡方程，取出薄板的一個微分塊的中面，如圖2所示，以微分塊中心的切向線為矩軸，可列出力矩平衡方程：

其中：Fsρ為橫向剪力，將內(nèi)力公式（8）代入（9）式得到：

1.2 無孔圓板求解

假設(shè)無孔圓形薄板上表面受均布壓力q0，取出半徑為ρ的中間部分的薄板，由平衡條件可得：

將D（ρ）代入薄板控制方程（10）中得到：

將（13）式進(jìn)行無量綱化，令 r=ρ/a得

注意到上述方程（14）是一個二階變系數(shù)常微分方程，系數(shù)項與r的級數(shù)項有關(guān)，因此考慮用級數(shù)展開法求解，右端項采用泰勒級數(shù)展開：

經(jīng)過整理，（14）式按級數(shù)形式可以寫成

對于0≤r≤1均成立，所以r的各次項前的系數(shù)應(yīng)相等，即可求出各項系數(shù)的關(guān)系。接下來從數(shù)學(xué)上證明此級數(shù)具有收斂性。

根據(jù)方程（17）可以求得如下結(jié)果：

綜上，該項級數(shù)是絕對收斂的。

而關(guān)于各個系數(shù)的求解，需要考慮邊界條件，圓板周邊設(shè)為固支邊，應(yīng)該滿足邊界位移和轉(zhuǎn)角為零：

選取各項參數(shù)如下：

將上述參數(shù)代入方程，求解級數(shù)項中的系數(shù)。

當(dāng)m=0，楊氏模量為常數(shù)，E=E0，即退化到了普通的圓形薄板，在經(jīng)典彈性力學(xué)中有精確解：

代入數(shù)值求解得到

利用本文的方法求解出的結(jié)果與之完全一致。

1.3 中心含圓孔圓板求解

假設(shè)在該圓板的中心位置有一個與圓板同心的圓孔，內(nèi)徑為b，圓板半徑仍為a，厚度為δ，上表面受橫向均布載荷q0，圓板仍會發(fā)生軸對稱彎曲?？刂品匠膛c之前的一致：

考慮與無孔圓板相同的邊界條件，外邊界固支，內(nèi)邊界自由。內(nèi)外兩個邊界，四個邊界條件為：

取出半徑為ρ的中心部分，內(nèi)邊界自由，無剪力作用，根據(jù)z方向的力平衡方程有

（23）式右邊項展開為如下形式：

考慮到當(dāng)下的研究對象是環(huán)形區(qū)域，因此對于上述控制方程的求解與無孔圓板不同，采用洛朗級數(shù)求解，不僅包括正次項系數(shù)，也包括負(fù)次項系數(shù)。將轉(zhuǎn)角函數(shù)設(shè)為如下形式：

將（24）式代入方程（23）得：

略去零系數(shù)項，轉(zhuǎn)角函數(shù)（24）可以重新寫為

則撓度函數(shù)為

注意到其中需要求解的未知系數(shù)僅為a1，c0兩項，而條件中含有四個邊界條件（20）-（21）式，分別為：

顯然方程數(shù)大于未知數(shù)個數(shù)，需要對解進(jìn)行優(yōu)化，本文采用最小二乘法進(jìn)行優(yōu)化。若僅將a1，c0代入上述邊界條件求解，顯然很難完全滿足邊界條件，大大降低了求解精度，因此結(jié)合（26）式對所有系數(shù)進(jìn)行優(yōu)化。假設(shè)未知系數(shù)為a-1～an，c0。首先由（29b-29d）式求得：

代入（26）式即可得到關(guān)于系數(shù)a2～an的n+1個方程，通過數(shù)學(xué)軟件Maple的幫助使用最小二乘法求得系數(shù) a2～an，然后由（29a-29d）求得 a-1、a0、a1和 c0。 代入數(shù)值求解：

用上一節(jié)的方法同樣可以驗證級數(shù)收斂，由圖3（b）可見，取級數(shù)展開的前十項系數(shù)即可。

2 有限元求解

本文研究的是面內(nèi)變剛度的FGM軸對稱圓形薄板受軸對稱載荷的情況，因此本節(jié)用軸對稱單元進(jìn)行有限元編程計算。在實(shí)際求解過程中，只需要沿著半徑方向劃分單元，求得半徑上各節(jié)點(diǎn)的撓度和轉(zhuǎn)角值即可代表整個圓板的變形。在經(jīng)典薄板小撓度彎曲理論中，忽略剪切變形，所以采用的單元為平面三節(jié)點(diǎn)環(huán)形單元。為與第一節(jié)結(jié)果對比，材料參數(shù)的變化函數(shù)一致，如第一節(jié)（11）式所示。

采用圓柱坐標(biāo)系，半徑方向為ρ坐標(biāo)，軸向為z坐標(biāo)，由于軸對稱無需考慮角度坐標(biāo)φ。

設(shè)w為z方向的位移 （撓度），a為圓板半徑，圓板受均布載荷q0。將各參數(shù)無量綱化，令：

現(xiàn)在來考察位移模式的收斂性：

（1）完備性：薄板的位移由中面位移決定，位移W中含有的常數(shù)項代表了中面沿z軸的剛體位移，一次項代表了繞半徑方向的剛體轉(zhuǎn)動，由于軸對稱，不考慮與角度相關(guān)的轉(zhuǎn)動項。

（2）協(xié)調(diào)性：板殼問題的泛函包含撓度的二階導(dǎo)數(shù)，因此插值函數(shù)應(yīng)具有C1連續(xù)性，即位移W及其一階導(dǎo)數(shù)在單元交界處連續(xù)。這里的環(huán)形單元，兩個單元交界為前一單元的終點(diǎn)和后一單元的起點(diǎn)，徑向坐標(biāo)一致，顯然求得的位移W和一階導(dǎo)數(shù)均為常數(shù)，完全滿足位移連續(xù)條件。綜上可知，該位移模式能夠滿足有限元收斂條件。

將單元的三節(jié)點(diǎn)分別代入計算：

由于小撓度彎曲理論中所有變量都是撓度的函數(shù)，所以只需要得到W的形式即可求出所有變量。由上三式將W重新表示為：

按照基爾霍夫薄板小撓度彎曲假設(shè)，在軸對稱圓柱坐標(biāo)系下，εz、ερz、εzφ、ερφ均為零，應(yīng)變?yōu)?/p>

由本構(gòu)方程可得應(yīng)力為：

單元剛度矩陣為：

單元載荷矩陣為：

將單元剛度矩陣和單元載荷矩陣組裝為總剛度矩陣和總載荷矩陣后，還需要設(shè)置邊界條件。實(shí)心圓板外邊界固支，因此設(shè)置最后一個單元的終點(diǎn)撓度、轉(zhuǎn)角為零，并考慮到對稱性，設(shè)置圓心處的轉(zhuǎn)角為零。按照這個條件，采用置大數(shù)法，將總剛度矩陣中對應(yīng)的對角線項設(shè)為大數(shù)，如1e30，從而在求解時使對應(yīng)的自由度值為零。空心圓板外邊界固支，內(nèi)邊界自由，只需設(shè)置最后一個單元的終點(diǎn)撓度、轉(zhuǎn)角為零，對應(yīng)總剛度矩陣對角項置大數(shù)，內(nèi)邊界無需外加條件。

對于以上各式的積分，采用六點(diǎn)高斯積分求解以保證準(zhǔn)確性，然后采用高斯消去法求解線性方程組。解得各節(jié)點(diǎn)的位移值W、θr后，為與第一節(jié)結(jié)果對比，取撓度w=as·W，θ=-θr。

為了驗證本文所提出的有限元法的準(zhǔn)確性，在（34）式中設(shè)m=0，則FGM板退化為均質(zhì)圓板，該種情況已有經(jīng)典解，求解后對比撓度曲線如圖5所示，實(shí)線為經(jīng)典解，“■”代表有限元解，顯然結(jié)果完全吻合，證明了本文有限元方法的可信性。

3 討 論

3.1 有限元結(jié)果與理論計算結(jié)果比較

當(dāng)m=1時，為FGM材料圓板，圖6和圖7分別為實(shí)心圓板和含圓孔圓板的結(jié)果對比。實(shí)線為第一節(jié)理論求解的結(jié)果，實(shí)心方塊代表第二節(jié)有限元的計算結(jié)果。

由圖6可見，利用兩種方法求解的實(shí)心圓板的撓度曲線和轉(zhuǎn)角曲線完全吻合。圖7中，圓板中心含圓孔時，從內(nèi)邊界至外邊界劃分單元，內(nèi)邊界自由，不加另外的邊界條件，孔洞半徑和圓板尺寸均與第一節(jié)相同。兩種方法求得的曲線趨勢一致，在內(nèi)邊界處的誤差最大，隨著半徑的增大，誤差逐漸減小。在內(nèi)邊界r=0.3 m處，有限元的計算結(jié)果為撓度值3.03E-4m，轉(zhuǎn)角值為5.55E-4，而利用前文的理論求解得撓度值為3.113E-4m，轉(zhuǎn)角值為5.92E-4，撓度誤差約為2.73%，轉(zhuǎn)角誤差約為6.6%。李堯臣的論文[15]中計算了面內(nèi)變剛度矩形板，并利用商業(yè)有限元軟件計算，誤差約為10%，本文的有限元方法更加準(zhǔn)確，理論解與之的誤差可能源于1.3節(jié)中因超定方程采用了解的優(yōu)化，無法完全保證所有方程精確滿足。

3.2 含孔圓板與不含孔圓板比較

3.3 楊氏模量變化對薄板變形的影響

顯然孔周的彎矩變化與圓板楊氏模量的變化函數(shù)有關(guān)。上例中設(shè)定E=E0e-ρ/a，從圓心開始沿半徑向外楊氏模量遞減。改變楊氏模量的變化規(guī)律，可以有效減小薄板的變形和孔邊的彎矩。

圖9和圖10分別為在m取不同的值，從1變到-1時撓度和彎矩變化。由圖線可以明顯看到當(dāng)m從1開始不斷減小，撓度逐漸減小，薄板的彎曲變形程度變小，孔邊的彎矩逐漸減小，整板的彎矩都有明顯的減小。m=0時為均質(zhì)板的特殊情況，而m變?yōu)樨?fù)數(shù)，楊氏模量沿半徑方向增大，撓度和彎矩比均質(zhì)板的情況更進(jìn)一步減小。由此可見，適當(dāng)改變功能梯度材料的結(jié)構(gòu)變化規(guī)律，可有效改善孔邊及整個結(jié)構(gòu)的受力情況。

4 結(jié) 論

本文討論了楊氏模量隨半徑呈指數(shù)函數(shù)形式變化時，圓形薄板和中心含圓孔的圓形薄板受均布載荷發(fā)生軸對稱彎曲的級數(shù)解和有限元解。與前人的方法對比，本文所用的方法更具有普適性，對于任意楊氏模量的變化均可求解，且不只局限于面內(nèi)變參數(shù)的功能梯度材料，對于連續(xù)變厚度的平板亦能適用，求解過程更直接可行，計算結(jié)果的精確度也有一定保證。