薄壁箱梁空间计算的板梁单元法

薄壁箱梁空间计算的板梁单元法

0采用设置箱梁补偿原理构造箱梁板件由于薄板结构的结构通常必须承受双向弯曲和扭转，所以有必要产生双向弯曲、畸变、变形和弯曲滞后。用一般方法,如广义坐标法、弹性地基梁比拟法及一般的有限元法等计算薄壁箱梁结构挠曲、扭转、畸变和剪力滞后都是比较复杂和繁冗的。1949年符拉索夫提出了对平板围成的闭口截面杆件考虑截面外形轮廓线变形的约束扭转理论的广义坐标法,这也是当今计算闭口截面薄壁杆件约束扭转的主要方法。但这一方法要求箱梁截面双轴对称,并且不能考虑箱梁顶板与底板的悬臂作用。在70年代,我国著名桥梁专家李国豪院士提出了箱梁挠曲、扭转有限元法,该方法按扭转中心确定扭转翘曲变形沿箱梁横截面的分布,当箱梁横截面不是双轴对称时,这与A.A.乌曼斯基薄壁结构扭转理论不符,畸变中心计算列式的物理意义也不够明确。文献提出了单室箱梁扭转和截面畸变计算的弹性地基梁相似法,该方法提出了计算畸变的良好思路,但其截面畸变翘曲惯性矩JⅡ按扭转中心计算似乎不妥。另外多室箱梁剪切中心和畸变中心的计算都是比较烦琐的,给工程应用带来了诸多不便。本文根据箱梁力学行为为梁特征的特点,将围成闭口截面箱梁的板件视为板梁,即相当于箱梁梁段子单元。这种方法不管箱梁是单室还是多室都可用板梁子单元组拼形成箱梁梁段单元。本文提出的计算方法不需计算箱梁的剪切中心和畸变中心,且计算思路与方法较简洁,物理概念明确,自由度较少,计算精度较高,可用于大型多室箱梁结构的空间计算分析。1箱梁截面位移参数为了说明问题的方便,本节将以单室箱梁为例给出其位移模式。如图1(b)所示,对于单室箱梁,其截面位移参数为四个角点的纵向位移wul、wuR、wll、wlR,箱梁顶板、底板的横向位移uu、ul及左、右上角的竖向位移vl、vR,为了考虑箱梁顶板的剪力滞后,还须在箱梁顶板增加轴向非均匀变形f1sin((x+bu2)πbu)+f2sin((3x+bu2)πbu)[4‚5],这样整个箱梁截面具有特征的变形参数共有10个。若把箱梁梁段单元作为一个整体单元考虑,用10个截面位移参数来描述箱梁两个方向的挠曲、扭转、畸变和剪力滞后,这给选择位移函数带来一定的困难。由于组成薄壁箱梁的各板件共同受力的结果,箱梁的应力主要是以梁的应力特征反映出来,这也说明组成箱梁的各板件是以梁的力学行为为特征的。因此可将箱梁的各板件划分为板梁子单元。这样就很方便地用10个截面位移参数来描述箱梁两个方向的挠曲、扭转、畸变及剪力滞后。对于组成箱梁的各板梁子单元,其横向和竖向位移模式在纵向采用三次多项式、扭转变位和轴向均匀与非均匀纵向变位均采用一次多项式。图2给出了箱梁顶板子单元的坐标布置图。对顶板横向位移uuc=[Ν(z)]{δuu}竖向位移vuc=[Ν(z)]{δvu}纵向位移wuc=[Μ(z)]{δwu}+[Μ(z)]sin((x+bu2)πbu){δf1}+[Μ(z)]sin((3x+bu2)πbu){δf2}扭转位移φu=[Ν(z)]{δφu}}(1)式中,[N(z)]=[N1(z)N2(z)N3(z)N4(z)];Ν1(z)=1-3(zl)2+2(zl)3;Ν2(z)=z-2z2l+z3l2;Ν3(z)=3(zl)2-2(zl)3;Ν4(z)=-z2l+z3l2;[M(z)]=[M1(z)M2(z)];Μ1(z)=1-zl;Μ2(z)=zl;{δuu}e=[uuciuuci´uucjuucj´]Τ;{δvu}e=[vucivuci´vucjvucj´]Τ;{δwu}e=[wuciuucj]Τ;{δf1}=[f1if1j]Τ;{δf2}=[f2if2j]Τ;{δφu}e=[φuiφuj]Τ上述式中,{δuu}、{δvu}、{δwu}、{δφu}指箱梁顶板板梁子单元质心处横向变位、竖向变位、轴向均匀变位及扭转变位的节点位移参数,{δf1}、{δf2}为考虑箱梁顶板剪滞效应而增加的非均匀变形的正弦波峰值。对底板ulc=[Ν(z)]{δul}vlc=[Ν(z)]{δvl}wlc=[Μ(z)]{δwl}φl=[Μ(z)]{δφl}}(2){δul}=[ulciulci´ulcjulcj´]Τ{δvl}=[vlcivlci´vlcjvlcj´]Τ{δwl}=[wlciwlcj]Τ{δφl}=[φliφlj]Τ右腹板ufR=[Ν(z)]{δufR}vfR=[Ν(z)]{δvfR}wfR=[Μ(z)]{δwfR}φfR=[Μ(z)]{δφfR}}(3){δufR}=[ufRiufRi´ufRjufRj´]Τ{δvfR}=[vfRivfRi´vfRjvfRj´]Τ{δwfR}=[wfRiwfRj]Τ{δφfR}=[φfRiφfRj]Τ左腹板ufl=[Ν(z)]{δufl}Τvfl=[Ν(z)]{δvfl}Τwfl=[Μ(z)]{δwfl}Τφfl=[Μ(z)]{δφfl}Τ}(4){δufl}=[ufliufli´ufljuflj´]Τ{δvfl}=[vflivfli´vfljvflj´]Τ{δwfl}=[wfliwflj]Τ{δφfl}=[φfliφflj]Τ2[wlcvucvlc]本文第1节给出了箱梁截面位移参数及其板梁子单元的位移模式,为了组拼箱梁各板件的刚度矩阵,并以箱梁截面位移参数作为节点自由度,于是需要建立箱梁截面位移参数与其板梁子单元位移参数的关系。对于箱梁板件的板梁子单元,由梁单元的基本假定,不考虑板梁单元纵向纤维的挤压,可略去箱梁板件宽度与厚度方向纤维挤压变形,于是左腹板下部竖向位移为vlt=vl+tanα(uu-ul),右腹板下部竖向位移为vRt=vR+tanα(ul-uu)。根据组成箱梁各板件板梁子单元之间的变形协调关系可得到箱梁截面位移参数与其板梁子单元位移参数关系:对顶板有{du}=A{d}(5)对底板有{dl}=B{d}(6)对左腹板有{dfl}=C{d}(7)对右腹板有{dfR}=D{d}(8)式中{du}=[wucuucvucuuc′vuc′φuf1f2]T{dl}=[wlculcvlculc′vlc′φl]T{dfl}=[wfluflvflufl′vfl′φfl]T{dfR}=[wfRufRvfRufR′vfR′φfR]T{d}=[uuulvlvRwulwuRwllwlRf1f2]TA=[Aij]8×10A1,5=A1,6=12;A2,1=1;A3,3=A3,4=12;A4,5=-A4,6=1b;A5,5=A5,6=-A5,7=-A5,8=12h;A5,9=sinα1+sinα32h;A5,10=sinα2+sinα42h;A6,3=-A6,4=-1b;A7,9=A8,10=1B=[Bij]6×10B1,7=B1,8=12;B2,2=1;B3,3=A3,4=12;B4,7=-B4,8=1bl;B5,5=B5,6=-B5,7=-B5,8=12h;B5,9=sinα1+sinα32h;B5,10=sinα2+sinα42h;B6,1=-B6,2=-2btanα;B6,3=-B6,4=-1blC=[Cij]6×10C1,5=C1,7=12;C1,9=sinα12;C1,10=sinα22;C2,1=sinα;C2,3=cosα;C3,1=cos2α2cosα;C3,2=12cosα;C3,3=-sinα;C4,5=-C4,7=1hb;C4,9=sinα1hb;C4,10=sinα2hb;C5,5=12bcosα-sinαh;C5,6=-12bcosα;C5,7=12blcosα+sinαh;C5,8=-12blcosα;C5,9=-sinαsinα1h;C5,10=-sinαsinα2h;C6,1=-C6,2=1hD=[Dij]6×10D1,6=D1,8=12;D1,9=sinα32;D1,10=sinα42;D2,1=-sinα;D2,4=cosα;D3,1=-cos2α2cosα;D3,2=-12cosα;D3,4=sinα;D4,6=-D4,8=1hb;D4,9=sinα3hb;D4,10=sinα4hb;D5,5=3sin2α-cos2α2bcosα;D5,6=sinαh-3sin2α-cos2α2bcosα;D5,7=-12blcosα;D5,8=12blcosα-sinαh;D5,9=-sinαsinα3h;D5,10=-sinαsinα4h;D6,1=-D6,2=1h其余系数Ai,j、Bi,j、Ci,j、Di,j(i=1,2,…,6或8,j=1,2,…,10)为零。3[[[][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][][]][][]][]][][][][][][]][][][]][]][]]][]][]]dz]d[][]][]][]][][][][][][]][][][]][]][]][]][][][][]][]][][][][][][][][][][][][][][][][][][π=πu+πl+πfl+πfR-{F}eT[{di}T{dj}T]eT(9)式中,πu、πl、πfl及πfR分别为薄壁箱梁单元顶板板梁、底板板梁、左腹板板梁及右腹板板梁子单元的应变能。即πu=12∫ΩuE([Μ′(z)]{δwuc}+sin((x+bu2)πbu)[Μ′(z)]{δf1}+sin((3x+3bu2)πbu)[Μ′(z)]{δf2}-x[Ν″(z)]{δuu}-y[Ν″(z)]{δvu})2dΩ+12∫l0GJu[Μ′(z)]{δφu}[Μ′(z)]{δφu}dz=12[{δwu}Τ{δf1}Τ{δf2}Τ{δuu}Τ{δvu}Τ][ΚΤu][{δwu}Τ{δf1}Τ{δf2}Τ{δuu}Τ{δvu}Τ]Τ+12∫l0GJu{δφu}Τ[Μ′(z)]Τ[Μ′(z)]{δφu}dz=12{du}Τ[Κu]{du}(10)式中,[KΤu]=[KΤuij];KΤu1,1=EA∫l0[M′(z)]T[M′(z)]dz;ΚΤu2,1=∫ΩuE[Μ′(z)]Τ[Μ′(z)]sin((x+bu2)πbu)dΩ;ΚΤu2,2=∫ΩuE[Μ′(z)]Τ[Μ′(z)]sin2((x+bu2)πbu)dΩ;ΚΤu3,1=∫ΩuE[Μ′(z)]Τ[Μ′(z)]sin((3x+bu2)πbu)dΩ;ΚΤu3,2=∫ΩuE[Μ′(z)]Τ[Μ′(z)]sin((x+bu2)πbu)⋅sin((3x+bu2)πbu)dΩ；ΚΤu3,3=∫ΩuE[Μ′(z)]Τ[Μ′(z)]sin2((3x+bu2)πbu)dΩ；KΤu4,4=EIuy∫l0[N″(z)]T[N″(z)]dz;KΤu5,5=EIux∫l0[N″(z)]T[N″(z)]dz;而[Ku]是由[KΤu]和箱梁顶板扭转刚度矩阵∫l0GJu[M′(z)]T[M′(z)]dz叠加后再按自由度{du}编号重排后而得到的。πl=12∫l0EΙlx[Ν″(z)]{δul}[Ν″(z)]{δul}dz+12∫l0EΙly[Ν″(z)]{δvl}[Ν″(z)]{δvl}dz+12∫l0GJl[Μ′(z)]{δφl}[Μ′(z)]{δφl}dz+12∫l0EA[Μ′(z)]{δwlc}[Μ′(z)]{δwlc}dzπfl=12∫l0EΙflx[Ν″(z)]{δufl}[Ν″(z)]{δufl}dz+12∫l0EΙfly[Ν″(z)]{δvfl}[Ν″(z)]{δvfl}dz+12∫l0GJfl[Μ′(z)]{δφfl}[Μ′(z)]{δφfl}dz+12∫l0EA[Μ′(z)]{δwfl}[Μ′(z)]{δwfl}dzπfR=12∫l0EΙfRx[Ν″(z)]{δufR}[Ν″(z)]{δufR}dz+12∫l0EΙfRy[Ν″(z)]{δvfR}[Ν″(z)]{δvfR}dz+12∫l0GJfR[Μ′(z)]{δφfR}Τ[Μ′(z)]{δφfR}Τdz+12∫l0EA[Μ′(z)]{δwfR}Τ[Μ′(z)]{δwfR}Τdz其中,E、G分别为箱梁的弹性模量和剪切模量,Iux、Iuy、Ilx、Ily、Iflx、Ifly、IfRx、IfRy分别为箱梁顶板、底板、左腹板、右腹板其相应的局部坐标x轴和y轴的惯性矩,Ju、Jl、Jfl、JfR分别为其相应的扭转惯性矩。将πu、πl、πfl、πfR代入式(9),再根据变分原理有([A00A]Τ[Κu][A00A]+[B00B]Τ[Κl][B00B]+[C00C]Τ[Κfl][C00C]+[D00D]Τ[ΚfR][D00D])⋅{{di}{dj}}={F}e(11)式中,[Ku]、[Kl]、[Kfl]、[KfR]分别为箱梁顶板、底板、左腹板及右腹板子单元的单元刚度矩阵,除[Ku]外,这些子单元刚度矩阵与普通空间梁单元刚度矩阵形式相同。其刚度矩阵形式可参阅有关结构力学和有限元方面的书籍,但要注意绕各板局部坐标转角与横向、竖向位移一阶导数的关系及自由度排列序号。令{δ}=[{di}Τ{dj}Τ]Τ[Κe]=[A00A]Τ[Κu][A00A]+[B00B]Τ[Κl][B00B]+[C00C]Τ[Κfl][C00C]+[D00D]Τ[ΚfR][D00D]式(11)变为[K]e{δ}e={F}e(12)4横向部结构变位本文第3节式(12)是箱梁单元的平衡方程,而对整个箱梁结构计算就需要考虑横隔板的作用。箱梁结构中设置横隔板主要是抵抗箱梁横截面畸变。设箱梁畸变角为γ,由图1(b)变位模式有γ=uu-ulh-12bl[vR+vR+tanα(ul-uu)]-[vl+vl+tanα(uu-ul)]=(1h+tanαbl)(uu-ul)-1bl(vR-vl)(13)又设横隔板抗剪刚度为Rd(横隔板顶点产生单位侧向水平位移时作用于顶点的横向水平力),则产生畸变角γ的剪切力为Rdγh,如图3,横隔板剪切变形应变能πRd等于剪切力所做之功,则图3横隔板畸变示意πRd=12Rdh2γ2=12Rdh2((1h+tanα1bl)-1bl(vR-vl))2(14)对式(14)进行变分有δπRd=Rdh2((1h+tanα1bl)(uu-ul)-1bl(vR-vl))⋅((1h+tanα1bl)(δuu-δul)-1bl(δvR-δvl))(15)在对箱梁结构进行有限元离散时,尽可能将横隔板位置划分到箱梁截面节点位置,否则还需根据箱梁中板梁子单元的变位关系来确定横隔板位置处的uu、ul、vR、vl。为了使问题变得简单,设横隔板在箱梁i截面节点处,则式(15)可写成δπRd=Rdh2⋅((1h+tanα1bl)(uui-uli)-1bl(vRi-vli))⋅((1h+tanα1bl)(δuui-δuli)-1bl(δvRi-δvli))=Rdh2[δu