两相流大作业.doc

水水蒸汽两相相变界面的数值模拟 两相流动与热物理大作业姓名 张蛟龙_学号 201328013524021_班级 物理308_ 指导教师 刘捷_ 完成时间 _2014.5.8_ 水水蒸汽两相相变界面的数值模拟报告一文献综述 作为化石资源的替代产品，核能的高效，清洁一直备受青睐，然而光环之下，核废料的处理不禁让人黯然神伤。强致命性辐射，动辄千年的半衰期，惯用的办法只能是深埋，等待下一代的聪明才智。与此同时，核废料的利用和加速衰减一直是核能大国们的研究重点。欧洲的ADS系统第六代散裂靶模型计划的目标就是要验证高水平的核废料转换的可行性。散裂靶作为连接加速器和核废料的装置需要工作在高辐射和高热流密度的条件下，因此散裂靶的设计是ADS系统研制最有挑战的部分。由加速器产生的高能质子流轰击靶核产生中子作为外源中子驱动和维持次临界堆的运行。散裂靶在极小的空间内需承受极大的热负荷，质子束通道与靶核的自由面相邻更加剧了设计难度。受材料限制，流体的温度不能超过550度，因此必须保证流体维持在一定的流量。但同时又要考虑高流速带来的飞溅和回流造成的局部温度过高。这一装置在水作为散裂靶的实验中获得了成功。2 问题描述2.1.模型及尺寸 图1、欧洲液态金属散裂靶V0.10示意图1如图1所示的欧洲加速器驱动次临界堆（ADS）之无窗散裂靶示意图，液态铅铋合金从上方管间流下并汇合，形成两相界面，质子束由中间的真空管进入打在自由面上。此次模拟用的是水，详细物理背景见文献1。2.2. 控制方程 连续性方程 动量方程 能量方程 3 Openfoam求解 有关Openfoam的下载和安装在老师给的安装指导的推荐网站上有详细的操作，在此就不赘述。网址为：/download/ubuntu.php。3.1. OpenFoam求解简述 Openfoam是一款基于linex的开源可编程软件，其求解过程的关键是三个文件夹的设置，即0，constant和system。0文件夹里存放的是初始条件和边界条件设置文件；constant文件夹里存放的是网格文件，物性参数和求解器模型；system文件夹里存放的是求解过程控制，差分格式和代数方程求解器设置文件。以下就三个文件的设置展开简述初始条件、边界条件、物性参数，网格个数、疏密设置差分格式、界面捕获算法、气蚀模型等的选择和设置。3.2. 0文件夹 包含有5个文件，分别为alph-water，p_rgh，U，epsilon，k，详细设置见附录1，这里只着重强调在大作业完成过程中几个曾经连续考虑的点。 首先是参数的量纲设置。在Openfoam文件中常会见到这样一行代码：dimensions 0 0 0 0 0，这便是量纲，单位顺序依次是 质量，长度，时间，温度，物质的量，电流，光强。 其次是边界条件和初始条件的设置。在alph-water中，alpha代表水所占比例，参照userguide，1时表示全部为液相，0时表示全部为气相。初始内部场的设置均为1，即起始时刻，散裂靶内部充满水。水入口是边界类型为“定值”，即fixedValue；上下两个出口的边界类型为“进出口”，即inletOutlet。在p_rgh中初始内部场设置值为5330，为压力较大的出口的压力值，查阅资料推荐内部场使用大值以减小汽蚀。质子束通道理论上应为真空，但考虑到自由面上水的汽化，选定了水的饱和压力2330以减小汽化，保证系统正常工作。文献显示上下两个出口压力差选定在2000到3000之间，此次试验选用差值为3000，即下面出口值设定为5330，也就是前面提到的出场值。在U中，因为速度时矢量，在设定时需要根据网格的情况处理。初始内部场设定为uniform（0 0 0），即均匀的静止态。进水口的类型为fixedValue，且有一个1m/s的向下流的速度。两个出口的类型为pressureInletOutletVelocity，即由压力决定速度的边界条件。在epsilon和k中，参照userguide中给出的公式：；， ，其中D为等效直径，分别为散裂靶上口总的直径和质子入口的直径，计算得k=0.0025，=0.00327。参照苏军伟的书，水入口类型为fixedValue(定值)，两个出水口的类型为inletOutlet，而墙的类型为epsilonWallFunction和kqRWallFunction。3.3. constant文件夹 首先介绍网格绘制。网格绘制有两种方法，在大作业探索过程中都进行了尝试。第一，关于Gambit直接绘制三维网格然后导进Openfoam进行计算。网格数大概为3万多，导入的操作过程如下：在Openfoam中建立一个文件夹，注意此文件夹必须包括system这个文件夹，因为导入过程需要system中所含的controlDict文件，否则将无法导入，可以先“借用”其中任意一个例子的system，后期计算时再根据实际情况更改。将网格文件复制到该文件夹下，在控制终端输入命令： fluentMeshToFoam ni.msh -scale 0.001其中fluentMeshToFoam为导入命令，ni.msh为网格名称，-scale 0.001的作用是把Gambit中默认为毫米的单位改成Openfoam中默认的以米为单位。单位转换的操作对初学者尤为要注意，开始的很多次尝试失败的原因就是忽略了两个软件之间单位的不相容性。导入完成后会生成Polymesh的文件夹，所包含的正是全部的网格信息。第二，关于Gambit绘制二维网格后导入Openfoam进行旋转。注意在Gambit绘制的二维网格在保存时一定要保存成二维的格式。与三位的导入相同，导入命令是： fluentMeshToFoam ni.msh -scale 0.001同样生成constant文件夹然后是旋转成五度的楔形的操作，命令是： makeAxialMesh -axis ax -wedge frontAndBackPlanes其中makeAxialMesh是生成轴网格的命令，-axis ax表示制定的旋转轴是边界名称型为ax的边，-wedge frontAndBackPlanes表示绕轴的旋转面是名称为frontAndBackPlanes的面。当然完成这一操作的前提是系统装有makeAxialMesh的软件。可以在控制终端直接输入makeAxialMesh这个命令，然后按提示进行下载和安装就可以使用。命令完成后终端会提示进行网格清理，命令为： collapseEdges在执行这一命令前，需要在system文件夹中添加一个collapseDict;文件规定网格清理的条件。具体设置见附录。之后可以检查一下网格，命令为： checkMesh这三个操作完成会生成一个新的包含网格信息的文件夹，它就可以作为新的0文件夹供计算使用。 然后是其余的几个文件。包括transportProperties，turbulenceProperties，RASProperties，g。在物性参数设置transportProperties中，水的饱和压力设定为2370Pa，表面张力sigma为0.07。流动选择湍流模型中的不可压缩RAS模型，选用k-两方程模型进行求解。在重力场设置时，因为也是矢量，所以要根据导进去的模型的坐标关系进行适当设定。具体设置见附录。3.4. system文件夹 system文件夹中一共包含了三个文件： controlDict，fvSchemes和fvSolution分别在求解过程控制，离散格式，求解器选择中发挥作用。首先看一下controlDict，文件中选定的是interPhaseChangeFoam，即带有自由表面流求解器。参照数值传热资料的设定Co=0.5，为了保证物理过程的清晰和碍于Co数的限制，时间步长按老师要求取作1e-5。参照前期同学的经验，采用adjustabeRunTime，即可调步长，自动调节最后一次的时间步长，以便准确输出。开启adjsutTimeStep，由于最大Co数的限制，当程序计算值大于最大Co数时，能够根据最大Co数自动调节减小步长；同时开启runtimeModifiable,在solver允许调节步长时，可根据最大Co数自动调节加大时间步长。 下面看数值格式，也就是各个量的离散形式。时间(ddtSchemes)的离散采用一阶隐式Euler形式。梯度格式（gradSchemes）采用有限元高斯线性插值。对流项的离散（divSchemes）格式选取特别要注意，由数值传热学的知识可以知道，二阶迎风有着优良的迁移和守恒特性，所以速度的离散采用具有二阶精度的二阶迎风，其他参数则采用一阶迎风格式。扩散项（laplacianSchemes）采用二阶高斯守恒格式。表面插值（interpolationSchemes）采用中心差分格式。其他的参数设置见附录。关于界面捕获方法，此次模拟采用的是VOF方法，即Volume of fluidmethod来捕获界面。接下来是代数方程求解器的设定。参照Openfoam当中的例子，Alph-water的计算选用PBiCG求解器，即预条件双共轭梯度求解器，残差（tolerance）设定为1e-8，相对残差（relTol）设定为0，当程序运行结果满足两个条件之一时停止。计算采用的是改进版的simple算法pimple，压力场与速度场分别存放在两套不同的网格中，所以压力求解器选用的是代数多重网格求解器（GAMG）。参照孙军委的指导书，U，k，全部选用的是光滑求解器（smoothSolver）。在pimple算法设定里，压力修正次数（nCorrectors）设为2，不进行非正交修正（nNonOrthogonalCorrectors）。其他设置见附录。4 计算流程和结果4.1. 运行界面以下为在1.3826s时的运行程序界面：DILUPBiCG: Solving for alpha.water, Initial residual = 1.01747e-06, Final residual = 9.65653e-09, No Iterations 1Phase-1 volume fraction = 0.581635 Min(alpha1) = 1.59102e-05 Max(alpha1) = 1MULES: Correcting alpha.waterLiquid phase volume fraction = 0.581612 Min(alpha1) = 1.59102e-05 Max(alpha1) = 1GAMGPCG: Solving for p_rgh, Initial residual = 0.000606079, Final residual = 5.81781e-09, No Iterations 8GAMGPCG: Solving for p_rgh, Initial residual = 8.81889e-06, Final residual = 3.87137e-09, No Iterations 5smoothSolver: Solving for epsilon, Initial residual = 0.000916708, Final residual = 9.80608e-07, No Iterations 33smoothSolver: Solving for k, Initial residual = 0.000708262, Final residual = 9.89131e-07, No Iterations 122ExecutionTime = 4.41 s ClockTime = 4 sCourant Number mean: 0.12362 max: 0.370546deltaT = 0.0002Time = 1.3826 我们可以看到运行程序里依次是对各个变量的求解，残差计算，迭代次数以及Co数和时间步长。另外，我们可以发现程序中每次运算都会计算Liquid phase volume fraction 等，这便是界面捕获VOF方法的体现。4.2. 运行结果及物理分析 以在Openfoam中旋转生成的网格为例，流动在1.5s前已经达到稳定。 以下为alpha-water的瞬态和稳态图 由图，我们可以看出，模拟开始后首先在散裂靶的水进口末端和突扩部分的顶端由于压力小于水的饱和压力而产生了汽化现象。稳态时在质子束通道和突扩的上半部分形成了稳定的近似真空区。以下为p_rgh的瞬态和稳态图上面第一个图是流动起始时压力的分布，可以看到，质子束口压力随着水面的下降稳定在2330Pa的设定值，水入口压力在流场中是最高的。由第二张图可以看出，出现汽化的部位的压力为水的饱和压力。最后的稳态图与alpha-water的稳态图几乎一致。以下为U的瞬态和稳态图 由图可以看出，在起始时刻，由于出口压力较低，在重力作用下，突扩管道内水的流速急剧增大，流速增大导致压力减小，局部出现汽化区，压力下降并稳定在饱和压力，液体又再次流回汽化区，如此往复作用，最终在重力和压力的联合作用下行程了稳态。理论上过程是这样的，但在模拟的过程中，一直没有得到理想的稳态速度分布图，还需要继续修改一些参数。五体会启示 老师依然是那个风格老师，作业却是更有挑战性的作业。这学期一直在学数值传热加上一个有关fluent的谈论课，所以对fluent已经不是那么陌生。不过，让人没想到的是怀柔这一片科研净土竟给了我们机会接触一下开源的Openfoam。没听过倒不打紧，关键是得在linux上运行。大作业开始的前期日子，我们一直在积极探索如何安装虚拟机和ubuntu。计算机本来内存就已经很小，运行虚拟机简直有些让人抓狂，再加上老师说得巨大的计算量，我们开始考虑放弃虚拟机，改向直接装并行系统。无意中得知班里有平时玩linux的高手，请来不到一个小时完成安装。一下子感觉前些天的努力都白搭了，看来老师提醒我们的要学会利用资源是很有必要的。接下来便是绘制网格，前期的计划是在较为熟悉的Gambit里绘制三维图形然后导如Openfoam。考虑到计算量的问题，并没有绘制完整的三维图形，而是在原模型的基础上截了一个五度的楔形用于模拟。网格绘制难度并不大，之后参照Openfoam的userguide，网格的导入同样是比较容易。之后就进入到了艰难的初始条件设置和求解设置过程，由于参照Openfoam中的例题system的设置大同小异，所以主要心思就放在了初始条件设置上，但在一系列尝试都失败后转而怀疑是网格的原因，所以才有了重新学习在Openfoam中旋转二维网格。网格的绘制总是比较简单，信心满满以为这次差不多了，但运行结果还是老样子，运行百步后Co数太大导致时间步长太小超过了计算机的浮点数，系统自动溢出。无奈之下，向理化所的王超同学请教。他认为可能是求解器的设置有问题，于是在依据他的建议进行了求解器设置的更改后计算才慢慢有了起色。以上便是本学期大作业的一个概述。 依然觉得设置大作业很有必要，并且是加大难度很有必要。因为事实证明，只要问题可解，我们总是可以找到解决的方法，只不过是前期做的功课的长短问题。首先，经过此次散裂靶的模拟，我们认识了一款Openfoam软件，在与其斗争的过程中，接触到了很多内部参数的设置选择，这是在fluent等软件计算时所看不到的。这样便为我们深入探究每个参数和设定对计算的影响提供了很好的平台尝试机会。然后就是这个学习的过程，当接触到一个完全陌生的软件或者问题时先应该从哪里下手，如何快速的找到自己需要的元素并加以利用；当在探索过程遇到障碍时如何有条理地进行分析，找到症结所在。最后要感谢老师又给我们提供了一次“痛然后快乐”的机会，感谢王超，范俊辉和邹文杰同学的无私帮助，感谢朱鹏同学分享的苏军伟老师的书，感谢陪我一起探索的栾奕军和丁林超同学。参考文献1 A.G. Class, D. Angeli, A. Batta, etc. XT-ADS Windowless spallation target thermohydraulic design & experimental setup, Journal of Nuclear Materials, 415, 378-384, 2011.2 A. Batta, A.G. Class, Numerical investigations on Geometrical Designs of the Windowless XT ADS spallation target, Proceedings of ICAPP 2007, Nice, France, May 13-18, 2007.3 A. Batta, A.G. Class, Free surface modeling and simulation of the water experiment for the XT ADS spallation target, Proceedings of ICAPP 08, Anaheim, CA USA, June 8-12, 2008.4 F. Roelofs, N.B. Siccama, H. Jeanmart, K. van Tichelen, M. Dierckx, P. Schuurmans, Application of a controlled swirl in the XT-ADS spallation target, Proceedings of ICAPP 08, Anaheim, CA USA, June 8-12, 2008 说明：in1为质子束通道的口，in2为水入口，ax为轴， frontAndBackPlanes是最初导进去的二维的旋转面，frontAndBackPlanes_pos与frontAndBackPlanes_neg为旋转后生成的楔形的前后两个面。附录一；初始条件设置alpha-waterdimensions 0 0 0 0 0;internalField uniform 1;boundaryField in2 type fixedValue; value $internalField; in1 type inletOutlet; inletValue uniform 1; value $internalField; out type inletOutlet; inletValue uniform 1; value $internalField; wall type zeroGradient; ax type symmetryPlane; frontAndBackPlanes type empty; frontAndBackPlanes_pos type wedge; frontAndBackPlanes_neg type wedge; epsiloninternalField uniform 0.00327;boundaryField wall type epsilonWallFunction; value uniform 0.00327; Cmu 0.09; kappa 0.41; E 9.8; value uniform 0.00327; out type inletOutlet; inletValue uniform 0.00327; value uniform 0.00327; in2 type fixedValue; value uniform 0.00327; in1 type inletOutlet; inletValue uniform 0.00327; value uniform 0.00327; ax type symmetryPlane; frontAndBackPlanes type empty; frontAndBackPlanes_pos type wedge; frontAndBackPlanes_neg type wedge; kinternalField uniform 0.0025;boundaryField wall type kqRWallFunction; value uniform 0.0025; out type inletOutlet; inletValue uniform 0.0025; value uniform 0.0025; in2 type fixedValue; value uniform 0.0025; in1 type inletOutlet; inletValue uniform 0.0025; value uniform 0.0025; ax type symmetryPlane; frontAndBackPlanes type empty; frontAndBackPlanes_pos type wedge; frontAndBackPlanes_neg type wedge; p_rghinternalField uniform 5330boundaryField in2 type zeroGradient; in1 type totalPressure; p0 uniform 2330; U U; phi phi; rho rho; psi none; gamma 1; value $internalField; out type totalPressure; p0 uniform 5330; U U; phi phi; rho rho; psi none; gamma 1; value $internalField; wall type zeroGradient; ax type symmetryPlane; frontAndBackPlanes type empty; frontAndBackPlanes_pos type wedge; frontAndBackPlanes_neg type wedge; uinternalField uniform (0 0 0);boundaryField in2 type fixedValue; value uniform (0 -1 0); in1 type pressureInletOutletVelocity; phi phi; value $internalField; out type pressureInletOutletVelocity; phi phi; value $internalField; wall type fixedValue; value uniform (0 0 0); ax type symmetryPlane; frontAndBackPlanes type empty; frontAndBackPlanes_pos type wedge; frontAndBackPlanes_neg type wedge; 附录二：constant设置RASPropertiesRASModel kEpsilon;turbulence on;printCoeffs on;transportPropertiesphases (water vapour);phaseChangeTwoPhaseMixture SchnerrSauer;pSat pSat 1 -1 -2 0 0 2370;sigma sigma 1 0 -2 0 0 0 0 0.07; water transportModel Newtonian; nu nu 0 2 -1 0 0 0 0 9.9838e-07; rho rho 1 -3 0 0 0 0 0 998.12;vapour transportModel Newtonian; nu nu 0 2 -1 0 0 0 0 5.5526e-04; rho rho 1 -3 0 0 0 0 0 0.017529;KunzCoeffs UInf UInf 0 1 -1 0 0 0 0 1.0; tInf tInf 0 0 1 0 0 0 0 0.5; Cc Cc 0 0 0 0 0 0 0 10; Cv Cv 0 0 0 0 0 0 0 10;MerkleCoeffs UInf UInf 0 1 -1 0 0 0 0 1.0; tInf tInf 0 0 1 0 0 0 0 0.5; Cc Cc 0 0 0 0 0 0 0 80; Cv Cv 0 0 0 0 0 0 0 1e-03;SchnerrSauerCoeffs n n 0 -3 0 0 0 0 0 1.0e+12; dNuc dNuc 0 1 0 0 0 0 0 2.0e-06; Cc Cc 0 0 0 0 0 0 0 1; Cv Cv 0 0 0 0 0 0 0 1; turbulenceProperties simulationType RASModel;附录三：system设置controlDictapplication interPhaseChangeFoam;startTime 0;stopAt endTime;endTime 5;deltaT 5e-5;writeControl adjustableRunTime;writeInterval 0.001;purgeWrite 0;writeFormat ascii;writePrecision 6;writeCompression uncompressed;timeFormat general;runTimeModifiable yes;adjustTimeStep on;maxCo 0.5;maxAlphaCo0.5;maxDeltaT1; fvSchemesddtSchemes default Euler;gradSchemes default cellLimited Gauss linear 1.0;divSchemes default none; div(phi,alpha) Gauss u