1、水水蒸汽两相相变界面的数值模拟两相流动与热物理大作业姓名 张蛟龙_学号 201328013524021_班级 物理 308_指导教师 刘捷_ 完成时间 _2014.5.8_ 水水蒸汽两相相变界面的数值模拟报告一文献综述作为化石资源的替代产品，核能的高效，清洁一直备受青睐，然而光环之下，核废料的处理不禁让人黯然神伤。强致命性辐射，动辄千年的半衰期，惯用的办法只能是深埋，等待下一代的聪明才智。与此同时，核废料的利用和加速衰减一直是核能大国们的研究重点。欧洲的 ADS 系统第六代散裂靶模型计划的目标就是要验证高水平的核废料转换的可行性。散裂靶作为连接加速器和核废料的装置需要工作在高辐射和高热流密度的
2、条件下，因此散裂靶的设计是ADS 系统研制最有挑战的部分。由加速器产生的高能质子流轰击靶核产生中子作为外源中子驱动和维持次临界堆的运行。散裂靶在极小的空间内需承受极大的热负荷，质子束通道与靶核的自由面相邻更加剧了设计难度。受材料限制，流体的温度不能超过 550 度，因此必须保证流体维持在一定的流量。但同时又要考虑高流速带来的飞溅和回流造成的局部温度过高。这一装置在水作为散裂靶的实验中获得了成功。2问题描述2.1.模型及尺寸图 1、欧洲液态金属散裂靶 V0.10 示意图1如图 1 所示的欧洲加速器驱动次临界堆（ADS）之无窗散裂靶示意图，液态铅铋合金从上方管间流下并汇合，形成两相界面，质子束由中
3、间的真空管进入打在自由面上。此次模拟用的是水，详细物理背景见文献1。2.2. 控制方程 连续性方程动量方程 能量方程 3 Openfoam 求解有关 Openfoam 的下载和安装在老师给的安装指导的推荐网站上有详细的操作，在此就不赘述。网址为：http:/www.openfoam.org/download/ubuntu.php。3.1. OpenFoam 求解简述Openfoam 是一款基于 linex 的开源可编程软件，其求解过程的关键是三个文件夹的设置，即 0，constant 和 system。0 文件夹里存放的是初始条件和边界条件设置文件；constant 文件夹里存放的是网格文件，
4、物性参数和求解器模型；system 文件夹里存放的是求解过程控制，差分格式和代数方程求解器设置文件。以下就三个文件的设置展开简述初始条件、边界条件、物性参数，网格个数、疏密设置差分格式、界面捕获算法、气蚀模型等的选择和设置。3.2. 0 文件夹 包含有 5 个文件，分别为 alph-water，p_rgh，U ，epsilon，k，详细设置见附录 1，这里只着重强调在大作业完成过程中几个曾经连续考虑的点。首先是参数的量纲设置。在 Openfoam 文件中常会见到这样一行代码：dimensions 0 0 0 0 0，这便是量纲，单位顺序依次是 质量，长度，时间，温度，物质的量，电流，光强。其次
5、是边界条件和初始条件的设置。在 alph-water 中，alpha 代表水所占比0)(iiuxtibjiijii Fxpt jb jcijkijijjiuFQueuxet )21()21(例，参照 userguide，1 时表示全部为液相， 0 时表示全部为气相。初始内部场的设置均为 1，即起始时刻，散裂靶内部充满水。水入口是边界类型为“定值” ，即 fixedValue；上下两个出口的边界类型为 “进出口” ，即 inletOutlet。在p_rgh 中初始内部场设置值为 5330，为压力较大的出口的压力值，查阅资料推荐内部场使用大值以减小汽蚀。质子束通道理论上应为真空，但考虑到自由面上水
6、的汽化，选定了水的饱和压力 2330 以减小汽化，保证系统正常工作。文献显示上下两个出口压力差选定在 2000 到 3000 之间，此次试验选用差值为3000，即下面出口值设定为 5330，也就是前面提到的出场值。在 U 中，因为速度时矢量，在设定时需要根据网格的情况处理。初始内部场设定为 uniform（0 0 0） ，即均匀的静止态。进水口的类型为 fixedValue，且有一个 1m/s 的向下流的速度。两个出口的类型为 pressureInletOutletVelocity，即由压力决定速度的边界条件。在 epsilon 和 k 中，参照 userguide 中给出的公式：； ， ，其
7、中 D 为等效)(212 zyxUklkC5.17.0U0l07.直径， ， 分别为散裂靶上口总的直径和质子入口的直径，21D12计算得 k=0.0025， =0.00327。参照苏军伟的书，水入口类型为 fixedValue(定值)，两个出水口的类型为 inletOutlet，而墙的类型为 epsilonWallFunction 和kqRWallFunction。3.3. constant 文件夹首先介绍网格绘制。网格绘制有两种方法，在大作业探索过程中都进行了尝试。第一，关于 Gambit 直接绘制三维网格然后导进 Openfoam 进行计算。网格数大概为 3 万多，导入的操作过程如下：在
8、Openfoam 中建立一个文件夹，注意此文件夹必须包括 system 这个文件夹，因为导入过程需要 system 中所含的controlDict 文件，否则将无法导入，可以先“借用”其中任意一个例子的system，后期计算时再根据实际情况更改。将网格文件复制到该文件夹下，在控制终端输入命令：fluentMeshToFoam ni.msh -scale 0.001其中 fluentMeshToFoam 为导入命令，ni.msh 为网格名称， -scale 0.001 的作用是把 Gambit 中默认为毫米的单位改成 Openfoam 中默认的以米为单位。单位转换的操作对初学者尤为要注意，开始的
9、很多次尝试失败的原因就是忽略了两个软件之间单位的不相容性。导入完成后会生成 Polymesh 的文件夹，所包含的正是全部的网格信息。第二，关于 Gambit 绘制二维网格后导入 Openfoam 进行旋转。注意在 Gambit 绘制的二维网格在保存时一定要保存成二维的格式。与三位的导入相同，导入命令是：fluentMeshToFoam ni.msh -scale 0.001同样生成 constant 文件夹然后是旋转成五度的楔形的操作，命令是：makeAxialMesh -axis ax -wedge frontAndBackPlanes其中 makeAxialMesh 是生成轴网格的命令，-
10、axis ax 表示制定的旋转轴是边界名称型为 ax 的边， -wedge frontAndBackPlanes 表示绕轴的旋转面是名称为frontAndBackPlanes 的面。当然完成这一操作的前提是系统装有 makeAxialMesh的软件。可以在控制终端直接输入 makeAxialMesh 这个命令，然后按提示进行下载和安装就可以使用。命令完成后终端会提示进行网格清理，命令为：collapseEdges在执行这一命令前，需要在 system 文件夹中添加一个 collapseDict;文件规定网格清理的条件。具体设置见附录。之后可以检查一下网格，命令为：checkMesh这三个操作完
11、成会生成一个新的包含网格信息的文件夹，它就可以作为新的 0文件夹供计算使用。然后是其余的几个文件。包括transportProperties，turbulenceProperties，RASProperties，g。在物性参数设置transportProperties 中，水的饱和压力设定为 2370Pa，表面张力 sigma 为 0.07。流动选择湍流模型中的不可压缩 RAS 模型，选用 k- 两方程模型进行求解。在重力场设置时，因为也是矢量，所以要根据导进去的模型的坐标关系进行适当设定。具体设置见附录。3.4. system 文件夹system 文件夹中一共包含了三个文件： control
12、Dict，fvSchemes 和fvSolution 分别在求解过程控制，离散格式，求解器选择中发挥作用。首先看一下 controlDict，文件中选定的是 interPhaseChangeFoam，即带有自由表面流求解器。参照数值传热资料的设定 Co=0.5，为了保证物理过程的清晰和碍于 Co 数的限制，时间步长按老师要求取作 1e-5。参照前期同学的经验，采用adjustabeRunTime，即可调步长，自动调节最后一次的时间步长，以便准确输出。开启 adjsutTimeStep，由于最大 Co 数的限制，当程序计算值大于最大 Co数时，能够根据最大 Co 数自动调节减小步长；同时开启 r
13、untimeModifiable,在solver 允许调节步长时，可根据最大 Co 数自动调节加大时间步长。下面看数值格式，也就是各个量的离散形式。时间(ddtSchemes)的离散采用一阶隐式 Euler 形式。梯度格式（gradSchemes）采用有限元高斯线性插值。对流项的离散（divSchemes）格式选取特别要注意，由数值传热学的知识可以知道，二阶迎风有着优良的迁移和守恒特性，所以速度的离散采用具有二阶精度的二阶迎风，其他参数则采用一阶迎风格式。扩散项（laplacianSchemes）采用二阶高斯守恒格式。表面插值（interpolationSchemes）采用中心差分格式。其他的
14、参数设置见附录。关于界面捕获方法，此次模拟采用的是 VOF 方法，即Volume of fluidmethod 来捕获界面。接下来是代数方程求解器的设定。参照Openfoam 当中的例子， Alph-water 的计算选用 PBiCG 求解器，即预条件双共轭梯度求解器，残差（tolerance）设定为 1e-8，相对残差（ relTol）设定为 0，当程序运行结果满足两个条件之一时停止。计算采用的是改进版的 simple 算法pimple，压力场与速度场分别存放在两套不同的网格中，所以压力求解器选用的是代数多重网格求解器（GAMG） 。参照孙军委的指导书，U ，k， 全部选用的是光滑求解器（s
15、moothSolver） 。在 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) =
16、 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-0
17、6, 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 ClockTim
18、e = 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 的瞬态和稳态图由图，我们可以看出，模拟开始后首先在散裂靶的水进口末端和突扩部
19、分的顶端由于压力小于水的饱和压力而产生了汽化现象。稳态时在质子束通道和突扩的上半部分形成了稳定的近似真空区。以下为 p_rgh 的瞬态和稳态图上面第一个图是流动起始时压力的分布，可以看到，质子束口压力随着水面的下降稳定在 2330Pa 的设定值，水入口压力在流场中是最高的。由第二张图可以看出，出现汽化的部位的压力为水的饱和压力。最后的稳态图与 alpha-water 的稳态图几乎一致。以下为 U 的瞬态和稳态图由图可以看出，在起始时刻，由于出口压力较低，在重力作用下，突扩管道内水的流速急剧增大，流速增大导致压力减小，局部出现汽化区，压力下降并稳定在饱和压力，液体又再次流回汽化区，如此往复作用，
20、最终在重力和压力的联合作用下行程了稳态。理论上过程是这样的，但在模拟的过程中，一直没有得到理想的稳态速度分布图，还需要继续修改一些参数。五体会启示老师依然是那个风格老师，作业却是更有挑战性的作业。这学期一直在学数值传热加上一个有关 fluent 的谈论课，所以对 fluent 已经不是那么陌生。不过，让人没想到的是怀柔这一片科研净土竟给了我们机会接触一下开源的Openfoam。没听过倒不打紧，关键是得在 linux 上运行。大作业开始的前期日子，我们一直在积极探索如何安装虚拟机和 ubuntu。计算机本来内存就已经很小，运行虚拟机简直有些让人抓狂，再加上老师说得巨大的计算量，我们开始考虑放弃虚
21、拟机，改向直接装并行系统。无意中得知班里有平时玩 linux 的高手，请来不到一个小时完成安装。一下子感觉前些天的努力都白搭了，看来老师提醒我们的要学会利用资源是很有必要的。接下来便是绘制网格，前期的计划是在较为熟悉的 Gambit 里绘制三维图形然后导如 Openfoam。考虑到计算量的问题，并没有绘制完整的三维图形，而是在原模型的基础上截了一个五度的楔形用于模拟。网格绘制难度并不大，之后参照 Openfoam 的 userguide，网格的导入同样是比较容易。之后就进入到了艰难的初始条件设置和求解设置过程，由于参照 Openfoam 中的例题 system 的设置大同小异，所以主要心思就放
22、在了初始条件设置上，但在一系列尝试都失败后转而怀疑是网格的原因，所以才有了重新学习在 Openfoam 中旋转二维网格。网格的绘制总是比较简单，信心满满以为这次差不多了，但运行结果还是老样子，运行百步后 Co 数太大导致时间步长太小超过了计算机的浮点数，系统自动溢出。无奈之下，向理化所的王超同学请教。他认为可能是求解器的设置有问题，于是在依据他的建议进行了求解器设置的更改后计算才慢慢有了起色。以上便是本学期大作业的一个概述。依然觉得设置大作业很有必要，并且是加大难度很有必要。因为事实证明，只要问题可解，我们总是可以找到解决的方法，只不过是前期做的功课的长短问题。首先，经过此次散裂靶的模拟，我们
23、认识了一款 Openfoam 软件，在与其斗争的过程中，接触到了很多内部参数的设置选择，这是在 fluent 等软件计算时所看不到的。这样便为我们深入探究每个参数和设定对计算的影响提供了很好的平台尝试机会。然后就是这个学习的过程，当接触到一个完全陌生的软件或者问题时先应该从哪里下手，如何快速的找到自己需要的元素并加以利用；当在探索过程遇到障碍时如何有条理地进行分析，找到症结所在。最后要感谢老师又给我们提供了一次“痛然后快乐”的机会，感谢王超，范俊辉和邹文杰同学的无私帮助，感谢朱鹏同学分享的苏军伟老师的书，感谢陪我一起探索的栾奕军和丁林超同学。参考文献1 A.G. Class, D. Angel
24、i, A. Batta, etc. XT-ADS Windowless spallation target thermohydraulic design internalField uniform 1;boundaryFieldin2type fixedValue;value $internalField;in1type inletOutlet;inletValue uniform 1;value $internalField;outtype inletOutlet;inletValue uniform 1;value $internalField;walltype zeroGradient;
25、axtype symmetryPlane;frontAndBackPlanestype empty;frontAndBackPlanes_postype wedge;frontAndBackPlanes_negtype wedge;epsiloninternalField uniform 0.00327;boundaryFieldwalltype epsilonWallFunction;value uniform 0.00327;Cmu 0.09;kappa 0.41;E 9.8;value uniform 0.00327;outtype inletOutlet;inletValue unif
26、orm 0.00327;value uniform 0.00327;in2type fixedValue;value uniform 0.00327;in1type inletOutlet;inletValue uniform 0.00327;value uniform 0.00327;axtype symmetryPlane;frontAndBackPlanestype empty;frontAndBackPlanes_postype wedge;frontAndBackPlanes_negtype wedge;kinternalField uniform 0.0025;boundaryFi
27、eldwalltype kqRWallFunction;value uniform 0.0025;outtype inletOutlet;inletValue uniform 0.0025;value uniform 0.0025;in2type fixedValue;value uniform 0.0025;in1type inletOutlet;inletValue uniform 0.0025;value uniform 0.0025;axtype symmetryPlane;frontAndBackPlanestype empty;frontAndBackPlanes_postype
28、wedge;frontAndBackPlanes_negtype wedge;p_rghinternalField uniform 5330boundaryFieldin2type zeroGradient;in1type totalPressure;p0 uniform 2330;U U;phi phi;rho rho;psi none;gamma 1;value $internalField;outtype totalPressure;p0 uniform 5330;U U;phi phi;rho rho;psi none;gamma 1;value $internalField;wall
29、type zeroGradient;axtype symmetryPlane;frontAndBackPlanestype empty;frontAndBackPlanes_postype wedge;frontAndBackPlanes_negtype wedge;uinternalField uniform (0 0 0);boundaryFieldin2type fixedValue;value uniform (0 -1 0);in1type pressureInletOutletVelocity;phi phi; value $internalField;outtype pressu
30、reInletOutletVelocity;phi phi;value $internalField;walltype fixedValue;value uniform (0 0 0);axtype symmetryPlane;frontAndBackPlanestype empty;frontAndBackPlanes_postype wedge;frontAndBackPlanes_negtype wedge;附录二：constant 设置RASPropertiesRASModel kEpsilon;turbulence on;printCoeffs on;transportPropert
31、iesphases (water vapour);phaseChangeTwoPhaseMixture SchnerrSauer;pSat pSat 1 -1 -2 0 0 2370;sigma sigma 1 0 -2 0 0 0 0 0.07; watertransportModel Newtonian;nu nu 0 2 -1 0 0 0 0 9.9838e-07;rho rho 1 -3 0 0 0 0 0 998.12;vapourtransportModel Newtonian;nu nu 0 2 -1 0 0 0 0 5.5526e-04;rho rho 1 -3 0 0 0 0
32、 0 0.017529;KunzCoeffsUInf 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;MerkleCoeffsUInf 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;SchnerrSauerCoeffsn n 0 -3 0 0 0 0 0 1.0e+12;dNuc dNuc
33、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; turbulencePropertiessimulationType RASModel;附录三：system 设置controlDictapplication interPhaseChangeFoam;startTime 0;stopAt endTime;endTime 5;deltaT 5e-5;writeControl adjustableRunTime;writeInterval 0.001;purgeWrite 0;writeFormat ascii;
34、writePrecision 6;writeCompression uncompressed;timeFormat general;runTimeModifiable yes;adjustTimeStep on;maxCo 0.5;maxAlphaCo 0.5;maxDeltaT 1;fvSchemesddtSchemesdefault Euler;gradSchemes default cellLimited Gauss linear 1.0;divSchemesdefault none;div(phi,alpha) Gauss upwind; div(phirb,alpha) Gauss
35、upwind ;div(rhoPhi,U) Gauss linearUpwind grad(U);div(phi,k) Gauss upwind;div(phi,epsilon) Gauss upwind;div(muEff*dev(T(grad(U) Gauss linear;laplacianSchemesdefault Gauss linear corrected;interpolationSchemesdefault linear;snGradSchemesdefault corrected;fluxRequireddefault no;p_rgh;pcorr;alpha.water;
36、fvSolutionsolvers“alpha.water.*“cAlpha 0;nAlphaCorr 1;nAlphaSubCycles 1;MULESCorr yes; nLimiterIter 5; solver PBiCG;preconditioner DILU;tolerance 1e-8;relTol 0;maxIter 10;minIter 1;pcorrsolver PCG; preconditionerpreconditioner GAMG;tolerance 1e-4;relTol 0;smoother DICGaussSeidel;nPreSweeps 0;nPostSw
37、eeps 2;nFinestSweeps 2;cacheAgglomeration true;nCellsInCoarsestLevel 10;agglomerator faceAreaPair;mergeLevels 1;tolerance 1e-4; relTol 0;maxIter 100;p_rgh$pcorr;tolerance 1e-8;relTol 0;p_rghFinal$p_rgh;tolerance 1e-8;relTol 0;“(U|k|epsilon)“solver smoothSolver;smoother GaussSeidel;tolerance 1e-6;relTol 0;nSweeps 1;“(U|k|epsilon)Final“$U;relTol 0;PIMPLEmomentumPredictor no;nOuterCorrectors 1; nCorrectors 2; nNonOrthogonalCorrectors 0;nAlphaCorr 1;nAlphaSubCycles 3;cAlpha 1;correctPhi yes;relaxationFactors“ .*“ 1;
