计算机模拟实验第6章 式(6-7-14)在离散化的有限体积V上在时间(，1十△)范围内积分有 [品os)av+e)aw]山=[n.rgaw+s(aw]u （6-7-15) 控制方程的离散化就是对式(67-15)中的各个项进行离散。式(6-7-15)中的瞬态项 采用欧拉隐式时间差分方案离散化，可表示为： -(). (6-7-16) 这一离散化方案在时间上是一阶准确的，并且保证了有界性和无条件稳定性。 式(6-7-15)中对流项的离散化是先用高斯定理将体积分转换为面积分，再进行数值 离散，对流项的离散表达式为： .(p90)dW=2S。·(pU)。≈2s。·(oU)g=∑F。(6-7-17) 式中，F=S。·()。代表通过面b的质量通量。为了计算对流项的离散表达式，需要给 出有限控制体积界面处的物理量中。，该值可通过与界面相邻的两个有限控制体积中对应 的物理量插值求出。常用的插值方法有中心插值、上风插值和复合插值等，不同的插值方 法具有不同的精度和稳定性。 式(6-7-15)中扩散项的离散化类似于对流项的离散化，扩散项的离散表达式为： v.(pr7)dW=2S·(pP)≈∑(p).Se·()(6-7-18) 式中，（）。可认为是线性变化的，也可通过插值求出。在两相流体流动中需要考虑两相 界面方向与单元网格边界方向的差异，该差异代表了NS方程中的黏度，因此在两相流体 流动模拟中有： (T)=%(T)+(1-)(T) (6-7-19) 式(6-7-15)中的广义源项S()包括所有未被表述为对流，扩散或时间量的项。广义 源项可为·的一般函数。在单元网格有限体积上线性化并积分后的离散表达式为： S(p)dV≈MVrp+MprV (6-7-20) 式中，M和M,为线性化过程中生成的两个系数。 (3)离散方程组的求解。 控制方程离散化后，在每个有限控制体积上都会得到一组线性代数方程，虽然每个方 程的具体形式都不同，但均可写成一般形式为： argrut+Saxgyut =c (6-7-21) 式(6-7-21)可用矩阵形式表示为： (6-7-22) 式中，A表示对应于每个有限体积的稀疏矩阵，包括对角线上的元素系数ar和非对角线 上的元素系数aN:中为求解变量·的矩阵：C为包含源项的矩阵。 线性代数方程组(6-7-22)中稀疏矩阵A通常无法直接求出其逆矩阵A1,因此一般 应用选代方法来求出中的近似解。迭代方法从一个初始的猜测值开始，通过方程逐步迭 -345
— ３４５ — 式(６Ｇ７Ｇ１４)在离散化的有限体积VP 上在时间(t,t＋Δt)范围内积分有: ∫ t＋Δt t ∂ ∂t∫VP (ρϕ)dV ＋∫VP Ñ(ρUϕ)dV é ë ê ê ù û ú údt＝∫ t＋Δt t ∫VP Ñ(Γ Ñϕ)dV ＋∫VP S(ϕ)dV é ë ê ê ù û ú údt (６Ｇ７Ｇ１５) 控制方程的离散化就是对式(６Ｇ７Ｇ１５)中的各个项进行离散.式(６Ｇ７Ｇ１５)中的瞬态项 采用欧拉隐式时间差分方案离散化,可表示为: ∂ ∂t∫VP (ρϕ)dV ≈ (ρPϕPVP )(t＋Δt)－ (ρPϕPVP )t Δt (６Ｇ７Ｇ１６) 这一离散化方案在时间上是一阶准确的,并且保证了有界性和无条件稳定性. 式(６Ｇ７Ｇ１５)中对流项的离散化是先用高斯定理将体积分转换为面积分,再进行数值 离散,对流项的离散表达式为: ∫VP Ñ(ρϕU)dV ＝∑b Sb(ρϕU)b ≈ ∑b Sb(ρU)bφb ＝∑b Fϕb (６Ｇ７Ｇ１７) 式中,F＝Sb(ρU)b 代表通过面b 的质量通量.为了计算对流项的离散表达式,需要给 出有限控制体积界面处的物理量ϕb,该值可通过与界面相邻的两个有限控制体积中对应 的物理量插值求出.常用的插值方法有中心插值、上风插值和复合插值等,不同的插值方 法具有不同的精度和稳定性. 式(６Ｇ７Ｇ１５)中扩散项的离散化类似于对流项的离散化,扩散项的离散表达式为: ∫VP Ñ(ρΓÑϕ)dV ＝∑b Sb(ρΓÑϕ)b ≈ ∑b (ρΓ)bSb(Ñϕ)b (６Ｇ７Ｇ１８) 式中,(Γ)b 可认为是线性变化的,也可通过插值求出.在两相流体流动中需要考虑两相 界面方向与单元网格边界方向的差异,该差异代表了 NＧS方程中的黏度,因此在两相流体 流动模拟中有: (Γ)b＝ηb (Γ)P b ＋(１－ηb)(Γ)N b (６Ｇ７Ｇ１９) 式(６Ｇ７Ｇ１５)中的广义源项S(ϕ)包括所有未被表述为对流、扩散或时间量的项.广义 源项可为ϕ 的一般函数.在单元网格有限体积上线性化并积分后的离散表达式为: ∫VP S(ϕ)dV ≈ M１VP ＋M２ϕPVP (６Ｇ７Ｇ２０) 式中,M１和 M２为线性化过程中生成的两个系数. (３)离散方程组的求解. 控制方程离散化后,在每个有限控制体积上都会得到一组线性代数方程,虽然每个方 程的具体形式都不同,但均可写成一般形式为: aPϕP(t＋Δt)＋ ∑N aNϕN(t＋Δt)＝c (６Ｇ７Ｇ２１) 式(６Ｇ７Ｇ２１)可用矩阵形式表示为: AΦ＝C (６Ｇ７Ｇ２２) 式中,A 表示对应于每个有限体积的稀疏矩阵,包括对角线上的元素系数aP 和非对角线 上的元素系数aN ;Φ 为求解变量ϕ 的矩阵;C 为包含源项的矩阵. 线性代数方程组(６Ｇ７Ｇ２２)中稀疏矩阵A 通常无法直接求出其逆矩阵A－１,因此一般 应用迭代方法来求出Φ 的近似解.迭代方法从一个初始的猜测值开始,通过方程逐步迭
物理实验教程 一近代物理实验 代，直到两个解之间的差异达到规定精度为止。 三、液滴撞击壁面问题的OpenFOAM模拟方法 OpenFOAM是一个基于C十十程序设计语言编写的面向对象的开源CFD软件库】 般需要在Linux系统下运行。OpenFOAM的主要优点有：源代码开源，研究者可根据 自己的需要在开源协议的基础上自主修改使用这些代码：软件架构设计优越，采用面向对 象编程技术，研究者可根据需要方便植入新的模型：模型设计、网格划分和导人十分方便： 从底层支持并行计算，计算效率高：配套的辅助工具库多，如swak4Foam和Paraview,可 方便实现计算数据的后处理和可视化。 用OpenFOAM软件模拟液滴撞击壁面的过程，可方便地完成计算区域网格划分、控 方程的有限体积法离散化、复杂代数方程组求解等一系列过程。模拟的只体实现方法 主要句括下列五个步骤： (1)物理模型建立。 液滴撞击壁面问题可简化为一个二维物理模型：单个液滴在重力的作用下从一定高 度自然下落撞击水平固体壁面，即研究液滴撞击水平壁面过程的运动状态及液滴自由界 面的变化规律。 (2)计算风域设定与网格划分 计算区域可设定为一个长方形，其大小可根据实际模拟问题确定，一般应不小于液滴 大小的10倍。长方形的四条边为计算区域的边界，底部边界为固体壁面边界，其余边界 为大气边界。基于边界规则，整个计算区域可划分规则网格，网格数量应足够多，特别是 固体壁面边界的网格需要加密，以便能够精确地描述液滴撞击壁面过程的运动状态。根 据实际模拟液滴的大小，液滴初始形状设置为计算区域中一个圆形的区域，液滴形状在模 拟计算过程中会实时变化。 (3)边界条件与控制参数设置 固体壁面边界可设置为无滑移，润湿性边界条件，其余边界可设置为大气边界条件。 控制参数需要根据大气和液滴物性以及实际模拟条件设置。 （4)求解器选取。 为了能够模拟两相流体流动，可选取OpenFOAM软件库中自带的interFOAM求解 器，开展液滴撞击壁面过程的模拟计算。 (5)计算结果可视化处理。 根据模拟计算结果，可绘制液滴运动和撞击壁面过程中的液体分布、压力场和速度场 图像及相关曲线，用于分析液滴运动的变化特征 【实验内容】 一、基础性实验内容 L.液滴撞击壁面模拟的网格参数优化 选取液体为水、气体为空气，设置液滴初始直径为2.68mm,初始速度为0.678m/s, 接触角为60°。改变划分网格的大小，多次模拟计算，并通过后处理软件展示液滴撞击壁 面过程，分析网格大小对计算结果，计算时间的影响，优化出最佳的网格划分方案。 -346
— ３４６ — 代,直到两个解之间的差异达到规定精度为止. 三、液滴撞击壁面问题的 OpenFOAM 模拟方法 OpenFOAM 是一个基于 C＋＋程序设计语言编写的面向对象的开源 CFD 软件库, 一般需要在 Linux系统下运行.OpenFOAM 的主要优点有:源代码开源,研究者可根据 自己的需要在开源协议的基础上自主修改使用这些代码;软件架构设计优越,采用面向对 象编程技术,研究者可根据需要方便植入新的模型;模型设计、网格划分和导入十分方便; 从底层支持并行计算,计算效率高;配套的辅助工具库多,如swak４Foam 和 Paraview,可 方便实现计算数据的后处理和可视化. 用 OpenFOAM 软件模拟液滴撞击壁面的过程,可方便地完成计算区域网格划分、控 制方程的有限体积法离散化、复杂代数方程组求解等一系列过程.模拟的具体实现方法 主要包括下列五个步骤: (１)物理模型建立. 液滴撞击壁面问题可简化为一个二维物理模型:单个液滴在重力的作用下从一定高 度自然下落撞击水平固体壁面,即研究液滴撞击水平壁面过程的运动状态及液滴自由界 面的变化规律. (２)计算区域设定与网格划分. 计算区域可设定为一个长方形,其大小可根据实际模拟问题确定,一般应不小于液滴 大小的１０倍.长方形的四条边为计算区域的边界,底部边界为固体壁面边界,其余边界 为大气边界.基于边界规则,整个计算区域可划分规则网格,网格数量应足够多,特别是 固体壁面边界的网格需要加密,以便能够精确地描述液滴撞击壁面过程的运动状态.根 据实际模拟液滴的大小,液滴初始形状设置为计算区域中一个圆形的区域,液滴形状在模 拟计算过程中会实时变化. (３)边界条件与控制参数设置. 固体壁面边界可设置为无滑移、润湿性边界条件,其余边界可设置为大气边界条件. 控制参数需要根据大气和液滴物性以及实际模拟条件设置. (４)求解器选取. 为了能够模拟两相流体流动,可选取 OpenFOAM 软件库中自带的interFOAM 求解 器,开展液滴撞击壁面过程的模拟计算. (５)计算结果可视化处理. 根据模拟计算结果,可绘制液滴运动和撞击壁面过程中的液体分布、压力场和速度场 图像及相关曲线,用于分析液滴运动的变化特征. 【实验内容】 一、基础性实验内容 １ 液滴撞击壁面模拟的网格参数优化 选取液体为水、气体为空气,设置液滴初始直径为２６８mm,初始速度为０６７８m/s, 接触角为６０°.改变划分网格的大小,多次模拟计算,并通过后处理软件展示液滴撞击壁 面过程,分析网格大小对计算结果、计算时间的影响,优化出最佳的网格划分方案
计算机模拟实验第6章 2。液滴撞击壁面运动特征的模拟分析 在最佳的网格划分方案条件下，改变液滴的大小和初速度以及液体的黏度等影响参 数，模拟计算液滴撞击壁面过程，并分析不同参数条件下液滴运动的特征。 二、设计性实验内容 在结合现有模拟计算条件的基础上，探究固体壁面粗糙度和表面形状，壁面润湿性对 液滴撞击壁面运动状态的影响，并分析液滴运动规律。 【注意事项】 (I)OpenFOAM软件一般在Linux系统环境中运行，因此需要先学习一些基本的 Linux操作命令。 (2)模拟计算时应根据宙诺数不同调整收敛误差的判据值 (3)划分网格时可在液滴运动路径范围和固体壁面处加密网格，以便兼颜计算效率 和计算精度的提高。 【思考与讨论】 (1)阐述计算流体力学的基本方法。 (2)液滴撞击壁面问题的控制方程如何用有限体积法离散化？ (3)根据所模拟的流体流动问题，如何在OpenFOAM软件库中选用合适的求解器 (4)试结合模拟实验结果分析VOF法追踪流体界面有哪些特点。 (⑤)试结合模拟实验结果分析液滴撞击壁面现象有哪些特征。 【参考文献】 [1]王先智.物理流体力学[M门.北京：清华大学出版社，2018 [2]吴小胜，黄晓鹏.计算流体力学基础[M们.北京：北京理工大学出版社，2021， [3] VVERSTEEG H K.MALALASEKERA W.An Introduction to Computationa Fluid Dynamics:The Finite Volume Method[M].北京：世界图书出版公司北 克公司.2010 [4]黄先北，郭墙.OpenFOAM从入门到精通[M们.北京：中国水利水电出版 牡.2021. [5]法德尔·穆卡莱德，卢卡·曼加尼，马尔万·达尔维什.计算流体力学中的有限 体积法OpenFOAM和Matlab高级导论[M们.李敏译.北京：中国水利水电出 版社，2020. [6]MOUKALLED F.LUCA MANGANI,MARWAN DARWISH.The Finite Volume Method in Computational Fluid Dynamies:An Advanced Introduction with Open- FOAM and Matlab[M].Switzerland:Springer International Publishing,2016. -347
— ３４７ — ２ 液滴撞击壁面运动特征的模拟分析 在最佳的网格划分方案条件下,改变液滴的大小和初速度以及液体的黏度等影响参 数,模拟计算液滴撞击壁面过程,并分析不同参数条件下液滴运动的特征. 二、设计性实验内容 在结合现有模拟计算条件的基础上,探究固体壁面粗糙度和表面形状、壁面润湿性对 液滴撞击壁面运动状态的影响,并分析液滴运动规律. 【注意事项】 (１)OpenFOAM 软件一般在 Linux系统环境中运行,因此需要先学习一些基本的 Linux操作命令. (２)模拟计算时应根据雷诺数不同调整收敛误差的判据值. (３)划分网格时可在液滴运动路径范围和固体壁面处加密网格,以便兼顾计算效率 和计算精度的提高. 【思考与讨论】 (１)阐述计算流体力学的基本方法. (２)液滴撞击壁面问题的控制方程如何用有限体积法离散化? (３)根据所模拟的流体流动问题,如何在 OpenFOAM 软件库中选用合适的求解器? (４)试结合模拟实验结果分析 VOF法追踪流体界面有哪些特点. (５)试结合模拟实验结果分析液滴撞击壁面现象有哪些特征. 【参考文献】 [１] 王先智．物理流体力学[M]．北京:清华大学出版社,２０１８． [２] 吴小胜,黄晓鹏．计算流体力学基础[M]．北京:北京理工大学出版社,２０２１． [３] VVERSTEEG H K,MALALASEKERA W．AnIntroductiontoComputational FluidDynamics:TheFiniteVolumeMethod[M]．北京:世界图书出版公司北 京公司,２０１０． [４] 黄 先 北,郭 嫱．OpenFOAM 从 入 门 到 精 通 [M]．北 京:中 国 水 利 水 电 出 版 社,２０２１． [５] 法德尔穆卡莱德,卢卡曼加尼,马尔万达尔维什．计算流体力学中的有限 体积法 OpenFOAM 和 Matlab高级导论[M]．李敏译．北京:中国水利水电出 版社,２０２０． [６] MOUKALLEDF,LUCA MANGANI,MARWAN DARWISH．TheFiniteVolume MethodinComputationalFluidDynamics:AnAdvancedIntroductionwithOpenＧ FOAM® andMatlab® [M]．Switzerland:SpringerInternationalPublishing,２０１６．
附 录 附录1中英文对照的常用物理量与单位 附表11中英文对照的常用物理量与单位 序号 物理量 单位 名称 符号 备注 The metre is the length of the path travelled by light in vacuu 长度 米 length uting a time interal of 1/200 792 458 of n vorond 17th CGPM (1983.Resolution 1:CR.97) 质量 千克 The kilogram is the unit of the maof the in- kg ernational prototype of the kilogram mass kilogram 3rd CGPM (1901.CR.70) The sccond is the duration of 9 192 631 770 periods of the radintion responding to the transition hetween the two hyperfine levels o 时间 秒 the ground state of the caesium 133 atom ti I3h CCPM (1867/68 Roio LCR 103) (Added by CIPM in 1997) The ampere is that constant current which.if maintained in two straight aralle nductorfneenghbeciruar 电流 安培 croonplaced 1 metre apart in vacuum.would produce A clectric current ampere between these conductors a force equal to 2 X 10-?newtons pe metre of length. 9th CGPM (1948.Resolution 2 and 7.CR.70) The kevin unit of thermodynamic temrerature.is the fraction 1/2731 13th CGPM (1967/68.Resolution 4CR.104) 热力学温度 开尔文 water having the isoto composition do K by the following am 0.00 155 76 mole of mole ofH.of 052 mole ofper mole of. (Added by CIPM in 205) 348
— ３４８ — 附录１ 中英文对照的常用物理量与单位 附表１Ｇ１ 中英文对照的常用物理量与单位 序号 物理量 单位 名称 单位 符号 备 注 １ 长度 length 米 meter m Themetreisthelengthofthepathtravelledbylightinvacuum duringatimeintervalof１/２９９７９２４５８ofasecond． １７thCGPM (１９８３,Resolution１;CR,９７) ２ 质量 mass 千克 kilogram kg Thekilogramistheunitofmass;itisequaltothemassoftheinＧ ternationalprototypeofthekilogram． ３rdCGPM (１９０１,CR,７０) ３ 时间 time 秒 second s Thesecondisthedurationof９１９２６３１７７０periodsoftheradiation correspondingtothetransitionbetweenthetwohyperfinelevelsof thegroundstateofthecaesium１３３atom． １３thCGPM (１９６７/６８,Resolution１;CR,１０３) Thisdefinitionreferstoacaesiumatomatrestatatemperatureof０K． (AddedbyCIPMin１９９７) ４ 电流 electriccurrent 安培 ampere A Theampereisthatconstantcurrentwhich,ifmaintainedintwo straightparallelconductorsofinfinitelength,ofnegligiblecircular crossＧsection,andplaced１ metreapartinvacuum,wouldproduce betweentheseconductorsaforceequalto２×１０－７ newtonsper metreoflength． ９thCGPM (１９４８,Resolution２and７;CR,７０) ５ 热力学温度 thermodynamic temperature 开尔文 kelvin K Thekelvin,unitofthermodynamictemperature,isthefraction１/２７３．１６ ofthethermodynamictemperatureofthetriplepointofwater． １３thCGPM (１９６７/６８,Resolution４;CR,１０４) ThisdefinitionreferstowaterhavingtheisotopiccompositiondeＧ fined exactly by the following amount of substance ratios: ０．０００１５５７６moleof２Hpermoleof１H,０．０００３７９９moleof１７O permoleof１６O,and０．００２００５２moleof１８Opermoleof１６O． (AddedbyCIPMin２００５)
附录 续表 序号 物理量 单位单位 名称符号 备注 ust be s 物质的量 may amount of substance particles n this definition.is understood that unbound aoms of carbor A 1980) The candela is the luminous intensity.in agiven direction.of 发光强度 坎德拉 souree that emits monochromatic radiation of frequeney5401 ed luminous intensity candela herz andtht hasaradiant intensity in that direction of1/683wat per steradian. 16th CGPM(1979.Resolution 3:CR.100) 平面角 弧度 rad m.m- ang 立体角 球面度 m2 m- solid angle steradian 顿率 赫芝 frequency hertz 力 牛顿 N k.ms- 压强 12 帕斯卡 PaN·m-2:kg·m-1·s-2 pressure pascal 能量 丝耳 13 N.mikg.m2.s-2 ioule 14 功 Js,kgm2s pow 电量 库仑 15 ·A electric charge coulomb 由压 伏特 16 J·C-1:kg·m2.8-3,A-1 electric potential yolt 电 VA1kgm21,A- ohm 电容 18 法拉第 F capacitance farad 349-
— ３４９ — 续表 序号 物理量 单位 名称 单位 符号 备 注 ６ 物质的量 amountofsubstance 摩尔 mole mol Themoleistheamountofsubstanceofasystem whichcontainsas manyelementaryentitiesasthereareatomsin０．０１２kilogramof carbon１２;itssymbolismol． Whenthemoleisused,theelementaryentitiesmustbespecified andmaybeatoms,molecules,ions,electrons,otherparticles,or specifiedgroupsofsuchparticles． １４thCGPM (１９７１,Resolution３;CR,７８) Inthisdefinition,itisunderstoodthatunboundatomsofcarbon １２,atrestandintheirgroundstate,arereferredto． (AddedbyCIPMin１９８０) ７ 发光强度 luminousintensity 坎德拉 candela cd Thecandelaistheluminousintensity,inagivendirection,ofa sourcethatemitsmonochromaticradiationoffrequency５４０×１０１２ hertzandthathasaradiantintensityinthatdirectionof１/６８３watt persteradian． １６thCGPM (１９７９,Resolution３;CR,１００) ８ 平面角 angle 弧度 radian rad mm－１ ９ 立体角 solidangle 球面度 steradian sr m２m－２ １０ 频率 frequency 赫兹 hertz Hz s－１ １１ 力 force 牛顿 newton N kgms－２ １２ 压强 pressure 帕斯卡 pascal Pa Nm－２;kgm－１s－２ １３ 能量 energy 焦耳 joule J Nm;kgm２s－２ １４ 功率 power 瓦特 watt W Js－１;kgm２s－３ １５ 电量 electriccharge 库仑 coulomb C sA １６ 电压 electricpotential 伏特 volt V JC－１;kgm２s－３A－１ １７ 电阻 resistance 欧姆 ohm Ω VA－１;kgm２s－３A－２ １８ 电容 capacitance 法拉第 farad F CV－１;kg－１m－２s４A２