自学教程：C++ twoSymm函数代码示例

void Giesekus_Cross::correct(){    // Velocity gradient tensor    volTensorField L = fvc::grad(U());    // Convected derivate term    volTensorField C = tau_ & L;    // Twice the rate of deformation tensor    volSymmTensorField twoD = twoSymm(L);    // Two phase transport properties treatment    volScalarField alpha1f = min(max(alpha(), scalar(0)), scalar(1));    volScalarField lambda =  alpha1f*lambda1_ + (scalar(1) - alpha1f)*lambda2_;    volScalarField etaP =  alpha1f*etaP1_ + (scalar(1) - alpha1f)*etaP2_;    volScalarField alpha =  alpha1f*alpha1_ + (scalar(1) - alpha1f)*alpha2_;     // Stress transport equation    tmp<fvSymmTensorMatrix> tauEqn    (        fvm::ddt(lambda, tau_)        + lambda * fvm::div(phi(), tau_)        ==        etaP * twoD        + lambda * twoSymm( C )        - (alpha * lambda / mu()) * ( tau_ & tau_)        - fvm::Sp( scalar(1), tau_ )    );    Info << "test 0" << endl;    tauEqn().relax();    solve(tauEqn);}

void DeardorffDiffStress<BasicTurbulenceModel>::correct(){    if (!this->turbulence_)    {        return;    }    // Local references    const alphaField& alpha = this->alpha_;    const rhoField& rho = this->rho_;    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;    const volVectorField& U = this->U_;    volSymmTensorField& R = this->R_;    fv::options& fvOptions(fv::options::New(this->mesh_));    ReynoldsStress<LESModel<BasicTurbulenceModel>>::correct();    tmp<volTensorField> tgradU(fvc::grad(U));    const volTensorField& gradU = tgradU();    volSymmTensorField D(symm(gradU));    volSymmTensorField P(-twoSymm(R & gradU));    volScalarField k(this->k());    tmp<fvSymmTensorMatrix> REqn    (        fvm::ddt(alpha, rho, R)      + fvm::div(alphaRhoPhi, R)      - fvm::laplacian(I*this->nu() + Cs_*(k/this->epsilon())*R, R)      + fvm::Sp(Cm_*alpha*rho*sqrt(k)/this->delta(), R)     ==        alpha*rho*P      + (4.0/5.0)*alpha*rho*k*D      - ((2.0/3.0)*(1.0 - Cm_/this->Ce_)*I)*(alpha*rho*this->epsilon())      + fvOptions(alpha, rho, R)    );    REqn.ref().relax();    fvOptions.constrain(REqn.ref());    REqn.ref().solve();    fvOptions.correct(R);    this->boundNormalStress(R);    correctNut();}

tmp<volSymmTensorField> LienCubicKE::devReff() const{    return tmp<volSymmTensorField>           (               new volSymmTensorField               (                   IOobject                   (                       "devRhoReff",                       runTime_.timeName(),                       mesh_,                       IOobject::NO_READ,                       IOobject::NO_WRITE                   ),                   -nuEff()*dev(twoSymm(fvc::grad(U_))) + nonlinearStress_               )           );}

tmp<volSymmTensorField> laminar::devRhoReff() const{    return tmp<volSymmTensorField>    (        new volSymmTensorField        (            IOobject            (                "devRhoReff",                runTime_.timeName(),                U_.db(),                IOobject::NO_READ,                IOobject::NO_WRITE            ),           -mu()*dev(twoSymm(fvc::grad(U_)))        )    );}

tmp<volSymmTensorField> GenSGSStress::devBeff() const{    return tmp<volSymmTensorField>    (        new volSymmTensorField        (            IOobject            (                "devRhoReff",                runTime_.timeName(),                U_.db(),                IOobject::NO_READ,                IOobject::NO_WRITE            ),            B_ - nu()*dev(twoSymm(fvc::grad(U())))        )    );}

tmp<volSymmTensorField> tatunRNG::devRhoReff() const{    return tmp<volSymmTensorField>    (        new volSymmTensorField        (            IOobject            (                "devRhoReff",                runTime_.timeName(),                mesh_,                IOobject::NO_READ,                IOobject::NO_WRITE            ),           -muEff()*dev(twoSymm(fvc::grad(U_)))        )    );}

void kEqn<BasicTurbulenceModel>::correct(){    if (!this->turbulence_)    {        return;    }    // Local references    const alphaField& alpha = this->alpha_;    const rhoField& rho = this->rho_;    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;    const volVectorField& U = this->U_;    volScalarField& nut = this->nut_;    fv::options& fvOptions(fv::options::New(this->mesh_));    LESeddyViscosity<BasicTurbulenceModel>::correct();    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));    tmp<volTensorField> tgradU(fvc::grad(U));    volScalarField G(this->GName(), nut*(tgradU() && dev(twoSymm(tgradU()))));    tgradU.clear();    tmp<fvScalarMatrix> kEqn    (        fvm::ddt(alpha, rho, k_)      + fvm::div(alphaRhoPhi, k_)      - fvm::laplacian(alpha*rho*DkEff(), k_)     ==        alpha*rho*G      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)      - fvm::Sp(this->Ce_*alpha*rho*sqrt(k_)/this->delta(), k_)      + kSource()      + fvOptions(alpha, rho, k_)    );    kEqn.ref().relax();    fvOptions.constrain(kEqn.ref());    solve(kEqn);    fvOptions.correct(k_);    bound(k_, this->kMin_);    correctNut();}

tmp<volSymmTensorField> kEpsilon::R() const{    return tmp<volSymmTensorField>    (        new volSymmTensorField        (            IOobject            (                "R",                runTime_.timeName(),                mesh_,                IOobject::NO_READ,                IOobject::NO_WRITE            ),            ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),            k_.boundaryField().types()        )    );}

void Foam::Maxwell::correct(){    // Velocity gradient tensor    volTensorField L = fvc::grad(U());    // Twice the rate of deformation tensor    volSymmTensorField twoD = twoSymm(L);     // Stress transport equation    fvSymmTensorMatrix tauEqn    (        fvm::ddt(tau_)     ==        etaP_/lambda_*twoD      - fvm::Sp( 1/lambda_, tau_ )    );    tauEqn.relax();    tauEqn.solve();}

Foam::tmp<Foam::volSymmTensorField>Foam::eddyViscosity<BasicTurbulenceModel>::devRhoReff() const{    return tmp<volSymmTensorField>    (        new volSymmTensorField        (            IOobject            (                IOobject::groupName("devRhoReff", this->U_.group()),                this->runTime_.timeName(),                this->mesh_,                IOobject::NO_READ,                IOobject::NO_WRITE            ),            (-(this->alpha_*this->rho_*this->nuEff()))           *dev(twoSymm(fvc::grad(this->U_)))        )    );}

Foam::tmp<Foam::volSymmTensorField>Foam::eddyViscosity<BasicTurbulenceModel>::R() const{    tmp<volScalarField> tk(k());    // Get list of patchField type names from k    wordList patchFieldTypes(tk().boundaryField().types());    // For k patchField types which do not have an equivalent for symmTensor    // set to calculated    forAll(patchFieldTypes, i)    {        if        (           !fvPatchField<symmTensor>::patchConstructorTablePtr_                ->found(patchFieldTypes[i])        )        {            patchFieldTypes[i] = calculatedFvPatchField<symmTensor>::typeName;        }    }    return tmp<volSymmTensorField>    (        new volSymmTensorField        (            IOobject            (                IOobject::groupName("R", this->U_.group()),                this->runTime_.timeName(),                this->mesh_,                IOobject::NO_READ,                IOobject::NO_WRITE,                false            ),            ((2.0/3.0)*I)*tk() - (nut_)*dev(twoSymm(fvc::grad(this->U_))),            patchFieldTypes        )    );}

tmp<volSymmTensorField>continuousGasKEpsilon<BasicTurbulenceModel>::R() const{    tmp<volScalarField> tk(this->k());    return tmp<volSymmTensorField>    (        new volSymmTensorField        (            IOobject            (                IOobject::groupName("R", this->U_.group()),                this->runTime_.timeName(),                this->mesh_,                IOobject::NO_READ,                IOobject::NO_WRITE            ),            ((2.0/3.0)*I)*tk() - (nutEff_)*dev(twoSymm(fvc::grad(this->U_))),            tk().boundaryField().types()        )    );}

Foam::tmp<Foam::volSymmTensorField>Foam::eddyViscosity<BasicTurbulenceModel>::R() const{    tmp<volScalarField> tk(k());    return tmp<volSymmTensorField>    (        new volSymmTensorField        (            IOobject            (                IOobject::groupName("R", this->U_.group()),                this->runTime_.timeName(),                this->mesh_,                IOobject::NO_READ,                IOobject::NO_WRITE            ),            ((2.0/3.0)*I)*tk() - (nut_)*dev(twoSymm(fvc::grad(this->U_))),            tk().boundaryField().types()        )    );}

Foam::tmp<Foam::volSymmTensorField>Foam::rasVofForceAndTorqueFunctionObject::devRhoReff() const{    const fvMesh& mesh =        time_.lookupObject<fvMesh>(regionName_);    if (mesh.foundObject<incompressible::RASModel>("RASProperties"))    {        const incompressible::RASModel& ras =            mesh.lookupObject<incompressible::RASModel>("RASProperties");        return ras.devReff();    }    else if (mesh.foundObject<incompressible::LESModel>("LESProperties"))    {        const incompressible::LESModel& les =            mesh.lookupObject<incompressible::LESModel>("LESProperties");        return les.devBeff();    }    else if (mesh.foundObject<twoPhaseMixture>("transportProperties"))    {        const twoPhaseMixture& twoPhaseProperties =            mesh.lookupObject<twoPhaseMixture>("transportProperties");        const volVectorField& U = mesh.lookupObject<volVectorField>("U");        return twoPhaseProperties.nu()*dev(twoSymm(fvc::grad(U)));    }    else    {        FatalErrorIn("floatingBody::devRhoReff()")            << "No valid model for viscous stress calculation."            << exit(FatalError);        return volSymmTensorField::null();    }}

tmp<volSymmTensorField> GenEddyVisc::devBeff() const{    return -nuEff()*dev(twoSymm(fvc::grad(U())));}

tmp<volSymmTensorField> GenEddyVisc::B() const{    return ((2.0/3.0)*I)*k() - nuSgs_*twoSymm(fvc::grad(U()));}

//- evaluate gradients needed for optimisationvoid surfaceOptimizer::evaluateGradients(    const scalar& K,    vector& gradF,    tensor& gradGradF) const{    gradF = vector::zero;    gradGradF = tensor::zero;    tensor gradGradLt(tensor::zero);    gradGradLt.xx() = 4.0;    gradGradLt.yy() = 4.0;    forAll(trias_, triI)    {        const point& p0 = pts_[trias_[triI][0]];        const point& p1 = pts_[trias_[triI][1]];        const point& p2 = pts_[trias_[triI][2]];        if( magSqr(p1 - p2) < VSMALL ) continue;        const scalar LSqrTri        (            magSqr(p0 - p1) +            magSqr(p2 - p0)        );        const scalar Atri =            0.5 *            (                (p1.x() - p0.x()) * (p2.y() - p0.y()) -                (p2.x() - p0.x()) * (p1.y() - p0.y())            );        const scalar stab = sqrt(sqr(Atri) + K);        const scalar Astab = Foam::max(ROOTVSMALL, 0.5 * (Atri + stab));        const vector gradAtri        (            0.5 * (p1.y() - p2.y()),            0.5 * (p2.x() - p1.x()),            0.0        );        const vector gradAstab = 0.5 * (gradAtri + Atri * gradAtri / stab);        const tensor gradGradAstab =            0.5 *            (                (gradAtri * gradAtri) / stab -                sqr(Atri) * (gradAtri * gradAtri) / pow(stab, 3)            );        const vector gradLt(4.0 * p0 - 2.0 * p1 - 2.0 * p2);        //- calculate the gradient        const scalar sqrAstab = sqr(Astab);        gradF += gradLt / Astab - (LSqrTri * gradAstab) / sqrAstab;        //- calculate the second gradient        gradGradF +=            gradGradLt / Astab -            twoSymm(gradLt * gradAstab) / sqrAstab -            gradGradAstab * LSqrTri / sqrAstab +            2.0 * LSqrTri * (gradAstab * gradAstab) / (sqrAstab * Astab);    }    //- stabilise diagonal terms    if( mag(gradGradF.xx()) < VSMALL ) gradGradF.xx() = VSMALL;    if( mag(gradGradF.yy()) < VSMALL ) gradGradF.yy() = VSMALL;}

void PDRkEpsilon::correct(){    if (!turbulence_)    {        // Re-calculate viscosity        nut_ = Cmu_*sqr(k_)/epsilon_;        nut_.correctBoundaryConditions();        // Re-calculate thermal diffusivity        //***HGWalphat_ = mut_/Prt_;        // alphat_.correctBoundaryConditions();        return;    }    RASModel::correct();    volScalarField divU(fvc::div(phi_/fvc::interpolate(rho_)));    if (mesh_.moving())    {        divU += fvc::div(mesh_.phi());    }    tmp<volTensorField> tgradU = fvc::grad(U_);    volScalarField G(GName(), rho_*nut_*(tgradU() && dev(twoSymm(tgradU()))));    tgradU.clear();    // Update espsilon and G at the wall    epsilon_.boundaryFieldRef().updateCoeffs();    // Add the blockage generation term so that it is included consistently    // in both the k and epsilon equations    const volScalarField& betav =        U_.db().lookupObject<volScalarField>("betav");    const volScalarField& Lobs =        U_.db().lookupObject<volScalarField>("Lobs");    const PDRDragModel& drag =        U_.db().lookupObject<PDRDragModel>("PDRDragModel");    volScalarField GR(drag.Gk());    volScalarField LI        (C4_*(Lobs + dimensionedScalar(dimLength, rootVSmall)));    // Dissipation equation    tmp<fvScalarMatrix> epsEqn    (        betav*fvm::ddt(rho_, epsilon_)      + fvm::div(phi_, epsilon_)      - fvm::laplacian(rho_*DepsilonEff(), epsilon_)     ==        C1_*betav*G*epsilon_/k_      + 1.5*pow(Cmu_, 3.0/4.0)*GR*sqrt(k_)/LI      - fvm::SuSp(((2.0/3.0)*C1_)*betav*rho_*divU, epsilon_)      - fvm::Sp(C2_*betav*rho_*epsilon_/k_, epsilon_)    );    epsEqn.ref().relax();    epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());    solve(epsEqn);    bound(epsilon_, epsilonMin_);    // Turbulent kinetic energy equation    tmp<fvScalarMatrix> kEqn    (        betav*fvm::ddt(rho_, k_)      + fvm::div(phi_, k_)      - fvm::laplacian(rho_*DkEff(), k_)     ==        betav*G + GR      - fvm::SuSp((2.0/3.0)*betav*rho_*divU, k_)      - fvm::Sp(betav*rho_*epsilon_/k_, k_)    );    kEqn.ref().relax();    solve(kEqn);    bound(k_, kMin_);    // Re-calculate viscosity    nut_ = Cmu_*sqr(k_)/epsilon_;    nut_.correctBoundaryConditions();    // Re-calculate thermal diffusivity    //***HGWalphat_ = mut_/Prt_;    // alphat_.correctBoundaryConditions();}

tmp<volSymmTensorField> laminar::devBeff() const{    return -nu()*dev(twoSymm(fvc::grad(U())));}

void LRR<BasicTurbulenceModel>::correct(){    if (!this->turbulence_)    {        return;    }    // Local references    const alphaField& alpha = this->alpha_;    const rhoField& rho = this->rho_;    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;    const volVectorField& U = this->U_;    volSymmTensorField& R = this->R_;    ReynoldsStress<RASModel<BasicTurbulenceModel> >::correct();    tmp<volTensorField> tgradU(fvc::grad(U));    const volTensorField& gradU = tgradU();    volSymmTensorField P(-twoSymm(R & gradU));    volScalarField G(this->GName(), 0.5*mag(tr(P)));    // Update epsilon and G at the wall    epsilon_.boundaryField().updateCoeffs();    // Dissipation equation    tmp<fvScalarMatrix> epsEqn    (        fvm::ddt(alpha, rho, epsilon_)      + fvm::div(alphaRhoPhi, epsilon_)      - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)     ==        Ceps1_*alpha*rho*G*epsilon_/k_      - fvm::Sp(Ceps2_*alpha*rho*epsilon_/k_, epsilon_)    );    epsEqn().relax();    epsEqn().boundaryManipulate(epsilon_.boundaryField());    solve(epsEqn);    bound(epsilon_, this->epsilonMin_);    // Correct the trace of the tensorial production to be consistent    // with the near-wall generation from the wall-functions    const fvPatchList& patches = this->mesh_.boundary();    forAll(patches, patchi)    {        const fvPatch& curPatch = patches[patchi];        if (isA<wallFvPatch>(curPatch))        {            forAll(curPatch, facei)            {                label faceCelli = curPatch.faceCells()[facei];                P[faceCelli] *= min                (                    G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL),                    1.0                );            }        }    }

tmp<volSymmTensorField> GenEddyVisc::B() const{    return ((2.0/3.0)*I)*k() - (muSgs_/rho())*dev(twoSymm(fvc::grad(U())));}

void kEpsilon::correct(){    // Bound in case of topological change    // HJ, 22/Aug/2007    if (mesh_.changing())    {        bound(k_, k0_);        bound(epsilon_, epsilon0_);    }    if (!turbulence_)    {        // Re-calculate viscosity        mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);        mut_.correctBoundaryConditions();        // Re-calculate thermal diffusivity        alphat_ = mut_/Prt_;        alphat_.correctBoundaryConditions();        return;    }    RASModel::correct();    volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));    if (mesh_.moving())    {        divU += fvc::div(mesh_.phi());    }    tmp<volTensorField> tgradU = fvc::grad(U_);    volScalarField G("RASModel::G", mut_*(tgradU() && dev(twoSymm(tgradU()))));    tgradU.clear();    // Update epsilon and G at the wall    epsilon_.boundaryField().updateCoeffs();    // Dissipation equation    tmp<fvScalarMatrix> epsEqn    (        fvm::ddt(rho_, epsilon_)      + fvm::div(phi_, epsilon_)      - fvm::laplacian(DepsilonEff(), epsilon_)     ==        C1_*G*epsilon_/k_      - fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_)      - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)    );    epsEqn().relax();    epsEqn().boundaryManipulate(epsilon_.boundaryField());    solve(epsEqn);    bound(epsilon_, epsilon0_);    // Turbulent kinetic energy equation    tmp<fvScalarMatrix> kEqn    (        fvm::ddt(rho_, k_)      + fvm::div(phi_, k_)      - fvm::laplacian(DkEff(), k_)     ==        G      - fvm::SuSp((2.0/3.0)*rho_*divU, k_)      - fvm::Sp(rho_*epsilon_/k_, k_)    );    kEqn().relax();    solve(kEqn);    bound(k_, k0_);    // Re-calculate viscosity    mut_ = rho_*Cmu_*sqr(k_)/epsilon_;    mut_.correctBoundaryConditions();    // Re-calculate thermal diffusivity    alphat_ = mut_/Prt_;    alphat_.correctBoundaryConditions();}

void LaunderSharmaKE::correct(){    // Bound in case of topological change    // HJ, 22/Aug/2007    if (mesh_.changing())    {        bound(k_, k0_);        bound(epsilon_, epsilon0_);    }    if (!turbulence_)    {        // Re-calculate viscosity        mut_ == rho_*Cmu_*fMu()*sqr(k_)/(epsilon_ + epsilonSmall_);        mut_ == min(mut_, muRatio()*mu());        // Re-calculate thermal diffusivity        alphat_ = mut_/Prt_;        alphat_.correctBoundaryConditions();        return;    }    RASModel::correct();    // Calculate parameters and coefficients for Launder-Sharma low-Reynolds    // number model    volScalarField E = 2.0*mu()*mut_*fvc::magSqrGradGrad(U_)/rho_;    volScalarField D = 2.0*mu()*magSqr(fvc::grad(sqrt(k_)))/rho_;    volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));    if (mesh_.moving())    {        divU += fvc::div(mesh_.phi());    }    tmp<volTensorField> tgradU = fvc::grad(U_);    volScalarField G("RASModel::G", mut_*(tgradU() && dev(twoSymm(tgradU()))));    tgradU.clear();    // Dissipation equation    tmp<fvScalarMatrix> epsEqn    (        fvm::ddt(rho_, epsilon_)      + fvm::div(phi_, epsilon_)      - fvm::laplacian(DepsilonEff(), epsilon_)     ==        C1_*G*epsilon_/k_ + fvm::SuSp((C3_ - 2.0/3.0*C1_)*rho_*divU, epsilon_)      - fvm::Sp(C2_*f2()*rho_*epsilon_/k_, epsilon_)    //+ 0.75*1.5*flameKproduction*epsilon_/k_      + E    );    epsEqn().relax();    solve(epsEqn);    bound(epsilon_, epsilon0_);    // Turbulent kinetic energy equation    tmp<fvScalarMatrix> kEqn    (        fvm::ddt(rho_, k_)      + fvm::div(phi_, k_)      - fvm::laplacian(DkEff(), k_)     ==        G - fvm::SuSp(2.0/3.0*rho_*divU, k_)      - fvm::Sp(rho_*(epsilon_ + D)/k_, k_)    //+ flameKproduction    );    kEqn().relax();    solve(kEqn);    bound(k_, k0_);    // Re-calculate viscosity    mut_ == Cmu_*fMu()*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);    mut_ == min(mut_, muRatio()*mu());    // Re-calculate thermal diffusivity    alphat_ = mut_/Prt_;    alphat_.correctBoundaryConditions();}

// Berechnung des ExtraspannungstensorsFoam::tmp<Foam::volSymmTensorField> Foam::power::devRhoReff() const{    if (obr_.foundObject<compressible::RASModel>("RASProperties"))    {        const compressible::RASModel& ras            = obr_.lookupObject<compressible::RASModel>("RASProperties");        return ras.devRhoReff();    }    else if (obr_.foundObject<incompressible::RASModel>("RASProperties"))    {        const incompressible::RASModel& ras            = obr_.lookupObject<incompressible::RASModel>("RASProperties");        return rho()*ras.devReff();        Info << "XXXXXXXXXXXXX"<< nl;    }    else if (obr_.foundObject<compressible::LESModel>("LESProperties"))    {        const compressible::LESModel& les =        obr_.lookupObject<compressible::LESModel>("LESProperties");        return les.devRhoReff();    }    else if (obr_.foundObject<incompressible::LESModel>("LESProperties"))    {        const incompressible::LESModel& les            = obr_.lookupObject<incompressible::LESModel>("LESProperties");        return rho()*les.devReff();    }    else if (obr_.foundObject<fluidThermo>("thermophysicalProperties"))    {        const fluidThermo& thermo =             obr_.lookupObject<fluidThermo>("thermophysicalProperties");        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);        return -thermo.mu()*dev(twoSymm(fvc::grad(U)));    }    else if    (        obr_.foundObject<transportModel>("transportProperties")    )    {        const transportModel& laminarT =            obr_.lookupObject<transportModel>("transportProperties");        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);        return -rho()*laminarT.nu()*dev(twoSymm(fvc::grad(U)));    }    else if (obr_.foundObject<dictionary>("transportProperties"))    {        const dictionary& transportProperties =             obr_.lookupObject<dictionary>("transportProperties");        dimensionedScalar nu(transportProperties.lookup("nu"));        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);        return -rho()*nu*dev(twoSymm(fvc::grad(U)));    }    else    {        FatalErrorIn("power::devRhoReff()")            << "No valid model for viscous stress calculation"            << exit(FatalError);        return volSymmTensorField::null();    }}

    // Reynolds stress equation    tmp<fvSymmTensorMatrix> REqn    (        fvm::ddt(alpha, rho, R)      + fvm::div(alphaRhoPhi, R)      - fvm::laplacian(alpha*rho*DREff(), R)      + fvm::Sp(((C1_/2)*epsilon_ + (C1s_/2)*G)*alpha*rho/k_, R)     ==        alpha*rho*P      - ((1.0/3.0)*I)*(((2.0 - C1_)*epsilon_ - C1s_*G)*alpha*rho)      + (C2_*(alpha*rho*epsilon_))*dev(innerSqr(b))      + alpha*rho*k_       *(            (C3_ - C3s_*mag(b))*dev(S)          + C4_*dev(twoSymm(b&S))          + C5_*twoSymm(b&Omega)        )    );    REqn().relax();    solve(REqn);    this->boundNormalStress(R);    k_ = 0.5*tr(R);    correctNut();    // Correct wall shear-stresses when applying wall-functions    this->correctWallShearStress(R);

void mixtureKEpsilon<BasicTurbulenceModel>::correct(){    const transportModel& gas = this->transport();    const twoPhaseSystem& fluid = gas.fluid();    // Only solve the mixture turbulence for the gas-phase    if (&gas != &fluid.phase1())    {        // This is the liquid phase but check the model for the gas-phase        // is consistent        this->liquidTurbulence();        return;    }    if (!this->turbulence_)    {        return;    }    // Initialise the mixture fields if they have not yet been constructed    initMixtureFields();    // Local references to gas-phase properties    const surfaceScalarField& phig = this->phi_;    const volVectorField& Ug = this->U_;    const volScalarField& alphag = this->alpha_;    volScalarField& kg = this->k_;    volScalarField& epsilong = this->epsilon_;    volScalarField& nutg = this->nut_;    // Local references to liquid-phase properties    mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence =        this->liquidTurbulence();    const surfaceScalarField& phil = liquidTurbulence.phi_;    const volVectorField& Ul = liquidTurbulence.U_;    const volScalarField& alphal = liquidTurbulence.alpha_;    volScalarField& kl = liquidTurbulence.k_;    volScalarField& epsilonl = liquidTurbulence.epsilon_;    volScalarField& nutl = liquidTurbulence.nut_;    // Local references to mixture properties    volScalarField& rhom = rhom_();    volScalarField& km = km_();    volScalarField& epsilonm = epsilonm_();    eddyViscosity<RASModel<BasicTurbulenceModel> >::correct();    // Update the effective mixture density    rhom = this->rhom();    // Mixture flux    surfaceScalarField phim("phim", mixFlux(phil, phig));    // Mixture velocity divergence    volScalarField divUm    (        mixU        (            fvc::div(fvc::absolute(phil, Ul)),            fvc::div(fvc::absolute(phig, Ug))        )    );    tmp<volScalarField> Gc;    {        tmp<volTensorField> tgradUl = fvc::grad(Ul);        Gc = tmp<volScalarField>        (            new volScalarField            (                this->GName(),                nutl*(tgradUl() && dev(twoSymm(tgradUl())))            )        );        tgradUl.clear();        // Update k, epsilon and G at the wall        kl.boundaryField().updateCoeffs();        epsilonl.boundaryField().updateCoeffs();        Gc().checkOut();    }    tmp<volScalarField> Gd;    {        tmp<volTensorField> tgradUg = fvc::grad(Ug);        Gd = tmp<volScalarField>        (            new volScalarField            (                this->GName(),                nutg*(tgradUg() && dev(twoSymm(tgradUg())))            )        );        tgradUg.clear();        // Update k, epsilon and G at the wall        kg.boundaryField().updateCoeffs();        epsilong.boundaryField().updateCoeffs();//.........这里部分代码省略.........

void tatunRNG::correct(){    if (!turbulence_)    {        // Re-calculate viscosity        mut_ = rho_*Cmu_*sqr(k_)/epsilon_;        mut_.correctBoundaryConditions();        // Re-calculate thermal diffusivity        alphat_ = mut_/Prt_;        alphat_.correctBoundaryConditions();        return;    }    RASModel::correct();    volScalarField divU(fvc::div(phi_/fvc::interpolate(rho_)));    if (mesh_.moving())    {        divU += fvc::div(mesh_.phi());    }    tmp<volTensorField> tgradU = fvc::grad(U_);    volScalarField S2((tgradU() && dev(twoSymm(tgradU()))));    tgradU.clear();    //volScalarField G(type() + ".G", mut_*S2);    volScalarField G(GName(), mut_*S2); // changed from 2.1.1 --> 2.2.2    volScalarField eta(sqrt(mag(S2))*k_/epsilon_);    volScalarField eta3(eta*sqr(eta));    volScalarField R    (        ((eta*(-eta/eta0_ + scalar(1)))/(beta_*eta3 + scalar(1)))    );    // Update epsilon and G at the wall    epsilon_.boundaryField().updateCoeffs();    // Dissipation equation    tmp<fvScalarMatrix> epsEqn    (        fvm::ddt(rho_, epsilon_)      + fvm::div(phi_, epsilon_)      - fvm::laplacian(DepsilonEff(), epsilon_)      ==        (C1_ - R)*G*epsilon_/k_      - fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_)      - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)    );    epsEqn().relax();    epsEqn().boundaryManipulate(epsilon_.boundaryField());    solve(epsEqn);    bound(epsilon_, epsilonMin_);    volScalarField YMperk = 2.0*rho_*epsilon_/(soundSpeed_*soundSpeed_);    // Turbulent kinetic energy equation    tmp<fvScalarMatrix> kEqn    (        fvm::ddt(rho_, k_)      + fvm::div(phi_, k_)      - fvm::laplacian(DkEff(), k_)      ==        G - fvm::SuSp(2.0/3.0*rho_*divU, k_)      - fvm::Sp(rho_*epsilon_/k_, k_)      - fvm::Sp(YMperk, k_)    );    kEqn().relax();    solve(kEqn);    bound(k_, kMin_);    // Re-calculate viscosity    mut_ = rho_*Cmu_*sqr(k_)/epsilon_;    mut_.correctBoundaryConditions();    // Re-calculate thermal diffusivity    alphat_ = mut_/Prt_;    alphat_.correctBoundaryConditions();}

void RNGkEpsilon<BasicTurbulenceModel>::correct(){    if (!this->turbulence_)    {        return;    }    // Local references    const alphaField& alpha = this->alpha_;    const rhoField& rho = this->rho_;    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;    const volVectorField& U = this->U_;    volScalarField& nut = this->nut_;    eddyViscosity<RASModel<BasicTurbulenceModel> >::correct();    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));    tmp<volTensorField> tgradU = fvc::grad(U);    volScalarField S2((tgradU() && dev(twoSymm(tgradU()))));    tgradU.clear();    volScalarField G(this->GName(), nut*S2);    volScalarField eta(sqrt(mag(S2))*k_/epsilon_);    volScalarField eta3(eta*sqr(eta));    volScalarField R    (        ((eta*(-eta/eta0_ + scalar(1)))/(beta_*eta3 + scalar(1)))    );    // Update epsilon and G at the wall    epsilon_.boundaryField().updateCoeffs();    // Dissipation equation    tmp<fvScalarMatrix> epsEqn    (        fvm::ddt(alpha, rho, epsilon_)      + fvm::div(alphaRhoPhi, epsilon_)      - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)     ==        (C1_ - R)*alpha*rho*G*epsilon_/k_      - fvm::SuSp(((2.0/3.0)*C1_ + C3_)*alpha*rho*divU, epsilon_)      - fvm::Sp(C2_*alpha*rho*epsilon_/k_, epsilon_)      + epsilonSource()    );    epsEqn().relax();    epsEqn().boundaryManipulate(epsilon_.boundaryField());    solve(epsEqn);    bound(epsilon_, this->epsilonMin_);    // Turbulent kinetic energy equation    tmp<fvScalarMatrix> kEqn    (        fvm::ddt(alpha, rho, k_)      + fvm::div(alphaRhoPhi, k_)      - fvm::laplacian(alpha*rho*DkEff(), k_)     ==        alpha*rho*G      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)      - fvm::Sp(alpha*rho*epsilon_/k_, k_)      + kSource()    );    kEqn().relax();    solve(kEqn);    bound(k_, this->kMin_);    correctNut();}