C++ parcelType::f方法代码示例

void Foam::LambertWall<CloudType>::evaluatePendularWall
(
    typename CloudType::parcelType& p,
    const point& site,
    const WallSiteData<vector>& data,
    scalar pREff
) const
{
    const scalar& st = this->surfaceTension();
    const scalar& ca = this->contactAngle();
    const scalar& lf = this->liqFrac();
    const scalar& vis = this->viscosity();
    const scalar& ms = this->minSep();

    scalar Vtot = lf*(p.Vliq());

    vector r_PW = p.position() - site;

    vector U_PW = p.U() - data.wallData();

    scalar r_PW_mag = mag(r_PW);

    scalar normalOverlapMag = pREff - r_PW_mag;

    scalar S = -normalOverlapMag;

    vector rHat_PW = r_PW/(r_PW_mag + VSMALL);


    // Normal force
    scalar capMag =
        4*mathematical::pi*pREff*st*cos(ca)/
        (1+max(S, 0)*sqrt(mathematical::pi*pREff/Vtot));

    //Info << "the value of capMag is " << capMag << endl;
    //Info << " the value of overlapMag S is " << S << endl;
    //Info << " the volume of Vtot is " << Vtot << endl;

    scalar Svis = max(pREff*ms, S);

    scalar etaN = 6*mathematical::pi*vis*pREff*pREff/Svis;

    vector fN_PW = (-capMag - etaN*(U_PW & rHat_PW)) * rHat_PW;

    p.f() += fN_PW;

    vector UT_PW = U_PW - (U_PW & rHat_PW)*rHat_PW
                  - ((pREff*p.omega()) ^ rHat_PW);

    scalar etaT =
        6*mathematical::pi*vis*pREff*(8./15.*log(pREff/Svis) + 0.9588);

    vector fT_PW = -etaT * UT_PW;

    p.f() += fT_PW;

    p.torque() += (pREff*-rHat_PW) ^ fT_PW;

}

void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
(
    typename CloudType::parcelType& pA,
    typename CloudType::parcelType& pB
) const
{
    vector r_AB = (pA.position() - pB.position());

    scalar dAEff = pA.d();

    if (useEquivalentSize_)
    {
        dAEff *= cbrt(pA.nParticle()*volumeFactor_);
    }

    scalar dBEff = pB.d();

    if (useEquivalentSize_)
    {
        dBEff *= cbrt(pB.nParticle()*volumeFactor_);
    }

    scalar r_AB_mag = mag(r_AB);

    scalar normalOverlapMag = 0.5*(dAEff + dBEff) - r_AB_mag;

    if (normalOverlapMag > 0)
    {
        //Particles in collision

        vector rHat_AB = r_AB/(r_AB_mag + VSMALL);

        vector U_AB = pA.U() - pB.U();

        // Effective radius
        scalar R = 0.5*dAEff*dBEff/(dAEff + dBEff);

        // Effective mass
        scalar M = pA.mass()*pB.mass()/(pA.mass() + pB.mass());

        scalar kN = (4.0/3.0)*sqrt(R)*Estar_;

        scalar etaN = alpha_*sqrt(M*kN)*pow025(normalOverlapMag);

        // Normal force
        vector fN_AB =
            rHat_AB
           *(kN*pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));

        // Cohesion force
        if (cohesion_)
        {
            fN_AB +=
                -cohesionEnergyDensity_
                *overlapArea(dAEff/2.0, dBEff/2.0, r_AB_mag)
                *rHat_AB;
        }

        pA.f() += fN_AB;
        pB.f() += -fN_AB;

        vector USlip_AB =
            U_AB - (U_AB & rHat_AB)*rHat_AB
          + (pA.omega() ^ (dAEff/2*-rHat_AB))
          - (pB.omega() ^ (dBEff/2*rHat_AB));

        scalar deltaT = this->owner().mesh().time().deltaTValue();

        vector& tangentialOverlap_AB =
            pA.collisionRecords().matchPairRecord
            (
                pB.origProc(),
                pB.origId()
            ).collisionData();

        vector& tangentialOverlap_BA =
            pB.collisionRecords().matchPairRecord
            (
                pA.origProc(),
                pA.origId()
            ).collisionData();

        vector deltaTangentialOverlap_AB = USlip_AB*deltaT;

        tangentialOverlap_AB += deltaTangentialOverlap_AB;
        tangentialOverlap_BA += -deltaTangentialOverlap_AB;

        scalar tangentialOverlapMag = mag(tangentialOverlap_AB);

        if (tangentialOverlapMag > VSMALL)
        {
            scalar kT = 8.0*sqrt(R*normalOverlapMag)*Gstar_;

            scalar etaT = etaN;

            // Tangential force
            vector fT_AB;

            if (kT*tangentialOverlapMag > mu_*mag(fN_AB))
            {
//.........这里部分代码省略.........

void Foam::WallLocalSpringSliderDashpot<CloudType>::evaluateWall
(
    typename CloudType::parcelType& p,
    const point& site,
    const WallSiteData<vector>& data,
    scalar pREff
) const
{
    // wall patch index
    label wPI = patchMap_[data.patchIndex()];

    // data for this patch
    scalar Estar = Estar_[wPI];
    scalar Gstar = Gstar_[wPI];
    scalar alpha = alpha_[wPI];
    scalar b = b_[wPI];
    scalar mu = mu_[wPI];

    vector r_PW = p.position() - site;

    vector U_PW = p.U() - data.wallData();

    scalar normalOverlapMag = max(pREff - mag(r_PW), 0.0);

    vector rHat_PW = r_PW/(mag(r_PW) + VSMALL);

    scalar kN = (4.0/3.0)*sqrt(pREff)*Estar;

    scalar etaN = alpha*sqrt(p.mass()*kN)*pow025(normalOverlapMag);

    vector fN_PW =
        rHat_PW
       *(kN*pow(normalOverlapMag, b) - etaN*(U_PW & rHat_PW));

    p.f() += fN_PW;

    vector USlip_PW =
        U_PW - (U_PW & rHat_PW)*rHat_PW
      + (p.omega() ^ (pREff*-rHat_PW));

    scalar deltaT = this->owner().mesh().time().deltaTValue();

    vector& tangentialOverlap_PW =
        p.collisionRecords().matchWallRecord(-r_PW, pREff).collisionData();

    tangentialOverlap_PW += USlip_PW*deltaT;

    scalar tangentialOverlapMag = mag(tangentialOverlap_PW);

    if (tangentialOverlapMag > VSMALL)
    {
        scalar kT = 8.0*sqrt(pREff*normalOverlapMag)*Gstar;

        scalar etaT = etaN;

        // Tangential force
        vector fT_PW;

        if (kT*tangentialOverlapMag > mu*mag(fN_PW))
        {
            // Tangential force greater than sliding friction,
            // particle slips

            fT_PW = -mu*mag(fN_PW)*USlip_PW/mag(USlip_PW);

            tangentialOverlap_PW = vector::zero;
        }
        else
        {
            fT_PW =
                -kT*tangentialOverlapMag
               *tangentialOverlap_PW/tangentialOverlapMag
              - etaT*USlip_PW;
        }

        p.f() += fT_PW;

        p.torque() += (pREff*-rHat_PW) ^ fT_PW;
    }
}