蘇海東，韓陸超，頡志強

(長江科學院 材料與結構研究所，武漢 430010)

1 研究背景

梁板殼結構分析通?；谝恍┍粚嶒烌炞C的基本變形假設(已歸納為具有很高模擬精度的近似理論)[1]：對于經典的薄梁板殼，有Euler-Bernoulli梁理論和Kirchhoff-Love板殼理論，要求原垂直于中面的橫截面在變形過程中始終保持為平面且仍垂直于中面，忽略橫向剪應變；對于需要考慮橫向剪應變的厚梁板殼，采用Timoshenko梁理論和Reissner-Mindlin板殼理論，只要求橫截面保持為平面。當然，前者也可當作后者在厚度變薄時的特殊情況。此外，沿厚度方向產生的直接應變可忽略。

從薄梁板殼假設出發(fā)，前人推導了不同于實體分析的控制方程[1]，最終都歸結為關于撓度的4階微分方程，由此帶來在有限元數(shù)值分析中要求插值函數(shù)滿足C1連續(xù)性(即導數(shù)連續(xù))的難題，至今都未能得到很好的解決[2-3]。

進一步對于曲梁和曲殼而言，由于存在初始曲率，其控制方程的推導本身就比較復雜(特別是曲殼[4])，因而一般不采用由控制方程建立積分弱形式的標準方式，而是用直梁或平板單元近似地模擬曲梁或曲殼，容易產生幾何誤差，進而帶來力學分析上的誤差。

筆者基于獨立覆蓋流形法和分區(qū)級數(shù)解思想，提出梁板殼結構數(shù)值計算的新方法[5-7]，針對一般的中面參數(shù)方程的幾何描述，將中面的精確幾何變化在應變計算中反映為中面局部坐標系的方向余弦關于整體坐標的導數(shù)運算，首次實現(xiàn)了精確幾何的曲梁和曲殼分析，但目前僅限于厚梁板殼假設[6-7]，還未在經典的薄梁板殼假設中實現(xiàn)，這是該方法的基礎理論部分尚未完成的工作。另外，對于各種中面的曲線或曲面描述，其相關的幾何公式都需要進行一遍復雜繁瑣的人工推導。

本文針對上述問題開展研究。首先簡要介紹獨立覆蓋流形法及分區(qū)級數(shù)解的思想，并討論其近似函數(shù)的C1連續(xù)性，然后以平面薄曲梁為例詳細給出計算公式，并簡述薄曲殼的計算過程，最后通過典型算例驗證方法的有效性。

2 獨立覆蓋流形法及分區(qū)級數(shù)解簡介

1991年，數(shù)值流形方法[8-9](以下簡稱流形法)的發(fā)明者石根華將現(xiàn)代數(shù)學的流形思想引入數(shù)值計算：如圖1(a)所示，一系列相互重疊的覆蓋用于分割求解域(要求覆蓋整體上包含求解域)，使整體分布復雜的物理場在每個覆蓋區(qū)域內簡單化，以利于局部近似函數(shù)的逼近。

從2011年開始，筆者在流形法中首次引入“獨立覆蓋”[10]，即覆蓋的獨立區(qū)域(非重疊部分)，如圖1(b)所示，每個覆蓋中的“藍色”區(qū)域就是獨立覆蓋，其近似函數(shù)通常采用多項式等完備級數(shù)。每個覆蓋中的“黃色”區(qū)域是該覆蓋與周邊覆蓋相重疊的部分，在此重疊區(qū)域內，采用單位分解函數(shù)φi(在流形法中又稱為權函數(shù))將各覆蓋中的級數(shù)Vi聯(lián)系起來[11]，以保證整體近似函數(shù)V的連續(xù)性,即

(1)

式中n為重疊的覆蓋數(shù)。

筆者進而提出獨立覆蓋流形法[12]，強調獨立覆蓋是計算分析的主要對象：如圖1(c)所示，獨立覆蓋占據一個覆蓋的主要面積，重點研究獨立覆蓋內的級數(shù)，覆蓋之間的重疊區(qū)域形成面積盡可能小(可只占覆蓋面積的1%)的窄“條形”，僅起到連接各覆蓋級數(shù)的作用，并以有限單元的線性形函數(shù)作為權函數(shù)，比如兩個覆蓋之間采用一維有限元的形函數(shù)[2]，即:

(2)

式中ξ為在條形厚度方向以條形中點為原點的局部坐標。這樣，由各覆蓋分區(qū)的級數(shù)聯(lián)合形成的整體近似函數(shù)，通過加權殘值法(如伽遼金法)逼近真實解，稱為基于流形思想的“分區(qū)級數(shù)解”(這是獨立覆蓋流形法的內核)，可作為微分方程的數(shù)值解答。

圖1 流形覆蓋示意圖Fig.1 Manifold covers

筆者在文獻[12]中討論了該方法的收斂性，指出：隨著各分區(qū)級數(shù)的階次升高，整體近似函數(shù)逐步逼近真實解；不僅場函數(shù)本身是收斂的，其導數(shù)也是收斂的；收斂性不受分區(qū)形狀的影響，也不受各覆蓋在條形重疊區(qū)域是否錯位連接的影響。因此發(fā)展出了基于任意網格(任意形狀、任意連接和任意加密)的數(shù)值計算[13]，以及無需人工參與的自動計算[14]。

采用分區(qū)級數(shù)解的一大優(yōu)勢是，可事先構造出滿足物理場整體或局部特性的級數(shù)形式。筆者提出的梁板殼數(shù)值計算新方法，對于多項式級數(shù)中包含厚度方向坐標的各項，僅在模擬板殼的薄膜位移(或梁的軸線向位移)時保留至厚度方向坐標的1階項，就可將厚梁板殼的基本變形假設固化在近似函數(shù)中[5-7]。同時在文獻[5]的直梁分析中，嘗試將所有含厚度方向坐標的級數(shù)項去掉，實現(xiàn)了細長梁(以下也稱為薄梁)的Euler-Bernoulli梁理論，對于其所要求的C1連續(xù)性，當時的解釋是“基于導數(shù)的收斂性，近似函數(shù)的導數(shù)趨向于連續(xù)的真實解”，本文補充說明如下。

圖2 2個覆蓋及其條形連接Fig.2 Two covers and a stripbetween them

在條形重疊區(qū)域內有

V=φ1V1+φ2V2。

(3)

不難驗證C0連續(xù)性。條形區(qū)域關于x方向的一階導數(shù)為

在完全收斂的情況下，在條形區(qū)域內由于真實解V*的連續(xù)性，有V1=V2→V*。即使在未完全收斂的情況下，考慮到條形區(qū)域很窄，則在達到一定階次時，條形處的V1和V2非常接近，即此處V1≈V2，因此以上情況都有

(5)

3 薄曲梁的分區(qū)級數(shù)表達及應變矩陣計算

如圖3所示，在整體直角坐標系x-y中，由曲線參數(shù)方程描述中面坐標(x0,y0)，建立沿中面曲線變化的局部正交坐標系xi-yi[6]，其中xi沿梁的軸線方向(中面曲線的切向)，yi沿厚度方向(中面曲線的法向)，xi軸在整體坐標系下的方向余弦為cosαxi、cosβxi，yi表示梁上的點到中面的距離，yi軸的方向余弦為cosαyi、cosβyi。

圖3 整體直角坐標系下的曲梁及其局部坐標系Fig.3 A curved beam and its local coordinates underthe global rectangular coordinate system

根據梁保持平截面的假設，局部坐標系下的兩個方向位移為[2]

(6)

其中，對于滿足Euler-Bernoulli理論的薄曲梁，有[4]

(7)

式中：一維多項式f1(xi)和f3(xi)分別表示梁中面的軸線方向和厚度方向的位移;f2(xi)=θ(xi)表示截面轉角；R為中面點處的曲率半徑；h為弧長微分系數(shù)。

因此只需考慮f1(xi)和f3(xi)，如圖4所示，對于多項式級數(shù)，在軸線和厚度方向去掉所有關于yi的項，只讓關于xi的一維級數(shù)(圖4中的虛線左側)參與計算(式(8)),就能通過式(6)和式(7)實現(xiàn)薄梁的變形假設，而無需推導薄梁控制方程。

(8)

式(6)—式(8)寫成矩陣形式,則有

(9)

分別記cosαxi、cosβxi為cosα、cosβ，同時考慮兩個局部坐標軸的垂直關系，將局部坐標下的位移轉換到整體坐標下的位移，即

(10)

并記其中的函數(shù)公式(對于第q項)為

二維應變子矩陣(對于第q項)為

(12)

其中，由鏈式法則，各項對整體坐標的偏導數(shù)為

(13)

形成3×3數(shù)組D。其中的數(shù)組運算(式(14)中j=1～3)為

(14)

然后D數(shù)組的每一項再與C數(shù)組的每一項按式(14)的方式進行組合。最后依次放入式(12)的對應位置，形成應變子矩陣Bq，并代入剛度矩陣公式計算單元剛度矩陣[6]。

對于連接覆蓋1和覆蓋2的條形區(qū)域，求關于整體坐標x或y的偏導數(shù)時，Bq矩陣每一項系數(shù)c還要與權系數(shù)φi(i=1,2)進行組合，形成應變子矩陣Biq，則有

(15)

其他關于積分方法等內容參見文獻[6]。另外，取R→∞和h=1可退化到計算細長直梁。

4 曲梁的幾何計算公式

(1)給定中面曲線的參數(shù)方程為

(16)

式中t為參數(shù)。取局部坐標xi=t，yi=r，r為曲梁上的各點到中面的距離。

(2)關于t的1—3階導數(shù)為：

(3)弧長的微分系數(shù)h表達式為

(18)

(4)曲率半徑R為

(19)

(5)方向余弦為:

(20)

(6)局部坐標關于整體坐標的導數(shù)。曲梁上的任一點為

用隱函數(shù)求導法分別對x和y求偏導，并求解方程組，令g1=x′t-r(cosβ)′t，g2=y′t+r(cosα)′t，g=g1cosα+g2cosβ，則：

(22)

綜上所述，在給定中面參數(shù)方程后，人工推導僅限于式(17)中相對簡單的各階求導，其他運算都由程序自動完成。因此，本文方法解決了以往在幾何公式推導上的難題，對于任意給定的中面參數(shù)方程都具有通用性。

另外，本文參考文獻[13]，采用邊界條的方式嚴格施加邊界條件。如式(23)所示的在邊界xi=xc處的邊界覆蓋級數(shù)(一般取為2階多項式即可，此處l為邊界條的寬度)，若施加簡支條件，則僅需令系數(shù)a0=b0=0，其余各項在xi=xc處自動為0，從而嚴格滿足邊界條件。若在xi處施加固支條件，考慮到式(7)中的f3(xi)的導數(shù)對ui的影響，還需令b1=0。

(23)

5 薄曲殼計算簡述

如圖5所示的曲殼中面，xi和yi為正交曲線坐標，zi表示殼體上的點(整體坐標為(x,y,z))到中面的距離[7]。給定中面參數(shù)方程為

(24)

圖5 整體直角坐標系下的曲殼及其局部坐標Fig.5 A curved shell and its local coordinates underthe global rectangular coordinate system

按薄板殼的Kirchhoff-Love理論，定義局部坐標系下位移覆蓋函數(shù)為

(26)

其中，對于正交的局部坐標系，兩個方向的轉角為[4]

式中：h1、h2分別為xi、yi兩個局部坐標的拉梅系數(shù)；R1、R2分別為兩個局部坐標方向的曲率半徑(對于正交的局部坐標系，本文為兩個主曲率半徑)[4]，其表達式為

(28)

其中：

(29)

將局部坐標系下的位移轉化到整體坐標系下的位移，并記其中的函數(shù)公式(對于第q項)為

(30)

式中：fq=fq(xi,yi)為關于(xi,yi)的第q項多項式；cosαxi、cosβxi、cosγxi分別是xi軸關于整體坐標x、y、z的方向余弦，其他以此類推。

薄曲殼計算的關鍵是三維應變子矩陣Bq[7]。整體的計算思路與上兩節(jié)的薄曲梁相似，只是過程更復雜一些：同樣要拆解成數(shù)組運算并進行組合，最后放入Bq矩陣的相應位置；在幾何運算中，需要事先推導中面方程關于局部坐標的1—3階導數(shù)公式(曲率半徑的微分需要3階導數(shù)[15])，之后的運算都由程序自動完成。除了式(28)和式(29)以外，其他幾何量的表達式見文獻[7]。

本文參考文獻[13]，將曲梁和曲殼的計算程序改寫成標準的二維獨立覆蓋及條形計算方式，而不采用以往基于強制約束形成獨立覆蓋及條形的方式[5-7]，這樣可以大幅減少計算量，并能實現(xiàn)任意網格劃分，以及通過邊界條嚴格施加邊界條件。

6 算 例

采用薄梁板殼假設，將文獻[6]和文獻[7]的曲梁和曲殼主要算例重新計算一遍。

6.1 常曲率的曲梁算例——圓形曲梁

如圖6所示，梁的矩形截面厚度為0.01 m。左、右兩端采用邊界條方式施加固端約束。法向集中力F=1 kN作用于梁的中點。彈性模量E=105kN/m2，泊松比μ=0。荷載作用點處的位移及與細密網格的有限元結果對比見表1。

圖6 圓形曲梁及其變形Fig.6 Circular curved beam and its deformation

表1 荷載作用點處的位移比較(圓形曲梁)Table 1 Comparison of the displacements of loadingaction point of the circular-curved beam

6.2 變曲率的曲梁算例——橢圓形曲梁

設橢圓的長軸a=2 m，短軸b=1 m，梁的矩形截面厚度為0.01 m。彈性模量E=108kN/m2，泊松比μ=0。如圖7所示，只用1個獨立覆蓋，左端全約束，右端分別承受水平力和豎向力F=1 kN。荷載作用點處的位移及其與細密網格的有限元結果對比見表2。

圖7 橢圓形曲梁及其變形Fig.7 Ellipse-curved beam and its deformation

表2 荷載作用點處的位移比較(橢圓形曲梁)Table 2 Comparison of the displacements of loadingaction point of an ellipse-curved beam

6.3 球面殼算例

計算如圖8(a)所示的第1象限內的球面殼，半徑為1 m，θ=π/4～π/2，φ=0～π/2。殼體厚度為0.01 m，左端(φ=0)和右端(φ=π/2)為固支約束，球面殼上表面承受法向均布壓力P=10 kN/m2。彈性模量E=106kN/m2，泊松比μ=0。圖8(b)顯示了用于計算球面殼參考解的劃分為5 000個平板單元的有限元網格。

圖8 球坐標及第1象限內的球面殼有限元網格Fig.8 Spherical coordinates and finite element meshesof a spherical shell in the first quadrant

本文在φ-θ的投影平面上分別采用3種網格：網格a——1個獨立覆蓋網格；網格b——均勻劃分為5×3(φ向×θ向)網格；網格c——均勻劃分為9×5網格。計算得到的各方向最大位移見表3，出現(xiàn)在(θ=π/2，φ=π/4)點處，均與細密網格的有限元解非常接近。以其中的網格a為例，本文方法的自由度為417，而文獻[7]采用厚殼假設時為669。

表3 球殼最大位移計算結果比較Table 3 Comparison of the maximum displacements ofthe spherical shell

將球殼兩端改為簡支約束，采用網格a計算，則9階情況下的最大位移u=v=-0.654 cm，w=-0.321 cm，而細密網格的有限元解為u=v=-0.658 cm，w=-0.320 cm。

程序也可用于計算平板。采用網格c，平板范圍為x=(0～π/2)m，y=(π/4～π/2)m，考慮四周簡支約束，其他條件不變。按薄板理論，采用4階多項式計算得到中點位移為0.460 2 m，與有限元參考解0.460 5 m非常接近。

與文獻[6]和文獻[7]相比，上述各算例在兩種變形假設下的計算結果幾乎一致，本文基于薄梁板殼假設可節(jié)省自由度約30%(對于平板彎曲算例，還可去掉所有關于薄膜位移的自由度)，但單元計算要復雜一些。

7 結 語

本文在精確幾何曲梁和曲殼的前期研究基礎上，完成了基于Euler-Bernoulli梁理論和Kirchhoff-Love板殼理論的薄曲梁和曲殼分析，不僅補上了梁板殼數(shù)值計算新方法在基礎理論上的最后一塊拼圖，而且解決了幾何公式推導復雜的問題(人工只需推導中面方程關于參數(shù)坐標的最多3階導數(shù))。以下再次總結新方法的特點和優(yōu)勢：

(1)采用實體計算模式，只需設定各分區(qū)的多項式級數(shù)中含厚度方向坐標的級數(shù)項不參與計算(或僅在薄膜位移中保留厚度方向坐標的1階項)，就能模擬薄或厚梁板殼的2種基本假設，而無需梁板殼控制方程及推導相應的數(shù)值計算公式。

(2)在收斂性方面具備C1連續(xù)，且用厚梁板殼假設計算薄梁板殼時無閉鎖現(xiàn)象，也沒有通常實體計算中厚度遠小于其他尺度導致的數(shù)值病態(tài)。

(3)對于任意給定中面參數(shù)方程和厚度分布的曲梁和曲殼(目前對于薄曲殼假設情況，中面參數(shù)方程暫定為正交曲線坐標描述)，兩種基本假設情況均實現(xiàn)了精確幾何描述下的力學分析，避免了由直梁或平板單元近似模擬曲梁或曲殼造成的幾何誤差。

因此，基于獨立覆蓋“分區(qū)級數(shù)解”的梁板殼數(shù)值計算新方法不僅避免了現(xiàn)有方法的諸多難題，而且展現(xiàn)了在精確幾何曲梁和曲殼分析方面的獨特優(yōu)勢。近期將進行梁板殼與一般實體單元連接以及梁板殼相互連接的算例研究，然后再結合獨立覆蓋流形法在任意網格劃分上的優(yōu)勢開展梁板殼自適應分析以及自動計算研究，并體現(xiàn)其高階級數(shù)收斂快、自由度少的優(yōu)勢。

本文研究的另一個重要意義在于：對于有限元法而言，構造C1連續(xù)的近似函數(shù)一直存在比較大的困難，而獨立覆蓋流形法具有導數(shù)自然收斂的特性，無需特殊操作就可實現(xiàn)C1連續(xù)，這為4階微分方程的數(shù)值求解打下了很好的基礎。

致謝：感謝石根華先生的指導！