CFD中文网

@青山 我做的不涉及三相，所以只有两相。然后你这个说只能在alphat里算感觉没道理啊，那如果是层流，这个输入都没有，就不允许闪蒸了嘛，应该不会把

@青山 既然壁面没有加热，那么就和一般的walls绝热条件应该是一样的设置。至于发生相变的相和之后的相，这个我就不太确定了，可能是phaseProperties里？

@李东岳 谢谢李老师，我把sqr函数看成sqrt了，现在我明白了，是能对上的

@李东岳 谢谢李老师，但是还是有一点问题。
Re的计算有：phasePair.C

Foam::tmp<Foam::volScalarField> Foam::phasePair::Re() const
{
    return magUr()*dispersed().d()/continuous().thermo().nu();
}

CdRe的计算有
SchillerNaumann.C

Foam::tmp<Foam::volScalarField> Foam::dragModels::SchillerNaumann::CdRe() const
{
    volScalarField Re(pair_.Re());

    return
        neg(Re - 1000)*24*(1.0 + 0.15*pow(Re, 0.687))
      + pos0(Re - 1000)*0.44*max(Re, residualRe_);
}

最后，K的计算有
dragModel.C

Foam::tmp<Foam::volScalarField> Foam::dragModel::Ki() const
{
    return
        0.75
       *CdRe()
       *swarmCorrection_->Cs()
       *pair_.continuous().rho()
       *pair_.continuous().thermo().nu()
       /sqr(pair_.dispersed().d());
}


Foam::tmp<Foam::volScalarField> Foam::dragModel::K() const
{
    return max(pair_.dispersed(), pair_.dispersed().residualAlpha())*Ki();
}

综上，Kf的表达式为：

比较您的文章中的

可以发现液滴直径的参与方式上不同，这个是什么情况，您了解吗

如题，检查曳力模型代码
位于dragModel.C

Foam::tmp<Foam::volScalarField> Foam::dragModel::Ki() const
{
    return
        0.75
       *CdRe()
       *swarmCorrection_->Cs()
       *pair_.continuous().rho()
       *pair_.continuous().thermo().nu()
       /sqr(pair_.dispersed().d());
}

两个问题，一是为什么有粘度，文献中基本是不考虑粘度的把，第二个问题是为什么液滴直径要开根号？
这里以李东岳老师的教程为例，文献形式一般为

检查Cd的计算，与文献中一致。检查到方程的部分，也是直接乘以的速度。
前后都没有问题，所以这里为什么会有粘度和被开了根号的直径呢？

@青山 他案例是这样的，但是实际上没有影响，界面内的相变也是正常算的。你初始的空泡份额和边界的别设为0，设置为一个很小的值，然后给液相热源，应该是可以发生闪蒸的

各位老师和学者好，我最近在尝试使用multiphaseEulerFoam来计算水的闪蒸过程。但是遇到了一些问题。对于纯水的两相流动（没有空气），应该是调用ThermalPhaseChangePhaseSystem

type thermalPhaseChangeMultiphaseSystem;

phases (steam water);

phaseChange on;

steam
{
    type            purePhaseModel;
    diameterModel   isothermal;
    isothermalCoeffs
    {
        d0              3e-3;
        p0              1e5;
    }

    residualAlpha   1e-3;
}

water
{
    type            purePhaseModel;
    diameterModel   constant;
    constantCoeffs
    {
        d               1e-4;
    }

    residualAlpha   1e-3;
}

在此基础上，后续的换热模型选用了RanzMarshall，但是我发现根本不会发生闪蒸（我在能量方程里对water相直接补充了热量，保证了温度上升）

in EEqns.H

forAll(fluid.anisothermalPhases(), anisothermalPhasei)
    {
        phaseModel& phase = fluid.anisothermalPhases()[anisothermalPhasei];

        const volScalarField& alpha = phase;
        tmp<volScalarField> tRho = phase.rho();
        const volScalarField& rho = tRho();
        tmp<volVectorField> tU = phase.U();
        const volVectorField& U = tU();
        
        volScalarField q = alpha*rho*(U&g)-alpha*rho*(U&g);

        if(phase.name() == "water")
        {
        q[93] = 5.0E8;
        q[94] = 5.0E8;
        q[95] = 5.0E8;
        q[96] = 5.0E8;
        q[97] = 5.0E8;
        }
···········

不会发生相变的原因确定于代码中，对于质量源项的计算有

ThermalPhaseChangePhaseSystem.C

        // Interfacial mass transfer update
        {
            volScalarField& dmdtf(*this->dmdtfs_[pair]);
            volScalarField& Tf(*this->Tfs_[pair]);

            const volScalarField Tsat(saturationModelIter()->Tsat(thermo1.p()));

            const volScalarField L
            (
                volatile_ != "none"
              ? this->Li(pair, volatile_, dmdtf, Tsat, latentHeatScheme::upwind)
              : this->L(pair, dmdtf, Tsat, latentHeatScheme::upwind)
            );

            volScalarField H1(this->heatTransferModels_[pair].first()->K(0));
            volScalarField H2(this->heatTransferModels_[pair].second()->K(0));

            volScalarField dmdtfNew((H1*(Tsat - T1) + H2*(Tsat - T2))/L);

其中，换热系数H1与2的计算
RanzMarshall.C

Foam::tmp<Foam::volScalarField>
        
        Foam::heatTransferModels::RanzMarshall::K(const scalar residualAlpha) const
        
        {
        
            volScalarField Nu(2 + 0.6*sqrt(pair_.Re())*cbrt(pair_.Pr()));
        
        
        
            return
        
                6
        
               *max(pair_.dispersed(), residualAlpha)
        
               *pair_.continuous().thermo().kappa()
        
               *Nu
        
               /sqr(pair_.dispersed().d());
        
        }

可以发现，当离散相为水蒸气（初始条件下为0），连续相为水时，输入的残差为0时，

            volScalarField H1(this->heatTransferModels_[pair].first()->K(0));
            volScalarField H2(this->heatTransferModels_[pair].second()->K(0));

这个H1和H2必然为0，因此计算所得到的相变量dmdtfNew也必然为0，无论液体的温度高过饱和温度多少。

请问有老师碰到过这种情况吗，任何点拨都感激不尽

@李东岳 我目前感觉是在压力方程中进行了一下处理，但是我不太确定。这一项和原始动量方程，以及和连续性方程，能量方程的迭代修正方式。在压力方程中有个这个

这里面的dmdt就是实际上的源项，但是我目前还没看太明白，他怎么迭代回去到能量和连续性方程，以及这一项怎么和速度联系的。 最近还有一些其他的事情，这块代码啃的比较慢

@Tens dgdt是可压缩性项，而非源项项目，对应在李东岳老师的说明里是这个部分

首先，我们假设一个速度恒定，不考虑压力和体积力，但考虑相变的两相流动过程，此时，k相的控制方程有：


其中的Sk表示k相的源项。我们知道，对于一个小网格内，Sk如果较大，而从气到液体或者液体到气体的相变过程势必会在网格上引起额外的对流通量，因为相变会使得体积发生变化。
在OpenFOAM中所使用的欧拉多相流求解器multiphaseEulerFoam中，相分数方程并未做额外处理，但是在能量方程中，却多了一项，具体的有：
in file AnisothermalPhaseModel.C

···
    const tmp<volScalarField> tcontErr(this->continuityError());
    const volScalarField& contErr(tcontErr());
···
    tmp<fvScalarMatrix> tEEqn
    (
        fvm::ddt(alpha, rho, he)
      + fvm::div(alphaRhoPhi, he)
      - fvm::Sp(contErr, he)

      + fvc::ddt(alpha, rho, K) + fvc::div(alphaRhoPhi, K)
      - contErr*K

      + this->divq(he)
     ==
        alpha*this->Qdot()
    );

其中，continuityError()为
in file MovingPhaseModel.C

template<class BasePhaseModel>
void Foam::MovingPhaseModel<BasePhaseModel>::correctContinuityError
(
    const volScalarField& source
)
{
    volScalarField& rho = this->thermoRef().rho();

    continuityError_ = fvc::ddt(*this, rho) + fvc::div(alphaRhoPhi_) - source;
}

不难看出，当相变/源项存在时，为了防止体积溢出的部分的能量影响网格内能量的计算，添加了一个额外的连续性修正项contErr来消除溢出项的影响（实际上动量方程中也添加了该项）。但是，这样的方式显然会破坏全场内的能量与质量守恒，两相条件蒸发，网格内产生的膨胀体积应当输运到其他网格中，但是这样的处理方式显然没有对该过程进行考虑。
因为膨胀引起网格面额外通量的分析和计算过程本身十分复杂，考虑起来非常困难，因此只能牺牲一些守恒性换来鲁棒性。但是无论在openfoamwiki中的理论手册、李东岳老师还是其他商业软件中都没有把这一项写在理论说明的文档里。想请教一下各位老师和学者，这样的处理方式是两相求解领域的共识吗？

@upc_ngh 有一个更简单的方法，我最近发现的，你索引的时候不要用

···lookupObject···

用

···lookupObjectRef···

用这个函数查找的结果可以直接修改并传值回去

@李东岳 谢谢李老师的回复，这个问题我最后发现是一个很简单的错误导致的。因为decompsePar没有把边界场进行拆分 ，导致各个处理器文件夹下的场定义缺失所以发生了错误

为了实现一些特殊需求，所以需要两个边界相互交互。但是修改后的边界条件仅在串行条件下可以运行。

下述代码为例：

//~in patch
label patchi = this->patch().boundaryMesh().findPatchID(patchName_);
const volScalarField& paraGet = this->db().objectRegistry::lookupObject<volScalarField>("paraA");//这一步可以正常运行
const scalarField& para = paraGet.boundaryField()[patchi];//这一步OK
Info << para << endl;//这一步在并行条件下会报错，也就是引用的对象在使用的时候会报错

在此基础上，我修改了代码，发现如果不是引用，直接取值，也会报错，因此可以断定是取边界值的函数在并行条件下有问题。

报错代码

//~in patch
label patchi = this->patch().boundaryMesh().findPatchID(patchName_);
const volScalarField& paraGet = this->db().objectRegistry::lookupObject<volScalarField>("paraA");//这一步可以正常运行
const scalarField para = paraGet.boundaryField()[patchi];////这一步在并行条件下会报错

最后，也尝试过从网格转移获取边界值，但是并行时也会在取值时报错。

//~in patch
label patchi = this->patch().boundaryMesh().findPatchID(patchName_);
const volScalarField& paraGet = this->db().objectRegistry::lookupObject<volScalarField>("paraA");//这一步可以正常运行
const fvMesh& mesh = patch().boundaryMesh().mesh();//正常运行
const fvPatchField<scalar>& para = mesh.boundary()[patchi].lookupPatchField<volScalarField, scalar>("paraA");//这一步可以正常运行
Info << para << endl;//这一步在并行条件下会报错，也就是引用的对象在使用的时候会报错

所以想向各位老师和同僚请教下怎么处理这种问题

@wangfei9088 好的，谢谢 那看来要用OpenFOAM做一些对流扩散问题的研究还有很多需要修正的地方

一般而言，组分方程中的拉普拉斯项，其广义扩散系数为D，即一种流体在另外一种流体中的扩散系数。然而，在OpenFOAM中的组分方程中，没有一个使用的是D，均使用的是alpha（热扩散系数）和mu（粘度，或者叫动量扩散系数）。

源代码中：
在chtMultiRegionFoam的组分方程中

            fvScalarMatrix YiEqn
            (
                fvm::ddt(rho, Yi)
              + mvConvection->fvmDiv(phi, Yi)
              - fvm::laplacian(turbulence.muEff(), Yi)//注意，这里调用了粘度
             ==
                reaction.R(Yi)
              + fvOptions(rho, Yi)
            );

在fireFoam的组分方程中

fvScalarMatrix YiEqn
            (
                fvm::ddt(rho, Yi)
              + mvConvection->fvmDiv(phi, Yi)
              - fvm::laplacian(turbulence->alphaEff(), Yi)//注意，此处为热扩散系数
              ==
                parcels.SYi(i, Yi)
              + surfaceFilm.Srho(i)
              + combustion->R(Yi)
              + fvOptions(rho, Yi)
            );

那么，为什么好像从来没有人提到或者询问这个问题（甚至可以说是很根本的问题），还是说我对代码的理解有误。望各位老师和同僚赐教

@upc_ngh 好的，非常感谢

@upc_ngh 从const转换到非const是怎么做的 ，我现在只能通过mesh.lookupObject<volScalarField>获得场值。

@upc_ngh 您好，我刚好最近也遇到这个问题，想要把src文件夹下的值传递到application中，请问您具体是怎么实现和解决的呢，望赐教，感激不尽

@bestucan 好的，谢谢老师，我试试

@bestucan 谢谢老师的回复，我这里关于第③点的说法可能不太准确。我只输入rho时（也就是Air源项为0），此时Air的质量是增加的，且水蒸气的总体质量也在变小，但是总体的质量分数是守恒为1的 所以还是比较奇怪吧，他这里怎么处理的代码我还没找到，好难挖，OpenFOAM这种嵌套