热-固耦合（精选7篇）
热-固耦合 篇1
1 柴油机气缸盖模型的建立
为了提高柴油机功率,使其缸内平均有效压力不断地提高,但是与此同时气缸盖所承受的机械负荷和热负荷也变得越来越严重,使气缸盖的寿命大大降低。为此本文对气缸盖进行了热-固耦合分析。
1.1 柴油机气缸盖的边界条件确定
在对气缸盖进行有限元分析时不仅要确定机械负荷作用时的边界条件,还需确定热边界条件。
(1)约束边界条件:
在施加位移约束时,由于气缸沿中心面基本对称,故在气缸体、气缸盖和气缸垫的对称面上施加对称约束,在气缸体底面施加全约束。
(2)载荷边界条件:
在预紧工况下,采用杆单元模拟螺栓、刚性单元模拟螺栓连接,通过在杆单元节点上施加负温度让杆单元收缩,模拟螺栓预紧力。
(3)接触边界条件:
建立缸盖与气缸垫、气缸垫与气缸体顶面的接触,计算程序根据接触状态自动在接触面上建立接触单元,传递作用力,完成接触模拟。
在稳定工况下,气缸盖燃气侧表面0.5 mm深处温度波动值不超过±1.2 ℃。因此,在对缸盖进行热强度分析时,可将其作为稳态进行分析。
1.2 柴油机气缸盖模型的建立
对气缸盖模型进行适当的简化,同时为了在较小的有限元计算规模基础上提高计算精度,必须对几何模型进行分区处理,在不同的分区采用不同的单元类型及单元尺寸。
缸盖的实体模型如图1所示,缸盖有限元网格模型如图2所示。
2 柴油机缸盖温度场的计算
图3为缸盖温度场云图。由图3可以看出高温区出现在底板火力面上,气门座孔与喷油器座孔之间和进排气门之间的鼻梁区域温度最高,达到327 ℃。
3 柴油机气缸盖应力计算
本文仿真计算了柴油机在机械负荷下缸盖的机械应力和热负荷下缸盖的热-机械耦合应力,同时对气缸盖关心区域进行采点比较分析,区域的选取主要考虑热负荷对整个气缸盖负荷的影响。图4为所选取的11个区域,这些区域主要是气缸盖火力面、喷油器轴线以及气缸盖底板水腔面等。
3.1 机械负荷的计算
预紧工况下缸盖只受到螺栓预紧力作用。经过有限元计算后,图5显示了预紧工况下整个缸盖等效应力的分布情况。
图6为预紧工况下考察区域应力曲线图。由图6可以看出,在4、7、8号区域等效应力最大。
3.2 热负荷计算
图7为热负荷等效应力云图,图8为热负荷热应力曲线图。由图7及图8可以看出由热负荷引起的热应力主要出现在喷油器座及进气门座靠近水腔侧的3、6、7号区域。热应力是由于温差引起的,火力面温度虽然最高,但其温差并不大,而3、6、7号区域温度虽不是最高,但由于和水腔紧贴,出现了较大的温差,因此成为最大热应力区域。
3.3 热-固耦合计算
在边界条件的加载中,不仅要考虑约束边界条件和接触边界条件,还要考虑热边界条件。图9显示了在耦合预紧工况下缸盖等效应力分布情况,图10为耦合预紧工况下考察区域应力曲线。
4 结论
通过以上计算分析可以发现气缸盖各考察区域的应力变化趋势不同。在预紧工况下气缸盖受螺栓预紧力,最大等效应力出现在气缸盖上的4、7、8号区域;当考虑热负荷作用时,气缸盖受到热-固耦合共同作用,气缸盖内部拉压情况变得更加复杂,各部分应力不同程度地改变,而且这种改变并不是简单的线性叠加,呈现明显的非线性。通过仿真计算表明高应力区域在不同负荷与工况下是不同的,因此在设计气缸盖时,不能简单考虑某一工况下单一负荷的作用,必须全面考虑各种不同负荷的耦合作用,这样才能使气缸盖在使用过程中不易损坏。
摘要：在建立气缸盖实体模型的基础上,利用相关边界条件计算了气缸盖上各区域的换热系数,并对气缸盖进行了温度场仿真计算;分析了在多种工况下气缸盖的应力,确定了结构上的薄弱区域。认为当气缸盖受到热-固耦合共同作用时,气缸盖内部拉压情况变得更加复杂,各部分应力呈现明显的非线性。
关键词：气缸盖,温度场,热-固耦合,应力
参考文献
[1]沈红斌.495汽油机缸盖的有限元结构强度计算分析[D].天津:天津大学,2003:21-22.
[2]贾明,解茂昭.均质压燃柴油机燃烧与排放的多区模型模拟[J].燃烧科学与技术,2005,37(5):263-264.
[3]张骜.柴油机气缸盖热负荷与排气温度的相关性研究[J].武汉船舶职业技术学院学报,2003(4):14-15.
[4]代秀红.Z6110柴油机缸盖结构强度[D].大连:大连理工大学,2003:32-34.
热-固耦合 篇2
液压滑阀应用领域极为广泛, 在使用过程中经常会出现磨损, 在中高压系统使用过程中, 阀芯的磨损情况更为严重, 甚至出现阀芯卡紧的现象, 其原因是当环境温度变化时, 液压滑阀阀芯、阀套之间的滑动副间隙发生变化[1]。液压滑阀卡紧故障的诊断主要通过液压系统的压力、流量检测来实现, 需要依赖专门的传感器件[2]。在液压阀的设计过程中, 阀芯和阀套理论上是完全同心的, 不管它在多大的压力下工作, 移动阀心所需的力只需克服黏性摩擦力, 数值上应该是很小的, 但在实际中并非如此, 特别是在中高压系统中, 阀开口处的节流作用产生温升, 造成阀内温度分布不均匀, 阀内局部受热, 阀芯阀套发生热膨胀, 使阀套与阀心间的配合间隙减小, 出现阀芯阀套间的磨损、卡死、卡滞、压力损失增加或泄漏加剧等现象。当发生卡死现象时, 系统便不能正常工作[3]。
很多学者对液压滑阀进行了研究, 也取得了相应的成果[4-5], 主要手段是计算流体力学 (CFD) 和实验, 或是两者结合, 但研究内容多涉及流场结构、压损、液动力以及气穴噪声, 也有部分文献考虑了液压黏性热的影响[6]。本文主要利用流-固-热耦合来研究液压滑阀工作过程中液压阀内部流场的变化情况以及液压阀表面温度场和应力场的分布情况, 为减少滑阀设计中阀芯卡死现象奠定理论基础。
1 理论分析与建模
滑阀主要由阀体和阀心等组成, 一般阀体上均有压力平衡槽。滑阀内有多个阀腔, 由于相似性, 用其中的一个阀腔完全能说明阀内的流体流动状态。阀腔内的流体流动和传热的现象十分复杂, 但都离不开3个物理规律的支配, 即质量守恒、动量守恒和能量守恒[7?8]。
质量守恒方程为
动量守恒方程为
能量守恒方程为
式中, 为散度;U为速度矢量;p为流体压力;μ为流体动力黏度;T为流体温度;cp为流体的质量定压热容;λ为热导率;F为作用在流体上的质量力;q为流体所吸收的热量;Φ为能量耗散函数;ρ为流体密度;ε1为流体的变形张量, 代表流体克服黏性所消耗的机械能, 它将不可逆转地转化为热而耗散掉。
耦合分析选取了温度场的第三类边界条件和耦合特有的边界条件。第三类边界条件为固体与流体因温度差而发生了对流换热, 此时固体表面的热流密度与温度差成正比, 即
式中, n为换热表面的外法线;h为传热系数;Tw为边界面温度;Tf为流体温度。
由于流体和固体遵循不同的控制方程, 所以固体和流体交换壁面上温度和热流密度必须满足连续性边界条件[9], 即
式中, Ts为固体温度;λf、λs分别为流体和固体的热导率。
其中, , 需要通过耦合迭代求解耦合面的温度和热流密度。
初始条件是过程开始时, 物体整个区域中所具有的温度
本文计算中采用的流体与固体的材料属性参数分别如表1和表2所示。
进行稳态计算分析时, 本文对流体作了如下假设:流体为不可压缩恒定流动的牛顿流体, 在模型中的流动状态主要为紊流, 采用k-ε紊流模型, 流体计算模型如图1a所示, 阀口局部细化如图1b所示。
2 结果分析
2.1 速度场分析
为了得到液压油黏性加热温升的大小和位置变化趋势, 本文采用ICEM和FLUENT软件计算了流场, 得到图2所示液压油速度场的分布情况。液压油从过流面A0流出节流槽口时, 由于过流面面积减小, 液压油流动速度提高, 此区域为B1区;经过过流面A1时开始形成高速射流, 在过流面Ak附近处达到最大流速, 沿着液压油流动方向形成一个明显的液压油高速射流区域, 即B2区;随后液压油冲入阀腔内B3区, 此时, 随着过流面面积的增大, 液压油流动速度逐步减小。
图3、图4所示分别为液压滑阀槽口深度和宽度对液压滑阀流场的影响情况。图3中射流角 (流体射流方向与水平之间的夹角) 是随着槽口深度的增大而不断增大的, 最高射流区域会随着槽口深度的增大而减小, 并逐步往槽口出流口处移动。最高射流区域的大小会随着槽口宽度的增大而增大, 这主要是因为过流面积增大了。图4则描述了最高流速与槽口深度、宽度和阀芯开口度之间的关系。可以看出, 最高流速是随着阀芯开口度的增大而先减小后增大的, 主要和过流面积有关, 即上文的阀口迁移现象造成的。
2.2 阀芯表面温度场分析
本文采用ANSYS WORKBENCH软件进行阀芯表面温度场分布计算, 油液初始温度设为50℃并忽略温度变化对液体黏度的影响。
图5为阀芯温度场分布示意图, 可以看出, 阀芯中间段及槽口附近是明显的高温区, 阀芯两端是低温区, 右边的低温区明显比左边的低温区大, 这是因为右边是油液上流, 温度较低;左边是油液下流, 油液经过槽口, 速度增大, 黏性加热显著, 产生大量热量。热量一部分随着油液流向出口, 一部分则通过槽口表面传向阀芯。同时, 在槽口内部存在一个相对低温区, 如图5所示。
图6描述了槽口宽度和深度对阀芯表面温度场的影响, 可以看出, 阀芯表面的高温区域主要集中在槽口附近, 并且是随着开口度的增大逐渐往槽口入口处移动的, 这是随高速射流主要区域的变化而产生的。最高温度区域是随着槽口宽度的增大而增大, 随着槽口深度的增大而从槽口中心逐渐往出口方向移动的, 这同样是和高速射流区域的变化相关的。
阀芯表面的最高温度数值与槽口深度、宽度以及阀芯开口度关系如图7、图8所示, 可以看出, 阀芯表面的最高温度随着开口度的增大而先增大后减小, 随槽口宽度增大而减小。
.H=1.0mm, 阀芯2.H=1.5mm, 阀芯3.H=2.0mm, 阀芯.H=1.0mm, 阀套5.H=1.5mm, 阀套6.H=2.0mm, 阀套
1.B=1.0mm, 阀芯2.B=1.5mm, 阀芯3.B=2.0mm, 阀芯4.B=1.0mm, 阀套5.B=1.5mm, 阀套6.B=2.0mm, 阀套
将阀芯表面的温度场耦合到固体结构分析中, 计算出阀芯和阀套的最大变形量, 如图9、图10所示。
由上文分析可知, 阀芯和阀套受热膨胀的径向变形量不同, 由图9可知, 阀芯变形量受到槽口深度的影响较大, 在其他条件不变的情况下, 液压阀不同槽口深度下槽口越深, 发生摩擦卡紧的机率越大;由图10可知, 阀套变形量受到槽口宽度的影响较大。液压阀不同槽口宽度下槽口越宽, 阀套径向变形量越大, 发生摩擦卡紧的机率越大。
1.H=1.0mm，阀芯2.H=1.5mm，阀芯3.H=2.0mm，阀芯4.H=1.0mm，阀套5.H=1.5mm，阀套6.H=2.0mm，阀套
1.B=1.0mm, 阀芯2.B=1.5mm, 阀芯3.B=2.0mm, 阀芯4.B=1.0mm, 阀套5.B=1.5mm, 阀套6.B=2.0mm, 阀套
3 结语
综上所述, 液压滑阀在节流口处出现高温膨胀使得阀芯磨损甚至卡紧的趋势与试验结果基本吻合。液压滑阀在节流口处出现高速射流, 引起液压油温度升高, 节流口处局部温升过高, 并且不同的结构有不同的结果, 所以, 在液压滑阀的工程设计和应用时, 应综合考虑液压黏性温升的影响, 在满足工况需要的情况下, 尽量选择黏度系数小的液压油, 并且在相同节流面积下, 选择节流槽口宽度相对大的液压阀。总之, 应综合考虑液压黏性温升对液压阀的影响, 使液压阀在满足工况的同时, 还能满足在高温环境下均正常工作的要求。
摘要：针对液压滑阀在中高压系统的使用过程中, 黏性热效应使得油流温升显著, 阀芯受热膨胀而出现磨损甚至卡紧, 导致液压滑阀失效的现象, 建立了液压滑阀的计算流体力学三维模型, 从不同开口度、不同槽口深度和宽度三方面进行流场和温度场的耦合分析, 得到液压滑阀的最高流速和最高温度大小和区域的分布情况, 以及阀芯和阀套的径向变形量, 为液压滑阀卡紧机理分析提供了强有力的支撑和参考。
关键词：滑阀,温度场,径向变形,流-固-热耦合
参考文献
[1]訚耀保, 肖其新, 闫世敏.温度对电液伺服阀的影响分析[J].流体传动与控制, 2008 (6) :23-25.Yin Yaobao, Xiao Qixin, Yan Shimin.The Influence of Temperature on the Electro-hydraulic Servo Valve[J].Fluid Drive and Control, 2008 (6) :23-25.
[2]李伟, 李伟波, 吴根茂, 等.滑阀机构液压卡紧无传感器诊断方法研究[J].中国机械工程, 1999, 10 (4) :436-438.Li Wei, Li Weibo, Wu Genmao, et al.The Research of the Slide Valve Hydraulic Clamping Sensorless Method[J].China Mechanical Engineering, 1999, 10 (4) :436-438.
[3]刘晓红, 柯坚, 刘桓龙.液压滑阀径向间隙温度场的CFD研究[J].机械工程学报, 2006, 42 (S) :231-234.Liu Xiaohong, Ke Jian, Liu Huanlong.The Radial Clearance CFD Study of Temperature Field of Hydraulic Slide Valve[J].Journal of Mechanical Engineering, 2006, 42 (S) :231-234.
[4]高殿荣, 王益群.液压锥阀流场的有限元法解析[J].机床与液压, 2000 (2) :12-14.Gao Dianrong, Wang Yiqun.The FEM Analysis of the Hydraulic Cone Valve Flow[J].Hydromechatronics Engineering, 2000 (2) :12-14.
[5]王林翔, 陈鹰, 路甬祥.液压阀道内三维流体流动的数值分析[J].中国机械工程, 1999, 10 (2) :127-129.Wang Linxiang, Chen Yin, Lu Yongxiang.The Three Dimensional Numerical Analysis of the Hydraulic Valve Passage for the Fluid Flow[J].China Mechanical Engineering, 1999, 10 (2) :127-129.
[6]Eckert E R G, Faghri M.Viscous Heating of High Prandtl Number Fluids with Temperature-dependent Viscosity[J].Heat Mass Transfer, 1986, 29 (8) :1177-1183.
[7]Ju Y.Thermal Conductionand Viscous Heating in Microscale Coquette Flows[J].Heat Transfer, 2000, 122:817-818.
[8]Deng Jian, Shao Xueming, Fu Xin, et al.Evaluation of the Viscous Heating Induced Jam Fault of Valve Spool by Fluid-structure Coupled Simulations[J].Energy Conversion and Management, 2009, 50 (4) :947-954.
基于岩体地下流固耦合理论的研究 篇3
关键词：裂隙岩体,三场耦合,流固耦合,数学模型
0前言
对裂隙岩体而言, 其工程特性决定于裂隙岩体所处的地质环境, 即裂隙岩体存在的地下水渗流场、应力场和温度场环境等等, 这些环境因素相互作用和影响, 使得裂隙岩体时时处于这些因素构成的动态平衡体系中。因而研究裂隙岩体的工程特性, 是一个涉及到诸如岩体力学、土力学、工程地质学、地下水力学、工程热力学等多学科的问题。与此同时, 一方面, 随着科学技术的发展及人类认识和利用自然能力的提高, 上述各学科自身都在不同程度地丰富和发展;另一方面, 针对岩土工程的实际需要, 上述各学科之间的相互交叉和渗透又日益突出。基于这一情况, 开展裂隙岩体地下水渗流场、应力场与温度场三场之间耦合体系作用机理的研究, 在促使地下水力学、岩体力学、工程热力学及工程地质学等各学科之间相互交叉和渗透方面, 具有重要的理论意义;从工程实践的角度出发, 大型土建工程的修筑, 使工程设计和建设者都面临许多新问题, 同时土建工程体自身会产生诸多不良反应 (包括自身稳定性、持久性和对环境产生的不良影响等) , 这时用传统的单个学科来解决以上新问题已显得无能为力, 即使勉为其难, 效果往往也不十分理想。这就需要相关多个学科的相互结合、补充和渗透, 以对问题的解决起到事半功倍的效果。进行裂隙岩体地下水渗流场、应力场与温度场三场耦合作用的研究, 是试图建立一种基于多学科综合应用基础上的广义理论体系, 并应用于工程体修筑中不良地质灾害 (包括涌水、岩爆和岩体失稳等) 的防治研究, 因而这项研究又具有重要的工程实际价值。
1三场耦合作用研究的模式
裂隙岩体地下水渗流场、应力场与温度场耦合作用是一个相对复杂的研究问题, 主要表现在裂隙岩体所处地质环境的各个组成部分, 即地下水渗流场、应力场与温度场自身都随时间、空间发生变化。与此同时, 以上各个组成部分之间的耦合作用处于一种复杂动态变化过程之中。为了研究的可行性, 笔者从简化场耦合作用具体过程的角度出发, 基于裂隙岩体对应场性能等效的原则, 通过裂隙岩体结构特征的量化研究, 对由裂隙结构面网络离散化所引起裂隙岩体的渗透性能、力学性能及热物理性能的非连续性特点, 分别进行了等效连续化处理。在此基础上, 结合研究工作的目的, 本文提出了基于深层地下水资源量评价预测目的基础上的裂隙岩体地下水渗流、应力、温度三场之间“非完全”耦合作用研究模式, 如图1所示。同时, 针对大型地下工程中所涉及地质灾害的成灾预测预报及防治需要, 给出了裂隙岩体地下水渗流场、应力场与温度场“完全”耦合作用研究模式见图1。
2孔隙介质流固耦合
在连续介质力学框架内建立了Eu ler型多相流体运移和变形孔隙介质耦合的理论模型。在这些研究中, 孔隙介质都被当作具有同一渗透率单一类型孔隙连续分布的介质, 即具有单孔隙度/单渗透率的介质。考虑有效应力影响及孔隙压力的一般应力应变关系可表示为:
undefined
式中, a为压力比系数, 可取α=K/H, K为固体介质的有效体积模量, H为比奥常数, 不计体力和惯性力时的方程为:
σij, j=0
应变-位移关系为:
undefined
可得固体变心控制方程为:
Gui, jj+ (λ+Guk, ki) +αP, i=0
式中, G是剪切模量;λ是拉梅常数。
对于流体相, 达西速度为:
undefined
式中, k是渗透率, μ是流动粘度。
3裂隙岩体流固耦合研究
自从1954年12月法国的M alpaset拱坝失事以后, 裂隙岩体渗流问题日益受到人们的重视, 许多人通过试验及理论分析研究了单一裂隙渗流与应力的关系。如Lou is根据一些钻孔压力试验结果给出以下经验公式有:
kf=k0e-ασ
式中, k0为初始渗透系数;σ为法向应力;α为参数。
Jones提出碳酸盐类岩石裂隙渗透系数的经验公式为:
(kfA) =Q0 (pc-pf) -n
式中, A为过水面积; pc为总压力;pf为内部孔隙压力; n为常数。
Liu J根据裂隙岩体的试验结果, 得出如下公式:
undefined
式中, R m 为岩体分类指数;Δεi为应变;φundefined为初始孔隙度; i= 1, 2, 3。
4裂隙岩体渗流与应力耦合分析的数学模型
目前, 对裂隙岩体渗流的数学处理主要有以下几种模型:
(1) 连续介质数学模型。该模型的优点是可以直接应用较成熟的孔隙介质饱和渗流分析方法来求解裂隙岩体饱和渗流问题。不足之处在于把裂隙网络等效为连续介质, 不能很好地刻画出裂隙的特殊导水作用, 故拟真性不好。另外, 该模型的适用范围很受限制, 不是所有岩体都可等效为连续介质。
(2) 离散裂隙网络模型。该模型的优点是较好地描述了裂隙岩体的非均匀各向异性, 故当岩块很致密, 确实可忽略其渗透性时, 具有拟真性好、精度高的优点。其不足之处在于研究域内全部裂隙的几何参数很难获得。另外, 模型的计算量较大。
(3) 双重介质模型。该模型认为流体不仅通过岩体的裂隙面流动, 而且还通过结构体流动。显然, 它能较为全面地反映裂隙岩体的渗流特征。据岩体不同的渗透率又可分为双孔隙度/单渗透率模型和双孔隙度/双渗透率模型。
(4) 离散介质-连续介质耦合模型。该模型将岩体中的裂隙根据其迹长和开度等划分为主干裂隙和次要裂隙。该模型的耦合条件是:裂隙流体压力作为连续介质域的压力边界;离散介质域与连续介质域间的流体交换作为离散介质域的流量边界。该模型综合了连续介质模型和离散介质模型的优点, 即能反映裂隙特殊的导流作用, 又能体现岩块的贮藏作用, 故很好地解决了精度与可操作性之间的矛盾。但由于描述连续介质域的流体运动和描述离散介质域流体运动方程不同, 给数学处理带来困难。另外, 该模型同样存在双重介质模型中的流量交换的问题。
5结论及意义
(1) 对孔隙介质的流固耦合机理目前主要集中在宏观研究上, 对微观机理研究较少, 故有必要加强。
(2) 目前, 国内外对裂隙单裂隙渗流规律主要精力放在如何修正立方体定律上, 但此定律中一些参数难以确定, 为便于工程应用, 应研究渗透性与岩体的分类指标的相关性上, 这样也可充分利用日益先进的测井技术成果。
(3) 许多流固耦合问题与温度场有关, 目前国外许多学者在这方面做了不少工作, 但国内这方面研究工作开展得还较少, 应加强。
(4) 尽管人们已经认识到孔隙、裂隙的形态分布规律均具有一定的分形特征, 但在已有的研究中, 真正利用分形几何等非线性理论建模的仍较少, 应加强分形等非线性科学在流固耦合问题中应用的研究。
(5) 混合物理论在孔隙固体流固耦合问题中已得到很好地利用, 今后应加强利用混合物理论的基本观点和方法来研究双重介质渗流, 热-液-力三场耦合等问题。
(6) 较符合实际的双重介质模型和离散介质-连续介质模型的拟真性, 在一定程度上取决于流体交换项确定的准确性, 但现有的交换项公式精度仍不高, 如何确定一个简单而又精确的流体交换项公式有待于进一步研究。
参考文献
[1]刘建军, 薛强.岩土热-流-固耦合理论及在采矿工程中的应用[J].武汉工业学院学报, 2004, 23 (3) .
[2]Desai C S.Introduction to the Finite Element Method[M].New York:Van Nostrand Reinhold Co, 1972.
[3]仵彦卿, 张倬元.岩体水力学导论[M].成都:西南交通大学出版社, 1994.
[4]黄涛, 陈一立.工程岩体地下水渗流-应力-温度耦合作用数学模型研究[J].西南交通大学学报, 1999, 34 (1) :11-15.
[5]刘建军, 刘先贵, 胡雅衽.低渗透储层流固耦合渗流规律研究[J].岩石力学与工程学报, 2002, 21 (1) :46-51.
[6]冉启全, 顾小芸.油藏渗流与应力耦合分析[J].岩土工程学报, 1998, 20 (2) :69-73.
[7]赵阳升.矿山岩石流体力学[M].北京:煤炭工业出版社, 1994.
[8]柴军瑞.混凝土坝渗流场与稳定温度场耦合分析的数学模型[J].水力发电学报, 2000, (1) :27-35.
[9]南建林.混凝土的温度一应力耦合本构关系[J].清华大学学报, 1997, 37 (6) :87-90.
热-固耦合 篇4
基坑工程的设计，首先要确定围护墙上的水土压力分布，渗流条件下的基坑围护墙上水土压力计算一直是深受重视的问题，引起了广泛的讨论。魏汝龙[1,2]、沈珠江[3]和李广信[4]分析了水土压力合算与水土压力分算的条件，认为水土合算方法与土力学基本原理相背，计算中低估了主动状态的水压力作用而高估了被动状态的水压力作用，偏于不安全;水土分算方法概念清楚，但在实际应用中有效应力指标难以确定，也无法考虑土体在不排水剪切时产生的超静孔压的影响。另外，地下水位较高、水量丰富的地区，基坑开挖坑内外水位差较大，基坑围护墙的水土压力分布与稳定性不可避免地受到地下水渗流的影响。Tanaka[5]分析了二维渗流条件下带支撑围护结构的涌土破坏及围护稳定性降低问题;Benmebarek[6]分析了深基坑围护墙周围地下水绕流引起的土体渗透破坏问题，指出了土体强度参数和围护墙摩擦等对渗透的影响;Schfer[7]分析了地下连续墙施工过程的水土压力分布模式，对不同基坑排水情况对土体参数和库伦土压力的影响，指出了在围护结构两侧存在较大水位差时考虑渗流效应的重要性及其计算方法。
一些学者对不同的计算方法做了较全面的对比，如陈愈炯，温彦锋[8]分析了基坑周边土体内孔隙压力随工程进展而变化的情况，叙述了不同工况下的水土压力分算和合算方法;王钊[9]，李广信[10]等总结分析了在静水压力、稳定渗流和超静水压力作用下，挡土结构上水土压力的计算方法，通过实例说明了采取水土分算及考虑地下水的渗流作用和超静孔压力作用对挡土结构上总压力计算的影响;董诚[11]等分析了传统水土压力计算方法的土层实用性，将传统计算结果与PLAXIS的结果对比，指出PLAXIS计算的合理性，但是分析仅限于简单土层还难以考虑基坑工程的施工工况等因素。
基坑内外水头差使坑外的地下水绕过围护墙向坑内渗流，导致土体中孔隙水压力和有效应力的改变以及基坑周围土体变形和渗透特性的改变，而土体固结变形和渗透特性的改变又通过渗透动水压力影响土体的应力、应变状态，此即为基坑周围土体中的渗流固结耦合效应。由基坑围护墙与土体的作用效应可知，土体的渗流边界改变及流—固耦合效应必然影响围护墙上的水土压力分布，本文针对地下水绕过围护墙渗流情况，分析了传统的水土压力分算、合算及考虑土体渗流—固结变形方法计算土压力的区别，并利用实测数据进行对比，期望改进软粘土和粉土地层中基坑围护墙上水土压力的计算方法。
1 不同计算方法简介
当前工程界对于基坑工程中的水土压力计算方法的讨论，多集中在透水性弱的粉土和粘性土中围护墙上的水土压力分布问题，概括起来，主要有“水土压力分算法”、“水土压力合算法”和根据基坑周围流网，考虑渗透力影响的水—土压力分算法，下面对它们作简单介绍。
1.1 水土压力分算
水土压力分算法从土的有效应力理论出发，采用有效强度指标c′和φ′，将围护墙上的水压力与土压力分开计算，其公式为:
围护墙上的水、土压力的总和计算公式为:
式中:p′a,p分别为主动、被动土压力;k′a,p分别为主动、被动土压力系数;σ′z为竖向有效应力，不考虑渗流效应，σ′z等于γ′z;u0为静水压力。
对于砂土等透水性强的土层，采用水土压力分算的方法计算没有异议，其假设符合土力学的基本原理。
计算粘性土特别是软粘土的土压力时，有学者建议使用不排水强度指标，如Peck[12]等建议，对于粘土和粉土采用不排水剪切指标计算主动土压力和被动土压力，公式分别为:
式中:在水下γ为饱和容重;c为不排水抗剪强度。
通过分析可知，不排水情况下φu=0，水的侧压力系数、土的主动和被动土压力系数都是1.0，水土分算和合算的公式表达形式相同。计算的主动土压力偏大而被动土压力偏小，在基坑开挖较深而c值不大的情况下，由于主动土压力与被动土压力的增加斜率相等，基坑围护墙难以平衡，设计较为保守。
1.2 水土压力合算
水土压力合算采用总应力强度指标计算土压力，把围护墙上的水土压力合并计算。采用的公式为:
式中:γ的值在水下取饱和容重;ka,p分别为与总应力指标对应的主动、被动土压力系数。
工程中，由于有效应力指标难以测定及基坑开挖产生的负超静孔压难以准确计算，工程师们采用三轴固结不排水强度或固结快剪强度指标将水压力与土压力合并计算，不考虑孔隙压力的变化与作用，但其在理论上存在缺陷。
1.3 考虑渗流的水土压力计算
龚晓南[13]等阐述了渗流对土压力的影响，指出了不同土层条件计算的水土压力分布也不同，利用土层中流线的水流速度相等原理计算水头损失对土压力的影响，取如图1所示的流线计算坑外地下水绕过墙趾向坑内的渗流并计算水土压力的变化。
不考虑渗流水头损失，围护墙上的水压力分布线为AE与FD，假设围护墙不透水，取紧贴墙壁的绕流曲线计算水力坡降和水头损失，则围护墙上的水压力分布变为AE′与FD′，阴影三角形△AB′C为作用在围护墙上的净水压力分布。
对于多层土的情况，根据基坑围护墙附近流网特征，利用紧贴墙边的流线流速相等原理[14]，流经各土层的水力坡降i与渗透系数K成反比，由基坑内外的水头损失H1(为各层的总和)，可计算各土层中的水力坡降。公式为:
式中:v为渗流流速;K1,K2……Ki为各层土的渗透系数;i1,i2……ii为各层土的水力坡降;hi为各土层的厚度。
由于基坑外侧的地下水从墙底绕流进入坑内，主动侧地下水向下渗流引起土体自重应力增大，在计算过程中计入渗透力的作用，得出黏性土水压力与有效土压力的计算公式，围护墙上的总水平荷载为水压力与有效土压力的合力，计算公式见文献[13]。
2 不同水土压力计算方法的比较
基坑周围土体为分布层数较多的土层时，由于各层土体的强度和渗透特性不同，基坑围护墙上的水土压力分布变得较为复杂。为了对比水土合算、分算与考虑渗流情况的水土压力计算方法的不同，选取文献[14]的基坑实例进行分析，研究不同算法对水土压力分布的影响。
工程算例:中国平安金融大厦地处浦东小陆家嘴金融贸易中心区，大厦塔楼39层，两侧裙房4层，主体结构均设置3层地下室，基坑面积约为18000m2，基坑开挖深度主楼区域约为17.90m，裙楼区域约为16.90m，属超大型深基坑。基坑采用地下连续墙作为围护体，坑内设置三道钢筋混凝土围檩和支撑。场地地下水属潜水类型，地下水静止水位埋深约为0.50～2.19m，计算时采用1.0m。基坑土层参数如表1所示。
由前面的分析可知，水土压力分算时，对于水压力分布，有考虑渗流的影响与不考虑渗流的影响两种算法，两种方法计算的水压力分布如图2所示。
图2中，分布1的提出是因为淤泥质土等的渗透系数较上部土层较小，在淤泥质土层中水头损失很大，故文献[12]提出了简化的三角形分布;分布2是根据公式(5)计算而得，在淤泥质土层中，水头损失较大，围护墙主动侧上的水压力显著减小，在墙底端水压力小于按静水压力的计算值，在坑内被动侧的水压力大于按静水压力的计算值，渗流效应可以改善围护墙的水压力分布;分布3是将分布2的水压力在坑内与坑外相抵消，与分布1相比，考虑渗流效应计算的水土压力明显减小。
利用水土合算、水土分算及考虑渗流效应的分算方法计算围护墙上的总压力，开挖到坑底时的水土压力合力分布如图3所示。
从图中可以看出，水土压力合算，采用固结不排水或快剪强度指标，土体φ不等于零，将主动侧、被动侧水压力分别乘以土压力系数Ka、Kp，则计算所得的主动侧总水土压力变小，被动侧水土压力变大，使基坑围护墙的计算偏于不安全;不考虑渗流的水土压力分算法，将静水压力作用于围护墙上，计算的水土压力值在主动侧大于水土合算的结果，被动侧小于合算的结果;考虑渗流效应的水土压力算法，由于坑外侧的地下水从墙底绕流进入坑内，主动侧地下水向下渗流引起土体自重应力增大，被动侧地下水向上渗流引起土体自重应力减小，故计算的结果介于水土合算与水土分算之间，其分析方法可以反映基坑土层对土压力的影响。
以上分析的考虑渗流效应的简化计算方法，能较好地反映围护墙的水土压力分布，但是这种简化方法考虑的因素较少，难以获得较准确的计算结果。基坑开挖地下水渗流的过程，土体固结变形和围护墙位移等都对土压力产生影响，我们有必要分析土体渗流固结的耦合效应对土压力的影响，把土体流—固耦合效应下的土压力与传统计算方法及实测结果进行对比，发现他们的区别，找出符合实际的基坑围护墙水土压力分布的计算方法。
3 考虑流—固耦合的水土压力计算
基坑土体开挖可简化为卸荷条件下卸除侧向约束，同时存在地下水渗流固结的变形过程，土体中地下水渗流与固结变形都对基坑围护墙上土压力产生影响，利用岩土工程有限元软件可建立分析模型，设立特定水力边界条件，通过数值方法模拟流—固耦合效应来考虑渗流对土压力的影响，计算围护墙上的水土压力。
3.1 PLAXIS流—固耦合的土压力计算
PLAXIS中流—固耦合计算[15]，使用Biot理论的固结基本方程，渗流问题采用达西(Darcy)定理，并且假设土体骨架弹性变形和基于小应变理论。以增量形式表示的Biot固结本构方程为:
式中:为有效应力增量;为应变增量，且;M是材料刚度矩阵。
程序通过有限单元离散，建立刚度矩阵、耦合矩阵和荷载增量矢量的关系，利用连续方程引入渗流问题和边界条件，求得渗透力从而得到应力与应变，然后利用耦合反解孔隙率和渗透矩阵，再用新的渗透特性求解渗流和应力应变，经过反复迭代求得满足精度的解。
程序中采用水土分算的方法，通过输入地下水水头执行地下水渗流程序进行计算，利用单元应力点上的压力水头求得孔隙水压力，将围护墙与土体接触界面上的有效压力与孔隙水压力值相加，得到基坑围护墙上总的水土压力分布。
3.2 基坑模型的建立
PLAXIS建模计算，将文献[14]中的基坑围护参数输入程序，土体的强度指标、弹性模量、泊松比等参数根据原工程资料确定，将土层模型定义为不排水状态，将地下连续墙和界面定义为非多孔状态，建立基坑计算模型。其中，c′、φ′采用三轴固结不排水剪切试验测定的ccu、φcu指标由上海地区经验公式[16]估算而得，估算φ′的经验公式(纯数值统计关系)为:
式中:ccu为三轴固结不排水剪切粘聚力(kPa);φcu为三轴固结不排水剪切内摩擦角(°)。
3.3 基坑围护墙上的土压力
中国平安金融大厦基坑土层分四层开挖施工，第一层至第三层分别开挖到各道钢筋混凝土围檩及支撑底标高位置，第四层开挖到地下室底板标高位置。基坑设置了2个主监测断面，埋设土压力盒测量开挖过程中围护墙上的水土压力合力。现在，我们利用PLAXIS程序并考虑流—固耦合分析，将计算的围护墙上土压力值与传统的水土分算、水土合算及现场土压力实测值进行对比，研究各工况下不同计算方法的计算结果与实测水土总压力的差别。计算中，考虑到开挖时间较短，并且为了将有限元计算的水土总压力结果与常规方法的比较，围护墙上的土压力包含了超静孔压和负孔压力，为水土总压力。
基坑开挖到第一道支撑底(开挖深度1.8m)，围护墙上的水土压力实测值与计算值的分布如图4所示。
图4表明，在主动侧，采用水土合算的土压力计算值与实测值相差不大，而采用水土分算和PLAXIS计算的主动土压力明显大于实测值，在被动侧，水土分算和水土合算计算的朗肯被动土压力均远大于土压力实测值，而PLAXIS计算的土压力则与实测值很好地吻合。其原因可分析为，第一层土开挖深度较浅，围护墙的侧移很小，被动侧土体远没有达到被动状态，土压力没有发挥到被动土压力的极限水平，PLAXIS程序分析，基坑围护墙水平位移较小，所得土压力符合实际情况。上海市基坑工程设计规程，针对这种情况建议将水土分算的土压力乘以0.5倍的折减系数，即采用降低的被动土压力，文献[14]的计算值与实测值也较吻合。
基坑开挖到第二道支撑底(开挖深度7.5m)，围护墙上的水土压力实测值与计算值的分布如图5所示。
从图中可以看出，在主动侧，水土合算的土压力计算值仍接近于实测值，支撑标高附近PLAXIS计算的主动侧土压力有突出的趋势，而在较深的部位，PLAXS计算值略小于水土分算的结果，证明此时渗流固结计算使水压力减小;在被动侧，采用水土分算的朗肯被动土压力和PLAXIS计算值与实测值均较接近，并且PLAXIS计算值与实测值吻合的更好，而采用水土合算的朗肯被动土压力计算值仍较实测值大得多。
开挖到第三道支撑底(开挖深度13.15m)，围护墙上的水土压力实测值与计算值的分布如图6所示。
在该工况下，主动侧水土合算的土压力计算值更接近于实测值，PLAXIS的计算值与实测值均有变大的趋势，表明基坑支撑约束了围护墙的水平位移，使水土压力大于主动土压力;被动侧的土压力，水土分算的朗肯被动土压力和PLAXIS计算值仍与实测值较接近，表明开挖第三层土时，围护墙有一定的水平位移，被动侧的土压力有变成被动土压力的趋势。
开挖到坑底标高时，围护墙上的水土压力实测值与计算值的分布如图7所示。
主动侧最下面测点的土压力实测值出现增大的反常结果，其原因尚不清楚。PLAXIS计算的结果也出现增大的现象，并且在上端的支撑附近位置，土压力的值比水土合算与分算的结果都要大，主动侧的土压力两端增大，呈“R”形分布，表明多道支撑的基坑主动侧土压力有呈“R”形分布的趋势;在被动侧，PLAXIS计算的水土压力分布较其他方法计算的结果更接近于实测值，表明随着开挖深度的变大，考虑土层渗流固结变形和围护墙位移作用计算的土压力较符合实际的土压力分布情况。
将图4至图7各个工况下的土压力实测值与计算值汇总可以发现，传统的水土压力分算与水土压力合算在围护墙主动侧的计算值，在基坑不同开挖工况下均相等，其计算的主动土压力不能反映渗流条件及开挖工况的影响;在被动侧，开挖深度较小时的水土合算结果偏大，与实测值存在较大差距。PLAXIS计算的水土压力与实测值能较好地吻合，说明考虑渗流固结耦合的算法能综合反映土层的渗流固结及围护墙位移变形与水土压力的相互作用，其分析结果合理可靠。
4 结论
本文总结了在基坑内外存在水位差，土层中地下水有渗流条件下围护墙上的水土压力计算方法，通过工程算例分析了不同计算方法的区别，可得以下结论。
(1)基坑周围土层为分层土时，根据地下水绕墙渗流考虑水头损失和渗透力的水土压力计算，能反映不同土层对水压力和土压力的影响，降低了围护墙主动侧的水压力，为基坑设计优化提供条件。
(2)对粉土和黏性土采用的“水土压力合算法”，计算的墙后水土压力与实测值较为吻合，但是，计算的墙前水土压力偏大。当基坑开挖较浅、围护墙变形较小时，水土合算的极限状态的土压力与基坑开挖的实测土压力差距较大，在软土基坑中，水土压力合算确定的嵌固深度将偏于不安全。
(3)利用岩土有限元PLAXIS程序可以考虑基坑土体的流—固耦合效应，根据基坑开挖过程中地下水渗流边界和围护墙变形计算的水土压力分布，在被动侧和实测值较好地吻合，随着开挖深度的增大，主动侧计算结果与实测值逐渐吻合，其与水土合算的值相比偏于安全。
(4)基坑围护墙上的水土压力计算，应仔细分析土层分布、围护墙的受力条件以及开挖渗流过程中围护墙的位移变形，利用岩土工程有限元程序，通过流—固耦合来考虑渗流对土压力的影响，可以反映不同条件的水土压力分布情况，计算结果合理可信，为复杂条件下基坑围护墙的水土压力计算提供了新途径。
摘要：深基坑开挖基坑内外产生水位差,地下水绕过围护墙底部向基坑内渗流影响围护墙上的水土压力分布。本文总结了基坑围护墙上水土压力的计算方法,针对基坑周围土体成层分布的特点,分析了传统的水土压力分算、合算以及考虑渗流影响的水土压力分算方法。利用岩土工程有限元方法,通过渗流—应变耦合来计算渗流条件下的水土压力,结合工程算例,比较了传统水土压力计算方法与有限元方法的计算结果以及与实测土压力的差别。发现考虑固结变形的水土压力计算值与实测值能很好地吻合,而传统的水土压力计算结果与实测值偏差较大,其计算的假定条件存在误差以及没能考虑渗流和不同开挖工况的影响。考虑渗流—固结耦合的有限元方法能综合考虑不同土层和边界条件以及基坑围护墙的受力条件与变形的影响,其计算的水土压力分布是合理的。
热-固耦合 篇5
微流体系统中大都要求有微泵。其中,薄膜往复振动式无阀微泵是以扩散/收缩单元为阀门,通过泵膜的振动驱动流体,其制造简单,可以驱动一些非均相的流体,是众多研究的焦点。根据驱动原理不同,无阀微泵可划分为电磁、静电、压电、形状记忆合金等多种类型。其中,静电无阀微泵具有较低功耗,易于与IC工艺集成兼容等优点,而得到广泛关注。
由于微泵的尺度很小,薄膜往复振动式无阀微泵的绝大部分物理量难以通过实验测量,特别是瞬态量。因此,在微泵的研制阶段,利用数值方法对动态性能进行仿真、预测显得很重要。它不仅可降低费用,而且能更好地了解微泵的工作原理以及可能的潜在缺陷。但因微泵中柔性泵膜、电驱动和流体相互耦合,加之流体在扩散/收缩单元中的不同流动方向表现出的不同压力特性,增加了对整个泵的仿真难度和复杂性。早期的研究采用低阶集总参数模型[1,2]和等效电网络模型[3]等简化方法,忽略了参数的空间分布特性。目前大部分的研究没有考虑泵的流固耦合效应,也无法建立电驱动信号量与流体动力学量(比如流量)之间的关系。或是研究驱动器的动态特性[4,5,6,7];或是简化流体场,根据扩散/收缩单元的压力损失系数,将泵膜的运动与泵腔内外压力以及扩散/收缩单元流量联系起来,研究泵膜的动态特性[8,9,10,11];或是将驱动器的位移输出,作为流体场的移动壁面问题来研究泵流体的动力学特性[12,13,14]。实际上,泵膜的动态特性与流体黏滞损失的非线性和不稳定性等是相互影响的。为了理解整个微泵的动态特性,模型必须能有效地描述这些耦合效应[15,16]。本文以静电无阀微泵为对象,建立静电-结构-流体全耦合的3维模型,利用数值方法,仿真并分析各个变量的非线性动态特性。
1 模型的建立
1.1 控制方程、边界条件以及初始条件
图1所示为静电无阀微泵的结构示意图。对于泵送流体,采用任意拉格朗日-欧拉(ALE)描述比较方便。某一时刻,在参考系中
式中,
流体域Ωf的不可压缩流动动量守恒定律以及连续方程可以表示为
式中,ρf、f、σ分别为流体密度、体积力向量和Cauchy应力张量;[0,T]为所考察的时间间隔。
牛顿流体的本构方程为
σ=-pI+2μ∇s
式中,I为单位张量;μ为流体的黏度;∇s
在泵膜流固耦合界面Γf-s处,边界条件如下:
式中,
式(5)为流固界面Γf-s处的无滑移条件;在流固界面Γf-s上,参考系的速度
对于泵膜,可以利用标准的Lagrangian描述建立其运动学方程,即
式中,ρs为固体的密度。
泵膜流固耦合边界Γf-s满足:
边界条件式(9)、式(10)分别对应于式(5)和式(7)。
如果没有包含自由电荷,则描述静电场的Poisson方程可以表示为
∇·ε∇ϕ=0 (11)
E=(E1,E2,E3)=-∇ϕ
式中,ϕ为静电势;ε为介电常数张量;E为静电场电场强度。
Maxwell静电应力张量TE由下式计算:
式中,FE为静电载荷。
各个场变量的初始条件都设置为零。为了计算易于收敛,将一个很小的初值输入到模型中。
1.2 模型参数及材料属性
图1所示静电微泵的泵膜厚度为10μm,泵腔深度为100μm,泵膜半径为2000μm;扩散/收缩单元的几何结构相同,长度为1000μm,最小端面宽度尺寸为80μm,扩张角为10°;连接的进出口管道长度皆为2000μm,宽度为200μm。流体介质为去离子水,物性参数如表1所示。
假设层流流动,考虑对称性,在CFD-ACE中建立微泵的一半模型。利用结构化网格划分流体,单元数为52 880个,泵膜划分为4974个壳单元。使用基于压力修正的SIMPLE-C 算法对流动以及能量进行顺序积分,获得流场的解。流体域的对流和扩散采用一阶迎风格式,时间积分采用了Crank-Nicolson格式。使用基于有限体积法的计算格式求解静电场,这样可以处理不同介质介电常数问题。利用一阶壳单元来表达泵膜,并进行大变形和接触等几何非线性分析。考虑到静电、结构以及流动耦合求解的需要,利用网格重划分工具,使用预测-校正方法对局部变形进行连续修正,将流动速度和结构速度联系起来,实现流固耦合分析。流体和结构求解器的耦合迭代频率设置为fc=1,使得流体场和结构场信息能及时得到交换和更新。当量纲一残差小于10-4时,认为计算是收敛的。施加的周期性驱动电压V=200(1-cos(20πt))(泵腔底电压为零),在Pentium-3.0GHz的PC机中,完成3.5个周期的求解,共350个子步数,求解时间约为59h。
2 结果及讨论
2.1 泵膜的动力学特性
图2所示为计算得到的泵膜中心点垂直方向上的位移随时间变化曲线,图中最大变形量为28.433μm。考虑流固耦合效应后,在流体阻尼作用下,泵膜无法回复至初始零位置,0.653μm为新的平衡位置。稳态情况下,泵膜位移的变化滞后于输入电压信号的变化,一个周期的滞后时间约为0.01s。
2.2 流场的瞬态特性
实际上,泵流体动态特性更为设计者所关注,为了更好地诠释泵膜振动与流体特性之间的耦合关系,用图3所示的泵膜流固耦合界面上中心点处的垂直方向位移与泵腔内流体压力来显示其对应关系,从图中可以明显地看出,流体压力的变化与泵膜的位移存在着密切关系。在静电力的作用下,泵膜下拉变形的速度大于回弹的速度,所以在流体压出泵腔的过程中其压力变化较大,最高压力为945.1Pa;而在流体吸入泵腔过程中,其压力变化相对较小,最大负压为-558.0Pa;整个泵送过程中,流体压力呈现非对称性变化。图3还表明,泵初始工作点从零点开始,经过第一个周期的瞬态过程后,第二个周期的轨迹与第三个周期的轨迹已经重合,表明泵已经进入稳定工作状态。流体压力的波动周期同样滞后于驱动电压的变化(约0.01s)。
图4所示为进出口单元处的瞬态体积流量变化曲线,从图中可以看出,在泵膜下拉(往固定电极)运动阶段(泵送阶段),出口单元的瞬态流量大于进口单元的瞬态流量;而在泵膜回弹阶段(泵吸入阶段),进口单元的流量则大于出口单元的流量。最大的瞬态流量差发生在泵腔内压力达到极大值时刻,而在泵腔内压力与进出口压力相差不大的时间段,进出口单元的流量相差不大。表明要使扩散/收缩单元发挥“整流”作用,其两端的压力差必须达到一定值,并且压力差要尽可能地大。对稳态情况下的泵净流量进行积分计算,得到的净输出流量为10.764nL/s。
图5所示为微泵横截面方向上流体压力的分布情况,由图可知泵腔内的压力几乎分布一致,说明在平面布置的无阀微泵中,采用Reynolds方程来描述泵腔内流体的瞬态动力学特性具有一定的可行性。另外,微泵所连接的进出口管道对微泵的工作性能也有影响。
考虑了流固耦合特性后,模型不仅能够仿真泵流场的动态特性,还能够仿真泵膜的应力特性变化。图6所示为t=28ms时刻的泵腔内流体压力分布与泵膜应力分布情况,以及这一时刻的扩散/收缩单元内流体的速度场分布。此时,泵腔内流体压力达到最大值945.3Pa,泵膜中心挠度为13.20μm;泵膜边缘处VonMises应力为28.39MPa。微泵腔内的流速较小,扩张管/收缩管颈部的平均流速达到最大,计算得到的平均流速为0.0867m/s,雷诺数为7.71,远小于宏观条件下通常认定的临界雷诺数(2000左右),这也证实了前面假设层流流动模式的正确性。另外,图6b和图6c表明,收缩/扩散管内未出现流体与固体壁面分离流动现象。
3 结论
(1)扩散/收缩单元的最大流量差都发生在泵腔内压力与进出口压力相差最大的时刻,表明要使得扩散单元发挥“整流”作用,其单元两端的压力差必须达到一定值。
(2)在仿真的每个时刻,泵腔内流体压力分布几乎一致,表明用Reynolds方程来描述泵腔内流体动力学特性具有可行性。
(3)仿真得到的最大雷诺数远小于通常认定的临界雷诺数,表明微泵流体具有层流特性。
(4)泵腔内流体动态特性与泵膜的运动存在着密切关系,如果忽略了流固耦合效应,简单地将微泵的驱动器与腔内的流场分离处理,必然歪曲了泵流体场的动态特性。泵膜变形、流体压力和进出口流速等响应与输入电压信号之间存在着滞后现象。对微泵的3D全耦合仿真,能够获得驱动电信号与流体流量、泵送压力等输出变量之间的直接关系,有利于从整体意义上实现微泵的优化设计。
发动机流固耦合传热自动优化分析 篇6
关键词：内燃机,水套,计算流体力学,流固耦合,自动优化
0概述
发动机冷却系统不仅要在高转速高负荷时带走热量,控制发动机温度在合理范围,还要满足低温低速低负荷时的热机需求[1]。发动机水套应保证在最小容积下带走最多的热量。单纯水套的计算流体力学(computational fluid dynamic,CFD)分析已成为水套设计的通用手段[2,3],该分析可通过对比冷却液的流速、传热系数来选出最优水套,但不能反应水套流域的增减对结构温度的影响。流固耦合技术的出现正好弥补了此缺点。间接的流固耦合需要流体计算与结构温度计算互为边界,并需要反复多轮的相互映射边界,才能得到准确结果,其计算繁琐且周期很长[4]。直接的流固耦合计算同时建立流体和固体的网格模型进行传热计算,结构的温度分布直接反映了水套的优劣。在简化计算边界保证精度的同时大大缩短了计算周期。多学科优化技术在仿真领域越来越受重视,一维计算的参数优化及外形的拓扑优化应用较多,而流域的CFD自动优化分析因为其复杂性及对计算机硬件需求较大,仍处于起步阶段。而借助自动优化技术得到最优化设计对于经验积累较欠缺的中国自主发动机开发非常重要。
应用商业流体软件对某款四缸发动机缸套、冷却液、机体进行稳态流固耦合传热分析,并与试验值对比验证其有效性;应用Java软件将计算流程自动化编程,应用网格变形软件对结构进行参数化变形;最后,在多学科优化软件的平台下进行试验设计,基于试验设计的结果进行优化分析,得到最优水套设计。
1 理论基础
1.1 传热模型
发动机传热过程从开始传热到热平衡一般超过10min,缸内燃气的瞬态效应可以忽略。而本文重点关注发动机在全速全负荷工况下达到热平衡后的温度,因此发动机传热过程可简化为稳态传热问题。高温燃烧的气体经过燃烧室壁面传到机体,机体的高温侧热量通过热传导传到低温侧,并通过水套中的低温冷却液带走热量,最终把发动机温度控制在合理范围内。平衡后,单位面积的换热量应相等,燃气侧和冷却水侧与壁面的换热可采用第三类热边界条件,如式(1)所示。
式中,h为对流传热系数;λ为结构热导率;T为结构温度;n为结构厚度;Tw为壁面温度;Tf为壁面附近流体的温度[5]。
1.2 优化方法与优化流程
优化是在总结一定样本结果的基础上,得到目标值与变量的规律并寻优。一般用试验设计(design of experiment,DOE)的方法来设计样本点。常用的DOE方法有全因子设计、部分因子设计、拉丁超力方设计、田口法等。优化的算法有梯度优化、直接搜索与全局优化算法等[6]。试验设计方法与优化算法的选择均需要权衡变量数量、优化目标及时间限制。
优化技术是将传统的水套试错改进方法进行参数化与自动化,计算网格根据定义的变量范围自动变形,不需要设计的二次介入,变形的网格自动输入给流体分析软件进行流固耦合分析,得到流场和温度场结果输入给优化软件进行参数敏感性分析,根据需要建立响应面模型,设立优化目标进行优化分析,将计算的最优方案输出给设计进行设计修改,修改的模型还需要进行流固耦合分析验证是否达到了预期的优化目标。达到最优化则设计完成,否则仍需要通过修改变量、增加样本点、改变试验设计方法或优化方法等措施重新优化。优化流程如图1所示。
2 整机与水套流固耦合分析
2.1 流固耦合传热分析模型
基于STAR-CCM+软件建立了某四缸汽油机的整机温度-水套耦合分析模型。模型包括缸体、缸盖、气缸垫、气门、气门导管、垫圈等,如图2所示。
计算模型将冷却液流经的水套区域封闭,形成单独的流体计算域;燃气侧相关的面区域保留完整以便施加燃气对流换热边界条件。结构本体进行了适当的简化,没有考虑油道与附件。
2.2 边界条件
发动机在标定转速全负荷运转时整机结构的热负荷最大,对水套的冷却能力要求最高,此时缸盖最高温度限值为260℃,缸体最高温限值为220℃。整机传热的流动换热边界均取自标定转速全负荷工况。整机传热主要包括燃烧释放热量、冷却液带走热量、外部空气带走热量、活塞摩擦生热等。燃烧热源通过对缸内燃烧过程进行三维瞬态CFD模拟得到,并对缸内0°~720°曲轴转角瞬时燃气温度和传热系数进行平均后映射到流固耦合传热分析模型中的燃气侧面网格上[7]。图3为映射后的气道与燃烧室侧的平均温度与平均传热系数。由于冷却液与整机的对流换热同时计算,冷却液侧的边界只需要给定水泵入口的流量120L/min,入口温度90℃,水套出口outflow。外壁面设定为与空气的强制对流,传热系数80W/(m2·K),空气温度50℃。未考虑活塞摩擦生热。
2.3 计算结果分析与验证
缸体缸盖的耦合计算结果包括了水套的流动特性与整机的温度分布。图4为缸体与缸盖水套的传热系数分布云图。该水套缸体传热系数整体偏小,1#缸到4#缸依次递减,从底部到顶部基本均匀分布,并没有在顶部火力岸附近加强换热。缸盖部分整体传热系数高于缸体,排气门鼻梁区传热系数从1#缸到4#缸逐步加强。图5为缸体缸盖的结构温度分布。缸盖最高温区域出现在排气门鼻梁区中间,最高温度不超过208℃,远低于缸盖温度限值260℃。缸体最高温区域出现在各缸之间的鼻梁区上方,温度不超过190℃,低于缸体最高温度限值220℃。
为了考察流体软件计算结构温度的准确度,将流体软件STAR-CCM+直接耦合计算的缸体缸盖温度值与结构软件ABAQUS计算的温度值进行了对比。由于试验限制,将缸体缸间温度与试验值进行了对比。缸盖没有相对应的实测值,缸盖的取点位置分布在排气门鼻梁区中间往下偏移5mm位置,缸体测点位置相对于各缸之间的鼻梁区上顶面向下偏移5 mm。表1为各测点的计算值与实测值。STAR-CCM+计算的缸体温度与ABAQUS计算值及实测值比较接近,其分布趋势有所偏差,最高温度出现在1/2缸间,而不是3/4缸间。影响因素一是没有考虑活塞的摩擦,二是缸垫直接与缸体缸盖传热,没有考虑实际的接触热阻,导致缸盖的温度分布对缸体有一定的影响。STAR-CCM+计算的缸盖温度比ABAQUS计算的温度偏小5~8℃。STAR-CCM+计算值与ABAQUS计算值及实测值总体吻合得较好,基本能反映发动机本体的温度水平。
3 优化分析
由于四缸机水套容积大导致整机温升慢。但计算与实测的缸体与缸盖温度远小于温度限值,具备一定的优化潜力。基于缸盖水套模型复杂,容积不大,而缸体水套容积为同排量相似发动机缸体水套容积的两倍,且缸体水套相对简单。因此此次优化只针对缸体水套。变形软件采用Meshwork,优化软件采用Optimus。
3.1 优化变量设置
缸体水套优化基于四缸机的原缸体模型。为了说明优化变量,过1#缸、2#缸体水套中间位置截取剖面A-A,截面位置及变形前后的水套剖面如图6所示。其中,实线为原缸体水套;虚线为变形后的缸体水套。优化变量主要为进排气外侧水套往内侧移动以便减薄水套厚度的变形量W;水套两侧底部向上移动以便减少水套的深度的变形量D1;水套各缸之间的中间冷却区域向上移动以便减少缸间的冷却深度的变形量D2。
3.2 DOE分析
变量W、D1与D2的变化范围如表2所示。为了减少计算量,DOE计算时每次只变化一个变量,其他两个变量为0。输出值为缸体、缸盖的最高温度与水套压损。
3.3 优化结果分析
图7为Optimus优化模型图。在Optimus平台下调用Morpher变形运算及将变形后的模型输出给STAR-CCM+自动进行整机流固耦合传热分析,最后将计算结果输出给Optimus进行统计分析。图8为缸体最高温度随着变量的变化曲线。缸间变薄的变量W对温度影响最大,7mm之前缸体水套变薄,流速加快,传热加强,缸体最高温度最多降低12℃,但继续减薄温度又会有所上升。深度变浅的变量D1与D2对温度的影响小于3℃,而且没有明显的变化规律。
由于优化的目标是在保证缸体与缸盖的冷却前提下尽可能地减小缸体水套容积,水套容积很难作为CFD计算的输出值,因此本文没有进行自动寻优,而是根据计算结果和工艺限制,选取了W、D1与D2的最终值为7、40、40mm。根据设计部门提供的最终优化模型进行流固耦合验算,最高温仍出现在各缸之间的鼻梁区上顶面,缸体的温度最高值由原方案的196.6℃降低到185.0℃,缸体在各缸之间的鼻梁区上顶面向下偏移5mm的位置测的缸间温度相对于原方案降低约4~5℃。缸盖由于结构与热负荷均未改动,排气门鼻梁区的最高温计算值与原方案相差小于1℃。水套的压损由原方案的18.1kPa增加到20.4kPa,压损增加可接受。缸体水套容积由原来的2.3L下降到1.0L,容积减小相当可观。优化分析对于挖掘这种设计方案的潜力效果非常明显。
4 试验验证
图9为优化水套的砂芯快速成型件。发动机缸体框架整体结构基本不变,只减小水套容积的情况下,在-20℃的低温暖机温升试验时,暖机时间相对于原方案提高了15%,水套容积的减小使得发动机热机性能得到了较大提升。标定转速全负荷工况下缸体在各缸之间的鼻梁区上顶面向下偏移5mm的位置测的缸间温度平均降低2~3℃,实测缸间温度降低幅度稍小于仿真值,但优化水套有利于减小缸间温度的结论得到了验证。
5 结论
(1)整机流场-温度场直接耦合计算在得到流场信息的同时能够得到结构的温度分布,通过直接评价结构温度来评价水套更加直观和准确。
(2)流体软件计算得到的整机温度与结构软件计算结果及实测值基本吻合,易于实现计算自动化与优化分析。
(3)基于Optimus平台调用Morpher网格变形与流体软件计算,可以实现多样本点分析,并能基于样本结果寻优,从而简化仿真分析流程及缩短分析时间。
参考文献
[1]HEYWOOD J B,WELLING O Z.Trends in performance characteristics of modern automobile SI and diesel engines[C]//SAE Paper.Cambridge,Massachusetts,USA,2009,2009-01-1892.
[2]STEPHEN S,EDWIN I,JUN X.Engine knock toughness improvement through water jacket optimization[C]//SAE Paper.Pittsburgh,Pennsylvania,USA,2003,2003-01-3259.
[3]郑清平,张惠明,黎苏.基于三维CFD技术的发动机冷却水套流动分析[J].内燃机工程,2009,30(6):37-40ZhENG Q P,ZhANG H M,LI S.Flow analysis in cooling water jacket of engine based on three dimensional CFD technology[J].Chinese Internal Combustion Engine Engineering,2009,30(6):37-40.
[4]李迎,俞小莉,陈红岩,等.发动机冷却系统流固耦合稳态传热三维数值仿真[J].内燃机学报,2007,25(3):252-257.LI Y,YU X L,CHENG H Y,et al.3Dsimulation of steady heat transfer of fluid solid coupled system in engine coolant system[J].Transactions of CSICE,2007,25(3):252-257.
[5]杨世铭,陶文铨.传热学[M].4版.北京:高等教育出版社,2006:46-47.
[6]赖宇阳.Isight参数优化理论与实例详解[M].北京:北京航空航天大学出版社,2012:94-95.
降雨条件下的边坡流固耦合分析 篇7
关键词：ABAQUS,耦合原理,降雨入渗理论,边坡稳定性
0前言
降雨是诱发滑坡的重要因素, 全世界每年都会由于降雨发生大量的滑坡事件, 给人们的生产生活带来巨大的影响, 受到各国政府与机构的高度重视。1982年美国加利福尼亚旧金山湾地区34 h内降雨616 mm, 在10个县内诱发了数千处滑坡、泥石流, 造成25人死亡, 6千6百万美元的直接经济损失。2012年10月我国云南省彝良县发生山体滑坡, 滑坡体积约4.5万m3, 造成了19人遇难、1人受伤, 损坏教室3间、村民房屋9间。滑坡体堵塞镇河, 形成蓄水量约1万m3的小型堰塞湖。
土体作为一种多孔介质, 雨水渗入土体实际上就是一个渗流场与应力场的耦合过程。一方面渗流场的变化会对土骨架造成一定的影响, 孔隙水压力的改变影响了作用在土骨架上的有效应力, 土骨架发生变形, 从而造成边坡内的应力场与位移场的变化;另一方面边坡内应力场发生变化, 土骨架体积变形导致土体的孔隙率发生变化, 进而改变土体的渗透系数导致边坡内的渗流场发生改变。因此, 研究降雨条件下的边坡流固耦合分析十分重要的。
1 流固耦合分析
1.1 有效应力原理
土体作为一种多孔可压缩材料, 可以模拟成多相体, 满足有效应力原理。材料孔隙率为孔隙体积和总体积的比值, 即。
在ABAQUS中通常使用孔隙比e, 并不采用孔隙率, 孔隙比与孔隙率的换算关系为:
饱和度s定义为水的体积与孔隙总体积的比值为:
土体为饱和时s=1, 土体干燥时s=0。
土体中一点的应力状态服从有效应力原理:即总应力是由平均孔隙水压力乘以系数x和作用于土骨架上的有效应力组成的。
由于实验研究数据较少, ABAQUS将x取为饱和度s。
1.2 应力平衡方程和渗流连续方程
1) 平衡方程。土体应力平衡方程采用虚功原理来表示, 即某一时刻t土体的虚功与作用在该土体上作用力 (体力和面力) 产生的虚功相等, 即
2) 渗流连续方程。渗流连续方程是由同一时间内流入土体的水量等于土体体积变化量这一连续条件来建立的, 即
式中, vw为流体的平均速度;n为边界S外法线, 该方程对ρw0进行了归一化。
在非饱和的土体中渗透系数的表达公式为:
式中, 为非饱和土的渗透系数;ks为考虑饱和度对渗透系数影响的一个折减系数, ABAQUS中ks=s3;k为饱和度的渗透系数。
2 降雨入渗分析
降雨入渗实际上就是水分进入土壤并在非饱和区与饱和区运动的过程。根据研究的结论, 含水量的剖面可以分为4个区:饱和区、含水量有明显下降的过渡区、含水量变化不大的传导区和含水量迅速降到初始值的湿润区, 每个区都有自己的特征。
降雨的初期, 在雨水施加于土体表面后的很短的时间内, 地表的含水量很快的由初始值增大到最大值。随着降雨入渗的进行, 水会不断的向边坡内部下移, 含水量的分布曲线由比较陡直逐渐变得相对平缓。
对于边坡而言, 实际入渗流量为qs (t) , 且垂直于坡面方向, 可以得到降雨强度和实际入渗流量的关系为:
qn (t) 为法线方向的降雨强度, Rn (t) 为法线方向的入渗率。
3 数值算例分析
以某一均质边坡作为研究对象, 该边坡土体饱和时的渗透系数为0.018 m/h, 弹性模量和泊松比分别为E=10MPa、V=0.3, 初始孔隙比为e=1, 凝聚力c=15 k Pa, 摩擦角φ'=30°。坡角为40°, 初始地下水位位于坡脚处, 边坡尺寸如图2所示。
计算采用的边界条件为:底面为零流量边界, 左右两侧为定水头边界, 水头都为10 m, 上部降雨入渗雨量大小作为定流量的边界。模拟降雨72 h, 0~24 h内降雨强度由0增大到20 mm/h, 24~48 h内保持20 mm/h降雨强度不变, 48~72 h内降雨强度由20 mm/h减小到0 mm/h。
为了获得初始孔压和饱和度的分布情况, 计算过程中先不考虑降雨的影响进行稳态渗流计算, 然后考虑降雨入渗情况下的瞬态分析, 在给定降雨强度的条件下获得不同时刻的孔隙压力、边坡位移、塑性区等分布情况。
对比图3~4可以发现随着降雨入渗的进行, 非饱和土体中的孔隙水压力不断地上升, 孔压由-200 k Pa变化到了-43.75 k Pa。非饱和区范围减小, 基质吸力也相应的减小, 导致土体抗剪强度降低, 影响边坡的稳定性。由图4~5可知降雨减小以后, 随着时间的延长, 孔隙水压力又下降到了-122.7k Pa, 土体浅层的基质吸力有所增加。
由图6~7可知最大的水平位移发生在坡脚处, 为10.25 mm;最大的沉降发生在边坡的中部, 为4.651 mm。在坡顶的位置, 当降雨入渗后, 基质吸力降低, 孔压增大, 有效应力有所减小, 出现了卸载回弹的现象所以最大的沉降没有出现在坡顶位置。由图8的位移矢量图可知, 由于降雨的入渗, 边坡有滑动变形的趋势, 可以推断边坡的稳定性是有所降低的。
由图9~10可见降雨入渗条件下, 在坡脚的浅层首先出现塑性区随后沿着坡面向上发展。由图11知大约在t=15 h在坡脚开始出现塑性区, 随后降雨强度达到最大, 塑性区也随时间快速增加, 大约到了50 h后, 降雨强度开始减小, 塑性区不再扩展, 等效塑性应变也保持不变。
图12为边坡坡面左上角顶点的场变量随水平位移变化的曲线。以位移的拐点作为评价标准, 可以看出在降雨前边坡的安全系数Fs为1.51, 但是受到降雨的影响以后安全系数降Fs低为1.38, 从数值上可以看出由于降雨入渗的作用导致边坡的稳定性降低。
4 结论
考虑降雨条件下的边坡流固耦合分析, 得出了降雨过程中孔压、位移、塑性区的变化规律。随着降雨的进行土体的含水量和孔压不断地上升导致基质吸力减小或消失, 土体抗剪强度降低影响边坡的稳定性。最大水平位移出现在坡脚处, 塑性区也是形成于坡脚并沿着坡面开展的, 所以降雨条件下边坡的浅层最可能出现失稳现象。
参考文献
[1]王金昌, 陈页开.ABAQUS在土木工程中的应用[M].杭州:浙江大学出版社, 2006.
[2]费康, 张建伟.ABAQUS在岩土工程中的应用[M].北京:中国水利水电出版社, 2009.
[3]张芳枝, 梁志松, 周秋娟.非饱和土性状及其边坡稳定性[M].北京:中国水利水电出版社, 2011.
[4]雷志栋, 杨诗秀, 谢森传.土壤水动力学[M].北京:清华大学出版社, 1988.
[5]王永义, 王专翠, 胡以高.降雨入渗补给规律分析[J].地下水, 1998, 15 (2) :74-75.
[6]Zienkiewicz O C, Humpheson C and Lewis R W.Associated and non-associated visco-plasticity and plasticity in soil mechanics[J].Geotechnique, 1975, 25 (4) :671-689.
[7]刘金龙, 栾茂田, 赵少飞, 等.关于强度折减有限元方法中边坡失稳判据的讨论[J].岩土力学, 2005, 27 (8) :1345-1348.