使用B樣條無單元法進行梯形蓋板受力分析

李 彬

(中鐵第四勘察設計院集團有限公司,武漢 430063)

鋼筋混凝土梯形蓋板大量應用于既有鐵路斜交涵洞中，板厚與跨度比一般在0.1～0.25，屬中厚板，在既有線改造或增建二線工程中，需要一種快速有效的方法對此類結構進行受力分析與評估。有限元法[9]是分析中厚板的有效方法，它通用性強，精度高，但是由于有限元法一般采用分片低階多項式插值來逼近問題的解函數，在邊界上一般只能達到C0連續，從而不可避免地造成了自由度多、工作量大的局面。許多學者將半解析法用于中厚板分析[10-11]，取得了較好的效果。本文介紹了利用B樣條無單元方法進行涵洞出入口鋼筋混凝土梯形蓋板受力分析的研究及其Matlab[12]實現。

1 基本方程

圖1 梯形蓋板示意

如圖1所示，考慮一個在x-y平面內的等厚梯形蓋板i-j-k-1,其中a為梯形蓋板短邊長、b為蓋板長邊長，L為計算跨度，t為厚度。i-j、k-1為自由邊，i-1、j-k為簡支邊。

板的位移包括了3個分量：第一個是沿z方向的位移w，第二個是沿y軸方向的轉角θx，第3個是沿x軸正方的轉角θy。采用中厚板理論進行分析[9]，面外剪應變及面內應變可分別表述為

γ=γxz

γyz=?w?x

?w?y-θx

θy≡w-θ(1)

ε=εx

εy

γxy=-z??x0

0??y

??y??x·θx

θy≡-zLθ(2)

設板厚為t，楊氏模量為E，剪切模量為G，泊松比為ν，則其應變能可表示為

int=1(Lθ)TD·((Lθ)·dΩ+

其中

D=Et312(1-ν2)α=56Gt

采用變分原理[13]，式(3)的變分形式為

2 樣條無單元法剛度列式的推導

將整個區域映射為在自然坐標(ξ,η)上的單位正方形，在此區域內對豎向位移采用雙四次B樣條進行擬合，對轉角采用雙三次B樣條進行擬合，近似公式如下

w=∑q,rBq,4(ξ)·Br,4(η)·a(q,r)=

[B4(ξ)?B4(η)]·{a}(5)

θx=∑q,rBq,3(ξ)·Br,3(η)·bx(q,r)=

[B3(ξ)?B3(η)]·{bx}(6)

θy=∑q,rBq,3(ξ)·Br,3(η)·by(q,r)=

[B3(ξ)?B3(η)]·{by}(7)

其中，Bk,d為d次B樣條的基函數；

Bd(ξ)=[B0,d(ξ),B1,d(ξ),B2,d(ξ),…,Bn,d(ξ)]；

Bd(ξ)?Bd(η)為Bd(ξ)與Bd(η)的張量積；

{a}、{bx}、{by}分別為對w、θx、θy進行擬合的系數向量

{a}=[…,a(q,r),…]T

{bx}=[…,bx(q,r),…]T

{by}=[…,by(q,r),…]T

下文中，記

{C}為進行受力分析時的基本未知量，相當于有限元法中的位移。

函數f關于(x,y)和(ξ,η)的導數的變換關系如下

?f?ξ

?f?η=?x?ξ?y?ξ

?x?η?y?η?f?x

?f?y=[J]?f?x

?f?y(8)

?f?x

?f?η(9)

按式(2)有

Lθ=??x0

0??y

??y??x·θx

θy=?θx?x

?θy?y

?θx?y+?θy?x(10)

將式(9)代入上式

?θx?η

?θy?ξ

?θy?η(11)

使用式(6)、(7)的θx、θy插值公式，則上式成為

??ξ[B3(ξ)?B3(η)]·{bx}

??η[B3(ξ)?B3(η)]·{bx}

??ξ[B3(ξ)?B3(η)]·{by}

??η[B3(ξ)?B3(η)]·{by}=

bx

by=

[Bb]·{C}(12)

按式(1)、(9)有

w-θ=?w?x

?w?y-θx

?w?η-θx

θy(13)

使用式(5)、(6)、(7)的w、θx、θy近似公式，則上式成為

w-θ

??η[B4(ξ)?B4(η)]·{a} -

[B3(ξ)?B3(η)]·{bx}

[B3(ξ)?B3(η)]·{by} =

B3(ξ)?B3(η)0

0B3(ξ)?B3(η) ×

bx

bx

by-

0B3(ξ)?B3(η)0

00B3(ξ)?B3(η)·a

bx

by=

0B3(ξ)?B3(η)0

00B3(ξ)?B3(η)·{C}=

[Bg]·{C}(14)

利用式(12)、(14)，式(4)可用離散形式表達為

(C)=δ(C)T·[K]·(C)(15)

其中[K]為剛度矩陣

3 等效荷載列式

對板上所施加的均布荷載q，其等效荷載列式按以下方式計算，列出其荷載勢能的變分式

δB4(ξ)?B4(η)00·a

bx

by·dΩ=

等效荷載列式可表達為

4 邊界條件

以在ξ=0邊上施加軟簡支(SS1)邊界條件為例，采用罰函數在整個能量泛函中加入一個修正項

penalty=12ρΓ(18)

其中為Γ簡支邊，對上式進行變分

δpenalty=ρ·δwξ=0Twξ=0dΓ=

B4ξξ=0?B4(η)·{a}dΓ=

B4ξξ=0?B4(η)00dΓ·(C) =

δ(C)T·Kpenalty·(C)(18)

求解時最后的剛度矩陣應為[K]+[Kpenalty]

5 算法的Matlab實現及驗證

本方法中B樣條節點在ξ、η方向上均勻分布，在ξ方向上有Kx+1列，在η方向上有Ky+1列，為減少節點數量，在首尾處的節點重復d次(d為B樣條次數)，即為所謂的開放均勻B樣條，其系數向量的長度為(Kx+d-1)×(Ky+d-1)，在豎向位移、轉角分別采用雙四次B樣條和雙三次B樣條進行擬合的情況下，系統總的自由度數即的長度為

(Kx+3)×(Ky+3)+2×(Kx+2)×(Ky+2)

各節點坐標

(ξi,ηj)=(i/Kx,j/Ky)…i∈[0:Kx]j∈[0:Ky]

在樣條節點所構成的網格背景上采用3點Gauss公式進行數值積分，以下將Gauss點稱作樣條的配點，利用Matlab樣條工具箱中的spcol函數計算配點處樣條的各階導數值，計算張量積則使用kron函數。

由于樣條無單元法具有算法簡潔、后處理方便的特點，借助于Matlab樣條工具箱的強大威力，開發的斜交梯形板分析程序僅用了不多300行的代碼，開發效率高，后期維護方便。

為驗證該算法的可靠性，計算了承受均布荷載的四邊簡支方板在不同厚度情況下的中心撓度與彎矩，并與文獻[14]提出的RPAQ單元及解析解作了比較，結果見表1，表中q為均布荷載，a為方板邊長，wc為方板中心點處的撓度，Mc為方板中心點處的彎矩，D為彎曲剛度。

表1 不同厚度四邊簡支方板中心撓度與中心彎矩

計算時取μ=0.3，采用12×12分劃的樣條進行插值，RPAQ單元采用16×16的單元劃分方式。從計算結果可以看出，本文算法不僅有著非常好的精度，而且沒有剪切閉鎖現象。

6 結論

本文提出了用于梯形蓋板分析的B樣條無單元法，詳列了相關剛度矩陣推導的過程，說明了生成荷載列式與施加邊界條件的方法。采用與文獻中類似方法相比，對位移采用比轉角高一次的B樣條進行插值，不僅提高了精度，而且避免了板變薄時出現的剪切閉鎖現象。

借助于Matlab及其Spline Toolbox編制了梯形蓋板專用分析及評估程序，擺脫了此類結構分析完全依賴于大型通用有限元軟件的局面，具有效率高、后處理方便的優點。通過對區域映射公式進行簡單調整，該方法就可以推廣到任意四邊形的中厚板及薄板分析，具有一定的通用性。

[1] H.Antes. Bicubic fundamental splines in plate bending[J]. Int. J. Numer. Meth. Engng， 1974,8(3):503-511.

[2] Atluri S N， Zhu T L. New Concepts in Meshless Methods[J]. Int. J.Numer. Meth. Engng， 2000,47(1-3):537-556.

[3] Ahmed A. Shabana， Ashraf M. Hamed. Use of B-Spline in the Finite Element Analysis: Comparison With ANCF Geometry[J]. J.Comput.Nonlinear Dynam. Volume 7， Issue 1， 011008， 2012.

[4] 石鐘慈.樣條有限元[J].計算數學,1979(1)：50-72.

[5] 何廣乾，周潤珍，林春哲.樣條函數法在解板殼問題中的應用[J].建筑結構學報，1981，2(2).

[6] 朱明杈.解任意四邊形板彎曲問題的樣條有限元法[J].計算教學，1987(1)：23-42.

[7] Fan， S. C. and Cheung， Y. K.， Analysis of shallow shells by spline finite strip method[J]， Engineering Structures， 1983,5(4):255-263.

[8] 劉效堯.樣條單元族及其在橋梁結構中的應用[J].土木工程學報,1983(2)：10-20.

[9] O. C. Zienkiewicz， R. L. Taylor. The Finite Element Method for Solid and Structural Mechanics[M]. Sixth Edition (Volume 2) Butterworth-Heinemann， 2005.

[10] 趙艷林.中厚板彎曲分析的樣條函數方法[J].廣西大學學報：自然科學版，1993(3):52-57.

[11] 沈小璞，何沛祥.樣條有限元法解中厚度矩形板彎曲[J].安徽建筑，1996(S1).

[12] Eva P?rt-Enander. Anders Sj?berg.MATLAB5手冊[M].北京：機械工業出版社，2000.

[13] 胡海昌.彈性力學的變分原理及其應用[M].北京：科學出版杜，1981.

[14] 龍馭球，陳曉明，岑松.一個不閉鎖和抗畸變的四邊形厚板元[J].計算力學學報,2005(8):385-391.