有限元 4-薄板彎曲問題

1、土木工程學院有限單元法有限單元法P-1/64有限單元法有限單元法第4章 彈性薄板彎曲問題的有限元法土木工程學院有限單元法有限單元法P-2/64薄板是一種常見的工程構件形式機械、航空和土建工程應用廣泛特殊形式小撓度薄板薄板彎曲問題在理論上和應用上都具有重要意義，并有專門著作加以論述（如楊耀乾平板理論）。象其它彈性力學問題一樣，用微分方程、差分法等經典方法所能求解的薄板問題很有限，一般只能解決等厚、小孔口、支承情況較簡單的單跨板。故工程設計中以往多采用簡化、近似、圖表等方法來解決板的設計問題。 土木工程學院有限單元法有限單元法P-3/64板殼結構板殼結構 板和殼是指厚度比其他尺寸要小得多的平面或曲

2、面構件。由于它的這種幾何特點，三維單元并不適合用來分析它們的變形。因為三維單元在三個方向的尺寸應盡量接近，否則求解精度由于“剪切自鎖”(shear locking)或系統(tǒng)矩陣病態(tài)而大大降低，甚至得到錯誤的結果。所以必須采用很細密的網格來適應板和殼的幾何特征，但是這將導致有限元模型的自由度瘋狂地增長。土木工程學院有限單元法有限單元法P-4/64平板上所承受的荷載分類平板上所承受的荷載通常有兩種:1. 面內拉壓荷載面內拉壓荷載。 由面內拉壓剛度承擔, 屬平面應力問題。2. 垂直于板的法向荷載垂直于板的法向荷載。 彎扭變形為主，具有梁的受力特征，即常說的彎曲問題。平板在垂直于板面的荷載作用下產生撓度

3、w。 板的面內變形與彎曲變形相互獨立，一般單獨推出各自的單元，然后再組裝在一起。平板殼體單元平面應力單元彎曲板單元土木工程學院有限單元法有限單元法P-5/64 當最大撓度w遠小于t 時, 稱為小撓度問題（or剛性板）當最大撓度w與t相差不大時,稱為大撓度問題（or柔性板）大撓度問題與小撓度問題51tw551tw5tw工程定義: 為剛性板；為柔性板；為絕對柔性板。土木工程學院有限單元法有限單元法P-6/64薄板與中厚板在板的分析中，常取板的中性面為xoy平面(如圖)。平板結構按其厚度t與短邊a的比值大小而分為:中厚板（Middle-Thick plate）和薄板（Thin plate)兩種。當

4、時稱為薄板薄板薄板寬度與厚度比值寬度與厚度比值在在15以上以上1at土木工程學院有限單元法有限單元法P-7/64中厚板薄板膜厚度與邊長之比較小很小非常小理論Mindlin 板理論Reissner 板理論Kirchhoff 板理論膜理論特點考慮橫向剪切變形忽略橫向剪切變形忽略彎曲變形與梁比較深梁細長梁弦土木工程學院有限單元法有限單元法P-8/64 當最大撓度w遠小于t 時, 稱為小撓度問題（or剛性板）當最大撓度w與t相差不大時,稱為大撓度問題（or柔性板）51tw551tw5tw工程定義: 為剛性板；為柔性板；為絕對柔性板。土木工程學院有限單元法有限單元法P-9/64一、基本假定一、基本假定、

5、略去垂直于中面的法向應力。（ ），并以中面上沿Z方向的撓度w代表板的撓度。、變形前垂直中面的任意直線，變形后仍保持為垂直中面的直線。（法向假定, ）、板彎曲時，中面不產生平面內應力。(中面中性層假定， ) 上述假定常稱為薄板小撓度問題假定（or 柯克霍夫Kirchhoof假定）。符合上述假定的平板即為剛性板剛性板。4.1 基本理論基本理論0z0zx0zy,0 xyyx和梁對比理解土木工程學院有限單元法有限單元法P-10/64二、基本方程二、基本方程 以上述假定為基礎，板分析中常用撓度w(x,y)作為基本未知量。、幾何方程（應變、幾何方程（應變撓度關系）撓度關系）彈性曲面沿x, y 方向的傾角從

6、中面取出一微小矩形ABCD，如圖所示，設其邊長為dx, dy，變形后彎曲成曲面ABCD土木工程學院有限單元法有限單元法P-11/64設A點撓度為w, 則沿x方向傾角(繞y軸)xwyB點撓度 dxxww沿y方向傾角(繞x軸) ywxD點撓度 dyyww土木工程學院有限單元法有限單元法P-12/64 沿x, y 方向位移作平行于xoz平面,設中面上點A到A1的距離為z，變形后，A點有撓度w, 同時發(fā)生彎曲，曲面沿x方向的傾角為 , 根據法線假定,則A1點沿x方向的位移: (負號為方向與x相反) 同理取 yoz平面得，zywvzxwuxw /土木工程學院有限單元法有限單元法P-13/64 z平面的應

7、變分量和曲、扭率基本假定，由于 , 故板內任意點的應變與平面問題相同: 0zzxzyxvyuyvxuxyyx yxwzywzxwzxyyx222222將u,v代入此為z平面的應變撓度幾何方程。 (4-1-2) 土木工程學院有限單元法有限單元法P-14/64上式中的 , , 為曲面在x,y方向的曲、扭率。這三個參數稱作彈性曲面的彎扭變形分量，它們完全確定了薄板內各點的應變分量。記為： 22xw22ywyxw2 z所以, chi yxwywxw222222凱，器，凱，器，西西土木工程學院有限單元法有限單元法P-15/64、物理方程（應力、物理方程（應力撓度關系）撓度關系）由于忽略z 對變形的影響,

8、 因此z平面的應力應變關系具有與平面問題相同的形式: xyxyxyyyxxEEE121122將(4-1-2)代入 022222222222111DyxwEzxwywEzywxwEzxyyx土木工程學院有限單元法有限單元法P-16/64或簡寫為： 0Dz式中彈性矩陣（平面應力問題）: 21000101120ED(4-1-4) 土木工程學院有限單元法有限單元法P-17/64板與梁的比較M土木工程學院有限單元法有限單元法P-18/64、內力方程（內力、內力方程（內力撓度關系）撓度關系）從板內取微元體tdxdy,由其上正應力x,y和剪應力xy ，可在截面上合成合力矩: Mx（yoz面上由x產生的繞y軸

9、彎矩） My（xoz面上由y產生的繞x軸彎矩）扭矩: Mxy（由剪應力產生） 土木工程學院有限單元法有限單元法P-19/64應力符號：正面正向為正xyzMxMxMyMyMxyMxyMyxMyxO 彎矩和扭矩Mx , My , Mxy , Myx 按右手螺旋法則用雙箭頭矢量來表示力偶，如圖所示（圖中所示各力偶的方向均為正）土木工程學院有限單元法有限單元法P-20/64t/2t/2zdzxyzxMxyMxyMyxMyOxyyx11假定Mx, My ,Mxy分別表示單位寬度上的內力矩。于是，內力矩陣： yxwywxwDt2222203212 dzDzdzzMMMMttttxyyx222202土木工程

10、學院有限單元法有限單元法P-21/64簡寫成 0312DtM (4-1-5) 前面推導了(4-1-4) Mzt312由此可見，平板上、下表面處的應力最大： Mttz226 0Dz(4-1-4) 比較(4-1-4)和(4-1-5)可得用內力矩表示的平板應力: 土木工程學院有限單元法有限單元法P-22/644、板的平衡方程),()2(4422444yxqywyxwxwD薄板彎曲的彈性曲面方程：)1 (1223EtD式中： 以上是薄板彎曲問題中的基本公式，從中可見其撓度w是彎曲問題中的基本未知函數。且由于忽略了z方向的變化，因此它只是x，y 的函數: w=w(x, y)。若w w已知，則位移、內力、

11、應力均可按上述相應公式求出。在經典解析法中，w(x, y)常設為三角級數形式。土木工程學院有限單元法有限單元法P-23/64在經典解析法中，w(x, y)常設為三角級數形式。例如，四邊簡支矩形板的w(x, y)設為: (納維爾解)11sinsin,mnmnbynaxmAyxw式中Amn為待定系數。假定荷載 11sinsin,mnmnbynaxmqyxq則可得位移函數: bynaxmbnamqDyxwmnsinsin1,222224土木工程學院有限單元法有限單元法P-24/644.2 有限元分析方法有限元分析方法 一、矩形單元的典型形式一、矩形單元的典型形式 將圖示矩形薄板沿x,y方向劃分成若干

12、小矩形(常取等分)。從中取出一小矩形(單元)，共有四個結點，此時不能象在平面問題中一樣，將結點此時不能象在平面問題中一樣，將結點視為視為“鉸鉸”，而是，而是“剛性的剛性的”，即每個結點有三，即每個結點有三個位移分量：個位移分量： 土木工程學院有限單元法有限單元法P-25/64即結點i的位移方向傾角上節(jié)為沿軸轉角繞方向傾角上節(jié)為沿軸轉角繞撓度xyyx wyx iyixiiixwywwwd 4 , 3 , 2 , 1i同理,相應的結點力 ）軸力偶（上節(jié)中的繞）軸力偶（上節(jié)中的繞豎向力xyMyMxyixiziiMMFF方向和正負號問題？右手螺旋法則土木工程學院有限單元法有限單元法P-26/64符號重

13、新定義是為了有限元表示的方便，由此得單元結位移向量 TyxyxTewwddd44411141節(jié)點力 TeFFF41TyxzyxzMMFMMF444111土木工程學院有限單元法有限單元法P-27/64二、 位移模式（函數）、位移模式的選取、位移模式的選取插值多項式取為： 29283726524321,xyyxxyxyxyxyxw312311310 xyyxy (4-2-1) 在上式中,前10項取到了三次項的全部,最后兩項則是從五個四次項 中選用了兩個。沒選 是因為它沒有多一項與其配對，沒選 ,它們在邊界上給出的撓度函數是四次的，比 和 要高一次,較之更難滿足邊界的協調條件。 432234yxyy

14、xyxx22yx44,yxyx33xy土木工程學院有限單元法有限單元法P-28/64、位移模式的檢驗、位移模式的檢驗（三個基本要求： 剛體位移，常應變，盡可能的邊界協調） 前三項含單元的剛體位移狀態(tài)：第一項1 與坐標x, y無關, 表示z方向的撓度是常量, 剛體移動表示剛體轉動第三項第二項32xxyyywxw 二次項代表均勻變形狀態(tài):曲率 4222xw6222yw扭率 522yxw土木工程學院有限單元法有限單元法P-29/64 能保證相鄰單元在公共邊界上撓度的連續(xù)性 以單元的邊界為例，在此邊界上y=-b，為常量，代入位移模式 (4-2-1)，可知邊界上的撓度w是x的三次函數，合并整理后可得：3

15、4232121xcxcxccw兩個端點共有4個邊界條件(結點1,2的撓度w1 , w2 ,和轉角y1,y2 ）。利用他們可唯一確定四個常數C1 C4。因為相鄰單元在結點1, 2的w, y對應相同，則兩個單元依據四個條件得到的C1 C4 亦相同，即兩單元在邊界具有同一撓度函數w。xyzbbaa土木工程學院有限單元法有限單元法P-30/64 不能保證相鄰單元在公共邊界上法線轉角的連續(xù)性 仍以1-2邊界為例，將y=-b代入后，此時342321xdxdxddx 但對x來講，1, 2結點只能提供2個已知條件，不能完全確定上式中的四個待定常數，故邊界的法線轉角不能保證連續(xù)性。 因此，這種單元是非協調元，但

16、可以驗證這種非協調元是能通過分片試驗的。因而當單元劃分不斷縮小時，計算結果仍能收斂于精確解。 分片試驗可以看作：在常應變情況下，位移的不協調部分對勢能有無貢獻 土木工程學院有限單元法有限單元法P-31/64三、形函數和形函數矩陣三、形函數和形函數矩陣 分別將單元結點1, 2, 3, 4的坐標值代入(4-2-1) ，并事先求出 ,便可得到各結點的位移值。一共可得12個關于i 的方程組，聯立求解可得i 。 如果采用局部坐標（，）并用單元節(jié)點位移值來表示位移函數，則 ywxxwy4, 3 , 2, 1),(iyiyixixiiiNNwNw土木工程學院有限單元法有限單元法P-32/64形函數矩陣: x

17、NxNyNyNNNNNNNNyyyxyx/4141444111 edNw (4-2-2) 4, 3 , 2, 1iyiyixixiiixyNyNwyN4, 3 , 2, 1iyiyixixiiiyxNxNwxN寫成矩陣形式: 土木工程學院有限單元法有限單元法P-33/64式中形函數:2221181iiiiiN211181iiixibN211181iiiyiaN(i=1 2 3 4) (4-2-3)在上面的推導中，我們仍然選用了局部坐標(無因次坐標)。局部坐標與整體坐標的關系為: 01xxa01yybz 土木工程學院有限單元法有限單元法P-34/64四、單元的幾何矩陣四、單元的幾何矩陣B和內力矩

18、陣和內力矩陣S1幾何矩陣B由前可知 ,將（4-2-2）的w代入,得到幾何矩陣： z yxNyxNyNyNxNxNxNxNByyyyx421224221224221221221222（4-2-4） B以子塊形式表示： B=B1 B2 B3 B4。 edBzz則應變可以表示為：土木工程學院有限單元法有限單元法P-35/64土木工程學院有限單元法有限單元法P-36/64土木工程學院有限單元法有限單元法P-37/642.內力矩陣內力矩陣S由基本方程（4-1-5） 可得到： eedSdBDF0（4-2-6） S稱為內力矩陣，把單元的四個結點坐標分別代入(4-2-4)，求得B后，即可獲得S ,內力矩陣 S

19、的顯式如下： DtF123其中薄板彎曲問題的彈性矩陣：21000101)1 (12230EtD2100010112ED內力向量，和結點力向量不同 xyyxMMMF土木工程學院有限單元法有限單元法P-38/64)1 (1223Et土木工程學院有限單元法有限單元法P-39/64土木工程學院有限單元法有限單元法P-40/64五、單元剛度矩陣五、單元剛度矩陣由一般公式得: aabbTdxdyBDBK 0將幾何矩陣B和彈性矩陣D的表達式代入，積分可得薄板彎曲問題矩形單元的單元剛度矩陣的顯式公式： 土木工程學院有限單元法有限單元法P-41/643659111416192124810131518201111

20、07161512212017365192114162418201315121201716151236591124810111107365241000000000000030kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkabDK土木工程學院有限單元法有限單元法P-42/64土木工程學院有限單元法有限單元法P-43/64土木工程學院有限單元法有限單元法P-44/64六、荷載等效變換六、荷載等效變換由荷載等效變換的一般公式可得 aabbTdxdyyxFyxNR,),(1法向均布荷載q 1111,),(ddFyxNa

21、bTTababababqabR331331331331TqF00,代入上述公式得：土木工程學院有限單元法有限單元法P-45/642單元內部點受法向集中力(或力偶)ababababPR22228如：單元中心點受法向集中力 aabbPPTdxdyyxPyxNR,),(),(),(PNPyxNTPPTPPTPP00 代入上述公式可得：土木工程學院有限單元法有限單元法P-46/64七、位移邊界條件七、位移邊界條件 對稱、固定邊和簡支邊上支點的已知位移條件如下：對稱軸: 法線轉角=0固定邊: 撓度=0 (或已知值) 邊線轉角=0 (或已知值) 法線轉角=0 (或已知值) 簡支邊: 撓度=0 (或已知值)

22、 邊線轉角=0 (或已知值)自由邊上節(jié)點的撓度、邊線和法線轉角均為特定參數，同內部節(jié)點一樣。與板鉸接的固定立柱，其節(jié)點撓度w = 0，也可以是已知值。 土木工程學院有限單元法有限單元法P-47/64八、計算例題例題1 計算圖示四邊固定方板。 方板的邊長為l，厚度為t，彈性模型量為E，波松比=0.3，全板承受均布法向荷載，求薄板中的撓度和內力。 單元劃分：單元劃分： 為了說明解題方法，采用最簡單的網格22，即把方板分成四個矩形單元。由于對稱性，只需計算一個單元，例如，計算圖中有陰影的單元，單元的節(jié)點編號為，。土木工程學院有限單元法有限單元法P-48/64此時，單元的a, b是:a=b=l/4 計

23、算節(jié)點荷載計算節(jié)點荷載: 由前面的均布荷載計算公式得：TllllllllqlR 21 12 12 121922邊界條件：邊界條件：邊界23和34為固定邊，因此節(jié)點2, 3, 4的撓度、邊線和法線轉角均為零。邊界12和14為對稱軸，因此x1 =0、y1 =0。于是，在4個節(jié)點12個位移分量中，只有一個待求的未知量w1 。土木工程學院有限單元法有限單元法P-49/64解出 16)681(15815821201120qlwlDwklD04100148. 0Dqlw 323009158. 0)1 (12EtEtD其中： 結構的代數方程組：結構的代數方程組：這是一個單元的計算題目，單元剛度矩陣在此處即為

24、總剛度矩陣。引入支承條件后，在總剛度矩陣中只取第一行、列元素，在方程組右端項中只保留第一個元素。于是結構的代數方程為：土木工程學院有限單元法有限單元法P-50/64內力內力: 利用式(4-2-6)可求得單元4個節(jié)點力矩為：位移向量僅有w1土木工程學院有限單元法有限單元法P-51/64由表看出，網格越密，計算結果越接近于精確答案。還可看出，位移的精度一般比內力的精度高，這是因為在位移法中，位移是由基本方程直接求出的，而內力則是根據位移間接求出的。 土木工程學院有限單元法有限單元法P-52/644.3 薄板有限元程序設計薄板有限元程序設計 一、總框圖一、總框圖 根據彎曲板有限元分析方法的解題過程，

25、可寫出其總框圖如下：土木工程學院有限單元法有限單元法P-53/64二、子框圖二、子框圖、單元坐標結點編號及單剛形式為了取撓度向下為正，又能與前述坐標系統(tǒng)統(tǒng)一，特將前述坐標前翻180 xzy土木工程學院有限單元法有限單元法P-54/64程序中單剛數組為 DK(20, 20), 子程序: Subroutine DG(A, B, E, T, U)為其形成單剛的子程序。 為了能適用板的彈性分析,程序采用了應力元和彎曲元的組合形式，即每個結點考慮5個位移分量: u, v, w, x , y , 前2個為平面應力問題的未知量,后3個為彎曲板的結點未知量。當只作彈性分析時，平面應力元和彎曲元是非藕連的，即單

26、剛的兩個副塊恒為0, 單剛的形式為：土木工程學院有限單元法有限單元法P-55/64、自動形成單元編號信息（ 單元信息數組：IB）。、結點定位向量。、形成荷載列向量。（ a. 結點力； b. 非結點力(只考慮均布力)）、總剛，Subroutine ZG(M, N, LD, A, B, E, T,U)、解方程。 FJZG( ), HUD( )、算單元力。 Subroutine DYL( ) 、算等效結點力。 、彈性矩陣D 。10、幾何矩陣B。 土木工程學院有限單元法有限單元法P-56/64三、輸入數據說明、總信息。共11個(見程序)、結點約束信息數組 JBJB(I,1) i結點的結點號JB(I,2

27、) i結點的約束分量號(15)結點約束信息應根據支承條件或對稱條件決定，如算例中所給出的四邊簡支方板，承受滿布均布力，此時可只取板的1/4作為分析對象，如下圖只取右上角1/4板，采用66網絡，則每個單元的邊長為1米(a=0.5, b=0.5)。土木工程學院有限單元法有限單元法P-57/64設結點編號如圖示: 在y=0的邊界上(1-7結點): 撓度 w=0 (第3個分量) 繞y軸轉角y =0 (第5個分量) 同理,在平行于y軸的x=6m邊界上: w=0, x=0 (3,4分量) 土木工程學院有限單元法有限單元法P-58/64在對稱軸x=0 邊界上 u=0, y =0在對稱軸y=6m 邊界上 v=0, x =0中點(43點)除W外, 其余均=0同時，