两相流大作业.doc

水—水蒸汽两相相变界面的数值模拟 ——两相流动与热物理大作业姓名 张蛟龙_______学号 201328013524021__班级 物理308_____ 指导教师 刘捷__ 完成时间 _2014.5.8_ 水—水蒸汽两相相变界面的数值模拟报告一．文献综述 作为化石资源的替代产品，核能的高效，清洁一直备受青睐，然而光环之下，核废料的处理不禁让人黯然神伤强致命性辐射，动辄千年的半衰期，惯用的办法只能是深埋，等待下一代的聪明才智与此同时，核废料的利用和加速衰减一直是核能大国们的研究重点欧洲的ADS系统第六代散裂靶模型计划的目标就是要验证高水平的核废料转换的可行性散裂靶作为连接加速器和核废料的装置需要工作在高辐射和高热流密度的条件下，因此散裂靶的设计是ADS系统研制最有挑战的部分由加速器产生的高能质子流轰击靶核产生中子作为外源中子驱动和维持次临界堆的运行散裂靶在极小的空间内需承受极大的热负荷，质子束通道与靶核的自由面相邻更加剧了设计难度。

受材料限制，流体的温度不能超过550度，因此必须保证流体维持在一定的流量但同时又要考虑高流速带来的飞溅和回流造成的局部温度过高这一装置在水作为散裂靶的实验中获得了成功2． 问题描述2.1.模型及尺寸 图1、欧洲液态金属散裂靶V0.10示意图[1]如图1所示的欧洲加速器驱动次临界堆（ADS）之无窗散裂靶示意图，液态铅铋合金从上方管间流下并汇合，形成两相界面，质子束由中间的真空管进入打在自由面上此次模拟用的是水，详细物理背景见文献[1]2.2. 控制方程 连续性方程 动量方程 能量方程 3． Openfoam求解 有关Openfoam的下载和安装在老师给的安装指导的推荐网站上有详细的操作，在此就不赘述网址为：http://www.openfoam.org/download/ubuntu.php3.1. OpenFoam求解简述 Openfoam是一款基于linex的开源可编程软件，其求解过程的关键是三个文件夹的设置，即0，constant和system0文件夹里存放的是初始条件和边界条件设置文件；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和kqRWallFunction3.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。