深水桥墩墩-水耦合作用计算方法对比研究

深水桥墩墩-水耦合作用计算方法对比研究

1动水压力计算模式中国许多桥梁位于水深。地震激励下的地板和水坝的相互作用会对桥墩产生动水压（桩水耦合）。这种功能不容忽视。它改变了码头结构的动态特征，改善了结构的动态响应。国内、外抗震设计规范对动水压力的计算还不完善,中国、印度、日本和欧洲等规范对作用在桥墩上的动水压力的计算作了规定,通过对比发现各国规范对于墩-水耦合振动的计算有较大的差异。我国《公路工程抗震设计规范》(JTJ044-89)和《铁路工程抗震设计规范》(GB50111-2006)也给出了当水深超过5m时,动水压力的近似计算公式,但仅限于圆形桥墩和矩形桥墩,且《公路桥梁抗震设计细则》(JTG/TB02-01-2008)1.0.2条指出“本细则主要适用于单跨跨径不超过150m的混凝土梁桥”,而目前我国已建、在建或拟建的深水桥梁大多是单跨跨径超过150m的长大桥梁。因此需要对深水桥梁的动水压力计算模式进行深入研究。目前动水压力计算主要分为解析法、数值分析法和半解析半数值分析法3类。其中数值分析法分为位移(结构)～位移(流体)、位移(结构)～压力(流体)2种模式;半解析半数值法对流体动水压力以解析解或半解析解表示,并将结构大幅简化后使用有限元软件计算或编程求解。国外学者Hisao.Goto等对圆形墩柱动水压力的解析解及影响因素进行了研究,赖伟等采用辐射波浪理论推导了圆形桥墩外域水的附加动水压力,刘振宇、李乔等推导了圆形和矩形空心墩内域水体附加动水压力的计算解析式,杨万理等对地震作用下桥墩非竖直棱面动水应力进行了分析。本文将对目前应用较广的3种计算模式[Morison公式、基于辐射波浪法的半解析半数值解法(本文简称辐射波浪法)和基于流体声单元的有限元数值解法(本文简称流体声单元法)]理论进行探讨,并建立、改进相应的计算模型,通过对常用深水桥墩动力特性以及动水压力的对比分析,得出它们的优缺点及适用范围。23几种计算模式的基本理论2.1内载荷的刚度Morison公式是1950年由Morison等提出的计算墩柱上波浪力的方法。该方法假定当墩柱的尺寸与波长相比较小时,墩柱并不影响波动场,且要求墩柱直径D与波长L之比D/L≤0.15时才较适用。设圆形墩柱固结于水底,水为理想可压缩流体、无漩运动,流体质点与墩柱节点沿X方向的绝对位移分别为uX、ut。则单位高度墩柱上受到的沿X方向的波浪力由惯性力和阻尼力2部分组成:F(X,Ζ,t)=ρπD24üX+CmρπD24(üX-üt)+12CDρD(˙uX-˙u)|˙uX-˙ut|(1)如果是在静水中,可得:F(X,Ζ,t)=-CmρπD24üt-12CDρD˙ut|˙ut|(2)式中,Cm、CD分别为惯性力系数和阻尼系数,通过试验或经验来选取,常取Cm=1.0,CD=1.2;D为桥墩直径;ρ为流体密度。2.2弹性振动引起的动水应力水体初始是静止的,圆形墩柱的运动产生一个向四周辐射的波。设辐射波是小波幅的线性波,并且波浪与墩柱的作用也是线性的。根据拉普拉斯方程及边界条件,求出圆形墩柱刚体运动和弹性振动引起的速度势,并利用速度势与动水应力的关系得到由刚体运动和弹性振动引起的动水应力。动水应力沿着圆柱表面积分,可得刚体运动和弹性振动引起的作用于各水下节点的动水压力,分别为:Ρout(1)jj=ω2u(1)p0eiwtQout(1)jj=-ü(1)pΜout(1)jj-˙u(1)pCout(1)jj(3)Ρout(2)ij=ω2u(2)p0(Ζi)eiwtQout(2)ij=-ü(2)pΜout(2)ij-˙u(2)pCout(2)ij(4)式中,相关变量含义及参数取值见参考文献。2.3品种c:压力+速度+振动假设流体是可压缩的,无粘、无漩、无热交换的均匀流,变形小,流速远小于流体中的声速。根据流体动量方程和连续方程可以得到声波动方程。将声波动方程用散度和梯度表达,经过Galerkin方法离散,乘上压力的变分并在流体域内积分,可得:∫v1c2δΡ∂2Ρ∂t2dv+∫v({L}ΤδΡ)({L}Ρ)dv=∫S{n}ΤδΡ({L}Ρ)dS(5)式中,c为流体介质声速;P为流体压力;t为时间。考虑流体法向压力梯度和墩柱法向加速度之间的关系{n}⋅{∇Ρ}=-ρ0{n}⋅∂2{ue}∂t2,并考虑流固耦合界面的阻尼引起的能量损耗,可得:[ΜΡe]{⋅⋅Ρe}+[CΡe]{˙Ρe}+[ΚΡe]{Ρe}+ρ0[Re]Τ{üe}=0(6)流固耦合界面考虑流体施加给墩柱的压力,墩柱的振动方程为:[Μe]{üe}+[Ce]{˙u}+[Κe]{ue}={Fe}+{Fpre}(7)式中,{Pe}为节点压力;{ue}为桥墩节点位移矢量,需要联立方程(6)、(7)求解,可得:[[Μe][Μfs][Μpe]]{{üe}{⋅⋅Ρe}}+[[Ce][Cpe]]{{˙ue}{˙Ρe}}+[[Κe][Κfs][Κpe]]{{ue}{Ρe}}={{Fe}{0}}(8)式中,[Me]为墩柱质量矩阵;[Ce]为墩柱阻尼矩阵;[Ke]为墩柱刚度矩阵;[Fe]为外部荷载矩阵;{Fpre}为流体压力荷载矩阵;[Mfs]=ρ0[Re]T;[Kfs]=-[Re]。对于非流固耦合界面的流体和固体,按照传统的波动和振动方程来求解。3数值计算模型3.1Morison公式计算模型在式(2)中,墩柱绝对位移ut=us+ug,us为墩柱相对位移,ug为地面运动位移,并令Μw=CmρπD24‚,将式(2)作为外力施加在墩柱基本运动方程(9)的左侧,整理后可得式(10)。Μüs+Cu˙s+Κus=-Μüg(t)(9)(Μ+Μw)üs+Cu˙s+Κus=-üg(t)(Μ+Μw)+F(t)(10)式中,F(t)=-ΜD(u˙g+u˙s)|u˙g+u˙s|,作为外力施加在墩柱上。据此建立桥墩离散后的计算模型[见图1(b)]。图1中mu、cu、ku分别模拟墩柱上部结构的质量、阻尼与刚度。动力计算中桥墩节点速度u˙si未知,在每一载荷步中需先设定i节点时刻t的速度初始值等于上一时刻t-Δt的速度值进行迭代计算,即u˙si(t)=u˙si(t-Δt),直至2次计算的节点速度之差在误差控制范围之内,本文对此问题编程求解。3.2辐射波浪法计算模型将式(3)、式(4)所求解的刚体运动所致动水压力,以及弹性振动所致动水压力施加在墩柱基本运动方程式(9)左侧,并考虑到up(1)=ug,up(2)=us,整理后可得:(Μ+Μ(2)out)üs+(C+C(2)out)u˙s+Κus=-(Μ+Μ(2)out)üg+F(t)(11)式中,F(t)=-(Μ(1)out-Μ(2)out)üg-C(1)outu˙g作为外力施加在墩柱上。则可以建立桥墩离散后的计算模型[见图1(c)]。3.3流体声单元法计算模型及改进模型在ANSYS软件中流体声单元Fluid30、Fluid130能够模拟流体介质,并能与其它单元配合模拟流体和固体的相互作用。Fluid30单元有8个节点,每个节点有4个自由度:X、Y、Z方向的位移以及压力P。Fluid130单元用来对Fluid30单元建立的流体域进行封装,模拟无限流体域对声波的吸收,使到达模型边界的压力波最大程度地被吸收,最小程度地返回到流体域中。Fluid130要求建立的流体域为球形流体域,这要求墩柱完全包含在流体域中,并且球形流体域的半径R=0.2c/f(c为流体中的声速,f为压力波的主频)。通常情况下球的半径较大,将引起建模、网格划分的困难,导致计算效率低甚至无法计算。因此,需要对这种有限元建模方式改进。Fluid30单元中mu用来指定边界导纳(吸收系数),通常由试验确定,取值为0～1,mu=0表示没有波被吸收,mu=1表示波被完全吸收。根据桥墩截面形状,将流体域取为圆柱体形,令最外层Fluid30单元的吸收系数为1,模拟Fluid130单元对声波的吸收。基于这种流体域边界条件简化方法,得到改进的墩-水耦合有限元计算模型[见图1(d)]。4基于深水桥的动态特性和动态响应分析以墩高及入水深度均为60m、桥墩的直径1～20m的圆形实心墩作为实例进行对比分析。4.1种方法计算结果的比较设桥墩在空气中的第i阶自振频率为ωai,在水中的第i阶自振频率为ωwi,定义decωi=100%×(ωai-ωwi)/ωai为桥墩第i阶自振频率的降低率。由结果可知:3种方法计算的1阶自振频率降低率(见图2)与2阶自振频率降低率相似;Morison公式计算的降低率恒定在10.4%左右,不随桥墩直径发生变化;辐射波浪法和流体声单元法计算的频率降低率随桥墩直径的增大而降低并趋于一致,稳定在6%左右。总体上看,Morison公式计算结果偏大,辐射波浪法略大于流体声单元法。对直径较大的桥墩,Morison公式高估了动水压力的作用。4.2下挖深度和高建立高．墩高60m、直径10m的桥墩节点位移地震荷载可视为不同频率和振幅的简谐荷载的叠加。对上述桥墩在X方向施加简谐荷载üg=5sin(2πωt),ω=3Hz,取步长0.01s,计算时长TE=2.0s。由于桥墩顶部位移X方向分量(UX)相对较大;桥墩底部截面受力复杂,在地震中容易破坏。因此对墩顶UX、墩底反力X方向的分量(FX)及绕Y轴的弯矩(MY)进行研究。设桥墩第i截面节点t时刻的位移为uit,整个荷载作用期间峰值位移uitmax=max(uit),t=0,…,TE。当不考虑水体时节点的峰值位移为ui-airtmax,考虑水体时峰值位移为ui-wattmax,水体引起的桥墩节点的峰值位移增量ui-incretmax=100%×(ui-wattmax-ui-airtmax)/ui-airtmax。与此类似,桥墩第i截面的峰值反力增量Fi-incretmax=100%×(Fi-wattmax-Fi-airtmax)/Fi-airtmax;峰值弯矩增量Mi-incretmax=100%×(Mi-wattmax-Mi-airtmax)/Mi-airtmax。动力响应计算结果见图3～5。从图3可见:对于墩高60m、直径10m的圆形桥墩,考虑动水压力时,桥墩不同高度处各节点的位移峰值都比不考虑动水压力影响时大;考虑动水压力时,3种计算模式得到的各节点位移相差不大。分析发现动水压力对该桥墩不同高度处FX峰值、MY峰值的影响与UX峰值相似。从图4可见,在桥墩直径较小时(10m以下),Morison公式与辐射波浪法、流体声单元法墩顶UX增量差值在5%之内,在桥墩直径10m以上时,差值不断增大;辐射波浪法与流体声单元法UX增量差值不大,在2%左右。从图5可见,Morison公式与辐射波浪法的FX增量比较接近,桥墩直径增大时,Morison公式仍高估动水作用;流体声单元法计算的结果比前两者都小,与辐射波浪法相比,差值在7%左右。分析还发现,在桥墩直径小于10m时,Morison公式与辐射波浪法、流体声单元法MY增量差值在5%之内,随桥墩直径增大,Morison公式计算结果越趋于偏大;辐射波浪法计算结果介于Morison公式和流体声单元法之间,在直径较大时,与流体声单元法的差值约为3%,与Morison公式的差值约为4%。总体上看:各模式计算的动水压力对桥墩动力反应的影响随桥墩直径的变化有相似的变化规律;Morison公式计算结果较大,高估了动水压力的作用,流体声单元法较小,辐射波浪法介于两者之间。4.3mopson公式及流体声单元法Morison公式推导过程清晰明确,动水压力表达式简洁易于计算。考虑动水压力的墩柱模型建模简单,虽然动力响应计算中每一载荷步因节点速度未知,需要迭代计算,但合理的计算方案能将迭代次数控制在数步之内,所以计算效率仍较高。Morison公式忽略了墩柱本身对于流场的影响,墩柱直径较小时,对流场影响不大;墩柱直径较大时,对流场的影响不可忽略,此时Morison公式将高估动水压力对墩柱的影响。通过选择不同截面的阻力系数和惯性力系数,Morison公式能用于圆形、正方形(或带倒角的正方形)、三角形(或带倒角的三角形)截面的实心墩柱计算,但不能计算空心墩柱内域水产生的动水压力。阻力系数和惯性力系数根据试验或经验确定,这2个系数的准确度对计算结果的精度产生直接影响。辐射波浪法在求半解析解时力学概念清晰,推理严谨,求解弹性振动引起的动水压力时将墩柱离散为节点,形成半解析半数值解法。动水压力表达式繁杂,需要借助数学软件编程求解。考虑动水压力的墩柱模型建模简单,同Morison公式建立的墩柱模型相似,并且不用迭代计算,计算效率高。该方法精度较好,目前仅适用于简单截面,如圆形实心墩、圆形空心墩及矩形空心墩。流体声单元法是采用位移(结构)～压力(流体)格式的有限元方法。Fluid30单元要求与之耦合的墩柱为三维实体单元,且流体域通常较大,因此墩-水耦合有限元模型包含的单元和节点数量大,每一载荷步中需在流场内部、结构场内部及流场和结构场之间迭代计算直至收敛,计算效率较低。合理的流体域范围、网格划分能平衡计算精度和计算效率的需求。该方法适用于各类截面的实心、空心墩柱。5mopson公式的适用性(1)通过对流体声单元法中流体域边界条件的改进,可简化流体域建模,提高该方法的计算效率。(2)Morison公式、辐