C++ TrackData::cloud方法代码示例

bool Foam::CollidingParcel<ParcelType>::move
(
    TrackData& td,
    const scalar trackTime
)
{
    typename TrackData::cloudType::parcelType& p =
        static_cast<typename TrackData::cloudType::parcelType&>(*this);

    td.switchProcessor = false;
    td.keepParticle = true;

    const polyMesh& mesh = td.cloud().pMesh();
    const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
    const scalarField& V = mesh.cellVolumes();

    switch (td.part())
    {
        case TrackData::tpVelocityHalfStep:
        {
            // First and last leapfrog velocity adjust part, required
            // before and after tracking and force calculation

            p.U() += 0.5*trackTime*p.f()/p.mass();

            p.angularMomentum() += 0.5*trackTime*p.torque();

            break;
        }

        case TrackData::tpLinearTrack:
        {
            scalar tEnd = (1.0 - p.stepFraction())*trackTime;
            const scalar dtMax = tEnd;

            while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
            {
                // Apply correction to position for reduced-D cases
                meshTools::constrainToMeshCentre(mesh, p.position());

                // Set the Lagrangian time-step
                scalar dt = min(dtMax, tEnd);

                // Remember which cell the parcel is in since this
                // will change if a face is hit
                const label cellI = p.cell();

                const scalar magU = mag(p.U());
                if (p.active() && magU > ROOTVSMALL)
                {
                    const scalar d = dt*magU;
                    const scalar maxCo = td.cloud().solution().maxCo();
                    const scalar dCorr = min(d, maxCo*cbrt(V[cellI]));
                    dt *=
                        dCorr/d
                       *p.trackToFace(p.position() + dCorr*p.U()/magU, td);
                }

                tEnd -= dt;
                p.stepFraction() = 1.0 - tEnd/trackTime;

                // Avoid problems with extremely small timesteps
                if (dt > ROOTVSMALL)
                {
                    // Update cell based properties
                    p.setCellValues(td, dt, cellI);

                    if (td.cloud().solution().cellValueSourceCorrection())
                    {
                        p.cellValueSourceCorrection(td, dt, cellI);
                    }

                    p.calc(td, dt, cellI);
                }

                if (p.onBoundary() && td.keepParticle)
                {
                    if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
                    {
                        td.switchProcessor = true;
                    }
                }

                p.age() += dt;
            }

            break;
        }

        case TrackData::tpRotationalTrack:
        {
            notImplemented("TrackData::tpRotationalTrack");

            break;
        }

        default:
        {
            FatalErrorIn
            (
//.........这里部分代码省略.........

void Foam::DsmcParcel<ParcelType>::hitWallPatch
(
    const wallPolyPatch& wpp,
    TrackData& td,
    const tetIndices& tetIs
)
{
    label wppIndex = wpp.index();

    label wppLocalFace = wpp.whichFace(this->face());

    const scalar fA = mag(wpp.faceAreas()[wppLocalFace]);

    const scalar deltaT = td.cloud().pMesh().time().deltaTValue();

    const constantProperties& constProps(td.cloud().constProps(typeId_));

    scalar m = constProps.mass();

    vector nw = wpp.faceAreas()[wppLocalFace];
    nw /= mag(nw);

    scalar U_dot_nw = U_ & nw;

    vector Ut = U_ - U_dot_nw*nw;

    scalar invMagUnfA = 1/max(mag(U_dot_nw)*fA, VSMALL);

    td.cloud().rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA;

    td.cloud().rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA;

    td.cloud().linearKEBF()[wppIndex][wppLocalFace] +=
        0.5*m*(U_ & U_)*invMagUnfA;

    td.cloud().internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA;

    td.cloud().iDofBF()[wppIndex][wppLocalFace] +=
        constProps.internalDegreesOfFreedom()*invMagUnfA;

    td.cloud().momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA;

    // pre-interaction energy
    scalar preIE = 0.5*m*(U_ & U_) + Ei_;

    // pre-interaction momentum
    vector preIMom = m*U_;

    td.cloud().wallInteraction().correct
    (
        static_cast<DsmcParcel<ParcelType> &>(*this),
        wpp
    );

    U_dot_nw = U_ & nw;

    Ut = U_ - U_dot_nw*nw;

    invMagUnfA = 1/max(mag(U_dot_nw)*fA, VSMALL);

    td.cloud().rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA;

    td.cloud().rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA;

    td.cloud().linearKEBF()[wppIndex][wppLocalFace] +=
        0.5*m*(U_ & U_)*invMagUnfA;

    td.cloud().internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA;

    td.cloud().iDofBF()[wppIndex][wppLocalFace] +=
        constProps.internalDegreesOfFreedom()*invMagUnfA;

    td.cloud().momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA;

    // post-interaction energy
    scalar postIE = 0.5*m*(U_ & U_) + Ei_;

    // post-interaction momentum
    vector postIMom = m*U_;

    scalar deltaQ = td.cloud().nParticle()*(preIE - postIE)/(deltaT*fA);

    vector deltaFD = td.cloud().nParticle()*(preIMom - postIMom)/(deltaT*fA);

    td.cloud().qBF()[wppIndex][wppLocalFace] += deltaQ;

    td.cloud().fDBF()[wppIndex][wppLocalFace] += deltaFD;

}

void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
{
    if (!this->active())
    {
        return;
    }

    // Retrieve the film model from the owner database
    const regionModels::surfaceFilmModels::surfaceFilmModel& filmModel =
        this->owner().mesh().time().objectRegistry::template lookupObject
        <regionModels::surfaceFilmModels::surfaceFilmModel>
    (
        "surfaceFilmProperties"
    );

    if (!filmModel.active())
    {
        return;
    }

    const labelList& filmPatches = filmModel.intCoupledPatchIDs();
    const labelList& primaryPatches = filmModel.primaryPatchIDs();

    const fvMesh& mesh = this->owner().mesh();
    const polyBoundaryMesh& pbm = mesh.boundaryMesh();

    forAll(filmPatches, i)
    {
        const label filmPatchi = filmPatches[i];
        const label primaryPatchi = primaryPatches[i];

        const labelList& injectorCellsPatch = pbm[primaryPatchi].faceCells();

        cacheFilmFields(filmPatchi, primaryPatchi, filmModel);

        const vectorField& Cf = mesh.C().boundaryField()[primaryPatchi];
        const vectorField& Sf = mesh.Sf().boundaryField()[primaryPatchi];
        const scalarField& magSf = mesh.magSf().boundaryField()[primaryPatchi];

        forAll(injectorCellsPatch, j)
        {
            if (diameterParcelPatch_[j] > 0)
            {
                const label celli = injectorCellsPatch[j];

                // The position could bein any tet of the decomposed cell,
                // so arbitrarily choose the first face of the cell as the
                // tetFace and the first point on the face after the base
                // point as the tetPt.  The tracking will pick the cell
                // consistent with the motion in the first tracking step.
                const label tetFacei = this->owner().mesh().cells()[celli][0];
                const label tetPti = 1;

//                const point& pos = this->owner().mesh().C()[celli];

                const scalar offset =
                    max
                    (
                        diameterParcelPatch_[j],
                        deltaFilmPatch_[primaryPatchi][j]
                    );
                const point pos = Cf[j] - 1.1*offset*Sf[j]/magSf[j];

                // Create a new parcel
                parcelType* pPtr =
                    new parcelType
                (
                    this->owner().pMesh(),
                    pos,
                    celli,
                    tetFacei,
                    tetPti
                );

                // Check/set new parcel thermo properties
                td.cloud().setParcelThermoProperties(*pPtr, 0.0);

                setParcelProperties(*pPtr, j);

                if (pPtr->nParticle() > 0.001)
                {
                    // Check new parcel properties
                    //                td.cloud().checkParcelProperties(*pPtr, 0.0, true);
                    td.cloud().checkParcelProperties(*pPtr, 0.0, false);

                    // Add the new parcel to the cloud
                    td.cloud().addParticle(pPtr);

                    nParcelsInjected_++;
                }
                else
                {
                    // TODO: cache mass and re-distribute?
                    delete pPtr;
                }
            }
        }
    }
}

void Foam::SurfaceFilmModel<CloudType>::inject(TrackData& td)
{
    typedef regionModels::surfaceFilmModels::surfaceFilmModel filmModelType;

    bool modelPresent =
        this->owner().db().objectRegistry::foundObject<filmModelType>
        (
            "surfaceFilmProperties"
        );

    if (!modelPresent)
    {
        return;
    }

    // Retrieve the film model from the owner database
    const filmModelType& filmModel =
        this->owner().db().objectRegistry::lookupObject<filmModelType>
        (
            "surfaceFilmProperties"
        );

    if (!filmModel.active())
    {
        return;
    }

    const labelList& filmPatches = filmModel.intCoupledPatchIDs();
    const labelList& primaryPatches = filmModel.primaryPatchIDs();

    const fvMesh& mesh = owner_.mesh();
    const polyBoundaryMesh& pbm = mesh.boundaryMesh();

    forAll(filmPatches, i)
    {
        const label filmPatchI = filmPatches[i];
        const label primaryPatchI = primaryPatches[i];

        const directMappedPatchBase& mapPatch =
            filmModel.mappedPatches()[filmPatchI];
        const mapDistribute& distMap = mapPatch.map();

        injectorCellsPatch_ = pbm[primaryPatchI].faceCells();

        cacheFilmFields(filmPatchI, primaryPatchI, distMap, filmModel);

        const vectorField& Cf = mesh.C().boundaryField()[primaryPatchI];
        const vectorField& Sf = mesh.Sf().boundaryField()[primaryPatchI];
        const scalarField& magSf = mesh.magSf().boundaryField()[primaryPatchI];

        forAll(injectorCellsPatch_, j)
        {
            if (diameterParcelPatch_[j] > 0)
            {
//                const point& pos = this->owner().mesh().C()[cellI];
                const scalar offset =
                    max
                    (
                        diameterParcelPatch_[j],
                        deltaFilmPatch_[primaryPatchI][j]
                    );
                const point pos = Cf[j] - 1.1*offset*Sf[j]/magSf[j];
                const label cellI = injectorCellsPatch_[j];

                // Create a new parcel
                typename CloudType::parcelType* pPtr =
                    new typename CloudType::parcelType(td.cloud(), pos, cellI);
                setParcelProperties(*pPtr, j);

                // Check new parcel properties
//                td.cloud().checkParcelProperties(*pPtr, 0.0, true);
                td.cloud().checkParcelProperties(*pPtr, 0.0, false);

                // Add the new parcel to the cloud
                td.cloud().addParticle(pPtr);

                nParcelsInjected_++;
            }
        }
    }
}

bool Foam::KinematicParcel<ParcelType>::move
(
    TrackData& td,
    const scalar trackTime
)
{
    typename TrackData::cloudType::parcelType& p =
        static_cast<typename TrackData::cloudType::parcelType&>(*this);

    td.switchProcessor = false;
    td.keepParticle = true;

    const polyMesh& mesh = td.cloud().pMesh();
    const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
    const scalarField& V = mesh.cellVolumes();
    const scalar maxCo = td.cloud().solution().maxCo();

    scalar tEnd = (1.0 - p.stepFraction())*trackTime;
    const scalar dtMax = tEnd;

    while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
    {
        // Apply correction to position for reduced-D cases
        meshTools::constrainToMeshCentre(mesh, p.position());

        // Set the Lagrangian time-step
        scalar dt = min(dtMax, tEnd);

        // Remember which cell the parcel is in since this will change if
        // a face is hit
        const label cellI = p.cell();

        const scalar magU = mag(U_);
        if (p.active() && magU > ROOTVSMALL)
        {
            const scalar d = dt*magU;
            const scalar dCorr = min(d, maxCo*cbrt(V[cellI]));
            dt *=
                dCorr/d
               *p.trackToFace(p.position() + dCorr*U_/magU, td);
        }

        tEnd -= dt;
        p.stepFraction() = 1.0 - tEnd/trackTime;

        // Avoid problems with extremely small timesteps
        if (dt > ROOTVSMALL)
        {
            // Update cell based properties
            p.setCellValues(td, dt, cellI);

            if (td.cloud().solution().cellValueSourceCorrection())
            {
                p.cellValueSourceCorrection(td, dt, cellI);
            }

            p.calc(td, dt, cellI);
        }

        if (p.onBoundary() && td.keepParticle)
        {
            if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
            {
                td.switchProcessor = true;
            }
        }

        p.age() += dt;

        td.cloud().functions().postMove(p, cellI, dt);
    }

    return td.keepParticle;
}

bool Foam::KinematicParcel<ParcelType>::move
(
    TrackData& td,
    const scalar trackTime
)
{
    typename TrackData::cloudType::parcelType& p =
        static_cast<typename TrackData::cloudType::parcelType&>(*this);

    td.switchProcessor = false;
    td.keepParticle = true;

    const polyMesh& mesh = td.cloud().pMesh();
    const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
    const scalarField& cellLengthScale = td.cloud().cellLengthScale();
    const scalar maxCo = td.cloud().solution().maxCo();

    scalar tEnd = (1.0 - p.stepFraction())*trackTime;
    scalar dtMax = trackTime;
    if (td.cloud().solution().transient())
    {
        dtMax *= maxCo;
    }

    bool tracking = true;
    label nTrackingStalled = 0;

    while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
    {
        // Apply correction to position for reduced-D cases
        meshTools::constrainToMeshCentre(mesh, p.position());

        const point start(p.position());

        // Set the Lagrangian time-step
        scalar dt = min(dtMax, tEnd);

        // Cache the parcel current cell as this will change if a face is hit
        const label cellI = p.cell();

        const scalar magU = mag(U_);
        if (p.active() && tracking && (magU > ROOTVSMALL))
        {
            const scalar d = dt*magU;
            const scalar dCorr = min(d, maxCo*cellLengthScale[cellI]);
            dt *=
                dCorr/d
               *p.trackToFace(p.position() + dCorr*U_/magU, td);
        }

        tEnd -= dt;

        scalar newStepFraction = 1.0 - tEnd/trackTime;

        if (tracking)
        {
            if
            (
                mag(p.stepFraction() - newStepFraction)
              < particle::minStepFractionTol
            )
            {
                nTrackingStalled++;

                if (nTrackingStalled > maxTrackAttempts)
                {
                    tracking = false;
                }
            }
            else
            {
                nTrackingStalled = 0;
            }
        }

        p.stepFraction() = newStepFraction;

        // Avoid problems with extremely small timesteps
        if (dt > ROOTVSMALL)
        {
            // Update cell based properties
            p.setCellValues(td, dt, cellI);

            if (td.cloud().solution().cellValueSourceCorrection())
            {
                p.cellValueSourceCorrection(td, dt, cellI);
            }

            p.calc(td, dt, cellI);
        }

        if (p.onBoundary() && td.keepParticle)
        {
            if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
            {
                td.switchProcessor = true;
            }
        }

        p.age() += dt;
//.........这里部分代码省略.........

void Foam::InjectionModel<CloudType>::inject(TrackData& td)
{
    if (!active())
    {
        return;
    }

    const scalar time = owner_.db().time().value();
    const scalar carrierDt = owner_.db().time().deltaTValue();
    const polyMesh& mesh = owner_.mesh();

    // Prepare for next time step
    label parcelsAdded = 0;
    scalar massAdded = 0.0;
    label newParcels = 0;
    scalar newVolume = 0.0;

    prepareForNextTimeStep(time, newParcels, newVolume);

    // Duration of injection period during this timestep
    const scalar deltaT =
        max(0.0, min(carrierDt, min(time - SOI_, timeEnd() - time0_)));

    // Pad injection time if injection starts during this timestep
    const scalar padTime = max(0.0, SOI_ - time0_);

    // Introduce new parcels linearly across carrier phase timestep
    for (label parcelI=0; parcelI<newParcels; parcelI++)
    {
        if (validInjection(parcelI))
        {
            // Calculate the pseudo time of injection for parcel 'parcelI'
            scalar timeInj = time0_ + padTime + deltaT*parcelI/newParcels;

            // Determine the injection position and owner cell
            label cellI = -1;
            vector pos = vector::zero;
            setPositionAndCell(parcelI, newParcels, timeInj, pos, cellI);

            if (cellI > -1)
            {
                // Lagrangian timestep
                scalar dt = time - timeInj;

                // Apply corrections to position for 2-D cases
                meshTools::constrainToMeshCentre(mesh, pos);

                // Create a new parcel
                parcelType* pPtr = new parcelType(td.cloud(), pos, cellI);

                // Assign new parcel properties in injection model
                setProperties(parcelI, newParcels, timeInj, *pPtr);

                // Check new parcel properties
                td.cloud().checkParcelProperties(*pPtr, dt, fullyDescribed());

                // Apply correction to velocity for 2-D cases
                meshTools::constrainDirection
                (
                    mesh,
                    mesh.solutionD(),
                    pPtr->U()
                );

                // Number of particles per parcel
                pPtr->nParticle() =
                    setNumberOfParticles
                    (
                        newParcels,
                        newVolume,
                        pPtr->d(),
                        pPtr->rho()
                    );

                // Add the new parcel
                td.cloud().addParticle(pPtr);

                massAdded += pPtr->nParticle()*pPtr->mass();
                parcelsAdded++;
            }
        }
    }

    postInjectCheck(parcelsAdded, massAdded);
}

Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
(
    TrackData& td,
    const scalar dt,
    const label cellI,
    const scalar Re,
    const scalar Pr,
    const scalar kappa,
    const scalar NCpW,
    const scalar Sh,
    scalar& dhsTrans,
    scalar& Sph
)
{
    if (!td.cloud().heatTransfer().active())
    {
        return T_;
    }

    const scalar d = this->d();
    const scalar rho = this->rho();

    // Calc heat transfer coefficient
    scalar htc = td.cloud().heatTransfer().htc(d, Re, Pr, kappa, NCpW);

    if (mag(htc) < ROOTVSMALL && !td.cloud().radiation())
    {
        return
            max
            (
                T_ + dt*Sh/(this->volume(d)*rho*Cp_),
                td.cloud().constProps().TMin()
            );
    }

    htc = max(htc, ROOTVSMALL);
    const scalar As = this->areaS(d);

    scalar ap = Tc_ + Sh/As/htc;
    scalar bp = 6.0*(Sh/As + htc*(Tc_ - T_));
    if (td.cloud().radiation())
    {
        tetIndices tetIs = this->currentTetIndices();
        const scalar Gc = td.GInterp().interpolate(this->position(), tetIs);
        const scalar sigma = physicoChemical::sigma.value();
        const scalar epsilon = td.cloud().constProps().epsilon0();

        ap = (ap + epsilon*Gc/(4.0*htc))/(1.0 + epsilon*sigma*pow3(T_)/htc);
        bp += 6.0*(epsilon*(Gc/4.0 - sigma*pow4(T_)));
    }
    bp /= rho*d*Cp_*(ap - T_) + ROOTVSMALL;

    // Integrate to find the new parcel temperature
    IntegrationScheme<scalar>::integrationResult Tres =
        td.cloud().TIntegrator().integrate(T_, dt, ap*bp, bp);

    scalar Tnew =
        min
        (
            max
            (
                Tres.value(),
                td.cloud().constProps().TMin()
            ),
            td.cloud().constProps().TMax()
        );

    Sph = dt*htc*As;

    dhsTrans += Sph*(Tres.average() - Tc_);

    return Tnew;
}

void Foam::ThermoParcel<ParcelType>::calc
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
    // Define local properties at beginning of time step
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    const scalar np0 = this->nParticle_;
    const scalar mass0 = this->mass();


    // Calc surface values
    // ~~~~~~~~~~~~~~~~~~~
    scalar Ts, rhos, mus, Pr, kappas;
    calcSurfaceValues(td, cellI, this->T_, Ts, rhos, mus, Pr, kappas);

    // Reynolds number
    scalar Re = this->Re(this->U_, this->d_, rhos, mus);


    // Sources
    // ~~~~~~~

    // Explicit momentum source for particle
    vector Su = vector::zero;

    // Linearised momentum source coefficient
    scalar Spu = 0.0;

    // Momentum transfer from the particle to the carrier phase
    vector dUTrans = vector::zero;

    // Explicit enthalpy source for particle
    scalar Sh = 0.0;

    // Linearised enthalpy source coefficient
    scalar Sph = 0.0;

    // Sensible enthalpy transfer from the particle to the carrier phase
    scalar dhsTrans = 0.0;


    // Heat transfer
    // ~~~~~~~~~~~~~

    // Sum Ni*Cpi*Wi of emission species
    scalar NCpW = 0.0;

    // Calculate new particle temperature
    this->T_ =
        this->calcHeatTransfer
        (
            td,
            dt,
            cellI,
            Re,
            Pr,
            kappas,
            NCpW,
            Sh,
            dhsTrans,
            Sph
        );


    // Motion
    // ~~~~~~

    // Calculate new particle velocity
    this->U_ =
        this->calcVelocity(td, dt, cellI, Re, mus, mass0, Su, dUTrans, Spu);


    //  Accumulate carrier phase source terms
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if (td.cloud().solution().coupled())
    {
        // Update momentum transfer
        td.cloud().UTrans()[cellI] += np0*dUTrans;

        // Update momentum transfer coefficient
        td.cloud().UCoeff()[cellI] += np0*Spu;

        // Update sensible enthalpy transfer
        td.cloud().hsTrans()[cellI] += np0*dhsTrans;

        // Update sensible enthalpy coefficient
        td.cloud().hsCoeff()[cellI] += np0*Sph;

        // Update radiation fields
        if (td.cloud().radiation())
        {
            const scalar ap = this->areaP();
            const scalar T4 = pow4(this->T_);
            td.cloud().radAreaP()[cellI] += dt*np0*ap;
            td.cloud().radT4()[cellI] += dt*np0*T4;
            td.cloud().radAreaPT4()[cellI] += dt*np0*ap*T4;
        }
//.........这里部分代码省略.........

void Foam::ReactingMultiphaseParcel<ParcelType>::calc
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
    const CompositionModel<reactingCloudType>& composition =
        td.cloud().composition();


    // Define local properties at beginning of timestep
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    const scalar np0 = this->nParticle_;
    const scalar d0 = this->d_;
    const vector& U0 = this->U_;
    const scalar T0 = this->T_;
    const scalar mass0 = this->mass();

    const scalar pc = this->pc_;

    const scalarField& YMix = this->Y_;
    const label idG = composition.idGas();
    const label idL = composition.idLiquid();
    const label idS = composition.idSolid();


    // Calc surface values
    scalar Ts, rhos, mus, Prs, kappas;
    this->calcSurfaceValues(td, cellI, T0, Ts, rhos, mus, Prs, kappas);
    scalar Res = this->Re(U0, d0, rhos, mus);


    // Sources
    //~~~~~~~~

    // Explicit momentum source for particle
    vector Su = vector::zero;

    // Linearised momentum source coefficient
    scalar Spu = 0.0;

    // Momentum transfer from the particle to the carrier phase
    vector dUTrans = vector::zero;

    // Explicit enthalpy source for particle
    scalar Sh = 0.0;

    // Linearised enthalpy source coefficient
    scalar Sph = 0.0;

    // Sensible enthalpy transfer from the particle to the carrier phase
    scalar dhsTrans = 0.0;


    // 1. Compute models that contribute to mass transfer - U, T held constant
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    // Phase change in liquid phase
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    // Mass transfer due to phase change
    scalarField dMassPC(YLiquid_.size(), 0.0);

    // Molar flux of species emitted from the particle (kmol/m^2/s)
    scalar Ne = 0.0;

    // Sum Ni*Cpi*Wi of emission species
    scalar NCpW = 0.0;

    // Surface concentrations of emitted species
    scalarField Cs(composition.carrier().species().size(), 0.0);

    // Calc mass and enthalpy transfer due to phase change
    this->calcPhaseChange
    (
        td,
        dt,
        cellI,
        Res,
        Prs,
        Ts,
        mus/rhos,
        d0,
        T0,
        mass0,
        idL,
        YMix[LIQ],
        YLiquid_,
        dMassPC,
        Sh,
        Ne,
        NCpW,
        Cs
    );


    // Devolatilisation
//.........这里部分代码省略.........