CFD中文网

余

目前我也在研究这个，不知道你还在不在研究这个问题了

余

模型背景

我正在使用 OpenFOAM 的 COxidationKineticDiffusionLimitedRate 模型：

配置如下：（Sb, C1, C2, E），温度约在 1250 K；

从源码可知 dmC 的计算包含 Ap·rhoc、RR·Tc·YO₂/WO₂、D₀·Rₖ/(D₀+Rₖ)，其组合单位似乎为 J/m，但源码最终把输出 dmC 当作 碳质量消耗 kg 来处理；

#include "COxidationKineticDiffusionLimitedRate.H"
#include "mathematicalConstants.H"
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::COxidationKineticDiffusionLimitedRate<CloudType>::
COxidationKineticDiffusionLimitedRate
(
    const dictionary& dict,
    CloudType& owner
)
:
    SurfaceReactionModel<CloudType>(dict, owner, typeName),
    Sb_(this->coeffDict().getScalar("Sb")),
    C1_(this->coeffDict().getScalar("C1")),
    C2_(this->coeffDict().getScalar("C2")),
    E_(this->coeffDict().getScalar("E")),
    CsLocalId_(-1),
    O2GlobalId_(owner.composition().carrierId("O2")),
    CO2GlobalId_(owner.composition().carrierId("CO2")),
    WC_(0.0),
    WO2_(0.0),
    HcCO2_(0.0)
{
    // Determine Cs ids
    label idSolid = owner.composition().idSolid();
    CsLocalId_ = owner.composition().localId(idSolid, "C");
 
    // Set local copies of thermo properties
    WO2_ = owner.thermo().carrier().W(O2GlobalId_);
    const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
    WC_ = WCO2 - WO2_;
 
    HcCO2_ = owner.thermo().carrier().Hc(CO2GlobalId_);
 
    const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
    const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
    Info<< "    C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
}
 
 
template<class CloudType>
Foam::COxidationKineticDiffusionLimitedRate<CloudType>::
COxidationKineticDiffusionLimitedRate
(
    const COxidationKineticDiffusionLimitedRate<CloudType>& srm
)
:
    SurfaceReactionModel<CloudType>(srm),
    Sb_(srm.Sb_),
    C1_(srm.C1_),
    C2_(srm.C2_),
    E_(srm.E_),
    CsLocalId_(srm.CsLocalId_),
    O2GlobalId_(srm.O2GlobalId_),
    CO2GlobalId_(srm.CO2GlobalId_),
    WC_(srm.WC_),
    WO2_(srm.WO2_),
    HcCO2_(srm.HcCO2_)
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
(
    const scalar dt,
    const scalar Re,
    const scalar nu,
    const label celli,
    const scalar d,
    const scalar T,
    const scalar Tc,
    const scalar pc,
    const scalar rhoc,
    const scalar mass,
    const scalarField& YGas,
    const scalarField& YLiquid,
    const scalarField& YSolid,
    const scalarField& YMixture,
    const scalar N,
    scalarField& dMassGas,
    scalarField& dMassLiquid,
    scalarField& dMassSolid,
    scalarField& dMassSRCarrier
) const
{
    // Fraction of remaining combustible material
    const label idSolid = CloudType::parcelType::SLD;
    const scalar fComb = YMixture[idSolid]*YSolid[CsLocalId_];
 
    // Surface combustion active combustible fraction is consumed
    if (fComb < SMALL)
    {
        return 0.0;
    }
 
    const SLGThermo& thermo = this->owner().thermo();
 
    // Local mass fraction of O2 in the carrier phase
    const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[celli];
 
    // Diffusion rate coefficient
    const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);
 
    // Kinetic rate
    const scalar Rk = C2_*exp(-E_/(RR*Tc));
 
    // Particle surface area
    const scalar Ap = constant::mathematical::pi*sqr(d);
 
    // Change in C mass [kg]
    scalar dmC = Ap*rhoc*RR*Tc*YO2/WO2_*D0*Rk/(D0 + Rk)*dt;
 
    // Limit mass transfer by availability of C
    dmC = min(mass*fComb, dmC);
 
    // Molar consumption
    const scalar dOmega = dmC/WC_;
 
    // Change in O2 mass [kg]
    const scalar dmO2 = dOmega*Sb_*WO2_;
 
    // Mass of newly created CO2 [kg]
    const scalar dmCO2 = dOmega*(WC_ + Sb_*WO2_);
 
    // Update local particle C mass
    dMassSolid[CsLocalId_] += dOmega*WC_;
 
    // Update carrier O2 and CO2 mass
    dMassSRCarrier[O2GlobalId_] -= dmO2;
    dMassSRCarrier[CO2GlobalId_] += dmCO2;
 
    const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
 
    // carrier sensible enthalpy exchange handled via change in mass
 
    // Heat of reaction [J]
    return dmC*HsC - dmCO2*HcCO2_;
}

COxidationKineticDiffusionLimitedRate 是 OpenFOAM 中用于描述 焦炭表面碳 (C(s)) 与氧气 (O₂) 反应生成 CO₂ 的动力学–扩散受限表面反应模型。
- 该模型只支持单一路径：C(s) + Sb·O₂ → CO₂；无法处理生成 CO 或 H₂ 的气化反应（如 char + CO₂ 或 char + H₂O）。(OpenFOAM)
社区用户曾指出，若使用 tutorial 中默认参数（Sb=1、C1=5e‑12、C2=0.002、E=7.9e7），在合理温度（如 Tc = 1250 K）下计算 Rk → 0，从而导致 dmC/dt → 0，看似模型不起作用。(CFD Online)

核心关注点

1. `dmC` 的单位与物理意义：

公式如下（简化版本）：

dmC = Ap * rhoc * RR * Tc * YO2 / WO2_ * (D0 * Rk / (D0 + Rk)) * dt;

各项单位：

Ap * rhoc → kg/m
RR*Tc*YO2/WO2_ → J/kg 或 J·kmol/kg²
D0 * Rk / (D0 + Rk) → 1/s
dt → s

组合后似乎是 J/m（能量密度单位），而非质量单位（kg），令人难以直接信服 dmC 实为 kg 碳质量损耗。

社区质疑这一点，但也指出源码后续使用 dmC/WC 计算摩尔变化，再映射为气相 O₂ 与 CO₂ 质量，最终模型确实是按质量守恒执行。(CFD Online)

2. 模型常参数来源与适用性

默认参数值（tutorial 案例中提供）在合理温度下几乎不启动：
- 当 Tc = 1250 K，Rk ≈ 0，导致 dmC/dt = 0，即不发生碳消耗；
- 需要极高温度（如 Tc → 50000 K）才会有所反应，这显然不现实。(CFD Online)
参数来源据称为 Baum & Street (1971) 研究，但社区提问第二点是：
- 确切常数来源不明；
- 模型中使用 Tc（载体温度）而非 Tp（碳颗粒表面温度），与文献描述不一致。(CFD Online)

总结

我对以下三点非常困惑，想向大家询问：

单位体系是否正确？
虽然中间结果为 J/m，但 dmC 最终确实代表 kg 碳损耗？有哪些地方负责单位换算以确保质量守恒？
为何采用载体温度 Tc 而非颗粒表面温度 Tp？
文献模型多数以 Tp（carbon particle temperature）为基础，但代码却用 Tc，是否存在误用或简化？
默认参数值是否合理
常数如 C2=0.002, E=7.9e7 是否适用？

余

补充一下：现行dev版本的OpenFOAM只有单向耦合和双向耦合，关于四向耦合即颗粒颗粒之间的pair collision似乎没有找到对应的Lagrangian Models，等待后续更新吧！

余

@李东岳李老师，我其实没太懂新老版本两个版本之间的粒子追踪算法的区别，尤其是到NCC界面或者AMI界面时会有什么质的区别

余

@李东岳感谢李老师，我准备转向org版本进行研究，采用全新求解器重新实操一下

余

自问自达一下，我找到了！在dev版本中确实已经集成发布了,

目前，field-based Lagrangian库实现了一些与固体、恒定密度粒子动力学相关的云。粒子尚未具有热力学特性，热传递、质量传递和相变过程尚未实现。开发团队正在努力添加这些功能。

已经实现的新功能包括：

二阶属性演化和二阶（抛物线）运动
注入粒子属性的泛化初始化，包括分布和属性之间的函数关系
内部和外部校正迭代
连续更新的平均属性
基于单元点平滑的载体源累积

兼容性和过渡

传统的particle-based Lagrangian库将继续保留，直到新的field-based库被认为功能足够完善。两个库通过大小写区分：
基于粒子的库使用小写的"lagrangian"命名对象、文件和目录
基于字段的库使用大写的"Lagrangian"

可用的教程和测试案例

目前有有限的教程和测试案例可用：

$FOAM_TUTORIALS/incompressibleFluid/TJunction教程已更新，使用基于新field-based Lagrangian库的函数对象。
在$WM_PROJECT_DIR/test/Lagrangian中添加了一些测试，以便方便地验证代码的行为。

余

我最近研究OpenFOAM文档时，发现了一个由Will Bainbridge于2025年1月7日编写的关于全新"field-based Lagrangian"系统的技术文档。这似乎是对传统particle-based Lagrangian库的彻底重构和革命性改进。通过详细阅读，我注意到这个新系统引入了大量令人兴奋的变化和功能增强。
新系统的主要创新点包括但不限于：
1.全新的场域-粒子耦合架构：从传统的particle-based转向更灵活、更一致的field-based架构。
2.改进的源项处理：引入了精细的SpSchemes配置，允许更精确地控制源项的隐式/显式处理
3.扩展的数值方案：包括新的插值格式、离散化方法和限制器。
……
……

这些创新看起来非常有前景且对我很有用，但我在最新的OpenFOAM-dev版本的tutorials目录中未能找到任何使用这一新系统的示例案例。我想问的是：

1.这个field-based Lagrangian系统目前处于什么开发阶段？它是否已经在某个特定的OpenFOAM版本或分支中实现？或者还只是计划中的功能？

2.是否有任何可用的教程案例、示例代码或测试案例来演示这一系统的完整功能？特别是展示新的SpSchemes、云属性定义、物理模型和求解器设置的例子。

有人知道相关信息或者有关指南文章可以参考吗？

余

我在使用MPPICDyMFoam模拟颗粒两相流时遇到了一个问题。计算到中间某一步时，颗粒被添加到流场后，程序似乎陷入了无限循环状态，无法继续计算到下一个时间步。（版本号of2406,采用并行32核）

具体表现为：程序成功完成了PIMPLE迭代计算（包括压力方程和速度方程的求解，残差也正常收敛），并显示了连续性误差和执行时间。随后程序进入下一个时间步（Time = 0.206549），正确创建了AMI接口，并开始求解颗粒云（Solving 3-D cloud kinematicCloud）。颗粒注入器也成功添加了29个新的颗粒包（parcels），系统中当前共有57958个颗粒包。

在这一步后，程序就停滞不前，没有任何错误提示，也没有继续进行下一步计算。日志显示的关键参数如下：

最大Courant数约为20
系统中的最大颗粒体积分数为0.153386
系统中的当前质量为0.0159749
平均每个parcel中的实际颗粒数为227453
我怀疑这可能与以下因素有关：

颗粒浓度在某些区域过高导致MPPIC模型难以收敛？
Courant数过大引起数值不稳定？
AMI接口权重不完全匹配（日志中显示有些权重不是1）？
颗粒-颗粒或颗粒-壁面相互作用的处理问题？
颗粒丢失没有被捕捉到，陷入无限循环？

我已尝试观察计算过程中的内存和CPU使用情况，没有发现明显资源耗尽的现象。请问有谁遇到过类似问题或有解决方案建议？是否需要调整特定的模型参数或时间步长控制策略？

Courant Number mean: 0.0772217 max: 20.0091
deltaT = 4.24916e-05
Time = 0.206506

AMI: Creating AMI for source:AMI1 and target:AMI1_slave
AMI: Patch source faces: 240
AMI: Patch target faces: 238
AMI: Patch source sum(weights) min:0.978615 max:1 average:0.998005
AMI: Patch target sum(weights) min:0.934139 max:1 average:0.992771
AMI: Creating AMI for source:AMI2 and target:AMI2_slave
AMI: Patch source faces: 10508
AMI: Patch target faces: 9877
AMI: Patch source sum(weights) min:0.99992 max:1.00033 average:1.00004
AMI: Patch target sum(weights) min:0.999899 max:1.00024 average:1.00004
AMI: Creating AMI for source:AMI3 and target:AMI3_slave
AMI: Patch source faces: 3438
AMI: Patch target faces: 3416
AMI: Patch source sum(weights) min:0.981432 max:1 average:0.999092
AMI: Patch target sum(weights) min:0.981187 max:1 average:0.999103

Solving3-D cloud kinematicCloud

Cloud: kinematicCloud injector: model1
Added 29 new parcels

Cloud: kinematicCloud
Current number of parcels = 57958
Current mass in system = 0.0159749
Linear momentum = (0.000108739 -0.000681819 -0.133831)
|Linear momentum| = 0.133833
Linear kinetic energy = 1.59464
Average particle per parcel = 227453
Injector model1:
- parcels added = 59501
- mass introduced = 0.0163974
Surface film:
- parcels absorbed = 0
- mass absorbed = 0
- parcels ejected = 0
Parcel fate: system (number, mass)
- escape = 1686, 0.000461168
Parcel fate: patch (walls|rotatewall|wall-blade|inlet-airbig|inlet-airsmall) (number, mass)
- escape = 0, 0
- stick = 0, 0
Parcel fate: patch inlet-par (number, mass)
- escape = 0, 0
- stick = 0, 0
Parcel fate: patch outletcoarse (number, mass)
- escape = 0, 0
- stick = 0, 0
Parcel fate: patch outletfine (number, mass)
- escape = 1543, 0.000422535
- stick = 0, 0
Parcel fate: patch (AMI1|AMI1_slave|AMI2|AMI2_slave|AMI3|AMI3_slave) (number, mass)
- escape = 0, 0
- stick = 0, 0
Min cell volume fraction = 0
Max cell volume fraction = 0.153386
Min dense number of parcels = 1

PIMPLE: iteration 1
smoothSolver: Solving for U.airx, Initial residual = 0.000518979, Final residual = 1.41353e-05, No Iterations 1
smoothSolver: Solving for U.airy, Initial residual = 0.000566608, Final residual = 1.61927e-05, No Iterations 1
smoothSolver: Solving for U.airz, Initial residual = 0.000457069, Final residual = 1.44488e-05, No Iterations 1
GAMG: Solving for p, Initial residual = 0.165607, Final residual = 0.00108703, No Iterations 2
GAMG: Solving for p, Initial residual = 0.011492, Final residual = 0.000100387, No Iterations 4
GAMG: Solving for p, Initial residual = 0.00196178, Final residual = 1.24262e-05, No Iterations 5
time step continuity errors : sum local = 2.99461e-08, global = -3.17012e-10, cumulative = 4.8858e-07
GAMG: Solving for p, Initial residual = 0.110414, Final residual = 0.000698195, No Iterations 2
GAMG: Solving for p, Initial residual = 0.00769959, Final residual = 5.49637e-05, No Iterations 5
GAMG: Solving for p, Initial residual = 0.00131523, Final residual = 7.52728e-08, No Iterations 15
time step continuity errors : sum local = 1.22513e-08, global = -1.95327e-10, cumulative = 4.88385e-07
PIMPLE: iteration 2
smoothSolver: Solving for U.airx, Initial residual = 0.000406817, Final residual = 3.40436e-06, No Iterations 3
smoothSolver: Solving for U.airy, Initial residual = 0.000448481, Final residual = 4.08898e-06, No Iterations 3
smoothSolver: Solving for U.airz, Initial residual = 0.000370139, Final residual = 4.3305e-06, No Iterations 3
GAMG: Solving for p, Initial residual = 0.0605394, Final residual = 0.000293902, No Iterations 3
GAMG: Solving for p, Initial residual = 0.00610706, Final residual = 3.71595e-05, No Iterations 4
GAMG: Solving for p, Initial residual = 0.000910397, Final residual = 6.90977e-06, No Iterations 4
time step continuity errors : sum local = 2.28203e-08, global = -2.98262e-10, cumulative = 4.88087e-07
GAMG: Solving for p, Initial residual = 0.00789404, Final residual = 6.29906e-05, No Iterations 3
GAMG: Solving for p, Initial residual = 0.00173749, Final residual = 1.36789e-05, No Iterations 3
GAMG: Solving for p, Initial residual = 0.000208345, Final residual = 7.866e-08, No Iterations 10
time step continuity errors : sum local = 1.2266e-08, global = -1.9661e-10, cumulative = 4.8789e-07
smoothSolver: Solving for k.air, Initial residual = 0.00056457, Final residual = 3.34182e-06, No Iterations 3
ExecutionTime = 2340.61 s ClockTime = 2357 s

Courant Number mean: 0.077187 max: 20.0223
deltaT = 4.24444e-05
Time = 0.206549

AMI: Creating AMI for source:AMI1 and target:AMI1_slave
AMI: Patch source faces: 240
AMI: Patch target faces: 238
AMI: Patch source sum(weights) min:0.978606 max:1 average:0.998006
AMI: Patch target sum(weights) min:0.93414 max:1 average:0.992772
AMI: Creating AMI for source:AMI2 and target:AMI2_slave
AMI: Patch source faces: 10508
AMI: Patch target faces: 9877
AMI: Patch source sum(weights) min:0.999918 max:1.00033 average:1.00004
AMI: Patch target sum(weights) min:0.999928 max:1.00024 average:1.00004
AMI: Creating AMI for source:AMI3 and target:AMI3_slave
AMI: Patch source faces: 3438
AMI: Patch target faces: 3416
AMI: Patch source sum(weights) min:0.981457 max:1 average:0.99909
AMI: Patch target sum(weights) min:0.981189 max:1 average:0.999102

Solving3-D cloud kinematicCloud

Cloud: kinematicCloud injector: model1
Added 29 new parcels

余

我也遇到一模一样的情况，请问解决了吗

余

目前用Gerris和Basilisk模拟气泡和雾化研究的人越来越多，可是中国网站上的资料太少了…………

余

http://basilisk.fr/src/examples/bubble.c
希望对你有用

余

请问大家，明年想去国外联培一短时间，目前从事油气多相浮射流的研究，很想找一个有做PBM传统的课题组，求介绍，谢谢！

余

@Mikasa 同志你会吗，能指点迷津吗，谢谢~~

余

气泡被水带着一直往下跑，而不会上升…………

余

余

VOF模型中使用上下周期性边界时，应该加一个反向的重力加速度，不然水就自由落体了，请问这个体力（bodyforce）的源项该怎么加，谢谢!，（不想用udf文件)

替代文字

余

@yhdthu 我的算例很简单

blocks
(
    hex (0 1 2 3 4 5 6 7) (160 320 160) simpleGrading (1 1 1)
);

余

@yhdthu 在边界处局部加密吗

CFD中文网

余正东

帖子

模型背景

核心关注点

1. `dmC` 的单位与物理意义：

2. 模型常参数来源与适用性

总结

已经实现的新功能包括：

兼容性和过渡

可用的教程和测试案例

CFD中文网

余正东

帖子

模型背景

核心关注点

1. dmC 的单位与物理意义：

2. 模型常参数来源与适用性

总结

已经实现的新功能包括：

兼容性和过渡

可用的教程和测试案例

1. `dmC` 的单位与物理意义：