C++ volScalarField::write方法代码示例

void thixotropicViscosity::correct
(
    const volScalarField& p,
    const volScalarField& T
)
{
    const kinematicSingleLayer& film = filmType<kinematicSingleLayer>();

    const volVectorField& U = film.U();
    const volVectorField& Uw = film.Uw();
    const volScalarField& delta = film.delta();
    const volScalarField& deltaRho = film.deltaRho();
    const surfaceScalarField& phi = film.phi();
    const volScalarField& alpha = film.alpha();
    const Time& runTime = this->film().regionMesh().time();

    // Shear rate
    const volScalarField gDot
    (
        "gDot",
        alpha*mag(U - Uw)/(delta + film.deltaSmall())
    );

    if (debug && runTime.writeTime())
    {
        gDot.write();
    }

    const dimensionedScalar deltaRho0
    (
        "deltaRho0",
        deltaRho.dimensions(),
        rootVSmall
    );

    const surfaceScalarField phiU(phi/fvc::interpolate(deltaRho + deltaRho0));

    const dimensionedScalar c0("c0", dimless/dimTime, rootVSmall);
    const volScalarField coeff("coeff", -c_*pow(gDot, d_) + c0);

    fvScalarMatrix lambdaEqn
    (
        fvm::ddt(lambda_)
      + fvm::div(phiU, lambda_)
      - fvm::Sp(fvc::div(phiU), lambda_)
      ==
        a_*pow((1 - lambda_), b_)
      + fvm::SuSp(coeff, lambda_)

        // Include the effect of the impinging droplets added with lambda = 0
      - fvm::Sp
        (
            max
            (
               -film.rhoSp(),
                dimensionedScalar(film.rhoSp().dimensions(), 0)
            )/(deltaRho + deltaRho0),
            lambda_
        )
    );

    lambdaEqn.relax();
    lambdaEqn.solve();

    lambda_.min(1);
    lambda_.max(0);

    mu_ = muInf_/(sqr(1 - K_*lambda_) + rootVSmall);
    mu_.correctBoundaryConditions();
}

int main(int argc, char *argv[])
{
    timeSelector::addOptions();

    #include "setRootCase.H"
    #include "createTime.H"

    instantList timeDirs = timeSelector::select0(runTime, args);

    #include "createMesh.H"
    #include "createFields.H"

    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);

        Info<< "Time = " << runTime.timeName() << endl;

        // Cache the turbulence fields

        Info<< "\nRetrieving field k from turbulence model" << endl;
        const volScalarField k = RASModel->k();

        Info<< "\nRetrieving field epsilon from turbulence model" << endl;
        const volScalarField epsilon = RASModel->epsilon();

        Info<< "\nRetrieving field R from turbulence model" << endl;
        const volSymmTensorField R = RASModel->R();

        // Check availability of tubulence fields

        if (!IOobject("k", runTime.timeName(), mesh).headerOk())
        {
            Info<< "\nWriting turbulence field k" << endl;
            k.write();
        }
        else
        {
            Info<< "\nTurbulence k field already exists" << endl;
        }

        if (!IOobject("epsilon", runTime.timeName(), mesh).headerOk())
        {
            Info<< "\nWriting turbulence field epsilon" << endl;
            epsilon.write();
        }
        else
        {
            Info<< "\nTurbulence epsilon field already exists" << endl;
        }

        if (!IOobject("R", runTime.timeName(), mesh).headerOk())
        {
            Info<< "\nWriting turbulence field R" << endl;
            R.write();
        }
        else
        {
            Info<< "\nTurbulence R field already exists" << endl;
        }

        if (!IOobject("omega", runTime.timeName(), mesh).headerOk())
        {
            const scalar Cmu = 0.09;

            Info<< "creating omega" << endl;
            volScalarField omega
            (
                IOobject
                (
                    "omega",
                    runTime.timeName(),
                    mesh
                ),
                epsilon/(Cmu*k),
                epsilon.boundaryField().types()
            );
            Info<< "\nWriting turbulence field omega" << endl;
            omega.write();
        }
        else
        {
            Info<< "\nTurbulence omega field already exists" << endl;
        }
    }

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

    // Update omega and G at the wall
    omega_.boundaryFieldRef().updateCoeffs();

    const volScalarField CDkOmega
        ( "CD",
        (2*this->alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
        );
    
    const volScalarField F1("F1", this->F1(CDkOmega));

    {
        volScalarField gamma(this->gamma(F1));
        volScalarField beta(this->beta(F1));

        // Turbulent frequency equation
        tmp<fvScalarMatrix> omegaEqn
        (
            fvm::ddt(alpha, rho, omega_)
          + fvm::div(alphaRhoPhi, omega_)
          - fvm::laplacian(alpha*rho*this->DomegaEff(F1), omega_)
         ==
            alpha*rho*gamma*S*W
          - fvm::Sp(alpha*rho*beta*omega_, omega_)
          - fvm::SuSp
            (
                alpha*rho*(F1 - scalar(1))*CDkOmega/omega_,
                omega_
            )
        );

        omegaEqn.ref().relax();

        omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());

        solve(omegaEqn);
        bound(omega_, this->omegaMin_);
    }

    // Turbulent kinetic energy equation
    const volScalarField FonLim(
        "FonLim",
        min( max(sqr(this->y_)*S/this->nu() / (
            2.2*ReThetacLim_) - 1., 0.), 3.)
    );
    const volScalarField PkLim(
        "PkLim",
        5*Ck_ * max(gammaInt()-0.2,0.) * (1-gammaInt()) * FonLim * 
        max(3*CSEP_*this->nu() - this->nut_, 0.*this->nut_) * S * W
    );

    tmp<fvScalarMatrix> kEqn
    (
        fvm::ddt(alpha, rho, k_)
      + fvm::div(alphaRhoPhi, k_)
      - fvm::laplacian(alpha*rho*this->DkEff(F1), k_)
     ==
        alpha*rho*(G*gammaInt() + PkLim)
        - fvm::Sp(max(gammaInt(),scalar(0.1)) * alpha*rho*this->betaStar_*omega_, k_)
    );

    kEqn.ref().relax();
    solve(kEqn);
    bound(k_, this->kMin_);

    this->correctNut(S2, this->F23());

   // Intermittency equation (2)
    volScalarField Pgamma1 = Flength_ * S * gammaInt_ * Fonset(S);
    volScalarField Pgamma2 = ca2_ * W * gammaInt_ * Fturb();
    tmp<fvScalarMatrix> gammaEqn
        (
            fvm::ddt(alpha, rho, gammaInt_)
            + fvm::div(alphaRhoPhi, gammaInt_)
            - fvm::laplacian(alpha*rho*this->DgammaEff(), gammaInt_)
            ==
            alpha*rho*Pgamma1 - fvm::Sp(alpha*rho*Pgamma1, gammaInt_) +
            alpha*rho*Pgamma2 - fvm::Sp(alpha*rho*ce2_*Pgamma2, gammaInt_)
        ); 
    
    gammaEqn.ref().relax();
    solve(gammaEqn);

    bound(gammaInt_,scalar(0));

    if (debug && this->runTime_.outputTime()) {
        S.write();
        W.write();
        F1.write();
        CDkOmega.write();
        const volScalarField Pgamma("Pgamma", Pgamma1*(scalar(1)-gammaInt_));
        Pgamma.write();
        const volScalarField Egamma("Egamma", Pgamma2*(scalar(1)-ce2_*gammaInt_));
        Egamma.write();
        FonLim.write();
        PkLim.write();
        Fonset(S)().write();
        FPG()().write();
    }
}