拉格朗日求解器代码求助

[李东岳]

$$\begin{array}{}\text{(1)}& r\frac{{\mathrm{\partial }}^{2}r}{\mathrm{\partial }{t}^2}+1.5\frac{\mathrm{\partial }r}{\mathrm{\partial }t}=\frac{1}{\rho }(p_{g0}{\left(r_0/r\right)}^{3\gamma }-p{\rho }_{\mathrm{c}}-\frac{2\sigma }{r}-\frac{4\mu }{r}\frac{\mathrm{\partial }r}{\mathrm{\partial }t})+\frac{({\mathbf{U}}_{\mathrm{c}}-\mathbf{U}{)}^2}{4}\end{array}$$

我把这个方程简化成这个行不行

$$\begin{array}{}\text{(2)}& r\frac{{\mathrm{\partial }}^{2}r}{\mathrm{\partial }{t}^2}+1.5\frac{\mathrm{\partial }r}{\mathrm{\partial }t}=A+\frac{({\mathbf{U}}_{\mathrm{c}}-\mathbf{U}{)}^2}{4}\end{array}$$

$A=0$

[youhaoyu]

@李东岳 这个把右侧括号部分考虑成0的话可能不大妥吧，因为我是个动态过程，气泡的入口和出口的压力差相差了两个数量级，还是挺大的。要看这个气泡慢慢变大过程。
如果你把A考虑成0，方程无法体现气泡内外压力差和表面张力，和我的模型可能也有一些问题。
方程的解法其实也不是很大问题，主要是这些变量如何去申明定义太恶心了。东岳老师您这种就简化的模型不大适合我目前的课题。 最开始我也想过，直接使用最简单的理想方程PV=NRT这种去修改，但是我觉得这种方程使用起来不是特别好

[李东岳]

我只是提供一个框架，而不是实际把r计算出来。后续你要在框架上自己进行修改，比如把A改成真实的。

[youhaoyu]

@李东岳 我的模型大概是1米高的管道，下方10000大气压，上方600Pa（液体是钢液，密度7000），底部通入气泡这种模型，对于气泡和液体之间压强差这个就是我们主要考虑的因素。这种来说A感觉不能直接考虑为0，所以我一直在看如何直接植入这个完整的方程

[杨英狄]

粒子性质的定义不仅在KinematicParcel.C和KinematicParcel.H，在其他文件中也涉及相关代码。你可以正则搜索一下，看看那些文件涉及粒子性质的声明和定义。

[youhaoyu]

@杨英狄 谢谢，这几天我我看了看，是我有些地方没有定义完成，我看您之前也在往DPMFoam添加RP方程，添加成功了吗？我最近添加编译成功了，库和求解器都成功了，但是计算出来结果颗粒直径没有变化。（您好久没登陆了，我还在尝试联系您，没成功）

[杨英狄]

@youhaoyu 如果计算过程中发现粒径没有变化，可能是添加的RP方程计算程序没有被调用。参考原生的lagrange库，粒子速度的更新函数calcVelocity（），是在函数calc()里面被调用的。你可以看看你的RP方程函数是否在calc()里面调用了。

[youhaoyu]

@杨英狄 是的，我就是这么编写的，我写了个calcDiameter，但是就是没发现哪儿有问题，现在我在一点点改看看哪儿有问题。
但是现在麻烦来了，修改了代码之后，编译能成功，但是使用求解器时候会报错，什么浮点数异常什么的，我初始条件都没改，哎麻烦死了。能不能给个联系方式，想咨询下您问题，因为刚好看到您也修改过

[mingyang]

@youhaoyu 你好，能学习一下您的植入方法吗，我最近也被困在这了，期待你的回复，感谢

[youhaoyu]

@mingyang 您好可以加我qq 2269609762一起探讨，我也没做成功

[李东岳]

我在几年前做了个一拉格朗日粒子直径变化的代码。我现在已经记不起来太多了。能想起来的是：

1. 粒子直径跟流场的湍流动能耗散率有关系
2. 植入的是一个跟PBM有关的粒径变化方程，考虑了coalescense和破碎
3. 代码可以编译并且测试过可以体现粒径变化

可以把下面的MomentumCloud.C：替换到OpenFOAM原来的代码然后编译。搞明白之后看看思路然后适配到你们的算法上。

MomentumCloud.C：

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.
 
    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
#include "MomentumCloud.H"
#include "integrationScheme.H"
#include "interpolation.H"
#include "subCycleTime.H"
 
#include "InjectionModelList.H"
#include "DispersionModel.H"
#include "PatchInteractionModel.H"
#include "StochasticCollisionModel.H"
#include "SurfaceFilmModel.H"
 
// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::setModels()
{
    dispersionModel_.reset
    (
        DispersionModel<MomentumCloud<CloudType>>::New
        (
            subModelProperties_,
            *this
        ).ptr()
    );
 
    patchInteractionModel_.reset
    (
        PatchInteractionModel<MomentumCloud<CloudType>>::New
        (
            subModelProperties_,
            *this
        ).ptr()
    );
 
    stochasticCollisionModel_.reset
    (
        StochasticCollisionModel<MomentumCloud<CloudType>>::New
        (
            subModelProperties_,
            *this
        ).ptr()
    );
 
    surfaceFilmModel_.reset
    (
        SurfaceFilmModel<MomentumCloud<CloudType>>::New
        (
            subModelProperties_,
            *this
        ).ptr()
    );
 
    UIntegrator_.reset
    (
        integrationScheme::New
        (
            "U",
            solution_.integrationSchemes()
        ).ptr()
    );
}
 
 
template<class CloudType>
template<class TrackCloudType>
void Foam::MomentumCloud<CloudType>::solve
(
    TrackCloudType& cloud,
    typename parcelType::trackingData& td
)
{
    if (solution_.steadyState())
    {
        cloud.storeState();
 
        cloud.preEvolve();
 
        evolveCloud(cloud, td);
 
        if (solution_.coupled())
        {
            cloud.relaxSources(cloud.cloudCopy());
        }
    }
    else
    {
        /////////////////////////////////
        Info<< "Dyfluid: Evolve function in MomentumCloud.C" << nl;
       
        List<DynamicList<label>> PartLabelPre(U_.size());
        List<DynamicList<label>> PartLabelPost(U_.size());
 
        DynamicList<label> cellWithPartPre;
        DynamicList<label> cellWithPartPost;
 
        pNumber_.primitiveFieldRef() = 0.0;
        forAllIter(typename MomentumCloud<CloudType>, *this, iter)
        {
            parcelType& p = iter();
            pNumber_[p.cell()] += p.nParticle();
            PartLabelPre[p.cell()].append(p.origId());                                           
        }
 
        forAll(U_, cell)
        {
            if (PartLabelPre[cell].size() != 0)
            {
                cellWithPartPre.append(cell);
            }
        }
 
        //Info<< cellWithPartPre << nl;
        /////////////////////////////////
    	
        cloud.preEvolve();
 
        evolveCloud(cloud, td);
 
        if (solution_.coupled())
        {
            cloud.scaleSources();
        }
 
        /////////////////////////////////
        forAllIter(typename MomentumCloud<CloudType>, *this, iter)
        {
            parcelType& p = iter();
            PartLabelPost[p.cell()].append(p.origId());                                           
        }
 
        const volScalarField& epsi = U_.mesh().lookupObject<volScalarField>("epsilon.water");
        scalar C1 = 0.00481;
        scalar C2 = 0.08;
        scalar sigma = 0.072;
        scalar rhoc = 998.0;
        scalar D1 = 0.88;
        scalar D2 = 9e6;
        scalar muc = 1e-3;
        const scalar dt = this->mesh().time().deltaTValue();
        
        forAll(U_, cell)
        {
            // if there are particles in some cell
            if (PartLabelPost[cell].size() != 0)
            {
                cellWithPartPost.append(cell);
        
                scalar allVolume = 0.0;
                scalar allnParticle = 0.0;
                scalar allnParticleNew = 0.0;
                scalar diamAve = 0.0;
                //Info<< "Cell [" << cell << "] has " << PartLabelPost[cell].size() 
                //    << " particles" << nl;
 
                // loop all over particles
                forAllIter(typename MomentumCloud<CloudType>, *this, iter)
                {
                    parcelType& p = iter();
 
                    // loop particles label in this cell
                    forAll(PartLabelPost[cell], i)
                    {
                        // if this particle belongs to this cell
                        if (p.origId() == PartLabelPost[cell](i))
                        {
                            allVolume += p.nParticle()*M_PI/6.0*pow3(p.d());
                            allnParticle += p.nParticle();
                            diamAve = pow(allVolume/allnParticle*6.0/M_PI, 1.0/3.0);
                            //Info<< " Labels are " << PartLabelPost[cell](i)
                            //    << ", its nParticle is " << p.nParticle() 
                            //    << ", volume is " << allVolume << nl;
                        }
                    }
                }
 
                const scalar d = diamAve;
 
                // cell manipulation, calculate average diameter
                if (d > SMALL)
                {
                    const scalar epsilon = epsi[cell];
                    scalar gN = 
                        allnParticle*C1*pow(epsilon, 1.0/3.0)/pow(d, 2.0/3.0)
                       *exp(-C2*sigma/(rhoc*pow(epsilon, 2.0/3.0)*pow(d, 5.0/3.0)));
 
                    scalar aSqrN = 
                        sqr(allnParticle)*D1*pow(epsilon, 1.0/3.0)*4.0*sqrt(2.0)*pow(d, 7.0/6.0)
                       *exp(-D2*muc*rhoc*epsilon/sqr(sigma)*pow4(d)/16.0);
                
                    allnParticleNew = allnParticle + (gN - aSqrN)*dt;
                
                    diamAve = diamAve*pow(allnParticle/allnParticleNew, 1.0/3.0); 
                    //Info<< "gN = " << gN << nl;
                    //Info<< "aSqrN = " << aSqrN << nl;
                    //Info<< "Particle increased by " << allnParticleNew - allnParticle << nl;
                    //Info<< "diamAve is " << diamAve << nl;
                }
            
                scalar newVolume = 0.0;
 
                forAllIter(typename MomentumCloud<CloudType>, *this, iter)
                {
                    parcelType& p = iter();
 
                    forAll(PartLabelPost[cell], i)
                    {
                        if (p.origId() == PartLabelPost[cell](i))
                        {
                            p.nParticle() = allnParticleNew/PartLabelPost[cell].size();
                            p.d() = diamAve;
                            newVolume += p.nParticle()*M_PI/6.0*pow3(p.d());
                            //Info<< " after brecoa, volume is " << newVolume
                            //    << "p.nParticle is " << p.nParticle() << nl;
                        }
                    }
                }
            }
        }
        /////////////////////////////////
    }
 
    cloud.info();
 
    cloud.postEvolve();
 
    if (solution_.steadyState())
    {
        cloud.restoreState();
    }
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::buildCellOccupancy()
{
    if (cellOccupancyPtr_.empty())
    {
        cellOccupancyPtr_.reset
        (
            new List<DynamicList<parcelType*>>(this->mesh().nCells())
        );
    }
    else if (cellOccupancyPtr_().size() != this->mesh().nCells())
    {
        // If the size of the mesh has changed, reset the
        // cellOccupancy size
 
        cellOccupancyPtr_().setSize(this->mesh().nCells());
    }
 
    List<DynamicList<parcelType*>>& cellOccupancy = cellOccupancyPtr_();
 
    forAll(cellOccupancy, cO)
    {
        cellOccupancy[cO].clear();
    }
 
    forAllIter(typename MomentumCloud<CloudType>, *this, iter)
    {
        cellOccupancy[iter().cell()].append(&iter());
    }
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::updateCellOccupancy()
{
    // Only build the cellOccupancy if the pointer is set, i.e. it has
    // been requested before.
 
    if (cellOccupancyPtr_.valid())
    {
        buildCellOccupancy();
    }
}
 
 
template<class CloudType>
template<class TrackCloudType>
void Foam::MomentumCloud<CloudType>::evolveCloud
(
    TrackCloudType& cloud,
    typename parcelType::trackingData& td
)
{
    if (solution_.coupled())
    {
        cloud.resetSourceTerms();
    }
 
    if (solution_.transient())
    {
        label preInjectionSize = this->size();
 
        this->surfaceFilm().inject(cloud);
 
        // Update the cellOccupancy if the size of the cloud has changed
        // during the injection.
        if (preInjectionSize != this->size())
        {
            updateCellOccupancy();
            preInjectionSize = this->size();
        }
 
        injectors_.inject(cloud, td);
 
 
        // Assume that motion will update the cellOccupancy as necessary
        // before it is required.
        cloud.motion(cloud, td);
 
        stochasticCollision().update(td, solution_.trackTime());
    }
    else
    {
//        this->surfaceFilm().injectSteadyState(cloud);
 
        injectors_.injectSteadyState(cloud, td, solution_.trackTime());
 
        CloudType::move(cloud, td, solution_.trackTime());
    }
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::postEvolve()
{
    Info<< endl;
 
    if (debug)
    {
        this->writePositions();
    }
 
    this->dispersion().cacheFields(false);
 
    forces_.cacheFields(false);
 
    functions_.postEvolve();
 
    solution_.nextIter();
 
    if (this->db().time().writeTime())
    {
        outputProperties_.writeObject
        (
            IOstream::ASCII,
            IOstream::currentVersion,
            this->db().time().writeCompression(),
            true
        );
    }
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::cloudReset(MomentumCloud<CloudType>& c)
{
    CloudType::cloudReset(c);
 
    rndGen_ = c.rndGen_;
 
    forces_.transfer(c.forces_);
 
    functions_.transfer(c.functions_);
 
    injectors_.transfer(c.injectors_);
 
    dispersionModel_.reset(c.dispersionModel_.ptr());
    patchInteractionModel_.reset(c.patchInteractionModel_.ptr());
    stochasticCollisionModel_.reset(c.stochasticCollisionModel_.ptr());
    surfaceFilmModel_.reset(c.surfaceFilmModel_.ptr());
 
    UIntegrator_.reset(c.UIntegrator_.ptr());
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::MomentumCloud<CloudType>::MomentumCloud
(
    const word& cloudName,
    const volScalarField& rho,
    const volVectorField& U,
    const volScalarField& mu,
    const dimensionedVector& g,
    const bool readFields
)
:
    CloudType(cloudName, rho, U, mu, g, false),
    cloudCopyPtr_(nullptr),
    particleProperties_
    (
        IOobject
        (
            cloudName + "Properties",
            this->mesh().time().constant(),
            this->mesh(),
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    ),
    outputProperties_
    (
        IOobject
        (
            cloudName + "OutputProperties",
            this->mesh().time().timeName(),
            "uniform"/cloud::prefix/cloudName,
            this->mesh(),
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        )
    ),
    solution_(this->mesh(), particleProperties_.subDict("solution")),
    constProps_(particleProperties_),
    subModelProperties_
    (
        particleProperties_.subOrEmptyDict("subModels", true)
    ),
    rndGen_(0),
    cellOccupancyPtr_(),
    cellLengthScale_(mag(cbrt(this->mesh().V()))),
    pNumber_
    (
        IOobject
        (
            "pNumbers",
            this->mesh().time().timeName(),
            this->mesh(),
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        this->mesh(),
        dimensionedScalar(dimless, 0.0)
    ),
    rho_(rho),
    U_(U),
    mu_(mu),
    g_(g),
    pAmbient_(0.0),
    forces_
    (
        *this,
        this->mesh(),
        subModelProperties_.subOrEmptyDict
        (
            "particleForces",
            true
        ),
        true
    ),
    functions_
    (
        *this,
        particleProperties_.subOrEmptyDict("cloudFunctions"),
        true
    ),
    injectors_
    (
        subModelProperties_.subOrEmptyDict("injectionModels"),
        *this
    ),
    dispersionModel_(nullptr),
    patchInteractionModel_(nullptr),
    stochasticCollisionModel_(nullptr),
    surfaceFilmModel_(nullptr),
    UIntegrator_(nullptr),
    UTrans_
    (
        new volVectorField::Internal
        (
            IOobject
            (
                this->name() + ":UTrans",
                this->db().time().timeName(),
                this->db(),
                IOobject::READ_IF_PRESENT,
                IOobject::AUTO_WRITE
            ),
            this->mesh(),
            dimensionedVector(dimMass*dimVelocity, Zero)
        )
    ),
    UCoeff_
    (
        new volScalarField::Internal
        (
            IOobject
            (
                this->name() + ":UCoeff",
                this->db().time().timeName(),
                this->db(),
                IOobject::READ_IF_PRESENT,
                IOobject::AUTO_WRITE
            ),
            this->mesh(),
            dimensionedScalar( dimMass, 0)
        )
    )
{
    setModels();
 
    if (readFields)
    {
        parcelType::readFields(*this);
        this->deleteLostParticles();
    }
 
    if (solution_.resetSourcesOnStartup())
    {
        resetSourceTerms();
    }
}
 
 
template<class CloudType>
Foam::MomentumCloud<CloudType>::MomentumCloud
(
    const word& cloudName,
    const volScalarField& rho,
    const volVectorField& U,
    const dimensionedVector& g,
    const fluidThermo& carrierThermo,
    const bool readFields
)
:
    MomentumCloud(cloudName, rho, U, carrierThermo.mu(), g, readFields)
{}
 
 
template<class CloudType>
Foam::MomentumCloud<CloudType>::MomentumCloud
(
    MomentumCloud<CloudType>& c,
    const word& name
)
:
    CloudType(c, name),
    cloudCopyPtr_(nullptr),
    particleProperties_(c.particleProperties_),
    outputProperties_(c.outputProperties_),
    solution_(c.solution_),
    constProps_(c.constProps_),
    subModelProperties_(c.subModelProperties_),
    rndGen_(c.rndGen_),
    cellOccupancyPtr_(nullptr),
    cellLengthScale_(c.cellLengthScale_),
    pNumber_(c.pNumber_),
    rho_(c.rho_),
    U_(c.U_),
    mu_(c.mu_),
    g_(c.g_),
    pAmbient_(c.pAmbient_),
    forces_(c.forces_),
    functions_(c.functions_),
    injectors_(c.injectors_),
    dispersionModel_(c.dispersionModel_->clone()),
    patchInteractionModel_(c.patchInteractionModel_->clone()),
    stochasticCollisionModel_(c.stochasticCollisionModel_->clone()),
    surfaceFilmModel_(c.surfaceFilmModel_->clone()),
    UIntegrator_(c.UIntegrator_->clone()),
    UTrans_
    (
        new volVectorField::Internal
        (
            IOobject
            (
                this->name() + ":UTrans",
                this->db().time().timeName(),
                this->db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            c.UTrans_()
        )
    ),
    UCoeff_
    (
        new volScalarField::Internal
        (
            IOobject
            (
                name + ":UCoeff",
                this->db().time().timeName(),
                this->db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            c.UCoeff_()
        )
    )
{}
 
 
template<class CloudType>
Foam::MomentumCloud<CloudType>::MomentumCloud
(
    const fvMesh& mesh,
    const word& name,
    const MomentumCloud<CloudType>& c
)
:
    CloudType(mesh, name, c),
    cloudCopyPtr_(nullptr),
    particleProperties_
    (
        IOobject
        (
            name + "Properties",
            mesh.time().constant(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        )
    ),
    outputProperties_
    (
        IOobject
        (
            name + "OutputProperties",
            this->mesh().time().timeName(),
            "uniform"/cloud::prefix/name,
            this->mesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        )
    ),
    solution_(mesh),
    constProps_(),
    subModelProperties_(dictionary::null),
    rndGen_(0),
    cellOccupancyPtr_(nullptr),
    cellLengthScale_(c.cellLengthScale_),
    pNumber_(c.pNumber_),
    rho_(c.rho_),
    U_(c.U_),
    mu_(c.mu_),
    g_(c.g_),
    pAmbient_(c.pAmbient_),
    forces_(*this, mesh),
    functions_(*this),
    injectors_(*this),
    dispersionModel_(nullptr),
    patchInteractionModel_(nullptr),
    stochasticCollisionModel_(nullptr),
    surfaceFilmModel_(nullptr),
    UIntegrator_(nullptr),
    UTrans_(nullptr),
    UCoeff_(nullptr)
{}
 
 
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
template<class CloudType>
Foam::MomentumCloud<CloudType>::~MomentumCloud()
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::setParcelThermoProperties
(
    parcelType& parcel,
    const scalar lagrangianDt
)
{
    parcel.rho() = constProps_.rho0();
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::checkParcelProperties
(
    parcelType& parcel,
    const scalar lagrangianDt,
    const bool fullyDescribed
)
{
    const scalar carrierDt = this->mesh().time().deltaTValue();
    parcel.stepFraction() = (carrierDt - lagrangianDt)/carrierDt;
 
    if (parcel.typeId() == -1)
    {
        parcel.typeId() = constProps_.parcelTypeId();
    }
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::storeState()
{
    cloudCopyPtr_.reset
    (
        static_cast<MomentumCloud<CloudType>*>
        (
            clone(this->name() + "Copy").ptr()
        )
    );
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::restoreState()
{
    cloudReset(cloudCopyPtr_());
    cloudCopyPtr_.clear();
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::resetSourceTerms()
{
    UTransRef().field() = Zero;
    UCoeffRef().field() = 0.0;
}
 
 
template<class CloudType>
template<class Type>
void Foam::MomentumCloud<CloudType>::relax
(
    DimensionedField<Type, volMesh>& field,
    const DimensionedField<Type, volMesh>& field0,
    const word& name
) const
{
    const scalar coeff = solution_.relaxCoeff(name);
    field = field0 + coeff*(field - field0);
}
 
 
template<class CloudType>
template<class Type>
void Foam::MomentumCloud<CloudType>::scale
(
    DimensionedField<Type, volMesh>& field,
    const word& name
) const
{
    const scalar coeff = solution_.relaxCoeff(name);
    field *= coeff;
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::relaxSources
(
    const MomentumCloud<CloudType>& cloudOldTime
)
{
    this->relax(UTrans_(), cloudOldTime.UTrans_(), "U");
    this->relax(UCoeff_(), cloudOldTime.UCoeff_(), "U");
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::scaleSources()
{
    this->scale(UTrans_(), "U");
    this->scale(UCoeff_(), "U");
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::preEvolve()
{
    // force calculation of mesh dimensions - needed for parallel runs
    // with topology change due to lazy evaluation of valid mesh dimensions
    label nGeometricD = this->mesh().nGeometricD();
 
    Info<< "\nSolving " << nGeometricD << "-D cloud " << this->name() << endl;
 
    this->dispersion().cacheFields(true);
    forces_.cacheFields(true);
    updateCellOccupancy();
 
    pAmbient_ = constProps_.dict().template
        lookupOrDefault<scalar>("pAmbient", pAmbient_);
 
    functions_.preEvolve();
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::evolve()
{
    if (solution_.canEvolve())
    {
        typename parcelType::trackingData td(*this);
 
        solve(*this, td);
    }
}
 
 
template<class CloudType>
template<class TrackCloudType>
void Foam::MomentumCloud<CloudType>::motion
(
    TrackCloudType& cloud,
    typename parcelType::trackingData& td
)
{
    CloudType::move(cloud, td, solution_.trackTime());
 
    updateCellOccupancy();
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::patchData
(
    const parcelType& p,
    const polyPatch& pp,
    vector& nw,
    vector& Up
) const
{
    p.patchData(nw, Up);
    Up /= p.mesh().time().deltaTValue();
 
    // If this is a wall patch, then there may be a non-zero tangential velocity
    // component; the lid velocity in a lid-driven cavity case, for example. We
    // want the particle to interact with this velocity, so we look it up in the
    // velocity field and use it to set the wall-tangential component.
    if (!this->mesh().moving() && isA<wallPolyPatch>(pp))
    {
        const label patchi = pp.index();
        const label patchFacei = pp.whichFace(p.face());
 
        // We only want to use the boundary condition value only if it is set
        // by the boundary condition. If the boundary values are extrapolated
        // (e.g., slip conditions) then they represent the motion of the fluid
        // just inside the domain rather than that of the wall itself.
        if (U_.boundaryField()[patchi].fixesValue())
        {
            const vector Uw1 = U_.boundaryField()[patchi][patchFacei];
            const vector& Uw0 =
                U_.oldTime().boundaryField()[patchi][patchFacei];
 
            const scalar f = p.currentTimeFraction();
 
            const vector Uw = Uw0 + f*(Uw1 - Uw0);
 
            const tensor nnw = nw*nw;
 
            Up = (nnw & Up) + Uw - (nnw & Uw);
        }
    }
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::updateMesh()
{
    updateCellOccupancy();
    injectors_.updateMesh();
    cellLengthScale_ = mag(cbrt(this->mesh().V()));
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::autoMap(const mapPolyMesh& mapper)
{
    Cloud<parcelType>::autoMap(mapper);
 
    updateMesh();
}
 
 
template<class CloudType>
void Foam::MomentumCloud<CloudType>::info()
{
    vector linearMomentum = linearMomentumOfSystem();
    reduce(linearMomentum, sumOp<vector>());
 
    scalar linearKineticEnergy = linearKineticEnergyOfSystem();
    reduce(linearKineticEnergy, sumOp<scalar>());
 
    Info<< "Cloud: " << this->name() << nl
        << "    Current number of parcels       = "
        << returnReduce(this->size(), sumOp<label>()) << nl
        << "    Current mass in system          = "
        << returnReduce(massInSystem(), sumOp<scalar>()) << nl
        << "    Linear momentum                 = "
        << linearMomentum << nl
        << "   |Linear momentum|                = "
        << mag(linearMomentum) << nl
        << "    Linear kinetic energy           = "
        << linearKineticEnergy << nl;
 
    injectors_.info(Info);
    this->surfaceFilm().info(Info);
    this->patchInteraction().info(Info);
}
 
 
// ************************************************************************* //