基于精確幾何的曲梁分析新方法

, 　,

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

1　研究背景

梁、板、殼是工程結(jié)構(gòu)的常見形式，其中，曲梁和曲殼由于其優(yōu)美的流線造型和良好的受力特性，在工程中有著廣泛的應(yīng)用。

目前，以有限元法為代表的數(shù)值計(jì)算方法是梁板殼力學(xué)分析的主要方法[1]。梁板殼分析一般基于2個(gè)基本假設(shè)：①對(duì)于原先垂直于中面的橫截面，在變形過(guò)程中始終保持為平面；②沿厚度方向產(chǎn)生的直接應(yīng)變可以忽略。對(duì)于細(xì)長(zhǎng)梁和薄板殼，可進(jìn)一步忽略剪切變形，在基本假設(shè)①的基礎(chǔ)上考慮“變形后的截面仍垂直于中面”的“直法線”假設(shè)。具體分析中，大都以中面位移為未知量，考慮位移沿厚度方向的線性變化，由此推導(dǎo)出不同于實(shí)體分析的控制方程，一般要求的C1連續(xù)性給構(gòu)造近似函數(shù)帶來(lái)一些困難；另一種方式是位移和轉(zhuǎn)角獨(dú)立插值，或直接從實(shí)體單元退化得到梁板殼單元，僅要求C0連續(xù)，并能考慮厚度較大時(shí)的剪切變形，但在求解細(xì)長(zhǎng)梁或薄板殼時(shí)易出現(xiàn)剪切自鎖。

對(duì)于曲梁和曲殼而言，由于存在初始曲率，其受力分析比直梁和平板要復(fù)雜得多。參考文獻(xiàn)[1]—文獻(xiàn)[4]，曲梁和曲殼的有限元分析主要有2種途徑：

(1)一種簡(jiǎn)單方式就是用直梁或平板單元近似地模擬曲梁或曲殼。隨著單元數(shù)量的增加，在幾何上和變形上都會(huì)逼近真實(shí)的曲梁或曲殼。這種方式固然簡(jiǎn)單，但為保證幾何相似性會(huì)使用較多單元，而且難以很好地反映曲梁或曲殼實(shí)際存在的一系列耦合問(wèn)題，造成求解過(guò)程中收斂慢和計(jì)算誤差大。

(2)另一途徑是直接構(gòu)造曲梁和曲殼單元。但由于曲率的存在，需要對(duì)任意形狀(常曲率、變曲率)下的應(yīng)變-位移關(guān)系進(jìn)行精確描述，還要考慮軸力、彎矩、剪力、扭矩等作用的耦合，甚至要涉及三維空間的非歐幾何，幾何方程、物理方程和平衡微分方程都比較復(fù)雜，還存在不同假設(shè)或簡(jiǎn)化條件下的殼體理論，因此，曲梁和曲殼單元的研究難度遠(yuǎn)大于直梁和平板單元。

近年來(lái)興起的等幾何分析方法[5]，幾何模型和分析模型采用統(tǒng)一的NURBS(Non-Uniform Rational B-Splines，非均勻有理B樣條)描述方式，實(shí)現(xiàn)了幾何模型的保形性，計(jì)算精度高，適合于對(duì)幾何精確性要求高的曲梁和曲殼分析。但等幾何分析方法也存在基函數(shù)相對(duì)復(fù)雜、本質(zhì)邊界條件不易施加、提高階次也不能完全避免薄殼計(jì)算中的剪切自鎖等問(wèn)題[6]。

考慮到曲梁和曲殼有著相似的變形特征，曲梁可看作曲殼在空間上的退化形式，因此，本文僅討論二維曲梁，但本質(zhì)上是著眼于曲梁和曲殼分析的一般性方法。

筆者在文獻(xiàn)[7]中基于“獨(dú)立覆蓋流形法”，提出梁的獨(dú)立覆蓋分析方法，采用實(shí)體計(jì)算退化方式，使某些多項(xiàng)式基函數(shù)不參與計(jì)算，就能模擬梁的基本假設(shè)，而且僅要求近似函數(shù)的C0連續(xù)并從根本上解決了剪切自鎖等問(wèn)題。在此基礎(chǔ)上，本文將這種方法推廣到二維曲梁，直接在精確的中面曲線參數(shù)方程基礎(chǔ)上計(jì)算曲梁的變形，為三維曲梁和下一步的曲殼分析探索新的路徑。為簡(jiǎn)便起見，本文以單位寬度的矩形截面梁為例，且暫時(shí)只考慮中面用一條曲線描述的情況。

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

圖1　任意形狀的獨(dú)立覆蓋及其條形連接Fig.1　Independent coverswith arbitrary shapeconnected through strips

圖2　曲梁的獨(dú)立覆蓋劃分及其條形連接Fig.2A curved beam divided by independent coversconnected through strips

2012年，作者在石根華博士的建議下，在流形法中首次引入“獨(dú)立覆蓋”[10-11]，即單位分解函數(shù)φi=1、近似函數(shù)V就是給定級(jí)數(shù)Vi的覆蓋區(qū)域，獨(dú)立覆蓋之間為窄條形的覆蓋重疊區(qū)域，其單位分解函數(shù)取為有限元線性形函數(shù)以實(shí)現(xiàn)覆蓋函數(shù)的線性過(guò)渡。這樣，每個(gè)覆蓋就包含1個(gè)獨(dú)立覆蓋及其條形重疊區(qū)域。首先研究了矩形獨(dú)立覆蓋，2013年提出任意形狀的獨(dú)立覆蓋[12]，如圖1所示，3個(gè)任意形狀的獨(dú)立覆蓋用條形(包括條形交叉處的過(guò)渡單元)連接，這樣，覆蓋的劃分不僅能適應(yīng)求解域的物理邊界，還能嚴(yán)格施加本質(zhì)邊界條件。比如，將圖2所示的細(xì)長(zhǎng)曲梁劃分為5個(gè)獨(dú)立覆蓋，其間用條形相連，條形上的單位分解函數(shù)就是一維有限單元的形函數(shù)，同時(shí)，在梁的端部通過(guò)條形與定義了邊界條件的結(jié)點(diǎn)相連。

2015年筆者在文獻(xiàn)[13]中正式命名“基于獨(dú)立覆蓋的數(shù)值流形方法”，簡(jiǎn)稱“獨(dú)立覆蓋流形法”，其中，獨(dú)立覆蓋和條形都是基本的計(jì)算單元(或稱為覆蓋網(wǎng)格)，但強(qiáng)調(diào)獨(dú)立覆蓋的主體作用。

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

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

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

與文獻(xiàn)[7]類似，將梁作為二維實(shí)體考慮，設(shè)局部坐標(biāo)xi和yi方向的位移分別為ui和vi，參考文獻(xiàn)[1]，對(duì)應(yīng)于文獻(xiàn)[7]中的Timoshenko直梁，設(shè)曲梁各點(diǎn)在局部坐標(biāo)下的位移為

(1)

梁中面的軸向和厚度方向位移分別設(shè)為f1(xi)和f3(xi)，假設(shè)xi=s，表明梁中面位移僅僅隨弧長(zhǎng)s變化，與yi沒(méi)有關(guān)系。根據(jù)基本假設(shè)②——沿厚度yi方向的應(yīng)變可忽略，因而對(duì)于固定了xi坐標(biāo)的梁上任一點(diǎn)，有vi(xi)=f3(xi)，再根據(jù)基本假設(shè)①，用yif2(xi)描述軸向纖維沿厚度方向呈線性變化的伸長(zhǎng)量(相對(duì)于中面位移f1(xi))，其中f2(xi)代表獨(dú)立定義的隨xi=s變化的截面轉(zhuǎn)角θ。

考慮關(guān)于局部坐標(biāo)的完全多項(xiàng)式，則有

(2)

(a) 軸線方向

(b) 厚度方向

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

圖5　曲梁的獨(dú)立覆蓋及其連接Fig.5　Independent covers and their connections oncurved beams

4　引入局部坐標(biāo)的曲梁分析公式

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

(3)

在獨(dú)立覆蓋1和獨(dú)立覆蓋2的條形重疊區(qū)域，比如圖5的A點(diǎn)和B點(diǎn)之間，有

(4)

設(shè)

Tink=wiLiNink。

(5)

(6)

(7)

上式：

(8)

其中：

(9)

與文獻(xiàn)[7]相比，考慮到αi隨坐標(biāo)的變化，因此式(7)多了局部坐標(biāo)系的方向余弦關(guān)于整體坐標(biāo)的導(dǎo)數(shù)。

如圖3所示，由簡(jiǎn)單的幾何關(guān)系，得到各點(diǎn)的整體直角坐標(biāo)與局部坐標(biāo)的轉(zhuǎn)換公式為

(10)

式(10)再結(jié)合中面(x0i,y0i)的曲線表達(dá)式，就可以將式(9)、式(8)乃至式(7)計(jì)算出來(lái)。

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

(11)

式中：Ω為計(jì)算單元的區(qū)域；i=1，2；j=1，2；q，t依次對(duì)第i個(gè)和第j個(gè)多項(xiàng)式覆蓋函數(shù)中參與計(jì)算的所有項(xiàng)進(jìn)行循環(huán)；D為平面應(yīng)力的彈性矩陣，其表達(dá)式為

(12)

式中:E為彈性模量;μ為泊松比。

集中荷載的子向量為

(13)

式中Fx和Fy分別為水平向和豎向的集中荷載。分布荷載也有類似的采用積分形式的公式。

將上述單元?jiǎng)偠染仃嚭蛦卧奢d向量分別組集成整體剛度矩陣和整體荷載向量，然后加上邊界條件，就可以求解線性方程組得到多項(xiàng)式的系數(shù)，過(guò)程與文獻(xiàn)[7]類似。

5　積分方式及具體計(jì)算過(guò)程

?ΩF(x(φ,r),y(φ,r))dxdy=?ΩF(φ,r)dsdr=

?ΩF(φ,r)s′(φ)dφdr。

(14)

式中，ds=s′(φ)dφ。

令

(15)

式中：r1(φ)和r2(φ)分別表示沿厚度方向積分的下限和上限(即-h/2到h/2，h為厚度，可以考慮h隨φ的變化)。采用一維Gauss數(shù)值積分[2]計(jì)算,即

(16)

積分點(diǎn)的個(gè)數(shù)nl與覆蓋函數(shù)關(guān)于xi=φ的階次有關(guān)，Hl為積分權(quán)值。對(duì)于固定的積分點(diǎn)φl(shuí)，式(15)的F1(φl(shuí))一般也采用一維Gauss積分，由于式(1)中的覆蓋函數(shù)只出現(xiàn)關(guān)于yi=r的常數(shù)項(xiàng)和一次式，因此，只用2點(diǎn)Gauss積分就可以滿足積分精度要求。更進(jìn)一步有可能在式(15)中用解析積分的方式消掉r，推導(dǎo)出簡(jiǎn)化的F1(φ)表達(dá)式，這樣就只需一次一維Gauss積分。

綜上所述，只需在積分點(diǎn)(包括沿厚度方向和沿軸線方向)上計(jì)算被積函數(shù)值，再乘以相應(yīng)的積分權(quán)值，就可以得到式(11)的剛度子矩陣。在具體計(jì)算過(guò)程中，關(guān)鍵是式(9)中的局部坐標(biāo)和式(7)中的方向余弦關(guān)于整體坐標(biāo)的導(dǎo)數(shù)計(jì)算，以下給出計(jì)算步驟：

(1)列出中面曲線的參數(shù)方程,即

(17)

(2)根據(jù)式(17)，求出積分點(diǎn)處的φ。

(5)求α的正弦和余弦關(guān)于x或y的導(dǎo)數(shù),即:

(18)

圖6　一段圓形曲梁及其變形(一端固定，另一端受水平力作用)Fig.6　A circular curved beam and its deformation(with one end constrained and another subjected tohorizontal load)

6　圓形曲梁算例

考慮中面曲線的曲率不變，如圖6所示，計(jì)算第1象限的1/4圓形曲梁。按上節(jié)的5個(gè)步驟計(jì)算如下：

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

式中:(x0,y0)為圓心坐標(biāo);R為半徑，得到弧長(zhǎng)的微分ds=Rdφ。

設(shè)圓的半徑為1 m，梁的截面厚度為0.01 m。只用1個(gè)獨(dú)立覆蓋，左端通過(guò)一個(gè)窄條與固定端相連，右端承受水平力F=1 kN。梁的彈性模量E=108kN/m2，泊松比μ=0。

軸向采用6點(diǎn)Gauss積分，計(jì)算得到荷載作用點(diǎn)處的位移為：覆蓋函數(shù)取5階多項(xiàng)式(總共17個(gè)自由度)時(shí)，u=-0.093 9 m，v=-0.059 8 m；取6階和7階多項(xiàng)式(分別為20個(gè)和23個(gè)自由度)時(shí)，u=-0.094 1 m，v=-0.059 8 m。

將文獻(xiàn)[14]的理論解用關(guān)于φ和r的多項(xiàng)式展開(考慮R=1，s=Rφ=φ)，得

(19)

可見，局部坐標(biāo)xi=φ，yi=r的位移uφ和ur在中面處隨參數(shù)方程的參數(shù)φ而變化，uφ沿厚度方向呈線性變化(關(guān)于r2項(xiàng)的系數(shù)很小)，ur沿厚度方向的變化可忽略，這與梁的基本假設(shè)一致。在底部φ=0處，理論解的位移與本文計(jì)算結(jié)果非常接近，不僅如此，計(jì)算得出的獨(dú)立覆蓋上關(guān)于φ和r多項(xiàng)式的各項(xiàng)系數(shù)，也與式(19)比較接近。

如圖7(a)所示，只用1個(gè)獨(dú)立覆蓋，左端和右端分別通過(guò)窄條與固定端相連，在梁的中點(diǎn)作用合力為F=1 kN的法向集中力。梁的彈性模量E=105kN/m2，泊松比μ=0。

軸向采用10點(diǎn)Gauss積分，覆蓋函數(shù)取10階和11階多項(xiàng)式(分別為32個(gè)和35個(gè)自由度)，計(jì)算得到荷載作用點(diǎn)處的位移u=v=-0.107 m，與細(xì)密網(wǎng)格的有限元解u=v=-0.112 m的相對(duì)誤差<5%。但進(jìn)一步升階時(shí)計(jì)算不收斂，有可能是積分精度不夠。重新劃分為5個(gè)獨(dú)立覆蓋，如圖7(b)所示，各獨(dú)立覆蓋只需采用4階多項(xiàng)式(總自由度為70個(gè))，軸向采用6點(diǎn)Gauss積分，得到u=v=-0.110 m。因此，實(shí)際計(jì)算中不一定要用很高的多項(xiàng)式階次，在控制好計(jì)算誤差的前提下，可采用升階和加密覆蓋的2種方式來(lái)提高精度。

圖7　圓形曲梁及其變形(兩端固定，中點(diǎn)受集中力)Fig.7　A circular curved beam and its deformation(with two ends constrained and the middle pointsubjected to concentrated load)

值得關(guān)注的是，當(dāng)采用1個(gè)獨(dú)立覆蓋時(shí)，獨(dú)立覆蓋的單元弧長(zhǎng)與截面厚度之比超過(guò)150∶1，但未出現(xiàn)有限元法的Timoshenko梁?jiǎn)卧谇蠼饧?xì)長(zhǎng)梁時(shí)的剪切自鎖現(xiàn)象。

圖8　橢圓形曲梁及其變形(一端固定，另一端受集中力作用)Fig.8　An ellipse curved beam and its deformation(with one end constrained and another subjected toconcentrated load)

7　橢圓形曲梁算例

考慮中面曲線的曲率變化，如圖8所示，計(jì)算第1象限的1/4橢圓形曲梁。具體計(jì)算過(guò)程如下：

(5)按式(18)求α的正弦和余弦關(guān)于x或y的導(dǎo)數(shù)。

設(shè)橢圓的半長(zhǎng)軸a=2 m，半短軸b=1 m，梁的截面厚度為0.01 m。只用1個(gè)獨(dú)立覆蓋，左端通過(guò)一個(gè)窄條與固定端相連，右端分別承受水平力和豎向力F=1 kN。梁的彈性模量E=108kN/m2，泊松比μ=0。

采用劃分100個(gè)單元的有限元結(jié)果作為對(duì)比，本文用10階多項(xiàng)式(自由度為32)的計(jì)算結(jié)果如下：水平力作用下，有限元為u=-0.168 5 m，v=-0.223 5 m，本文為u=-0.168 3 m，v=-0.223 0 m，如圖7(a)所示；豎向力作用下，有限元為u=-0.223 5 m，v=-0.326 0 m，本文為u=-0.223 0 m，v=-0.325 2 m，如圖7(b)所示?？梢钥闯?，各截面上的軸力、彎矩、剪力等各種內(nèi)力荷載組合作用下的變形都能夠反映出來(lái)。

類似地，關(guān)于r的表達(dá)式為低階有理多項(xiàng)式，可用解析積分方式，但公式推導(dǎo)相對(duì)于圓形曲梁更為復(fù)雜。

8　結(jié)　語(yǔ)

本文采用獨(dú)立覆蓋流形法，在前期研究的關(guān)于直梁的獨(dú)立覆蓋分析方法基礎(chǔ)上，只需借助隨中面參數(shù)方程變化的局部坐標(biāo)系，并考慮該坐標(biāo)系的局部坐標(biāo)以及方向余弦關(guān)于整體坐標(biāo)的導(dǎo)數(shù)，就能實(shí)現(xiàn)精確幾何描述下的曲梁分析。文中給出常曲率的一段圓形曲梁和變曲率的一段橢圓形曲梁2個(gè)算例，驗(yàn)證了方法的可行性。

本文方法的特點(diǎn)是：完全采用實(shí)體分析模式，只需使多項(xiàng)式覆蓋函數(shù)中的某些項(xiàng)不參與計(jì)算，就能模擬梁的基本假設(shè)，相比較而言，通常的曲梁?jiǎn)卧獦?gòu)造比較復(fù)雜，需要推導(dǎo)曲梁控制方程和相應(yīng)的離散公式，考慮曲率對(duì)幾何方程、平衡方程的影響以及各種內(nèi)力耦合；同時(shí)，本文方法還保持了梁的獨(dú)立覆蓋分析方法的一些基本優(yōu)點(diǎn)，如僅要求近似函數(shù)的C0連續(xù)性、無(wú)剪切自鎖、不存在由于厚度遠(yuǎn)小于其他尺度而導(dǎo)致的數(shù)值病態(tài)、與實(shí)體連接方便等；在實(shí)現(xiàn)幾何保形性方面，這是除了等幾何分析方法之外的另一途徑，與之相比，本文方法具有多項(xiàng)式基函數(shù)相對(duì)簡(jiǎn)單、嚴(yán)格施加本質(zhì)邊界條件、完全消除剪切自鎖等優(yōu)點(diǎn)。

正如引言所述，本文不限于討論二維曲梁，而是希望找到解決曲梁和曲殼類問(wèn)題的新途徑。后續(xù)研究將緊緊圍繞這一主題，探討新方法對(duì)各種曲梁和曲殼的適應(yīng)性：參考文獻(xiàn)[7]，考慮到截面轉(zhuǎn)角θ=dv/dxi，模擬細(xì)長(zhǎng)曲梁的“直法線假設(shè)”；考慮中面為不同曲線參數(shù)方程的連接，并給出多種厚度以及變截面算例；最近完成了母線為曲線描述的旋轉(zhuǎn)殼分析，將另文介紹；三維曲梁和曲殼分析的思路與本文一致，因此將按空間曲線或曲面參數(shù)方程建立隨中面變化的曲線坐標(biāo)系，并計(jì)算該坐標(biāo)系的曲線坐標(biāo)和方向余弦關(guān)于整體坐標(biāo)的導(dǎo)數(shù)，最終在三維空間內(nèi)實(shí)現(xiàn)精確幾何描述下的曲梁和曲殼分析。

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

參考文獻(xiàn)：

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

[2]王勖成．有限單元法[M]．北京：清華大學(xué)出版社，2003.

[3]潘科琪，劉錦陽(yáng). 柔性曲梁多體系統(tǒng)的研究現(xiàn)狀和展望[J]. 力學(xué)進(jìn)展，2011，41(6)：711-721.

[4]趙躍宇，康厚軍，馮銳，等. 曲線梁研究進(jìn)展[J]. 力學(xué)進(jìn)展，2006，36(2)：170-186.

[5]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.

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

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

[8]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.

[9]BABUSKA I, MELENK J M. The Partition of Unity Method[J]. International Journal for Numerical Methods in Engineering, 1997, 40:727-758.

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

[11] 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.

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

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

[14] 王敏中，王煒，武際可．彈性力學(xué)教程(修訂版)[M]．北京：北京大學(xué)出版社，2011.