采用獨(dú)立覆蓋流形法分析精確幾何描述的曲殼

,　,,　

(長江科學(xué)院 材料與結(jié)構(gòu)研究所，武漢　430010)

1　研究背景

曲殼(也稱為殼)不僅造型優(yōu)美，而且是一種非常經(jīng)濟(jì)的承載結(jié)構(gòu)，在工程中得到廣泛應(yīng)用。曲殼的力學(xué)分析，除了極少數(shù)的簡單情況有解析解答外，絕大多數(shù)都采用以有限元法[1]為主的數(shù)值計(jì)算方式。

理論上，曲殼可以采用由控制方程建立積分弱形式的標(biāo)準(zhǔn)方式來建立數(shù)值分析表達(dá)式。然而，曲殼控制方程(包括平衡微分方程、物理方程和幾何方程)的推導(dǎo)本身就比較復(fù)雜，涉及曲面微分幾何的很多內(nèi)容，而且基于不同假設(shè)會(huì)得到不同方程，因此鮮有直接應(yīng)用曲殼控制方程形成數(shù)值計(jì)算過程的嘗試[2]，僅有適用于特殊簡化理論的扁殼單元等極少類型。

在有限元分析中，應(yīng)用最多的是由平板單元拼接成曲殼的方式[1]。理論上，足夠多的平板可以逼近連續(xù)曲面的精確幾何形狀，并模擬其真實(shí)的力學(xué)行為。這種幾何精確性對(duì)曲殼非常重要，然而目前一般依賴人為判斷，容易產(chǎn)生幾何誤差進(jìn)而帶來力學(xué)分析上的誤差。即使對(duì)于平板單元本身而言，由薄板彎曲的4階微分方程得到的積分弱形式，要求近似函數(shù)具有C1連續(xù)性，這一直以來都是研究難點(diǎn)和熱點(diǎn)，直至現(xiàn)在還不斷有新的研究成果，由此產(chǎn)生了各種協(xié)調(diào)的或非協(xié)調(diào)的甚至是基于多場分析的平板單元，但大都有各自的適用范圍，給用戶的選擇造成困難。

采用實(shí)體計(jì)算的退化方式是另一種有效途徑[2]。它應(yīng)用殼體的Reissner-Mindlin假設(shè)，可以考慮橫向剪切變形，而且近似函數(shù)只需要C0連續(xù)性，更重要的是，可以在有限單元的各邊增加結(jié)點(diǎn)以模擬曲面幾何。但其缺陷是容易出現(xiàn)剪切自鎖和薄膜自鎖，雖然可以通過減縮積分等方式進(jìn)行處理，但可能帶來矩陣的奇異性等其他問題。

近年來出現(xiàn)的等幾何分析方法[3]，強(qiáng)調(diào)幾何的精確性，適用于對(duì)幾何要求很嚴(yán)的曲殼分析，但其基函數(shù)較為復(fù)雜，對(duì)復(fù)雜求解域的適應(yīng)性不強(qiáng)，而且難以完全避免自鎖等問題[4]。

我們?cè)谇捌谘芯恐谢讵?dú)立覆蓋流形法，采用二維實(shí)體分析方式，提出了直梁和曲梁的獨(dú)立覆蓋分析方法[5-6]。其中，只需借助隨中面參數(shù)方程變化的局部坐標(biāo)系，并考慮該坐標(biāo)系的局部坐標(biāo)以及方向余弦關(guān)于整體坐標(biāo)的導(dǎo)數(shù)，就能實(shí)現(xiàn)精確幾何描述下的曲梁分析。文獻(xiàn)[6]給出常曲率的圓形曲梁和變曲率的橢圓形曲梁算例，驗(yàn)證了該方法的可行性。本文將此方法推廣到三維板殼(包括三維曲殼和三維平板)分析，從而形成一整套梁板殼分析的新方法。

2　獨(dú)立覆蓋流形法簡介

圖1　獨(dú)立覆蓋及其條形連接Fig.1　Independent covers connected through strips

2012年，筆者在石根華博士的建議下，在流形法中首次引入“獨(dú)立覆蓋”[9-10]，即單位分解函數(shù)φi=1、近似函數(shù)V就是給定級(jí)數(shù)Vi的覆蓋區(qū)域，獨(dú)立覆蓋之間為窄條形的覆蓋重疊區(qū)域，其單位分解函數(shù)取為有限元的線性形函數(shù)以實(shí)現(xiàn)覆蓋函數(shù)的線性過渡。每個(gè)覆蓋包含一個(gè)獨(dú)立覆蓋及其周邊的條形重疊區(qū)域。如圖1(a)所示，首先研究了矩形獨(dú)立覆蓋系統(tǒng)，圖中的大矩形為獨(dú)立覆蓋；2013年提出任意形狀的獨(dú)立覆蓋[11]，如圖1(b)所示，A，B，C，D，E這5個(gè)任意形狀的獨(dú)立覆蓋用條形連接，為方便計(jì)算，不規(guī)則的條形以及條形交叉處都劃分為三角形，這樣，覆蓋的劃分不僅能適應(yīng)求解域的物理邊界，還能嚴(yán)格施加本質(zhì)邊界條件[12]。

2015年筆者在文獻(xiàn)[13]中正式命名“基于獨(dú)立覆蓋的數(shù)值流形方法”，簡稱“獨(dú)立覆蓋流形法”，其中，獨(dú)立覆蓋和條形都是基本的計(jì)算單元(稱為覆蓋網(wǎng)格)，但強(qiáng)調(diào)獨(dú)立覆蓋的主體作用。獨(dú)立覆蓋流形法的收斂性是在各覆蓋區(qū)域內(nèi)采用的完備級(jí)數(shù)(如多項(xiàng)式級(jí)數(shù))逼近真實(shí)解的基礎(chǔ)上建立的，不僅包括物理場本身的收斂，而且包括高階導(dǎo)數(shù)的收斂。覆蓋網(wǎng)格具有任意形狀、任意連接、任意加密的特性。

3　獨(dú)立覆蓋上的曲殼覆蓋函數(shù)

覆蓋函數(shù)的多樣性是獨(dú)立覆蓋流形法的顯著優(yōu)勢?？梢酝ㄟ^選擇覆蓋函數(shù)以反映物理場的局部特性或整體特性，前者如文獻(xiàn)[10]的裂紋分析，在裂紋尖端附近選擇解析級(jí)數(shù)作為覆蓋函數(shù)從而加快收斂，后者在文獻(xiàn)[5]、文獻(xiàn)[6]以及本文所討論的梁板殼分析中，選擇多項(xiàng)式覆蓋函數(shù)的合適形式以反映梁板殼的整體變形特征。

如圖2所示的整體直角坐標(biāo)系x-y-z下，在第i個(gè)獨(dú)立覆蓋中，中面坐標(biāo)(x0i,y0i,z0i)由殼體曲面參數(shù)方程描述，以(x0i,y0i,z0i)為原點(diǎn)建立沿中面變化的局部正交坐標(biāo)系xi-yi-zi，其中xi和yi沿中面內(nèi)的曲線坐標(biāo)方向(切向)，zi沿厚度方向(即中面的法向)。xi和yi考慮為曲線坐標(biāo)，zi表示殼體上的點(diǎn)到中面的距離。若中面由多個(gè)曲面連接而成，相同曲面的覆蓋可以合并使用同一曲面定義的曲線坐標(biāo)。

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

值得注意的是，一般的多項(xiàng)式級(jí)數(shù)表達(dá)都是關(guān)于整體坐標(biāo)系的，本文將其推廣到隨中面變化的局部坐標(biāo)系下的多項(xiàng)式級(jí)數(shù)表達(dá)。將殼體作為三維實(shí)體考慮，設(shè)局部坐標(biāo)xi，yi，zi方向的位移分別為ui,vi,wi，參考文獻(xiàn)[1]，根據(jù)殼體的Reissner-Mindlin假設(shè)，設(shè)殼體上的各點(diǎn)在局部坐標(biāo)下的位移為

(1)

式中：殼體中面位移僅僅隨曲線坐標(biāo)xi和yi變化，與zi沒有關(guān)系，因此3個(gè)方向的中面位移分別設(shè)為f1(xi,yi),f2(xi,yi),f3(xi,yi)；沿厚度zi方向的應(yīng)變可忽略，因而對(duì)于固定了xi和yi坐標(biāo)的殼體上的任一點(diǎn)，有wi(xi,yi)=f3(xi,yi)；再根據(jù)平截面假設(shè)，用zif4(xi,yi)和zif5(xi,yi)分別描述xi向和yi向纖維沿厚度方向呈線性變化的伸長量(分別相對(duì)于中面位移f1(xi,yi)和f2(xi,yi))，其中f4(xi,yi)和f5(xi,yi)分別表示獨(dú)立定義的隨中面坐標(biāo)變化的截面轉(zhuǎn)角θxi和θyi。

關(guān)于局部坐標(biāo)的完全多項(xiàng)式為

(2)

如圖3所示的多項(xiàng)式排序，顯示了從0階到2階的完全多項(xiàng)式。應(yīng)用式(2)構(gòu)造的局部近似函數(shù)，采用實(shí)體計(jì)算模式，只需考慮虛線左側(cè)的各項(xiàng)，而使虛線右側(cè)的各項(xiàng)不參與計(jì)算，就能得到式(1)，實(shí)現(xiàn)殼體的Reissner-Mindlin假設(shè)。另外，當(dāng)3個(gè)方向的位移取同階多項(xiàng)式時(shí)，wi=f3(xi,yi)比θxi=f4(xi,yi)

(a)xi和yi方向

(b)zi方向

圖3多項(xiàng)式排列
Fig.3Polynomialorder

和θyi=f5(xi,yi)高1階，從而避免了剪切自鎖問題。在xi和yi方向上，定義了關(guān)于xi和yi的完全多項(xiàng)式以充分反映薄膜變形，從式(1)可以看出，純彎模式在中面(zi=0)處不會(huì)產(chǎn)生附加的薄膜變形，因此避免了薄膜自鎖問題。

4　曲殼分析公式

如圖2所示，在第i個(gè)獨(dú)立覆蓋內(nèi)，引入轉(zhuǎn)換矩陣Li(各項(xiàng)為局部坐標(biāo)軸的方向余弦)，將曲殼各點(diǎn)的局部位移轉(zhuǎn)換到整體坐標(biāo)系下的位移，即

(3)

式中：cosαxi，cosβxi，cosγxi分別是xi軸關(guān)于整體坐標(biāo)x，y，z的方向余弦，其他依此類推。

如果網(wǎng)格在xi-yi上投影為矩形單元，可以采用矩形有限單元的形函數(shù)作為單位分解函數(shù)，并應(yīng)用強(qiáng)制約束關(guān)系[9-10]實(shí)現(xiàn)獨(dú)立覆蓋和條形。比如，在圖1(a)中將結(jié)點(diǎn)2、結(jié)點(diǎn)7、結(jié)點(diǎn)8的所有自由度約束到結(jié)點(diǎn)1，實(shí)現(xiàn)獨(dú)立覆蓋1；將結(jié)點(diǎn)4、結(jié)點(diǎn)9、結(jié)點(diǎn)10約束到結(jié)點(diǎn)3，實(shí)現(xiàn)獨(dú)立覆蓋3；兩者之間的條形為線性連接。

對(duì)于任意形狀覆蓋，根據(jù)文獻(xiàn)[11]，在圖1(b)的三角形單元1-6-7中，將結(jié)點(diǎn)6、結(jié)點(diǎn)7約束到結(jié)點(diǎn)1，實(shí)現(xiàn)獨(dú)立覆蓋A；獨(dú)立覆蓋B類似；A和B之間的條形由2個(gè)三角形單元組成，由結(jié)點(diǎn)7約束到結(jié)點(diǎn)1，結(jié)點(diǎn)8約束到結(jié)點(diǎn)2，實(shí)現(xiàn)A和B之間的線性過渡。因此所有的計(jì)算單元(包括獨(dú)立覆蓋及條形)都可以采用同樣的計(jì)算公式(暫不考慮強(qiáng)制約束)，不失一般性，以下公式都在有限單元內(nèi)采用高階的常規(guī)流形法公式進(jìn)行推導(dǎo)，位移表示為

(4)

式中：對(duì)于三角形單元，m=3；對(duì)于矩形單元，m=4；Wi對(duì)應(yīng)于第i個(gè)覆蓋(結(jié)點(diǎn))的單位分解函數(shù)，取為二維有限單元的形函數(shù)，可由單元內(nèi)各點(diǎn)在xi-yi上的投影坐標(biāo)確定[1]。設(shè)

(5)

(6)

(7)

其中,

(8)

fiq,y和fiq,z依此類推。

其中：

式(7)第3個(gè)等號(hào)后面的第2項(xiàng)矩陣為局部坐標(biāo)系的方向余弦關(guān)于整體坐標(biāo)的導(dǎo)數(shù)計(jì)算，對(duì)于平板分析，由于方向余弦不變，因此無需考慮。

如圖2所示，考慮zi軸上的殼體點(diǎn)(x,y,z)與中面點(diǎn)(x0i,y0i,z0i)構(gòu)成的矢量，由簡單的幾何關(guān)系，得到殼體點(diǎn)的整體直角坐標(biāo)與局部坐標(biāo)的轉(zhuǎn)換公式為

(10)

以式(10)再結(jié)合中面點(diǎn)(x0i,y0i,z0i)的曲面參數(shù)方程表達(dá)式，就可以將式(9)、式(8)乃至式(7)具體地計(jì)算出來，見下節(jié)中關(guān)于幾何的具體計(jì)算步驟。

按最小勢能原理推導(dǎo)單元?jiǎng)偠染仃嚨淖泳仃嚍?/p>

(11)

式中：Ω為計(jì)算單元的三維空間區(qū)域；i=1,2,…,m；j=1,2,…,m；q，t分別表示對(duì)第i和第j覆蓋所對(duì)應(yīng)的多項(xiàng)式覆蓋函數(shù)中參與計(jì)算的所有項(xiàng)進(jìn)行循環(huán)；D為空間彈性矩陣[1]。

集中荷載和分布荷載的子向量分別為：

(12)

(13)

式中：Fx,Fy,Fz分別為整體坐標(biāo)3個(gè)方向的集中荷載；qx，qy，qz為分布荷載集度；A為分布荷載作用面積。

將上述單元?jiǎng)偠染仃嚭蛦卧奢d向量分別組集成整體剛度矩陣和整體荷載向量，然后參考文獻(xiàn)[12]加上邊界條件，就可以求解線性方程組得到多項(xiàng)式的系數(shù)，再根據(jù)式(4)和式(7)得到整體坐標(biāo)系下的位移和應(yīng)變。

5　積分方式及幾何計(jì)算公式

本文僅考慮中面為一個(gè)曲面的情況，只定義一個(gè)隨中面變化的正交局部坐標(biāo)系。設(shè)中面的曲面參數(shù)方程為

(14)

將曲線坐標(biāo)設(shè)為xi=ξ，yi=η，厚度方向zi=r，且在中面處r=0，則局部坐標(biāo)表示為ξ-η-r。另外，αxi，βxi，γxi等分別用αx，βx，γx等表示。

?ΩF[x(ξ,η,r),y(ξ,η,r),z(ξ,η,r)]dxdydz=

?ΩF(ξ,η,r)h1h2h3dξdηdr。

(15)

式中：h1，h2，h3分別為3個(gè)局部坐標(biāo)的拉梅系數(shù)，其中，在zi坐標(biāo)軸上用r表示點(diǎn)到中面的距離，則h3=1。

令

(16)

式中r1(ξ,η)和r2(ξ,η)分別表示沿厚度方向積分的下限和上限(即-h/2到h/2，h為厚度，可以考慮h隨ξ-η的變化)。將單元各邊與單元形心相連形成三角形，采用Hammer數(shù)值積分[1]計(jì)算，即

(17)

綜上所述，只需在積分點(diǎn)(包括在中面投影區(qū)域內(nèi)和沿厚度方向)上計(jì)算被積函數(shù)值，再乘以相應(yīng)的積分權(quán)值，就可以得到式(11)的剛度子矩陣。式(13)關(guān)于分布荷載的計(jì)算類似。

從新媒體在2016年美國總統(tǒng)大選及在特朗普陣營發(fā)揮的作用看，特朗普現(xiàn)象不應(yīng)僅解讀為民粹主義本身的再度崛起，而是民粹主義結(jié)合新媒體取得成功的例證。民粹主義的出現(xiàn)既有歷史文化原因亦有受政治機(jī)會(huì)主義者煽動(dòng)利用因素。如今，民粹隨經(jīng)濟(jì)危機(jī)再度回歸，公眾不再跟隨精英的觀點(diǎn)成為2016年大選的背景（Gelman&Azari 2017），這一態(tài)勢為特朗普成功利用并進(jìn)一步塑造。民粹和民主的界限本身也頗為微妙，“服務(wù)人民”即便未曾出自每個(gè)美國總統(tǒng)競選者之口，也毫無疑問是不可否認(rèn)的命題。呼喚“多數(shù)人的統(tǒng)治”本身既是政治理念又是政治策略，這賦予民粹主義在美國政治生態(tài)中頑強(qiáng)的生命力。

在關(guān)于幾何的具體計(jì)算過程中，關(guān)鍵是局部坐標(biāo)和方向余弦關(guān)于整體坐標(biāo)的導(dǎo)數(shù)計(jì)算，以下給出計(jì)算步驟：

(1)根據(jù)中面的曲面參數(shù)方程，得到xi軸、yi軸、zi軸的基向量分別為(式(18)和式(19)暫用x，y，z分別代替xi，yi，zi：x=x0i(ξ,η);y=y0i(ξ,η);z=z0i(ξ,η)):

(18)

(2)由式(18)求出積分點(diǎn)(ξj,ηj)處xi軸、yi軸、zi軸關(guān)于整體坐標(biāo)的方向余弦為：

(19)

其中，

L=(y,ξz,η-z,ξy,η)2+(z,ξx,η-x,ξz,η)2+(x,ξy,η-y,ξx,η)2

(3)根據(jù)式(10)可得

(20)

對(duì)于確定的某種曲面，以上推導(dǎo)公式的過程是一次性的。

6　球面殼算例

如圖4(a)所示的球坐標(biāo)系，按第5節(jié)的4個(gè)步驟計(jì)算如下：

圖4　球坐標(biāo)及第1象限內(nèi)的球面殼有限元網(wǎng)格Fig.4　Spherical coordinates and a sphericalshell in the first quadrant

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

式中：(x0,y0,z0)為球心坐標(biāo)；R為半徑；ξ=θ；η=φ。因此有：

然后得到拉梅系數(shù)h1=R，h2=Rsinθ。

(2)αx，βx，γx，αy，βy，γy，αz，βz，γz的余弦值可表示為：

(3)由式(10)和球殼中面參數(shù)方程得

將上式對(duì)x求導(dǎo)，解方程求出：

同理分別對(duì)y和z求導(dǎo)，得到：

(4)將各方向余弦分別對(duì)x求導(dǎo)，得到

-ξ,xsinθcosφ-η,xcosθsinφ。

其他方向余弦的導(dǎo)數(shù)依此類推。

對(duì)y,z求導(dǎo)，只需將上式中的x改為y或者z。

計(jì)算如圖4(a)所示的第1象限內(nèi)的球面殼，半徑R=1 m，θ=π/4～π/2，φ=0～π/2。殼體厚度為0.01 m，左端(φ=0)和右端(φ=π/2)為固支約束，球面殼上表面承受法向均布?jí)毫=10 kPa。彈性模量E=106kPa，泊松比μ=0。圖4(b)顯示了用于計(jì)算球面殼參考解的有限元細(xì)密網(wǎng)格，由5 000個(gè)平板單元拼接而成。

如圖5所示，在φ-θ的投影平面上，流形法計(jì)算分別采用3種矩形有限元網(wǎng)格：網(wǎng)格1——只用1個(gè)獨(dú)立覆蓋；網(wǎng)格2——5×3(φ向×θ向)個(gè)獨(dú)立覆蓋；網(wǎng)格3——9×5個(gè)獨(dú)立覆蓋。采用文獻(xiàn)[12]的邊界條件施加方法：在網(wǎng)格1的左、右兩端各增加1個(gè)條形，與固定約束的結(jié)點(diǎn)相連；在網(wǎng)格2和網(wǎng)格3的左、右端各個(gè)獨(dú)立覆蓋中，直接設(shè)定覆蓋函數(shù)使之事先滿足邊界條件。

圖5　φ-θ投影平面上的覆蓋網(wǎng)格Fig.5　Cover meshes in φ-θ projection plane

計(jì)算得到的各方向最大位移見圖6，出現(xiàn)在(θ=π/2，φ=π/4)點(diǎn)處：當(dāng)采用網(wǎng)格1(1個(gè)獨(dú)立覆蓋)時(shí)，必須取到9階多項(xiàng)式覆蓋函數(shù)時(shí)才能收斂，x向(由對(duì)稱性，x向和y向相等)和z向最大位移分別為0.583和0.275(單位10-2m，下同)；網(wǎng)格2(5×3個(gè)獨(dú)立覆蓋)，5階多項(xiàng)式的結(jié)果為0.579和0.272，6階多項(xiàng)式為0.581和0.273；網(wǎng)格3(9×5個(gè)獨(dú)立覆蓋)，4階多項(xiàng)式的結(jié)果為0.580和0.272，5階多項(xiàng)式為0.581和0.273。而圖4(b)細(xì)密網(wǎng)格的有限元參考解為0.581和0.274，與3種網(wǎng)格下的收斂解幾乎一致?？梢?，不管采用粗網(wǎng)格配以較高階次的多項(xiàng)式，還是采用細(xì)網(wǎng)格配以較低階次的多項(xiàng)式，都有很好的收斂性。在控制好計(jì)算誤差的前提下，可采用升階和加密覆蓋的2種方式來提高精度。

圖6　球面殼的最大位移Fig.6　Maximum displacements of the spherical shell

關(guān)于ξ-η平面上的三角形Hammer積分點(diǎn)數(shù)，對(duì)網(wǎng)格1的試算表明：4階及以下的多項(xiàng)式覆蓋函數(shù)情況可用7點(diǎn)積分；5階和6階為12點(diǎn)積分；7階及以上至少需16點(diǎn)積分。因此，較高的多項(xiàng)式階次雖然可用較少的自由度達(dá)成收斂，但剛度矩陣積分的計(jì)算量也很大。

當(dāng)去掉式(7)等號(hào)右邊的第2項(xiàng)矩陣，就可以進(jìn)行平板計(jì)算。如圖5所示的平板，x=0～π/2(單位:m)，y=π/4～π/2(單位:m)，同樣的左、右兩端固支約束和均布荷載情況下，網(wǎng)格1(1個(gè)獨(dú)立覆蓋)采用4階多項(xiàng)式覆蓋函數(shù)計(jì)算得到中點(diǎn)處的最大位移為1.903(單位10-2m，下同)，與細(xì)密網(wǎng)格的有限元解相同。若再考慮上、下端面的固支約束，則6階多項(xiàng)式計(jì)算得到中點(diǎn)處的最大位移為0.115 6；網(wǎng)格3(9×5個(gè)獨(dú)立覆蓋)采用4階多項(xiàng)式得到的最大位移為0.115 8；而細(xì)密網(wǎng)格的有限元參考解同樣為0.115 8?？梢?，三維平板計(jì)算也有很好的收斂性。

值得注意的是，當(dāng)采用1個(gè)獨(dú)立覆蓋時(shí)，單元最大邊長與截面厚度之比超過150∶1，但未出現(xiàn)有限元法基于實(shí)體退化方式求解薄板殼時(shí)的自鎖現(xiàn)象。

從圖3的多項(xiàng)式選擇來看，目前采用的方式是在某一階次的完全多項(xiàng)式中去除關(guān)于zi的二次以上的高階項(xiàng)。事實(shí)上，在以彎曲為主的板殼分析中，可以減少圖3(a)中用于模擬薄膜變形的關(guān)于xi和yi的多項(xiàng)式階次，比如上面的平板算例屬于純彎變形，可以將關(guān)于xi和yi的多項(xiàng)式完全去掉，使自由度數(shù)大幅下降。進(jìn)一步可以參考文獻(xiàn)[5]，基于物理場導(dǎo)數(shù)的收斂性，考慮截面轉(zhuǎn)角θxi=?wi/?xi和θyi=?wi/?yi以模擬薄板的Kirchhoff-Love假設(shè)，從而再去掉圖3(a)中關(guān)于zi的自由度。

限于研究時(shí)間較短，本文僅給出了常曲率的球面殼算例，未給出變曲率的曲殼算例。但從前文的理論推導(dǎo)，以及文獻(xiàn)[6]給出的橢圓形曲梁算例來看，本文方法應(yīng)用于變曲率的曲殼分析也無問題。另外，對(duì)于形狀更為復(fù)雜的曲殼，可采用任意形狀的獨(dú)立覆蓋。

7　結(jié)　語

本文采用獨(dú)立覆蓋流形法，在前期研究的關(guān)于直梁和曲梁的獨(dú)立覆蓋分析方法基礎(chǔ)上，只需借助隨中面參數(shù)方程變化的局部坐標(biāo)系，并考慮該坐標(biāo)系的局部坐標(biāo)以及方向余弦關(guān)于整體坐標(biāo)的導(dǎo)數(shù)，就能實(shí)現(xiàn)精確幾何描述下的曲殼分析，而忽略坐標(biāo)系的變化即可分析平板。

文獻(xiàn)[5]的直梁、文獻(xiàn)[6]的曲梁，以及本文的曲殼和平板分析，初步構(gòu)成了一套基于獨(dú)立覆蓋流形法的梁板殼分析新方法。其特點(diǎn)是：

(1)完全采用實(shí)體分析模式，只需使多項(xiàng)式覆蓋函數(shù)中的某些項(xiàng)不參與計(jì)算，就能模擬梁板殼變形的基本假設(shè)，從而避免了推導(dǎo)梁板殼控制方程以及數(shù)值分析公式的過程，特別對(duì)于曲梁和曲殼而言，該過程的復(fù)雜性體現(xiàn)在應(yīng)用曲面微分幾何推導(dǎo)不同假設(shè)下的控制方程，以及由積分弱形式推導(dǎo)數(shù)值計(jì)算的離散公式。另外，與實(shí)體單元的連接也無需特殊處理。

(2)不管是模擬細(xì)長的Euler-Bernoulli梁或薄板的Kirchhoff-Love假設(shè)，還是考慮橫向剪切的Timoshenko梁或Reissner-Mindlin板殼假設(shè)，都僅要求近似函數(shù)的C0連續(xù)性，其中，前者需要考慮獨(dú)立覆蓋流形法的導(dǎo)數(shù)收斂性(薄板殼的相關(guān)公式及算例將在后續(xù)給出)。

(3)在數(shù)值計(jì)算的穩(wěn)定性方面，雖然是實(shí)體分析模式，但不存在有限元實(shí)體計(jì)算中由于厚度遠(yuǎn)小于其他尺度而導(dǎo)致的數(shù)值病態(tài)，也無任何自鎖現(xiàn)象。

(4)在實(shí)現(xiàn)幾何保形性方面，避免了通常由直梁或平板單元構(gòu)造曲梁或曲殼所造成的幾何誤差。這是除了等幾何分析方法之外的另一途徑，與之相比，本文方法具有多項(xiàng)式基函數(shù)相對(duì)簡單、嚴(yán)格施加本質(zhì)邊界條件、完全消除自鎖、覆蓋加密方便等優(yōu)點(diǎn)。

因此目前看來，這種新方法可以解決現(xiàn)有梁板殼分析存在的主要基礎(chǔ)性問題。究其原因，以有限元法的自鎖問題為例，有限元法基于插值方式構(gòu)造近似函數(shù)，當(dāng)采用實(shí)體退化方式分析梁板殼時(shí)，有限單元的插值函數(shù)并不是總能夠描述正確的彎曲變形[2]：剪切自鎖源于在純彎模式下產(chǎn)生了偽橫向剪切；薄膜自鎖也是在純彎模式下不能表達(dá)中面的無伸縮變形。而獨(dú)立覆蓋流形法是基于分區(qū)(覆蓋)的直接逼近方法，所選用的多項(xiàng)式函數(shù)能夠準(zhǔn)確地表達(dá)梁板殼的彎曲變形和薄膜變形。

當(dāng)然，目前還未就新方法的計(jì)算效率進(jìn)行分析比較。新方法可以采用粗網(wǎng)格配以高階覆蓋函數(shù)，或者細(xì)網(wǎng)格配以較低階的覆蓋函數(shù)，甚至是結(jié)合獨(dú)立覆蓋流形法的精度分析實(shí)現(xiàn)覆蓋任意加密的自適應(yīng)計(jì)算，需要研究哪種方式(或自適應(yīng)途徑)的計(jì)算量更小，并通過優(yōu)化算法進(jìn)一步減少計(jì)算量。最近，筆者已完成了母線為曲線描述的旋轉(zhuǎn)殼分析，正在考慮中面為不同曲線的連接；下一步要通過更多算例探討新方法對(duì)各種梁板殼(包括三維曲梁)分析的適應(yīng)性；將來進(jìn)行梁板殼大位移的幾何非線性分析，材料非線性分析，以及動(dòng)力分析等。

致謝:感謝石根華博士的指導(dǎo)！

參考文獻(xiàn)：

[1]ZIENKIEWICZ O C , TAYLOR R L. 有限元方法[M].5版. 莊茁，岑松,譯. 北京: 清華大學(xué)出版社, 2006.

[2]BELYTSCHKO T, LIU W K, MORAN B.連續(xù)體和結(jié)構(gòu)的非線性有限元 [M]．莊茁,譯．北京：清華大學(xué)出版社，2002.

[3]HUGHES T J R，COTTRELL J A，BAZILEVS Y. Isogeometric Analysis：CAD，F(xiàn)inite Elements，NURBS，Exact Geometry and Mesh Refinement[J]. Computer Methods in Applied Mechanics and Engineering，2005,194(39/40/41)：4135-4195.

[4]李新康. 層合結(jié)構(gòu)等幾何分析研究[D]. 杭州：浙江大學(xué)，2015.

[5]蘇海東，頡志強(qiáng). 梁的獨(dú)立覆蓋分析方法[J]. 長江科學(xué)院院報(bào), 2018,35(4):143-150.

[6]蘇海東，周朝，頡志強(qiáng). 基于精確幾何的曲梁分析新方法[J]. 長江科學(xué)院院報(bào), 2018,35(4):151-157,166.

[7]SHI G H. Manifold Method of Material Analysis[C]∥U.S. Army Research Office. Transactions of the Ninth Army Conference on Applied Mathematics and Computing, Minneapolis, Minnesota, U.S.A., June 18-21, 1991: 51-76.

[8]BABUSKA I,MELENK J M.The Partition of Unity Method[J]. International Journal for Numerical Methods in Engineering, 2015, 40(4):727-758.

[9]祁勇峰, 蘇海東, 崔建華．部分重疊覆蓋的數(shù)值流形方法初步研究[J].長江科學(xué)院院報(bào), 2013,30(1):65-70.

[10] SU Hai-dong, QI Yong-feng, GONG Ya-qi,etal. Preliminary Research of Numerical Manifold Method Based on Partly Overlapping Rectangular Covers[C]∥DDA Commission of International Society for Rock Mechanics. Proceedings of the 11th International Conference on Analysis of Discontinuous Deformation (ICADD11), Fukuoka, Japan, August 27-29, 2013, London: Taylor & Francis Group, 2013: 341-347.

[11] 蘇海東, 祁勇峰, 龔亞琦, 等. 任意形狀覆蓋的數(shù)值流形方法初步研究[J]. 長江科學(xué)院院報(bào), 2013, 30(12): 91-96.

[12] 蘇海東，頡志強(qiáng)．獨(dú)立覆蓋流形法的本質(zhì)邊界條件施加方法[J]. 長江科學(xué)院院報(bào), 2017,34(12)：140-146.

[13] 蘇海東，頡志強(qiáng)，龔亞琦, 等．基于獨(dú)立覆蓋的流形法收斂性及覆蓋網(wǎng)格特性[J]. 長江科學(xué)院院報(bào)，2016, 33(2)：131-136.