弹性圆柱薄壳稳定性问题的不完全双二次非协调板元解

收稿日期:2004-10-11弹性圆柱薄壳稳定性问题的不完全双二次非协调板元解金坚明,高忠社,肖　强(西北师范大学数学与信息科学学院,甘肃兰州　730070)摘　要:　介绍了弹性圆柱薄壳稳定性问题,给出了不完全双二次非协调板元,并用不完全双二次非协调板元求解弹性圆柱薄壳特征值问题.关键词:　弹性圆柱薄壳;稳定性问题;总势能;不完全双二次非协调板元中图分类号:　O 357;O 241 文献标识码:　A 文章编号:100420366(2005)022*******I nco m plete Bi quadra ti c Nonconfor m i n g Pla te Ele m en t Soluti on to theProble m of St abil ity of an Ela sti c C i rcular Cyl i n dr i ca l ShellJ I N J ian 2m ing ,GAO Zhong 2she ,X I A O Q iang(Colleg e of M athe m atics and Inf or m ation S cience ,N orthw est N or m al U niversity ,L anzhou 730070,China )Abstract :　T he p roble m of stability of an elastic circular cylindrical shell w as introduced si m p ly ,and then the incomp lete biquadratic nonconfor m ing p late ele m en t w as given ,and the eigenvalue p roble m of elastic circular cylindrical shell w as s o lved by the incomp lete biquadratic nonconfo r m ing p late ele m ent .Key words :　elastic circular cylindrical shell ;the p roble m of stability ;total potential energy ;incomp lete biq 2uadratic nonconfor m ing p late ele m en t图1　圆柱薄壳受力绕中心轴的扭矩作用1　圆柱壳体的稳定性理论设有一各向同性,均匀材料的圆柱薄壳.它的厚度t 远小于其半径r 和长度L ,壳体同时承受轴向力p ,横向静水压力q 和绕中心轴之扭矩M k 的作用,如图1所示.由于薄壳厚度相当小,显然在丧失稳定之前,圆柱壳体内应力皆为沿厚度不变之均匀薄膜应力,即应力分量可表达为Ρx =p2Πrt=常量,(1)Ρy =-q rt=常量,(2)Σx y =Mk2Πr 2t=常量,(3)其中:p >0表示拉力,q >0为静水外压力.当壳体开始丧失稳定时,壳体的总势能U 应为U =U u +U c ,(4)第17卷　第2期2005年6月 甘肃科学学报Journal of Gansu Sciences V ol .17　N o .2Jun .2005其中:U u 为弯曲应变势能,U c 为薄膜应变势能.如果我们忽略曲率之高阶影响,则内力矩和中曲面的曲率变化关系和平板一样,因而U u 应当与平板之弯曲应变势能相同,即[1]U u =D2∫L 0∫2Πr{[52Ξ5x 2+52Ξ5y 2]2+2(1-Λ)[(52Ξ5x 5y )2-52Ξ5x 252Ξ5y 2]}d x d y ,(5)薄膜应变势能为U c =t2E∫L∫2Πr{[52F x 2+52F y 2]2+2(1-Λ)[(52F 5x 5y )2-52F 5x 252F5y 2]}d x d y ,(6)其中:D -E t312(1-Λ2)圆柱薄壳刚度;E ,Λ,F 材料的弹性模量,泊松比,应力函数1最后,得到总应变势能的表达式 U =U u +U c =D 2∫L 0∫2Πr 0{[52Ξ5x 2+52Ξ5y 2]2+2(1-Λ)[(52Ξ5x 5y )2-52Ξ5x 252Ξ5y 2]}d x d y +t 2E ∫L 0∫2Πr 0{[52F 5x 2+52F 5y 2]2+2(1+Λ)[(52F 5x 5y )2-52F 5x 252F5y 2]}d x d y ,(7)其中应力函数F 与挠曲函数Ξ之间,按应变连续条件存在下面关系:4F =-E r ・52Ξ5x2.(8)在上述受载状态下,应力函数之通解可设(F )通=c 1x 2+c 2x y +c 3y 2,(9)其中:c 1,c 2,c 3为应力状态系数.根据应力函数和应力分量的微分关系,不难看出它的物理意义.c 1=12Ρy =-qr2t,(10)c 2=-Σx y =-Mk2Πr 2t,(11)c 3=12Ρx =p 4Πrt,(12)设V 为圆柱壳体系统的总势能,则V =U -W ,(13)其中W =Ρx t ∫L 0∫2Πr0(-12)(5Ξ5x)2d x d y +q ∫L 0∫2Πr0r2(5Ξ5y )2d x d y +Σx y t∫L 0∫2Πr 0[-(5Ξ5x )・(5Ξ5y)]d x d y +W 0,(14)此处W 0在相应载荷情况下为一常量,,系统处于随遇平衡状态,它的总势能应具有极值,即∆V =0.(15)2　仅受轴向均布压力作用下圆柱薄壳的稳定性:(Ρy =Σxy =0,Ρx ≠0)此时,式(10)～(12)为c 1=c 2=0.c 2=12Ρx =84Πrt,于是 (F )通=c 3y 2.则有U c =t 2E ∫L 0∫2Πr 0{[2c 3]2}d x d y =t 2E ∫L 0∫2Πr 0(4Πrt )2d x d y =8t E Π2r 2t2L ・2Πt =16L E Πrt ,于是U =U u +U c =D2∫L 0∫2Πr{[52Ξx 2+52Ξy 2]2+2(1-Λ)[(52Ξx y )2-52Ξx 252Ξy2]}d x d y +16L E Πrt .这时W 为2 甘肃科学学报 2005年　第2期W =Ρx t∫L 0∫2Πr(-12)(5Ξ5x)2d x d y +W 0.最后弹性圆柱薄壳体仅受轴向均布压力作用时系统的总势能为 V =D 2∫L 0∫2Πr 0{[52Ξ5x 2+52Ξ5y 2]2+2(1-Λ)[(52Ξ5x 5y )2-52Ξ5x 252Ξ5y2]}d x d y +16L E Πrt -Ρx t ∫L 0∫2Πr0(-12)(5Ξ5x)2d x d y -W 0.(16)此时圆柱薄壳仅受轴向均布压力作用,且两端简支时问题的弱解,就是求解Ξ∈H 。

2(A ),使得a (Ξ,Ξ3)=b (Ξ,Ξ3),ΠΞ3∈H 。

2(A ),(17)其中a (Ξ,Ξ3)=D ∫∫A{∃Ξ∃Ξ3+(1-Λ)[252Ξ5x 5y 52Ξ35x 5y -52Ξ5x 252Ξ35y 2-52Ξ5y 252Ξ35x 2]}d x d y ,(18)b (Ξ,Ξ)=-Ρx∫∫A5Ξ5x 5Ξ35xd x d y ,A ={(x ,y ) 0≤x ≤L ,0≤y ≤2Πr },两端简支即:Ξ 5A =0,M n 5A =0,5A 是弹性圆柱薄壳x =0,x =L 时两端边界,n 为边界外向法失,M k 为弯矩.通过与文献[2]类同的分析可知此问题的弱解是存在且唯一的.3　8个自由度不完全双二次非协调板元对给定的区域A 作矩形剖分,e (A 1,A 2,A 3,A 4)为其上任一单元,对给定的Ξe (x ,y )求其上的不完全双二次插值∏Ξe=a 1+a 2x +a 3y +a 4x2+a 5x y +a 6y 2+a 7x 2y +a 8x y 2,(19)满足 Ξ(A i )=(∏Ξe )(A i ),5Ξ5n i (B i )=(55n i ∏Ξe )(B i ),i =1,2,3,4.(20)其中:B i 为e 的各边中点,n οi 为e 在B i 处的单位外法向量,引进局部坐标系Ν=(x -x c ) L 1,Γ=(y -y c ) L 2,(x c ,y c )是e 的形心坐标,L 1,L 2为e 的平行于x ,y 轴边的半长,e 在局部坐标系中的象e δ是边长为2的标准正方形,经过坐标变换可得∏Ξe=∑4i =1Ξ(Ai)N i (Ν,Γ)+∑i =1,3L2(5Ξn i)(B i )M i (Ν,Γ)+∑i =2,4L1(5Ξn i)(B i )M i (Ν,Γ),(21)其中 N 1(Ν,Γ)=14(1+ΝΓ-Ν2Γ-ΝΓ2), N 2(Ν,Γ)=14(1-ΝΓ-Ν2Γ+ΝΓ2), N 3(Ν,Γ)=14(1+ΝΓ+Ν2Γ+ΝΓ2), N 4(Ν,Γ)=14(1-ΝΓ+Ν2Γ-ΝΓ2), M 1(Ν,Γ)=14(-1-2Γ+Γ2+2Ν2Γ), M 2(Ν,Γ)=14(-1+2Ν+Ν2-2ΝΓ2), M 3(Ν,Γ)=14(-1+2Γ+Γ2-2Ν2Γ), M 4(Ν,Γ)=14(-1-2Ν+Ν2+2ΝΓ2).(22)3第17卷 金坚明等:弹性圆柱薄壳稳定性问题的不完全双二次非协调板元解 以插值式(21)作为位移模式的有限元称为不完全双二次元.由于∏Ξe 在相邻单元的共同边界上不连续,所以作为板元是非协调的.但由文献[3～6]可知它能通过分片检验,且收敛速度跟M o rl oy 元相同.4　有限元法第2部分已说明对各向同性,均匀材料的厚度远小于其半径r 和长度L 的圆柱薄壳,当仅受轴向均布压力作用,且两端简支时,问题的弱解存在且唯一.文中第3部分又详细介绍了不完全双二次非协调板元.现在我们就来讨论上述问题的不完全双二次非协调板元解.令{∆}e i =(Ξ(A 1),Ξ(A 2),Ξ(A 3),Ξ(A 4),L 2(5Ξ5n 1)(B 1),L 1(5Ξ5n 2)(B 2),L 2(5Ξ5n 3)(B 3),L 1(5Ξ5n 4)(B 4))T ,[N ]e i =14(N 1(Ν,Γ),N 2(Ν,Γ),N 3(Ν,Γ),N 4(Ν,Γ),M 1(Ν,Γ),M 2(Ν,Γ),M 3(Ν,Γ),M 4(Ν,Γ)).(23)从而可得Ξ=[N ]e i {∆}e i ,Ξ3=[N ]e i {∆3}e i.(24)这样广义应变向量为{ς1}e i =-52Ξ5x 2-52Ξ5y 2-252Ξ5x 5y=-52[N ]e i 5x 2-52[N ]e i 5y 2-252[N ]e i5x 5y{∆}e i =[B 1]e i {∆}e i.(25){ς2}e i =[5Ξ5x ]=[5[N ]e i5x]{∆}e i =[B 2]e i {∆}e i.(26)这里[B 1]e i 为3×8矩阵,[B 2]e i 为1×8矩阵.下面引出单元刚度矩阵[D ]e i =E t312(1-Λ2)1Λ0Λ1001-Λ2,{M }e i =[D ]e i {ς1}e i ,(27)这样定解问题的虚功方程可表示为单元虚功的叠加∑N Ei =1∫∫e i{ς31}T e i{M }e i d x d y =∑N Ei =1-Ρx∫∫e i{ς32}Te i {ς2}e i d x d y ,即∑N Ei =1∫∫e i{ς31}T e i[D ]e i {ς1}e i d x d y =∑N Ei =1-Ρx∫∫e i{ς32}Te i {ς2}e i d x d y ,(28)其中:N E 是矩阵A 剖分成的矩阵单元个数,将(25)代入上式左边可得 ∫∫e i{ς31}Te i [D ]e i {ς1}e i d x d y =∫∫e i([B 1]e i {∆3}e i )T [D ]e i ([B 1]e i {∆}e i )d x d y ={∆3}e i (∫∫e i([B 1]T e i [D ]e i [B 1]e i d x d y ){∆}e i ,令[K ]e i =∫∫e i[B 1]T e i [D ]e i [B 1]e i d x d y ,(29)这是一个8×8的方阵,称为单元刚度矩阵.此外-Ρx∫∫e i{ς32}Te i [ς2]e i d x d y =-Ρx∫∫e i{∆3}T e i ([B 2]Te i [B 2]e i {∆}e i )d x d y ={∆3}T e i (∫∫e i-Ρx [B 2]T e i [B 2]e i d x d y ){∆}e i ,4 甘肃科学学报 2005年　第2期令[F ]e i =∫∫e i-Ρx [B 2]T e i [B 2]e i d x d y ,(30)这是一个8×8的方阵,称为单元载荷矩阵.要得到一个单元刚度矩阵,单元载荷矩阵,只需要计算出 ∫∫e i[5[N ]e i 5x ]T [5[N ]e i 5x ]d x d y ,∫∫ei[52[N ]e i 5x 2]T [52[N ]e i5x 2]d x d y , ∫∫ei[52[N ]e i 5x 5y ]T [52[N ]e i5x 5y ]d x d y ,∫∫e i[52[N ]e i 5y 2]T [52[N ]e i5y 2]d x d y 即可,经过具体计算得到∫∫e i[5[N ]e i 5x ]T [5[N ]e i 5x]d x d y =L 1L 2720L 21　176-16-56-104-460-48　160　48-16　176-104-56-160　48　160-48-56-104　176-16　160　48-160-48-104-56-16　176　160-48-160　48-460-160　160　160　320 0-320 0-48　48　48-48 0-96 0　576　160　160-160-160-320 0　320 0　48-48-48　48 0　576 0-96,(31)∫∫e i[52[N ]e i 5x 2]T [52[N ]e i 5x2]d x d y =L 1L 2720L 21L 22　240　240-240-240-480 0　480 0 240　240-240-240-480 0　480 0　-240-240　240　240　480 0-480 0　-240-240　240　240　480 0-480 0　-480-480　480　480　960 0-960 0　0　0　0　0　0720　0720　　480　480-480-480-960 0　960 0　0　0　0　0　0720　0720,(32)∫∫e i[52[N ]e i 5y 2]T [52[N ]e i 5y 2]d x d y=L 1L 2720L 21L 22　240-240-240　240 0　480 0-480　-240　240　240-240 0-480 0　480　-240　240　240-240 0-480 0　480 240-240-240　240 0　480 0-480 0 0 0 0　720 0　720 0　480-480　480　480 0　960 0-960 0 0 0 0　720 0　720 0-480　480　480-480 0-960 0　960,(33)∫∫e i[52[N ]e i 5x 5y ]T [52[N ]e i 5x 5y ]d x d y =L 1L 2720L 21L 22　480 0-480 0-480　480　480-480 0　480 0-480-480-480　480　480-480 0　480 0　480-480-480　480 0-480 0　480　480　480-480-480-480-480　480　480　960 0-960 0　480-480-480　480 0　960 0　960　480　480-480-480-960 0　960 0-480　480　480-480 0　960 0　960,(34)5第17卷 金坚明等:弹性圆柱薄壳稳定性问题的不完全双二次非协调板元解 把单元刚度矩阵,单元载荷矩阵代入虚功方程式(28)得到∑N Ei =1{∆3}Te i[K ]e i {∆}e i =∑N Ei =1{∆3}T e i [F ]e i {∆}e i ,则有{∆3}T [K ]{∆}={∆3}T [F ]e i {∆},(35)式中:{∆}是由Ξ(x ,y )在节点处的函数值Ξ(A i )及在各剖分单元,各边中点外法向导数值L 1(5Ξ5n i)(B i )(这里的n οi 方向与x 轴方向相同),L 2(5Ξ5n i)(B i )(这里的n οi 方向与y 轴方向相同)所组成的向量,组成的顺序与e i 单元中相仿.[K ]称为总刚度矩阵,它由单元刚度矩阵元素块[K ]e i ,按其脚标对号入座形成.[F ]称为总载荷矩阵,它由单元载荷矩阵元素块[F ]e i ,按其脚标对号入座形成.在边界5A 上,自然边界条件不必列出,对强加的边界条件,在x =0,L 处当位移为零,转角为零时只需在总刚度矩阵及总载荷矩阵中划去边界节点函数值及边界单元的边界中点外法向导数值所涉及的行列即可,这时式(35)成{∆υ3}T [K ]{∆υ}={∆υ3}T [F ]{∆υ},(36)由虚功原理可得([K]-[F ]){∆υ}=0,(37)这样通过求线性齐次方程组(37)的非零解,令det ([K]-[F ])=0,就可求出临界力Ρx .5　计算实例例　两端简支圆柱薄壳受均匀侧压的屈曲分析.几何数据如下:纵向长度L =200c m ,r =100c m ,厚度t =1c m .圆柱薄壳的材质是铝合金,弹性模量E =71GPa ,Λ=0.33,求轴向均布压力的临界特征值.解　已知当圆柱薄壳的厚度远小于半径r 和长度L 时,且仅受轴向均布压力作用,此时压应力临界特征值无量纲系数公式和压应力临界特征公式分别为Ρ3x=13(1-Λ2),(38)Ρx =E3(1-Λ2)(t r)2,(39)其中:Ρ3x 是压应力临界特征值无量纲系数,由于Λ=0.33得Ρ3x =0.61,代入式(39)得Ρx =0.43GPa .当网格剖分为10×20用上面介绍的不完全双二次非协调板元求解可得Ρx =・0.435GPa 二者相差十分小.参考文献:[1]　铁摩辛柯S ,沃诺斯基S .板壳理论[M ].板壳理论翻译组1北京:科学出版社,1977.[2]　金坚明.薄板特征值问题误差分析[J ].兰州大学学报,2002,38(6):27232.[3]　Stumm el F .Basic Pompactness P roperties of N onconfor m ing and H ybrid F inite E le m ent Spaces [J ].Rairo A nalN um er ,1980,4(1):812115.[4]　C iarlet P G .The F inite E le m ent M ethod for E lli p tic P roble m [M ].N orth Holland ,Am sterda m ,1978.[5]　(日)鹫津久一郎.弹性和塑性力学中的变分法[M ].卢文达,黄择方,卢鼎霍译.北京:科学出版社,1984.[6]　石钟慈.关于不完全双二次非协调板元的收敛性[J ].计算数学,1986,8(1):53262.作者简介:金坚明,(19432)男,浙江省萧山人,1967年毕业于兰州大学数学系,现任西北师范大学数学与信息科学学院教授、研究生导师,主要研究方向为样条函数、固体力学有限元、小波分析.6 甘肃科学学报 2005年　第2期。