中厚板理論的適用范圍和精確程度的研究＊

龍述堯，姜 琛

（湖南大學(xué) 汽車車身先進(jìn)設(shè)計(jì)制造國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410082）

板和殼在生產(chǎn)建設(shè)中被廣泛使用，隨著科學(xué)技術(shù)的發(fā)展，在某些領(lǐng)域中，除了傳統(tǒng)的薄板（殼）外，還使用了各種各樣的特殊板殼結(jié)構(gòu).例如，氣冷核反應(yīng)堆的預(yù)應(yīng)力混凝土壓力容器，航空航天中的各種夾層板（殼）和復(fù)合材料板（殼）等.在板殼問(wèn)題中，工程上廣泛使用的是Kirchhoff－Love薄板理論，Re－issner中厚板理論[1]和 Mindlin中厚板理論[2].在工程實(shí)際問(wèn)題中，由于結(jié)構(gòu)的幾何形狀、邊界條件和荷載的復(fù)雜性，使用解析方法求解幾乎不可能，大多使用各種數(shù)值方法進(jìn)行分析，其中廣泛使用的是有限元法.但是有限元法需域內(nèi)劃分單元、應(yīng)力解精度差、計(jì)算成本高，所以在20世紀(jì)70年代，計(jì)算機(jī)發(fā)展未能滿足有限元法的需求時(shí)，許多研究者開始研究邊界積分方程法 （Boundary Integration Equation Method，BIEM），提出了一種只需在求解域邊界上劃分單元的新方法，由 C.A.Brebbia[3]首先命名為邊界元法 （Boundary Element Method，BEM）.之后，邊界元法被應(yīng)用到各個(gè)領(lǐng)域，其中F.Vander Weeёn[4]和 S.Y.Long，C.A.Brebbia，J.C.F.Telles[5]研究了邊界元法在Reissner理論中厚板中的應(yīng)用.

本文采用Reissner中厚板理論邊界元程序、Mindlin中厚板理論有限元程序和三維彈性力學(xué)有限元程序計(jì)算各種厚跨比下四邊固支和四邊簡(jiǎn)支矩形板，并進(jìn)行對(duì)比，來(lái)說(shuō)明三維理論、中厚板理論和薄板理論分析板時(shí)厚跨比的適用范圍以及精確程度，為采用何種理論來(lái)分析板提供理論依據(jù)和參考；同時(shí)也驗(yàn)證了采用縮減積分的Mindlin中厚板有限元法消除了剪切鎖死現(xiàn)象.

1 Reissner中厚板理論

Reissner中厚板理論是考慮了橫向剪切變形和擠壓變形的理論.考慮笛卡爾坐標(biāo)系（x1，x2，x3），其中x1，x2軸在板中面內(nèi)，x3則垂直于板中面，如圖1所示.在板的上下表面切應(yīng)力為零而正應(yīng)力為分布載荷集度，即x3＝±h／2，σ33＝±q／2.

圖1 中厚板幾何形狀和坐標(biāo)系Fig.1 Geometry shape and coordinate of moderately thick plate

各應(yīng)力在厚度方向上的變化關(guān)系[1]如下：

式中：i，j＝1，2；Mij為單位長(zhǎng)度彎矩；Qi為剪力.為了求解Reissner中厚板控制方程的基本解，在此假定體力在厚度方向上也按應(yīng)力一樣變化，即

由以上假定，Reissner中厚板的控制方程變?yōu)椋?/p>

其中的彎矩、扭矩和剪力用轉(zhuǎn)角和位移表示為：

式中：D＝Eh3／[12（1－ν2）]為板的彎曲剛度；ν為泊松比；是厚度參數(shù).式（4）中轉(zhuǎn)角和位移為沿板厚的平均值，即：

式中：vi和v3為坐標(biāo)軸方向的位移分量.

為了簡(jiǎn)便，將βi（i＝1，2）和w用ui（i＝1，2，3）表示，則對(duì)于區(qū)域Ω的邊界Γ，位移和廣義面力的邊界條件可以表示為：

其中S＝Su∪St，ti為廣義面力.假定邊界上某點(diǎn)外法線方向余弦為ni＝cos（n，xi），則廣義面力為

2 Reissner中厚板邊界元基本解和邊界積分方程

基本解可由式（3）中體力項(xiàng)定義為一塊無(wú)限大中厚板作用單位橫向集中力或單位彎矩求得，具體推導(dǎo)見(jiàn)文獻(xiàn)[3].僅列出關(guān)于位移和面力的基本解：

上兩式中：

其中K0（z）和K1（z）為修正貝塞爾函數(shù).

使用Betti功互等定理、格林函數(shù)或加權(quán)殘值法，可以得到Reissner中厚板的邊界積分方程：

式中：cij是一個(gè)與邊界形狀有關(guān)的系數(shù)，對(duì)于光滑邊界cij＝δij／2，對(duì)于有角點(diǎn)和尖點(diǎn)的邊界法線不連續(xù)的邊界，可以解析計(jì)算cij，但更方便的是間接計(jì)算cij的值.當(dāng)整個(gè)邊界無(wú)外力和載荷作用時(shí)，式（9）變?yōu)椋?/p>

上式的非零解為板的3種獨(dú)立剛體位移的任意組合，把這些位移代入式（10）有

當(dāng)均布荷載情況下，可以應(yīng)用散度定理把式（9）中的域積分簡(jiǎn)化為邊界積分.

數(shù)值方法求解式（9）后，未知邊界位移和邊界力可以全部求出.令cij＝δij可以求出內(nèi)部點(diǎn)的位移，而內(nèi)部點(diǎn)的內(nèi)力解，需要對(duì)式（9）關(guān)于源點(diǎn)ξ求微分，有如下形式：

式中：i，j≠3.P的分量為Pij＝Mij，P3j＝Qj.其中wij，Dijk和Sijk見(jiàn)文獻(xiàn)[6]中附錄Ⅰ.

3 數(shù)值計(jì)算

把所考慮的板的邊界劃分為M個(gè)單元，以每個(gè)節(jié)點(diǎn)為源點(diǎn)并在M個(gè)單元上積分，離散式（9），并將其中域積分轉(zhuǎn)換為邊界積分，則可得線性方程組

式中：H，G為影響矩陣；D為對(duì)應(yīng)于載荷q的邊界積分列向量.代入已知的邊界條件后，將其中的已知項(xiàng)移到方程等號(hào)右邊與D向量合并，未知項(xiàng)移到等號(hào)左邊，得到標(biāo)準(zhǔn)代數(shù)方程組

采用標(biāo)準(zhǔn)高斯消去法求解上述方程組，即可求得未知量.之后可以利用式（9）計(jì)算內(nèi)部任意一點(diǎn)的位移，利用式（12）求內(nèi)部點(diǎn)的內(nèi)力.

由于式中的奇異性，在計(jì)算奇異積分時(shí)，采用了三次多項(xiàng)式非線性變換[7].此種變換可以在奇異點(diǎn)附近集中較多的高斯點(diǎn)，因此對(duì)奇異積分采用標(biāo)準(zhǔn)高斯積分能達(dá)到相當(dāng)高的精度.

對(duì)于存在角點(diǎn)和尖點(diǎn)等邊界法線不連續(xù)的情況，采用雙節(jié)點(diǎn)或部分不連續(xù)單元.雙節(jié)點(diǎn)單元中2個(gè)節(jié)點(diǎn)有不能在同一方向上給定位移的限制，而不連續(xù)單元沒(méi)有這種限制.

根據(jù)前述離散化的公式，編制了Reissner中厚板邊界元的FORTRAN程序[8].

4 算 例

采用了2個(gè)具有代表性的算例來(lái)研究薄板和中厚板理論適用范圍和精確程度.算例為：1）受均布荷載作用的四邊固支矩形板，如圖2所示；2）受均布荷載作用的四邊簡(jiǎn)支矩形板.兩算例的長(zhǎng)、寬取a＝b＝1，橫向荷載集度q＝1，板的彎曲剛度D＝1，泊松比ν＝0.3.

4.1 受均布荷載作用的四邊固支矩形板

利用本文編制的Reissner中厚板邊界元程序計(jì)算了在不同厚跨比（h／a＝0.001，0.005，0.01，0.05，0.1，0.2，0.5）下的板中心點(diǎn)的撓度w，板中心點(diǎn)下表面應(yīng)力σx和板固支邊中點(diǎn)的彎矩Mx.并與利用有限元軟件ABAQUS中對(duì)應(yīng)的三維有限元解和板殼有限元解，還有部分薄板解析解進(jìn)行對(duì)比.

利用ABAQUS的三維8節(jié)點(diǎn)單元、二維4節(jié)點(diǎn)板殼單元和2節(jié)點(diǎn)邊界單元進(jìn)行分析計(jì)算，其中以三維8節(jié)點(diǎn)有限元分析的結(jié)果作為基準(zhǔn)，以說(shuō)明二維4節(jié)點(diǎn)板殼單元和2節(jié)點(diǎn)邊界單元計(jì)算結(jié)果的準(zhǔn)確性.

圖2 四邊固支矩形板Fig.2 Four edges clamped rectangular plate

在三維有限元、板殼有限元和中厚板邊界元計(jì)算中，取中面中心點(diǎn)的撓度和中心點(diǎn)板的下表面點(diǎn)的應(yīng)力σx作為比較對(duì)象.因?yàn)槿S單元只能得到節(jié)點(diǎn)的位移和應(yīng)力，要想得到截面上的內(nèi)力時(shí)需要進(jìn)行合成計(jì)算，為此取中心點(diǎn)板的下表面點(diǎn)的應(yīng)力σx作為比較對(duì)象.板殼有限元分析結(jié)果既可輸出內(nèi)力，也可輸出應(yīng)力，但由于邊界元只能輸出內(nèi)力，所以需要利用彎矩Mx與應(yīng)力σx的關(guān)系求得：

式中：z是應(yīng)力計(jì)算點(diǎn)到中面的距離.

為了更好、更精確地對(duì)比各方法的計(jì)算結(jié)果，對(duì)Reissner中厚板邊界元模型每邊劃分16個(gè)線性單元，整個(gè)板共64個(gè)單元；對(duì)板殼有限元模型同樣每邊劃分16個(gè)線性單元，共16×16個(gè)單元；三維有限元模型則在板面內(nèi)每邊劃分16個(gè)線性單元，厚度方向上網(wǎng)格的劃分是依照三維單元盡量要成正六面體的準(zhǔn)則來(lái)進(jìn)行的，這樣能保證作為基準(zhǔn)的三維實(shí)體有限元結(jié)果的精度.

圖3畫出了各種厚跨比下板中點(diǎn)的撓度.從圖3中可以看出在所有厚跨比情況下，3種方法的計(jì)算結(jié)果十分吻合，說(shuō)明了考慮剪切變形的Reissner中厚板理論是對(duì)三維彈性力學(xué)的較精確的簡(jiǎn)化.從圖中還可以看出，當(dāng)厚跨比大于0.2時(shí)，與三維單元的結(jié)果誤差明顯增大.

圖3 各厚跨比情況下四邊固支板中點(diǎn)撓度Fig.3 Central deflections of clamped plate for several thickness－to－span ratios

從圖4可看到，在厚跨比較?。╤／a＜0.2）時(shí)，邊界元、板殼有限元和三維有限元計(jì)算出的σx十分接近，但是在h／a＞0.2后，誤差越來(lái)越大.但是邊界元解精度比板殼有限元的解精度高，證明了邊界元法在理論上比有限元法應(yīng)力（內(nèi)力）的精度高.綜合各種厚跨比撓度和應(yīng)力的誤差，中厚板與厚板的厚跨比分界線一般取0.2為宜，厚跨比超過(guò)了0.2，無(wú)論是撓度還是應(yīng)力的誤差都會(huì)急劇增加.同時(shí)，從圖3和圖4還可以看出：利用縮減積分技術(shù)的Mindlin板殼有限元，在計(jì)算非常薄的板時(shí)大大緩解了剪切鎖死現(xiàn)象.

圖4 各厚跨比情況下四邊固支板中點(diǎn)下表面的σx誤差Fig.4 Central errors ofσxof clamped plate for several thickness－to－span ratios

圖5還給出了固支邊中點(diǎn)處的彎矩Mx.在計(jì)算彎矩時(shí)，三維實(shí)體單元選取20節(jié)點(diǎn)二次單元，因?yàn)槎螁卧奈灰剖嵌畏植嫉?，但是?yīng)力是線性分布的，所以二次單元積分出的彎矩應(yīng)該更精確.從圖5可以看出邊界元解比板殼有限元的解精確得多，而且在薄板階段和解析解十分靠近，充分說(shuō)明了中厚板邊界元解對(duì)薄板有很好的精度.邊界元解和有限元三維實(shí)體單元解有差異，應(yīng)當(dāng)是來(lái)源于彎矩計(jì)算時(shí)積分的誤差，因?yàn)槿S單元只能輸出應(yīng)力，想獲得彎矩還需對(duì)應(yīng)力進(jìn)行積分，而節(jié)點(diǎn)并不是在高斯積分點(diǎn)上，只能使用梯形法近似積分，致使精度較差.

圖5 各厚跨比下四邊固支板邊中點(diǎn)處的MxFig.5 CentralMxof clamped plate for several thikness－to－span ratios

4.2 受均布載荷作用的四邊簡(jiǎn)支矩形板

如圖6所示四邊簡(jiǎn)支矩形板，所有數(shù)據(jù)同例4.1.

圖6 四邊簡(jiǎn)支矩形板Fig.6 Simple－supported rectangular plate

板殼有限元仍然采用雙線性板殼單元，而中厚板邊界元采用2節(jié)點(diǎn)線性單元.另外單元?jiǎng)澐譁?zhǔn)則也跟算例4.1相同，邊界元在每邊上劃分16個(gè)單元，整板共64個(gè)單元；板殼有限元也在每邊劃分16個(gè)單元.

圖7 各厚跨比下四邊簡(jiǎn)支板中點(diǎn)撓度Fig.7 Central deflection of simple－supported plate for several thickness－to－span ratios

從圖7中可以看出，3種單元的四邊簡(jiǎn)支板中心的撓度的解十分接近.在厚板階段（h／a＞0.1），圖中清楚地表明了本文所用中厚板邊界元法比板殼有限元法更接近三維有限元的解，說(shuō)明了邊界元法精度比有限元高.

圖8 各厚跨比下四邊簡(jiǎn)支板中心處彎矩的誤差Fig.8 Central errors ofMxof simple－supported plate for several thickness－to－span ratios

從圖8中可以看出，無(wú)論是邊界元還是板殼有限元的解都十分接近三維有限元的解.從圖中可以看到，在薄板階段（h／a＜0.1），邊界元解和薄板理論十分接近，和三維有限元解、板殼有限元解有些許誤差，因?yàn)榘鍤び邢拊獾膬?nèi)力精度要低于邊界元解.

5 結(jié) 論

本文將Reissner中厚板理論運(yùn)用到邊界元法中，并通過(guò)數(shù)值方法將其編制成FORTRAN程序，通過(guò)幾個(gè)典型算例的計(jì)算，證明了邊界元法分析中厚板的程序的正確性.另外，通過(guò)與三維有限元、板殼有限元的對(duì)比，證明了Reissner理論是對(duì)三維彈性力學(xué)的正確簡(jiǎn)化，同時(shí)，中厚板邊界元法比板殼有限元法在撓度、內(nèi)力、應(yīng)力上具有更高的精度，更高的計(jì)算效率.因此，使用邊界元方法對(duì)中厚板進(jìn)行數(shù)值計(jì)算，是十分合適和有效的.

另外，從兩個(gè)算例可以得到結(jié)論，薄板和中厚板的分界線一般取厚跨比等于0.1，此后薄板理論解和中厚板理論解的誤差越來(lái)越大；而中厚板和厚板的分界線應(yīng)當(dāng)在厚跨比等于0.2處，因?yàn)榇撕笾泻癜謇碚摻饩腿S彈性力學(xué)解有相當(dāng)大的誤差.同時(shí)說(shuō)明了中厚板理論適用于薄板情況，并驗(yàn)證了利用縮減積分技術(shù)的Mindlin板殼有限元，在計(jì)算非常薄的板時(shí)大大緩解了剪切鎖死現(xiàn)象.

[1] REISSNER E.The effect of transverse shear deformation on the bending of elastic plate[J].J Appl Mech，1945，12：69－77.

[2] MINDLIN R D.Influence of rotatory inertia and shear on flexural motions of isotropic elastic plates[J].Journal of Applied Mechanics，1951，18：31－38.

[3] BREBBIA C A.The boundary element method for engineers[M].London：Pentech Press，1978：159－211.

[4] VANDER WEE?N F.Application of the boundary integral equation method to Reissner's plate model[J].International Journal for Numerical Methods in Engineering，1982，18（1）：1－10.

[5] LONG S Y，BREBBIA C A，TELLS J C F.Boundary element bending analysis of moderately thick plates[J].Engineering Analysis，1988，5（2）：64－74.

[6] KARAM V J，TELLES J C F.On boundary elements for Reissner’s plate theory[J].Engineering Analysis，1988，5（1）：21－28.

[7] TELLES J C F.A self－adaptive coordinate transformation for efficient numerical evaluation of general boundary element integral[J].Int J Numer Engineering，1987，24：959－973.

[8] 龍述堯.邊界單元法概論[M].北京：中國(guó)科學(xué)文化出版社，2002：214－253.LONG Shu－yao.Introduction of boundary elements method[M].Beijing：China Scientific ＆ Cultural Publishing House，2002：214－253.（In Chinese）