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

void SHF::breakupParcel
(
    parcel& p,
    const scalar deltaT,
    const vector& vel,
    const liquidMixture& 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 acceleration = p.Urel(vel)/p.tMom();
    vector trajectory = p.U()/mag(p.U());

    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;

    bool bag = false;
    bool multimode = false;
    bool shear = false;
    bool success = false;


    //  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_));

    if(weGas > weC && weGas < weB)
    {
        bag = true;
    }

    if(weGas >= weB && weGas <= weMM)
    {
        multimode = true;
    }

    if(weGas > weMM)
    {
        shear = true;
    }

    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;
//.........这里部分代码省略.........

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;
//.........这里部分代码省略.........