数值模拟实验（精选12篇）
数值模拟实验 篇1
1 概述
气浮工作台通过压缩空气在重物和支撑表面间形成空气薄膜,产生法向的推力支承重物,其优点主要是克服了传统工作台响应滞后大、行程短、存在摩擦等不足,可在大规模集成电路制造、微机械加工、二维精细测量上广泛应用,代表着超精密、长行程的工作台的发展方向[1,2]。
雷诺方程可以描述流体在气膜间隙内的流动,但其形式较为复杂,多数情况下只能通过数值计算来求近似解。在流体内部流场的数值计算中,常用的方法有有限差分法、有限元法、有限边界元法等[3]。运用有限差分法推导的方程可保证满足守恒定律,且可获得较高的计算效率和计算精度,是目前流动问题数值计算中广泛应用的方法。
本文基于FLUENT软件,对气浮工作台内部流体的流场进行数值计算,给出了节流器的结构参数、供气压力等因素对承载力、静态刚度的影响关系,并进行气浮工作台实验,由此得到其静态性能。实验结果和理论计算之间有着较好的一致性,从而说明数值模拟的正确性和可行性。
2 气浮工作台数学模型
气浮工作台具有的承载能力、刚度等性能与气膜间隙内的气体压力分布有很大的关系。确定压力分布的方程可由气体力学的基本方程——气体运动方程、连续性方程、气体的状态方程和能量方程导出[4,5]。
气浮工作台数学模型如图1所示,依照这个模型推导出气膜间隙内气体压力分布的方程。
如图1建立坐标系,略去气体运动方程——Navier-Stokes方程中的外力项,可以给出下列方程式:
式(1)中,u,v,w是x,y,z方向气体的速度,μ是气体的动力粘度,ρ是密度。
为了说明气浮工作台各方面的性能,需要知道其承载能力(W)以及静态刚度(K)。其中承载能力是将气膜内的压力分布进行积分得到表达式为:
气浮工作台静态刚度K=,这里采用差分法代替微分法计算,精度满足工程需要,其表达式为;
3 FLUENT软件的应用
FLUENT软件是将求解区域划分成一系列有单个节点作代表的控制体积,利用控制方程与控制体积的积分得出具有守恒特性的离散方程[6,7],方程系数物理意义明确,是目前流场和热传导问题的数值计算中应用最广泛的方法。
FLUENT和其他有限元分析软件一样,包括前处理器、求解器和后处理器,FLUENT中各模块之间的关联如图2所示。
本文采用FLUENT软件进行模拟计算,需先用GAMBIT软件进行前处理,最后进行计算结果的后处理,将结果以直观的形式输出[8,9]。FLUENT软件的使用步骤如图3所示。
4 数值仿真结果及分析
设润滑气体为常温空气,动力粘度μ=1.833×10-10kgs/cm2,绝热指数k=1.4,密度ρ=1.226×10-16kg/cm3,气体常数R=29.27 m/K,常温T=288 K,喷嘴流量系数c0=0.85。在此前提下,将对节流孔径的大小、供气压力的大小等因素对静态性能的影响关系分别进行分析。
4.1 计算承载能力
供气压力分别是0.5 MPa,0.6 MPa和0.7 MPa时,相对的节流孔直径是0.3 mm,不同气膜厚度下气浮工作台承载能力如图4。
4.2 刚度的计算
依据气浮工作台承载能力在不同的压力下的分析得知,其约为在供气压力不一样的情况下的刚度,在节流孔直径为0.3mm,供气压力分别是0.5MPa,0.6MPa和0.7MPa时,气膜厚度变化导致其刚度变化由图5表示。
从以上结果可知:
(1)相应参数不变时,气浮工作台的承载能力随供气压力变大;
(2)相应参数不变时,供气的压力越大,刚度也就越大。
5 实验分析
气浮工作台性能实验的装置主要包括搭建气浮实验平台、涡流位移的传感器以及LMS SCA-DAS数据采集前端、涡流位移的测试元、砝码、千分表以及微动台。
实验软件部分采用LMS Test.Lab 9A中的Signature Testing模块,实验方法如图6所示,设计的气浮工作台受到砝码重力的作用,气膜的间隙产生变化,此时传感器的侧头和位置测试的距离发生改变,传感器信号随之改变,这个信号经过涡流位移的测试单元之后经过信号放大,通过数据采集前端而记录下来。接下来加载砝码,就可以取得不同信号的数据,同时利用处理方法就可以获得承载能力、气膜间隙以及静刚度的关系。
理论分析曲线和实验测试结果对比分析如图7所示。
理论分析的结果与实验测试得到数据基本一致。实验过程中,气浮工作台测试系统达到了预期的效果,实验结论与仿真结果相符,证明模型建立正确,同时也说明仿真结果的正确。
6 结论
给出的理论模型适用于气浮工作台的性能研究。应用FLUENT软件对气浮工作台的三维流场模拟计算,分析了节流孔孔径的大小、气源供气压力的大小等因素对其性能的影响关系,计算结果稳定性高、一致性好。在实验平台上进行了气体实验,并将实验测试结果与数值仿真结果进行了比较,两者有较好的吻合性,说明了该数值模拟方法具有一定的可行性,也为进一步改进气浮工作台的设计、提高其性能提供了重要依据。
摘要：运用计算流体力学FLUENT软件数值模拟,分析研究影响气浮工作台静态性能有关因素,并得到相应的变化曲线。进行气体静压实验,从实验结果可知,实验与数值模拟计算的结果有着较好的一致性,证明了该方法在气体润滑领域应用的可行性,也为进一步改进气浮工作台的结构设计,提升其性能提供了良好的依据。
关键词：气浮工作台,FLUENT,实验
参考文献
[1]刘丹,程兆谷,高海军,等.步进扫描投影光刻机工件台和掩膜台的进展[J].激光与光电子学进展,2003,40(5):14-20.
[2]刘强,张从鹏.直线电机驱动的H型气浮导轨运动平台[J].光学精密工程,2007,15(10):1540-1546.
[3]马明建,辛世界,王辉林,等.静压式平面空气轴承压力场的有限元分析[J].农业机械学报,1994(4):91-97.
[4]刘暾,刘育华,陈世杰.静压气体润滑[M].哈尔滨:哈尔滨工业大学出版社,1990.
[5]刘燕霞,朱宇姝,汪先明,等.气膜润滑及其应用[J].江西科学,2002,20(3):174-178.
[6]王福军.计算流体动力学分析-CFD软件原理与应用[M].北京:清华大学出版社,2004.
[7]韩占忠,王敬,兰小平.FLUENT流体工程仿真计算实例与应用[M].北京:北京理工大学出版社,2004.
[8]孙雅洲,卢泽生,饶河清.基于FLUENT软件的多孔质静压轴承静态特性的仿真与实验研究[J].机床与液压,2007,35(3):170-172.
[9]李兴华.空气静压轴承的静特性计算及其实际应用[J].同济大学学报,1999,27(1):69-73.
数值模拟实验 篇2
以一种新型的.Boussinesq型方程为控制方程组,采用五阶Runge-Kutta-England格式离散时间积分,采用七点差分格式离散空间导数,并通过采用恰当的出流边界条件,从而建立了非线性波传播的新型数值模拟模型.通过对均匀水深水域内波浪传播的数值模拟说明,模型能较好地模拟大水深水域和强非线性波的传播.通过设置不同的入射波参数来进行潜堤地形上波浪传播的物理模型实验,并将数值解与物理模型实验结果进行了比较.
作 者：张洪生 冯文静 王亚玲 吴中 缪国平ZHANG Hong-sheng FENG Wen-jing WANG Ya-ling WU Zhong MIAO Guo-ping 作者单位：张洪生,冯文静,缪国平,ZHANG Hong-sheng,FENG Wen-jing,MIAO Guo-ping(上海交通大学,船舶海洋与建筑工程学院,上海,30)
王亚玲,WANG Ya-ling(上海海事大学,交通运输学院,上海,35)
吴中,WU Zhong(河海大学,交通与海洋工程学院,江苏,南京,210024)
尿素热解制取氨气数值模拟 篇3
关键词：尿素溶液；热分解；NH3；数值模拟
中图分类号：X701 文献标识码：A 文章编号：1006-8937（2016）21-0172-03
1 概 述
SCR烟气脱硝技术是世界上最成熟的脱硝技术之一，其在我国已得到广泛应用。目前SCR系统还原剂氨气主要来源于液氨、氨水和尿素[1]，液氨为危险化学品，目前其安全隐患问题日益受到大家的重视，特别在一些城市热电，距离城市近，一旦出现氨泄漏将会对附近居民生活造成重大影响；而氨水由于其浓度低，从而降低了其危险性，但其耗量将大大增加，运输成本高；尿素作为一种无危险的绿色肥料，利用其热解制氨具有与液氨相同的脱硝性能，且便于运输、存储和使用，因而越来越多的城市电厂倾向于采用尿素热解制氨技术[2-4]。
尿素热解制氨技术是通过把质量浓度低于50%的尿素溶液在热解装置中雾化，蒸发后热解生成氨气。Tokmakov等[5]认为单独尿素分解的产物最有可能是NH3与HNCO。Chen等[6]通过热重分析-质谱联用技术研究了尿素的热解，发现尿素在熔点（132 ℃）之前已经开始分解，但分解量很少。Schaber等[7-8]报导了在温度高于413 K时，尿素由熔融态蒸发为气态，且当温度高于425 K时尿素热分解为NH3和HNCO。吕洪坤[9]等在一管式石英反应器上实验研究了尿素溶液的高温热分解特性以及添加Na2CO3后对相关因素的影响，保持很高的尿素有效分解率时所能达到的HNCO水解率很低，Na2CO3可以有效地促进 HNCO的水解。Gentemann等[10]在800～1 300 K的温度范围内研究了尿素溶液的热分解，讨论了温度、氧含量对NH3、CO2生成的影响。本文对尿素热分解的进行机制进行了分析并建立了一个管式尿素热解反应器模型，通过数值计算的方法分析热解温度、加热风量、液滴粒径等对尿素热解转化率的影响，为工程实际应用提供理论指导。
2 模拟计算对象
管式热解反应器示意图，如图1所示。反应器直径为1 500 mm，高温空气从入口到热解炉出口整体长度为6 000 mm。在1 000 mm轴线中心处设置一支尿素溶液喷嘴，尿素溶液通过压缩空气雾化后喷入反应器，雾化喷嘴喷射角为90 °，流量为0.035 kg/s，根据不同工况调节反应入口空气温度、流速以及尿素溶液雾化粒径。
3 模型选择及设定
反应器内尿素溶液雾化热解过程涉及到湍流流动、气液两相流、传热传质、液滴蒸发、尿素热解以及化学反应动力学等多方面，是一个极其复杂的物理、化学反应过程。本文针对流动的湍流模型选用标准湍流模型；传热模型选用P-1辐射模型。
3.1 离散相模型
对于尿素水溶液在气相中的两相流动，采用离散相模型DPM（Discrete Phase Model），即采用拉格朗日坐标系下跟踪液滴相，采用欧拉坐标系处理气相。同时由于颗粒的喷射角度还是其喷出时间都是随机分布的，认为尿素水溶液液滴在反应器内的运动负荷随机轨道模型，并耦合了两相间的相互作用，考虑了动量、质量和热量。
对于尿素水溶液液滴，采用多组分颗粒（multicomponent）模型，尿素溶液与水溶液按照50%配比，颗粒温度为50 ℃，防止尿素水溶液结晶。同时雾化模型选择solid cone类型，喷射雾化角为90 °，流量密度根据边界条件确定。
3.2 液滴蒸发模型
尿素水溶液液滴浓度为50%，尿素浓度较高，此时处理溶液雾化蒸发时，不能完全处理为纯水的喷射蒸发，应考虑尿素溶液的蒸发。由于对颗粒采用多组分颗粒（multicomponent）模型，可分别设置尿素和水的蒸发参数，其中尿素溶液设定其汽化潜热为1 398 KJ/kg，蒸发温度为420 K，沸腾点为483 K，而水溶液汽化潜热为2 263 KJ/kg，蒸发温度为284 K，沸腾点为373 K，两者混合雾化颗粒蒸发平衡采用拉乌尔定律。
对于雾化尿素液滴蒸发过程的计算，主要是通过野地的加热、蒸发、沸腾过程的模拟来分别考虑。并考虑采用准稳态模型，液滴的加热、蒸发、沸腾过程的质量和能量平衡方程如下：
3.3 尿素热分解模型
尿素热分解路径为固态/液态尿素先蒸发为气相NH2CONH2，然后气相尿素分解为NH3和HNCO，模型示意图，如图2所示[11，12]。
尿素的热分解受限于动力学参数，因此尿素会在熔融液态保持一段时间，而气态尿素在高温环境中并不稳定，根据文献[5]尿素热解最的产物最有可能是NH3与HNCO，且该反映为一个快速反应过程，此后HNCO进一步水解生成NH3和CO2，其被认为气相均反应[11]。
3.4 动力学模型
本文采用二步总包反应模型，尿素水溶液首先在高温空气中雾化蒸发生成尿素蒸汽和水蒸汽，随后经热分解分解成NH3与HNCO，而后HNCO进一步水解生成NH3和CO2，动力学方程式及参数，见表1[13，14]。
4 模拟结果与分析
4.1 尿素水溶液热解数值计算分析
尿素水溶液热解模拟结果，如图3所示，为入口空气温度为873 K，流速0.5 m/s，液滴粒径为100微米条件下尿素水溶液热解模拟结果。
从图可见，尿素水溶液喷入反应器后被迅速加热，由于尿素蒸发温度为420 K，高于水分的沸腾温度，因此水分首先从液滴中蒸发，并随着水分的蒸发，液滴表面的尿素浓度越来越高，待液滴中水分几乎蒸发完全后，此时尿素才开始逐步蒸发热解。同时随着尿素水溶液蒸发热解，空气温度逐步降低，氨气浓度不断增加，出口烟气温度降低至673 K左右。分析HNCO浓度分布可见，在高温段中尿素热解生成的HNCO随后与水蒸气发生水解反应，并在出口处基本水解完成，完全转化成氨气。
4.2 温度对热解效率的影响
入口空气流速0.5 m/s，尿素水溶液液滴粒径100 μm，分析了573～1 073 K温度区间内不同温度工况下尿素热解制氨的转化率的影响，模拟图，如图4所示。
由图4可以看出，尿素水溶液液滴在热空气流场中停留时间越长，其NH3转化率不断增加，这是因为随着尿素液滴的停留时间的增加，尿素热解越彻底，更加有利于尿素热解。
同时从图4还可以看出，随着温度的升高，尿素热解产物NH3转化率增大。前期随着温度的升高，NH3转化率大幅增加，当温度达到873 K以上时，尿素水溶液液滴在热空气中停留时间为10 s时，NH3转化率基本已经达到100%， 此后随着温度的升高尿素水溶液在更短的停留时间内就能够达到100%NH3转化率。这是由于尿素热解反应推进率常数随着温度升高而增大[15]，温度越高，尿素分解越彻底，NH3转化率越高，可见温度是尿素热解的一个关键的因素。尿素热解是一个吸热反应，温度的越高，其反应越剧烈，反应速率越快，所需的反应时间也越短，这也就解释了在温度高于873K时，随着温度的升高，尿素水溶液在更短的停留时间内就能够达到彻底转化。
4.3 空气流速对热解效率的影响
空气流速决定了进入反应器的空气流量，空气流速的变化其首先影响液滴在反应器内的停留时间，其次作为热源，空气流量大小影响着热量的供给。本模拟研究了空气温度873 K时，尿素水溶液液滴粒径100 μm条件下，空气流速在0.25～1.5 m/s区间内空气流速对热解效率的影响。
模拟分析空气流速对热解效率的影响，如图5所示。
随着空气流速的增加，在0.25～0.5 m/s区间内，尿素热解效率快速增加，而此后尿素热解效率基本不怎么变化。其原因应当是：在0.25～0.5 m/s区间内，由于空气流量低，导致其热量供给不足，尿素热解得到充足的热量，从而对热解效率影响较大；而此后随着空气流速增大，热量供给增大，且热量的增大抵消了其停留时间变短的影响，热解效率基本不变。
分析计算得到，空气流速在0.45 m/s时，空气流量达到当尿素完全热解后氨气浓度为5%，在空气温度873 K时，其热解效率基本彻底，由此可知，烟气温度达到873 K以上，在确保氨气浓度低于5%时，空气流量大小对尿素热解效率基本无影响。
4.4 雾滴颗粒粒径对热解效率的影响
在空气温度873 K，空气流速0.5 m/s工况下，分析颗粒粒径对热解效率的影响，模拟结果，如图6所示。
图6表明了随着雾滴粒径的增大，前期热解效果显著降低，尿素水溶液达到同等热解效率所需的停留时间将增大。雾滴粒径的增大加大了雾滴蒸发所需的时间，使得雾滴不能快速蒸发，同时由于蒸发吸热，在雾滴周围形成一个局部低温区，不利于尿素的热解，从而使得前期尿素热解缓慢，热解所需停留时间增大。从图中可知，当粒径大于250 μm后，现有反应器的停留时间将无法满足尿素颗粒完全热解。
5 结 语
①对尿素水溶液雾滴的蒸发热解过程进行模拟分析，发现由于水与尿素的蒸发温度不同，前期主要为水分蒸发，并随着水分的蒸发尿素水溶液浓度逐渐增大，待液滴中水分几乎蒸发完全后，此时尿素才开始逐步蒸发热解。
②热解温度对尿素热解效率有显著影响，随着温度的升高，NH3转化率热解效率增大，当温度达到873 K时，NH3转化率基本已经达到100%，此后随着温度的升高，尿素水溶液达到彻底热解的停留时间可减少，即高温度下所需的停留时间更短。
③空气流量决定了反应器内的热量供给，过低的空气流量将导致热解效率降低，同时过低的空气流量将无法保证氨气浓度低于5%，烟气温度在873 K时，在确保氨气浓度低于5%的烟气流量条件下，空气流量大小对尿素热解效率基本无影响。
④雾滴颗粒粒径的增大，使得其蒸发热解所需的时间增加，且雾滴局部温度脚底，尿素热解将受阻，要实现尿素完全热解所需停留时间将增长，反应器的尺寸将加长，设备投资增大。
参考文献：
[1] 郭伟，催宁.尿素热解制氨SCR脱硝技术在电厂的应用与优化[J].锅炉 技术，2012，43（3）：77-80.
[2] 喻小伟，李宇春，蒋娅，等.尿素热解研究及其在脱硝中的应用[J].热力 发电，2012，41（1）：1-5.
[3] 杜成章，刘诚.尿素热解和水解技术在锅炉烟气脱硝工程中的应用[J]. 华北电力技术，2010（6）：39-41.
数值模拟实验 篇4
随着现代工业的发展,直径在0.1～1.0mm之间的小孔结构广泛存在于航空航天、武器装备、汽车、电子、船舶等领域中。小孔结构材料一般是难加工材料,并且小孔的数量少则几十个、多则上万个,其加工质量、精度、效率对产品的性能、质量和成本有很大影响。传统的机械切削加工由于刀具刚性的影响及存在切削力,被加工零件易产生变形,不适合用于难加工材料的小孔批量加工。特种加工技术中的激光加工、电火花加工、电解加工可用于小孔加工。但是激光加工、电火花加工都属于热加工,加工表面存在热影响区和再铸层,影响工件的使用寿命和安全可靠性,因此需要后续光整处理,导致加工周期长、成本高[1,2]。电解加工小孔包括成形金属管电解(shaped tube electrolytic machining, STEM)加工、毛细管电解(capillary drilling, CD)加工、电射流电解(electro-stream drilling, ESD)加工和喷射电解(jet electrolytic drilling, JED)加工[3,4,5,6,7,8]。前三种方法都需要制作成形阴极,并且加工中阴极需要进给,以进入工件中。这些方法存在阴极制作困难,加工中容易出现短路而烧伤工件等缺陷。喷射电解加工小孔方法避免了上述缺陷,在加工过程中金属阴极仅对电解液束辉光充负电,使电解液束起到工具阴极的作用,打入工件,使工件按电解液束的形状阳极溶解作用而成形,从而可实现小孔高质量、高精度的加工目标[9]。
本文针对喷射电解加工小孔方法,基于电场理论,建立了二维数学模型,并利用该模型进行了数值模拟和工艺实验验证。在此基础上,通过工艺对比实验研究了喷射孔径、加工电压和电流方式对加工质量的影响,总结了喷射电解加工小孔的典型工艺特征。
1 二维数学模型
1.1理论假设
根据喷射电解加工小孔中加工间隙固定不变且无成形阴极的特点,提出如下假设:①加工过程及整个加工间隙中,电解液各向同性,即各处电导率和电流效率都相同,均分别保持为常数;②加工过程中,阳极/电解液界面和阴极/电解液界面上的电极电位均分别视为常数;③阳极电化学溶解服从法拉第定律;④加工间隙中的电场视为无源稳恒电流场。
1.2数学模型
根据上述理论假设,针对喷射电解加工小孔的电流场建立了简化模型,如图1所示。阴极喷嘴的孔直径为2R;喷嘴与工件表面距离固定,即加工间隙为L;直流稳压电源正极接工件、负极与喷嘴相连。
由电场理论可知,喷射电解加工小孔区电场的电位ϕ分布符合拉普拉斯方程,即
根据法拉第定律,图1所示的阳极溶解速度vn可表示为
式中,η为电流效率;ω为体积电化学当量;ia为电流密度。
图1中z向的边界移动速度vz可表示为
式中,α为阳极溶解速度与z向边界移动速度的夹角。
由几何关系可知:
式中,za为喷射电解加工小孔二维轮廓曲线函数;r为曲线函数变量。
电流密度与电场强度的关系如下:
式中,κe为电解液电导率。
阳极某点电场强度等于该点电位梯度,即
归纳以上关系式可得如下描述喷射电解加工小孔间隙中电位分布的数学模型方程:
阳极表面边界条件为
阴极表面边界条件为
根据上述二维数学模型,由加工电场电位的分布可求得加工区各点的电场强度,再计算出各点的电流密度,进而可得到各点的阳极溶解速度,最后可得出加工后的二维轮廓。
2 数值模拟与实验验证
2.1工艺参数
实验选用YJ63型稳压电源,输出电压为0～300V,额定电流为5A。阴极喷嘴的孔径有三种尺寸,分别为1.5mm、1.0mm和0.4mm。电解液选用NaNO3溶液,NaNO3的质量分数为18%,该电解液的电导率为12.2S/m。电解液进液压力为1.5MPa,加工间隙为2mm。工件为0.5mm厚的不锈钢片(1Cr18Ni9Ti),该材料的体积电化学当量ω为2.1×10-9m3/(A·min),电流效率η约为0.6。
2.2数值模拟结果及验证
根据上述二维数学模型,对喷射电解加工小孔过程进行了数值模拟与实验验证,图2是数值模拟结果与实验结果的对比图。实验结果是利用ADE公司的MicroXAM3D Profiler形貌仪测量所得。加工参数包括:喷嘴直径0.4mm,加工电压200V,加工时间分别为5s和30s。由图2可看出,喷射电解加工小孔数值模拟结果与实验结果比较吻合,径向误差较小,不超过0.01mm,孔在深度方向上误差较大,最大处在孔的底部中心地带,误差值达0.03mm,实验结果的二维整体轮廓略大于模拟结果。出现这一现象的原因是:在200V高电压作用下,工件阳极除了电化学溶解去除材料以外,阴阳两极之间产生的辉光放电也有助于去除材料。在喷射电解加工过程中,阴阳两极之间高速流动的电解液束周围的气压低于大气压,从而形成围绕电解液束的低压环形区域,而电解加工中阴极析出氢气,阳极溶解产生电解产物有时还包括氧气、二氧化氮等气体,因此在这个低压区域中,存在一定的气流。低压区的气流在高电压的作用下形成辉光放电,在阴阳两极之间生成放电通道,通道内的负极性粒子会撞击阳极表面,这一作用有助于去除阳极材料。另外辉光放电形成的加工区域温度梯度也有助于阳极氧化溶解。
从整体轮廓对比来判断,根据本文建立的二维数学模型可以较好地模拟实验加工结果。由于难以模拟入口处的杂散腐蚀,本文建立的数学模型适用于直径在1.0mm以下小孔的喷射电解加工的数值模拟。
3 实验结果与分析
在喷射电解加工小孔过程中,阴阳两极相对固定,没有进给。影响小孔加工质量的重要因素有:喷嘴孔径、加工电压和电流方式。因此,针对这三个重要因素对加工质量的影响进行了工艺对比实验和分析。
3.1喷嘴孔径对加工质量的影响
由喷射电解加工小孔原理可知,喷嘴孔径直接决定了喷射液束的粗细,也就决定了加工孔的大小。图3对比了采用三种不同孔径喷嘴进行喷射电解加工小孔的入口形貌。
从图3a可以看出,当喷嘴直径为1.5mm时,喷射电解加工的小孔入口周围存在明显的过切现象,即图中的黑色区域,该区域的面积是加工孔入口面积的两倍以上,另外,入口外围的杂散腐蚀非常严重。
对比图3a与图3b可以看出,当喷嘴直径由1.5mm降至1.0mm后,喷射电解加工孔入口周围的过切现象有所减轻,过切区域的面积与加工孔入口面积相当,孔的加工质量没有明显改善。当喷嘴直径降至0.4mm时,如图3c所示,小孔的入口过切现象已大为减轻,过切区域的面积仅为小孔入口面积的1/5,杂散腐蚀基本去除,孔的圆度也大为改善。
由图3可以看出,在相同的进液压力(1.5MPa)和加工间隙(2mm)下,随着喷嘴直径的减小,加工孔的质量也随着提高,在加工直径为1.0mm以下的小孔时,可以有效减少入口的杂散腐蚀,提高孔的加工精度。
3.2加工电压对加工质量的影响
喷射电解加工小孔阴阳两极之间施加的电压大小对加工效率有决定性作用,常用材料去除率Rm来表示加工效率,即
式中,Δmi为第i次实验所得的加工前后的工件质量差;n为实验重复次数;t为加工时间,每次实验加工时间相同。
图4显示了加工电压对喷射电解加工小孔效率的影响。喷嘴直径为0.4mm,加工时间统一为30s,每组实验重复5次。由图4可以看出,加工电压在120V以下时,喷射电解加工小孔的材料去除率在3mg/min以下;当加工电压增加至120V以上时,材料去除率有明显上升的趋势。这说明,当加工电压在120V以上时,喷射电解加工小孔的电流效率较高,并伴随辉光放电效应,使得加工效率得到跃升。因此,喷射电解加工小孔采用的加工电压应在120V以上,才能取得更高的加工效率。
3.3电流方式对加工质量的影响
喷射液束电解加工小孔采用的电流具有直流电流和脉冲电流两种方式。图5是利用MicroXAM 3D Profiler表面形貌仪拍摄的形貌图,对比了两种不同电流方式对喷射电解加工小孔的影响。喷嘴直径为0.4mm,直流电源电压为200V;脉冲电源电压为400V,占空比为50%,频率为2kHz。
从图5可以看出,脉冲喷射电解加工的孔入口轮廓好于直流电解加工的孔入口轮廓,孔的锥度也明显小于直流电解孔的锥度。这是因为,脉冲电流电解加工是以周期间歇供电代替传统的连续供电,使工件阳极发生周期断续的电化学阳极溶解,它可以充分利用脉冲间隙的断电间隙去极化、散热等作用,更新加工间隙的电场和流场,并有助于及时排除电解产物。因此,与直流喷射电解加工相比,脉冲喷射电解加工能够获得更高的加工精度。
4 结论
(1)基于电场理论,针对喷射电解加工小孔方法,建立了二维数学模型,并通过了数值模拟和实验验证,证实了该模型可用于模拟预测实验结果,该模型适用于直径1.0mm以下的小孔加工。
(2)通过不同孔径喷嘴进行了喷射电解加工小孔对比实验,研究了喷嘴孔径对加工质量的影响,实验结果表明在进液压力和加工间隙一定的条件下,喷嘴孔径越小,喷射电解加工小孔质量越高。该方法在加工直径1.0mm以下的小孔时,可以有效减少入口的杂散腐蚀,提高孔的加工精度。
(3)研究了加工电压和电流方式对喷射电解加工小孔质量的影响,实验表明,喷射电解加工小孔所施加的电压应在120V以上;与直流电源相比,使用脉冲电源能有效减小孔的锥度、提高加工精度。
参考文献
[1]郭文渊,王茂才,张晓兵.镍基超合金激光打孔再铸层及其控制研究进展[J].激光,2003,24(4):1-3.
[2]夏劲武,徐家文,赵建社.电火花加工表面质量的研究及进展[J].电加工与模具,2008(6):11-15.
[3]Chryssolouris G,Wollowitz M,Sun N P.Electro-chemical Hole Making[J].CIRP Annals-Manufac-turing Technology,1984,33(1):99-104.
[4]Zhu Di,Xu Huiyu.Improvement of Electrochemi-cal Machining Accuracy by Using Dual Pole Tool[J].Journal of Materials Processing Technology,2002,129(1/3):15-18.
[5]Kozak J,Rajurkar K P,Balkrishna R.Study ofElectrochemical Jet Machining Processes[J].Transactions of the ASME,Journal of Manufactur-ing Science and Engineering,1996,118(4):490-498.
[6]Li Yong,Zheng Yunfei,Yang Guang,et al.Local-ized Electrochemical Micromachining with Gap Con-trol[J].Sensors and Actuators A:Physical,2003,108(1/3):144-148.
[7]Sharma S,Jain V K,Shekhar R.ElectrochemicalDrilling of Inconel Superalloy with Acidified SodiumChloride Electrolyte[J].The International Journalof Advanced Manufacturing Technology,2002,19(7):492-500.
[8]Datta M.Microfabrication by Electrochemical MetalRemoval[J].IBM Journal of Research and Develop-ment,1998,42(5):665-669.
数值模拟实验 篇5
局部扰动对主坑道爆炸波发展的数值模拟与实验研究
在地下建筑物,如隧道、地下储油库、人防工程、地下物资仓库等里面,由主通道旁结分支通道是最常见的一种布置形式,是一种典型的复杂受限空间结构.一旦有可燃气体发生爆炸燃烧,爆炸压力波和火焰的传播将受众多因素的影响,其中局部扰动的影响是主要因素之一.本文通过实验和数值模拟的`方法研究了油气混合物在该复杂受限空间中由弱点火引起爆炸燃烧的发展过程,湍流强度经旁接分支坑道后在主通道中的变化,以及爆炸压力波和火焰经局部扰动后的变化过程;并将数值模拟结果与实验结果进行了对比和综合分析,得到了与地下受限空间安全相关的重要结论.湍流强度是复杂受限空间中可燃气体爆炸燃烧发展过程的主要影响因素之一,局部扰动将增强爆炸流场的湍流强度,加速燃烧化学反应,能量的释放量和速率大大提高.这些能量的快速加入促进了高峰值压力波的形成,火焰也被加速,爆炸从此由弱转强,出现跃升.研究结果对地下受限空间爆炸过程的进一步研究以及爆炸灾害的预防都有参考价值.
作 者：蒋新生 杜扬 袁占国 高建丰 钱海兵 JIANG Xin-sheng DU Yang YUAN Zhan-guo GAO Jian-feng QIAN Hai-bing 作者单位：解放军后勤工程学院供油系,重庆,400016刊 名：安全与环境学报 ISTIC PKU英文刊名：JOURNAL OF SAFETY AND ENVIRONMENT年，卷(期)：5(5)分类号：X932关键词：安全工程 分支通道 爆炸 湍流 燃烧
城市隧道施工的三维数值模拟分析 篇6
【关键词】城市隧道；数值模拟；地表变形；CRD工法
0.引言
目前，随着岩土工程数值方法和计算机技术的快速发展，复杂定解条件问题的处理才能成为可能，数值方法被人们广泛用来进行求解隧道施工过程中遇到的各种问题的最有效的通用方法。
目前对于CRD法在城市浅埋隧道施工中的运用及其对地表变形和衬砌受力的内在因素和机理方面的研究较少。本文结合龙岩至厦门高速铁路（简称龙厦铁路）石桥头隧道进口段工程， 对CRD法施工进行数值模拟，重点研究隧道施工引起的地表变形沉降规律、衬砌受力特点、衬砌受力与围岩变形的关系，最后通过现场监测结果验证数值模拟的合理性。
1.工程概况
龙厦高速铁路石桥头隧道位于福建省龙岩市城区，隧道整体由西北折向东南方向，为电气化双线隧道，线间距4~4.4m，行车速度120km/h，位于R=1000m的右偏曲线上，隧道净高11.2m，净宽12.8m，全长1586m，是龙厦铁路重点控制性工程。地质调查和钻探揭示，进口段DK2+450~DK3+021长571m，为Ⅴ级围岩。且该段隧道埋深浅，开挖跨度大，围岩稳定条件差，因此，控制围岩变形及地表沉降是该段隧道施工中的关键技术问题，因此本文以此段为依据进行数值模拟分析，来研究隧道施工引起的地表变形规律。
2.隧道施工数值模型的建立
由于进口段地质条件较差，选取隧道埋深约17m建立计算模型。隧道净空断面为（12.8m×11.2m），根据隧道开挖的影响范围，参考已有的计算理论和经验，模型计算范围为：中心线两侧各取50m，竖向选取地表以下70m，为了节约时间，隧道纵向取三个开挖步长即9m。在本文分析计算中，假定隧道围岩按弹塑性材料考虑，破坏准则采用摩尔-库伦准则。结合该工程，在模拟隧道施工过程中，超前管棚预加固，采用改善围岩参数的等效方法来模拟；钢拱架加喷射混凝土初期支护，采用壳体结构单元模拟；围岩、加固圈和二次衬砌，采用实体单元模拟。整个地层数值模型计算区域为100m×9m×70m，共划分网格单元13620个，节点23360个。
模型计算边界条件为：左右边界施加水平方向的约束，底部边界施加固定约束，顶部边界为自由面。模拟计算中，土层力学参数和围岩支护参数，按石桥头隧道地质勘查报告和规范选取，其值见表1。
表1模型參数表
3.CRD工法施工过程的数值模拟
建立三维网格模型→初始条件模拟分析(自重应力和顶面的均布荷载) →加固圈加固处理(超前管棚加固) →开挖左上导洞模拟计算、初期支护模拟计算→开挖右上导洞模拟计算、初期支护模拟计算→开挖左下导洞模拟计算、初期支护模拟计算→开挖右下导洞模拟计算、初期支护模拟计算→二次衬砌施加模拟计算→最后对模拟计算结果进行分析。
注意的是在开挖前必须采取一些加固措施(如超前管棚等)进行预加固地层，每部开挖之后要及时施作初期支护，最后拆除临时支护并施作二次衬砌。
通过对CRD工法在城市浅埋隧道施工进行三维数值模拟，得到各部开挖支护和二次支护的竖向位移及影响范围分布(结合tecplot10.3后期处理软件得图1)、围岩应力分布(结合tecplot10.3后期处理软件得图2)、衬砌受力分布云图(如图3)、开挖结束后的变形矢量及塑性状态分布(如图4、图5)。
第一部开挖支护第二部开挖支护
第三部开挖支护 第四部开挖支护
二次支护
图1各部开挖支护和二次支护竖向位移及影响范围
表2模拟开挖地表下沉统计/mm
由图3和表2模拟计算结果表明：
3.1采用CRD法对隧道进行施工，隧道上方出现了变形沉降，沉降影响范围随着开挖位置变化而变化，沉降影响范围首先产生并偏向于开挖一则，随着开挖的继续，影响范围也随之变化，最后影响区域几乎对称位于隧道中心线上方两侧各26m左右，最大沉降值为16.2mm，位于隧道拱顶正上方。因此在对地表沉降进行监控量测时，基准点要埋设远离隧道中心线26m以外。
3.2隧道各部开挖进行支护后，沉降虽然继续增加，但数值不大，并且主要影响范围明显缩小，由此知道在隧道开挖过程中，隧道毛洞开挖结束时，初期支护和临时支护要及时跟进，对稳定隧道开挖面极为重要。
3.3隧道不同部位的开挖对地表沉降的影响各异，左上导洞开挖沉降比例较大，因土体开挖打破了原有土体三向应力平衡，地应力第一次将重新分布，必然引起周围土体的位移和变形，由于第一步开挖非常关键，在进行开挖前有必要采取相应措施（如超前管棚加固），并及时配合临时支护，尽量控制因开挖引起的变形位移。上部导洞开挖沉降量占总沉降量多达70%，由于隧道上部是主要支撑拱顶的关键部分，而这部分土体被开挖，上部围岩此时处于临空状态，内部应力大量释放，使得拱部土体受到了较大的拉应力作用，其上部土体岩层变形较大，导致地表沉降量也较大，因此在该阶段施工过程中支护要及时，并且拱顶部位钢拱架喷混凝土的固定及连接要严密，超前小导管注浆及时跟进，以保证开挖掌子面的整体稳定。随后进行的下部开挖，地表沉降进一步增大，但相对上部而言，该部位开挖对地表沉降影响明显较小，但在施工过程中也要严格按照设计方案施工。最后二次衬砌的沉降量占总沉降4%左右，在隧道施工过程中，要求初次衬砌结束后，并且在隧道地表沉降和洞内收敛基本稳定时，再进行二次衬砌施工。
图2a围岩水平应力分布图图2b围岩竖向应力分布图
图3a衬砌竖向应力云图图3b衬砌水平应力云图
图4 变形矢量图图5 塑性状态图
由图2到图5模拟计算结果表明：
3.4衬砌周围土层应力主要为压应力，最大压应力分布在隧道顶部和底部区域，且顶部土层压应力分布与地表沉降类似为抛物线，仅在隧道两侧的局部区域有较小的拉应力。
3.5在隧道周围土体和地表荷载作用下，衬砌应力以压应力为主，局部为拉应力。在垂直方向上受力最大值达4.06MPa，主要位于隧道拱腰处，且应力比较集中，因此在施工过程中，拱腰处的钢支撑接头需重点处理，做好锁脚加强支护，防止应力集中对衬砌造成变形开裂，降低衬砌结构的整体刚度，从而影响整个围岩的稳定。隧道两腰上半部分区域承受拉应力，最大值约为1.28MPa，拉应力的产生与地层应力分布、隧道断面形状及拱顶变形等诸多因素有关，且有扩大趋势，并向腰部，顶部集中。拱顶部位受到水平应力集中影响，最大值达0.72MPa，这主要是由于拱顶部位在开挖过程中产生较大沉降有关，且该部位在开挖时作为主要受力部位承受垂直和水平两个方向的地应力。因此施工时，要保证拱顶处的固定连接，接头的连续性会对衬砌变形和结构整体沉降产生较大影响，再配合必要的加固措施以保证开挖面的整体稳定。
3.6隧道顶部受力复杂，应力集中，顶部以上区域变形较大，是隧道开挖造成地面沉降主要区域；在隧道坑壁及向外延伸线上，主要受压缩剪切，是剪切破坏的主要区域，如果长期不支护，可能会出现流动现象，发生剪切破坏，造成隧道失稳而影响地面变形过大。
4.变形监控量测
对于城市浅埋隧道施工而言，监控量测具有重要地位，它是信息化施工的基础，是确保隧道施工和周围环境安全的重要手段，是检验设计参数、地面稳定性、评价施工方法的主要依据。对于石桥头隧道工程，地表沉降使用高可靠性磁阻尼摆式补偿(补偿范围±15、、)的索佳SDL30电子水准仪，配合RAB码玻璃钢水准尺(标准差±1.0mm)进行量测；由于隧道上方多为民居低层建筑物，且建筑物周围空地较窄，因而采用通过测量建筑物基础相对沉陷的方法来确定建筑物的倾斜，因此量测仪器与地表沉降测量仪器相同，而对于一些高层建筑物，则采用可实现高精度，远距离无协作目标测距，反射片测距及棱镜测距的索佳SET230RK全站仪。由图6看出，隧道地表最终沉降的实测值与模拟值相比，误差在1.8~3.9mm之间，对于相同里程，拱顶最终沉降量一般比地表最终沉降量稍大，但最终沉降量都处于规范允许之內。地表房屋得到安全保障，隧道顺利贯通，工程工期明显缩短，不仅为之后铺轨作业提供了充裕时间，且为业主单位节省了工程经费，达到预期经济和社会效益。
图6实测地表和拱顶沉降曲线
5.结论
５．１本文所建立的三维有限差分模型，可有效的分析CRD法在城市浅埋隧道施工引起的地表沉降规律和衬砌受力特点。
５．２采用CRD法对城市浅埋隧道施工，沉降影响范围随着开挖位置变化而变化，沉降影响范围首先产生并偏向于开挖一则，随着开挖的继续，影响范围也随之变化，最后影响范围几乎对称的位于隧道中心线上方。
５．３采用CRD法对城市浅埋隧道施工，初期支护对地表沉降影响范围有明显控制效果。上部导坑开挖对整个地表变形起到关键作用，沉降量占总沉降量的70%左右。
５．４衬砌周围的土层应力主要为压应力，且顶部土层应力的分布与地表沉降类似，为抛物线，仅在隧道两侧的局部区域有较小的拉应力。而衬砌应力同样以压应力为主，局部区域为拉应力，高应力区集中在顶部和腰部。
５．５隧道顶部受力复杂，应力集中，顶部以上区域变形较大，是隧道开挖造成地表沉降主要区域。在隧道坑壁及向外延伸线上，主要受压缩剪切，是剪切破坏的主要区域，如果长期不支护，可能会出现流动现象，发生剪切破坏，造成隧道失稳而影响地面变形过大。
【参考文献】
［１］严绍洋，李亮辉，张春宇.公路隧道开挖渗流场的有限差分法分析[J].中外公路，2007，27：120-123．
［２］杨玉胜.红砂岩隧道开挖稳定性的三维数值分析[J].路基工程，2007，67-69．
数值模拟实验 篇7
关键词：稳定温度分层,CFD,分层流水槽
在广袤的自然水域中,温度分层流是一种常见的现象,经常出现在大气、海洋等自然现象中,也经常出现在能源、化工等工业过程里,这其中又以重力方向上温度不均匀引起的流体密度变化而产生的分层流最为常见。在污水排放口的附近、水池,以及江河、湖泊、海洋中都会产生温度分层。
温度分层对浮力射流的浮升运动影响非常显著,当浮力射流在温度分层环境进行数值模拟的过程中,稳定的温 度分层是 模拟精度 的重要保证。
然而,目前关于分层流的数值计算模拟研究,多数学者重点关注于密度分层与盐度分层,关于温度分层流方向的研究文献并不多见。高殿武[1]、A. Eidelman等[2]研究了热分层湍流计算模拟的数学模型; 金海生[3]、胡振红[4]等对存在温度和盐度梯度的分层流进行了数值模拟; Hayat T[5]进行了温度影响条件下Oldroyd-B黏性流体的流动分析; Rorai C[6]讨论了同时存在层流和湍流情况下的热分层运动; 姚志崇等[7]开展了分层流体中拖曳球体尾流及辐射内波试验研究。
然而,前几位学者在研究过程中,关注的重点均在于数学模型和物理模型的描述,在壁面条件上未给予足够关注。多数学者在壁面问题的处理上,采用了黏性无滑移的壁面条件,并假定没有质量或热量交换; 同时,对近壁网格节点应用壁函数方法来处理近壁黏性次层。
稳定温度分层的形成,热平衡是其前提条件,必须有持续的能量注入与交换,否则随着时间的推移, 温度分层将被破坏。以简化的海洋温度分层为例, 太阳的辐射能不断注入海洋,海水表层吸收太阳能转化为热能,在各种热交换的作用下达到动态热平衡,最终形成稳定温度分层。在真实的环境中,海洋相对静止,航行体相对运动; 而在数值模拟中,采用的是航行体相对静止,将航行体的运动转化为海洋中来流流动的处理方式。在这样的方式下,采用无滑移壁面的设定条件,流体与壁面将产生黏性作用; 采用绝热的壁面设定,将导致流体在流动过程中失去热量的交换通道并使得热量交换缺乏合理约束; 两者共同作用必将对温度场产生影响,不能维持来流的温度分层。因此,按先前的方法设定的壁面条件显然无法满足要求。
本文在借鉴前人分析思路和方法的基础上,结合温度分层的成因,重点对壁面条件的设定进行了研究,提出了以滑移壁面和给定壁面温度函数为核心,结合来流温度函数的边界条件设置方法,利用数值计算的方法成功实现了计算域内稳定温度分层的形成与保持,并通过与分层流水槽实验结果进行比较,证明了该方法的可行性。
1物理模型与数学模型
研究对象如图1所示,采用笛卡尔坐标系,坐标原点位于来流入口上部中心点,长方体计算域,- z轴方向为环境流体方向,z区间定义为0 ~ 8 000 mm,x轴为计算域宽度方向,x区间定义为 - 1 000 ~+ 1 000 mm,- y轴为计算域深度方向,y区间定义为0 ~ 1 200 mm。来流沿 - z轴方向流入并流出,并沿 - y轴方向上 有近似线 性的温度 梯度分布。
实际海洋环境中流体密度的垂向分布多是非线性的,在近表水面可近似为线性分布。使用笛卡尔直角坐标系,在计算中忽略黏性耗散项,有温度梯度分层流的运动规律理论模型为:
连续性方程:
式( 1) 中: ρ 代表流体密度,ui代表着三个方向的速度分量。
动量方程:
式( 2) 中: p代表静压; μ 代表黏性系数; gi代表重力加速度分量。
能量方程:
式( 3) 中: T代表流体温度; k代表热传导系数; cp代表定压比热。
状态方程:
采用改进型的k-ε 模型来模拟分层流运动过程中流动与传热的强耦合问题。改进型的k-ε 模型可由瞬态的N-S方程组导出,它与标准的k-ε 模型具有相似的方程形式。将水视为不可压缩流体,改进型的k-ε 模型中湍流动能k的方程为:
耗散率 ε 的方程为:
湍动能产生项:
式中:,S为平均应变率; ρ 代表流体密度,ui代表着三个方向的速度分量;为湍流黏性系数; Gk和Gb分别代表平均速度梯度引起的湍动能和浮力作用引起的湍动能的影响; Prt代表普朗特常数,在改进型的k-ε 模型中,代表热膨胀系数;C3ε代表浮力对耗散率的影响,,ν 代表与重力方向平行的速度分量,μ 代表与重力方向垂直的速度分量。C1ε、C2、σk、σε均为常量,在CFX中,默认有C1ε= 1. 44,C2= 1. 9,σk= 1. 0,σε= 1. 2,其中,σk、σε分别代表湍动能普朗特常数和耗散率普朗特常数。
2数值模拟方法
使用CFX软件求解N-S( Navier-Stock) 方程,采用了改进型的有限体积法,利用全隐式多网格耦合求解技术进行数值求解,利用Ansys CFX求解时采用SIM2PLEC算法和一阶迎风格式,考虑重力影响, 所有方程均要求同时达到10- 4精度以上。
网格划分:
为保证计算精度与计算速度,针对模型特点,使用ICEM CFD为计算域生成结构化六面体网格,当网格数量达到309万个时,网格数量满足网格无关性要求。
入口边界条件:
u、w、v分别为计算域入口边界处x、y、- z方向上的速度分量,T为计算域入口的边界温度,来流流速为v0,来流温度函数为Tw,有:
壁面边界条件:
壁面设定为滑移壁面,壁面y方向温度分布定为来流温度函数Tw;
出口边界条件:
设定为流量边界条件,通过计算求得流量参数q使得qin= qout。
时间步长选为物理时间步长:。
3温度分层流的实验模拟
温度分层流的模拟实验是在分层流水槽中进行,如图2所示。分层流实验水槽的长宽高分别为8 m × 1 m × 1 m,水深为0. 85 m。水槽以三角铁为框架,壁面采用可视效果较好的钢化玻璃,并用玻璃胶进行密封粘接。
实际海洋中的温度分层是由于太阳辐射造成的,为更加逼真模拟实际环境,实验采用40盏功率为275 W的辐射灯分两排架设于离水面0. 3 m的高度处对水照射加热 ,如图3所示。 通过辐射传热的方法,形成温度分层。实验中,初始水温为27. 3 ℃ ,辐射灯照 射2 h后停止。利用K型热电偶测量分层流体的温度,热电偶按照预定的深度布置在垂直插入水 中的测量 杆上,共24个测点。使用安捷伦数据采集器记录热电偶的测量数据,测点的深度位置以及通过测点测得温度详见表1。
4模拟结果分析与对比
将实验数据进行线性拟合,得到来流温度函数Tw,代入模型进行计算,得到如下仿真结果,如图4所示。
在仿真计算的结果中,在任意的一个XY截面上,取x = - 900、0、700 mm时沿y轴方向上的温度分布进行对比,从图5可以看出,温度分布曲线完全重合,由此可见,在X轴方向上,温度分布均匀。
沿Z轴方向上取任意的三个XY截面,在截面上取沿Y轴方向上的温度分布进行对比,从图6可以看出,温度分布曲线完全重合,由此可见,在Z轴方向上,温度分布完全相同。采用原先的方法模拟的流动入口、中部和出口的温度分布如图7所示,随流体向下游流动,温度分层向深度扩散,上游与下游的温度分层不能保持稳定。采用本文提出的方法后,之前方法中温度分层不能维持的情况被成功消除,有效维持了温度分层的稳定。
任取一截面的温度分布与分层流水槽实验结果进行对比,从图8可见,温度分布几乎完全相同,在允许误差的范围内,该种方法与实际情况相吻合。
5结论
本文建立了描述温度分层流运动的质量、动量和能量守恒方程,采用改进型的k-ε 湍流模型模拟了温度分层流运动的流动与传热强耦合过程,得到的主要结论有:
运用有限体积法,提出了以滑移壁面和给定壁面温度函数为核心,结合来流温度函数的边界条件设置,对现有温度分层模拟方法进行了改进,成功实现了计算域内稳定温度分层的形成与保持。通过与前人方法相对比,证明了本文提出的温度分层流数值模拟方法的有效性。
采用辐射加热灯加热水槽流体表面,成功实现了流体的温度分层。数值模拟与分层流水槽试验结果对比,证明了本文数值模拟方法的可行性与准确性。
数值模拟实验 篇8
飞机结冰是指水分在机体表面凝冻成 冰层的现 象。结冰不仅会增加飞机的飞行重量,同时也改变了飞机的气动外 形,影响气流的流动情况,造成阻力增大,升力减小,操纵性、稳定性恶化[1,2]。目前应用最为广泛的热气防除冰系统在发动机引气不足的条件下会受到一定制约且能耗较大;而电热防除冰系统的能耗达到26kW;电脉冲除冰系统的能耗为3kW[3]。近十年来国内外提出了一种应用于飞机的 新型除冰 技术———压电除冰,但目前还处于实验研究阶段[4],本文将探 讨这一除 冰技术,通过数值模拟与实验论证其可行性与节能性。
1压电除冰原理
飞机机械除冰方法是通过机械力破坏冰和黏附 层之间的黏结正应力或者剪应力,使冰破裂以致从黏附层脱落。通过人们研究发现,冰从物体表面剥离所需的拉压正应力是剪切应力的几十倍[4],根据这一特性,利用压电材 料的逆压 电效应制 作压电驱动器来破坏 机翼蒙皮 和冰层的 弱剪应力,除冰效果 将更佳。
2压电除冰的有限元分析
为了验证理论研究的可行性并指导后续实验,采用有限元分析软件ANSYS进行数值模拟。通过模态分析与谐响应分析求解出研究模型不同振型下冰和蒙皮之间的剪切应力,利用计算结果获得压电驱动器的大概安装位置与实验激振频率。可分为无压电片安装时的数值模拟与压电片安装时的数值模拟。
2.1无压电片的试验蒙皮数值模拟
本文采用低速翼型NACA0030进行仿真模拟。由于实际情况下冰层主要集中在机翼前缘,为了研究问题的方便,对模型进行简化,只取机翼前缘部分作为研究对象。整个试验蒙皮的尺寸分别为:展长0.3m,翼型厚度0.3m,蒙皮厚2.0mm,其物性参数为:密度ρAl=2.78×103kg/m3,弹性模量EAl=7.05×1010Pa,泊松比λAl=0.33。采用Solid45划分该试验蒙皮,并约束断面处节点的全部自由度。
建模完成后对该试验蒙皮进行模态分析,可以得到机翼的多个固有频率与振型。当驱动器的激振频率达到构件的固 有频率时,振动的振幅将达到最大,但不同固有频率下激振的 振幅也不同,找到合适的共振频率和振型就成了首要任务,因为这是作为压电除冰系统设计的基础。本计算取前10阶模态进行分析,并观察其对应的振型与固有频率,图1为第六阶 振型示例。根据图1可以得到整个构件的最大形变位于:x=0.09,y=0.11,z=±0.3和x=0.18,y=0.148,z=±0.3处,故将驱动用的压电片粘贴在该位置再做进一步的研究。
2.2带压电片时的试验蒙皮数值模拟
根据2.1所述,选择大小为0.05m×0.05m×0.003m、材料为PZT-8的压电片两块粘在试验蒙皮形变最大处,并建立相应的压电驱动分析模型。蒙皮采用Solid45,压电片采用耦合单元Solid5,同样约束蒙皮断面后断点节点所有自由度。
模型建完后再次进行模态分析验证粘贴好压 电片后其 振型变化的差异,从求解结果来看,尽管频率有所变化,但形变位置变化不大。由此,可考虑进一步对压电片施加 电压,并设定激振频率范围进行谐响应分析。确定扫频范围0~6000 Hz,压电驱动电压为320V,求解可得到不同频率下的应力分布情况。当激振频率为1600Hz时,所读取的节点XZ方向的最大应力为3.42 MPa,YZ方向的最大剪应力为2.34 MPa,XY方向的正应力为4.74 MPa,其等效应力为13.5 MPa。在分析中通常选取2 MPa作为剪切黏附强度[3],通过该计算结果可认为利用逆压电效应进行除冰是可行的。因此,将搭建实验台进一步论证。
3压电除冰实验
3.1实验设备
压电除冰实验台除 如上设计 的蒙皮实 验件、压电驱 动器外,还需要配备信号发生器、提高输出功率的功率放大器、压电传感器、采集图像的示波器与制冰冰箱。
3.2实验过程及结果
按照ANSYS仿真求解的结果,把压电片粘贴在机翼的内表面,并根据3.1所述的设备完成实验台的搭建。采用人工结冰方法完成蒙皮结冰,再进行除冰实验。使用信号发生器输出正弦波交流电,首先利用较低的电压改变信号发生器的输出频率,在数值模拟得到的多个共振频率附近,通过观察示波器 波形变化来寻找实 验机翼的 共振频率。 发现当激 振频率达 到1530Hz、4320Hz、4830Hz、5920Hz时,构件的振幅将相对达到最大。选定1530Hz作为除冰激振频率,并调整幅值电压至650V,输出峰值电流为0.1A,观察到试验蒙皮上的冰块可成功去除。与计算所述1600Hz能去除冰层的结论相比,其驱动频率误差小于10%;同时所耗功率仅为32.5 W,远远小于电热除冰系统与电脉冲除冰系统所耗的电能。
4结语与展望
数值模拟实验 篇9
国内大多数换热器流速偏低, 其主要原因是设计保守和用户在传热能力不能满足要求时加大换热面积或者增加换热台数, 造成流速更低, 容易结垢, 运行效率更低, 甚至进一步恶化[1]。污垢使换热器的运行效率下降50%, 这种低流速的换热器数量较多, 实际生产中对污垢自动清洗技术的企盼也较迫切[2]。换热管内插物是一种经试验验证能有效提高传热效率的方法, 其中行之有效的2种内插物是自转式强化传热转子与螺旋扭带, 如图1～5所示[3]。
转子与扭带在管内流体的推力作用下旋转, 增加流体湍流程度, 起到了强化传热的作用。为了更好地对比研究两种换热元件的换热性能优劣, 参照实际工况, 首先从数值模拟的角度, 采用CFD流体分析软件Fluent 6.0对不同导程的转子与螺旋扭带进行传热性能的对比研究, 并引入光管进行参照。
表1所示为转子与扭带参数表。
1 模拟部分
由于换热管的L/D较大, 因此为了有利于分析, 对比实际工况, 适当缩短换热管长度, 取2m, 流体流速适当降低, 设管程热流体进口速度1.9m/s, 壳程冷流体入口速度1.5m/s, 方式为逆流传热。其中热流体进口温度60℃, 冷流体进口温度为20℃。为了突出有无转子和扭带的效果, 尽量减小其他因素的干扰, 将内管管壁设为一理想薄片, 其热阻始终为0。
1.1 控制方程
模拟包括流体/固体间换热, 换热类型包括热传导、对流换热。求解的能量方程如下:
undefined
式中:keff—有效导热率, keff=kt+k;undefined—组分j的扩散通量;Sh—化学反应热和其他体积热源。
在式 (1) 中, 右边分别表示由于热传导、组分扩散、粘性耗散而引起的热量转移。
其中, undefined, 对于不可压流体, undefined。
式中:mj—组分j的质量分数, 组分j的焓定义为undefined, 其中Tref=298.15K[4]。
由于介质为水, 组分扩散, 粘性耗散等引起的热量可以忽略, 则固体计算域的能量方程可简化为[5]:
undefined
1.2 边界条件
进口边界:采用速度进口形式, 管程入口速度为1.9/s, 壳程入口速度为1.5m/s;
出口条件:壳程与另一端管程出口均采用outflow形式;
流畅区域:管内壁与转子或扭带外表面的之间的空间, 壳程与管程之间的空间;
壁面条件:管壁壁面采用无滑移固定粗糙壁面, 管程、壳程内管壁耦合换热, 管程、壳程外管壁分别为耦合换热和绝热。
1.3 网格划分
为了避免划分过程中出现最小负体积, 采用四面体网格划分, 总网格个数为203万。
1.4 模拟结果分析
为了定性研究转子与扭带在换热过程中的性能优劣, 在其他条件不变的情况下, 该研究中将导程为100mm、200mm、400mm的组合转子以及100mm+200mm交替排列、200mm+400mm交替排列的组合转子与导程200mm的扭带进行实验对比, 引入光管作为参照, 图6所示为转子与扭带在换热管内工作情况。
为了更直观了解其在管程中的状态, 将截图75%透视, 图7所示为数值模拟光管、各导程转子、扭带在强化传热过程中横截面温度分布云图。热流体自上而下进入内管 (管程) , 冷流体自下而上进入外管 (壳程) , 冷热流体在逆流过程中实现热交换, 在转子或扭带的作用下达到不同程度的强化传热效果。
由图7可以得出数值模拟结论:
(1) 由于光管内流体扰流较小, 隔热边界层较厚, 换热不均匀, 换热系数较低。
(2) S=100转子:导程比较小, 旋转速度较快, 破坏管内壁流体边界层, 起到很好的强化传热的效果, 但是由于阻力较大, 泵功率一定的情况下, 管程流体流量小, 在管外流体换热量一定的条件下, 管内流体进出口温差与流体流量成反比, 因此温差较大。
(3) S=100+200转子:换热均匀, 换热系数较高。转子阻力减小, 管内流体流量比导程100转子大, 而管外冷流体得到充分换热, 比较适合工业应用。
(4) S=200转子:热流体很快从60℃降到45℃左右, 冷流体温度变化均匀, 沿管长方向温度均匀增加, 可见S=200转子换热效果也很明显。
(5) S=200+400转子:管内流体换热均匀, 但温度从60℃降到45℃所用时间比较长, 分析原因:阻力小, 因此流量比S=200大;对边界层的扰动不如S=200明显。壳程流体与S=200比温度无明显变化。
(6) S=400转子:冷热流体两端温差比前几种转子的小, 换热系数不高, 但强于光管。
(7) S=200扭带:换热均匀, 管外冷流体两端温差与S=400转子几乎相同, 但由于扭带运转具有整体性, 对隔热边界层破坏不明显, 并且流体经过扭带的径向边缘时, 与垂直管束的传热方向的夹角较大, 根据场协同原理, 换热效果明显不如导程S=100、S=100+200、S=200、S=200+400转子, 其换热效果与S=400基本相同。
2 实验部分
为了更好地检验模拟结果, 将以上几种换热元件做实验并与模拟结果进行对比。实验装置由动力系统 (电机和多级泵) 、换热系统、加热系统、数据采集系统、测试系统以及控制系统等组成, 实验系统采用计算机自动控制和自动数据采集。应用基于模糊自适应方法的自动控制系统, 保证实验段进口流量和温度的稳定可靠, 可以实现实验段进口流量±2%和进口温度±0.5℃的控制精度。壳程套管尺寸为ϕ47×3.5, 管程换热管直径ϕ25×1.5。实验装置密封性能良好, 管外介质为冷水, 管内介质是热水。实验段、所有管线及储液箱等采用高发泡聚乙烯材料保温, 以减少热损失、提高实验精度。
2.1 数据处理
该实验用总传热热阻减去壳程对流换热热阻以及管壁热阻来替代, 直接用威尔逊图解法计算, 这样能避免相对较大的误差, 即:
undefined
其中, undefined
式中:K—总传热系数, undefined;A—换热面积;Rw—管壁热传导热阻, undefined;下标i—管程;下标o—壳程;ΔTm—对数平均温差;Qi—总传热率;C—威尔逊图解法的斜率。
根据Dirker、Meyer等人的研究, 将m定义为0.8。根据Dittus-Boelter公式, 将n定义为0.4[5]。
2.2 实验结果
图8所示为光管、各导程转子、螺旋扭带随管程流体流速的增加, 换热系数的变化曲线, 引入流速1.9m/s时数值模拟的结果作为比较。
注:数值模拟传热系数的计算结果没有考虑内管管材热阻, 因此由实验所得传热系数稍大。
2.3 实验结果分析
由图8实验结果可知:
(1) 无论转子还是螺旋扭带, 与光管相比都达到了强化传热的目的;
(2) 转子强化传热的效果普遍好于螺旋扭带的强化传热效果;
(3) 流速为1.9m/s时的模拟结果与实验结果基本吻合 (若考虑内管管材热阻情况下) ;
(4) S=100 、S=100+200转子对层流边界扰动最明显, 传热系数最高。
3 结论
通过数值模拟和实验结果, 总结出自转式转子、螺旋扭带与光管的换热性优劣, 结果如下:
K100转子≥K100+200转子>K200转子>K200+400转子>K400转子>K扭带≈K400转子>K光管
由此可知, 管内插强化传热元件中, 螺旋转子破坏传热边界层的能力突出, 强化传热效果好, 其中, 叶片螺旋角的大小, 即导程是影响转子强化传热能力的重要因素, 导程越大, 对流体扰动越明显, 但是随之阻力也增大, 耗能增多。
综上分析, 从提高换热系数、降低阻力的角度考虑, 转子导程为100+200的经济效益最优, 在电厂冷凝器、化工厂换热器等换热设备的强化传热、节能减排方面具有较高的应用价值。上述结果对螺旋转子在工业应用中具有指导意义。
参考文献
[1]俞天兰, 蒋少青, 刘桂英, 等.我国冷冻水冷器设备污垢影响的能量损失研究[J].低温工程, 2002, (3) :34-36.
[2]蒋少青, 俞天兰, 彭德其, 等.外螺旋塑料管自动清洗及其传热强化技术的研究[J].冷饮与速冻食品工业, 2005, 11 (1) :15-18.
[3]洪蒙纳, 邓先和.管壳式换热器管程强化传热研究进展[J].广东化工, 2005, (3) :41-42.
[4]王瑞金, 张凯, 王刚.Fluent技术基础应用实例[M].北京, 清华大学出版社, 2007.
[5]李月, 丁玉梅, 杨为民.内置洁能芯换热管强化传热数值模拟研究[J].石油机械, 2009, 37 (10) :14-17.
数值模拟实验 篇10
磁性液体是由单畴铁磁性颗粒构成的高稳定胶状悬浮液。它是由近似球形、直径约10nm的铁磁性或亚铁磁性微粒(微粒表面包覆着直径约为2nm且具有长链分子结构的单分子层表面活性剂或分散剂)高度弥散于选定的基载液中,所构成的具有高稳定性的固液两相胶体溶液。无论是在重力场,还是电磁场作用下,磁性液体都能长期稳定地存在,是一种极其稳定的实用新型纳米功能材料。
在传感研究与应用方面,磁性液体广泛应用于航空、航天、宇航站等尖端国防军事领域,以解决各种情况复杂、条件恶劣环境下的检测问题[1,2,3,4,5,6,7,8,9,10,11]。由于国外对磁性液体系列传感器的研究绝大部分集中在军事领域,以解决苛刻条件下的测试问题[5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23],所以,公开报道材料很少。在可查阅的公开文献中,美国、英国、印度和罗马尼亚等国进行磁性液体及其传感研究主要集中在动力传感及其在航空装置中的应用[7]、惯性传感测量装置[10],以及磁性液体的电磁学特性对电量和非电量传感器性能的影响[15]。我国对磁性液体水平传感器的研究尚处于探索起步阶段,对于如何选择磁性液体传感器的设计参数、如何量化影响其灵敏度等输出特性的因素等方面的问题,仍缺乏系统的理论分析与研究。
本文在研究磁性液体水平传感器检测原理的基础上,建立数学模型,利用有限元分析软件对磁性液体水平传感器的工作状态进行模拟,研究影响其灵敏度等输出特性的因素和规律,并与实验结果进行对比和分析。
1 数学模型
图1所示为磁性液体水平传感器的基本结构形式。它由内径为r并充有一定体积磁性液体的非磁性圆柱形容器,以及在其上均匀缠绕的三组线圈构成。三组线圈的匝数、长度、绕制方式相同,两个感应线圈S1与S2的中心距离为G,对称分布在激励线圈S0两侧。按图示方法连接外电路:S0外接电压源以提供激励电压Uin;S1与S2反相串接后,连接输出电路以提供输出感应电势差Uo。
假设,当容器处于水平位置时,L0、L1与L2分别为线圈S0、S1与S2的自感系数,M0为S1、S2与S0间的互感系数;当容器发生小角度倾斜时,S1与S2的自感系数分别变为L′1与L′2,S1、S2与S0间的互感系数变为M。
又有,当容器水平时,感应线圈S1与S2的总自感L为
当发生小角度倾斜时,感应线圈S1与S2的总自感变为L′,即
且有
设激励线圈S0通以电流I,则感应电动势为
又因为L1=L2,因此,当管状容器水平,线圈所包围的磁介质未发生变化时,输出感应电动势为零。
当容器相对于水平面发生微小倾角时,感应电动势为
则输出电动势为
因感应线圈S1和S2均为密绕螺线管,则
式中,μ为管状容器中磁性液体与空气的总磁化率;N为感应线圈的匝数;S为线圈的横截面积;d为线圈缠绕的线性长度。
当容器相对水平位置发生微小倾斜时,设ΔL为自感变化值,Δμ为磁导率变化值,则有
又令ΔM为互感变化值,则有
化简得
当容器相对水平位置发生小角度倾斜时,忽略Δμ的高阶无穷小;设I0sinω t为通过激励线圈S0的电流,当激励线圈与感应线圈的匝数、几何形状、尺寸都相同时,有L0=L,则输出电动势,即式(1)可以化简为
设k=Vm/V,其中,Vm为磁性液体体积,V为管状容器容积。假设容器中只存在着空气与磁性液体的混合物,那么器中磁性液体与空气的总磁导率可表示为
μ=kμm+(1-k)μa
式中,μm为磁性液体的磁导率;μa为空气磁导率。
设ΔVm为磁性液体体积变化值,则当ΔVm≈2rhd≈rGdθ时,两感应线圈内部磁液体积发生变化,则有
那么,在两个感应线圈所包围的范围内,参与作用的磁性液体的磁导率也发生变化,有
又由于,μm=μ0(1+χm)且μa≈μ0,其中,χm为磁性液体的磁化率,即磁化系数,μ0为真空磁导率,所以有
则式(1)简化为
式中,n为单位长度上激励线圈的匝数。
值得指出的是:由前面的推导过程可知,当注入的磁性液体体积是容器容积的一半时,输出电势差Uo最大,即可以测得的一维水平倾角最大。因此,磁性液体水平传感器水平倾角的最大测量范围与该传感器的结构参数有关。
输出电势差通过放大、整流、滤波后,式(3)可以表示为
Uo=PGrχmnNUinfinθ (4)
式中,P为常数;fin为激励信号的输入频率。
磁性液体水平传感器的灵敏度定义为
ζ=Uo/θ=PGrχmnNUinfin (5)
由式(4)易知,当水平倾角θ很小时,输出电势差Uo与水平倾角θ成正比。因此,被测姿态的一维水平倾角这个非电量,可以由电信号测量。
由式(5)易知,磁性液体水平传感器的灵敏度与在其中作为敏感元件的磁性液体的磁化参数,盛有磁性液体的非磁性容器的结构参数,容器上所缠绕的各组线圈的结构、尺寸和位置参数,以及激励信号的输入电压、输入频率参数有关。而且,当所测倾角很小、结构参数确定之后,该种传感器的灵敏度ζ分别与磁性液体的磁化系数χm、激励信号的输入电压Uin以及激励信号的输入频率fin成正比,而与被测倾角大小无关。
2 数值模拟
2.1 磁导率的处理
在传统的电磁学有限元分析和数值模拟中,通常只出现单纯的非磁性流体或稀土永磁体,它们的磁导率是一个定值。就磁性液体水平传感器而言,当被测姿态相对于水平面发生微小倾斜时,盛装在容器内的磁性液体将随着装置的倾斜而由高向低流动,同时,磁性液体具有磁性,而且其磁性的大小和磁性液体单位体积的浓度有关[12]。因此,当磁性液体发生流动时,即磁性液体水平传感器在不同倾斜位置时,感应线圈所包围的不同区域内的磁液体积会随着被测姿态的倾斜而改变,这可以看作是在相应区域参与作用的磁性物质的磁导率随着倾角的变化而改变。
因此,当被测姿态相对于水平面发生小角度倾斜时,由于角度变化不大、磁场变化也不大,可以认为B-H曲线是线性的。另外,考虑到磁性液体的连续性,磁性液体的位置是随着倾角的确定而确定的,可以理想地认为在确定的瞬间,磁性液体的磁导率是定值。如式(2)所示,这个瞬时磁导率显然是倾角θ的函数。
2.2 电磁耦合分析
虽然,就磁性液体水平传感器而言,激发磁性液体中磁性颗粒磁化的外磁场是激励线圈外接电源产生的瞬态时变电磁场,它在磁性液体内部会产生附加磁感应强度,从而造成磁性液体磁化不均匀。但是,从宏观上看,磁性液体的磁化是均匀的。在这种情况下,电场与磁场是一个整体,电磁场耦合且不可分割,所以其电场分量与磁场分量存在确定的依赖关系。
把磁性液体区域看作是由有限多个球形铁磁性颗粒与非磁性的基载液组成的区域,由于粒子运动的时间远远大于电磁作用的时间,可以认为磁性液体体系遵循静态磁场方程:
式中,H为磁场强度矢量;j为电流密度矢量;B为磁感应强度矢量;∇为哈密尔顿算子。
由于磁性液体中的磁性颗粒按照外磁场的方向磁化,所以磁化强度方向与磁场强度方向一致。对于磁性液体,引入表达式:
式中,M为磁化强度矢量;μr为相对磁导率。
为了描述B与H之间的关系,可以定义方程
B=μ0μrH=μH
式中,μ为总磁导率。
电磁场麦克斯韦方程组为
式中,D为电位移矢量;ρ为电荷密度;E为电场强度矢量。
引进标势φ与磁矢势A,由于在低频情况下,磁化是线性的,可以假设全部空间充满各向同性、线性介质,则有
D=εE
式中,ε为介电常量。
则原麦克斯韦方程组可以化简为
∇
∇
式(6)与式(7)即为电磁场耦合微分方程。
2.3 模型建立与求解
根据图1所示结构及磁性液体水平传感器的工作原理,建立有限元分析模型。设定各设计参数值,其中三组线圈分别缠绕三层,三层匝数共为3780,导线直径为0.05mm,磁液体积为容积的一半。容器内外直径分别为25mm、30mm,内壁总长210mm,可测倾角范围为0～0.11rad。通过分析可知,磁性液体水平传感器核心工作部分的结构可以由三维轴对称问题简化为四分之一结构的二维平面问题。然后,在ANSYS前处理器中建立基于磁性液体的水平传感器的结构模型。
模型建立后,再确定此磁性液体水平传感器模型结构的物理环境。该模型的物理环境主要包括单元的选择及其属性的确定,磁性液体、非磁性管状容器、线圈导线以及空气等非磁性材料性能数据的输入等几个方面。首先将选项框设为Electrical-Magnetic Nodal;采用PLANE53平面单元、INFIN110远场单元以及CIRCUI124电路单元,建立2-D瞬态电磁场分析的物理环境。材料性能数据可以从已经定义好的材料性能数据库中直接输入,也可以通过材料模型对话框输入。所选用的线圈导线材料为铜丝,电阻率为1.7×10-8Ω·m,模型中没有永磁材料。由于磁性液体的磁性能是磁性液体水平传感器感应信号变化的关键参数,所以不能像以往研究有永磁体材料存在时,认为磁性液体相对于空气的相对磁导率是1,而应将磁性液体相对磁导率设为略大于1的数值,例如MURX=1.3。另外,容器、线圈以及空气的相对磁导率均设为1。
物理环境设定后,针对模型的不同部分赋予相应的材料属性,然后进行合适的网格划分。接着对空气边界施加磁力线平行的边界条件,即定义矢量磁位Az=0,而模型内部的边界条件由ANSYS自动施加。最后,选择合适的求解器进行求解。经分析得到的磁力线分布如图2所示。
3 实验过程
适用于磁性液体水平传感器的磁性液体既要具有良好的流动性,又要具有较强的饱和磁化强度。经过前期实验的比较与分析,由于煤油基磁性液体的黏度较低,可以使磁液具有良好的流动性的同时,保持较高的磁化强度,因此,选择煤油作为基载液,Fe3O4作为固相磁性纳米颗粒,油酸作为包裹于纳米颗粒外的表面活性剂,制备煤油基磁性液体。制备出的磁液主要性能参数为:密度1.15×103kg/m3,黏度2.747Pa·s,磁化强度370.09×103A/m。
按照数值模拟的参数,设计实验模型,进行实验验证,如图3所示。
1.磁性液体水平传感器实验模型 2.信号发生器 3.电路系统 4.示波器
4 结果与分析
4.1 数值模拟与实验验证结果
通过有限元分析,图4、图5、图6是由数值模拟结果绘出的曲线图。
图7、图8、图9是由实验验证的结果绘制的曲线图。
另外,实验所用磁性液体经过稀释,可以改变其磁化系数。设稀释后与稀释前的磁性液体磁化系数之比为χ′m/χm。当输入电压确定时,灵敏度与磁性液体相对磁化系数χ′m/χm的关系曲线如图10所示,图示表明传感器模型的灵敏度与磁性液体的磁化系数呈线性关系,与理论分析式(5)的结果相一致。因此,在保证磁性液体具有良好的流动性的前提下,增加磁性液体的磁化系数,可以相应地增强传感器的灵敏度。
4.2 模拟与实验结果分析与对比
在相同的输入电压、输入频率参数条件下,将数值模拟结果与实验结果进行对比,结果如图11、图12及图13所示,由图可知,磁性液体水平传感器模型工作状态的数值模拟结果与其实验结果基本一致,变化规律相吻合。不论是数值模拟结果,还是实验结果,都符合式(4)、式(5)的线性关系。
经分析,值得指出的是:利用最小二乘法,得到实验结果各拟合直线的相关系数r的变化范围为0.97467<r<0.99969,由此分析此传感器线性度指标,取显著性水平α=0.05,实验组数为6,查得相关系数临界值rmin=0.811。由此可见,此传感器的输出电压与输入量即被测水平倾角之间具有较好的线性度。数值模拟结果与实验结果存在差异,经分析认为这主要是由实验模型制造精度不高、线圈缠绕部分制作不均匀、不对称等加工制作因素引起的。这些因素也是引起模型实验结果出现零位误差的原因,除提高加工精度外,还可以采用电路补偿的方式来消除这些误差。如何更好的消除这些差异,也是今后研究工作的一项重要内容。
5 结论
(1)磁性液体水平传感器的最大量程,与其结构参数有关。当结构参数确定之后,磁性液体灵敏度正比于磁性液体的磁化系数、激励信号的输入电压有效值、激励信号的输入频率,而与被测倾角大小无关。
(2)将磁性液体在某一瞬时的磁导率看作定值,可以简化磁性液体水平传感器模型工作状态的数值模拟复杂性。经实验验证,数值模拟结果与实验结果基本一致。
(3)当其他条件不变,随着输入电压的增加,磁性液体水平传感器的灵敏度线性增加,不受其他实验参数的影响;而输入频率信号的取值范围在500～2000Hz之间,磁性液体水平传感器灵敏度与输入频率信号的线性度更好。
数值模拟实验 篇11
关键词:土质边坡;渐进性破坏;数值模拟;物理相似模拟试验
边坡的破坏是一个渐进的过程,这个概念早在20世纪60年代,Skempton在研究超固结粘土边坡的稳定性时就提出了。他认为,土体的强度并非在整个滑裂面上同时发挥作用,而是当土中某一点的剪应力增加到超过土的强度时,该点发生剪切破坏。由于这种剪切破坏逐渐传递,使得渐进破坏面逐渐扩大,最后当滑动推力上升到超过滑裂面的抗剪强度时便发生坡体的整体滑动。所以,滑坡是一个土体随时间而缓慢变形并最终突然崩塌的过程。
对于一定的坡型和坡体材料,其坡脚的塑性区范围是一定的。而一旦出现局部塑性破坏区,该区域的承载能力会明显下降,坡体应力场会出现重新调整,并导致新的塑性破坏区的产生。新的塑性破坏区在形成过程中,释放残余能量,该残余能量将向塑性破坏区周围的弹性区扩展,使之受到其影响,而受到影响最大的是其临近区域。这些区域的某些范围可能由于叠加到相邻塑性破坏区所释放的能量,而超过其弹性强度,从而进入塑性破坏区。这样在不断的能量释放、转移、调整过程中,塑性破坏区范围不断扩大,破坏面不断延伸,释放的能量不断减小,当到达一定程度,弹性区在叠加了残余能量以后,所受到的应力仍然小于材料的弹性强度,则塑性破坏区不再扩展,破坏面不再扩大。但是该叠加应力若大于材料的弹性强度,则塑性区继续扩展,直至破坏面形成。本项研究基于渐进破坏理论,采用有限元数值计算分析方法,对土质边坡的渐进破坏过程进行模拟,并通过室内的物理相似模拟试验对其进行验证。
一、土质边坡渐进破坏分析程序设计
1.有限元建模
本程序采用有限元建模FORTRAN90编写程序。由一个主程序和30个子程序所构成。程序编写时,把主要部分划分成几个子块,每一子块完成一项工作,并将可能用得比较多的程序段,用一个个小的子程序编写出来,便于多次调用。主程序主要用于控制运算路径和进程,各个子程序分块编写,分级调用。这种模块化的结构体系便于程序的编写、修改、调试,是化繁为简的行之有效的方法。
程序编写根据有限元建模的一般步骤进行,将网格自动生成、单元刚度矩阵、整体刚度矩阵、荷载列向量、解方程出位移等编写成子程序,在主程序MAIN中顺次调用。经过程序的一次运算后,可以得到边坡模型中各单元形心应力及处理后的绕节点平均应力。
2.FORTRAN程序实现渐进破坏过程
上述一次程序运算完毕,求出边坡的弹性应力场,计算出节点位移、单元及节点应力。然后根据德鲁克-普拉格强度准则求出系统的超过强度条件的最不安全单元,标记下来,对该单元进行处理,将其强度降低,并将其过量应力按照弹塑性力学的原理分配到该单元节点上。对破坏单元进行处理通常是将弹性模量E降低,一般降低为原弹性模量的 左右。修改好破坏单元的参数后,一个新的模型诞生。该模型中,与前次破坏单元相关的节点获得了破坏单元卸载的过量载荷,使得这些单元较其他单元更有可能称为下一破坏单元。也就是说,与破坏单元相邻的单元是下一级破坏单元产生的区域,这与很多理论研究和实验结果都是一致的。
处理后的单元与其他单元一起进入程序第二次运算,得到新的破坏单元,对新的破坏单元给予同样的强度折减处理,再运算,如此循环,直至破坏单元的数量不再增加,说明破坏面延伸至此终止,程序运行结束。
这样,在载荷不变的情况下,坡体从第一批屈服破坏单元开始出现,直至最终滑动面生成,清晰的体现了渐进的思想。屈服破坏区域逐渐扩大,破坏单元逐渐增多,破坏滑动面逐渐“生成”,渐进的过程一目了然。
3.程序运算结果
经过本程序计算后,破坏单元呈“渐进”方式出现,即破坏单元数目越来越多,破坏区域越来越大。将破坏单元渐进出现的区域及次序用图1展示。图1中,标有数字的单元就是破坏单元,数字1、2、3表示的是单元破坏的次序。图中所示破坏单元构成的是一个破坏区域,而不是理论上所说的破坏面,这是因为程序计算所取精度的关系。
二、物理相似模拟试验
1.试验设计
根据边坡实际尺寸和模型架装置尺寸,选择模型边坡顶部宽度为0.4米,坡高0.65米,坡体底面宽度取1.05米,模型厚度为0.25米。
填筑模型所用材料全部采用现场采集的砂土,进行破碎、筛分处理后,分层填装。根据路基填筑施工中采用的分层高度,并考虑模型填筑的方便,模型填筑的分层高度取5cm。整个边坡模型高度为0.65m,共分13层进行填装。在每次分层填筑后,于适当位置放置自制观测元件,直至整个模型填筑完毕。自制测点交错布置,呈梅花桩形式排列,共布置8层24个测点。
模型填筑完成时,用相机记录各个自制测点的位置。模型填筑完毕一天后,进行加载。荷载施加完毕后待模型稳定时,每间隔一段时间再用直尺测一次位移变化,并用相机记录下此刻的测点位置,从而得到各测点的位移变化量。
2.试验结果分析
在模型填筑完毕,荷载加上去不多时,可以很清楚的发现,测点位置发生了改变。随着荷载施加,自制测点的位移开始发生比较大的变化,并在坡脚附近出现第一道裂纹。一旦裂纹出现,坡体的强度必然发生比较大的变化,可以看到,坡体里面的测点发生了很明显的位移变化。随着荷载施加,裂纹开始开展,慢慢“长成”为一条裂缝,见图2与图3。
可以看到,在裂缝的上下测点位移的变化情况有了很大的不同,裂缝以上的部分,测点位移发生了相当大的变化,可以看到,从第一幅到第三幅,位移变化有5cm之多,且对照三幅图斜坡面上竖向裂缝的位置,说明坡体沿着裂缝开始下滑了。而裂缝以下部分的测点位移基本没有什么变化,这也从一定程度上说明该试验坡体的破坏不是整体的坍塌过程,而是裂缝局部发展扩大的渐进过程。
需要说明的是,由于土质拌合均匀程度、夯实均匀程度以及实验条件等一些问题,使得裂缝不像理论值那样规整,但裂缝的形状已能在一定程度上说明问题。通过试验记录下来的数据与用相机拍摄的图片,证明了试验结果与本文数值模拟的结果一致,坡体在自重和外载作用下发生的渐进破坏是从坡脚开始,并逐渐发展形成最终的滑动面。
三、結语
边坡破坏不是瞬时发生的整体破坏,而是从坡体局部出现裂缝,并一步步“生长”,形成滑动面的渐进破坏过程。
经过自己编写程序的运算,可以追踪到边坡破坏单元渐进出现的过程,从而证明边坡破坏渐进理论的正确性。
数值模拟实验 篇12
近20年来,由植物油制备生物柴油作为石油燃料的替代物,已引起了世界各国的广泛关注。且主要以实验手段为主,研究各种参数(温度、压力、催化剂等)对转化率的影响,故关于其工艺、原料和产品后处理方面的报道颇多[1,2,3,4]。由于超临界环境下反应的复杂性,关于数值模拟结果方面的报道很少。反应温度影响反应速率和时间,准确模拟反应装置中的传热过程极其重要。功能强大的CFD商业软件用在传热方面的研究效果显著[5,6],故采用CFD软件对一套高温高压间歇式反应装置进行了数值模拟,希望能进一步了解反应过程并为此类反应装置的设计、运行和改造提供参考信息。
1实验装置介绍
反应装置平面图如图1所示,锡池尺寸为50×138×50mm。该装置中根据反应的需要设定锡池的温度,反应釜和反应液处于同一温度(环境温度),以锡池加热反应釜和反应液。经过一段时间传热后,三者温度达到平衡,制备生物柴油的反应在恒温下进行。
2传热理论分析
对于本文研究的实验装置锡池为热源(热体),反应釜与反应液为热阱(冷体)。所涉及的传热方式主要为热传导。描述热传导的机理及数学模型为傅里叶(Fourier)定律。控制方程为三维不可压缩流动的连续性方程和Navier-Stokes动量及能量方程:
式(5)中,流体的比能E定义为单位质量流体的内能和动能之和,E=i+(u2+v2+w2)/2。对于不可压缩流体,内能i=CVT,CV是定容比热容, ᐁ·u=0。于是能量方程可以写成下面的温度方程:
FLUENT采用有限体积法(FVM)求解传热用的是如下形式的能量方程和N-S方程:
式中,keff=kt+k,为有效导热系数,Tref=298.15 K。
在两相交界面处引入耦合边界条件:
这里假设区域Ⅱ为流体,区域Ⅰ为固体。
3计算与模拟
3.1计算模型
为合理设置边界条件,简化计算过程,对传热过程的整体模型作以下假设:
(1)锡池整个外表面不对空气辐射传热,即没有热量损失;
(2)反应液充满反应釜的整个腔体,固-液接触表面没有间隙;
(3)各流体为不可压缩流体,即流体密度不随时间变化。
采用Fluent6.2软件求解传热过程温度场,以反应装置的实际尺寸作为计算的物理模型,选择锡池外表面包含的区域为计算域。由于实验中采用的是超临界甲醇法来制备生物柴油,所以在此以甲醇作为模拟反应釜内的升温过程的流体,模拟的最高温度为573 K。网格划分如图2,锡池的节点数为298 089个,反应釜的节点数为100 214个,反应液的节点数为36 432个。
采用SIMPLE算法求解压力—速度耦合方程,动能和能量的离散格式均采用一阶迎风差分格式,反应釜和各流体壁面均采用标准壁面函数方法处理,材料相关参数来自文献[7,8]。
3.2计算结果与讨论
由非稳态解可知,在经过大约89s时间的传热后,锡池、反应釜和反应液三者温度达到平衡,此时沿中心轴所在对称面得到的温度场分布云图如图3所示。由于前3s时间所求解不稳定,所以截取了3s以后各时刻的锡液传热速率随时间的变化,如图4所示。从锡液传热速率的变化曲线可以看出,前30s内明显高于后段传热时间的传热速率,随着锡池、反应釜和反应液三者温差的降低,传热速率降低,传热过程满足热力学第一定律。将模拟的温升曲线和试验测得的温升曲线进行了对比,如图5,模拟的升温过程和实验数据的误差在8%以内。
甲醇的超临界温度为512.6 K,着重考察了反应釜腔体内流体在512.6 K温度附近的温度场分布,且采取了强化传热的措施,将锡池设为流动状态,图2中的B表面为流体入口,A表面为出口,尺寸减小到40 mm×128 mm×40 mm,以节省锡的用量和减小反应装置的体积,对比了两种传热过程。两种情形反应装置温度场分布对比方图如图6,沿反应液中轴线X-Y图分别为图7及图8所示。其中X-Y图所取的中轴线分别处于X方向15—82 mm和10—77 mm从仿真结果可以得到以下结论。
(1)原来的传热方式在18 s左右使流体温度达到临界温度,温度分布在(516.39—534.7) K,而后者只须8 s左右,温度分布在(520.69—521.69) K,有效地缩短了传热时间和改善了温度分布的不均匀性。与有关文献[9,10]报道的流体在超临界温度下呈现温度不均匀分布的特点相符。
(2)观察图6可以知道反应釜腔体内的流体温度要比左边反应釜壁面温度低,这是因为金属吸热快的缘故。(a)图中流体温度低于反应釜密封端的温度,(b)图中流体温度高于反应釜密封端的温度。由于两种情形传热时间的差别,而反应釜密封端金属较多,前者经过时间较长,密封端吸收了足够多的热量,温度升高较流体快;后者经过时间短,密封端吸收热量较少,温度低于流体。这样可以减少反应釜吸收的热量,减少锡液温度的波动,有利于控制反应温度。
4结论
对某制备生物柴油反应装置传热进行了数值模拟,得出了固-液耦合相各时刻的温度场分布,模拟结果与实验结果的吻合度较高,验证了模型的适应性。着重分析了甲醇在超临界温度下的温度场,模拟结果表明此时反应釜腔体内流体温差较大,即温度场分布不均匀。为了改善这种温度场分布的不均匀性,对该装置采取了强化传热措施的升温过程进行了仿真。对比两种传热方式后发现,传热方式对反应釜腔体内的流体温度影响极大。强化传热可以有效地改善温度分布的不均匀性并缩短传热时间,进而减少了反应釜吸收的热量,可以更充分地利用能量。
超临界状态下的反应机理和传热过程十分复杂,要深入研究超临界状态反应过程中的传热还需要创立新的计算方法,手动编程来实现能够真实模拟考虑反应吸热和放热特性的计算机算法。因此将这种自主开发的传热模型的计算法与CFD软件相结合,进一步提高模型的精确性、通用性及仿真度,以得到更多模拟数据来指导此类装置的设计、运行和改造,并在揭示反应机理方面充当更重要的角色, 是今后工作的主要研究方向。
参考文献
[1]Yin Jianzhong,Xiao Min,Song Jibin.Biodiesel from soy bean oil su-percritical methanol with co-solvent.Energy Conversion and Manage-ment,2008;49:908—912
[2]Kusidiana D,Saka S.Two-step preparation for catalyst free biodiesel fuel production hydrolysis and methyl esterification.Appl Biochem Biotech,2004;115(1—3):78l—791
[3]Shimada Y,Watanade Y,Samukawa T,et al.Conversion of vegetable oil to biodiesel using immobilized Candida Antarctica lipase.J Am Oil Chem Soc,1999;76:789—793
[4]Ma F,Clements L D,Hanna M A.Biodiesel fuel form animal fat.an-cillary studies on transesterification of beef tallow.Ind Eng ChemRes,1998;37(9):3768—3771
[5]Jayakumar J S,Mahajani S M,Mandal J C,et al.Experimental and CFD estimation of heat transfer in helically coiled heat exchangers.Chemical Engineering Research and Design,2008;86:221—232
[6]Jafari A,Tynjala T,Mouswi S M,et al.Simulation of heat transfer in a ferrofluid using computational.Heat Fluid Flow,2008;7:1—6
[7]李春胜,黄德彬.金属材料手册.北京:化学工业出版社,2004:56—61
[8]张书圣,温永红,丁采风,等.有机化学手册.北京:化学工业出版社,2006;6:499—521
[9]李会雄,孙树翁,郭斌,等.超临界水在倾斜光管中的传热不均匀性研究.工程热物理学报,2008;2(29):241—245