双向流固耦合（共8篇）
双向流固耦合 篇1
1 引言
流体动力学是研究流体平衡的条件及压强分布、流体运动规律、以及流体与固体之间的相互作用等, 研究结果对分析管道的振动及影响因素有重要意义。本文针对新疆某石化公司的10-K-302C离心式甲烷制冷压缩机自开机以来润滑油管线振动较大的问题, 通过对管内流体流动状态进行模拟分析, 得出了流体耦合前后动力特性的变化及管道振动的原因。
2 双向流固耦合分析原理
流固耦合要遵循质量守恒定律、动量守恒定律、能量守恒, 所以在流固耦合交界面处, 应满足流体域固体应力 (σ) 、位移 (d) 、温度 (T) 、热流量 (q) 等变量的相等或守恒, 即满足下面四个方程:
式中:df、ds分别为液体、固体的边界位移。σf、σs分别为液体、固体应力。
3 流体和管道的计算模型
就10-K-302C离心式甲烷制冷压缩机装置的润滑油管线位移较大现象, 通过分析润滑油耦合前后的动力学特性, 找出流体运动特性, 对寻找该管道振动原因有重要指导作用。出口管道的管路图如图1:选取润滑油在弯管中心轴线处的1、、2、、3、、4、点, 及在出口处5、为观测点。
用Workbench有限元分析软件对管道结构进行建模, 并利用workbench自带的填充功能生成润滑油模型。将实际的润滑油简化为不可压缩、无粘性、平均密度和压力在整个流域中均一的流体。分析时设定管道内润滑油的材料属性:质量864Kg/m3, 动力粘度为4.33cp, 比定压热容为2063J/ (Kg.K) 。建立润滑油模型, 为了更好地描述边界层处参数变化情况, 在网格划分时设定膨胀层。
4 流体耦合前后分析对比
4.1 流体耦合前后的压力和速度对比
润滑油耦合前后的速度流线图2和图3, 可以看出耦合前润滑油的最大流速达到了5.168m/s, 耦合后的最大速度为4.555m/s, 较耦合前变小, 但耦合前后流体速度均在弯头处较大,
耦合前后润滑油与管道接触壁面的压力云图4和图5。绝对压力均在入口处较大, 弯头处较其连接处的直管压力较大。耦合前润滑油壁面的最大绝对压力为772KPa, 最小绝对压力为759.9KPa, 压力波动值为1.58%, 压力波动较小。流固耦合后接触壁面的压力大小和分布与耦合前几乎相同。
4.2 耦合前后流体观测点随时间的变化对比
观察观测点在耦合前后轴向速度和压力随时间变化曲线 (图6、7) 。从图6、7可以看出, 润滑油在耦合前后各个弯头处的测试点的速度相差较大, 同时达到稳定的时间不同;耦合前测试点3、4、达到稳定时间较长, 出口处的稳定速度较入口变大;耦合后2、、3、、4、速度波动较大且要达到稳定的时间比较长;耦合前后弯管处的1、、2、、3、、4、的压力同时达到了最值;耦合后除了初始阶段压力值波动较大外, 稳定时间和稳定值都一样。
4.3 耦合前后流体中心轴线处的压力速度对比
耦合前后润滑油从入口到出口的中心轴线处的压力和速度变化情况如图8和9。可以看出耦合前后的流体压力和速度变化情况基本一致, 且都是在管道弯头处压力和速度变化都达到大。
5 结论
(1) 车间提供的出口管路的测试压力波动范围在660-680KPa之间, 软件模拟分析的流体的压力波动范围在660-672KPa之间, 与利用cfx计算结果的差值为1.17%, 说明了可以用双向流固耦合有限元模拟方法进行管道流体的分析研究。
(2) 在润滑油速度为1.93m/s的情况下流体耦合前后的流体压力变化很小, 说明流体的压力能在流固耦合后变化很小, 对管道产生的影响较小;耦合前, 润滑油最大流速为5.168m/s, 耦合后的最大速度为4.555m/s, 流体速度较耦合前减小了11.86%, 由耦合时能量守恒可知, 流体的这部分动能传递给了管道。
(3) 耦合时流体的速度波动比较大, 特别是观测点2、3、4处;耦合时流体速度达到稳定的时间较耦合前达到稳定的时间较长, 这样对管道的影响也会变得较大。耦合后最大速度比设定速度1.93m/s大了136%, 且在弯头处变化最大, 这样会对管道产生很大的冲击力, 造成管道的振动。
参考文献
[1]谢龙汉.ANSYS CFX流体分析及仿真[M].北京:电子工业出版, 2012, 1-20
[2]赵兴艳.CFD方法在流体机械设计中的应用[J].流体机械, 2000, 28 (3) ;22-25
[3]宋学官.ANSYS流固耦合分析与工程实例[M].北京:中国水利水电出版社, 2012, 3-13
双向流固耦合 篇2
应用流固耦合数值方法研究机翼的射流减振技术
采用流固耦合的数值计算方法研究了NACA 0015翼型在大迎角(15°～60°)范围的颤振,以及在翼型背部部分引入射流的.减振技术.在流体区域用高精度、高分辨率算法求解Farve平均的Navier-Stokes方程,在固体区域用四阶Runge-Kutta法求解振动方程,并且每一个时间步后都在两个区域之间传递边界条件.计算结果表明在翼型背部引入适当射流会降低翼型的振动,并提高升力.但如果引入射流的速度过高,会在叶背处形成新的分离流,升力反而会下降.
作 者：金琰 袁新 作者单位：清华大学热能工程系,北京,100084刊 名：空气动力学学报 ISTIC EI PKU英文刊名：ACTA AERODYNAMICA SINICA年，卷(期)：20(3)分类号：V211.3关键词：NACA 0015 翼型 流固耦合 气动弹性 颤振 射流 分离流
双向流固耦合 篇3
双向流道轴流泵装置由于具有占地面积小,运行管理方便等特点,广泛应用于我国排涝灌溉领域中。近年来,双向轴流泵装置的安全稳定运行越来越受到重视,由于双向流道轴流泵装置的可靠性很大程度上依赖于规范的结构设计,而结构设计的基础是正确掌握泵各部件应力分布及变形的情况。所以在研究泵装置时需要将流场计算与结构计算综合考虑。
流固耦合作用是自然界中客观存在的一种特殊现象,是指流体与固体之间的相互作用。在水力机械领域中,流固耦合理论首先应用于水轮机领域[1,2,3,4]。随着水泵向大型化发展,流固耦合理论也逐渐应用于泵领域。袁寿其等[5]及其他学者[6,7,8,9]分析了叶轮流固耦合作用对离心泵内部流场的影响。裴吉等[10]及其他学者[11,12,13]利用流固耦合理论对离心泵叶轮进行应力及应变分析,施卫东等[14]及其他学者[15,16]通过流固耦合模型对轴流泵叶轮的应力应变进行了数值研究。吴广宽等[17]基于瞬态流固耦合技术分析了混流式转轮叶片裂纹产生的原因。张德胜[18]采用双向耦合方法对蜗壳结构产生的振动位移和振动速度进行了数值模拟。但是,目前为止,在轴流泵领域中流固耦合方法主要应用于泵段的数值模拟,并没有考虑进出水流道的影响。而双向流道轴流泵装置中,双向进水流道由于其形状特征容易产生涡带,进水流态不良将直接导致叶轮室内的流态不良,发生脱流、漩涡等现象,使得叶轮产生强烈振动[19],严重时会导致叶轮叶片产生裂纹,影响泵装置的安全运行。预防叶轮叶片出现裂纹的一种措施是对叶片不同区域进行不同程度的加厚,而叶片加厚是以叶片的应力与变形分布为理论基础。
为了获得叶轮应力及应变的分布情况,本文首次采用双向同步求解的方法,在考虑双向进出水流道的情况下,对轴流泵内流场和叶轮结构响应进行联合求解,通过数值模拟可以获得在不同工况下叶轮叶片的应力及变形规律,不仅可以为叶轮叶片加厚提供理论指导,也可以为双向流道泵装置的优化设计与安全运行提供参考。
1 数值模拟
1.1 计算模型
本文以国内某大型双向流道泵站为研究对象,该泵站在运行过程中会产生一定程度的水力振动,从而产生一定程度的机械振动,使叶轮根部产生不同程度的裂纹。以下将以该泵站的模型泵装置为研究对象进行计算。
双向流道轴流泵装置主要包括双向进出水流道,轴流泵叶轮,轴流泵导叶。其中,叶轮叶片数为3,导叶叶片数为5,其结构如图1所示。模型泵装置设计参数为:流量Q=0.33 m3/s,转速n=1 550 r/min,扬程H=2.76 m。
1.2 网格划分
利用ICEM对流体域进行网格划分,其中叶轮水体与导叶水体采用结构网格,进出流道采用非结构网格。利用ANSYS对固体域进行网格划分,固体域只考虑叶轮结构模型。流体域网格与固体域网格如图2所示。水体域网格总数为4 030 944,固体域网格数为18 540。
1.3 边界设置
流体域计算中,进水流道进口采用质量流量,前池的表面设置为自由水面,自由水面对速度和湍动能均采用对称平面处理。出水流道出口采用固定总压,总压设定为101 325 Pa,采用自由出流边界条件。所有壁面为光滑壁面,采用无滑移边界条件。将叶轮与进水流道,叶轮与导叶的交界面设置为瞬态冻结转子。将导叶与出水流道的交界面设置为普通连接。将叶轮与流体的接触面设置为动网格。叶轮每旋转3°作为一个时间步,时间步长为0.000 32 s。
叶轮材料选用结构钢,叶轮材料特性参数如表1所示。定义叶轮轮毂圆柱面为固定约束,设置叶轮叶片为流固耦合作用面。固体域的时间步长与流体域的时间步长保持一致。
双向耦合计算的顺序是流体域计算与固体域计算同步进行,并通过耦合作用面进行数据传递,当固体域计算与流体域计算同时收敛时,进入下一步的耦合计算。具体计算步骤如图3所示。
为了保证双向流固耦合的精度,双向耦合的初始值为流体域非定常计算10个周期后的数据。
2 结果与分析
2.1 外特性结果分析
由于低扬程泵装置实验测试误差较大,本文中实验值为原型泵装置运行时的单机扬程值,流固耦合计算值是指模型泵进行双向流固耦合计算后的第二个周期的时均扬程。数值模拟中采用实际尺寸缩小10倍的模型泵,转速利用等扬程几何相似换算定律所得:
式中:H,HM为实际泵扬程,模型泵扬程;D2,D2M为实际泵叶轮出口直径,模型泵叶轮出口直径;n,nM为实际泵转速,模型泵转速;ηh,ηh M为实际泵水力效率,模型泵水力效率。
如图4所示,流固耦合计算值在小流量工况点略低于实验值,在其他工况点皆高于实验值。两条曲线变化趋势基本相同,在最高效率点(Q=0.396 m3/s),流固耦合计算值与实验值相对误差为6.2%,说明双向同步耦合计算结果是可靠的,也验证了所选湍流模型与网格类型的适用性。
2.2 叶轮等效应力及变形分布图
图5(a)和图5(b)分别为叶轮吸力面与叶轮压力面的变形分布。如图5所示,随着流量增加叶轮叶片的变形量逐渐减小,且叶轮吸力面的变形量大于叶轮压力面的变形量。叶轮压力面与吸力面的最大变形量都分布在靠近轮缘的进水边处,且叶轮压力面与吸力面的变形量由轮缘向轮毂方向逐渐减小。分析其原因主要是叶轮叶片厚度由轮缘到轮毂方向逐渐增大,叶片轮缘厚度较小,容易受应力影响产生变形。
图6(a)和图6(b)分别为叶轮吸力面与叶轮压力面的等效应力分布。如图所示,叶片吸力面的等效应力大于叶片压力面,且随着流量增大,叶片的等效应力逐渐减小。叶片等效应力由叶缘向轮毂方向逐渐增大,主要原因是由于叶片形状,越接近轮毂叶轮吸力面越容易出现脱流现象导致越接近轮毂叶片压力面与吸力面的压差越大。
2.3 等效应力时域图分析
表2为最大等效应力与最大变形量值随流量的变化。如表2所示,在Q=0.8Qopt处,出现最大值分别为137.6 MPa与0.086 771 mm。随着流量增加,最大等效应力与最大变形量逐渐减小。
由于在不同工况下的最大等效应力值相差较大,很难在一张图内表现出3个工况下的最大等效应力的变化趋势,故采用最大等效应力方差的形式来表现变化趋势:
式中:X为常数;E(X)为X的平均值;D(X)为X的方差。
图7为不同工况下,最大等效应力的方差随时间的变化。如图所示,在Q=0.8Qopt时,最大等效应力在一个叶轮周期内出现了2个极大值。在Q=1.0Qopt时,最大等效应力在一个叶轮旋转周期内只出现了1个极小值。在Q=1.1Qopt时,最大等效应力在一个叶轮旋转周期内并没有出现极值。可以看出,随着流量的增加,最大等效应力的变化频率逐渐减弱,3个工况下最大等效应力的相位互不相同。
由于在3种工况下,叶轮的最大等效应力都分布在轮毂附近,而轮毂与叶轮的交界面处容易产生应力集中。所以对叶轮前缘根部特征点P1,P2,以及叶轮尾缘根部特征点P3,P4进行分析研究,其中P1,P3位于叶片吸力面上,P2,P4位于叶片压力面上。特征点分布图如图8所示。
图9-图12为P1,P2,P3,P4处最大等效应力方差随时间的变化。如图所示,P1,P2,P3,P4处的等效应力脉动变化趋势基本相同。在Q=0.33 m3/s时,4个特征点处均有比较明显的等效应力脉动,但随着流量增加,4个特征点处的等效应力脉动都逐渐减小。4个特征点中,P1,P2处的等效应力脉动幅值远远大于P3,P4处的等效应力脉动幅值。P2处的等效应力脉动范围最大,为0~3 MPa,P4处的压力脉动范围最小仅为0~0.015 MPa。
图13为P1,P2,P3,P4处的最大等效应力随流量的变化。如图所示,随着流量增大,除了P4处的最大等效应力出现了先降低后上升的趋势,其余特征点处的最大等效应力皆逐渐下降。4个特征点中,P1,P2处的最大等效应力值远大于P3,P4处的最大等效应力值。在3个工况下,P2处的最大等效应力均大于其余特征点处,P4处的最大等效应力均小于其他特征点处。
通过以上比较可以看出,由于在叶轮旋转时,入流造成叶片前缘处的内部流动不稳定且叶片的最大等效应力区分布在叶轮轮毂处。所以P1,P2处的等效应力脉动和等效应力值远大于P3,P4处。以上计算分析结果对泵站在长期运行后叶片根部产生一定程度裂纹起到了很好的验证作用。
3 结语
(1)双向轴流泵装置中,叶轮所受轴向力随叶轮旋转呈周期性变化,其幅值受流量变化影响明显。
(2)双向轴流泵装置中叶轮叶片厚度对叶片的变形量起主导作用,最大变形区集中分布在叶轮前缘与轮缘的夹角处,且叶片变形量由轮缘到轮毂方向逐渐降低。叶轮形状特征导致叶片所受等效应力由轮缘到轮毂方向逐渐增大,最大应力区出现在轮毂附近。
(3)本文中,最大等效应力脉动随流量逐渐降低,在泵装置最高效率点Q=1.0Qopt时,并不是叶轮的最高效率点。进一步证明了对泵装置进行模拟计算时需要考虑进出水流道对泵装置的影响,才能准确预测泵装置的性能。
(4)在双向轴流泵装置中,叶轮进水边所受应力远远大于叶轮出水边所受应力。原因可能是双向进水流道形状导致叶轮进口处流态紊乱,入流造成叶轮前缘处的内部流道不稳定。
综上所述可以看出叶轮根部和叶轮进水边所受应力较大,需对该区域进行特别加厚。叶轮所受应力随流量增加有明显下降,在大流量工况下,为了提高叶轮做功效率可以适当降低叶轮叶片厚度。本文所得结论可以为双向轴流泵装置的结构优化与性能预测提供有效参考。
摘要：双向流道轴流泵装置广泛应用于农田灌溉以及抗洪排涝领域。双向进水流道由于形状特征容易导致叶轮产生振动,严重时叶轮叶片甚至会产生裂纹,影响泵装置的安全运行。故采用双向同步求解的方法,在考虑双向进出水流道的情况下,对轴流泵内流场和叶轮结构响应进行联合求解,分析叶轮表面的应力及变形分布,可以为叶轮进行结构改进提供理论指导。计算结果表明,轴向力随叶轮旋转产生周期性变化,叶轮最大等效应力主要分布在靠近轮毂处,最大变形量主要分布在叶轮前缘与轮缘夹角处。当流量增大时,叶轮最大等效应力与最大变形量逐渐减小。在此基础上,研究了位于叶轮前缘与叶轮后缘上的4个特征点所受等效应力的情况,通过对比发现,叶轮前缘所受等效应力远大于叶轮后缘。
发动机流固耦合传热自动优化分析 篇4
关键词：内燃机,水套,计算流体力学,流固耦合,自动优化
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.
双向流固耦合 篇5
由于爆炸冲击波产生的威力巨大, 并且可以通过计算炸药的各项参数得到产生爆炸冲击波的各项数据, 能够较好的得到控制和利用, 所以爆炸冲击问题一直是工程和军事领域的热点问题, 越来越受到科研工作者的重视。
北京理工大学的赵文杰[1]通过实验和理论分析, 对圆抛物面薄壳型雷达天线在爆炸冲击波作用下的变形特点及有关的变形特征参量进行了分析讨论, 并给出了工程计算式。中国矿业大学力学与建筑工程学院的李清[2]进行了爆炸载荷作用下动态裂纹扩展实验研究。洛阳水利工程技术研究所[3]做了冲击波对工程结构及装备的动载实验, 另外, 有些理论已在工程实际中得到了验证与运用, 如地下隧道试验[4]、地形对冲击波影响试验[5]及人防工事试验[6]。对装备部件进行爆炸冲击与毁伤模拟, 属于结构分析的范畴, 本文采用流固耦合算法对爆炸冲击问题进行仿真分析。
1 爆炸冲击波形成机理分析
炸药爆炸完成后, 产生的高温、高压爆轰产物像一个超音速推进的活塞, 把空气从原来的位置上迅速地排挤出去, 形成一个以超音速运动的、状态参数有突跃的压缩空气层, 即空气冲击波。在此过程中, 爆轰产物的压力在未降到空气初始状态压力前继续膨胀, 不断把能量传输给空气冲击波。当爆轰产物的压力下降到空气初始状态压力时, 空气冲击波形成, 此时爆轰产物的体积称为极限体积。由于惯性作用, 爆轰产物体积继续膨胀, 使压力低于周围空气压力, 所以在空气冲击波尾部会形成一个稀疏区。
爆炸冲击波的形成和传播如图1所示, T1、T2、T3、T4分别表示爆炸后不同瞬间冲击波的状态, 其中T3时刻冲击波开始与爆轰产物分离, T4时刻出现了稀疏区。
2 冲击波毁伤靶板的有限元分析
利用有限元方法对靶板和装备部件进行爆炸冲击毁伤仿真, 采用如图2所示的步骤。
2.1 冲击波模型
近距离爆炸时, 爆炸产物和爆炸冲击波共同对目标进行毁伤, 如果建立炸药模型使之产生冲击波, 则无法将爆炸产物的作用分离出来;另外, 目标与炸药距离扩大, 会使建模空间扩大, 给计算带来不便, 容易造成较大误差。为此, 本文中采用了对靶板直接施加冲击波载荷的方法来研究冲击波对靶板的毁伤效果。
实际情况下的冲击波载荷相当复杂, 一般和结构变形相耦合且不易精确测定, 为研究方便, 略去荷载中的次要因素, 对实际载荷进行理想化处理, 冲击载荷简化后的主要形式有:正弦脉冲、三角脉冲和矩形脉冲。针对爆炸冲击波的特点, 有些学者将爆炸冲击波简化为三角脉冲波的形式。Jones在讨论船舶板的抨击损坏时采用Ochi提出的三角形压力脉冲形式的抨击时间历程曲线[7]。
为使问题的分析过程得以简化, 本文将爆炸冲击波取为三角形脉冲载荷的形式, 载荷-时间历程曲线如图3所示。
如图3所示, 爆炸冲击波作用时间为t, 在t/2时, 冲击波峰值超压达到最大, 由于空气中爆炸冲击波的强度一般在80~130MPa, 所以峰值超压在此范围内取值, 本文选用峰值超压为80MPa和100MPa的冲击波对靶板进行冲击, 以研究其毁伤机理。
2.2 靶板材料模型
由于在战损仿真分析中, 一般将材料等效为LY-12铝合金, 为使仿真分析结果的一致性, 本文中的靶板材料也定为LY-12铝合金, 其材料参数如表1所示[8]:
LY-12铝合金为应变率无关材料, 因此, 可选择*MAT_PLASTIC_KINEMATIC材料模型进行有限元分析。
不同应变率下, LY-12铝合金的损伤应力介于0.55至0.56 GPa之间, 因此将损伤应力取为0.555GPa, 由弹性模量E=72GPa, 可计算出该材料的失效应变为0.77%。其他材料参数为:屈服极限γ=0.31Gpa, 硬化系数β=0.7。
2.3 冲击载荷定义和施加
冲击载荷的定义包括爆炸冲击波作用时间的定义和强度的定义, 且强度随时间变化而变化。为了对靶板施加爆炸冲击载荷, 需定义载荷—时间变量数组, 赋值情况见表2。冲击载荷定义完成后, 对靶板进行冲击载荷施加。
2.4 网格划分与网格密度
为确定合适的网格密度, 在不同的靶板网格密度下, 用强度100Mpa、作用时间0.00500s的爆炸冲击波对LY-12铝合金靶板 (靶板尺寸300×300×8 mm) 进行冲击实验, 结果如表3所示:
由表3可以看出, 随着靶板网格密度的增大, 计算时间显著增长, 但靶板的最大位移基本趋于稳定。为此, 在后续研究中选用了靶板长和宽均50等分, 厚度4等分的划分方法。
3 基于流固耦合算法的爆炸毁伤仿真
对装备部件进行毁伤仿真分析, 实质上就是利用有限元分析中的流固耦合方法, 分析爆炸复合物和冲击波对部件的损伤情况, 研究其损伤机理。多物质流固耦合算法的原理是, 对炸药及其他流体材料 (如空气、水、土壤等) 采用Euler算法, 对其他的结构采用Lagrange算法, 然后通过流固耦合方式来处理相互作用 (*CONSTRAINED_LA-GRANGE_IN_SOLID) , 该方法的优点是炸药和流体材料在Euler单元中流动, 不存在单元的畸变问题, 并且通过流固耦合方式来处理相互作用, 能方便地建立爆炸模型。
本文利用TNT炸药近距离爆炸进行毁伤分析, 将所有物体都选择为三维实体单元, 选定使用LS-DYNA中提供的SOLID164实体单元。
LS-DYNA专门为各种炸药提供了材料模型 (*MAT_HIGN_EXPLOSIVE_BURN) , 再结合JWL状态方程对它进行描述。TNT炸药的材料模型*MAT_HIGN_EXPLOSIVE_BURN中的相关参数为:密度R0=1.7g/cm3, 爆速D=0.753cm/μs, Chapman-Jouget压力PCJ=0.255×1011Pa。
JWL状态方程用于描述压力与体积应变之间的关系, 其形式为
对于TNT炸药, 在g-cm-μs单位制中, 以上方程的输入参数为:A=5.4094, B=0.093726, R1=4.5, R2=1.1, ω=0.35。
4 爆轰产物的破坏作用分析
大量实验表明, 距离爆炸中心10~15r0 (r0为装药半径) 的范围内, 目标受到爆炸产物和爆炸冲击波的复合作用, 而超过10~15r0范围时, 目标只受到爆炸冲击波的作用[9]。
基于上述流固耦合方法进行爆炸冲击仿真, 以装药密度为1.6g/cm3的TNT炸药为例, 爆炸完成瞬间爆轰产物的各项参数如下:压强=196MPa, 密度=2.13g/cm3, 温度=3350℃, 速度=1750m/s。
从这些数据可以看出, 爆炸产物的各项参数都变得很大, 爆炸发生瞬间, 目标在10~15 r0范围内时, 会受到很大的冲量, 导致目标严重毁伤。
从以上分析可以看出, 爆炸产物的破坏威力大, 但作用范围小, 只能对直接接触或近距离接触的目标发生作用。不过随着现代高精武器的出现, 能够使炸药在距目标距离很小时产生爆炸, 所以利用爆炸产物的巨大毁伤作用对目标进行毁伤打击, 已成为现实。所以研究炸药近距离爆炸时对目标的毁伤情况, 具有重要意义。
参考文献
[1]赵文杰, 蒋浩征, 王秀兰.爆炸冲击波对圆抛物面薄壳雷达天线毁伤效应研究[J].弹箭与制导学报, 2000, (1) :8～13.
[2]李清, 杨仁树等.爆炸载荷作用下动态裂纹扩展试验研究[J].岩土力学与工程学报, 2005 (8) , 24 (16) :2913～2916.
[3]张正宇, 赵根等.三峡三期碾压混凝土围堰拆除爆破研究[J].工程爆破, 2003 (3) , 9 (1) :1～8.
[4]唐海, 李海波等.地形地貌对爆破振动波传播的影响实验研究[J].岩石力学与工程学报, 2007 (9) , 26 (9) :1818～1823.
[5]郭涛, 徐全军等.钢筋混凝土地下人防工事爆破拆除[J].工程爆破, 2004 (3) , 10 (1) :35～37.
[6]姚熊亮.舰船结构振动冲击与噪声[M].国防工业出版社.2007 (2) .
[7]Jones H.Plastic Behavior on Ship Structures.Trans.SNAME.1976, 84:115～145.
[8]Finite Element Modeling of the RAH-66 ComancheHelicopter Tailcone Section Using PATRAN and DYTRAN.ADA392096, 2001.
双向流固耦合 篇6
在港区压载、洗舱含油污水处理回收装置中, 气携式液液旋流器是工艺流程中预处理工段的核心设备, 担负着将油污分层回收的功能。该装置集气浮罐和液液旋流器两种分离设备的优点于一体, 其油水分离效率可达90%以上, 且占地面积小, 同时也是未来港区含油污水处理工艺小型化、船载化发展方向重点研发的装备[1,2]。
气携式液液旋流器不是将气浮罐和旋流器的构型和功能简单地叠加, 而是要从根本上进行全新的构型设计, 以实现气、液、油三相间的分离, 该设备的详细构型图如图1所示[3]。其中注气腔与旋流器大椎段之间的微孔旋流管 (套筒) 是该设备最核心的部件, 空压机进气通过该部件的微孔进入旋流器锥段, 实现气浮功能, 同时利用注气腔内一定的正压, 确保锥段内油污水正常旋流分离, 不流出到注气腔形成倒灌[4,5]。微孔旋流套筒的构型设计和微孔孔径参数, 对设备整体分离效率、处理性能有决定性影响, 是工程技术人员研发的重点。
利用数值模拟分析, 如CFD软件Fluent建立各介质在分离设备的场中分布和运动轨迹, 分析各个构型参数和水利条件对分离效率的影响, 优选构形参数, 是近年来化工传质分离设备工程设计的一种趋势[6]。CFD模拟离散相在流场中的传质分离计算, 工程上一般采用一维计算方法, 且采用经验公式给定诸多参数[7]。油水两相旋流分离的数值模拟, 流涉、水油两相相互作用、油滴间相互作用及油滴流动的随机性等的行为, 多采用欧拉—拉格朗日方法, 该方法在Fluent中表现为对旋流器内的离散相进行模拟[8]。气携式液液旋流器由于其引入了气浮相, 旋流器内流场结构呈现三维不对称性, 由于气体上升流动较为复杂, 还伴有气泡破裂吸纳油滴的过程, 所以采用传统的计算方法无法准确预测上述复杂过程。近年来, 随着计算流体力学和数值传质分离学的发展, 将二者结合起来进行一体化求解的耦合计算方法为多相间的传质分离问题的解决开辟了新的方向, 并且已经得到了良好的应用[9]。
本文拟将这种方法应用到气携式液液分离器微孔旋流套筒内的多离散相分析中, 对整个套筒构型进行模拟研究, 其中考虑了固体壁面作用, 以及流场内部和油水传质过程。由于油滴粒径较小、浓度低, 水相流场与油滴运动是相互影响的, 因此在计算中只需给定整个计算域外边界的边界条件。同时对油滴采用相间耦合随机轨道模式进行预测, 考虑油滴运动与水体湍流的相互作用, 利用油滴的运动轨迹来计算整体分离效率[10]。该类方法基于传统经验公式, 且各相间计算不再孤立进行, 流场与固壁的模拟整合为耦合求解, 能更精确地反映流场特性, 是目前最流行的模拟方法。
2 研究设备参数和模型描述
2.1 工艺流程描述
本实验所采用的轻质分散相液液分离的气携式液液水力旋流器, 它将旋流与气浮原理结合起来, 克服了常规旋流器液滴易破碎和小油滴去除效果差的缺点。它的构型主要包括一个小锥段、一个带有微孔的旋流套筒和注气腔, 并通过旋流腔入口垫板和旋流腔上盖板将三者联合形成密闭空间。通过螺栓、法兰以及密封垫片构建的Ⅰ、Ⅱ、Ⅲ、Ⅳ4个腔体室, 原水腔、注气腔之间开有2个平行的切向入口, 油污水从此进入设备。油污水在旋流管内形成高速旋转的涡流, 油水之间密度有差从而在离心力作用下, 密度较小的油滴分散相向内运动, 形成泊核, 并最终螺旋向上通过溢流出水管排出。而水体密度较大且是连续相, 逐渐在外侧向下运动, 由底流出水管排出。同时空压机将压缩气体注入注气室, 并经过旋流管壁上的微孔进入旋流管内形成微气泡, 它们可以携带细小油滴发生气浮作用, 有效去除油滴, 提高分离效率。
2.2 旋流器构形参数
实验所用的旋流器, 微孔旋流管 (即核心部件, 图1中部件20) 由具有特殊微孔结构的微孔管加工而成, 微孔管是标准型号制品, 由聚乙烯粉末烧结加工而成。微孔管在注气后可产生大量均匀细小的气泡, 从而起到气浮效果。按工程经验, 实验设备采用的微孔直径范围在5~10μm之间, 大锥角α为20°, 小锥段锥角θ为1.5°时, 油水分离效果最佳[11]。
2.3 模型和边界条件
由于旋流器结构复杂, 某些局部细节在不影响计算结果的前提下为计算方便需进行适当简化, 如法兰端面密封圈突出部;旋流器流场网格布设模型尽量简化, 分为大锥段、小锥段、注气腔、溢流口、出水口等多个圆柱 (圆锥) 体;为了保证网格生成的质量, 在保证流量相等的前提下压缩空气入口由圆型改为方形便于计算。对流动变化大的局部区域进行适当加密, 以保证计算精度。分别对流体区, 不同介质所对应的固体区生成各自的体网格, 整个网格以六面体网格为主。网格数量约为100万个。划分后的网格如图2所示。
计算中认为流体在旋流器中的流动是三维不可压的粘性湍流流动, 湍流模型选用标准模型。求解方法为“Simple-C”方法, 计算采用CFD常用Fluent计算软件完成。计算中主要考虑了液液传质分离和气液分压平衡两种物理分离形式, 并耦合嵌入“ω-β”模型, 以体现气泡将油滴颗粒化并上升带走这一过程。“ω-β”湍流模型是CFD软件中一种标准模型, 它比标准模型修正了湍动黏度, 考虑了平均流动中的旋转及旋转流动, 很适合处理高应变率及流线弯曲程度较大的流动。它的边界条件设定如下。
进口边界:采用速度入口, 根据已知流量以及入口直径, 直接得到气相入口速度, 其余相应的湍流参数一并计算得出。根据油滴循环量确定固相油滴浓度, 其入口速度具体操作参数为:两侧入口流速为5m/s, 轴向入口流速为2.5m/s。
出口边界:出口按照湍流流动充分发展处理, 采用自由出流“outflow”。
固壁边界:壁面为无滑移边界条件, 选定默认壁面粗糙度, 其值为0.5。采用标准壁面函数法处理边界湍流。设定旋流器最外层壁面为绝热壁。水—油、气—水两个接触壁面, 均勾选“Coupled”耦合壁面边界条件。设计水温为15℃, 气体温度为25℃, 压力1.0MPa.G, 水力停留时间为可调参数, 将外边界设置为等温壁。为节省计算时间, 可以先求解流动初场, 此时的水—油壁面设置为绝热壁。当计算达到初步收敛时, 加入气—水耦合壁面, 最终达到收敛。这样做的好处是在计算网格量很大的情况下, 能够提高计算效率。如果开始就加入气—水耦合求解, 需要更加长的时间才能达到收敛, 计算效率会明显降低。
3 模型计算结果分析
3.1 结算结果分析
图3给出了数值计算迭代的步数与残差的关系和旋流器内流体迹线图。在设定各个相关参数, 残差数量级为10-3时, 程序计算到500步, 质量流率的监控曲线已趋近平衡, 说明收敛。一般而言各个参量的残差随着迭代步数的增加逐渐减小, 且变化曲线光滑, 即可认为模型的选择和初始操作参数设置合理。迭代收敛的步数还与残差等级精度有关, 但主要的评判标准是看进出口液体质量流率是否趋于稳定。
模型中液体在旋流器内是以螺旋流形式运动, 且从固壁向中心形成两种运动方向相反的螺旋线。外部液体向下流向底流口形成外旋流, 内部液体向上流向溢流口形成内旋流, 实现旋流分离。油相液滴与连续相液体介质在给定速度时由环向入口匀速进入旋流器, 油水两相的密度差, 决定了在离心力作用下较轻的油相向旋流器中心移动, 在轴线附近形成分散相的核心, 并向上从溢流口排出。
3.2 油相浓度分布
图4为水力旋流器内部沿轴向纵断面的油滴等浓度分布图 (以分流比F=10%为例) , 在已有模型下进行预测处理, 当油水混合液由入口进入旋流器后可以在其内部实现分离。初始阶段在器壁附近油相浓度几乎为零, 随着向中心移动油相浓度逐渐增高, 在溢流口附近和大小锥段中心处油相浓度最大。在底流出口处, 油相体积浓度减至0.1%~0.3%左右, 分离率85%。此外锥段下部靠近底流口处的迹线形状基本呈直线流束, 可说明底流口附近区域基本无螺旋流存在, 也说明底流口处没有分离作用。将分流比操作参数加大, 底流口附近的直线流束起始点随之逐渐远离底流口, 因此表明旋流器的油水分离主要是在旋流腔、大锥段和小锥段实现, 这符合旋流分离理论分析。当流体进入大锥段加速后, 由于流速急剧增大, 使流体形成高速旋转的螺旋流并产生较高的离心力, 分散油相在离心力的作用下分离出来聚集到中心处, 形成油芯并反向流向溢流低压区, 进而从溢流口排出。
此外在分流比为F=6%、F=10%、F=15%条件下, 旋流器纵断面油相浓度分布得出的分流比, 它对分离效率的影响很大。当F=10%时, 分离效果最好, 底流含油体积分数可降至0.5%以下, 分离效率在85%以上, 实际工业用旋流器分流比也大致在10%~12%这一区间, 也侧面证实了这一观点。F=6%时, 分离效果不错, 但由于分流比过小, 油滴大部分都从底流流出, 使底流含油浓度加大。F=15%时, 分离效果不好, 由于加大了分流比, 旋流器内部的流场发生了变化, 此时流体的迹线在接近小锥段尾部时几乎不再成螺旋形, 故而效果较差。
4 速度场分布特性
4.1 切向速度流场分布
在速度场方面, 对于切向速度、轴向速度和径向速度三个速度分量, 切向速度是它们最重要的单独分量, 它表明油滴受离心力大小, 且它的实际测量参数在三个分量中最大。因为它容易测定, 经多长期研究, 多数学者认为切向速度场可以分为内部的准强制涡和外部的准自由涡。本次模拟实验中就是要将旋流器内的双涡结构表现出来, 即中心处的强制涡和外部的自由涡。
图5中切向速度在径向方向上先增大后减小, 然后速度反向增大又减至零。同时壁面处速度不为零, 这与假定的边界与固壁条件相符。分析大锥段A (Z=0.01m) 和小锥段B (Z=0.04m) 处截面的切向速度矢量图, 看出各个截面的切向速度分布大体相似, 随着流体向下流动切向速度的数值不断衰减, 沿径向的读数也有相同规律, 都可表明中心处的强制涡和外部自由涡的存在, 以及反向速度证明了内螺旋流的存在。
4.2 轴向速度场
轴向速度场分两部分, 即占主要趋势的走向溢流的内旋流, 以及走向底流的外旋流, 它内部存在一个零轴向速度包络面。该包络面将流场分成内、外两个旋流区。外旋流区轴向速度指向底流, 在旋流器器壁附近轴向速度达到最大值, 随着半径减小轴向速度亦减小, 这一区域的流体边旋转边向底流流动, 最终从底流口排出。图6是x=0面的轴向速度矢量图, 图中可以看到内、外两个旋流区, 曲线清晰地反映出了油水分离情况。
4.3 径向速度场
在水力旋流器内部速度场中的径向速度是三个速度分量中数量级最小, 常见的LDA系统无法实现对径向速度的直接测量, 故而它较难测定。径向速度影响被分离介质径向迁移, 对旋流器内径向速度分布规律具有重要意义。
本实验的径向速度结果在图7中表示, 可以看出分散相油滴沿径向的滑移速度指向中心, 由器壁向中心处逐渐增大, 最大滑移速度出现在距离中心1cm左右, 此处离心加速度最大。这是由于离心分离过程导致油水两相的径向速度差异。而图中在中心处滑移速度几乎为零, 这是因为在中心处油滴集合成油芯, 基本上不存在相对滑移速度。
5 结论
从以上的模拟分析我们可以得出这样的认识:原型机对油水两相的分离效果最高可达85%以上, 进一步提高其分离效率主要是看微孔注气的气浮作用。如果对原型机进行注气, 由于空气密度为油、气、水三相中最小, 因此当气体会有相当一部分进入旋流器的核心处, 最终由溢流口排出, 如果注气量及气泡的粒径大小适当, 会更好地携带油滴并加速油滴的运移过程, 同时也将更多的油滴带入核心处。但如果注气量过大或气泡粒径过大, 它向中心运移的速度会明显加快, 其携带油滴的能力也会受到影响, 并且许多个这样的气泡会占据旋流器核心的大部分空间, 不利于油从溢流口的排出。
实验证实在设计气携式液液旋流器时, 在旋流器入口处注入一定的气体, 并将溢流管直径加大;在单体旋流器锥段开注气口, 并将溢流管直径加大;旋流腔用微孔材料代替刚性管, 并加上钢套使此处形成环形腔, 加大溢流管直径三种方式, 提升旋流器油水分离效果。同时旋流器的操作参数, 分流比在10%~12%区间, 将获得较好的油水分离效果, 且注气量不宜过大, 微孔的孔径也需要较为合适才可以形成较好的气浮分离效果。
摘要：指出了气携式液液旋流器是港区小型化、船载化洗舱、压仓含油废水预处理的关键设备。采用流固耦合数值模拟方法再现压缩机进气、油水旋流与设备本体的整体气浮—旋流—分离耦合计算模型, 对其核心部件微孔旋流套管的构型、孔径、内外压差及腔内流场分布进行了三维数值模拟, 得到了微孔旋流腔内汽、水、油三相的流场分布, 合理确定了注气腔气—水平衡分压、溢流比和微孔孔径, 并模拟了压缩机、潜水泵正常工况下进气、进水、出水、油污溢流的流量、流速波动范围, 为该设备的一体化构型设计和加工选材提供了参考数据。
双向流固耦合 篇7
关键词：月球车驱动轮,月壤,流固耦合,ABAQUS
1 引言
美国“机遇号” (Opportunity) 、“勇气号” (Spirit) 火星车以及苏联的Lunakhod月球车的成功向世人展示了自主移动机器人在星球探测中的重大作用。月球车也因此成为了世界上新一轮探月热潮中的亮点。
月球车是各种探测仪器的载体, 其基本功能是具有在未知的复杂路面行走的能力, 以满足科学探测考察的需要 [1]。月球表面的土壤通常较为松软, 在松软月壤中, 月球车车轮容易发生滑转, 从而造成牵引性能下降。正是由于月壤可通过能力差, 因此, 深入研究在月面特殊环境下的月球车牵引性能, 对于月球车的性能评估及面向月球的高通过性行走机构研制具有重要意义。2008年李建桥等应用离散元软件PFC2D对月球车轮在月壤上的牵引通过性能进行了模拟研究 [2];陶建国等建立了月球车刚性车轮在松软土壤上前进和转向的轮地作用力学模型, 并设计了一套车轮性能参数测试系统, 对所设计的一种刚性车轮在松软土壤上运动的力学特性进行了初步的实验测试 [3];2010年杨艳静等建立了模拟月壤的仿真模型, 对月球车刚性车轮和模拟月壤的相互作用进行了仿真分析和试验验证, 结果发现挂钩牵引力的计算结果在低滑转率下和试验结果复合较好, 而在高滑转率下, 计算结果偏大 [4]。
在以往研究中, 大滑转率下驱动轮在松散的月壤上通过时导致的大变形问题一致没有得到很好的解决, 所以我们采用流固耦合方法建立轮壤作用模型, 其中月壤为欧拉体, 模拟计算结果表明, 本流固耦合轮壤作用模型是有效的, 可以应用于月球车轮的设计、优化和性能评价等方面, 并对火星车等其他移动
机器人以及地面车辆在松软土壤中的轮地相互作用也具有参考价值。
2 轮壤有限元模型的建立
2.1 驱动轮有限元模型的建立
由于月球车车轮多为金属车轮, 所以车轮在模拟月壤表面行走时, 模拟月壤的变形远大于车轮的变形;但是车轮的变形和应力分布不是我们关注的主要问题。由于显式积分计算本身时间步长非常小, 加上大滑转率下的大变形问题, 导致计算周期特别长, 因此这里将车轮视为刚性体。建立的车轮外轮廓模型, 采用壳单元来划分网格。车轮与轮轴中心点之间采用刚性耦合, 车轮与轮轴中心点一起运动, 以此简化车轮模型。如图1所示。车轮直径为350mm, 轮宽150mm。
车轮材料为铝合金, 密度为0.0027g/mm3, 弹性模量为70GPa, 泊松比为0.3。
2.2 欧拉模型的建立
2.2.1欧拉域的选择
在传统的拉格朗日分析中, 节点是由材料确定的, 材料变形则单元也变形。拉格朗日单元通常是100%的单一材料, 因此材料边界和单元边界是一致的。当物质发生大变形、断裂、分离时, 有限元法中的网格发生畸变, 使计算中断, 数值模拟难以处理。相对地, 在欧拉分析中, 节点是空间固定的, 单元不会发生变形, 而材料在单元间流动。欧拉单元可能不会是100%的充满材料, 很多情况下可能是部分填充甚至是空的。因此, 欧拉材料的边界必须在每个增量步中进行计算, 通常和单元边界并不一致。欧拉网格通常是由简单的矩形单元组成, 为材料提供流动和变形的空间。一旦欧拉材料移动到欧拉网格以外, 它就不再参与到欧拉分析中了。
我们根据滑转率、驱动轮前进速度、沉陷量来确定月壤模型大小。月壤模型长1500mm、宽460mm、高100mm。欧拉域 要大于月 壤模型 , 欧拉域长1500mm、宽460mm、高150mm。考虑到节省计算资源, 欧拉网格在高度方向的上半部分网格划分较密, 下半部分网格划分较粗;欧拉域上表面与车轮直接接触部分网格较密, 而远离车轮部位网格稀疏。欧拉域为六面体网格, 单元属性为EC3D8R, 如图2所示。
2.2.2 月壤本构模型的选择
由于月球的低重力, 所以在力学分析时可认为月壤是摩擦-粘性土壤。
在实际工程问题分析中, 大多采用等向强化弹塑性本构关系描述月壤材料。到目前为止, 适用于该种月壤建模的本构模型主要有Mohr-Coulomb模型和修正的Drucker-Prager模型, 前者在描述岩土强度特性时的有效性几乎不受限制, 是一种比较通用的本构模型;而后者则在数值计算方面具有更高的效率 [5,6]。但后者的有效性有一定的限制, 已有人证明, 只有当土壤的内摩擦角小于22°时, 修正的Drucker-Prager模型才能较好地逼近Mohr-Coulomb模型, 也就是有效的, 而月壤的内摩擦角大于22°, 不能采用Drucker
Prager模型, 只能采用Mohr-Coulomb模型来建立月壤模型。Drucker-Prager准则表达式如下:
式中: I1为应力第一不变量; J2为应力偏量第二不变量; α、K为试验参数, 与材料性质及模型选择有关, 可以由粘聚力C和内摩擦角A确定。
2.2.3 月壤材料参数的确定
建立月壤Drucker-Prager模型所需的参数主要有 [7,8]:月壤的密度参数、弹性参数和MohrCoulomb塑性参数。其中弹性参数包括月壤的杨氏模量和泊松比, 而Mohr-Coulomb塑性参数又可分为Plasticity参数和Hardening参数, 前者包括偏心率、子午线偏心率、摩擦角和膨胀角;后者包括内聚力和热膨胀塑性应变。
取月壤的内聚力和内摩擦角分别为2kP a和20°;偏心率和子午线偏心率取软件默认值;由于月壤内聚力几乎与温度无关, 故取热膨胀塑性应变为0。
通过查阅大量与月壤性质相近的粘土参数, 大致确定了其弹性模量的波动范围, 最终取月壤弹性模量为10Mpa;月壤密度取为0.0054g/mm3;由前面的表述可知, 月壤属于粘性土壤, 故取月壤的泊松比为0.39。
膨胀角反映的是介质剪切变形引起的体积改变。由于月壤剪切变形对体积改变影响不大, 月壤膨胀角取为5°。
3 驱动轮与月壤相互作用模型的建立
默认的, 所有的欧拉单元都是初始被定义为“虚”材料。在定义初始条件时, 欧拉单元中月壤部分定义为充满月壤材料, 即月壤域内材料的体积分数为1, 剩余材料被忽略。在分析中, 材料根据定义的载荷发生变形, 因而体积分数被重新计算。
为了正确地模拟车轮和月壤的相互作用, 需要在有限元模型上施加相应的边界条件。对整个分析模型建立显式通用接触, 其支持欧拉体和拉格朗日体之间的接触关系。由于定义了欧拉体, 这个通用接触自动地被扩展为流固耦合接触, 即欧拉-拉格朗日接触。欧拉域底面约束速度v1, v2, v3均为零;欧拉域侧面约束为垂直于X轴的侧面约束v1, 其值为零, 垂直于Y轴的侧面约束v2为零;驱动轮轮心约束U2, URl, UR3均为零。对整个模型施加重力 (G) , 在驱动轮轮心施加驱动轮承载 (F) 及驱动轮速度 (线速度沿X方向和角速度沿Y方向) 等来模拟月球车轮在月壤上的行走过程。
设定计算类型、输出变量、计算时间等得到驱动轮与月壤相互作用的有限元模型, 见图3。
4 数值模拟结果分析
在驱动轮承受3kg载荷, 线速度0.02mm/ms, 滑转率为0.2下进行模拟计算。
图4为模拟驱动轮滚过后的模拟的月壤痕迹。
图5为挂钩牵引力随时间变化曲线;图6为驱动力矩随时间变化曲线;图7为沉陷量随时间变化曲线。
可以看出挂钩牵引力、驱动力矩、沉陷量在行走一段时间之后才能达到较为稳定的状态, 且有明显波动。波动现象是由于轮刺是间隔的, 在滚动中轮下土壤受力也有周期性变化, 事实上轮轴产生不同程度的起伏, 轮子的挂钩牵引力、驱动力矩、沉陷量也就有了周期性变化。
由这些图可以直观地观察月壤的变形情况, 并方便地查取相关变形数据。不难发现, 该月壤模型的变形情况和预期的结果是相符的, 与试验结果吻合较好, 这在一定程度上表明了该流固耦合方法建立的月球车轮-壤作用模型的有效性。
5 结论
仿真计算某一车轮在月壌上的滚动结果, 表明了流固耦合模型的有效性。由于月壤参数并不是一成不变的, 仿真分析时, 应根据不同的地点和不同的应用场合适当调整月壤模型的尺寸和相关材料参数, 以使其更加符合真实情况。
该仿真方法可以应用于月球车轮的设计、优化和性能评价等方面, 并对火星车等其他移动机器人以及地面车辆在松软土壤中的轮地相互作用的特性分析具一定的有参考价值。
参考文献
[1]邹猛, 李建桥, 李因武, 等.刚性轮-月壤相互作用预测模型及试验研究[J].农业工程学报, 2007, 23 (12) :119-123.
[2]李建桥, 邹猛, 贾阳, 等.月球车轮与月壤相互作用动力学模拟[J].农业机械学报, 2008, 39 (4) :1-4.
[3]陶建国, 胡明, 高海波, 等, 月球车刚性车轮与土壤相互作用的力学模型与测试[J].空间科学学报, 2008, 28 (4) :340-344.
[4]杨艳静, 向树红.月球车刚性车轮与模拟月壤相互作用有限元仿真和试验验证[J].强度与环境, 2010, 37 (1) :47-52.
[5]李杰, 庄继德, 赵旗.车辆地面力学弹塑性本构关系的研究[J].吉林工业大学学报, 1999, 29 (2) :1-7.
[6]郑永春, 欧阳自远, 王世杰, 等.月壤的物理和机械性质[J].矿物岩石, 2004, 24 (4) :14-19;
[7]朱向荣, 王金昌.ABAQUS软件中部分土模型简介及其工程应用[J].岩土力学, 2004, 25 (2) :144-148.
双向流固耦合 篇8
目前, 在工程学科中, 特别是流体动力学领域中, 越来越多的实际问题需要进行耦合场的模拟分析, 例如水轮机的叶片在水流中的变形情况, 风机的叶片在风场中的变形情况等。
因为受软件开发水平的限制, 很多软件只能完成单一物理场的模拟, 而不能够实现多物理场的耦合。为了能够实现对某一物理模型的多场耦合分析, 国内的一些高校和研究机构通常对相关软件进行二次开发, 例如对Ansys[1]和FLUENT进行程序接口的二次开发, 来解决不同软件之间的数据交换问题, 但这种方法不仅要求开发者具有相当高的编程水平, 同时也需要耗费大量的时间, 而且这些机构开发出来的程序也往往只适用于他们自己所研究的领域, 所以在推广上具有很大的局限性。
虽然实现流固耦合分析的软件很多, 方法也不少, 但由于受方方面面因素的制约, 国内在这方面的资料却很稀缺, 本文以一螺旋S型风力机的叶片在风场中旋转时的某一瞬时状态为例, 来介绍如何通过Ansys Workbench[2] (ANSYS-CFX) 来实现单向流固耦合 (FSI) 分析的方法。
1三维建模
目前比较流行的三维建模工具有Pro/E, Solidworks, CATIA以及Gambit, 这些三维建模软件其实原理都是相同的, 各有利弊。本文介绍的螺旋S型风轮[3]的三维模型是用Gambit依据其原型 (见图1) 并经过数值模拟优化后的结果而建立的 (见图2) 。Gambit所建的三维模型完全尊重风轮的原型, 这样的优点是, 模拟更加接近真实情况。
在三维建模软件中完成建模后, 有两种方法可以将模型导入Ansys Workbench中, 一种是直接导入, 比如Solidworks软件就可以直接和Ansys Workbench相连接, 实现数据共享, 即在Ansys Workbench的界面中可以对Solidworks所建的模型进行相关的操作;另一种方法是先用三维建模软件把模型保存为一定的文件格式, 比如用CATIA可以保存为*.catia文件格式, 用Gambit可以保存为*.sat文件格式, 然后再将保存后的文件导入Ansys Workbench。
由于Ansys Workbench软件是ANSYS公司2007年新推的ANSYS11.0版本所含带的, 所以与早期的ANSYS版本相比, 模型导入后, 一般不会发生线条开裂、形状变形等问题。
2网格划分
Ansys Workbench里面划分网格都是系统默认的单元, 不需要操作者去选择, 只要控制单元大小或者分网方法就行了, 也可以细化, 这比Ansys的经典界面方便很多, 而且复杂模型网格质量也比较好。对于那些熟悉Ansys APDL的人可以直接在Workbench插入Command, 这些Command其实就是经典Ansys的APDL, 但需要注意的是, 由于Workbench是调用Ansys的求解器, 调到Ansys的模型只有节点和单元, 所以用Command只能处理有限元的模型, 即只有对节点和单元操作的APDL才有效, 而对实体操作的APDL是没有作用的。
在划分网格前需要进行一项特别重要的工作:就是需要为三维模型建立其所处的流场 (流场的介质可以为水, 也可以是空气等) , 流场可以是长方形的、正方形的或者是管道等等。
在FLUENT中划分网格时, 往往需要对叶片周围流场的网格进行加密, 以获得更精确的结果, 但在CFX中, 则不需要进行这样的操作, 这是由CFX自身所带流体场的性质决定的。
对于叶片的网格, 需要对其局部的网格进行细化, 这样算出来的结果更精确。
叶片和流场的网格一般是一起被划分的, 但可以分开进行尺寸设置, 本文所举的例子是风场, 所以其尺寸远远大于叶片的尺寸, 这样设置的原因也是为了能够更加真实地模拟风场。
3在CFX中进行流场的计算
在CFX里的流场计算原理和在FLUENT里面的原理是一样的, 但也有如下一些区别。
(1) CFX的优点是物理模型丰富、功能强大、基于有限元的有限体积离散方法, 精度比较高, 收敛的速度较快, 但是计算速度慢。
(2) CFX比较适合求解高速流体以及多物理场问题, 尤其是气固、热方面的问题。
(3) CFX计算所占用的内存要比FLUENT大。
(4) FLUENT物理模型比较缺乏, 很多问题没有对应的模型, 比如多相流中每相不能是多组分, 对于湿空气和其他流体组成的多相流就不能算了。而且FLUENT的前处理器格式封闭, 只能适合于FLUENT, CFX的前处理器icem输出格式丰富。
(5) 在人工干预求解方面, FLUENT要比CFX更好一点, 但在算法方面, 普遍认为CFX要更先进一些。
(6) FLUENT计算速度快, 但收敛比较慢。
3.1 流场的设置
用CFX对流场进行数值模拟, 有多种方程模型可以选择, 本例所选的是k-ε方程模型。
流场的设置主要是对气体的性质进行设置, 例如温度、密度、传热系数等等。
3.2 边界条件的设置
这里与FLUENT计算模型有区别的是, 在FLUENT里面叶片部分在风场模型中是挖空的, 而在CFX里则不需要挖空, 把叶片三维模型直接放到风场中进行计算就可以了。
3.2.1 风场部分设置的边界条件
(1) 风速。风速加在进风面, 本例速度V=10 m/s。
(2) 风道的四周面的属性设置成壁, 这个和FLUENT里面是一样的, 即为Wall。
(3) 出口因为是与空气相同, 通常这里的压强都设置成相对大气的压强, 其值为0。
3.2.2 叶片部分设置的边界条件
(1) 本例中, 叶片在风场中模拟的转速为40 r/min, 约为4.2 rad/s。
(2) 叶片的固体属性。在CFX里, 叶片也被设置成壁, 这样就更能模拟出其固体的属性了。
从算得结果的旋转速度矢量图可以看出, 叶片是以轴为中心的逆时针转动, 叶片尖的最大旋转速度矢量值和FLUENT中的计算模型是一致的。
3.3 CFX的计算结果
在计算所得的结果中, 可以发现叶片迎风面的压强值比背风面的要大很多, 负压区都出现在背风面。同时也发现, 叶片附近的风速发生了变化, 最大风速值达到近14 m/s, 超过了设置的最大风速10 m/s。而且风在绕过叶片之后产生了漩涡, 运动变得异常复杂。
在获得CFX流场数值分析的结果后, 我们将叶片表面的压强分布以及最大压强值和FLUENT计算的结果进行了对比, 结果发现如下。
(1) 压强分布基本吻合。
(2) 压强峰值接近 (CFX中最大压强约80 Pa, FLUENT中最大压强约90 Pa) 。
从而也从另一侧面验证了计算结果的可靠性。
4在Ansys中进行静力分析
在Ansys里对其结构[4]进行静力分析的步骤大体如下。
(1) 定义单元属性。
(2) 定义材料性质, 比如材料的弹性模量E、密度ρ、泊松比γ。
(3) 划分网格。
(4) 添加边界条件, 即约束条件, 这里需要添加的约束条件有重力、旋转速度、位移和从CFX导进的瞬态压强。
(5) 结果分析。
从计算结果中, 可以发现在各种荷载的综合作用下, 叶片最容易变形, 变形最大的部位为叶片下端的尖角处, 这样的分析结果可以让我们对叶片结构的设计进行优化, 在后续的优化过程中对叶片加了一些筋板, 不仅提高了风机风能的转化效率, 同时也起到稳定加固叶片结构的作用。而优化方案的基础, 就是这里所进行的结构分析所提供的数据。
而在叶片的应力分布图上, 一方面可以读出最大应力值, 来分析材料的结构是否处在安全的范围之内, 此例中应力最大值大约为24 MPa, 而低碳钢的安全强度一般在200 MPa左右, 所以从材料的角度来看, 如果使用低碳钢, 那么可以认为材料符合强度要求;另一方面也可以看出结构什么部位的应力最大, 从而对该部位进行更加详细的分析, 进而对其进行改进、优化, 以使得其在结构上获得最优。
此外, 在Ansys Workbench中, 还可以进行模态分析 (即振动分析) , 对叶片的振行进行分析, 其分析方法与传统的Ansys是一样的。
5总 结
本方法解决了一个实现流固耦合[5]的最大难题, 就是如何把分布不均的瞬态压强值导进Ansys模块进行静力分析, 使其模拟的变形状态和应力值更加精确, 更加接近真实情况。整个操作流程可以用图3来表示。
此分析方法可以被应用于所有的单向流固耦合问题的分析, 比如水泵、压缩机、风扇、吹风机、涡轮、膨胀器、涡轮增压器和鼓风机等等。其在应用上具有普遍性, 操作上具有简易性。
参考文献
[1] Moaveni, S.有限元分析-ANSYS理论与应用[M].北京:电子工业出版社, 2008.
[2]李兵, 陈雪峰.ANSYS Workbench设计、仿真与优化[M].北京:清华大学出版社, 2008.
[3]何宗敬.Darrieus-Savonius组合风轮气动性能的研究[J].太阳能学报, 1992, 14 (4) :239-244.
[4]朱英辉, 刘志璋.升、阻力混合型垂直轴风力机结构与性能的分析[J].华东电力, 2008, 36 (7) :99-101.