求解加筋板的位移與應力的一種解析計算方法

王小明

(中國艦船研究設計中心 武漢 430064)

0 引 言

加筋板是船舶、鋼結構橋梁和海洋平臺的主要構件.在橋梁領域，由于其縱向和橫向加筋導致這兩個方向的等效剛度不相等，通常將縱橫加筋板稱為正交異性板.對于這類正交異性板的計算理論有三個體系：①主梁體系，由面板、縱筋和橫筋組成主向梁，抵抗整體變形；②梁板體系，由面板、縱筋和橫筋組成一個整體，面板作為縱筋和橫筋的共同翼板；③面板體系，面板視為支撐在縱筋和橫筋上的各向同性連續板[1-3].對于第一體系，確定了翼板寬度，其力學分析問題可以由初等梁理論解決；對于第三體系，可以采用薄板理論分析；對于第二體系，其存在梁板耦合的系統，計算比較復雜，正是本文研究的內容.

關于第二體系變形和應力的計算方法，Younseok等[4]提出了將加筋板等效成正交異性板的詳細計算方法.這樣的算法存在一定的局限性，要求加強筋比較小，且密集布置，載荷均布或者加載在加強筋位置上，不滿足這兩點，計算精度就比較差，并且無法計算加強筋上的應力和位移.Jeom等[5]提出了基于Huber理論的P-E法，將縱橫加筋板的彎曲計算分為兩個步驟：①假設橫向加強筋剛度無限大，按剛性支持的正交異性板計算；②計算橫向加強筋的影響，然后修正第一步的計算結果.Bhaskar等[6]采用經典薄板理論和歐拉-伯努利梁理論研究了各向同性加筋板和正交異性加筋板的彎曲和自由振動，該方法的計算結果與三維有限元結果吻合良好.采用一階剪切變形理論和能量變分法推導加筋板的自由振動方程是常用的方法，關于振動微分方程的解法，Naveed等[7]并沒有完全采用廣為熟悉的雙三角級數的方法，而是引入了鐵木辛科(Timoshenko)梁函數作為形函數；Hamedani等[8]則是把微分方程轉化為有限元格式求解.Shanmugam等[9]采用有限元計算和試驗驗證的方式研究了加筋板在面內載荷和橫向載荷作用下的變形和破壞，研究結果認為面內載荷增加，橫向承載力將下降.彭林欣等[10-11]基于移動最小二乘近似，給出了加肋板線性彎曲的無網格法.運用無網格解得出的結果與加肋板三點彎曲實驗的位移結果吻合良好，與ANSYS有限元求解結果誤差不超過5%.Sapountzakis等[12-14]提出了一種板筋隔離法計算加筋板彎曲問題，將板筋拆解后，獨立研究，認為板筋系統的內力連同真實外力看作板的外力，板筋內力看作筋的外力，應用板筋結合面上的位移連續條件，求解橫向位移.研究中解微分方程的方法轉化為有限元格式計算，算例中重點關注的是橫向位移和板筋內力分布，對應力計算結果沒有研究.Lin等[15-17]先后用經典薄板理論、歐拉-伯努利梁理論、一階剪切變形理論和鐵木辛科梁理論推導加筋板的振動微分方程，以梁板間的耦合力和力矩為中間未知量，采用位移協調條件求解了自由邊界條件下、固支、簡支邊界條件下板的速度導納.在求解中，沒有考慮到筋板之間的水平剪切力，認為板筋的彎曲中面都是自身的幾何中面，在加強筋尺寸比較小，分布稀疏時計算結果差別不大，若加強筋尺寸大，且分布密集，則誤差較大.其重點討論的是一根筋的情況，實際工程中，往往是多根筋的，因此板筋之間的耦合力和力矩在不同加強筋之間互相影響；研究中關于加筋板的位移應力并沒有討論.雖然有限元計算已經成為當前的主流算法，但是采用通用有限元軟件計算，計算結果與網格劃分、計算者的水平和經驗密切相關.工程中的設計過程，經常是方案不斷變更，方案變更導致有限元模型修改工作量巨大.所以，在船舶設計領域，大部分的強度和剛度計算依然是以經典理論公式計算為主.在這樣的背景下，研究加筋板的解析計算方法依然非常必要.

為了求解加筋板的彎曲問題，文中提出了一種針對多根加強筋構成的加筋板的“板筋分離法”，與文獻[12-14]不同的是所有求解過程均采用半解析算法.這種方法以板筋之間的作用力和力矩為中間未知變量，將板筋之間的力、力矩和真實外力看作加筋板中板的外力，將板筋之間的力和力矩看作筋的外力，分別獨立求解二者的彎曲位移，通過位移協調條件確定板筋之間的耦合力和力矩.最終確定加筋板的位移和應力.這種方法實質上考慮了板筋之間的剪切力和法向力，比采用正交異性板法精度高；避免了能量法中需要針對不同截面分別計算變形能的麻煩，因此比能量法具有更廣的應用范圍，可以應用于工字型鋼加筋板，倒T型鋼加筋板等.

1 基本微分方程

建立坐標系，見圖1.加筋板沿x軸方向長度為a，沿y軸方向寬度為b，板厚為t.加強筋為角鋼，厚度為ts，高度為h.加強筋平行于x軸方向布置.外力大小為F0，作用點的坐標為(x0，y0).加筋板橫向彎曲時，自由度可以用φx，φy和w表示.φx為xz平面內板的轉角，φy為yz平面內板的轉角，都以板中面轉向z軸的方向為正.w為板的垂向(z向)位移，u0p，v0p分別為板中面的水平(x向和y向)位移.板和筋是焊接在一起的，因此分界面上具有相同的變形.

圖1 加筋板與坐標系統

受文獻[12-13]的啟發，并考慮板筋之間剪切力引起的剛體位移，對于多根加強筋的情形，板的彎曲微分方程可表示為

(1)

(2)

(3)

(4)

(5)

式中：N為加強筋的總數；D=Et3/12(1-μ2)為板的抗彎剛度；C=κEt/2(1+μ)為抗剪剛度；κ=5/6為剪切修正系數；E，μ分別為板的彈性模量和泊松比；δ為狄拉克函數.

四邊簡支的加筋板，其邊界條件可以表示為

x=0或a，則

(6)

y=0或b，則

(7)

對于加強筋而言，將其看作鐵木辛科梁，其自由度為xz平面的彎曲轉角θ(x)，yz平面的扭轉轉角φs(x)和垂向(z向)位移u(x).u0s為梁中面的水平(x向)位移.

以加強筋形心點為原點建立坐標系，即C為角鋼橫截面的形心位置，見圖2.加筋板在橫向彎曲時，第i根加強筋受到板對它的力qi、qxi和轉矩Mi的作用.

圖2 加筋板中板梁的耦合力與力矩

圖2中的Mi，僅示意其大小分布，方向由板(筋)上的雙箭頭表示.所以第i根加強筋的控制微分方程為

(8)

(9)

(10)

(11)

簡支梁的邊界條件為

x=0或a，則

(12)

2 基本微分方程的解

由式(1)～(5)組成的方程組與邊界條件式(6)～(7)構成定解問題.可以采用雙傅里葉級數來求解這一定解問題.

(16)

(17)

式中：km=mπ/a,kn=nπ/b.式(13)～(17)已經滿足了邊界條件式(6)～(7)，則僅需要通過方程組確定傅里葉級數的系數.方程(1)～(3)與方程(4)～(5)是互相獨立的，可以分別求解.將式(13)～(15)代入方程式(1)～(3).

(18)

(19)

由式(18)可以分別求出傅里葉系數.

(20)

(21)

(22)

方程(4)和(5)也可以采用相同的方法計算.

(23)

(24)

(25)

(26)

對加強筋所在的梁也可以作類似的求解與推導.以第i根加強筋為例，方程式(8)～(11)與邊界條件式(12)組成的定解問題，解的形式可以寫成

(27)

(28)

(29)

(30)

式(27)～(30)代入式(8)～(11)，可以求得，

(31)

由式(31)求出傅里葉系數.

(33)

(34)

(35)

(36)

根據位移協調條件，位于加筋板板與加強筋分界面位置的質點，他們的位移應該是相同的.則

ui(x)=w(x,yi),φis(x)=φy(x,yi)

(37)

u0p(x,yi)-φx(x,yi)t/2=u0si(x)+θi(x)h/2

(38)

將式(20)～(22)代入式(13)～(17)，式(33)～(36)代入式(27)～(30)，可以分別用板筋間的耦合力表示廣義位移，聯合協調條件式(37)～式(38)，可以求出耦合力.用矩陣表示為

(39)

(40)

(41)

(42)

式(41)～(42)中的β1=kny1，βN=knyN，η=1-μ.

通過式(39)～(42)就可以求解出板與筋之間的耦合力和力矩，再代入式(20)～(22)，則可以確定加筋板的位移和轉角.確定位移之后，求板中應力或者筋的應力分布，可以根據幾何方程和本構關系求得，即下面的式(43)～(46).

板的應力可以表示為

(43)

(44)

(45)

第i根加強筋的應力為

(46)

式中：ω為扇形坐標.因為加強筋除了彎曲產生的彎曲正應力，還有加強筋的約束扭轉產生的扇形正應力.

3 算例與討論

簡支加筋板長a=2 m、寬b=1.5 m.加強筋為扁鋼，高度h=0.1 m、厚度ts=0.01 m，等間距分布，間距ls=0.5 m，分別布置在y1=0.5 m，y2=1 m的位置.集中力F0=2×104N，作用在板中心的位置，即x0=a/2，y0=b/2.求解加筋板的位移分布和應力分布，求解加強筋中的正應力分布.

為了驗證本文算法的正確與否，下面采用4種算法，分別命名為方法1～4.方法1采用文獻[10]中的基本原理推導控制微分方程，同時考慮加強筋的扭轉變形能，但是求解微分方程的方法不是采用文獻中的最小二乘無網格法，而是采用伽遼金算法.方法2采用ANSYS有限元算法，板采用shell181，共4 800個單元；加強筋采用beam188模擬，共160個單元.方法3 本文算法.方法4 本文方法的退化算法，并忽略整體彎曲中面的偏移構件中面的影響.板采用薄板模型，加強筋采用歐拉-伯努利梁模型，則式(1)～(5)退化成一個關于橫向位移w(x,y)的算式(47).

(47)

加強筋的彎曲方程退化為

(48)

邊界條件和方程的求解也可以仿照前面的方法來調整.

相應板的應力計算式，則簡化為

(49)

(50)

加強筋的應力計算式修改為

(51)

選取幾個典型部位的位移分布和應力分布比較上述4種計算方法的精準度.圖3為y=b/4，b/2和x=a/4，a/2位置的位移分布；圖4為y=b/4位置處板上表面的x向應力分布(以下簡稱x向應力分布)；圖5為x=a/4位置處板上表面的y向應力分布(以下簡稱y向應力分布).圖6為加強筋下邊緣直角點位置處的正應力分布(以下簡稱正應力分布)，因為外力和加強筋都是關于y=b/2對稱，所以兩根筋的正應力分布完全相同.

圖3 不同位置的位移分布

圖4 y=b/4位置處的x向應力分布

圖5 x=a/4位置處的y向應力分布

圖6 加強筋中的正應力分布

由圖6可見：前3種方法關于位移的計算結果都比較接近，僅在中心的位置有細微差別.板中應力分布，方法2和方法3數據點比較一致，差別較小，方法1和方法4有的位置與方法2吻合很好，有的位置差別明顯.

有限元法(方法2)是目前普遍認同的準確算法，以方法2的計算值為標準值，則本文方法(方法3)關于位移的最大計算誤差為3.5%，板的應力最大計算誤差為-13.5%，加強筋中應力最大計算誤差為-0.1%.以上僅是以圖3～6少數幾個特殊位置計算結果得到的誤差分析結論，但是并非僅這些點的位置吻合較好，其他位置吻合就吻合不好.通過測試比較，其他位置的誤差也在上述范圍之內.需要說明的是，不能試圖比較載荷作用位置(算例中的板中心點)附近的應力計算差別，因為在載荷作用位置附近，質點的應力是無限大，各種方法都不收斂，計算差別很大.這一點，與薄板彎曲的結論是類似的.

4種計算方法中，方法4數值處理上最為簡單，精度上在y=b/4位置處的x向應力、加強筋中正應力與文中方法上差別較大.因為這種方法忽略了彎曲中面偏離加強筋中面的影響，即忽略了板筋之間的剪切作用；還忽略了加強筋的約束扭轉.加強筋的高度往往比板厚大得多，板筋之間的剪切作用比較明顯，彎曲中面應該在靠近板的一側，忽略這些，將導致加強筋附近板的應力和加強筋中應力計算誤差較大.

在數值處理上，本文方法不需要采用諸如伽遼金法或者李茲法求解微分方程，僅需要求解矩陣的逆矩陣或者解線性代數方程組，編程相對簡單，計算速度明顯比伽遼金算法快.上述4種方法，除方法2是由ANSYS軟件默認的迭代次數和計算精度外，其他3中方法中的三角級數的迭代累加次數都不可能是無限次數的，都是通過編程控制的.雙傅立葉級數收斂較慢，編程中所有廣義位移表達式級數的累加次數m=n≥11之后則收斂，方法3中Dm矩陣各個元素是單傅立葉級數，收斂很快，n≥7則收斂.上述計算結果和誤差分析，都是在級數收斂這個前提下所完成的.當然，累加次數越多，精度越高，但是計算機的耗時和內存消耗量也是巨大的.

關于本文方法的適用范圍，需要說明兩點.

1)公式推導中，以集中力作為代表，引入狄拉克函數來表達的，對于多個集中力的情況，僅需同時列出多個狄拉克函數累加，可以相應推導計算公式；均布載荷也是普遍存在的載荷形式，則需要把推導過程中的Fmn作如下替換即可.

2)上述插圖和算例中是以扁鋼作為代表的，對于其他形式的加強筋，如工字鋼，方鋼，倒T型鋼等也同樣適用，僅需要注意計算加強筋扭轉慣性矩應該采用下面的公式.

這一點是用能量法所不能達到的，因為計算應變能的時候，不同的加強筋截面需要采用不同的計算方法，如工字鋼、倒T型鋼和方鋼就要分別采用不同的公式計算應變能.

4 結 論

在推導加筋板中的板平衡微分方程，采用了一階剪切變形理論，加強筋(梁)采用了鐵木辛科梁理論，提出了一種“板筋分離法”計算加筋板的橫向彎曲.板與筋之間的耦合力和力矩作為中間變量，加強筋受到板對其的作用力和作用力矩，板也受到加強筋的反作用，以這兩對中間變量代入各自的平衡方程中，求出各自的位移.板與筋共同區域交界面上，二者具有相同的橫向位移和轉角，這就是位移協調條件.通過位移協調條件求解出中間變量耦合力和力矩，最后求出位移表達式.

1)通過算例驗證了本文計算方法的正確性，與有限元法相比，位移的計算結果最大誤差為3.5%；板應力計算結果與有限元法大部分區域吻合較好，最大誤差-13.5%；加強筋中正應力的最大誤差為-0.1%.造成這些局部誤差的主要原因是加強筋的約束扭轉沒有完全考慮.

2)文中的計算方法適用性上，可以用于計算集中載荷或均布載荷作用下的加筋板，加強筋的形式可以是扁鋼，工字鋼，方鋼，倒T型鋼等.