精確幾何薄曲梁曲殼分析的分區(qū)級數(shù)解

蘇海東，韓陸超，頡志強(qiáng)

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

1 研究背景

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

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

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

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

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

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

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

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

(1)

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

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

(2)

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

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

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

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

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

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

V=φ1V1+φ2V2。

(3)

不難驗(yàn)證C0連續(xù)性。條形區(qū)域關(guān)于x方向的一階導(dǎo)數(shù)為

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

(5)

3 薄曲梁的分區(qū)級數(shù)表達(dá)及應(yīng)變矩陣計(jì)算

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

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

根據(jù)梁保持平截面的假設(shè)，局部坐標(biāo)系下的兩個(gè)方向位移為[2]

(6)

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

(7)

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

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

(8)

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

(9)

分別記cosαxi、cosβxi為cosα、cosβ，同時(shí)考慮兩個(gè)局部坐標(biāo)軸的垂直關(guān)系，將局部坐標(biāo)下的位移轉(zhuǎn)換到整體坐標(biāo)下的位移，即

(10)

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

二維應(yīng)變子矩陣(對于第q項(xiàng))為

(12)

其中，由鏈?zhǔn)椒▌t，各項(xiàng)對整體坐標(biāo)的偏導(dǎo)數(shù)為

(13)

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

(14)

然后D數(shù)組的每一項(xiàng)再與C數(shù)組的每一項(xiàng)按式(14)的方式進(jìn)行組合。最后依次放入式(12)的對應(yīng)位置，形成應(yīng)變子矩陣Bq，并代入剛度矩陣公式計(jì)算單元?jiǎng)偠染仃嘯6]。

對于連接覆蓋1和覆蓋2的條形區(qū)域，求關(guān)于整體坐標(biāo)x或y的偏導(dǎo)數(shù)時(shí)，Bq矩陣每一項(xiàng)系數(shù)c還要與權(quán)系數(shù)φi(i=1,2)進(jìn)行組合，形成應(yīng)變子矩陣Biq，則有

(15)

其他關(guān)于積分方法等內(nèi)容參見文獻(xiàn)[6]。另外，取R→∞和h=1可退化到計(jì)算細(xì)長直梁。

4 曲梁的幾何計(jì)算公式

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

(16)

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

(2)關(guān)于t的1—3階導(dǎo)數(shù)為：

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

(18)

(4)曲率半徑R為

(19)

(5)方向余弦為:

(20)

(6)局部坐標(biāo)關(guān)于整體坐標(biāo)的導(dǎo)數(shù)。曲梁上的任一點(diǎn)為

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

(22)

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

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

(23)

5 薄曲殼計(jì)算簡述

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

(24)

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

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

(26)

其中，對于正交的局部坐標(biāo)系，兩個(gè)方向的轉(zhuǎn)角為[4]

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

(28)

其中：

(29)

將局部坐標(biāo)系下的位移轉(zhuǎn)化到整體坐標(biāo)系下的位移，并記其中的函數(shù)公式(對于第q項(xiàng))為

(30)

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

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

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

6 算 例

采用薄梁板殼假設(shè)，將文獻(xiàn)[6]和文獻(xiàn)[7]的曲梁和曲殼主要算例重新計(jì)算一遍。

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

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

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

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

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

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

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

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

6.3 球面殼算例

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

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

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

表3 球殼最大位移計(jì)算結(jié)果比較Table 3 Comparison of the maximum displacements ofthe spherical shell

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

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

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

7 結(jié) 語

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

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

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

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

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

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

致謝：感謝石根華先生的指導(dǎo)！