直接数值模拟（精选9篇）
直接数值模拟 篇1
在油气田开发和生产过程中产生了大量的含油污水,如联合站产出水、洗井回水以及钻井废水等。与此同时,在油气田开发后期,油气井大量出水,地层能量衰竭。为了保持地层压力,需要大量注水。将油、气井产水经过处理回注地层,既可以保持地层压力,维持产能,又可以减少污水排放,形成产与注平衡的良性循环,具有十分重要的工业与环境保护意义。对含油污水进行处理后进行回注不但具有重要的社会意义,而且也具有巨大的经济意义。
由于过滤过程变化规律的复杂性,目前无论是过滤装置的设计还是运行,人们主要还是依靠实验得到的经验数据或经验准则关联式,因而利用污水过滤过程数学模型来研究污水过滤过程的变化规律具有特别重要的意义。
1 污水过滤过程的数学模型
污水的过滤过程的数学模型包括固体滤料相物料衡算方程、液体污水相物料衡算方程滤层过滤动力学方程和滤层压降计算方程四部分。1986年Homer等人[1,2]对均值滤床过滤进行了更加深入研究,提出精确反映滤床运行情况的数学模型表达式:
式中:
c—污水中悬浮物的体积分数;σ—滤层中悬浮物体积比沉积,cm3/cm3;
x—滤层深度,m;t—过滤时间,h;λ—过滤系数,1/cm;
v—污水过滤滤速,cm/s;εd—沉积物的空隙率;
ε'=ε0-σ为滤床过滤进行中的滤料实际孔隙率;
ε0—滤层滤料孔隙率。
要求解液固两相的物料衡算方程,还必需给出描述过滤系数值λ的过滤动力学方程。本文采用Ives提出的过滤动力学方程式[3],其形式如下:
式中:
σu—过滤饱和比沉积量,cm3/cm3;Bf—滤层中滤料的充填系数。
λ0—干净滤层过滤系数
滤层阻力计算方程采用文献[4]推荐的以下方程式:
式中:
H—水头损失,cm;μ—水的动力粘性系数,g/cm·s;
g—重力加速度,cm/s2;s0—清洁滤层内比表面积,cm2/cm3;
ρ—水的密度,g/cm3;ε0—清洁滤层孔隙率
要从污水过滤过程数学模型得出有意义的解,必需根据问题的实际背景给出其对应的初始条件及边界条件。对于固定的滤床过滤过程,其初始条件、边界条件如下:
式中:
λ0—干净滤层过滤系数,1/cm;c0—污水过滤前悬浮物体积分数;
2 特征线法求解数学模型
由方程的双曲特性可以采用特征线法来求解,为了采用特征线法对一般固定床污水的过滤过程的数学模型进行数值模拟,需要将方程在其特征线方向上离散。
方程组(1)、(2)的特征方向方程为:
特征线上的微分关系为:
将方程在特征线上离散得到离散方程如下:
同样地计算(7)、(8)中的系数,在点3的二次近似t3(2),x3(2);σ3(2),c3(2)可依次解下面方程组得到:
这个过程可继续,直到点3的坐标t3,x3和σ3,c3的前后两次近似相差为一个允许误差即可。
本文为固定滤床,取vs=0,这时借助Matlab编程语言,划分空间网格结点40个,语句约227条求的该问题解。
3 实例计算
为了验证上述过滤方程特征宪法离散数值模型,应用如下试验资料进行验证,污水动力粘度μ=9×10-3(g/cm.s),清洁滤层内表面积S0=75.1(cm2/cm3),污水密度ρ=1.0(g/cm3),重力加速度g=980(cm/s2),清洁滤层孔隙率ε0=0.34,滤层饱和比沉积量σu=0.12,干净滤层过滤系数λ0=0.03844(1/cm),滤层中滤料的充填系数Bf=5.7,滤层厚度l=100(cm),污水过滤前悬浮物体积分数c0=0.001,将这些数据代入特征线离散模型(8),(9),(10),(11)求解得出l=20(cm),l=40(cm),l=60(cm),l=80(cm)处的c/c0随过滤时间的穿透曲线图,见图1。
l=20(cm),l=40(cm),l=60(cm),l=80(cm)处薄层的体积比沉积量σ随过滤时间的曲线图,见图2。
由图1可知,在l=20,40,60,80(cm)处的c/c0随过滤时间的穿透曲线与试验所测基本吻合。
4 结论
直接分析均质滤料床过滤精确控制方程,合理利用方程的双曲特征,采用特征线计算法,消除了通常的向前差分格式引入的人工粘度造成的计算误差[5]。试验证明,计算值与实测值吻合良好,表明该计算法是可行的,计算结果是可靠的。并且通过上述分析,还可得到如下结论:
(1)特征线法是计算污水过滤过程数学模型的一种有效的迭代数值方法,对空间步长和时间步长要求不高,数值计算表明稳定性很好,其数值计算结果表明与实验数据非常吻合。当滤床移动,即滤床速度不为零时,从网格的几何形态与走向分析特征线法更具有优势。
(2)在本文有关过滤过程的数学模型中包含有许多经验系数,这些系数在不同情况下有不同的取值(如滤料种类、滤料粒径以及滤层厚度等),并且其中的一些系数需要依靠实验来得到,因此在应用模型对过滤过程进行模拟分析时要注意模型的适用范围。
(3)利用数学模型方法可以深化对污水过滤过程规律的认识,从而能够对污水的设计和运行起到一定的指导作用。
摘要：在应用固定滤床处理污水过程中,由于过滤过程变化规律的复杂性,人们主要还是依靠实验得到的经验数据或经验准则关联式。本文直接分析均质滤料床过滤精确控制方程,合理利用方程的双曲特征,采用特征线计算法进行了数值模拟分析,消除了通常的向前差分格式引入的人工粘度造成的计算误差。结果表明特征线法数值计算稳定,空间、时间步长选取范围较大,模拟结果合理,变化趋势与实验及有关文献理论结论定性一致。
关键词：含油污水,过滤,固定床,数值模拟,特征线法
参考文献
[1]张建锋,等.均质滤料滤层设计的优化研究.给水排水,2002,28,2.
[2]Ho rner R M W,Tarvis R J,M ack ie R I.Deepbed filt rat ion:a new look at the basic equat ions[J].W at Res,1986,20:215-220.
[3]Delo A K.Theo ry of sand filt rat ion[J].J sanit EngngD iv p rocAm SOC C IV Engr,1969,95:399-422.
[4]W righ t A M,KavanaughMC.D iscussion of theo ry of sand filt rat ion[J].J Engng D ivV p roc AmSOC C IV Engrs,1970,96:620-622.
[5]伯帜帜,等.均质滤料床直接过滤的数学模型及应用.给水排水,2000,26,9.
直接数值模拟 篇2
湍流射流的数值模拟
目的.用K-ε和K-W湍流模型及亚格子涡模型进行湍流射流的数值模拟.方法采用4步Runge-Kutta方法离散时间导数项,3阶ENO格式离散对流通量项,中心差分格式离散粘性通量项,数值求解Reynolds平均可压缩N-S方程.结果不同湍流模型都取得与实验比较一致的结果.结论用K-ε和K-W湍流模型封闭Reynolds平均可压缩N-S方程组能较正确地反映射流的湍流特性.
作 者：黄振宇 徐文灿 Huang Zhenyu Xu Wencan 作者单位：北京理工大学机电控制工程系,北京,100081刊 名：北京理工大学学报 ISTIC EI PKU英文刊名：JOURNAL OF BEIJING INSTITUTE OF TECHNOLOGY年，卷(期)：19(6)分类号：V411.3关键词：湍流模型 射流 基本无振荡格式(ENO)
直接数值模拟 篇3
在太阳能采暖通风系统中,空气在太阳辐射作用下于吸热壁上所形成的非稳态自然对流边界层流动是决定系统热性能的关键所在,因此,有必要透彻了解和掌握此类低普郎特数流体在恒定热通量作用下的非稳态自然对流边界层流动动力学行为[1～4].研究表明,所形成的自然对流边界层流动在达到稳态前将先历经起始和过渡两个阶段.而表征其瞬时流动特征的主要参数是壁温、热边界层厚度、内层和外层速度边界层厚度、边界层内竖直方向的最大速度以及边界层发展达到稳态所需的时间.最近,我们利用标度分析得到了均匀和线性分层低普郎特数(Pr<1)流体沿恒定热通量竖直平板的非稳态自然对流边界层流动在各个发展阶段的标度关系[1].本文则利用直接数值模拟来验证和定量化所得的线性分层低普郎特数流体沿恒定热通量竖直平板的非稳态自然对流边界层流动在各个发展阶段的标度关系.
1 标度关系
本文所考虑的是线性分层Pr<1流体沿恒定热通量竖直平板的非稳态自然对流边界层.设定施加于平板上的热通量为常量(即假定平板上的温度梯度为常数),流体的起始温度是线性分层的,即Ta(Y)=T0+SY,这里S为温度分层系数,它是个大于零的常量,而T0为Y=0处流体的起始温度.
假定流体遵从Boussinesq假设,且流动为二维的层流,则流体的非稳态瞬时流动可由Navier-Stokes方程和温度传递方程(文献[1]中的方程(2)～(5))来描述.依据文献[5],流动的特征长度Lc,特征速度Vc,特征时间tc和特征温度Tc可分别表示为
这里Pr=v/k为普郎特数,β，v和κ分别为流体的热膨胀系数、黏度和热导率、g为重力加速度、为无量纲温度分层系数.利用这些特征参数,可把表征自然对流边界层流动瞬时流动特征的主要参数(壁温Tw,热边界层厚度Δτ,内层速度边界层厚度Δvi,外层速度边界层厚度Δvo。以及边界层内竖直方向的最大速度Vm)在边界层的各个发展阶段的标度关系转化成如表1所示的无量纲形式(参见文献[1]的表2.需要指出的是由于文献[1]和本文所采用的特征参数不同,因而所得到的标度关系的无量纲形式也不相同).而边界层发展达到稳态所需的时间te的无量纲标度关系为,这也是边界层在过渡区内波动的周期tp的无量纲标度关系.
2 直接数值模拟方法
流动的控制方程可利用Lc,Vc,tc和Tc转换为如下的无量纲形式
这里x,y,u和v分别为水平方向和竖直方向的无量纲坐标和速度分量,τ,p和θ分别为无量纲时间、压力和温度.
直接数值模拟采用了有限体积法,具体详情参见文献[6,7].而本文所用的直接数值模拟程序的可靠性和精度已利用Sparrow和Gregg[8]的相关理论分析结果得到了检验.
3 直接数值模拟结果和讨论
为了验证和定量化上述各种标度关系,本文在0.01≤Pr≤0.5和0.2≤s≤5的范围内选择了有代表性的8个工况来进行直接数值模拟计算.其中工况1～4(对应于Pr=0.01,0.05,1和0.5,但维持s=1不变)用于验证标度关系对Pr的依赖性;工况5～8(对应于s=0.2,0.5,2和5,但维持Pr=0.1不变)用于验证标度关系对s的依赖性.
由于壁温,热边界层厚度和外层速度边界层厚度具有完全相同的标度关系,因而限于篇幅这里只给出壁温,内层速度边界层厚度和边界层内竖直方向的最大速度的直接数值模拟结果.
图1所示为壁温随时间变化的直接数值模拟结果.由图1(a)和1(b)可见,不同高度y处的结果在流动发展的各个阶段基本上都是相同的,表明壁温和相关的时间标度关系均与y无关,这是和标度分析结果相一致的.而图1(c)和1(d)则表明,在起始阶段初期,不同s值的结果均可用如下的定量化标度关系来描述
(a)和(b)为当Pr=0.1和s=1时不同y处的结果;(c)和(d)为当Pr=0.1时y=90处不同s值的结果;(e)和(f)为当s=1时y=90处不同Pr值的结果,其中粗实线代表定量化的标度关系式(5),粗点点线代表定量化的标度关系式(6),粗点划线代表定量化的标度关系式(7),粗划线代表标度关系θm,s=1.
这是和标度分析结果相吻合的;在起始阶段末端,不同s值的壁温是相同的,且发生在同一个时刻τw,e,其定量化的标度关系可表示成
这使得壁温θw,e与s之间的标度关系可定量化为
这也是和标度分析结果相一致的;在稳态阶段,所有不同s值的结果均趋近于θw,s=1,这正是标度分析所得出的结论.此外,结果还表明过渡区间内波动的周期τw,p与s之间的定量化标度关系可表示成
这同样是和标度分析结果相一致的.
图1(e)所示为不同Pr值在起始阶段的结果.由图可见,虽然在起始阶段初期不同Pr值的结果均重合于由式(5)表示的直线,表明θw～Pr-1/4的标度关系是成立的.可是,在起始阶段末端,不同Pr值的壁温并不一致,表明除了由标度分析预测的θw,e～Pr-1/4的标度关系之外,θw,e和Pr之间还存在更进一步的相关关系,虽然相对于θw,e～Pr-1/4的标度关系而言,这种更进一步的相关关系的尺度要小得多.类似地,图1(e)也表明不同Pr值的结果到达起始阶段末端的时间也有一定差异,这同样表明τw,e和Pr之间除了由标度分析预测的、由式(6)表示的标度关系τw,e～Pr-1/之外,还有更进一步的相关关系存在,虽然同样地这种更进一步的相关关系的尺度相比要小得多.然而,由图1(f)可见,不同Pr值的壁温在稳态阶段都最终趋近于θw,s=1,这是和标度分析的结果相一致的.而过渡区间内不同Pr值的波动周期τw,p也可由式(8)表示的定量化标度关系来表示,这也是和标度分析的结果相一致的.
为了充分体现在起始阶段末端θw,e和Pr之间的所有内在相关关系,可把由标度分析所得的、定量化的标度关系式(7)修改为
这里fθ(Pr)是Pr值的函数,其具体表达式可由图2(a)所示的直接数值模拟计算结果确定为
同样地,为了充分体现在起始阶段末端τw,e和Pr之间的所有内在相关关系,可把由标度分析所得的标度关系式(6)修改为
这些修正后的标度关系充分体现了Pr值对壁温的所有内在影响,这可从图2(b)所示的结果中得到验证.
图3所示为vm随时间变化的直接数值模拟结果.由图3(a)和3(b)可见,不用高度y处的结果在流动发展的各个阶段基本上都是相同的,表明vm和相关的时间标度关系均与y无关,这是和标度分析结果相一致的.而图3(c)和3(d)则表明,在起始阶段初期,不同s值的结果均可用如下的定量化标度关系来描述
这是和标度分析结果相吻合的;在起始阶段末端,不同s值的最大速度是相同的,表明vm与s值无关,其定量化的标度关系由直接数值模拟结果确定为
这是与标度分析结果相吻合的.结果也表明,不同s值的最大速度到达起始阶段末端的时间也是相同的,其定量化的标度关系可表示成
在稳态阶段,所有不同s值的结果均趋近于vm,s=0.322,这也是与标度分析的结果相一致的.此外,结果还表明过渡区间内波动的周期τm,p与s之间的定量化标度关系可表示成
这同样是和标度分析结果相一致的.
图3(e)所示为不同Pr值在起始阶段的结果.由图可见,在起始阶段初期不同Pr值的结果并不重合于由式(12)表示的直线,表明除了由标度分析预测的vm～Pr的标度关系之外,vm和Pr之间还存在着更进一步的相关关系.同样地,在起始阶段末端,不同Pr值的vm,e差异很大,也表明除了由标度分析预测的vm,e～Pr1/4的标度关系之外,vm,e和Pr之间还存在着更进一步的相关关系.但是,结果也表明不同Pr值的结果到达起始阶段末端的时间基本上是相同的,均可由式(14)表示,这与标度分析所得的标度关系Tm,e～Pr-1/2是相同的.此外,由图3(f)可见,不同Pr值的vm,s在稳态阶段都最终趋近于0.322,这也是和标度分析的结果相一致的.而过渡区间内不同Pr值的波动周期τm,p也可由式(15)来表示,这也是和标度分析的结果相一致的.
为了充分体现在起始阶段末端vm,e和Pr之间的所有内在相关关系,可把由标度分析所得的、定量化的标度关系式(13)修改为
(a)和(b)为当Pr=0.1和s=1时不同y处的结果;(c)和(d)为当Pr=0.1时y=90处不同s值的结果;(e)和(f)为当s=1时y=90处不同Pr值的结果,其中粗实线代表定量化的标度关系式(12),粗点点线代表定量化的标度关系式(14),粗点划线代表定量化的标度关系式(13),粗划线代表标度关系vm,s=0.322.
这里fv(Pr)是Pr值的函数,其具体表达式可由图4(a)所示的直接数值模拟计算结果确定为
类似地,可把标度关系式(12)中作如下的修正,以充分体现vm和Pr值之间的标度关系
这些修正后的标度关系充分体现了Pr值对最大速度的所有内在影响,这可从图4(b)所示的结果得到验证.
类似地,直接数值模拟结果表明(限于篇幅,这里未给出相关的图)内层速度边界层厚度δvi在起始阶段初期不同s值的结果均可用如下的定量化标度关系来描述
这是和标度分析结果相吻合的;在起始阶段末端,不同s值的δvi是相同的,表明δvi与s值无关,其定量化的标度关系由直接数值模拟结果确定为
这是与标度分析结果相吻合的.结果也表明,不同s值的δvi到达起始阶段末端的时间也是相同的,其定量化的标度关系可表示成
在稳态阶段,所有不同s值的结果均趋近于δvi,s=0.785,这也是与标度分析的结果相一致的.此外,结果还表明过渡区间内波动的周期τvi,p与s之间的定量化标度关系可表示成
这同样是和标度分析结果相一致的.
然而,在起始阶段初期不同Pr值的结果并不重合于由式(19)表示的直线,表明除了由标度分析预测的δvi～Pr1/2的标度关系之外,δvi和Pr之间还存在着更进一步的相关关系.同样地,在起始阶段末端,不同Pr值的δvi,e差异很大,也表明除了由标度分析预测的δvi,e～Pr1/4的标度关系之外,δvi,e和Pr之间还存在着更进一步的相关关系.但是,结果也表明不同Pr值的结果到达起始阶段末端的时间基本上是相同的,均可由式(21)表示,这与标度分析所得的标度关系τvi,e～Pr-1/2是相同的.此外,不同Pr值的δvi在稳态阶段都最终趋近于δvi,s=0.785,这也是和标度分析的结果相一致的.而过渡区间内不同Pr值的波动周期τm,p也可由式(22)表示的定量化标度关系来表示,这也是和标度分析的结果相一致的.
为了充分体现在起始阶段末端δvi,e和Pr之间的所有内在相关关系,可把由标度分析所得的、定量化的标度关系式(20)修改为
这里fvi(Pr)是Pr值的函数,其具体表达式由直接数值模拟结果确定为
类似地,可把标度关系式(19)中作如下的修正,以充分体现δvi和Pr值之间的标度关系
4 结论
在太阳能采暖通风系统中,空气在太阳辐射作用下于吸热壁上所形成的非稳态自然对流边界层流动是决定系统热性能的关键所在.自然对流边界层流动在达到稳态前将先历经起始和过渡两个阶段.而表征其瞬时流动特征的主要参数是壁温、热边界层厚度、内层和外层速度边界层厚度、边界层内竖直方向的最大速度以及边界层发展达到稳态所需的时间.
标度分析表明,对于线性分层的Pr<1流体而言,起始阶段边界层的发展是与高度无关而只依赖于时间的.而稳态时各项标度关系与处于起始阶段终结时刻的大不相同,且都与高度和时间无关,而只与Pr数和线性分层度相关.
直接数值模拟结果表明由标度分析所得的各项标度关系比较准确地揭示了流动的特性参数与流动的控制参数之间的内在联系.但并没有充分体现出对Pr数的所有依赖关系.这一不足可利用直接数值模拟结果得到弥补.
摘要：在太阳能采暖通风系统中，空气在太阳辐射作用下于吸热壁上所形成的非稳态自然对流边界层流动是决定系统热性能的关键所在．对于线性分层的低普郎特数(Pr＜1)流体而言，标度分析表明在起始阶段边界层的发展是与高度无关而只依赖于时间的，而处于稳态时，各项标度关系与处于起始阶段终结时刻的大不相同，且都与高度和时间无关，而只与Pr数和线性分层度相关．直接数值模拟计算结果表明由标度分析所得的各项标度关系揭示了流动的特性参数与控制参数之间的决定性内在联系．但同时也表明标度关系并没有充分体现出对Pr数的所有依赖关系，这一不足可利用直接数值模拟计算结果得到弥补．
关键词：非稳态自然对流,直接数值模拟,标度分析,低普朗特数流体,边界层,标度关系
参考文献
[1]刘滔，林文贤．低Pr数流体非稳态自然对流流动的标度分析．力学与实践，2007，29(1)：23～28(Liu Tao,Lin Wenxian. Scaling analysis of unsteady natural convection flow of a fluid with low Prandt(?) number.Mechanics in Engineering, 2007,29(1):23～28 (in Chinese))
[2]李光正．斜压流体产生涡旋的数值模拟．力学与实践，1997，19(6):39～42 (Li Guangzheng.Numerical simulation of swirl produced by baroclinic fluid.Mechanics in Engineer- ing,1997,19(6):39～42 (in Chinese))
[3]王明皓，符松，章光华．竖直平板间自然对流的湍流频谱特性．力学与实践，2004，26(1)：21～23(Wang Minghao,Fu Song, Zhang Guanghua.Spectrum characteristics of turbulent natural convection between two walls.Mechanics in Engi- neering,2004,26(1):21～23 (in Chinese))
[4]马丽娟，徐丰，胡非等．侧加热腔体内重力波演变过程的数值模拟．力学与实践，2006，28(5)：19～24(Ma Lijuan,xu Feng,Hu Fei,et al.Numerical simulation of the forma- tion and development of internal gravity wave in a differ- entially heated cavity.Mechanics in Engineering,2006, 28(5):19～24 (in Chinese))
[5] Prandtl L.Essentials of Fluid Dynamics.London:Blackie, 1952
[6] Lin W.Scaling Analysis and Direction Simulation of Un- steady Weak Fountain and Natural Convection Flow.[PhD Thesis].The University of Sydney,2000
[7] Lin W,Armfield SW,Patterson JC.Cooling of a Pr＜1 fluid in a rectangular container.J Fluid Mech,2007,574: 85～108
喷管分离流动的数值模拟 篇4
针对高空喷管地面试验时产生流动分离的现象,用Beam-Warming近似因式分解法求解薄层N-S方程,对分离流场进行数值模拟.预示了分离点的`位置以及喷管壁面压强分布等,可以为高空发动机地面试车提供参考,并为进一步研究分离的影响因素打下基础.
作 者：孙得川 李江 蔡体敏 王珏 王维彬 田昌义 Sun Dechuan Li Jiang Cai Timin Wang Jue Wang Weibin Tian Changyi 作者单位：孙得川,李江,蔡体敏,Sun Dechuan,Li Jiang,Cai Timin(西北工业大学航天工程学院,西安,710072)
王珏,王维彬,田昌义,Wang Jue,Wang Weibin,Tian Changyi(北京丰源机械研究所,北京,100076)
直接数值模拟 篇5
关于白云鄂博矿含碳球团直接还原过程的研究主要采用实验研究方法,采用与数值模拟相结合的研究较少,实验研究并不能获得实验过程中炉内球团的温度分布,无法准确分析结果。本文建立白云鄂博矿含碳球团直接还原过程的数学模型,用C++编程求解数学模型,利用实验数据来验证数学模型的正确性。通过正交试验分析炉气温度、反应时间、碳氧比、球团直径对金属化率的影响规律,从而为白云鄂博矿含碳球团直接还原工艺的优化研究提供理论支持。
1 模型的建立
1.1 物理模型
建立含碳球团的一维物理模型如图1,对模型作出合理的假设:(1)将球团看作多孔介质,且球团为规则的球形;(2)气体在球团内部只有扩散传质且热量通过导热传递;(3)忽略固—固反应,只考虑气-固直接还原反应(C与CO2,CO与Fe O、Fe2O3、Fe3O4);(4)球团内部同一位置,固体与气体温度相等;(5)不考虑球团二次氧化;(6)在反应过程中球团的大小不变;(7)在初始阶段气固分布均匀。
根据以上物理模型,建立白云鄂博矿含碳球团直接还原过程的数学模型。数学模型中导热系数与扩散系数及化学反应速率的计算方法参考文献[15,16,17],控制方程及边界条件如下。
1.1.1 能量方程
式(1)中:ρm为球团综合密度,kg·m-3;Cp,m球团综合热容,J·kg-1·K-1;T为温度,K;t为时间,s;r为半径,m;λeff为有效导热系数,W·m-1·K-1;Ri为化学反应i的反应速率,mol·m-3·s-1;ΔHi为i反应的化学反应热J·mol-1。
1.1.2 气相组分质量守恒方程
式(2)中:ε为球团孔隙率;ρg,k为气相组分k的质量密度,kg·m-3;Deff,k为气体组分k的有效扩散系数m2·s-1;Mk为物质的相对分子质量,kg·mol-1;Vi,k为反应i中K物质的化学计量数。
1.1.3 固相组分质量守恒方程
式(3)中:ρs,k为固相组分K的质量密度,kg·m-3。
初始条件:当t=0,0≤r≤r0时:ρs,k=ρs,k,0,ρg,k=ρg,k,0,T=T0。
边界条件:当r=0时,
式中:e为球体表面黑度;σ为波尔兹漫常数,5.678×10-8W·m-2·K-4;Tf为炉温,K;h为对流传质系数m·s-1;ρg,k,∞为主流气体中气相组分k的质量密度,kg·m-3。
2 实验
2.1 实验原料
实验所用原料是白云鄂博矿铁精粉,主要化学成分见如表1。含碳球团使用铁精粉、皂土和半焦与水混合均匀,用圆盘造球机造球。球团平均粒径约为12 mm,平均质量约为2.1 g,将造好的球团放在鼓风干燥箱105℃下干燥48 h后取出。
2.2 实验装置及方法
图2是实验装置示意图,实验所采用的加热炉为管式电阻炉(SKL16-BYL.包头云捷电炉厂),用控温系统用来控制炉温。通过电子天平连续称重球团质量的变化,同时将质量变化的数据传输给电脑,记录球团质量随时间的变化。在实验前,根据含碳球团的配比可知还原前球团的含碳量m1、球团中与铁结合的氧元素总质量Oi和铁精粉中总铁和总氧之比K。在实验过程中,当炉温上升至所需温度时,以4 L/min的流量向炉管内通入氩气保护;随后将坩埚放入吊篮内,同时用电脑采集含碳球团质量随时间变化的数据。实验结束后,将反应后的球团通过化学分析化验得出球团的残炭量m2,根据还原前含碳量m1与残炭量计算焦粉失重量Δm,如公式(4)所示。根据球团失重量Δw和焦粉的失重量Δm计算球团失氧量ΔO,之后根据公式(5)、式(6)计算还原度和金属化率[18]。
式中:R为球团的还原度;M为球团金属化率。
2.3 实验结果
图3是碳氧比为0.9的含碳球团在不同炉气温度下失重量随时间的变化规律。当温度为900℃时,在0~50 min之间质量从56 g缓慢减小到约为54.7 g;当温度为1 000℃时,在0~20 min之间质量从56 g较为缓慢减小到约为54.25 g,之后质量变化很小;当温度为1 100℃时,在0~10 min之间质量从56 g迅速减小到约为54.2 g,之后失重量变化很小可以认为球团质量没有变化,球内还原反应结束。分析不同炉气温度下球团质量随时间的变化可知,随着炉气温度的升高,球团质量减小速率加快,还原反应结束的时间缩短,这是因为炉气温度越高使球团温度升高,从而导致还原化学反应速率增加。
1为氩气瓶;2为电脑数据处理系统;3为天平支架;4为炉衬;5为电子天平;6为挂钩;7为炉管;8为吊篮;9-10为控温系统;11为进气管路
3 结果与讨论
3.1 数学模型的验证
采用三种实验工况来验证数学模型的正确性,分别在炉气温度为900℃、1 000℃、1 100℃研究了碳氧比为0.9的白云鄂博矿含碳球团金属化率随时间的变化规律。图4为碳氧比为0.9时白云鄂博矿含碳球团在炉气温度为900℃、1 000℃、1 100℃下金属化率的模拟值与实验值的比较,可以看出模拟值与实验值吻合较好,因此我们所建立的数学模型可以模拟白云鄂博矿含碳球团的直接还原过程。
3.2 白云鄂博矿含碳球团直接还原工艺参数优化
综上所述我们可以看出影响金属化率的主要因素有:炉气温度、反应时间、碳氧比、球团直径。我们通过正交试验的方法来模拟球团的直接还原过程,并通过极差分析寻找参数的最佳组合。每个因素考虑三个水平具体参数见表2,正交试验方案及结果见表3。
表4是极差分析结果,P是该水平下金属化率的平均值,平均值越高说明在该水平下金属化率越高还原越好。R是极差(平均最大值减去最小值),极差越大说明该水平对金属化率影响越大。从表中可以看出RA>RB>RC>RD,因此可以得出影响白云鄂博矿含碳球团金属化率因素大小排序:炉气温度>反应时间>碳氧比>球团直径,分析可得温度是影响化学反应速率的主要因素,而铁的生成量是由反应速率和时间共同决定的,因此二者是影响金属化率的主要因素。但是温度不仅影响反应速率而且间接影响铁的生成,因此温度的影响大于时间,而当时间和温度都达到足够大时碳的含量是决定金属化率大小的主要因素,可以看出我们所得出的结果是符合实际的。而最佳组合为A3B3C1D1即:炉气温度1 100℃、反应时间30 min、球团粒径12 mm、碳氧比0.9。
3.3 炉中球团内部温度的分布
图5是粒径为12 mm、碳氧比0.9的含碳球团在不同炉温制度下球表和球心处的温度随时间变化规律图。从图5中可以看出,球心和球表温度在炉温为200℃和600℃时温差不大可以认为是均匀的;当炉气温度变为1 100℃时,球表与球心出现温差且温差逐渐增大,在大约70 min时球内温度趋于均匀。对比不同温度球心与球表的温度可知,在低炉气温度下球团温度分布较为均匀,而当炉气温度为1 100℃球心与球表出现温差最大约为60℃。这是因为在低温时球团内部反应较为缓慢,吸热也较为缓慢对球团内部温度影响较小。当炉气温度达到1 100℃时球团内部各反应开始剧烈进行并且吸热。因此可以看出较高的炉温才可以使化学反应剧烈进行,而化学反应对球团内温度的分布影响很大。
3.4 不同炉气温度对金属化率的影响
根据3.2的结果,模拟研究两种工况的转底(板推式)还原炉对粒径为12 mm、碳氧比0.9的含碳球团金属化率的影响,为转底(板推式)还原炉实验提供数据参考。将加热段分为三段,工况一的炉气温度从第一段到第三段分别为:200℃、600℃、1 100℃;工况二的炉气温度从第一段到第三段分别为:500℃、1 100℃、500℃。两种工况冷却段温度为200℃,且四段的时间均为30 min。
图6是粒径为12 mm、碳氧比0.9的含碳球团在两种不同工况下金属化率随时间的变化规律。从图中可知,工况一在前60 min内金属化率变化很小可以认为是0,在60~69 min之间金属化率从0迅速提高到约为93%之后金属化率变化不大,可以认为球团内部还原反应已经结束;工况二在前30 min内金属化率变化很小可以认为是0,在30~40 min之间金属化率迅速从0增加到约为93%之后金属化率变化不大,可以认为球团内部还原反应已经结束。对比不同工况下金属化率随时间的变化规律可知,低温对球团金属化率作用不明显,高温段1 100℃对金属化率的影响很明显。这是因为在低温时三氧化二铁生成铁的量很少,而氧化铁还原为铁需要很高的能量,因此炉气温度较低时球团内部可以认为没有铁生成金属化率为0。而高温段使球团内部温度迅速提升导致各还原反应开始剧烈进行,金属化率也迅速提高。粒径为12 mm、碳氧比0.9的含碳球团在炉气温度为1 100℃时,各段反应时间可以控制在10 min左右。可以根据不同的炉气温度、碳氧比、球团粒径模拟白云鄂博矿含碳球团直接还原过程,通过模拟的结果来选择合适的反应时间,为实际生产提供理论支持。
4 结论
本文建立了白云鄂博矿含碳球团还原过程的一维非稳态数学模型,通过对比不同炉气温度下金属化率的模拟值与实验值,数学模型的准确性得到了验证。研究了含碳球团直接还原过程中球团内温度分布,通过正交试验研究影响金属化率的主要因素,以及不同炉气温度对含碳球团金属化率的影响,得到的主要结论如下。
(1)当炉气温度从900~1 100℃、碳氧比从0.9~1.1、反应时间从10~30 min、粒径从12~24mm时,影响白云鄂博矿含碳球团金属化率的主要因素排序为:炉气温度>反应时间>碳氧比>球团直径。金属化率最佳组合为:炉气温度1 100℃、反应时间30 min、球团粒径12 mm、碳氧比0.9。
(2)炉气温度为200℃、600℃时对含碳球团内部温度影响很小,球团温度趋于均匀;当温度等于1 100℃时球团球心与球表出现最高60℃的温差。因此可以看出较高的炉温才可以使化学反应剧烈进行,而化学反应对球团内温度的分布影响很大。
坡地建筑火灾数值模拟 篇6
关键词：坡地建筑,火灾,数值模拟,坡度,风速,宽高比
坡地建筑发生火灾的特点及规律不同于普通的平原上的建筑, 往往受到山地特有的气候、地形等因素的影响。例如:受山谷风 (如图1所示) 的影响, 白天和夜晚的风向不同, 对火灾的影响也会不同。使用FDS, 设置了15种研究工况, 对山坡上的建筑火灾进行数值模拟计算与分析。
1 模拟场景的设置
为研究各影响因素对山坡上的建筑火灾影响程度, 设置了一个标准场景作为研究的基准, 其他工况均在此标准场景的基础上稍作变化设定而成。为较为真实地模拟现实场景, 标准场景的建立按照实际建筑物、山坡等比例创建。考虑到数值计算结果受边界条件的影响, 研究对象四周留下充足的缓冲区域和边界。标准场景尺寸和建筑物关系按照图2的比例和尺寸进行设置。为了便于分析和讨论, 场景中设置一个山坡面, 一个坡底平面, 一个坡顶平面, 12栋建筑 (分别用大写英文字母标注为A至L) 。
1.1 标准场景设置
标准场景设置为长192 m、宽160 m、高72 m的区域。建筑物间距10m, 建筑物的尺寸为10m×10m×10m。建筑物距坡底的坡底缓冲区长30m, 建筑物距坡顶的坡顶缓冲区长30m, 最外侧建筑物缓冲区域55m。建筑物表面材质使用FDS属性库中的黄松木 (Yellow Pine) 。坡度选用中坡地25%的坡度, 整个山坡水平长度130m, 坡面长度134m, 垂直高度32m, 坡前缓冲区域长30m, 坡底缓冲区域长30m, 坡顶缓冲区长30m, 坡顶平面的缓冲区长32m。
为使场景设置更接近真实环境, 需要对风的属性进行相关设置。FDS中可以使用ATMOSPHERIC参数对风的属性进行设置, 此时FDS能模拟产生典型的大气效果并严格按照风廓线公式产生风, 如式 (1) 所示。
将z0 (大气廓线起源高度) 设置为10m, 同时将p (大气廓线指数) 设置为0.25即1/4分布。选用3m/s沿X轴正向吹入的谷风, 风将依次经过坡前、坡底、山坡面、坡顶、坡顶平台后离开。计算时间为1 800s。为使模拟场景中风场较为稳定, 标准场景设置了一个点火控制逻辑, 在模拟进行到300s时使着火点出现进行燃烧。火源设置在建筑物顶部, 火源尺寸为10m×10m, 火源功率为1 000kW/m2, 为非稳态增长t2火, 设置60s后达到峰值。
1.2 标准场景网格设置
由于计算主要产生在场景中央, 为提高计算效率, 标准场景网格按照内精外粗的划分标准, 对整个模拟区域的网格进行划分。计算区域是最为关注的部位, 因此将精度设置为1m;为了满足建模的需要和平滑坡面阶梯之间的缝隙, 山坡面区域精度设置为2m;空中区域相对而言对精度要求略低, 精度设置为4 m。网格总数量为321 584个。
1.3 标准场景探点设置
标准场景中主要关注建筑表面温度、建筑表面热辐射两个物理量, 因此设置了温度探测器、热流量探测器两种探点。由于场景设置的对称性, 只在建筑一侧位置设置相关探点:建筑表面使用温度探测器, 命名为T;建筑表面用热流量探测器, 命名为RHF;建筑表面探点以建筑物为参照物设置, 按照图3所示的前、右、后、左、顶的顺序依次编号为1~5, 按照建筑物名加探点名加编号的方式命名。例如, 建筑物A顶部的感温探点以A-T05命名, 建筑物F前面的热流量探点以F-RHF01命名。
2 工况的设置
2.1 坡度工况设置
坡度选用我国坡地类型中占较大比例、较适宜建设建筑物的平地、缓坡地、中坡地和陡坡地, 设置4种工况进行计算与分析, 详见表1。
2.2 风工况设置
为研究风对山坡上的建筑火灾影响程度, 选用山坡常见的两种风向:山风和谷风;另外选用三种常见的风速:3、6、9m/s。设置6种工况进行计算与分析, 详见表2 (风设置采用了1/4法则, 参照式 (6) ) 。
2.3 宽高比工况设置
为研究宽高比对山坡上的建筑火灾影响程度, 设建筑高度不变, 选用宽高比1∶1、1∶2和2∶1三种工况进行计算与分析, 详见表3。
2.4 建筑布局工况设置
为研究建筑布局对山坡上的建筑火灾影响程度, 按照建筑对齐分布和交错分布设置两种工况进行计算与分析, 如图4所示。风向为谷风, FDS中风向为X轴正向, 风速为3m/s;工况一为交错分布, 工况二为对齐分布。
3 数值模拟与分析
通过FDS的SmokeView对计算结果进行可视化观察与分析, 同时对各探点数据绘制曲线图, 主要对结果中的温度、速度、热释放速率、建筑表面温度、建筑表面热流量等物理量及对应的矢量图进行对比分析以寻找规律与特点。
3.1 标准工况温度分析
建筑G、F、B、C受火灾影响程度最大, 其中平均温度最高的是G建筑表面的322℃, 而相对最低的C建筑表面平均温度也达到81℃。受火灾影响程度最高的建筑表面依次为建筑G的前部、建筑F的后部、建筑B和C的左侧。对受火灾影响最大的G建筑的前部进行重点分析。
3.2 坡度工况分析
选取坡度工况中Y=80m处的温度切片对4种工况进行分析, 如图5所示。从图5中的等温线可以看出, 随着坡度的增大, 着火建筑F对其下风方向建筑G的危害程度逐渐增大。燃烧产生的火焰和热烟气距离建筑G的前部和顶部越来越近。从建筑G前部的表面升温曲线 (如图6所示) 也可以看出, 随着坡度的增大, 建筑表面温度越来越高, 升温速度越来越快。其中, 25%的坡度和40%的坡度工况中升温曲线较为接近, 40%的坡度工况中温度略高于25%的坡度工况, 其他工况中升温曲线层次较为分明。
3.3 风工况分析
图7为720s时位于模拟场景Y=80m处的温度切片, 可观察到随着风速的增大, 着火建筑F对其坡上建筑G的危害程度逐渐增大。图7中阴影面为焰心高温区域, 黑色虚线为指定的中间对比温度, 黑色实线为指定温度25.5℃的等值面。其中, 工况一的中间对比温度为32℃、工况二的中间对比温度为30℃、工况三的中间对比温度为33℃。可以看出, 等值面随着风速的增大而增大, 燃烧产生的火焰和热烟气距离建筑G表面越来越近, 温度越来越高, 覆盖的范围也越来越大, 对建筑G的危害程度也逐渐增大。
下风方向的建筑在风的作用下受火灾影响的程度较大。不同风速下建筑G的前部表面升温曲线, 如图8所示。由图8可以看出, 建筑G受火灾影响最为显著;随着风速的增大, 建筑表面升温速度逐渐增大, 三组升温曲线层次较为分明。
环境风对避开主风向的建筑起到减缓升温的作用。同风速下建筑B的左侧表面升温曲线, 如图9所示。由图9可以看出, 随着风速的不断加大, 位于着火建筑F右侧的建筑B表面升温速度呈现出减缓趋势。一方面由于环境风带来的冷空气对建筑表面产生冷却效果, 风速越大冷却的越快;另一方面, 随着风速的不断增大, 气流将建筑燃烧产生的热量带走的也越来越多。
3.4 宽高比工况分析
图10为不同宽高比下热释放速率分布图, 可观察到随着宽高比的增大, 着火建筑F对下风方向建筑G和建筑H的危害程度逐渐减小。工况一时建筑G的前部和顶部被建筑F燃烧的火焰覆盖;工况二时建筑F与建筑G之间的安全距离加大, 建筑G受建筑F火灾的影响程度减小;工况三时建筑G与建筑F之间的距离进一步加大, 受火灾影响程度更小。
3.5 建筑布局工况分析
不同建筑布局下B建筑表面温度分布, 如图11所示。由图11可以看出, B建筑受建筑布局变化影响较小, 这是由于建筑布局的改变对B建筑的位置、高度等参数几乎没有影响, 只有环境风在两种布局下流经建筑表面略有不同。对齐分布时升温曲线呈小幅震荡上升, 而交错分布时则较为平滑。
不同建筑布局下G建筑表面温度分布, 如图12所示。与建筑B的升温曲线相比可以看出, 两种建筑布局的变化对其表面升温的影响非常显著, 对齐分布时升温较为迅速, 而当建筑布局改为交错分布时升温速度明显变缓。一方面交错分布时建筑G较对齐分布时位置向右移动了10m, 远离了着火建筑;另一方面, 由于错开的空间使环境风对建筑的冷却作用更加显著。
4 结论与建议
(1) 坡度与火灾危险性呈正相关的关系。若修建在较陡的山坡上则应提高建筑材料耐火等级并适当增大建筑间防火间距。
(2) 风对下风方位的建筑火灾危险性呈正相关的关系, 建议坡地建筑选址时应根据当地风向分布概率, 使建筑位置尽量避开主风向, 降低火灾发生时的危害程度。谷风风向时的火灾危险性大于山风风向时的火灾危险性, 建议救援人员在风向变换时注意加强个人防护和调整救援时机。
(3) 建筑宽高比与火灾危险性呈反相关的关系, 建议坡地建筑防火间距应在各消防设计规范中选用最不利条件设置并适当增大防火间距。
(4) 对齐分布的建筑布局火灾危险性高于交错分布的建筑布局, 建议坡地建筑群采用交错分布的建筑布局。
参考文献
尾矿坝溃坝数值模拟研究 篇7
近年来随着社会对矿产需求的增加, 尾矿的总量也在急剧的增长, 尾矿库数目也随之增加, 堆坝高度也越来越高, 同时尾矿坝溃坝的风险也越来越高[1], 又加上服务年数超限的尾矿库越来越多, 导致尾矿库溃坝事故的数目逐年递增。又由于多数尾矿堆存物具有一定的毒害性[2], 这更加加剧了尾矿溃坝带来的危害, 现在已成为政府重点关注对象。
尾矿坝溃坝作为一种灾害事故具有瞬间爆发性和很强的破坏性, 要想得到真实可靠的溃坝现场监测数据是非常困难的, 并且需要很大的人力和财力的投入。为了得到更多的数据, 更好地研究分析尾矿坝溃坝相关机理, 不得不另辟蹊径。经过大量的研究与尝试后, 认为使用Fluent3D软件来模拟泥砂流的流动过程有一定的可行性, 它能为研究提供运动和动力方面的数据, 填补现场数据的不足。
1 尾矿坝基本概况
1.1 尾矿颗粒组成
尾矿的颗粒组成不仅决定了它的物理性质, 也决定了其力学性质, 如渗透性、压缩性和剪切强度等。为了得到准确的尾矿颗粒组成情况, 试验人员随机选取了5袋尾矿样并将其混合均匀, 进行细化处置后, 取10组试样, 使用美国生产的Microtrac S3500型激光颗粒分析仪进行颗粒分析。从颗粒测试结果可知, 10组全尾矿样均为尾粉土, 全尾矿颗粒粒径的中值粒径在0.016~0.088 mm之间, 平均值为0.044 mm。
1.2 物理性质试验测试
按照试验规范要求, 取了3组尾矿样, 分别对它们进行了比重、密度、含水率、塑限和液限等物理性质参数的测定, 尾矿物理性质试验测试结果如表1所示。
注:密度是将制备好的尾矿样按照15%含水率配制, 然后根据《土工试验规程》测定。
2 尾矿坝溃坝模拟试验
根据库区地形, 试验模型按1∶500比尺缩小, 则模型库区尺寸为3 m×2.8 m×0.5 m (长×宽×高) , 参考以往尾矿坝溃坝泥石流冲击距离, 该次试验下游冲沟设计为15 m, 模拟溃坝泥石流在下游冲击距离为7.5 km的整个情况。根据设计模型尺寸及其他条件综合考虑, 确定相似系数如表2所示[3]。
式中, Lm为模型尺寸;Lp为原型尺寸;tm为模型过程时间;tp为原型过程时间;λρ=1。
2.1 溃坝后泥深变化规律
通过对下游各过流断面的泥浆淹没高程进行整理, 得到了各断面泥深随时间变化曲线, 如图1所示。
2.2 溃坝后冲击力变化规律
不同尾矿坝高度溃决后, 泥浆在库区下游5 m处的冲击力过程情况, 如图2所示。总体上看冲击力曲线呈现前陡后缓的趋势, 意味着溃后泥砂流的冲击力在很短的时间内便能达到其峰值, 然后会逐渐降低。
3 尾矿坝溃坝数值模拟
把尾矿库溃坝所形成的泥流浆体看作是空气和泥砂流的两相流体, 可用欧拉-欧拉模型来描述其运动。在欧拉-欧拉模型中, 常见的模型有三种, 分别为:流体体积模型 (VOF) 、混合物模型以及欧拉模型。经过验证, 选用流体体积模型 (即VOF模型) 对尾矿坝溃坝泥砂流的运动过程进行数值模拟比较合适[4,5]。
为了方便模型的建立, 将实际情况进行了相应的简化, 比如沟谷的边坡角度都统一设置成一个大致的平均值, 模型高100 m, 库区模型尺寸为1 500 m×1 400 m×100 m, 整个下游流动沟槽长度为7 500 m。参照前面得到的尾矿物理力学性质以及原始边界条件的初始化, 开展溃坝泥砂流流动的数值模拟。
3.1 溃坝后泥深变化规律
通过图3可以知道, 不同断面上的泥深的变化都呈现前面陡峭后面平缓的趋势。就某一断面来讲当溃泥砂流到达后, 在很短的时间内达到最大值, 随后就不断降低, 出现明显的拖尾现象。其中转弯处有点异常, 主要是因为转弯处由于弯道效应而使得泥深要比直道的情况下高出很多。综合来看, 溃坝后下游将会出现的最大淹没深度是95 m左右, 距离坝址下游5 km范围内最后会形成约20 m的泥砂沉积深度, 也就是说下游5 km内的村庄基本都会被淹没掉。
3.2 溃坝后冲击力变化规律
由图4中可以知道, 100 m高尾矿库溃决后对下游造成的冲击力可以高达2.6 MPa, 泥砂流一旦到达某处则迅速产生很高的冲击力, 在很短的时间内达到峰值, 随着库区泥砂流来源的减少和流动距离的增加冲击压力将迅速的减小。而对于下游0.6 km处的米茂村来说, 受到的冲击力可以说是最大的, 高达2.5 MPa左右, 所以该村的防护工作应重点把握。
4 数值模拟与相似模拟试验结果对比
为了验证数值计算结果的可靠性, 将其与相似模拟试验所得到的结果进行对比分析, 下面就整体趋势、泥深变化和压力变化的情况进行简单的比较。
4.1 整体趋势比对
根据相似模拟试验结果比对分析得出:无论是淹没深度还是冲击力变化情况, 总体上都呈现前面迅速增长、峰值以后都逐渐减小, 并趋于平稳值的现象。都有明显的峰值, 并且在相同条件下, 两种计算结果得出的规律是类似的。
4.2 淹没深度结果比较
当坝高100 m的情况下, 根据溃后库区下游2.5 km处的淹没深度的峰值情况, 包括相似模拟和数值模拟两个结果。两者之间的差值在25%左右, 处于可以接受的范围。
4.3 冲击力结果比较
同样是抽取坝高100 m, 库区下游2.5 km处的冲击力数据进行比对, 在此处冲击力的峰值结果相差有点大, 数值计算结果为1.9 MPa, 相似模拟试验结果为6 MPa。根据相关工程实际经验认为此结果可以接受, 处于合理误差范围[6]。
4.4 原因分析
相似模拟试验与数值模拟的结果不同主要有以下几点原因: (1) 试验设备器材与数值模拟中选材之间有一定的差异, 不尽相同; (2) 相似模拟试验中相似材料选取与数值模拟中材料有一定的差异, 不尽相同; (3) 相似模拟试验设备中溃坝瞬间控制系统无法像数值模拟那样做进行完美模拟; (4) 相似模拟试验设备中应力采集系统存在一定的误差, 而数值模拟中采用力学模型以及算法也存在一定的缺陷。这些都是导致两者结果有差异的原因, 但是结果在工程要求的范围内。
5 结论
利用Fluent3D软件对坝高100 m的尾矿坝溃决过程进行数值模拟, 可以得出以下结论:
(1) 用Fluent3D软件对尾矿坝溃决过程的模拟结果与采用相似模拟试验得到的结果基本相同, 溃坝后泥砂流的泥深分布与冲击力变化曲线都具有相同的趋势, 并且各自的峰值的相似度较大, 均在工程计算允许的误差范围内。
(2) 证实了用Fluent3D软件来模拟尾矿坝溃坝后泥砂流流动的可行性, 进一步拓展了尾矿坝溃坝研究的方式和手段, 能为以后溃坝后泥砂流的研究提供一个参考。
(3) 解决了因尾矿坝溃坝具有瞬间性及极强破坏性而不易监测到溃坝现场数据的难题, 为进一步研究尾矿坝溃坝防护措施打下了基础。
参考文献
[1]史冬梅.泥浆流变模型的判别方法[J].西部探矿工程, 2005 (9) :141-142
[2]周娜娜, 汤亚飞, 梁震.磷矿浮选尾矿水污染物释放的实验[J].武汉工程大学学报, 2011, 33 (3)
[3]敬小非.尾矿坝溃坝泥石流流动特性及灾害防护措施研究[D].重庆:重庆大学, 2011
[4]魏文礼, 郭永涛.基于加权本质无振荡格式的二维溃坝水流数值模拟[J].水利学报, 2007, 5 (5) :596-600
[5]ZHOU Di-hua, SHEN H W, LAI Jihn-sung, et a1.Approximate Riemann solvers in FVM for 2D hydraulicshock wave modeling[J].J Hydraul Eng, 1996, l22 (12) :692-702
水管路水力特性数值模拟 篇8
关键词：管路,水力特性,数值模拟
1概述
无论是湍流流动状态还是层流流动状态, 由于粘性力的作用, 在流动的时候, 或多或少都会出现机械能的损失。在水流以层流状态流动时, 工质水与水管壁之间, 工质的不同流速的流层之间, 在粘性力的作用下, 所产生的阻力都是对工质本身在做负功。工质的机械能大量地被损失[1]。所以研究不同规格管路引起机械能损失的大小就显得格外重要。
目前国内外学者对水管路水力特性已经做了相关研究, 赵宝峰等学者主要是考虑到了突然扩大管的局部水头损失系数计算中沿程损失的影响这一因素, 从而修正了研究并结合试验对公式[2];贺益英等人通过试验研究了弯管的局部水头损失系数并分析输水管线中弯管局部水头损失对其附近的水力参数的影响[3];Valiantzas通过考虑局部水头损失和沿程水头损失的关系对威斯一达西公式进行了改进[4];茅泽育等通过试验对管道汇流口局部阻力进行了研究[5];宋长福等对有连接管的调压室水头损失系数通过模型试验进行了研究[6];王开等通过模型试验验证考察了不同计算公式的性能对倒虹吸有压弯管的局部水头损失系数计算方法的影响对其进行比较研究[7]。
汽水管网是将锅炉、汽轮机、换热器等热力设备予以连接并传输工质的管道及其附件的总称, 是整个蒸汽动力系统的重要组成部分[8]。汽水管网系统的应用及其广泛, 在火力发电厂, 核电站, 供热系统中均有大量的汽水管网, 而水管路是其中必要的组成部分。其中直管和弯头所占比重较大, 因此对直管与弯头引起机械能损失的大小就显得格外重要。
2数学模型
水力计算基本方程:
2.1质量守恒定律 (连续性方程)
式中:———流体密度, kg/m3;
t———时间, s;
u, v, w——速度矢量沿x、y、z方向的分量, m/s。
2.2动量守恒定律
实质为牛顿第二定律。对于符合牛顿切应力公式的流体, 可得动量守恒方程为:
式中:p———流体微元体上的压力, Pa;
Fx, Fy, Fz——作用在工质微元体上的体积力, N。
2.3能量守恒定律
它是由热力学第一定律推导出来的, 能量守恒的矢量形式:
式中:cp———定压比热容, k J/ (kg·K) ;
T———热力学温度, K;
K———导热系数, W/ (m·K) ;
ST———粘性耗散项, 流体的内热源并且在粘性力作用下工质机械能损失的部分。
2.4控制方程的通用形式为:
流体流动中, 各方程可采取相同的表达形式, 使用通用变量, 通用微分方程为:
式中:———通用变量, 可以代表u, v, w, t等求解变量;
———广义扩散系数;4
T———广义源项。
3网格无关性验证验证
以DN25, 3m的直管为例, 具体观察网格数量对所模拟结果的影响。如图1所示。正如图1所示, ft表示用经验公式所计算出的摩擦因子的经验值, f表示用软件模拟出来的经验值, 从图中可以看出在网格数为150000左右时, 两条线近乎重合, 而从整体来讲网格数从100000到350000的过程中两条线的差距都是非常小的, 结合图和所计算工况可看出最大的差距为在网格数为320000时, 经验值为0.00629825, 而计算得出的值为0.0062485, 计算出来的相对于经验值的误差仅为0.78%。所以在对直管进行网格划分的时候, 网格数可以不用划分过多, 这样的话可以节省计算时间提高工作效率。从图中可知对直管的网格划分网格数控制在150000左右即可。这样会使效率最高。
下面以DN50, 直管段管长3m为例, 观察弯头的网格数对局部阻力系数的影响。如图2所示。从图2可以观察到, 弯头的所划分的网格数量对最终结果的影响还是比较大的, 在这一点上, 与直管的网格无关性的结论不一样。其实从两个物理模型的形状就可看出直管在网格划分一步的网格质量较弯头来讲是比较容易保证的。所以说, 对弯头的网格划分难度较大, 同样具有较高的不确定性。弯头的物理模型包括了两段直管短和弯头区域, 若在流速方面划分同样的份数, 弯头的网格质量较低。
网格数越多的话效果就会越好。从图中可见, 网格数越多局部阻力系数的计算值就会越来越接近经验值。所以在对弯头进行网格划分保证网格数量尽量多。而又要兼顾工作效率, 所以网格数选到300000左右即可。
4管件模拟计算分析
4.1直管水力特性分析
一共选取了七种大小的管径, 分别是DN15, DN20, DN25, DN32, DN40, DN50, DN65。在每一组管径下, 在对七个不同流速下的状况进行研究, 而为了保证网格质量, 将各管径不同的直管按不同长度设定。
如图3所示, 可以清楚地观察出管径和流速分解对摩擦因子的具体影响。最上面一条线为管径15的直管, 在下面一条是管径为20的, 以此类推, 最下面的一条线是管径为65的直管。所以, 从这张图上可以看出随着管径的增加, 摩擦因子是在逐渐变小的。这样的变化趋势也同样佐证了经验值的正确性。若只关注一条曲线的话, 每条曲线的大致趋势也是基本一致的。就是随着流速的增加, 摩擦因子的数值是在逐渐减小的。与此同时还可以看出, 各条曲线都是在流速较低的时候斜率大, 随着流速的增加斜率是逐渐变缓的。各条曲线大体的趋势都是如此, 足以说明这是一个共性的问题。上述特点就为各条曲线的共性。通过观察单条曲线的走势也可得出:当速度很大的时候, 流速基本不会对摩擦因子的大小产生影响。流速越小, 流速对摩擦因子的影响就越大。
4.2弯头水力特性分析
在研究管径大小对局部阻力系数的影响时, 仍可以采用控制变量的方法。但正如上文讨论所提, 利用控制两个变量的方法进行研究时, 图像过于简单, 给出的信息量较少, 又不能与其他变量之间进行横向的比较。所以在对最后一个变量因素进行研究的时候直接采用控制一个变量的方法。
首先, 对D/R进行控制变量, 探究流速与管径之间的关系, 如图4所示。
在图4中控制D/R=1, 令管径和流速变化从而观察规律。从图中看出, 最上方的一条线为DN20, 然后是DN32, 而最后一根线代表DN65。很明显的发现随着弯头管径的增加, 局部阻力系数是在不断地减小的。这是我们需要根据图上给出的信息将流速和管径两个因素作比较, 观察那个因素对局部阻力系数的影响更强烈。仍旧按照上述方法进行分析发现, 四条线之间并没有发生互相交叉重叠的现象。也就是说DN20的弯头流速无论怎么改变其所对应的局部阻力系数永远大于DN32。其他几条线的关系以此类推, 从而得出结论:弯头的管径越小, 其所对应的局部阻力系数就越大。管径作为影响局部阻力系数的一项因素, 其影响的程度远远大于流速对其的影响程度。在本章图上同样可以得出流速的规律, 仅仅观察一条曲线就会发现上述规律, 在这里不再赘余。
5结论
5.1在利用FLUENT软件对直管和弯头划分网格和数值模拟计算时, 通过各自的网格无关性验证可以得出结论:对直管所划分的网格数的多少基本不会对最终计算的摩擦因子的数值产生重大影响;然而在弯头方面, 所划分的网格数越多, 计算工况所得出的局部阻力系数就越接近相应局部阻力系数的经验值。
5.2对于直管的摩擦因子进行研究时, 可以得出结论:随着工质在直管内流速的增大, 摩擦因子的数值是在不断减小的, 这种现象在流速较低的时候比较明显。随着流速的渐渐增加, 流速的变化对摩擦因子的影响就越小, 直到流速足够大的时候, 流速的变化将不会对摩擦因子产生任何改变。同时, 随着直管管径的增加, 摩擦因子的数值是在不断减小的。
5.3对于弯头的局部阻力系数进行研究时, 可以得出结论:管内工质流速, 弯头的几何因素即D/R以及弯头的直径均会对局部阻力系数产生影响。弯头的直径越大, 计算所得的局部阻力系数就越小;弯头的D/R越大, 局部阻力系数就也随之增大;工质在弯头内流动的速度越大, 局部阻力系数就越小。
参考文献
[1]赵秀凤.引水遂洞洞内消能水工模型试验的研究[D].长春:华北水力水电学院, 2006.
[2]赵宝峰, 金英子, 卢玉邦等.对突然扩大局部水头损失的初探[J].东北农业大学学报, 1997, 6:175-178.
[3]贺益英、赵懿瑶、孙淑卿等.弯管局部阻力系数的试验研究[J].水力学报, 2003, 23:54-58.
[4]VALIANTZAS J D, Modified Hazen Williams and Darc-Weisbach.Equations for friction and local head losses along irrigation laterals[J].Journal of Irrigation and Drainage Engineering, 2005, 131 (4) :342-350.
[5]茅泽育、赵凯、赵璇等.管道汇流口局部阻力试验研究[J].水力学报, 2007, 27:812-818.
[6]宋长福、夏小娟.有连接管的调压室水头损失系数研究[J].水力科技与经济, 2005, 10:604-606.
[7]王开, 傅旭东, 王光谦.倒虹吸有压弯管的局部水头损失系数计算方法比较[J].水资源与水工程学报, 2006, 16:1-4.
导弹头锥烧蚀数值模拟 篇9
流场和热环境的研究是导弹头锥烧蚀和温度场计算的基础, 国内外许多研究人员求解了抛物化的N-S方程。在结构烧蚀方面, 比较普通地采用气动化学烧蚀模型。王中原先生通过应用小扰动理论研究了弹翼间的气动干扰, 采用轴对称比拟的方法考虑了有攻角 (非对称) 情况, 引入热化学反应, 建立了高温下的烧蚀准则[1]。陈林泉先生以热化学烧蚀二方程模型为基础, 用烧蚀与传热耦合计算的方法, 研究了梯度功能材料的烧蚀及热防护机理[2]。本文在“飞机结构三维温度场分析”程序[3]中增加了烧蚀功能程序, 进行了导弹头锥温度场/烧蚀耦合计算。
对于高速飞行的导弹, 工程上常用热壁热流来描述空气动力加热量, 热壁热流是壁焓的函数, 壁焓又是壁面温度和压力的函数, 热壁热流和温度很难建立起显式的代数关系, 因此空气动力加热和结构温度场分析是相互影响的, 在理论上是一个流固耦合问题, 不能将两者割裂开来讨论。头锥一般都采用碳-碳材料, 这种材料的物性随温度变化是一个瞬态导热过程, 加之烧蚀产生的热效应和物面边界的移动, 使得问题更为复杂。因此, 对导弹头锥应该进行烧蚀和瞬态温度场的耦合计算。
飞行器结构大量使用杆、梁、肋、蒙皮等组合结构, 各结构件之间多采用铆接、螺接、胶接等连接型式。由于在连接面上存在不完全接触, 因此严格来讲接触面上的传热是一个传导、对流、辐射并存的导热方式, 接触面两端对应点上的温度可能截然不同, 这是由于不同结构件的接触“热阻”造成的。因此, 在接触面上连续体温度场控制方程不适用。若不考虑这种接触“热阻”的影响, 用普通的连续体温度场分析方法进行飞机结构温度场分析, 将会造成重大的误差, 随之带来的是在进一步进行热应力分析时, 将无法计及由于接触面两端的温差而产生的面内应力。“飞机结构三维温度场分析”程序在大量热阻试验的基础上建立了热阻参数库[4], 从而可以较准确地计算出飞行器的内部温度。
1气动加热简述
当飞行器高速穿越大气层时, 由于气体的黏性作用, 靠近弹体表面的一薄层气流速度减小, 在弹体表面为零。这一薄层称为动力附面层, 简称附面层。在附面层中, 气体动能变为热能, 气体温度上升, 一部分热量传给飞行器壳体, 形成了气动加热。
1.1气动加热工程算法
根据沿气流方向的坐标X计算雷诺数:
当Rex≤2×106为层流工况;当Rex>2×106为湍流工况。
式 (1) 中:
ρH为大气密度;VH为气流速度;Prw为普朗特数;μH为大气密度;aH为音速。
层流工况下热流计算公式为:
式 (2) 中:μ*、ρ*是T*下的干空气物理参数。T*由下式确定:
湍流工况下热流计算公式为:
1.2烧蚀工程算法
烧蚀式防热是依靠防热材料的质量损失来消耗吸收气动加热以达到防热目的, 烧蚀计算方法都是针对具体的材料而建立的。本文采用文献[5]给出的方法进行烧蚀量的计算。
式中Tw为壁面温度;m﹒c为质量损失率。
2烧蚀外形计算和有限元模型建立
烧蚀量简化为物面径向缩聚量, 在建立模型时, 根据烧蚀量的最大可解性, 找寻头锥体内一不动点 (即不可能烧蚀到的点) 作为头锥建模的球心, 在球心向物面作射线与物面相交, 从不动点到头锥顶端的线段8等份, 以等分点到不动点之距为半径画弧与各射线成交点, 从而形成平面域网格 (图1) , 平面网格经过360°旋转则成体元网格。
物面点的烧蚀都沿径向缩聚, 我们的坐标变换方法是:把垂直物面的烧蚀边界的移动速度Vi转换到有限元计算的“球缺”物面点上, 若从t0到t (时间变化为Δt) 时, 其烧蚀量为ΔS, 则原OA长度缩减为OA′长度 (图2) , 原网格节点坐标在下一步计算时则以OA′长度重新等距计算形成, 进而重新形成总刚。
3 耦合计算实施
把烧蚀率计算程序嵌入飞机结构三维温度场分析程序之中, 首先进行气动/温度场耦合计算, 所得温度作为已知量计算当前的烧蚀量, 如果计算的烧蚀量大于0, 则根据烧蚀量计算当前的节点坐标, 在下一步温度迭代计算重新形成刚度矩阵之前读取新的节点坐标, 从而达到变边界的目的。二者如此交替迭代, 实现了直接耦合计算。
4 结论
数值计算结果表明, 弹箭的飞行速度、飞行时间、弹箭外形结构、弹体材料、飞行大气条件等均对烧蚀过程有影响, 这些影响参数的不同组合均会引起不同的烧蚀现象。
图3给出了不动点到头锥顶端点共9个节点的时历程曲线, 从图中可以看出, 温度场面向头锥内部逐渐降低。随着来流速度的增加, 经过相同的飞行时间, 锥部的烧蚀后退深度逐渐增加, 在来流速度较小时, 烧蚀后退深度亦较小, 烧蚀现象不明显, 而来流速度较大时, 则发生明显的烧蚀。这是因为随着来流速度的增大, 弹体表面与空气的摩擦加剧, 产生的热量增加;另一方面, 由于摩擦产生的热量增加, 使得弹体表面的温度升高, 导致化学反应的速率增大, 产生大量的化学热。
图4给出了导弹物面点的时间历程曲线, 温度是由头锥顶端点向后逐渐降低。
图5给出了头锥在某时刻的温度云图, 图6给出了头锥剖面下的烧蚀量, 头锥烧蚀后退深度与时间均近似呈线性关系增加。
摘要：针对高速飞行导弹建立了气动加热及烧蚀的数学计算模型, 应用“飞机结构三维温度场分析”程序, 根据物体面表面径烧蚀量移动边界, 通过“热阻”处理不同结构件的接触问题, 成功地进行了烧蚀和温度场的直接耦合求解。计算结果表明了导弹头锥的烧蚀和温度分布规律。
关键词：温度场,有限元,接触热阻,烧蚀,耦合计算
参考文献
[1]王中原, 史金光.超高速飞行弹箭气动烧蚀数值模拟研究.南京理工大学学报 (自然科学版) , 2003;27 (5) :595—602
[2]陈林泉, 张雁, 王书贤, 等.梯度功能材料的传热与烧蚀计算.科学技术与工程, 2004;4 (8) :663—677
[3]杨志斌, 成竹, 任青梅, 等.飞机结构三维温度场分析程序.计算力学学报, 2004, 21 (4) :498—501
[4]任青梅.航空结构热传导特性参数库.中国飞机强度研究所报告, 623S-9505-13.1995