基于有限元分析板时厚跨比的选择

基于有限元分析板时厚跨比的选择

板和壳在生产和建设中得到广泛使用。随着科学技术的发展，除传统的薄膜（外壳）外，多个特殊的板壳结构也被用于。例如，气冷核反应堆的预测混凝土压力容器、各种中间板（外壳）和环保板（外壳）。在工程中，广泛应用了kirchhoff-hele-hele-hele-hele-hele-hele-hele-hele的薄板理论和薄板的薄板理论。在实际工程中，由于结构几何形状、边界条件和负荷的复杂性，使用分析方法几乎是不可能的。其中大多数使用各种数值方法进行分析。然而，由于有限的元法只需划分区域中的单元，所需分解的精度和计算成本较大，因此大多数研究人员开始研究边界积分方程（biem）。在边界分析方法的介绍下，边界分析单元被引入不同的区域，其中f.vadweec.a.brebas，j.c.f.teles在边界分析理论中研究了厚板的边界方法。本文采用Reissner中厚板理论边界元程序、Mindlin中厚板理论有限元程序和三维弹性力学有限元程序计算各种厚跨比下四边固支和四边简支矩形板,并进行对比,来说明三维理论、中厚板理论和薄板理论分析板时厚跨比的适用范围以及精确程度,为采用何种理论来分析板提供理论依据和参考;同时也验证了采用缩减积分的Mindlin中厚板有限元法消除了剪切锁死现象.1厚板控制方程Reissner中厚板理论是考虑了横向剪切变形和挤压变形的理论.考虑笛卡尔坐标系(x1,x2,x3),其中x1,x2轴在板中面内,x3则垂直于板中面,如图1所示.在板的上下表面切应力为零而正应力为分布载荷集度,即x3=±h/2,σ33=±q/2.各应力在厚度方向上的变化关系如下:σij=12x3h3Μij‚σi3=32h[1-(2x3h)2]Qi‚σ33=14(2x3h)[3-(2x3h)2]q.(1)σij=12x3h3Mij‚σi3=32h[1−(2x3h)2]Qi‚σ33=14(2x3h)[3−(2x3h)2]q.(1)式中:i,j=1,2;Mij为单位长度弯矩;Qi为剪力.为了求解Reissner中厚板控制方程的基本解,在此假定体力在厚度方向上也按应力一样变化,即Ρi=Bi(x1,x2)12x3h3‚Ρ3=B3(x1,x2)32h[1-(2x3h)2].(2)Pi=Bi(x1,x2)12x3h3‚P3=B3(x1,x2)32h[1−(2x3h)2].(2)由以上假定,Reissner中厚板的控制方程变为:{Μij,j-Qi=0,Qi,i+q=0.(3)其中的弯矩、扭矩和剪力用转角和位移表示为:Μij=D(1-ν)2(βi,j+βj,i+2ν1-νβk,kδij)+ν1-νqλ2δij,Qi=D(1-ν2)λ2(βi+w,i).(4)式中:D=Eh3/[12(1-ν2)]为板的弯曲刚度;ν为泊松比;λ=√10/h是厚度参数.式(4)中转角和位移为沿板厚的平均值,即:βi=∫h/2-h/212x3hvidx3,w=∫h/2-h/232h[1-(2x3h)2]v3dx3.(5)式中:vi和v3为坐标轴方向的位移分量.为了简便,将βi(i=1,2)和w用ui(i=1,2,3)表示,则对于区域Ω的边界Γ,位移和广义面力的边界条件可以表示为:ui=¯ui,在Su上;ti=¯ti,在St上.其中S=Su∪St,ti为广义面力.假定边界上某点外法线方向余弦为ni=cos(n,xi),则广义面力为ti=Mijnj,t3=Qini.(i,j=1,2)(6)2厚板边界积分基本解可由式(3)中体力项定义为一块无限大中厚板作用单位横向集中力或单位弯矩求得,具体推导见文献.仅列出关于位移和面力的基本解:u*ij=12πD{21-ν[B(z)δij-A(z)r,ir,j]-12δij(lnz-12)-12r,ir,j}‚u*i3=-u*3i=12πD12rr,i(lnz-12)‚u*33=-12πD1λ2[21-νlnz-14z2(lnz-1)].(7)t*ij=-14πr[(4A(z)+2zΚ1(z)+1-ν)(δijr,nnj)+(4A(z)+1-ν)r,jni-2(8A(z)+2zΚ1(z)+1-ν)r,ir,jr,n]‚t*i3=λ22π[B(z)ni-A(z)r,ir,n]‚t*3i=-1-ν8π[(21+ν1-νlnz-1)ni+2r,ir,n]‚t*33=-12πrr,n.(i,j=1,2)(8)上两式中:A(z)=Κ0(z)+2z[Κ1(z)-1z];B(z)=Κ0(z)+1z[Κ1(z)-1z];z=λr;r=[(x1-ξ1)2+(x2-ξ2)2]1/2;r,i=(xi-ξi)/r;r,n=r,ini.其中K0(z)和K1(z)为修正贝塞尔函数.使用Betti功互等定理、格林函数或加权残值法,可以得到Reissner中厚板的边界积分方程:cijuj+∫Γt*ijujdΓ=∫Γu*ijtjdΓ+∫Ω[u*i3-ν(1-ν)λ2u*ik,k]qdΩ‚i,j=1,2,3；k=1,2.(9)式中:cij是一个与边界形状有关的系数,对于光滑边界cij=δij/2,对于有角点和尖点的边界法线不连续的边界,可以解析计算cij,但更方便的是间接计算cij的值.当整个边界无外力和载荷作用时,式(9)变为:cijuj+∫Γt*ijujdΓ=0.(10)上式的非零解为板的3种独立刚体位移的任意组合,把这些位移代入式(10)有cia(ξ)=-∫Γ[t*ia(ξ,x)+(ξa-xa)t*i3(ξ,x)]dΓ,ci3(ξ)=-∫Γt*i3(ξ,x)dΓ,(11)i=1,2,3;α=1,2.当均布荷载情况下,可以应用散度定理把式(9)中的域积分简化为边界积分.数值方法求解式(9)后,未知边界位移和边界力可以全部求出.令cij=δij可以求出内部点的位移,而内部点的内力解,需要对式(9)关于源点ξ求微分,有如下形式:Ρij=ν(1-ν)λ2qδij+∫Γ[Dijktk-Sijkuk+qwij]dΓ.(12)式中:i,j≠3.P的分量为Pij=Mij,P3j=Qj.其中wij,Dijk和Sijk见文献中附录Ⅰ.3标准高斯积分法把所考虑的板的边界划分为M个单元,以每个节点为源点并在M个单元上积分,离散式(9),并将其中域积分转换为边界积分,则可得线性方程组HU=GT+D.(13)式中:H,G为影响矩阵;D为对应于载荷q的边界积分列向量.代入已知的边界条件后,将其中的已知项移到方程等号右边与D向量合并,未知项移到等号左边,得到标准代数方程组AX=B.(14)采用标准高斯消去法求解上述方程组,即可求得未知量.之后可以利用式(9)计算内部任意一点的位移,利用式(12)求内部点的内力.由于式中的奇异性,在计算奇异积分时,采用了三次多项式非线性变换.此种变换可以在奇异点附近集中较多的高斯点,因此对奇异积分采用标准高斯积分能达到相当高的精度.对于存在角点和尖点等边界法线不连续的情况,采用双节点或部分不连续单元.双节点单元中2个节点有不能在同一方向上给定位移的限制,而不连续单元没有这种限制.根据前述离散化的公式,编制了Reissner中厚板边界元的FORTRAN程序.4固支矩形板.采用了2个具有代表性的算例来研究薄板和中厚板理论适用范围和精确程度.算例为:1)受均布荷载作用的四边固支矩形板,如图2所示;2)受均布荷载作用的四边简支矩形板.两算例的长、宽取a=b=1,横向荷载集度q=1,板的弯曲刚度D=1,泊松比ν=0.3.4.1接触模型的建立利用本文编制的Reissner中厚板边界元程序计算了在不同厚跨比(h/a=0.001,0.005,0.01,0.05,0.1,0.2,0.5)下的板中心点的挠度w,板中心点下表面应力σx和板固支边中点的弯矩Mx.并与利用有限元软件ABAQUS中对应的三维有限元解和板壳有限元解,还有部分薄板解析解进行对比.利用ABAQUS的三维8节点单元、二维4节点板壳单元和2节点边界单元进行分析计算,其中以三维8节点有限元分析的结果作为基准,以说明二维4节点板壳单元和2节点边界单元计算结果的准确性.在三维有限元、板壳有限元和中厚板边界元计算中,取中面中心点的挠度和中心点板的下表面点的应力σx作为比较对象.因为三维单元只能得到节点的位移和应力,要想得到截面上的内力时需要进行合成计算,为此取中心点板的下表面点的应力σx作为比较对象.板壳有限元分析结果既可输出内力,也可输出应力,但由于边界元只能输出内力,所以需要利用弯矩Mx与应力σx的关系求得:σx=12Μxh3z.(15)式中:z是应力计算点到中面的距离.为了更好、更精确地对比各方法的计算结果,对Reissner中厚板边界元模型每边划分16个线性单元,整个板共64个单元;对板壳有限元模型同样每边划分16个线性单元,共16×16个单元;三维有限元模型则在板面内每边划分16个线性单元,厚度方向上网格的划分是依照三维单元尽量要成正六面体的准则来进行的,这样能保证作为基准的三维实体有限元结果的精度.图3画出了各种厚跨比下板中点的挠度.从图3中可以看出在所有厚跨比情况下,3种方法的计算结果十分吻合,说明了考虑剪切变形的Reissner中厚板理论是对三维弹性力学的较精确的简化.从图中还可以看出,当厚跨比大于0.2时,与三维单元的结果误差明显增大.从图4可看到,在厚跨比较小(h/a<0.2)时,边界元、板壳有限元和三维有限元计算出的σx十分接近,但是在h/a>0.2后,误差越来越大.但是边界元解精度比板壳有限元的解精度高,证明了边界元法在理论上比有限元法应力(内力)的精度高.综合各种厚跨比挠度和应力的误差,中厚板与厚板的厚跨比分界线一般取0.2为宜,厚跨比超过了0.2,无论是挠度还是应力的误差都会急剧增加.同时,从图3和图4还可以看出:利用缩减积分技术的Mindlin板壳有限元,在计算非常薄的板时大大缓解了剪切锁死现象.图5还给出了固支边中点处的弯矩Mx.在计算弯矩时,三维实体单元选取20节点二次单元,因为二次单元的位移是二次分布的,但是应力是线性分布的,所以二次单元积分出的弯矩应该更精确.从图5可以看出边界元解比板壳有限元的解精确得多,而且在薄板阶段和解析解十分靠近,充分说明了中厚板边界元解对薄板有很好的精度.边界元解和有限元三维实体单元解有差异,应当是来源于弯矩计算时积分的误差,因为三维单元只能输出应力,想获得弯矩还需对应力进行积分,而节点并不是在高斯积分点上,只能使用梯形法近似积分,致使精度较差.4.2单元划分及整板划分如图6所示四边简支矩形板,所有数据同例4.1.板壳有限元仍然采用双线性板壳单元,而中厚板边界元采用2节点线性单元.另外单元划分准则也跟算例4.1相同,边界元在每边上划分16个单元,整板共64个单元;板壳有限元也在每边划分16个单元.从图7中可以看出,3种单元的四边简支板中心的挠度的解十分接近.在厚板阶段(h/a>0.1),图中清楚地表明了本文所用中厚板边界元法比板壳有限元法更接近三维有限元的解,说明了边界元法精度比有限元高.从图8中可以看出,无论是边界元还是板壳有限元的解都十分接近三维有限元的解.从图中可以看到,在薄板阶段(h/a<0.1),边界元解和薄板理论十分接近,和三维有限元解、板壳有限元解有些许误差,因为板壳有限元解的内力精度要低于边界元解.5板壳有限元法本文将Reissner中厚板理论运用到边界元法中,并通过数值方法将其编制成FORTRAN程序,通过几个典型算例的计算,证明了边界元法分析中厚板的程序的正确性.另外,通过与三维有限元、板壳有限元的对比,证明了Reissner理论是对三维弹性力学