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

摘 要：利用Reissner中厚板理論邊界元、Mindlin中厚板理論有限元和三維彈性力學有限元分析了各種厚跨比的板的位移和應力，以確定利用三維理論、中厚板理論和薄板理論分析板時厚跨比的適用范圍以及精確程度，為分析板時選擇采用薄板理論、中厚板理論還是三維理論提供理論依據；同時也驗證了采用縮減積分的Mindlin中厚板理論有限元法緩解了剪切鎖死現象.

關鍵詞：邊界元法；有限元法；中厚板；薄板；厚跨比

中圖分類號：TU470.3 文獻標識碼：A

Research on the Applicable Range and Accuracy of Moderately Thick Plate Theory

LONG Shu-yao，JIANG Chen

(State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body， Hunan Univ， Changsha， Hunan 410082，China)

Abstract: The displacements and stresses for various thickness-to-span ratios were obtained using the boundary element method (BEM) for Reissner’s moderately thick plates， and the finite element method (FEM) for Mindlin’s moderately thick plates and three dimensional problems.The applicable range and accuracy among the three dimensional theories， the moderately thick plate theory and thin plate theory were defined. A theory basis and reference for choosing the thin plate theory， the thick plate theory or the three dimensional theory to analyze plates were provided. The results have shown that shear locking is deleted when the reduced integration is used in the finite element of plates and shells.

Key words:boundary element method; finite element method; moderately thick plate; thin plate; thickness-to-span ratio

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

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

1 Reissner中厚板理論

Reissner中厚板理論是考慮了橫向剪切變形和擠壓變形的理論.考慮笛卡爾坐標系(x1，x2，x3)，其中x1，x2軸在板中面內，x3則垂直于板中面，如圖1所示.在板的上下表面切應力為零而正應力為分布載荷集度，即x3=±h/2，σ33=±q/2

各應力在厚度方向上的變化關系[1]如下：

σij=12x3h3Mij，

σi3=32h1－2x3h2Qi，

σ33=142x3h3－2x3h2q.(1)

式中：i，j=1，2；Mij為單位長度彎矩；Qi為剪力.為了求解Reissner中厚板控制方程的基本解，在此假定體力在厚度方向上也按應力一樣變化，即

Pi=Bi(x1，x2)12x3h3，

P3=B3(x1，x2)32h1－2x3h2.(2)

由以上假定，Reissner中厚板的控制方程變為：

Mij，j－Qi=0，

Qi，i+q=0. (3)

其中的彎矩、扭矩和剪力用轉角和位移表示為：

Mij=D(1－ν)2βi，j+βj，i+2ν1－νβk，kδij+

ν1－νqλ2δij，

Qi=D1－ν2λ2(βi+w，i).(4)

式中：D=Eh3/［12(1－ν2)］為板的彎曲剛度；ν為泊松比；λ=10/h是厚度參數.式(4)中轉角和位移為沿板厚的平均值，即：

βi=∫h/2－h/212x3hvidx3，

w=∫h/2－h/232h1－2x3h2v3dx3.(5)

式中：vi和v3為坐標軸方向的位移分量.

為了簡便，將βi(i=1，2)和w用ui(i=1，2，3)表示，則對于區域Ω的邊界Γ，位移和廣義面力的邊界條件可以表示為：

ui=ui，在Su上； ti=ti，在St上.

其中S=Su∪St，ti為廣義面力.假定邊界上某點外法線方向余弦為ni=cos (n，xi)，則廣義面力為

ti=Mijnj，

t3=Qini.(i，j=1，2) (6)

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

分方程

基本解可由式(3)中體力項定義為一塊無限大中厚板作用單位橫向集中力或單位彎矩求得，具體推導見文獻[3].僅列出關于位移和面力的基本解：

u*ij=12πD{21－νB(z)δij－A(z)r，ir，j－

12δij(ln z－12)－12r，ir，j}，

u*i3=－u*3i=12πD12rr，i(ln z－12)，

u*33=－12πD1λ221－νln z－14z2(ln z－1).(7)

t*ij=－14πr［(4A(z)+2zK1(z)+1－ν)(δijr，nnj)+

(4A(z)+1－ν)r，jni－2(8A(z)+2zK1(z)+

1－ν)r，ir，jr，n］，

t*i3=λ22π［B(z)ni－A(z)r，ir，n］，

t*3i=－1－ν8π21+ν1－νln z－1ni+2r，ir，n，

t*33=－12πrr，n.(i，j=1，2)(8)

上兩式中：

A(z)=K0(z)+2zK1(z)－1z；

B(z)=K0(z)+1zK1(z)－1z；

z=λr；

r=［（x1－ξ1)2+(x2－ξ2)2］1/2；

r，i=(xi－ξi)/r；

r，n=r，ini.

其中K0(z)和K1(z)為修正貝塞爾函數.

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

cijuj+∫Γt*ijujdΓ=∫Γu*ijtjdΓ+

∫Ωu*i3－ν(1－ν)λ2u*ik，kqdΩ，

i，j=1，2，3；k=1，2. （9）

式中：cij是一個與邊界形狀有關的系數，對于光滑邊界cij=δij/2，對于有角點和尖點的邊界法線不連續的邊界，可以解析計算cij，但更方便的是間接計算cij的值.當整個邊界無外力和載荷作用時，式（9）變為：

cijuj+∫Γt*ijujdΓ=0. （10）

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

cia(ξ)=－∫Γ［t*ia(ξ，x)+(ξa－xa)t*i3(ξ，x)］dΓ，

ci3(ξ)=－∫Γt*i3(ξ，x)dΓ， （11）

i=1，2，3;α=1，2.

當均布荷載情況下，可以應用散度定理把式(9)中的域積分簡化為邊界積分.

數值方法求解式(9)后，未知邊界位移和邊界力可以全部求出.令cij=δij可以求出內部點的位移，而內部點的內力解，需要對式(9)關于源點ξ求微分，有如下形式：

Pij=ν(1－ν)λ2qδij+∫Γ［Dijktk－Sijkuk+qwij］dΓ. （12）

式中:i，j≠3.P的分量為Pij=Mij，P3j=Qj.其中wij，Dijk和Sijk見文獻[6]中附錄Ⅰ.

3 數值計算

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

HU=GT+D. （13）

式中：H，G為影響矩陣；D為對應于載荷q的邊界積分列向量.代入已知的邊界條件后，將其中的已知項移到方程等號右邊與D向量合并，未知項移到等號左邊，得到標準代數方程組

AX=B.（14）

采用標準高斯消去法求解上述方程組，即可求得未知量.之后可以利用式（9）計算內部任意一點的位移，利用式（12）求內部點的內力.

由于式中的奇異性，在計算奇異積分時，采用了三次多項式非線性變換[7].此種變換可以在奇異點附近集中較多的高斯點，因此對奇異積分采用標準高斯積分能達到相當高的精度.

對于存在角點和尖點等邊界法線不連續的情況，采用雙節點或部分不連續單元.雙節點單元中2個節點有不能在同一方向上給定位移的限制，而不連續單元沒有這種限制.

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

4 算 例

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

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

利用本文編制的Reissner中厚板邊界元程序計算了在不同厚跨比(h/a=0.001，0.005，0.01，0.05，0.1，0.2，0.5)下的板中心點的撓度w，板中心點下表面應力σx和板固支邊中點的彎矩Mx.并與利用有限元軟件ABAQUS中對應的三維有限元解和板殼有限元解，還有部分薄板解析解進行對比.

利用ABAQUS的三維8節點單元、二維4節點板殼單元和2節點邊界單元進行分析計算，其中以三維8節點有限元分析的結果作為基準，以說明二維4節點板殼單元和2節點邊界單元計算結果的準確性.

圖2 四邊固支矩形板

Fig.2 Four edges clamped rectangular plate

在三維有限元、板殼有限元和中厚板邊界元計算中，取中面中心點的撓度和中心點板的下表面點的應力σx作為比較對象.因為三維單元只能得到節點的位移和應力，要想得到截面上的內力時需要進行合成計算，為此取中心點板的下表面點的應力σx作為比較對象.板殼有限元分析結果既可輸出內力，也可輸出應力，但由于邊界元只能輸出內力，所以需要利用彎矩Mx與應力σx的關系求得：

σx=12Mxh3z.（15）

式中: z是應力計算點到中面的距離.

為了更好、更精確地對比各方法的計算結果，對Reissner中厚板邊界元模型每邊劃分16個線性單元，整個板共64個單元；對板殼有限元模型同樣每邊劃分16個線性單元，共16×16個單元；三維有限元模型則在板面內每邊劃分16個線性單元，厚度方向上網格的劃分是依照三維單元盡量要成正六面體的準則來進行的，這樣能保證作為基準的三維實體有限元結果的精度.

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

lg(h/a)

從圖4可看到，在厚跨比較小(h/a<0.2)時，邊界元、板殼有限元和三維有限元計算出的σx十分接近，但是在h/a>0.2后，誤差越來越大.但是邊界元解精度比板殼有限元的解精度高，證明了邊界元法在理論上比有限元法應力(內力)的精度高.綜合各種厚跨比撓度和應力的誤差，中厚板與厚板的厚跨比分界線一般取0.2為宜，厚跨比超過了0.2，無論是撓度還是應力的誤差都會急劇增加.同時，從圖3和圖4還可以看出：利用縮減積分技術的Mindlin板殼有限元，在計算非常薄的板時大大緩解了剪切鎖死現象.

lg(h/a)

圖5還給出了固支邊中點處的彎矩Mx.在計算彎矩時，三維實體單元選取20節點二次單元，因為二次單元的位移是二次分布的，但是應力是線性分布的，所以二次單元積分出的彎矩應該更精確.從圖5可以看出邊界元解比板殼有限元的解精確得多，而且在薄板階段和解析解十分靠近，充分說明了中厚板邊界元解對薄板有很好的精度.邊界元解和有限元三維實體單元解有差異，應當是來源于彎矩計算時積分的誤差，因為三維單元只能輸出應力，想獲得彎矩還需對應力進行積分，而節點并不是在高斯積分點上，只能使用梯形法近似積分，致使精度較差.

4.2 受均布載荷作用的四邊簡支矩形板

如圖6所示四邊簡支矩形板，所有數據同例4.1.

圖6 四邊簡支矩形板

Fig.6 Simple-supported rectangular plate

板殼有限元仍然采用雙線性板殼單元，而中厚板邊界元采用2節點線性單元.另外單元劃分準則也跟算例4.1相同，邊界元在每邊上劃分16個單元，整板共64個單元；板殼有限元也在每邊劃分16個單元.

lg(h/a)

圖7 各厚跨比下四邊簡支板中點撓度

Fig.7 Central deflection of simple-supported plate for several thickness-to-span ratios

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

lg(h/a)

圖8 各厚跨比下四邊簡支板中心處彎矩的誤差

Fig.8 Central errors of Mx of simple-supported plate

for several thickness-to-span ratios

從圖8中可以看出，無論是邊界元還是板殼有限元的解都十分接近三維有限元的解.從圖中可以看到，在薄板階段(h/a<0.1)，邊界元解和薄板理論十分接近，和三維有限元解、板殼有限元解有些許誤差，因為板殼有限元解的內力精度要低于邊界元解.

5 結 論

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

另外，從兩個算例可以得到結論，薄板和中厚板的分界線一般取厚跨比等于0.1，此后薄板理論解和中厚板理論解的誤差越來越大；而中厚板和厚板的分界線應當在厚跨比等于0.2處，因為此后中厚板理論解均和三維彈性力學解有相當大的誤差.同時說明了中厚板理論適用于薄板情況，并驗證了利用縮減積分技術的Mindlin板殼有限元，在計算非常薄的板時大大緩解了剪切鎖死現象.

參考文獻

［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]．北京：中國科學文化出版社，2002：214-253.

LONG Shu-yao. Introduction of boundary elements method [M]. Beijing:China Scientific Cultural Publishing House， 2002:214-253.（In Chinese)