C++ parcel::m方法代码示例

void Foam::SHF::breakupParcel
(
    parcel& p,
    const scalar deltaT,
    const vector& vel,
    const liquidMixtureProperties& fuels
) const
{
    label cellI = p.cell();
    scalar T = p.T();
    scalar pc = spray_.p()[cellI];

    scalar sigma = fuels.sigma(pc, T, p.X());
    scalar rhoLiquid = fuels.rho(pc, T, p.X());
    scalar muLiquid = fuels.mu(pc, T, p.X());
    scalar rhoGas = spray_.rho()[cellI];

    scalar weGas = p.We(vel, rhoGas, sigma);
    scalar weLiquid = p.We(vel, rhoLiquid, sigma);

    // correct the Reynolds number. Reitz is using radius instead of diameter

    scalar reLiquid   = p.Re(rhoLiquid, vel, muLiquid);
    scalar ohnesorge  = sqrt(weLiquid)/(reLiquid + VSMALL);

    vector vRel = p.Urel(vel);

    scalar weGasCorr = weGas/(1.0 + weCorrCoeff_*ohnesorge);

    // droplet deformation characteristic time

    scalar tChar = p.d()/mag(vRel)*sqrt(rhoLiquid/rhoGas);

    scalar tFirst = cInit_*tChar;

    scalar tSecond = 0;
    scalar tCharSecond = 0;


    //  updating the droplet characteristic time
    p.ct() += deltaT;

    if (weGas > weConst_)
    {
        if (weGas < weCrit1_)
        {
            tCharSecond = c1_*pow((weGas - weConst_),cExp1_);
        }
        else if (weGas >= weCrit1_ && weGas <= weCrit2_)
        {
            tCharSecond = c2_*pow((weGas - weConst_),cExp2_);
        }
        else
        {
            tCharSecond = c3_*pow((weGas - weConst_),cExp3_);
        }
    }

    scalar weC = weBuCrit_*(1.0+ohnCoeffCrit_*pow(ohnesorge, ohnExpCrit_));
    scalar weB = weBuBag_*(1.0+ohnCoeffBag_*pow(ohnesorge, ohnExpBag_));
    scalar weMM = weBuMM_*(1.0+ohnCoeffMM_*pow(ohnesorge, ohnExpMM_));

    bool bag = (weGas > weC && weGas < weB);

    bool multimode = (weGas >= weB && weGas <= weMM);

    bool shear = (weGas > weMM);

    tSecond = tCharSecond*tChar;

    scalar tBreakUP = tFirst + tSecond;
    if (p.ct() > tBreakUP)
    {

        scalar d32 =
            coeffD_*p.d()*pow(ohnesorge, onExpD_)*pow(weGasCorr, weExpD_);

        if (bag || multimode)
        {

            scalar d05 = d32Coeff_*d32;

            scalar x = 0.0;
            scalar y = 0.0;
            scalar d = 0.0;
            scalar px = 0.0;

            do
            {
                x = cDmaxBM_*rndGen_.sample01<scalar>();
                d = sqr(x)*d05;
                y = rndGen_.sample01<scalar>();

                px =
                    x
                   /(2.0*sqrt(constant::mathematical::twoPi)*sigma_)
                   *exp(-0.5*sqr((x-mu_)/sigma_));

            } while (y >= px);

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

void reitzKHRT::breakupParcel
(
    parcel& p,
    const scalar deltaT,
    const vector& vel,
    const liquidMixture& fuels
) const
{

    label celli = p.cell();
    scalar T = p.T();
    scalar r = 0.5*p.d();
    scalar pc = spray_.p()[celli];

    scalar sigma = fuels.sigma(pc, T, p.X());
    scalar rhoLiquid = fuels.rho(pc, T, p.X());
    scalar muLiquid = fuels.mu(pc, T, p.X());
    scalar rhoGas = spray_.rho()[celli];
    scalar Np = p.N(rhoLiquid);
    scalar semiMass = Np*pow(p.d(), 3.0);

    scalar weGas      = p.We(vel, rhoGas, sigma);
    scalar weLiquid   = p.We(vel, rhoLiquid, sigma);
    // correct the Reynolds number. Reitz is using radius instead of diameter
    scalar reLiquid   = 0.5*p.Re(rhoLiquid, vel, muLiquid);
    scalar ohnesorge  = sqrt(weLiquid)/(reLiquid + VSMALL);
    scalar taylor     = ohnesorge*sqrt(weGas);

    vector acceleration = p.Urel(vel)/p.tMom();
    vector trajectory = p.U()/mag(p.U());
    scalar gt = (g_ + acceleration) & trajectory;

    // frequency of the fastest growing KH-wave
    scalar omegaKH =
        (0.34 + 0.38*pow(weGas, 1.5))
       /((1 + ohnesorge)*(1 + 1.4*pow(taylor, 0.6)))
       *sqrt(sigma/(rhoLiquid*pow(r, 3)));

    // corresponding KH wave-length.
    scalar lambdaKH =
        9.02
       *r
       *(1.0 + 0.45*sqrt(ohnesorge))
       *(1.0 + 0.4*pow(taylor, 0.7))
       /pow(1.0 + 0.865*pow(weGas, 1.67), 0.6);

    // characteristic Kelvin-Helmholtz breakup time
    scalar tauKH = 3.726*b1_*r/(omegaKH*lambdaKH);

    // stable KH diameter
    scalar dc = 2.0*b0_*lambdaKH;

    // the frequency of the fastest growing RT wavelength.
    scalar helpVariable = mag(gt*(rhoLiquid - rhoGas));
    scalar omegaRT = sqrt
    (
        2.0*pow(helpVariable, 1.5)
       /(3.0*sqrt(3.0*sigma)*(rhoGas + rhoLiquid))
    );

    // RT wave number
    scalar KRT = sqrt(helpVariable/(3.0*sigma + VSMALL));

    // wavelength of the fastest growing RT frequency
    scalar lambdaRT = 2.0*mathematicalConstant::pi*cRT_/(KRT + VSMALL);

    // if lambdaRT < diameter, then RT waves are growing on the surface
    // and we start to keep track of how long they have been growing
    if ((p.ct() > 0) || (lambdaRT < p.d()))
    {
        p.ct() += deltaT;
    }

    // characteristic RT breakup time
    scalar tauRT = cTau_/(omegaRT + VSMALL);

    // check if we have RT breakup
    if ((p.ct() > tauRT) && (lambdaRT < p.d()))
    {
        // the RT breakup creates diameter/lambdaRT new droplets
        p.ct() = -GREAT;
        scalar multiplier = p.d()/lambdaRT;
        scalar nDrops = multiplier*Np;
        p.d() = cbrt(semiMass/nDrops);
    }
    // otherwise check for KH breakup
    else if (dc < p.d())
    {
        // no breakup below Weber = 12
        if (weGas > weberLimit_)
        {

            label injector = label(p.injector());
            scalar fraction = deltaT/tauKH;

            // reduce the diameter according to the rate-equation
            p.d() = (fraction*dc + p.d())/(1.0 + fraction);

            scalar ms = rhoLiquid*Np*pow3(dc)*mathematicalConstant::pi/6.0;
            p.ms() += ms;
//.........这里部分代码省略.........