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

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

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

    scalar rho   = fuels.rho(pc, T, p.X());
    scalar sigma = fuels.sigma(pc, T, p.X());
    scalar mu    = fuels.mu(pc, T, p.X());

    // inverse of characteristic viscous damping time
    scalar rtd = 0.5*Cmu_*mu/(rho*r2);

    // oscillation frequency (squared)
    scalar omega2 = Comega_*sigma/(rho*r3) - rtd*rtd;

    if (omega2 > 0)
    {
        scalar omega = sqrt(omega2);
        scalar romega = 1.0/omega;

        scalar rhog = spray_.rho()[p.cell()];
        scalar We = p.We(Ug, rhog, sigma);
        scalar Wetmp = We/WeCrit_;

        scalar y1 = p.dev() - Wetmp;
        scalar y2 = p.ddev()*romega;

        scalar a = sqrt(y1*y1 + y2*y2);

        // scotty we may have break-up
        if (a+Wetmp > 1.0)
        {
            scalar phic = y1/a;

            // constrain phic within -1 to 1
            phic = max(min(phic, 1), -1);

            scalar phit = acos(phic);
            scalar phi = phit;
            scalar quad = -y2/a;
            if (quad < 0)
            {
                phi = 2*mathematicalConstant::pi - phit;
            }

            scalar tb = 0;

            if (mag(p.dev()) < 1.0)
            {
                scalar theta = acos((1.0 - Wetmp)/a);

                if (theta < phi)
                {
                    if (2*mathematicalConstant::pi-theta >= phi)
                    {
                        theta = -theta;
                    }
                    theta += 2*mathematicalConstant::pi;
                }
                tb = (theta-phi)*romega;

                // breakup occurs
                if (deltaT > tb)
                {
                    p.dev() = 1.0;
                    p.ddev() = -a*omega*sin(omega*tb + phi);
                }
            }

            // update droplet size
            if (deltaT > tb)
            {
                scalar sqrtWe = AWe_*pow(We, 4.0) + 1.0;
                scalar Kbr = k1_*omega*sqrtWe;

                if (We > WeTransition_)
                {
                    sqrtWe = sqrt(We);
                    Kbr =k2_*omega*sqrtWe;
                }

                scalar rWetmp = 1.0/Wetmp;
                scalar cosdtbu = max(-1.0, min(1.0, 1.0-rWetmp));
                scalar dtbu = romega*acos(cosdtbu);
                scalar decay = exp(-Kbr*dtbu);

                scalar rNew = decay*r;
                if (rNew < r)
                {
                    p.d() = 2*rNew;
//.........这里部分代码省略.........

void reitzDiwakar::breakupParcel
(
    parcel& p,
    const scalar deltaT,
    const vector& vel,
    const liquidMixture& fuels
) const
{
    /*
        These are the default values for this model...
        static const scalar Cbag   = 6.0;
        static const scalar Cb     = 0.785;
        static const scalar Cstrip = 0.5;
        static const scalar Cs     = 10.0;
    */

    const PtrList<volScalarField>& Y = spray_.composition().Y();

    label Ns = Y.size();
    label cellI = p.cell();
    scalar pressure = spray_.p()[cellI];
    scalar temperature = spray_.T()[cellI];
    scalar Taverage = p.T() + (temperature - p.T())/3.0;

    scalar muAverage = 0.0;
    scalar Winv = 0.0;
    for(label i=0; i<Ns; i++)
    {
        Winv += Y[i][cellI]/spray_.gasProperties()[i].W();
        muAverage += Y[i][cellI]*spray_.gasProperties()[i].mu(Taverage);
    }
    scalar R = specie::RR()*Winv;

    // ideal gas law to evaluate density
    scalar rhoAverage = pressure/R/Taverage;
    scalar nuAverage = muAverage/rhoAverage;
    scalar sigma = fuels.sigma(pressure, p.T(), p.X());

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    //     The We and Re numbers are to be evaluated using the 1/3 rule.
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    scalar WeberNumber = p.We(vel, rhoAverage, sigma);
    scalar ReynoldsNumber = p.Re(vel, nuAverage);

    scalar sqRey = sqrt(ReynoldsNumber);

    if (WeberNumber > Cbag_)
    {
        if (WeberNumber > Cstrip_*sqRey)
        {
            scalar dStrip =
                pow(2.0*Cstrip_*sigma, 2.0)/
                (
                    rhoAverage
                  * pow(mag(p.Urel(vel)), 3.0)
                  * muAverage
                );

            scalar tauStrip =
                Cs_ * p.d()
              * sqrt
                (
                    fuels.rho(pressure, p.T(), p.X())
                    / rhoAverage
                )
              / mag(p.Urel(vel));

            scalar fraction = deltaT/tauStrip;

            // new droplet diameter, implicit calculation
            p.d() = (fraction*dStrip + p.d())/(1.0 + fraction);
        }
        else
        {
            scalar dBag =
                2.0 * Cbag_ * sigma
              / (
                  rhoAverage
                * pow(mag(p.Urel(vel)), 2.0)
                );

            scalar tauBag =
                Cb_ * p.d()
                * sqrt
                  (
                      fuels.rho(pressure, p.T(), p.X())
                    * p.d()
                    / sigma
                  );

            scalar fraction = deltaT/tauBag;

            // new droplet diameter, implicit calculation
            p.d() = (fraction*dBag + p.d())/(1.0 + fraction);
        }

    }

}

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

    const PtrList<volScalarField>& Y = spray_.composition().Y();

    label Ns = Y.size();
    label cellI = p.cell();
    scalar pressure = spray_.p()[cellI];
    scalar temperature = spray_.T()[cellI];
    scalar Taverage = p.T() + (temperature - p.T())/3.0;

    scalar Winv = 0.0;
    for(label i=0; i<Ns; i++)
    {
        Winv += Y[i][cellI]/spray_.gasProperties()[i].W();
    }
    scalar R = specie::RR*Winv;

    // ideal gas law to evaluate density
    scalar rhoAverage = pressure/R/Taverage;
    scalar sigma = fuels.sigma(pressure, p.T(), p.X());

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    //     The We and Re numbers are to be evaluated using the 1/3 rule.
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    scalar rhoFuel = fuels.rho(1.0e+5, p.T(), p.X());

    scalar U = mag(p.Urel(vel));

    const injectorType& it =
        spray_.injectors()[label(p.injector())].properties();

    vector itPosition(vector::zero);
    label nHoles = it.nHoles();
    if (nHoles > 1)
    {
        for(label i=0; i<nHoles;i++)
        {
            itPosition += it.position(i);
        }
        itPosition /= nHoles;
    }
    else
    {
        itPosition = it.position(0);
    }
//    const vector itPosition = it.position();


    scalar lBU = B_ * sqrt
    (
        rhoFuel * sigma * p.d() * cos(angle_*mathematicalConstant::pi/360.0)
      / sqr(rhoAverage*U)
    );

    scalar pWalk = mag(p.position() - itPosition);

    if(pWalk > lBU && p.liquidCore() == 1.0)
    {
        p.liquidCore() = 0.0;
    }
}

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 myLISA_3_InjPos::atomizeParcel
(
    parcel& p,
    const scalar deltaT,
    const vector& vel,
    const liquidMixture& fuels
) const
{


    const PtrList<volScalarField>& Y = spray_.composition().Y();

    label Ns = Y.size();
    label cellI = p.cell();
    scalar pressure = spray_.p()[cellI];
    scalar temperature = spray_.T()[cellI];
//--------------------------------------AL____101015--------------------------------//
//  scalar Taverage = p.T() + (temperature - p.T())/3.0;
    scalar Taverage = temperature;
//-----------------------------------------END--------------------------------------//
    scalar Winv = 0.0;

    for(label i=0; i<Ns; i++)
    {
        Winv += Y[i][cellI]/spray_.gasProperties()[i].W();
    }

    scalar R = specie::RR*Winv;

    // ideal gas law to evaluate density
    scalar rhoAverage = pressure/R/Taverage;
    //scalar nuAverage = muAverage/rhoAverage;
    scalar sigma = fuels.sigma(pressure, p.T(), p.X());

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    //     The We and Re numbers are to be evaluated using the 1/3 rule.
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    scalar WeberNumber = p.We(vel, rhoAverage, sigma);

    scalar tau = 0.0;
    scalar dL = 0.0;
    scalar k = 0.0;
    scalar muFuel = fuels.mu(pressure, p.T(), p.X());
    scalar rhoFuel = fuels.rho(1.0e+5, p.T(), p.X());
    scalar nuFuel = muFuel/rhoFuel;

    vector uDir = p.U()/mag(p.U());

    scalar uGas = mag(vel & uDir);
    vector Ug = uGas*uDir;

    /*
        TL
        It might be the relative velocity between Liquid and Gas, but I use the
        absolute velocity of the parcel as suggested by the authors
    */

//    scalar U = mag(p.Urel(vel));
    scalar U = mag(p.U());
    p.ct() += deltaT;

    scalar Q = rhoAverage/rhoFuel;

    const injectorType& it =
        spray_.injectors()[label(p.injector())].properties();

    if (it.nHoles() > 1)
    {
        Info << "Warning: This atomization model is not suitable for multihole injector." << endl
             << "Only the first hole will be used." << endl;
    }

    const vector direction = it.direction(0, spray_.runTime().value());
//--------------------------------CH 101108--------------------------------------------------//
//    const vector itPosition = it.position(0);
    const injectorModel& im = spray_.injection();
    const vector itPosition = it.position(0) + im.injDist(0)*direction/mag(direction);
//------------------------------------END----------------------------------------------------//
    scalar pWalk = mag(p.position() - itPosition);

//  Updating liquid sheet tickness... that is the droplet diameter

//    const vector direction = it.direction(0, spray_.runTime().value());

    scalar h = (p.position() - itPosition) & direction;

    scalar d = sqrt(sqr(pWalk)-sqr(h));

    scalar time = pWalk/mag(p.U());

    scalar elapsedTime = spray_.runTime().value();

    scalar massFlow = it.massFlowRate(max(0.0,elapsedTime-time));

    scalar hSheet = massFlow/(mathematicalConstant::pi*d*rhoFuel*mag(p.U()));

    p.d() = min(hSheet,p.d());

    if(WeberNumber > 27.0/16.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;
//.........这里部分代码省略.........