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

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;
    }
}

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

            p.d() = d;
            p.ct() = 0.0;
        }

        if (shear)
        {
            scalar dC = weConst_*sigma/(rhoGas*sqr(mag(vRel)));
            scalar d32Red = 4.0*(d32*dC)/(5.0*dC - d32);
            scalar initMass = p.m();

            scalar d05 = d32Coeff_*d32Red;

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

            do
            {

                x = cDmaxS_*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);

            p.d() = dC;
            p.m() = corePerc_*initMass;

            spray_.addParticle
            (
                new parcel
                (
                    p.mesh(),
                    p.position(),
                    p.cell(),
                    p.tetFace(),
                    p.tetPt(),
                    p.n(),
                    d,
                    p.T(),
                    (1.0 - corePerc_)*initMass,
                    0.0,
                    0.0,
                    0.0,
                    -GREAT,
                    p.tTurb(),
                    0.0,
                    scalar(p.injector()),
                    p.U(),
                    p.Uturb(),
                    p.X(),
                    p.fuelNames()
                )
            );

            p.ct() = 0.0;
        }
    }
}

// Return 'keepParcel'
bool reflectParcel::wallTreatment
(
    parcel& p,
    const label globalFacei
) const
{
    label patchi = p.patch(globalFacei);
    label facei = p.patchFace(patchi, globalFacei);

    const polyMesh& mesh = spray_.mesh();

    if (mesh_.boundaryMesh()[patchi].isWall())
    {
        // wallNormal defined to point outwards of domain
        vector Sf = mesh_.Sf().boundaryField()[patchi][facei];
        Sf /= mag(Sf);

        if (!mesh.moving())
        {
            // static mesh
            scalar Un = p.U() & Sf;

            if (Un > 0)
            {
                p.U() -= (1.0 + elasticity_)*Un*Sf;
            }

        }
        else
        {
            // moving mesh
            vector Ub1 = U_.boundaryField()[patchi][facei];
            vector Ub0 = U_.oldTime().boundaryField()[patchi][facei];

            scalar dt = spray_.runTime().deltaT().value();
            const vectorField& oldPoints = mesh.oldPoints();

            const vector& Cf1 = mesh.faceCentres()[globalFacei];

            vector Cf0 = mesh.faces()[globalFacei].centre(oldPoints);
            vector Cf = Cf0 + p.stepFraction()*(Cf1 - Cf0);
            vector Sf0 = mesh.faces()[globalFacei].normal(oldPoints);

            // for layer addition Sf0 = vector::zero and we use Sf
            if (mag(Sf0) > SMALL)
            {
                Sf0 /= mag(Sf0);
            }
            else
            {
                Sf0 = Sf;
            }

            scalar magSfDiff = mag(Sf - Sf0);

            vector Ub = Ub0 + p.stepFraction()*(Ub1 - Ub0);

            if (magSfDiff > SMALL)
            {
                // rotation + translation
                vector Sfp = Sf0 + p.stepFraction()*(Sf - Sf0);

                vector omega = Sf0 ^ Sf;
                scalar magOmega = mag(omega);
                omega /= magOmega+SMALL;

                scalar phiVel = ::asin(magOmega)/dt;

                scalar dist = (p.position() - Cf) & Sfp;
                vector pos = p.position() - dist*Sfp;
                vector vrot = phiVel*(omega ^ (pos - Cf));

                vector v = Ub + vrot;

                scalar Un = ((p.U() - v) & Sfp);

                if (Un > 0.0)
                {
                    p.U() -= (1.0 + elasticity_)*Un*Sfp;
                }
            }
            else
            {
                // translation
                vector Ur = p.U() - Ub;
                scalar Urn = Ur & Sf;
                /*
                if (mag(Ub-v) > SMALL)
                {
                    Info << "reflectParcel:: v = " << v
                        << ", Ub = " << Ub
                        << ", facei = " << facei
                        << ", patchi = " << patchi
                        << ", globalFacei = " << globalFacei
                        << ", name = " << mesh_.boundaryMesh()[patchi].name()
                        << endl;
                }
                    */
                if (Urn > 0.0)
//.........这里部分代码省略.........

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

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

            // Total number of parcels for the whole injection event
            label nParcels =
                spray_.injectors()[injector].properties()->nParcelsToInject
                (
                    spray_.injectors()[injector].properties()->tsoi(),
                    spray_.injectors()[injector].properties()->teoi()
                );

            scalar averageParcelMass =
                spray_.injectors()[injector].properties()->mass()/nParcels;

            if (p.ms()/averageParcelMass > msLimit_)
            {
                // set the initial ms value to -GREAT. This prevents
                // new droplets from being formed from the child droplet
                // from the KH instability

                // mass of stripped child parcel
                scalar mc = p.ms();
                // Prevent child parcel from taking too much mass
                if (mc > 0.5*p.m())
                {
                    mc = 0.5*p.m();
                }

                spray_.addParticle
                (
                    new parcel
                    (
                        spray_,
                        p.position(),
                        p.cell(),
                        p.n(),
                        dc,
                        p.T(),
                        mc,
                        0.0,
                        0.0,
                        0.0,
                        -GREAT,
                        p.tTurb(),
                        0.0,
                        p.injector(),
                        p.U(),
                        p.Uturb(),
                        p.X(),
                        p.fuelNames()
                    )
                );

                p.m() -= mc;
                p.ms() = 0.0;
            }
        }
    }
}