场初始化与 setFields:工程设置与诊断验证

给出 setFieldsDict 的 defaultFieldValues 与 regions 两层结构、常用几何选择器的坐标写法与命令行选项,并用球体体积与网格分辨率估算初始化区域的单元数与体积分数,说明 topoSet 与表达式初始化工具的适用位置。

OpenFOAMsetFields结果诊断与可信度验证工程设置与参数选择体积分数静水压setFieldsDictboxToCell工程

场初始化与 setFields:工程设置与诊断验证

setFields 只做一件事:按几何选择器把指定单元上的场值改写为给定值,其余单元保留默认值。它不支持时间推进,也不理解物理,因此配置的核心是"选择器坐标写得对不对"和"网格能不能分辨出这个几何"。初场无法被"验证为正确",只能被验证为与设计意图一致。可核对的量有三类:区域体积与体积分数的积分值、压力初场与静水压关系的一致性、以及被选中单元的几何分布与设计几何的偏差。

常用选择器的坐标写法

boxToCell 用两个对角点定义轴对齐长方体,写法是 box (x0 y0 z0) (x1 y1 z1);,两点不要求有序,但必须是同一坐标系下的绝对坐标。sphereToCell 用 centre 与 radius。cylinderToCell 用 p1、p2、radius。rotatedBoxToCell 额外给 origin、i、j、k 三个方向向量。面向已有集合的写法用 cellSet 或 cellZoneToCell,只给集合名即可,几何细节交给 topoSet 处理。

坐标与网格必须同源。blockMeshDict 里的 vertices 与 setFieldsDict 里的坐标都是绝对坐标,单位为米;若几何由 snappyHexMesh 从 STL 生成而 STL 以毫米为单位,转换时通常靠 scale 条目,忘记缩放会让选择器落在网格之外。

用体积估算初始化规模

球体区域的体积为

Vsphere=43πr3V_{\text{sphere}} = \frac{4}{3}\pi r^{3}

半径 r=0.05 mr = 0.05\ \mathrm{m} 时 V=43×3.14159×1.25×10−4=5.236×10−4 m3V = \frac{4}{3}\times 3.14159\times 1.25\times 10^{-4} = 5.236\times 10^{-4}\ \mathrm{m^3}。单元尺度 5 mm 的立方单元体积为 1.25×10−7 m31.25\times 10^{-7}\ \mathrm{m^3},因此被选中的单元数约为

Nset≈VregionVcell=5.236×10−41.25×10−7≈4189N_{\text{set}} \approx \frac{V_{\text{region}}}{V_{\text{cell}}} = \frac{5.236\times 10^{-4}}{1.25\times 10^{-7}} \approx 4189

若计算域是边长 0.5 m 的立方体,总单元数为 0.125/1.25×10−7=1.0×1060.125/1.25\times 10^{-7} = 1.0\times 10^{6},被初始化单元占 0.42%0.42\%。这个比例有实际意义:体积分数场的初始占比直接决定自由液面演化到准稳态所需的时间,占比越小,瞬态越长,需要的物理时间也越长。若把单元尺度加密到 2 mm,单元体积降到 8×10−9 m38\times 10^{-9}\ \mathrm{m^3},同一球体覆盖约 65450 个单元,几何边界更贴合但内存与时间步代价同步上升。

defaultFieldValues 与 regions 的分工

setFieldsDict 只有两个顶层条目。defaultFieldValues 为整个场设定基线值,可以省略;省略时以现有 0/ 文件中的 internalField 作为基线,这一点决定了重复运行 setFields 的行为。regions 是一个列表,每个元素是一组"选择器 + fieldValues",按出现顺序依次施加,后面的区域覆盖前面的结果。

defaultFieldValues
(
    volScalarFieldValue alpha.water 0
);

regions
(
    boxToCell
    {
        box (0 0 0) (0.146 0.292 0.146);
        fieldValues
        (
            volScalarFieldValue alpha.water 1
            volVectorFieldValue U (0 0 0)
        );
    }

sphereToCell
    {
        centre  (0.25 0.25 0.25);
        radius  0.05;
        fieldValues
        (
            volScalarFieldValue T 350
        );
    }
);

区域顺序很重要:若把球体放在盒子之后,球内单元的 alpha.water 不会被重置,只会被改写 T;反过来则会同时改写两者。需要"从大区域中挖掉小块"时,应先设大区域,再用小块覆盖。

命令行选项与执行顺序

setFields -dict system/setFieldsDict
setFields -time 0 -noFunctionObjects
setFields -latestTime
grep -n "boxToCell\|sphereToCell" system/setFieldsDict

setFields 默认读取 system/setFieldsDict,用 -dict 可指向其他路径。-time 指定作用于哪个时间目录,默认是最新的;-latestTime 与之等价但语义更明确。-noFunctionObjects 在 controlDict 里挂了昂贵函数对象时能显著缩短启动时间。执行顺序上,setFields 必须在 blockMesh 之后、求解器之前,且对应场文件必须已存在于目标时间目录,否则会报找不到场。

topoSet 与 setFields 的分工是:topoSet 生成并持久化单元集合(cellSet、faceSet、cellZone),适合被多个工具反复引用;setFields 只做一次性赋值,不留下集合文件。若某个几何区域需要在边界条件、源项与后处理中多处引用,先用 topoSet 建集合,再用 cellZoneToCell 或 cellSet 在 setFieldsDict 中引用。表达式型初始化(例如按解析函数赋速度剖面)超出 setFields 的能力范围,需要 funkySetFields 一类的外部工具。

故障模式与判定试验

现象 根因 判定试验
区域内的场完全没有变化 选择器坐标落在网格外,选中零个单元 用 postProcess -func 'writeCellCentres' 导出单元中心,核对坐标范围
Cannot find field alpha.water 目标场未在 0/ 中创建 ls 0/ 确认场文件存在且类名正确
重复运行 setFields 结果越改越乱 省略 defaultFieldValues,基线取自上一次的结果 补上显式默认值,或先 cp -r 0.orig 0
球体边界呈阶梯状且体积偏差大 单元尺度相对几何尺寸过粗 按 NsetN_{\text{set}} 估算,单元边长应小于半径的 1/5
STL 几何与选择器对不上 STL 单位与网格单位不一致,缺少缩放 比较 STL 包围盒与 checkMesh 报告的 bounding box
初始平均体积分数是设计值的 1.5 倍 选择器几何大于设计几何 用 writeCellCentres 导出包围盒,与设计尺寸比较
计算起步阶段出现 0.1 m/s0.1\ \mathrm{m/s} 的虚假速度 压力初场与重力不平衡 用 h=Δp/(ρg)h = \Delta p/(\rho g) 反算液深,与设计液深比较
被选中单元边界呈明显锯齿且体积偏小 单元尺度相对几何过粗 加密到 2 mm 重跑,积分值应趋近理论体积
平均值随统计窗口起点强烈变化 初场瞬态被计入统计 用两个不同起点各算一次平均,差值应小于 2%
封闭算例中 αˉ\bar\alpha 持续漂移 相方程非守恒或边界与初场不兼容 检查入口出口通量是否为零,观察 εα\varepsilon_\alpha 是否随步数线性增长

与解析解的对照

setFields 之后的场是否可解,可以用一步无对流、无扩散的极限来检验:关掉对流项与扩散项、只留时间导数与源项,运行一步,场的变化应精确等于源项乘 Δt\Delta t。对 α\alpha 这类被代数约束的场,更实用的对照是体积守恒:在无进出流的封闭算例中,αˉ\bar\alpha 应随时间保持不变,漂移量

εα=αˉ(t)−αˉ(0)αˉ(0)\varepsilon_{\alpha} = \frac{\bar\alpha(t) - \bar\alpha(0)}{\bar\alpha(0)}

应保持在 10−610^{-6} 量级。若 0.5 s 内漂移到 10−210^{-2},说明求解器的相方程存在非守恒项,或者是 setFields 写出的非均匀场与边界条件不兼容。

用体积积分核对相分布

体积分数场的全域积分给出相体积,与设计几何的理论体积直接可比:

αˉ=∑iαiVi∑iVi\bar{\alpha} = \frac{\sum_i \alpha_i V_i}{\sum_i V_i}

设计值为半径 0.05 m0.05\ \mathrm{m} 的球,理论体积 43πr3=5.236×10−4 m3\frac{4}{3}\pi r^{3} = 5.236\times 10^{-4}\ \mathrm{m^3};计算域边长 0.5 m0.5\ \mathrm{m} 时总体积 0.125 m30.125\ \mathrm{m^3},理论平均体积分数为 5.236×10−4/0.125=4.189×10−35.236\times 10^{-4}/0.125 = 4.189\times 10^{-3},即 0.419%0.419\%。用 5 mm 单元离散后,被标记单元数与理论值 41894189 的偏差主要来自阶梯化,通常落在 ±5%\pm 5\% 以内。若积分得到的 αˉ\bar\alpha 是理论值的 1.5 倍,说明 boxToCell 或 sphereToCell 的坐标范围比设计几何大了一圈,属于配置错误而非离散误差。

取值的命令行方式是

postProcess -func 'fieldMinMax(alpha.water)' -time 0
postProcess -func 'writeCellCentres' -time 0
foamToVTK -time 0 -fields '(alpha.water U p)'

writeCellCentres 生成 C、Cx、Cy、Cz 五个场,把它们与体积分数一起导出,可以在外部脚本里直接算加权积分与几何包围盒。

静水压与初场的一致性

分层流或多相算例的压力初场若与重力不平衡,计算一开始就会产生虚假加速度。静水压关系为

h=Δpρgh = \frac{\Delta p}{\rho g}

已知水柱底部与顶部压差 Δp=2452.5 Pa\Delta p = 2452.5\ \mathrm{Pa}、ρ=1000 kg/m3\rho = 1000\ \mathrm{kg/m^3}、g=9.81 m/s2g = 9.81\ \mathrm{m/s^2},则液柱高度 h=2452.5/(1000×9.81)=0.25 mh = 2452.5/(1000\times 9.81) = 0.25\ \mathrm{m}。反过来,若设计液深 0.25 m 而 0/p 的压力差只有 245 Pa,那就是少了一个数量级,自由液面会在最初几十步内明显下沉或上浮。

判据是:初始时刻的最大虚假速度应远小于特征速度。若特征速度 Uref=1.0 m/sU_{\text{ref}} = 1.0\ \mathrm{m/s},静压不平衡导致的初始速度幅度应控制在 10−3 m/s10^{-3}\ \mathrm{m/s} 以下;若 log.run 第一步就报告最大速度达到 0.1 m/s0.1\ \mathrm{m/s},说明压力初场与重力项不匹配。

几何边界的逐单元核对

选择器与网格的偏差只有导出单元中心才看得清。做法是把 Cx、Cy、Cz 与 alpha.water 一起读入,筛出 α>0.5\alpha > 0.5 的单元,计算它们的坐标包围盒与球心距离分布。设计球心 (0.25,0.25,0.25)(0.25, 0.25, 0.25)、半径 0.05 m0.05\ \mathrm{m},则被选中单元的球心距应全部落在 0.05 m0.05\ \mathrm{m} 加半个单元对角线以内。5 mm 单元的对角线为 5×3≈8.66 mm5\times\sqrt{3}\approx 8.66\ \mathrm{mm},因此球心距上限约 0.0543 m0.0543\ \mathrm{m};若出现 0.07 m0.07\ \mathrm{m} 的单元,说明选择器中心写错或存在第二个区域意外覆盖。

初场对统计量的污染窗口

初场是人为的,它携带的瞬态扰动需要时间衰减。对自由液面坍塌这类算例,前段统计量的偏差可以用"与最终准稳态值的相对差"来量化。设液面高度在 t=0.05 st = 0.05\ \mathrm{s} 时为 0.300 m0.300\ \mathrm{m},在 t=0.40 st = 0.40\ \mathrm{s} 后稳定在 0.146 m0.146\ \mathrm{m},则前者的相对偏差为 (0.300−0.146)/0.146=1.055(0.300-0.146)/0.146 = 1.055,即超过 100%。只有当液面高度在连续 0.10 s 内的变化小于 2% 时,才可以开始做时间平均;把统计窗口从 t=0t = 0 起算会把 100% 量级的初场偏差摊进平均值,得到一个既非初场也非稳态的数。

判定试验:分别用"从 t=0t = 0 起算"和"从 t=0.30 st = 0.30\ \mathrm{s} 起算"两个窗口计算平均液面高度,两者相差应小于 2%。若差异达到 20%,说明统计窗口仍被初场污染。

参考资料

  1. OpenFOAM Foundation, OpenFOAM User Guide, v14, Section “Mesh and field manipulation: setFields”.
  2. OpenFOAM Foundation, OpenFOAM User Guide, v14, Section “The setFieldsDict dictionary”.
  3. C. J. Greenshields, Notes on Computational Fluid Dynamics: General Principles, CFD Direct, 2022.
  4. Rusche, H., Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions, PhD thesis, Imperial College London, 2002.
  5. H. K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, 2nd ed., Pearson, 2007.
  6. H. Jasak, Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows, PhD thesis, Imperial College London, 1996.
  7. OpenFOAM Foundation, OpenFOAM User Guide, v14, Section “Post-processing function objects”.
  8. J. U. Brackbill, D. B. Kothe, C. Zemach, “A continuum method for modeling surface tension”, Journal of Computational Physics, 100(2):335–354, 1992.
  9. G. B. Whitham, Linear and Nonlinear Waves, Wiley-Interscience, 1974.