大涡数值模拟范文
大涡数值模拟范文（精选7篇）
大涡数值模拟 第1篇
关键词：大涡模拟,导叶,攻角,叶片尾流,水轮机
水轮机流道形状和水流运动的复杂性, 使得水力振动难以把握。导叶是水轮机里重要的水力部件, 它调节进入转轮液流的速度环量。导叶水力性能的优劣直接影响到转轮内的流态, 从而影响水轮机的效率和稳定性。活动导叶开度的变化水流角度发生变化, 必然对叶片压力面或吸力面进行冲击, 引起脱流漩涡降低效率, 甚至进一步发展形成叶道涡, 产生中高频的压力脉动, 对水轮机稳定工作构成威胁。本文拟采用大涡模拟方法 (LES) 对某一具体工程的活动导叶进行数值模拟, 观测导叶流道中速度、压力、湍动能的分布以及对叶片尾流的运动机理进行分析。望能为探索影响水力机组的水力振动因素、改进叶型设计和提高整机效率等提供参考。
1控制方程
在LES方法中对变量Φ使用滤波函数
undefined
式中, undefined为大尺度的平均分量;D是流动区域;x′是实际流动区域中的空间坐标;x是滤波后的大尺度空间坐标;G (x, x′) 是滤波函数, 它决定了所求解的涡的尺度。
用式 (1) 表示的滤波函数处理瞬时状态下的Navier-stokes方程及连续方程有
undefined
undefined
其中undefined为亚格子尺度应力。
根据Smagorinsk的基本SGS应力模型有
undefined
式中, μt为亚格子尺度的湍动粘度。
undefined
式中, Δi代表i轴方向网格尺寸, Cs是Smagorinsky常数。
取Cs=Cso (1-ey+/A+)
式中, y+是到壁面的最近距离, A+ 是半经验常数取25.0, Cso是Van Driest常数取0.1。
2数值模型及边界条件
在时间离散上, 采用二阶全隐式格式。在空间上采用有限体积法对控制方程进行离散, 对流项采用二阶迎风格式, 扩散项采用中心差分格式。每个时间步上, 应用SIMPLEC方法进行离散方程的求解。在本计算中, 时间步长均选为0.000 5 s。
进口面上给定速度ui和通过计算得到的湍流强度k与扩散率ε。由于入口边界实际速度分布难准确测定, 所以将入口区域适当延长, 按充分发展匀速流给定, 即取平均速度;出口面上按零梯度处理, 出口距离适当延长采用自由出流的边界条件;壁面边界处应用无滑移边界条件;平行进口水流方向两侧边取周期边界条件;不计重力, 入口处压力暂取为0.1 MPa。
对某实际工程中的导叶进行三维数值分析, 模型雷诺数为4.2106, 计算湍流强度k=0.339 3, 扩散率ε=0.736 4, 划成30万个网格, 并在页面处和壁面处局部加密网格。分别在攻角为0°, 14.81°, 29.62°三种有代表性的开度情况下, 观察导叶流道中的速度流场分布、湍动能分布和导叶表面的压力系数分布情况。
3计算结果分析
经计算得到模型各速度残差值在10-4左右摆动, 不大于10-3。计算的进出口流量差最大为410-7, 固所得结果收敛。最终模拟结果如图1-9所列。
由于速度较大, 从图1中可看出在导叶开度为零时, 叶片后方出现少量湍流涡, 其能量不大, 随流动直线向后扩散至转轮。最大湍流粘度2.876 Pas。
可以看到导叶两侧出现明显的负压区其最小负压为-187 969.4 Pa。最大压力出现在导叶顶端其值为185 120.9 Pa。
计算得流场中最大速度为27 m/s, 由于导叶尾端上部有切角, 故速度矢量图上低速区向上方移动。
攻角为14.81°时导叶头部后方分离的湍流能量较大, 最大亚格子湍流粘度4.55 087 Pas。尾部的尖端湍流能量区也较大, 尾端湍流脱落近似卡门涡的脱落。造成导叶尾端周期震动。
导叶正面形成高压区, 背面出现负压区。最大压力出现在导叶头部值为190 807 Pa, 负压区最小压力为-454 670.1 Pa。
导叶两侧速度比前后端增大, 导叶负压面起始处速度最大, 值为33.298 m/s 负压面出现明显回流现象, 其最大速度为-14.796 m/s。
随着攻角增大到29.62°时, 负压面上边界层分离点前移, 边界层提前脱落, 导叶头部后方湍流能量增大并向后传递, 导叶尾部漩涡湍动能值最大, 其最大亚格子湍流粘度为7.369 9 Pas。在尾流区内, 类似于卡门涡列的扰流漩涡现象明显, 并先按导叶方向脱落, 脱落后沿流道方向向后传递。
叶片正压面上形成一个高压力区, 最大压力出现在导叶头部顶端, 值为197 718.4 Pa, 叶片负压面的最低压力区最小压力为-1 099 080 Pa。
从速度图中可看出叶片后方回流现象严重, 由于湍流影响x方向流体最大速度43.79 m/s 导叶负压面下回流最大速度值为-22.15 m/s。
从模拟结果可以看出在叶片两侧及尾部, 叶片正压面的相对流速较大与叶片负压面回流漩涡间的速度差造成垂直于叶片表面方向的较大速度分量, 该速度分量的方向在叶片尾部周期交替变化, 形成周期性扰流。还可得出导叶正压面压力系数随攻角的增加而增大, 负压面压力系数随攻角的增加而减小, 并且负压面上回流现象随着攻角的增大越趋明显。随攻角增大导叶尾部湍流能量逐渐增大, 在大攻角时, 出现周期性类似卡门涡街的湍流涡脱落现象, 造成导叶周期性震动。
在大攻角时导叶尾部出现湍动能极大值、速度极小值、压力极小值。绝对速度最大值、负压力极大值, 尾部湍流能量最大值等均在最大攻角时出现。这是由于负压面上边界层分离点前移, 边界层提前脱落, 形成尾迹涡, 如果流入转轮流道, 会引起水力振动。
4结语
(1) 片机械在大流速高雷诺数下叶片绕流流态十分复杂, 大攻角时在叶片负压面存在明显回流, 漩涡脱落, 压力降低甚至造成空蚀现象。
(2) 叶攻角较大时, 在叶片尾部区域湍动能值较大并在尾流区存在类似卡门涡街的湍流现象, 其周期性间隔脱落是造成导叶尾部周期震动的原因。
(3) 湍流涡带进一步向下游发展, 造成转轮进口速度及压力分布的明显不均匀, 将直接影响转轮内部的流动状态。
(4) 由于导叶尾部切角, 随着导叶攻角增大, 导叶尾部负压区前移。且整个负压区的范围渐渐扩大, 负压面上速度最大值点也明显提前。
参考文献
[1]王文全, 张立翔, 闫妍, 等.混流式水轮机导叶叶道内湍流场的大涡模拟[J].热能动力工程, 2009, 24 (2) :257-260.
[2]曹昆鸟, 姚志民.水轮机原理及水力设计[M].北京:清华大学出版社, 1991.
[3]王福军.计算流体动力学分析[M].北京:清华大学出版社, 2004.
[4]杨昌兵, 杨永泉.对称导叶在不同攻角下的水力特性数值分析[J].四川大学学报, 2002 (6) :28-31.
[5]张立翔, 王文全, 姚激.混流式水轮机转轮叶片流激振动分析[J].工程力学, 2007, 24 (8) :143-150.
[6]刘春梅.混流式水轮机转轮内部流动的CFD分析:[硕士论文].四川:西华大学, 2006.
[7]张昌兵.水力机械叶片尾流机理研究[J].四川大学学报, 2005 (6) :15-19.
[8]肖若富, 王正伟, 罗永要.基于流固耦合的混流式水轮机转轮静应力特性分析[J].水利发电学报, 2007 (3) :120-123.
[9] Zhang Lixiang, Wang Wenquan.LES of turbulent flowin 3Dskewbladepassage of a reacting hydro turbine[J].ModemPhysics Letters B, 2005, 19 (28.29) :1487-1490.
大涡数值模拟 第2篇
R-T界面不稳定性及湍流混合的大涡模拟
采用大涡模拟(LES)方法研究R-T(Rayleigh-Taylor)不稳定性,显示了不稳定发展的小扰动阶段、变形阶段、规则非线性阶段、不规则非线性阶段和湍流混合阶段界面演化与特征,分析了湍流混合层厚度增长的规律,同时还研究了湍流二阶相关量如动量、质量和能量通量在大尺度和亚格子尺度所占分量,验证了LES方法是处理R-T界面不稳定性和湍流混合的`一种有效方法.
作 者：谢正桐 李家春 王丽丽 作者单位：中科院力学所工程科学部,中科院非线性力学国家重点实验室(LNM)北京100080;北京应用物理与计算数学研究所,北京100080刊 名：计算物理 ISTIC EI PKU英文刊名：CHINESE JOURNAL OF COMPUTATIONAL PHYSICS年，卷(期)：19(2)分类号：O357.5 O242.1关键词：界面不稳定性 湍流混合 大涡模拟
大涡数值模拟 第3篇
关键词：大涡模拟；谐波合成；数值模拟；脉动风场
中图分类号：O35 文献标识码：A
文章编号：1674-2974（2015）11-0064-08
大涡模拟在计算资源与计算效率方面是雷诺时均模拟（RANS）和直接模拟（DNS）的折中，集成了二者优势，因此在CFD模拟方面得到了广泛应用.但在实际运用过程中，大涡模拟入口边界处的脉动信息给定问题还没有得到完全解决[1].目前，许多学者直接将均匀来流赋给入口边界，没有考虑脉动信息的作用.而在实际风场中，来流脉动对结构的脉动风压和瞬时风荷载起着非常重要的作用[2]，有着不可忽略的地位.因此，非常有必要在大涡模拟过程中考虑脉动信息的影响.
针对大涡模拟脉动入口研究方法，可以将其分为三类[3]，第一类为“预前模拟法”，是运用一个预前模拟区域，将目标风场预先模拟出来，然后赋给主模拟区域的入口边界.朱伟亮等[2]在预前模拟区域用流场循环的方法生成了脉动风场，并对钝体平面屋盖的结构风压进行了模拟，这种脉动风场生成方法的不足之处是需要生成匹配的数据库，会消耗巨大的内存和计算时间.序列合成法的另一种储存形式是用空间序列储存，Chung[4]利用泰勒方法协调了预前模拟方法两个计算域之间的离散差异性，对槽道流进行了模拟，运用该方法的缺点是不同时刻的速度不能严格的由泰勒级数展开来获取，因此插值的效果不佳.第二类方法为涡方法，针对涡方法的研究，Mathey等[5]人用一系列的二维随机漩涡加载到平均风场中去从而得到脉动风场，这种方法生成的脉动风场虽然具有较好的空间相关性和湍动能信息，但不足之处是目标风场很难满足风场各向异性目标谱的要求.第三类方法为序列合成法，序列合成法是运用傅里叶变换的方法来生成脉动风速时程.2002年Hanna等[6]人提出了一种简单的各向异性湍流脉动生成方法，生成了随机正态序列数，它是假设在顺风方向和垂直方向的脉动均方根和平均速度有不同的比例，然后通过加权得到脉动信息，这种方法的不足是缺少了雷诺应力信息，同时湍流也会很快衰减.此后，2006年，Kondo[7]等提出了动态湍流输出方法，但这种方法不能满足目标谱的需求.2008年，Xie和Castro[8]用雷诺应力的垂向分布构造出了脉动速度分量，得到的脉动量满足湍流的基本特性并具有雷诺应力信息，但是生成的湍流在入口处并不满足连续性的需求，需要经过很长时间的发展才能形成真正的湍流.后来，Hemon and Santi[9]，Huang等[10-11]也针对大涡模拟的入口脉动边界做出了相关研究，得到了不错的计算结果.
本文借助于经典的谐波合成方法，在充分考虑湍流的平均风剖面、目标谱、空间相关性等指标下，合成一系列的随机序列时程数据，把生成的时程数据与大涡模拟的入口建立对应关系，基于对FLUENT软件平台进行二次开发，把生成的随机序列数赋给大涡模拟的入口边界，从而生成入口脉动信息.最后通过对比数值模型和风洞试验结果发现：
基于谐波合成的脉动信息生成方法能满足大涡模拟入口边界的要求，其优势在于目标谱及相关性参数能根据不同情况的需要进行随意调整，合成的时程数据可直接赋给入口边界，生成速度快，模拟时间长短可以随意控制，模拟的点都相互独立，可以很好地进行并行计算，大大提高计算效率.
1 数值模拟理论与方法
在充分考虑脉动风场的平均速度、功率谱、相关性、时间步长等参数后，基于谐波合成理论，用MATLAB对风速时程数据进行合成，生成满足目标风场的随机序列数.然后将生成的随机序列数与入口网格在时间和空间上建立对应关系，时间上要满足时间步长的一致，保证谐波合成风速与大涡模拟风速在时间上的同步.空间上要对入口网格坐标进行严格控制，将模拟点的速度时程赋给对应的入口网格中心点，整个过程用UDF程序编制，对FLUENT软件用户自定义模块进行加载，重新定义入口速度，从而实现大涡模拟入口脉动信息的输入.
2 大涡模拟入口脉动风场特性分析
为验证本文入口处脉动信息输入的正确性，建立了一空风洞数值模型，里面没有任何障碍物，具体尺寸为15 m×2.5 m×3.0 m （x×y×z），总网格数为180万.入口平均风采用指数率风剖面形式，α值为0.16，取基准风速为20 m/s，地表粗糙高度为0.05 m.入口脉动量由MATLAB程序模拟生成，然后通过UDF程序将这些时程数据赋给入口边界网格，待计算稳定后提取模型入口0.5 m和1 m高度处的风速时程并与输入的谐波合成法生成的数据序列对比，如图1～2所示.
由图1～2可以看出，模型入口监测值与输入的谐波合成法生成的数据序列基本一致，说明了入口边界风速时程输入方法的正确性.为了验证入口边界功率谱的正确性，将数值模型入口1 m高度处顺风方向功率谱模拟值与我国规范进行对比，如图3所示.
从图3可以看出，功率谱能较好的与我国规范采用的Kaimal谱吻合，满足脉动风场特性要求.
对流场入口中心0.5 m和1 m高度处的速度时程数据进行相关性分析，定义两点分别为A和B.图4给出了入口监测点A和B的互相关性和B点的自相关性检验结果，图4（a）中目标值由公式（11）给出.
图4表明了入口处脉动风场相关性呈指数率衰减，模拟值与理论目标值衰减趋势一致，体现出了模拟风场在入口处具有良好的相关性.
通过对上述平均风速、功率谱、相关性等参数进行分析，发现基于谐波合成法生成的脉动风场满足湍流的基本特性，可以在大涡模拟的入口进行无缝对接.
3 数值模型
为验证数值模型内部风场正确性，以TJ-2风洞试验数据[18]为参照对象，建立了两个数值模型，一为没有任何障碍物的空风洞，与第二节所用空风洞验证模型为同一计算模型，利用上节所提出的谐波合成方法来生成入口脉动信息，简称为“空风洞数值模型”；二为与风洞试验条件一致的“尖劈粗糙元数值模型”.下面分别对两种模型的尺寸，网格数量，计算参数与边界条件进行介绍.
3.1 模型尺寸与网格数量
3.1.1 空风洞数值模型
TJ-2风洞具体尺寸为15 m×2.5 m×3.0 m （长×高×宽），空风洞数值模拟采用与物理风洞一致的模型尺寸，为保证计算精度，采用全六面体网格，并对总网格数为90万、180万、400万3种网格进行了网格无关性测试，权衡计算资源和计算精度后，选取了180万网格数模型作为最终计算模型，网格在模型底部进行加密，增长因子为1.1，最底层网格高度为0.01 m，具体网格如图5所示.
3.1.2 尖劈粗糙元数值模型
尖劈粗糙元数值模型与物理风洞试验条件保持一致，布置如图6所示.与上一节类似，为保证精度，对尖劈粗糙元数值模型进行全场六面体网格划分，并对260万，450万，800万3种网格进行了网格无关性测试，无关性测试结果如表1所示.
通过比较模型内部测量中心处湍流度发现260万网格数模型较450万和800万的偏差较大，但450万和800万网格数的计算结果偏差幅度很小，权衡计算精度和计算资源的影响后，最终选取了450万网格数作为最终计算网格.整个网格在底部进行加密，最底层网格高度为0.005 m，延伸率为1.1，网格如图7所示.计算过程中，对空风洞数值模型和尖劈粗糙元数值模型的“预定模型中心”位置进行风速监测，具体位置如图6所示.
3.2 边界条件设置
空风洞数值模型的入口条件是基于对FLUENT软件平台进行二次开发，将生成的随机序列数赋给入口边界；而尖劈粗糙元数值模型则采用平均风作为入口边界，详细边界条件与计算参数如表2所示.
4 计算结果
空风洞数值模型采用超线程12核工作站进行计算，模拟150 s时长（60 000步）需耗时60 h，其计算结果如图8（a）所示.尖劈粗糙元模型同样采用超线程工作站进行计算，模拟同样时长需耗时124 h，计算结果如图8（b）所示.
通过对 “预定模型中心”不同高度处的风速监测，得到监测中心处的平均风剖面与湍流度信息，对其无量纲处理并与风洞试验结果进行对比，其结果如下图所示.通过比较图9和图10发现，在风剖面和湍流度方面，两种数值模型相对风洞试验值的偏差均较小，能满足工程精度需要，也再一次证明了谐波合成方法生成脉动风场的正确性.用最小二乘法对空风洞模型风剖面进行指数率拟合，得到α指数为0.161，非常接近风洞试验的拟合值0.162，也接近于我国规范给出的B类风场特性.
对“预定模型中心处”1 m高度处的时程数据进行分析，得到空风洞模型和尖劈粗糙元模型的顺风方向功率谱，将其进行拟合并与Kaimal谱进行对比，如图11所示.
通过对两种方法同一监测点顺风方向功率谱与Kaimal谱进行对比发现，数值模拟结果与风洞试验的风速功率谱在低频段（频率<1）吻合较好，但在高频段（频率>1）出现陡降，而工程领域最为关心的频率出现在惯性子区（0.2<频率<1），本文两种数值模拟方法均能较好的捕捉该区域的风谱，因此所生成的脉动风能较好的满足工程需要.目前，数值模拟的风速时程普遍存在风谱高频段衰减现象，是由于高频段代表风场中小涡的能量贡献，而LES中小尺度涡是用亚格子模型进行封闭，不能直接对小尺度涡进行求解，因此数值模拟不能捕捉到脉动风场的高频段风谱，通过改进湍流模型和加密网格数量可以改善风谱高频衰减问题.
在计算效率方面，将两种数值模型所需的网格数和耗时量作出对比，如表3所示.
由表3可以发现，在达到允许精度要求下，空风洞模型所需的网格数不到尖劈粗糙元模型的一半，计算耗时量方面也远远优于尖劈粗糙元模型.
此外，空风洞数值模型具有更好的适应性，特别是对不同脉动风场的模拟，空风洞模型只需对目标风场的特性参数进行调整，操作简单方便，可调性强.而尖劈粗糙元数值模型只能通过不断调整尖劈粗糙元的位置和数量来对目标风场进行试算，计算工作量大，消耗时间长.
5 结 论
本文基于谐波合成方法生成了随机序列数，通过对商业软件FLUENT进行二次开发，利用UDF程序把时程数据转化成大涡模拟的入口脉动信息，同时建立了与物理风洞试验条件一致的尖劈粗糙元数值模型，通过对比研究得到了以下成果和结论.
1）实现了随机序列数在大涡模拟入口处的无缝对接，开发了大涡模拟入口脉动速度输入模块，可以很好地进行并行计算，解决了大涡模拟入口处湍流信息给定问题.
2）通过比较平均风剖面、湍流度和功率谱等结果发现：基于谐波合成方法生成的大涡模拟脉动风场具有较高的精度，能满足工程需要，但风速功率谱在高频段会出现衰减现象.
3）通过对比空风洞数值模型和尖劈粗糙元模型的计算效率发现，空风洞数值模型在网格数量和计算耗时方面均要远远优于尖劈粗糙元数值模型.
4）相比尖劈粗糙元数值模型，基于谐波合成脉动生成方法对不同风场模拟只需调整风场特性参数，操作方便，模拟速度快，是一种极具前景的脉动生成方法，也可为大涡模拟入口处脉动生成方法精细化研究提供参考.
基于谐波合成法的脉动风场生成方法是谱合成方法的一种拓展，拥有许多优点的同时也存在一些不足，如生成的随机风速时程不满足N-S方程的连续性，风谱高频衰减等问题还有待进一步解决.
参考文献
[1] JIANG Wei-mei，MIAO Shi-guang. 30 years review and perspective of the research on the large eddy simulation and atmospheric boundary layer[J]. Advances in Nature Sciences， 2004， 14： 11-19.
[2] 朱伟亮，杨庆山. 基于 LES 模型的近地脉动风场数值模拟[J]. 工程力学，2010，27（9）：17-21.
ZHU Wei-liang， YANG Qing-shan. Large eddy simulation of near ground turbulent wind field [J]. Engineering Mechanics， 2010， 27 （9）： 17-21. （In Chinese）
[3] BABA-AHMADI M H， TABOR G. Inlet conditions for LES using mapping and feedback control[J]. Computer & Fluid， 2009， 38（6）： 1229-1311.
[4] CHUNG Y M， SUNG H J. Comparative study of inflow conditions for spatially evolving simulation[J]. AIAA Journal， 1997， 35（2）： 269-274.
[5] MATHEY F， COKLJAT D， BERTOGLIO J P， et al. Assessment of the vortex method for large eddy simulation inlet conditions[J]. Progress in Computational Fluid Dynamics， 2006， 6： 1-3.
[6] HANNA S R， TEHRANIAN S， CARISSIMO B， et al. Comparisons of model simulations with observations of mean flow and turbulence within simple obstacle arrays[J]. Atmospheric Environment， 2002， 36： 5067-5079.
[7] KONDO H， ASAHI K， TOMIZUKA T， et al. Numerical analysis of diffusion around a suspended expressway by a multi-scale CFD model[J]. Atmospheric Environment， 2006， 40： 2852-2859.
[8] XIE Z T， CASTRO I P. Efficient generation of inflow conditions for large eddy simulation of street-scale flows[J]. Flow Turbulence and Combustion， 2008， 81： 449-470.
[9] HEMON P， SANTI F. Simulation of a spatially correlated turbulent velocity field using biorthogonal decomposition[J]. Journal of Wind Engineering and Industrial Aerodynamics， 2007， 95 （1）： 21-29.
[10]HUANG S H， LI Q S， WU J R. A general inflow turbulence generator for large eddy simulation[J]. Journal of Wind Engineering and Industrial， 2010， 98： 600-617.
[11]卢春玲，李秋胜，黄生洪，等. 大跨度屋盖风荷载的大涡模拟研究[J]. 湖南大学学报：自然科学版，2010，37（10）：7-12.
LU Chun-ling， LI Qiu-sheng， HUANG Sheng-hong， et al. Large eddy simulation of wind loads on long-span roofs[J]. Journal of Hunan University ： Natural Sciences， 2010，37 （10）： 7-12. （In Chinese）
[12]DEARDOFF J W. The use of sub grid transport equation in a three dimensional model of atmospheric turbulence[J]. ASME Journal of Fluid Engineering， 1973， 95： 429-438.
火羽流轴向温度大涡模拟与实验比较 第4篇
1 经典理论模型
1.1 Zukoski模型
Zukoski进行了不同面积、不同形状的池火灾实验, 对理想羽流模型进行了一定的优化修正, 得出了计算池火灾羽流中心线轴向温升的计算公式, 如式 (1) 所示。
式中:ΔT0为火羽流轴向温升, K;T∞为环境温度, K;g为重力加速度, 取为9.81m/s2;Qc为对流热, kW, 取Qc=0.7Q;Z为距离火源的高度, m;cp为空气的定压比热容, kJ/ (kgK) ;ρ∞为环境空气密度, 取1.29kg/m3。
1.2 Heskestad模型
Heskestad引入“点火源”的概念, 并取消了理想羽流模型的无辐射假设, 总结分析了大量实验数据, 得出了计算池火灾羽流中心线轴向温升的计算公式, 如式 (2) 、式 (3) 所示。
式中:Z0为虚点源的高度, m;Q为火源的热释放速率, kW;D为池火源的直径, m。
2 大空间火羽流轴向温度实验与模拟
2.1 大空间火羽流轴向温度实验
实验在江苏省徐州市黄山工业园大型厂房内进行, 其平面图及剖面图如图1所示。厂房为单层钢结构, 东西长60m, 南北宽30m, 高9.2m, 建筑面积为1 800m2, 厂房内部室温为284K。
实验以柴油为火源, 进行了火源功率为640、1 080kW的两组实验。在火源上方4.5~7.5m处每隔0.5m设置一个K型热电偶, 连接Agilent34970A数据采集器, 以获得火羽流中心线轴向温度的数据。
由于池火灾在初始阶段火源功率快速增长, 火羽流中心线轴向温度不稳定, 随着火源功率的增加会达到稳定燃烧阶段, 火羽流中心线的温度会趋于稳定。因此, 选取300~700s的温度数据作为研究对象, 图2、图3为300~700s火源上方不同高度热电偶记录的数据。
2.2 基于大涡模拟的火羽流轴向温度模拟
2.2.1 模型建立
建立全尺寸火灾计算模型, 模型长度为60 m, 宽为30m, 高为9.2m, 建筑面积1 800m2。由于模拟精度与网格大小有着密不可分的关系, 因此将模型网格大小设置为0.5m0.5m0.5m, 以确保能精确模拟池火灾发生及发展的全过程, 准确测量火羽流中心线轴向温度。模型的计算条件:火源位于模型中心部位, 热释放速率为640、1 080kW, 环境温度设置为284K, 并在火源的上方4.5~7.5m每隔0.5m布置一个热电偶, 测量火羽流中心线轴向温度。
2.2.2 模拟结果分析
为了便于研究, 选取300~700s的火羽流中心线轴向温度为研究对象, 不同位置热电偶记录了火羽流中心线轴向温度数据。图4、图5为300~700s火源上方不同高度热电偶数据。
2.3 模拟数据与实验及理论计算对比
利用Zukoski、Heskestad经典羽流模型及式 (1) ~ (3) 计算大空间池火灾4.5~7.5m处羽流中心线轴向温度。假定环境空气密度为1.29kg/m3, 比热容为1.04kJ/ (kgK) , 环境温度为284K。基于大涡模拟的池火灾发展全过程模型中热电偶测量火源上方中心线轴向温度, 选取300~700s的温度数据作为研究对象, 将记录温度的平均值作为火羽流中心线轴向温度, 并与全尺寸实体池火灾数据及经典羽流模型的预测值进行分析比较, 如图6、图7所示。
由图可知, 基于大涡模拟的FDS模拟得到的火羽流中心线轴向温度与实验数据变化趋势基本吻合, 均随着高度的升高而降低。与经典羽流模型预测的数据相比, FDS模拟得到的火羽流中心线轴向温度在火源上方4.5~5.0m范围内吻合度最高, 相对误差最小, 随着高度的增加相对误差逐渐增大。同时, FDS模拟得到的火羽流中心线轴向温度比实验及经典羽流模型数据普遍偏大, 但能保持在一个较为合理的误差范围之内。
为了进一步分析基于大涡模拟的FDS技术在预测大空间池火灾羽流中心线轴向温度的准确性, 对模拟值和经典羽流模型值与实验值进行处理分析, 图8、图9是640、1 080kW火源中心线轴向温度相对误差比较图。
由图可知, 基于大涡模拟的FDS技术和经典羽流模型都能较好地预测大空间小功率池火灾中心线轴向温度, 其相对误差保持在0.249%~8.140%。在火源上方4.5m处, 模拟值的相对误差均为最小, 有1.64%、0.249%, 随着距离高度的增加, 其相对误差有增大趋势, 但最大值也只有4.64%、5.81%, 均在允许的相对误差范围之内。同时, 其相对误差的平均值仅为3.10%、3.52%, 说明与实验数据的吻合度较高, 能有效预测火羽流中心线轴向温度。与经典羽流模型的预测值相比, 基于大涡模拟的模拟值相对误差普遍偏小, 仅在个别点出现相对误差偏大的现象, 其相对误差的平均值最小, 说明基于大涡模拟的FDS技术在预测大空间小功率池火灾中心线轴向温度方面要优于经典羽流模型。
3 结论
(1) 基于大涡模拟的FDS模拟得到的火羽流中心线轴向温度与实验数据变化趋势基本吻合, 且与经典羽流模型预测的数据相比, FDS模拟得到的火羽流中心线轴向温度在火源上方4.5~5.0m范围内吻合度最高, 同时模拟值比实验及经典羽流模型数据普遍偏大, 但能保持在一个较为合理的误差范围之内, 能确保应用到工程中的可靠性、安全性。
(2) 基于大涡模拟的模拟值的相对误差比经典羽流模型的预测误差普遍偏小, 仅仅在个别点会出现相对误差偏大的现象, 同时其相对误差的平均值最小, 说明基于大涡模拟的FDS技术在预测大空间小功率池火灾中心线轴向温度方面要优于经典羽流模型。
摘要：为了验证基于大涡模拟的FDS技术研究火羽流中心线轴向温度的可靠性, 建立FDS模型对某实体厂房中640、1 080kW的池火灾进行数值模拟, 并与实验室开展的全尺寸池火灾数据及经典理论模型预测值比较。研究结果发现, 基于大涡模拟的FDS技术得到的火羽流中心线轴向温度与试验数据变化趋势基本吻合;与经典羽流模型预测的数据相比, FDS模拟得到的火羽流中心线轴向温度吻合度最高, 其相对误差值普遍偏小, 仅仅在个别点会出现相对误差偏大的现象。
关键词：大涡模拟,火羽流,轴向温度,池火,FDS
参考文献
大涡数值模拟 第5篇
现代桥梁的跨径日益增大,风荷载也渐渐成为此类大跨度桥梁的主要控制荷载。对于桥梁风工程领域,其研究仍是以风洞试验为主。由于风洞中基于相似理论的模型实验有其自身的弱点,例如缩尺模型雷诺数模拟的困难,实验周期长,资金消耗多,实验设备复杂以及流场可视化困难等。CFD技术的发展为桥梁风工程提供了一种风洞实验之外的研究手段,其优点在于投资少,周期短,节省人力物力,可控性强,在计算机上具有很好的重复性和可视性。同时还能实现风洞缩尺模型实验中难以实现的条件。因此,配合日益高效的计算设备,CFD技术已得到了非常广泛的应用,得到了迅速的发展[1,2]。
目前的CFD模拟的方法主要是基于N-S方程的流体计算方法和基于粒子运动的LB方法[3]。N-S方法中根据流场的求解方法和对于计算系统的要求大致可分为RANS,DNS和LES方法。简单而言,RANS方法计算量小,湍流模型多且适用性较差,只有流场的统计信息,流场的细节被全部抹去,多用于工程统计,很多数值研究者将其应用在桥梁风工程方面,如桥梁断面气动三分力及表面压力系数的识别,湍流模型用于不同桥梁断面形式的比选[4],采用动网格技术对桥梁断面的气动导数进行了计算等[5]。DNS方法计算量最大,能获得全部的流场信息,限于目前的计算条件多用于低雷诺数的科学研究,工程应用较少;而LES则介于两者之间,吸引了大量数值研究者对其进行科研和工程应用方面的研究,并被认为是今后用于科研及工程研究的比较理想的数值计算方法。
1 大涡模拟方法
现代湍流理论认为,湍动能传输链中大尺度脉动几乎包含了所用的湍动能,而小尺度脉动主要耗散湍动能。大涡模拟方法基于N-S方程,在湍流数值模拟中仅计算大尺度脉动,小尺度脉动与大尺度的相互作用采用特定的模型来封闭。由于放弃了大量的小尺度脉动的计算,湍流模拟的空间和时间步长都可以放大,从而有效的降低计算资源的消耗[6]。
通过滤波方法可以将小尺度脉动从湍流脉动中去掉,得到可解尺度湍流脉动控制方程,其中含有不封闭的亚格子应力项,需要构造模型封闭。常用的Smagorinsky涡粘模式形式简单,只需确定涡粘系数即可,计算稳定性也较好,但其缺点是耗散过大。此后发展了动态确定涡粘系数的动态Smagorinsky模式,可以适应实际湍流的局部结构,因而克服了涡粘模式耗散过大的问题[7]。
实际应用中,复杂湍流计算域网格的近壁面处理始终是数值模拟的关键问题。沿壁面法线以上的不同距离可将流动分为近壁区(内区)和核心区(外区)。对于核心区可直接采用大涡模型计算,对于近壁区又可分为粘性底层、过渡层和对数律层。粘性底层中粘性力在动量、质量交换中起主导作用;对数律层中湍流切应力占主导地位,流动处于充分发展湍流状态;过渡层介于两者之间,粘性力和湍流切应力相当。对于大涡模拟,较好的网格划分要求近壁面第一层网格处于粘性底层内[8]。
2 计算模型
桥梁风工程中研究的一项基本内容即得到气动三分力系数,进而用于分析桥梁结构在静风荷载作用下的应力与位移等。风洞试验中通过高频天平获得一定风速下的模型表面压力和摩擦力时程,时间平均并无量纲化后得到气动力系数,定义如下:。
式中FD、FL和M为风轴坐标系下的阻力、升力和升力矩,U为来流风速,ρ为空气密度,H、B、L分别为模型高度、宽度及轴向长度[4]。
基于Fluent计算平台,网格划分使用GAMBIT,计算模型基于苏通长江大桥主桥闭口箱梁断面,采用动力Smagorinsky-lilly模式对桥梁断面气动力系数进行识别预测。计算模型按风洞试验模型1/70的几何缩尺比建立,无栏杆等附属物,进行-7°、-3°、0°、+3°及+7°五个风攻角的预测,具体模型尺寸见图2。
三维计算域取为2 m6 m0.115 m(设为XYZ平面),模型中心离入口2 m,距出口4 m,模型特征高度H=0.057 3 m,阻塞度为2.9%(±7度时最大阻塞度为5.2%)。考虑到计算模型为二维节段模型,认为其展向流动是统计均匀的,在一定的展向长度及网格分辨率条件下,展向边界可取为对称边界,这里取展向为两倍的断面特征长度,即2H,展向划分为20个网格。空气密度及空气粘滞度取默认值,计算雷诺数为Re=4.6104,均匀风速入口,与实验风速相同,取U=12 m/s,出口采用压力边界条件,表压设为0,垂直流向的上下两个边界同样也设为对称边界,以加快计算收敛速度。
XY平面计算域网格划分如图3,为分辨近壁面边界层复杂湍流,对模型近壁面网格进行了加密,并通过计算调整使近壁面第一层网格的y+<5,约为1,即保证第一层网格处于粘性底层,其物理高度约在4.510-5 m左右。最终划分的体网格数约为95万,节点数约100万。时间分辨率根据0.01B/U=0.000 48 s(B为模型宽度,U为流向风速)及d/U≈0.000 12 s(d为模型壁面附近流向网格尺度),将时间步长设为0.000 1 s。
对于非定长不可压流动的计算,采用PISO算法可以加快收敛,并可以采用较大的松弛因子和较长的时间步,但大涡模拟方法对时间分辨率要求较高,在不增大时间步的情况下采用PISO并不能节省机时,故本文采用SIMPLEC算法进行计算,空间离散采用大涡方法默认的Bounded Central Differencing方法[8]。
本文在计算三维计算域的同时也进行了二维计算域的计算,其二维网格划分同三维计算域的XY平面。
3 数值结果
本文计算时未采用任何近壁面模型来处理壁面湍流,通过计算调整控制近壁面第一层网格在粘性底层以内,即满足无量纲参数y+<5的当量范围内,图4为通过计算调整得到的近壁面第一层网格y+值,基本满足上述要求。
数值计算推进一定时间步后趋于稳定,继续推进一定的时间步后得到稳定状态下的气动力时程,对其作统计平均并无量纲化后得到气动力系数。表1为0°攻角下的气动三分力系数预测结果,并列举了实验值、RANS[4]计算值等。
可以看到,0度攻角三维计算域的计算值均较实验值偏小,阻力系数与实验值的误差在19%左右,升力系数误差在15%左右,但升力矩系数相差较大。二维网格结果与三维网格结果及实验值相比误差均较大,图5给出了五个攻角的预测结果对比。
三维网格计算的阻力系数与实验值符合较好,在0攻角时略小,-7°时略大;升力系数在-7°时相差最大,误差达到13%;升力矩系数符合较差,在+7°时略大,-7°时略小,误差均在40%左右。二维网格预测结果通常较三维结果略大,其中阻力系数差别达到50%。
大涡模拟方法相比于雷诺时均方法,保留了大尺度脉动的全部信息,图6、图7为计算中的流场压力等值线图,可以看到比较丰富的流场信息。三维计算域在展向也具有一定的速度,图中给出的是流场中间剖面的压力等值线。二维计算域由于没有展向的流动,其表面压力漩涡比三维计算域的更明显。
4 计算域及网格尺度影响
大涡计算方法对模型近壁面网格要求较高,下文降低了桥梁模型近壁面网格密度,其计算域及展向网格尺度与上文模型相同,近壁面第一层网格高度增加了一倍以上,其y+值约为5到30之间,总的计算网格数约为39万,节点约41万。计算时在近壁面采用了近壁面平衡层模型进行处理,具体为Wenner等[9]建议的用幂函数作为平衡层的速度分布模型。
由表2可见,近壁面网格高度增大后,算例2的计算结果略差,阻力系数较小,升力系数绝对值较上文算例1要大,与实验值的误差均有所增加。
风洞试验中,模型的长宽比较大,用于减小实验时三维端部效应对二维绕流结果的影响,数值计算中,模型的轴向计算域长度及展向网格尺度同样是重要参数。前文的计算采用2H的轴向计算域长度及展向20个网格,计算结果与实验值相比略偏小,表3给出了轴向计算域长度及展向网格分辨率小于上文模型的计算情况。
可以看出,相同的网格尺度,随着展向计算域的增加,阻力系数和升力系数变化不大,升力矩系数较小,计算值都在零附近。同时,相同的展向长度,展向网格分辨率越高,阻力系数越小,升力系数变化不大,升力矩系数与实验值相比偏小,几个算例之间差别不大。即展向网格尺度越大,阻力计算结果越大,当计算域为二维情况时,阻力系数偏大,与实验值误差达到36.6%。而本文展向网格尺度最小的算例7,其阻力系数偏小,与实验值误差有19%。
算例4和算例7的展向长度不同,展向网格尺度相同,计算结果算例7稍好。同样,算例5和算例8的展向长度不同,展向网格尺度相同,其阻力系数结果与实验值符合最好。
综上,对于此类固定断面气动力平均结果,大涡模拟方法的二维计算的结果与实验值相比偏大较多,三维的计算结果在展向长度大于一个桥梁断面特征长度的情况下,保证一定的网格尺度,其计算结果与实验值相比能在20%以内,在工程上具有一定的参考价值。
5 结论
1)采用大涡模拟对闭合箱梁断面三分力系数的识别具有一定的工程意义。
2)与壁面垂直方向的近壁面网格分辨率至关重要,理论上完全分辨率近壁面湍流,其网格分辨率需要占到70%以上,合理的近壁处理仍是大涡模拟亟待解决的问题。
3)二维桥梁断面的绕流计算,展向计算尺度及展向网格尺度对预测结果有一定的影响。限于条件,对更长展向尺度在计算上的效果有待进一步研究。二维大涡计算对于快速预测所需结果,以及检查网格划分等条件的设定是否合理具有一定的参考意义。
大涡模拟的优势在于保留了过滤尺度以上的脉动信息,可以用于桥梁断面绕流的气动力谱分析和尾流特征的分析。同时,实际来流均有一定的湍流度,可以在来流中生成一定的目标湍流风谱以更符合实际情况。采用大涡数值方法对桥梁断面的气动力特性分析有待于进一步的研究。
摘要：以苏通长江大桥主桥闭口箱梁断面为例,采用大涡模拟方法的动态Smagorinsky模式,计算了雷诺数与风洞试验相当的断面气动三分力系数,并与实验值及雷诺时均方法的结果进行了对比。同时也进行了二维计算域的计算,并对不同三维计算域尺度及网格分辨率进行了计算对比,计算得到了不同计算域及网格对计算结果的影响情况。
关键词：闭口箱梁,大涡模拟,气动力系数,计算域
参考文献
[1]Davenport A G.Past,present and future of wind engineering.Journal of Wind Engineering and Industrial Aerodynamics,2002;90:1371—1380
[2]项海帆.现代桥梁抗风理论与实践.北京:人民交通出版社,2005
[3]刘天成.桥梁结构气动弹性数值计算的Lattice Boltzmann方法.同济大学博士学位论文,2007
[4]邓文.桥梁断面二维数值风洞的雷诺时均模拟方法.同济大学硕士学位论文,2006
[5]翟志轩.桥梁断面气动参数数值识别的基本问题.同济大学硕士学位论文,2006
[6]张兆顺,崔桂香,许春晓.湍流大涡数值模拟的理论和应用.北京:清华大学出版社,2008
[7]Sagaut P.Large eddy simulation for incompressible flows.Springer,2005
[8]Fluent6.3Documentation.Fluent Inc,2006
大涡数值模拟 第6篇
早期对开口火溢流的研究主要来自Yokoi,研究了不同开口时溢流的温度和速度分布及溢出火焰的高度。后来的学者对相关实验数据进行分析,提出预测溢流火焰高度的经验模型。近期国内外的研究仍然集中在溢流火焰高度、温度分布以及壁面热流密度分布等方面。
在实际的高层建筑中还存在一种特殊形式:开口所在的外墙可能与另一面毗邻单侧翼墙构成L型外立面、或者与毗邻双侧翼墙构成U型外立面。实际火灾案例表明,当外墙开口附近存在另一翼墙时,火焰溢流对外墙壁面的热危害可能更加严重。Oleszkiewicz在实验中也发现,毗邻翼墙能够很大程度改变溢出火焰的形态、增大火焰高度,并且能显著增加开口上方的热流密度。虽然近期国内外一些评价外墙保温材料防火性能的测试标准已经考虑了L型结构的影响,但其对溢出火焰行为特征的定量影响仍未见报道。对于U型结构,李建涛等人通过小尺寸试验发现该类结构能加速外墙可燃保温材料表面的竖向火蔓延,但该项工作中未涉对开口火焰溢流特征的分析。Lu等人也通过小尺寸试验研究了U型结构下溢流火焰的高度,并提出一个全局参量将此类工况下与平整壁面工况下的溢流火焰高度进行关联,该研究未对无翼墙、单侧翼墙和双侧翼墙条件下的火焰形态和外立面热流密度进行量化对比分析。
笔者在前人研究的基础上开展工作,重点分析不同开口高宽比、毗邻单侧翼墙(L型外立面)和双侧翼墙(U型外立面)条件下的开口溢出火焰形态、高度、外立面辐射热流密度变化规律,揭示外立面构造形式对开口火溢流的量化影响,为建筑外立面防火设计提供参考依据。
1 火灾场景与数值计算模型
1.1 火灾场景设计
研究使用美国标准与技术研究院开发的计算流体力学程序FDS 5.5开展数值模拟分析。FDS主要针对火灾驱动下的流体流动进行计算模拟,主要采用大涡模拟(LES)方法数值求解低速、热浮力动下流动的Navier-Stokes方程,其可靠性已通过大量试验验证。大涡模拟方法主要描述了气态燃料和燃烧产物与周围空气的湍流混合,其基本思想是混合产生的流体涡旋足够大,可对流体动力学方程给出足够精确的计算结果。
火灾场景设计主要考虑开口高宽比、建筑外立面毗邻单侧翼墙(对应L型外立面)和双侧翼墙(对应U型外立面)对溢出火焰特征的影响。燃烧室尺寸、翼墙和火源功率的设置参考了ISO 13785-2,燃烧室尺寸为4 m(长)×4m(宽)×2m(高),外立面高度为8m。图1给出了数值计算中使用的简化建筑模型。
火源类型为稳态正庚烷油池火,火源功率5.7 MW,俯视位置如图2所示。沿燃烧室外墙开口中心线竖向设置4个热流密度监测点,从距离燃烧室开口上方0.25m处开始布置,其他三个布置间距分别为0.75、1.00、2.00m,具体位置如图3 所示。燃烧室开口面积固定为2.4m2,通过改变其高度和宽度实现高宽比差异化,开口设置如图4所示。
1.2 数值计算模型设置
数值计算中选取合适的网格尺度对于计算结果的准确性至关重要。使用FDS时需要确保网格尺寸不大于0.1 D*,其中D*为表征火羽流的特征尺度,其关系式见式(1)。
式中:ρ∞为环境气体密度,kg/m3,取1.177 4kg/m3;Cp为空气比热容,取1 003J/(kg·K);T∞为环境温度,取300K;g为重力加速度,取9.81m/s2;Q为火源的热释放速率,kW。
场景的稳态火源功率为5.7 MW,因此采用更为精细的网格尺寸,将网格单元设置为长宽高均为0.1m的立方体,网格总数为1 600 000个。此外,数值模拟计算时考虑到计算区域边界的影响,对计算区域进行延展,燃烧室前方延展4m,其余各方向分别延展2m。模拟燃烧时间300s。
2 结果与讨论
开口火溢流属于自由扩散燃烧,火焰存在明显的震荡现象,火焰高度随时间不断变化。因此,参照火灾研究中较为通用的间歇率法确定火焰平均高度,即选择550℃作为火焰温度判据,火焰间歇率为0.5的位置所在高度即为火焰平均高度。
2.1 不同外立面形式的火焰宏观形态对比
首先分析不同外立面限制条件下的溢流火焰宏观形态,提取双侧翼墙(U型)、单侧翼墙(L型)和无翼墙下的火焰温度云图进行对比分析,分别考虑r为0.60、1.67两种典型情况。如图5、图6所示,双侧翼墙条件下的火焰高度明显高于单侧翼墙条件下的火焰,无翼墙条件下的火焰高度最低。这是由于翼墙限制了溢出火焰对周围空气的卷吸,从而使中性线温度升高,火焰高度变高。外立面构造形式从无翼墙向单侧翼墙和双侧翼墙的变化过程中,溢出火焰附壁程度逐渐增强,且在双侧翼墙条件下发生完全附壁燃烧。火焰附壁是建筑外立面抑制火焰卷吸,并在火焰四周产生非均一压力场所致,该现象有可能导致火灾通过建筑外立面加速蔓延。因此,对于具有双侧翼墙的U型外立面形式,火焰对空气卷吸的受限程度最高,导致火灾向建筑上层蔓延的风险更加严重。
此外,通过比较图5、图6还可以发现,开口高宽比对溢流火焰高度也有一定影响。与高宽比数值较小的扁平型开口相比,高宽比数值较大时火焰高度明显较高。对于U型外立面,这可能还与火焰对空气卷吸的限制程度有关,也即高宽比较大时,双侧翼墙的间距随之减小,因此对火焰卷吸具有更强的限制作用,从而导致该条件下的火焰更高。
2.2 单侧翼墙对溢出火焰高度和外立面辐射热流密度的影响
表1给出了在不同开口高宽比条件下,无翼墙和单侧翼墙(L型外立面)对溢出火焰高度影响的对比数据。从表1可以看出,单侧翼墙条件下的平均火焰高度比无翼墙限制条件下的火焰高出10%以上。这主要是因为毗邻翼墙可限制火溢流对周围空气的卷吸,随着进入溢流冷空气量的减少,其中心线上的温度衰减缓慢,因而具有更高的火焰高度。同时还可以看出,开口高宽比越大,单侧翼墙对溢出火焰高度的影响也更大。
表2给出无翼墙与单侧翼墙限制条件下开口上方壁面辐射热流密度数据对比。可以看出,单侧翼墙条件下开口上方的辐射热流密度明显高于无翼墙时的热流密度。r=1.67时,单侧翼墙条件下的壁面辐射热流密度高出48.79%~73.04%;r=0.60时,单侧翼墙条件下的壁面辐射热流密度高出52.06%~99.93%。对此,可结合火焰形态的分析给出解释,单侧翼墙在一定程度上限制了火焰对周围空气的卷吸,使其具有更高的火焰高度,同时火焰还向壁面发生倾斜。因此,壁面接收到的辐射热流密度更高。
2.3 双侧翼墙对溢出火焰高度和外立面辐射热流密度的影响
表3、表4分别给出了双侧翼墙(U型外立面)与无翼墙条件下溢出火焰高度与外立面辐射热流密度的对比数据。首先,双侧翼墙的溢出火焰高度显著增加,高出无翼墙的火焰高度60%以上,也远高于单侧翼墙产生的影响。
从表4的外立面辐射热流密度数据可以看出,r=1.67时,双侧翼墙条件下的壁面辐射热流密度比无翼墙时高出126.61%~1 625.60%;r=0.60时,双侧翼墙条件下的壁面辐射热流密度比无翼墙时的高124.67% ~498.45%。结合对火焰形态的分析可知,这是由于双侧翼墙进一步限制了溢出火焰对周围空气的卷吸,导致火焰只能从下面或者正面卷吸空气,导致火羽流竖向温度衰减更加缓慢,火焰高度也显著变高;此外,双侧翼墙条件下的溢流火焰均发生非常明显的附壁燃烧现象,火焰紧贴外立面向上蔓延,因此显著提高了外立面与火焰之间的辐射换热程度。辐射热流密度数据进一步说明,由双侧翼墙形成的U型外立面可加剧开口火溢流的发展,若不对其进行优化设计,则会对建筑外立面产生较为严重的火灾风险。
3 结论
(1)由于翼墙对火羽流卷吸空气具有限制作用,单侧翼墙条件下的平均火焰高度比无翼墙限制条件下的火焰高出10%以上,双侧翼墙限制条件下的溢出火焰高度显著增加,高出无翼墙条件下的火焰高度60%以上。
(2)外立面构造形式从无翼墙向单侧翼墙和双侧翼墙的变化过程中,火溢流的附壁燃烧程度逐渐增强,且在双侧翼墙条件下发生完全附壁燃烧。
(3)当r=1.67时,单侧翼墙条件下的壁面辐射热流密度比无翼墙时高48.79%~73.04%;双侧翼墙条件下的壁面辐射热流密度高出126.61%~1 625.60%。
(4)当r=0.60时,单侧翼墙条件下的外立面辐射热流密度比无翼墙时高出52.06%~99.93%;双侧翼墙外立面辐射热流密度高出124.67%~498.45%。
(5)双侧翼墙的开口火溢流具有更高的火灾危险性,增大了建筑外立面火灾竖向蔓延的可能性。
摘要：通过FDS数值模拟分析不同开口高宽比、L型外立面和U型外立面条件下的开口溢出火焰形态、火焰高度和外立面辐射热流密度变化规律。结果表明,单侧翼墙条件下的平均火焰高度比无翼墙限制条件下的火焰高出10%以上,双侧翼墙限制条件下的平均火焰高度比无翼墙条件下高出60%以上;溢出火焰的附壁燃烧程度在外立面构造形式从无翼墙、单侧翼墙、双侧翼墙的变化过程中逐渐增强,对外立面产生的辐射热流密度也随之明显升高,且火焰在双侧翼墙条件下完全发生附壁燃烧,对建筑外立面形成的火灾危害显著增大;开口高宽比越大,火焰高度越高。
大涡数值模拟 第7篇
关键词：混流式水轮机,全流道,CFD,三维非定常湍流计算,大涡模拟
0 引 言
水轮机稳定运行工况非稳定场分析计算主要是分析水轮机的水力稳定性,确定可能存在的引起水轮机运行不稳定的水力激励力。非稳定工况不稳定流场特性计算主要是研究水轮机开停机、调负荷等变工况过程中的流场变化特性。国内外电站水轮机运行的研究表明,水轮机开停机过程中所产生的水力激励力是造成叶片破坏和叶片裂纹的主要原因。因此分析开停机、变工况等过程的内部流动特性对于分析叶片破坏和叶片的裂纹具有重要的作用。目前可以采用国际上最新的CFD(计算流体动力学)技术对水电厂水轮机的开停机过程及变工况过程进行非稳定场分析[1]。
在对水轮机内部流动进行数值模拟计算时,全三维湍流N-S方程是目前描述流体流动的比较准确的数学物理方程。然而由于计算资源的限制,对水轮机内部流动这样复杂的湍流流动问题采用直接数值模拟(DNS)还是不现实的。目前在工程上求解全三维湍流N-S方程的方法主要有两类:①雷诺时均方程法RANS(Reynolds-averaged Navier-Stokes Equations);②大涡模拟LES(Large Eddy Simulation)。其中RANS方程是目前工程上主要采用的方法,同时随着计算机技术的不断发展,大涡模拟在工程上的应用日益广泛。
本文针对某水电站混流式原型水轮机HL100-WJ-75,应用大涡模拟动态亚格子湍流模型和滑移网格技术进行了全流道动静干扰的非定常湍流数值模拟,获得了一些重要的流动细节随时间的变化情况,计算结果更接近实际流场,为水轮机的水力设计和改型优化提供重要的参考依据。
本文原型水轮机的基本参数为:设计水头为170 m,额定转速为1 000 r/min,转轮直径为750 mm,叶片数17个,蜗壳进口直径650 mm,固定导叶数为8个,活动导叶数为12个,尾水管为弯锥型。
1 数值计算模型和方法
1.1 计算区域及网格
本文进行流动计算的区域为原型水轮机蜗壳进口到尾水管出口的整个流道,如图1所示。
针对计算域空间复杂的特点,计算网格采用以适应性强的非结构化四面体为主的混合网格。整个流道的网格节点总数为42.6万个,单元总数为202万个。其中蜗壳及导水机构单元数为53万个,转轮单元数为118万个,尾水管单元数为31万个。图2给出了混流式水轮机整体计算的网格模型。
1.2 控制方程和湍流模型
水轮机内部非定常流动的连续方程和N-S方程的雷诺平均形式为[2,3]:
式中:
亚格子雷诺应力τij定义式为[4]:
式中:δij为Kronecker常数;
动态模型系数C为[5]:
其中
上标^为采用试验过滤算子
式中:w(t)为指数权重函数。采用式(7)的平均方法可以避免动态模型系数C的剧烈波动和计算的不稳定,保持计算模型在空间的局部特性。这些特点有利于进行类似水轮机等具有复杂几何流道的流场计算。
1.3 计算方法和边界条件
采用贴体坐标下的有限体积法和非交错网格对控制方程在空间上进行离散,对控制方程中的源项和扩散项应用二阶中心格式,对控制方程中的对流项应用二阶迎风格式。在时间离散上,采用二阶全隐式格式。每个时间步上应用SIMPLEC方法进行离散方程的求解。
进口边界条件:在计算域的进口面上,给定速度值Uin和湍流强度。
出口边界条件:计算前出口速度和压力未知,采用自由出流边界条件。
壁面边界条件:在固壁处采用无滑移边界条件,近壁区采用标准壁面函数。
初始流场条件:在整个计算过程中,首先将转轮固定在某一个位置,按照标准k-ε模型进行三维定常湍流计算,然后将得到的定常流场结果作为大涡模拟的初始流场。
本计算中时间步长为0.001 s,当计算收敛后,时间步向前推进,同时转轮网格相应转动到新的位置,开始进行新时间步上的计算。
1.4 动静干扰(RSI)计算方法
动静部件耦合计算是通过在转动部件(rotor)和非转动部件(stator)之间采用特定的干涉面(interface)模式进行的,并采用滑移网格(Moving Mesh)模型。在混流式水轮机中,由于转轮的转动,分别与上游导叶和下游尾水管形成了两级动静干扰,采用滑移网格模型,在转轮进口前和转轮出口后形成两个网格滑移面,模拟动静干扰的非定常流场。
2 计算结果及分析
原型水轮机流场非定常CFD分析主要对偏工况下的小开度(导叶开度为25%)和大开度(导叶开度为100%)两种工况进行湍流计算。
2.1 引水部件流场随时间的变化
在不同工况下,对原型水轮机内的非定常流场进行了模拟,记录下水平截面上流场随时间的变化。图3为小开度工况不同时刻水轮机水平截面的压强分布图。可以看出,随着转轮的转动,引水部件中压强随时间变化而周期性地变化。
2.2 转轮流道内流场随时间的变化
从大开度工况转轮流道内不同时刻不同位置的横截面湍动黏度分布图(图4)可以看出由于在转轮的转动和导叶的动静干涉作用下,对水轮机转轮流道内水流的复杂湍流运动状态的影响。靠近转轮叶片进水边上部湍动黏度较大,且随时间的变化较为明显,说明转轮内局部湍流运动剧烈。图5为大开度工况不同时刻叶片工作面压强分布图;图6为大开度工况不同时刻叶片背面压强分布图。可以看出,叶片表面压强随时间的变化较为显著。
2.3 尾水管内流场随时间的变化
在不同的工况下,分别对原型水轮机的非定常流场进行模拟。以大开度工况不同时刻尾水管纵截面压强(图7)及湍动黏度分布图(图8)为例,可以看出,尾水管内压强随时间的变化较为显著。湍动黏度随时间的变化在尾水管进口处较为明显,可见转轮旋转对尾水管的动静干扰较强,使尾水管内的湍流运动较为复杂。
3 结 语
本文使用大涡模拟动态亚格子湍流模型和先进的滑移网格技术模拟水轮机中动静干扰的三维非定常流场,计算结果表明大涡模拟对复杂流场的仿真效果较为显著。随着计算机性能的提高,计算网格节点的增加,这一方法将更加准确。
需进一步开展的工作:①应用和发展其他高级湍流模型如分离涡模拟(DES)等,以更准确地计算水轮机中的非定常流动,达到对水轮机中的水流可以精细的模拟。②复杂湍流运动状态下涡旋产生的力学机理研究,这样才能对数值模拟出的流场作出正确合理的分析。
参考文献
[1]王正伟.流体机械基础[M].北京:清华大学出版社,2006.
[2]戴会超,槐文信,吴玉林,等.水利水电工程水流精细模拟理论与应用[M].北京:科学出版社,2006.
[3]王福军.计算流体动力学分析-CFD软件原理与应用[M].北京:清华大学出版社,2004.
[4]Smagorinsky J.General circulation experiments with the primitiveequations[J].Monthly Weather Review,1963,91(3):99-164.
[5]Germano M,Piomelli U,Moin P,et al.A dynamic subgrid-scale eddyviscosity model[J].Physics of Fluids A,1991,3(7):1 760-1 765.
大涡数值模拟范文
声明：除非特别标注，否则均为本站原创文章，转载时请以链接形式注明文章出处。如若本站内容侵犯了原著者的合法权益，可联系本站删除。
