在输出时间步的同时，输出颗粒所受drag force文件

李东岳

fvVectorMatrix UcEqn
(
    fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + continuousPhaseTurbulence->divDevRhoReff(Uc)
 ==
    (1.0/rhoc)*cloudSU
);

后面(1.0/rhoc)*cloudSU是曳力项，你可以尝试把它写出来

zhangxc0223

@qingdong_wang @东岳谢谢一楼和东岳老师的回复。还想请问东岳老师，这个曳力项，应该写在哪个文件里输出？

我发现求解器(比如DPMFoam)里面的creatFields.H文件，输出的都是流体相的场文件，比如Uc, p, phic之类的，并没有颗粒相的场文件。我也查找了一楼所说的src/lagrangian/basic目录下的文件，也没发现控制输出颗粒场的文件。

所以还想麻烦您进一步指明，这个曳力项应该写在哪里？

cxzhmg

@zhangxc0223 老铁实现了吗？如何不通过createFields添加新变量来输出变量？或者怎么添加不需要read的变量来输出？我试了no_read还是提示没有对应文件

zhangxc0223

@cxzhmg 并没有琢磨明白
creatFields能添加的应该都是流体相的场，离散的颗粒相应该实现不了

李东岳

如果你要欧拉场，可以看cloudVolSuSp怎么输出的，如果你要拉格朗日场，去KinematicParcelIO.C文件中，添加你的dragForce:

参考下面都是要输出的信息


template<class ParcelType>
template<class CloudType>
void Foam::KinematicParcel<ParcelType>::readFields(CloudType& c)
{
    bool write = c.size();

    ParcelType::readFields(c);

    IOField<label> active
    (
        c.fieldIOobject("active", IOobject::MUST_READ),
        write
    );
    c.checkFieldIOobject(c, active);

    IOField<label> typeId
    (
        c.fieldIOobject("typeId", IOobject::MUST_READ),
        write
    );
    c.checkFieldIOobject(c, typeId);

    IOField<scalar> nParticle
    (
        c.fieldIOobject("nParticle", IOobject::MUST_READ),
        write
    );
    c.checkFieldIOobject(c, nParticle);

    IOField<scalar> d
    (
        c.fieldIOobject("d", IOobject::MUST_READ),
        write
    );
    c.checkFieldIOobject(c, d);

    IOField<scalar> dTarget
    (
        c.fieldIOobject("dTarget", IOobject::MUST_READ),
        write
    );
    c.checkFieldIOobject(c, dTarget);

    IOField<vector> U
    (
        c.fieldIOobject("U", IOobject::MUST_READ),
        write
    );
    c.checkFieldIOobject(c, U);

    IOField<scalar> rho
    (
        c.fieldIOobject("rho", IOobject::MUST_READ),
        write
    );
    c.checkFieldIOobject(c, rho);

    IOField<scalar> age
    (
        c.fieldIOobject("age", IOobject::MUST_READ),
        write
    );
    c.checkFieldIOobject(c, age);

    IOField<scalar> tTurb
    (
        c.fieldIOobject("tTurb", IOobject::MUST_READ),
        write
    );
    c.checkFieldIOobject(c, tTurb);

    IOField<vector> UTurb
    (
        c.fieldIOobject("UTurb", IOobject::MUST_READ),
        write
    );
    c.checkFieldIOobject(c, UTurb);

    label i = 0;

    forAllIter(typename CloudType, c, iter)
    {
        KinematicParcel<ParcelType>& p = iter();

        p.active_ = active[i];
        p.typeId_ = typeId[i];
        p.nParticle_ = nParticle[i];
        p.d_ = d[i];
        p.dTarget_ = dTarget[i];
        p.U_ = U[i];
        p.rho_ = rho[i];
        p.age_ = age[i];
        p.tTurb_ = tTurb[i];
        p.UTurb_ = UTurb[i];

        i++;
    }
}

马乔

这个题我会
第一步，在ParticleForce类中添加虚函数

virtual forceSuSp calcDragForce
(
      const typename CloudType::parcelType& p,
      const typename CloudType::parcelType::trackingData& td,
      const scalar dt,
      const scalar mass,
      const scalar Re,
      const scalar muc
 ) const
{
     return forceSuSp(Zero);
}

第二步，在ErgunWenYuDragForce类添加此函数，并返回

...
{
    return this->calcCoupled
   (
       p,
       td,
       dt,
       mass,
       Re,
       muc
   );

第三步，在ParticleForceList类中添加函数

forceSuSp calcDragForce
(
      const typename CloudType::parcelType& p,
      const typename CloudType::parcelType::trackingData& td,
      const scalar dt,
      const scalar mass,
      const scalar Re,
      const scalar muc
 ) const
{
    forceSuSp value(Zero);

    if (calcCoupled_)
    {
        forAll(*this, i)
        {
            value += this->operator[](i).calcDragForce(p, td, dt, mass, Re, muc);
        }
    }

    return value;
}

当然也可以这么写

...
{
    forceSuSp value(Zero);

    if (calcCoupled_)
    {
        forAll(*this, i)
        {
            if(isA<ErgunWenYuDragForce>(this->operator[](i)))
                 value = this->operator[](i).calcDragForce(p, td, dt, mass, Re, muc);
        }
    }

    return value;
}

这么写的话，前面两步就可以不用做了
第四步，在KinematicParcel类的calcVelocity函数中添加

template<class ParcelType>
template<class TrackCloudType>
const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
(
    TrackCloudType& cloud,
    trackingData& td,
    const scalar dt,
    const scalar Re,
    const scalar mu,
    const scalar mass,
    const vector& Su,
    vector& dUTrans,
    vector& drag,    //添加返回参数
    scalar& Spu
) const
{
    ...
    const forceSuSp FDrag = forces.calcDragForce(p, ttd, dt, mass, Re, mu);
    ...
    drag = -(FDrag.Sp()*td.Uc() + FDrag.Su());
    ...
}

在calc函数中添加

template<class ParcelType>
template<class TrackCloudType>
void Foam::KinematicParcel<ParcelType>::calc
(
    TrackCloudType& cloud,
    trackingData& td,
    const scalar dt
)
{
    ...
    vector drag = Zero;

    // Calculate new particle velocity
    this->U_ =
        calcVelocity(cloud, td, dt, Re, td.muc(), mass0, Su, dUTrans, darg, Spu);

    if (cloud.solution().coupled())
    {
        // Update momentum transfer
        cloud.UTrans()[this->cell()] += np0*dUTrans;

        cloud.dragForce()[this->cell()] += np0*drag;

        // Update momentum transfer coefficient
        cloud.UCoeff()[this->cell()] += np0*Spu;
    }

第五步，在KinematicCloud类中添加变量

autoPtr<volVectorField::Internal> drag_;

构造中添加

    drag_
    (
        new volVectorField::Internal
        (
            IOobject
            (
                this->name() + ":drag",
                this->db().time().timeName(),
                this->db(),
                IOobject::READ_IF_PRESENT,
                IOobject::AUTO_WRITE
            ),
            mesh_,
            dimensionedVector(dimForce, Zero)
        )
    ),

以及函数

template<class CloudType>
inline Foam::DimensionedField<Foam::vector, Foam::volMesh>&
Foam::KinematicCloud<CloudType>::drag()
{
    return *drag_;
}


template<class CloudType>
inline const Foam::DimensionedField<Foam::vector, Foam::volMesh>&
Foam::KinematicCloud<CloudType>::drag() const
{
    return *drag_;
}

第六步，主程序添加

        auto cloudDrag(kinematicCloud.drag());
        volVectorField cloudVolDrag
        (
            IOobject
            (
                "cloudVolDrag",
                runTime.timeName(),
                mesh
            ),
            mesh,
            dimensionedVector(cloudDrag.dimensions()/dimVolume, Zero),
            zeroGradientFvPatchVectorField::typeName
        );

        cloudVolDrag.primitiveFieldRef() = -cloudDrag/mesh.V();
        cloudVolSUSu.correctBoundaryConditions();

大概是这样，没编译过，新鲜的。里面的量纲和曳力应该还需要调整。

...

其实大部分情况下，一般都只会用到曳力，所以可以直接用kinematicCloud.SU(Uc)就行了，那这样的话，前面所有的都不用做了

李东岳

是的，OpenFOAM这面很多问题问起来2行就可以，但是要实现起来需要很多流程 @马乔感谢大佬分享！！！

马乔

修正下

value += this->operator[](i).calcCoupled(p, td, dt, mass, Re, muc);

cloud.drag()[this->cell()] += np0*drag;

cloudVolDrag.correctBoundaryConditions();

zhangxc0223

非常感谢大佬们的对话，受益匪浅啊

hotairballon1997

@马乔想借这里请问您一下，如果我只想输出颗粒的Uc的话（写入到每个时间步的/lagrangian/kinematicCloud中），应该怎么操作？

我代码只有入门水平，实在看不太明白您的一系列操作修改，望您有空能指教

CFD中文网

在输出时间步的同时，输出颗粒所受drag force文件