diff --git a/applications/solvers/dfHighSpeedFoam/calculateR.H b/applications/solvers/dfHighSpeedFoam/calculateR.H new file mode 100644 index 000000000..651b00ad8 --- /dev/null +++ b/applications/solvers/dfHighSpeedFoam/calculateR.H @@ -0,0 +1,49 @@ +start = std::clock(); + +combustion->correct(); + +label flag_mpi_init; +MPI_Initialized(&flag_mpi_init); +if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM); +end = std::clock(); +time_monitor_chem += double(end - start) / double(CLOCKS_PER_SEC); + +volScalarField Yt(0.0*Y[0]); + +start = std::clock(); +forAll(Y, i) +{ + volScalarField& Yi = Y[i]; + + if (i != inertIndex) + { + rhoYi[i].ref() += chemistry->RR(i)*runTime.deltaT(); + Info <<"max reaction rate "<< Yi.name() << " is " << max(chemistry->RR(i)).value() << endl; + + Yi = rhoYi[i]/rho; + Yi.max(0.0); + + Yi.correctBoundaryConditions(); + rhoYi[i] = rho*Yi; + Yt += Yi; + } +} + +Y[inertIndex] = scalar(1) - Yt; +Y[inertIndex].max(0.0); +rhoYi[inertIndex] = rho*Y[inertIndex]; + +end = std::clock(); +time_monitor_Y += double(end - start) / double(CLOCKS_PER_SEC); + +chemistry->correctThermo(); +Info<< "min/max(T) = " + << min(T).value() << ", " << max(T).value() << endl; + +p.ref() = rho()/psi(); +p.correctBoundaryConditions(); +rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField(); + + +Info << " finish calculate Reaction" << nl << endl; + diff --git a/applications/solvers/dfHighSpeedFoam/createFields.H b/applications/solvers/dfHighSpeedFoam/createFields.H index ddffc838e..e01dab723 100644 --- a/applications/solvers/dfHighSpeedFoam/createFields.H +++ b/applications/solvers/dfHighSpeedFoam/createFields.H @@ -1,5 +1,44 @@ #include "createRDeltaT.H" +word ddtSchemes("Euler"); +if (mesh.schemesDict().readIfPresent("timeScheme", ddtSchemes)) +{ + if(ddtSchemes == "RK2SSP" || ddtSchemes == "RK3SSP") + { + Info<< "ddtSchemes: " << ddtSchemes << endl; + } +} + +word chemScheme("RR"); +mesh.schemesDict().readIfPresent("chemScheme", chemScheme); + +if ((chemScheme == "wrate") || (chemScheme == "RR")) +{ + Info<< "chemScheme: " << chemScheme << endl; +} +else +{ + FatalErrorInFunction + << "chemScheme: " << chemScheme + << " is not a valid choice. " + << "Options are: wrate, RR" + << abort(FatalError); +} + + +if((ddtSchemes != "RK2SSP") && (ddtSchemes != "RK3SSP")) +{ + if(chemScheme == "wrate") + { + FatalErrorInFunction + << "This combination is not a valid choice. " + << "If you want to use wrate for chemistry, please use RK2SSP or RK3SSP scheme." + << abort(FatalError); + } +} + + + Info<< "Reading thermophysical properties\n" << endl; // psiThermo* pThermo = new hePsiThermo(mesh, word::null); diff --git a/applications/solvers/dfHighSpeedFoam/createFields_rk2.H b/applications/solvers/dfHighSpeedFoam/createFields_rk2.H new file mode 100644 index 000000000..25bdeffc3 --- /dev/null +++ b/applications/solvers/dfHighSpeedFoam/createFields_rk2.H @@ -0,0 +1,27 @@ +volScalarField rho_save("rho_save",thermo.rho()); + +volVectorField rhoU_save("rhoU_save",rho*U); + +volScalarField rhoE_save("rhoE_save",rho*(ea + 0.5*magSqr(U))); + +PtrList rhoYi_save(nspecies); +forAll(rhoYi_save,i) +{ + rhoYi_save.set + ( + i, + new volScalarField + ( + IOobject + ( + "rhoYi_save" + Y[i].name(), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + rho*Y[i] + ) + ); +} + diff --git a/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C b/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C index 67254ec4d..39c195cdb 100644 --- a/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C +++ b/applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C @@ -76,6 +76,7 @@ int main(int argc, char *argv[]) #include "createTime.H" #include "createDynamicFvMesh.H" #include "createFields.H" + #include "createFields_rk2.H" #include "createTimeControls.H" double time_monitor_flow=0; @@ -97,6 +98,27 @@ int main(int argc, char *argv[]) scalar CoNum = 0.0; scalar meanCoNum = 0.0; + std::vector rkcoe1(3); + std::vector rkcoe2(3); + std::vector rkcoe3(3); + scalar rk; + label nrk=0; + + if(ddtSchemes == "RK2SSP") + { + rkcoe1[0]=1.0; rkcoe2[0]=0.0; rkcoe3[0]=1.0; + rkcoe1[1]=0.5; rkcoe2[1]=0.5; rkcoe3[1]=0.5; + rkcoe1[2]=0.0; rkcoe2[2]=0.0; rkcoe3[2]=0.0; + rk=2; + } + else if(ddtSchemes == "RK3SSP") + { + rkcoe1[0]=1.0; rkcoe2[0]=0.0; rkcoe3[0]=1.0; + rkcoe1[1]=0.75; rkcoe2[1]=0.25; rkcoe3[1]=0.25; + rkcoe1[2]=0.33; rkcoe2[2]=0.66; rkcoe3[2]=0.66; + rk=3; + } + Info<< "\nStarting time loop\n" << endl; while (runTime.run()) @@ -127,197 +149,89 @@ int main(int argc, char *argv[]) } - // --- Directed interpolation of primitive fields onto faces - - surfaceScalarField rho_pos(interpolate(rho, pos)); - surfaceScalarField rho_neg(interpolate(rho, neg)); + volScalarField rho_rhs("rho_rhs",rho_save/runTime.deltaT()); + volVectorField rhoU_rhs("rhoU_rhs",rhoU_save/runTime.deltaT()); + volScalarField rhoYi_rhs("rhoYi_rhs",rhoYi_save[0]/runTime.deltaT()); + volScalarField rhoE_rhs("rhoE_rhs",rhoE_save/runTime.deltaT()); - PtrList rhoYi_pos(nspecies); - PtrList rhoYi_neg(nspecies); - forAll(rhoYi_pos,i) + if ((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) { - rhoYi_pos.set - ( - i, - new surfaceScalarField - ( - IOobject - ( - "rhoYi_pos" + Y[i].name(), - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - interpolate(rhoYi[i], pos,"Yi") - ) - ); - } + for( nrk=0 ; nrkmuEff()); - volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U)))); + Info<< "Time = " << runTime.timeName() << nl << endl; - // --- Solve density - #include "rhoEqn.H" + #include "phiCal.H" - start = std::clock(); - // --- Solve momentum - #include "rhoUEqn.H" - end = std::clock(); - time_monitor_flow += double(end - start) / double(CLOCKS_PER_SEC); + // --- Solve density + #include "rhoEqn.H" + start = std::clock(); + // --- Solve momentum + #include "rhoUEqn.H" + end = std::clock(); + time_monitor_flow += double(end - start) / double(CLOCKS_PER_SEC); - // --- Solve species - #include "rhoYEqn.H" + // --- Solve species + #include "rhoYEqn.H" - // --- Solve energy - #include "rhoEEqn.H" + // --- Solve energy + #include "rhoEEqn.H" + } turbulence->correct(); diff --git a/applications/solvers/dfHighSpeedFoam/phiCal.H b/applications/solvers/dfHighSpeedFoam/phiCal.H new file mode 100644 index 000000000..e5f569ed8 --- /dev/null +++ b/applications/solvers/dfHighSpeedFoam/phiCal.H @@ -0,0 +1,45 @@ +phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg; + +PtrList phiYi(nspecies); +forAll(phiYi,i) +{ + phiYi.set + ( + i, + new surfaceScalarField + ( + IOobject + ( + "phiYi_" + Y[i].name(), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + aphiv_pos*rhoYi_pos[i] + aphiv_neg*rhoYi_neg[i] + ) + ); +} + +surfaceVectorField phiUp +( + (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg) + + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf() +); + +surfaceScalarField phiEp +( + "phiEp", + aphiv_pos*(rho_pos*(ea_pos + 0.5*magSqr(U_pos)) + p_pos) + + aphiv_neg*(rho_neg*(ea_neg + 0.5*magSqr(U_neg)) + p_neg) + + aSf*p_pos - aSf*p_neg +); + +// Make flux for pressure-work absolute +if (mesh.moving()) +{ + phiEp += mesh.phi()*(a_pos*p_pos + a_neg*p_neg); +} + +volScalarField muEff("muEff", turbulence->muEff()); +volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U)))); diff --git a/applications/solvers/dfHighSpeedFoam/preCal.H b/applications/solvers/dfHighSpeedFoam/preCal.H new file mode 100644 index 000000000..879170583 --- /dev/null +++ b/applications/solvers/dfHighSpeedFoam/preCal.H @@ -0,0 +1,118 @@ +// --- Directed interpolation of primitive fields onto faces +surfaceScalarField rho_pos(interpolate(rho, pos)); +surfaceScalarField rho_neg(interpolate(rho, neg)); + +PtrList rhoYi_pos(nspecies); +PtrList rhoYi_neg(nspecies); +forAll(rhoYi_pos,i) +{ + rhoYi_pos.set + ( + i, + new surfaceScalarField + ( + IOobject + ( + "rhoYi_pos" + Y[i].name(), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + interpolate(rhoYi[i], pos,"Yi") + ) + ); +} + +forAll(rhoYi_neg,i) +{ + rhoYi_neg.set + ( + i, + new surfaceScalarField + ( + IOobject + ( + "rhoYi_neg" + Y[i].name(), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + interpolate(rhoYi[i], neg,"Yi") + ) + ); +} + +surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name())); +surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name())); + +volScalarField rPsi("rPsi", 1.0/psi); +surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name())); +surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name())); + +surfaceScalarField ea_pos(interpolate(ea, pos, T.name())); +surfaceScalarField ea_neg(interpolate(ea, neg, T.name())); + +surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos); +surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg); + +surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos); +surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg); + +surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf()); +surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf()); + +// Make fluxes relative to mesh-motion +if (mesh.moving()) +{ + phiv_pos -= mesh.phi(); + phiv_neg -= mesh.phi(); +} + +volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi)); +surfaceScalarField cSf_pos +( + "cSf_pos", + interpolate(c, pos, T.name())*mesh.magSf() +); +surfaceScalarField cSf_neg +( + "cSf_neg", + interpolate(c, neg, T.name())*mesh.magSf() +); + +surfaceScalarField ap +( + "ap", + max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero) +); +surfaceScalarField am +( + "am", + min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero) +); + +surfaceScalarField a_pos("a_pos", ap/(ap - am)); + +surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap))); + +surfaceScalarField aSf("aSf", am*a_pos); + +if (fluxScheme == "Tadmor") +{ + aSf = -0.5*amaxSf; + a_pos = 0.5; +} + +surfaceScalarField a_neg("a_neg", 1.0 - a_pos); + +phiv_pos *= a_pos; +phiv_neg *= a_neg; + +surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf); +surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf); + +// Reuse amaxSf for the maximum positive and negative fluxes +// estimated by the central scheme +amaxSf = max(mag(aphiv_pos), mag(aphiv_neg)); \ No newline at end of file diff --git a/applications/solvers/dfHighSpeedFoam/rhoEEqn.H b/applications/solvers/dfHighSpeedFoam/rhoEEqn.H index 532d064d1..2b96e8862 100644 --- a/applications/solvers/dfHighSpeedFoam/rhoEEqn.H +++ b/applications/solvers/dfHighSpeedFoam/rhoEEqn.H @@ -8,12 +8,23 @@ surfaceScalarField sigmaDotU & (a_pos*U_pos + a_neg*U_neg) ); -solve -( - fvm::ddt(rhoE) - + fvc::div(phiEp) - - fvc::div(sigmaDotU) -); +if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) +{ + rhoE_rhs = -fvc::div(phiEp)+fvc::div(sigmaDotU); + + rhoE = rkcoe1[nrk]*rhoE_save + + rkcoe2[nrk]*rhoE + + rkcoe3[nrk]*rhoE_rhs*runTime.deltaT(); +} +else +{ + solve + ( + fvm::ddt(rhoE) + + fvc::div(phiEp) + - fvc::div(sigmaDotU) + ); +} ea = rhoE/rho - 0.5*magSqr(U); ea.correctBoundaryConditions(); @@ -25,21 +36,39 @@ rhoE.boundaryFieldRef() == rho.boundaryField()*(ea.boundaryField() + 0.5*magSqr( if (!inviscid) { - fvScalarMatrix eEqn - ( - fvm::ddt(rho, ea) - fvc::ddt(rho, ea) - // alpha in deepflame is considered to calculate h by default (kappa/Cp), so multiply gamma to correct alpha - - fvm::laplacian(turbulence->alphaEff()*thermo.gamma(), ea) - + diffAlphaD - == - fvc::div(hDiffCorrFlux) - ); + if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) + { + rhoE_rhs = fvc::laplacian(turbulence->alphaEff()*thermo.gamma(), ea) + - diffAlphaD + + fvc::div(hDiffCorrFlux); + + rhoE += rkcoe3[nrk]*rhoE_rhs*runTime.deltaT(); + + ea = rhoE/rho - 0.5*magSqr(U); + ea.correctBoundaryConditions(); + + ha = ea + p/rho; + chemistry->correctThermo(); // before this, we must update ha = e + p/rho + } + else + { + fvScalarMatrix eEqn + ( + fvm::ddt(rho, ea) - fvc::ddt(rho, ea) + // alpha in deepflame is considered to calculate h by default (kappa/Cp), so multiply gamma to correct alpha + - fvm::laplacian(turbulence->alphaEff()*thermo.gamma(), ea) + + diffAlphaD + == + fvc::div(hDiffCorrFlux) + ); - eEqn.solve("ea"); + eEqn.solve("ea"); - ha = ea + p/rho; - chemistry->correctThermo(); - rhoE = rho*(ea + 0.5*magSqr(U)); + ha = ea + p/rho; + chemistry->correctThermo(); + rhoE = rho*(ea + 0.5*magSqr(U)); + } + } Info<< "min/max(T) = " diff --git a/applications/solvers/dfHighSpeedFoam/rhoEqn.H b/applications/solvers/dfHighSpeedFoam/rhoEqn.H index 21be7d5be..d02a21314 100644 --- a/applications/solvers/dfHighSpeedFoam/rhoEqn.H +++ b/applications/solvers/dfHighSpeedFoam/rhoEqn.H @@ -1,5 +1,20 @@ -solve -( - fvm::ddt(rho) - + fvc::div(phi) -); +if ((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) +{ + rho_rhs = -fvc::div(phi); + + rho = rkcoe1[nrk]*rho_save + + rkcoe2[nrk]*rho + + rkcoe3[nrk]*rho_rhs*runTime.deltaT(); + + Info <<"in rk"<< nrk+1 << " finish calculate rho" << nl << endl; +} +else +{ + solve + ( + fvm::ddt(rho) + + fvc::div(phi) + ); +} + + diff --git a/applications/solvers/dfHighSpeedFoam/rhoUEqn.H b/applications/solvers/dfHighSpeedFoam/rhoUEqn.H index f00a6899c..df6903eb9 100644 --- a/applications/solvers/dfHighSpeedFoam/rhoUEqn.H +++ b/applications/solvers/dfHighSpeedFoam/rhoUEqn.H @@ -1,23 +1,50 @@ -//same with rhoReactingFoam -solve -( - fvm::ddt(rhoU) - + fvc::div(phiUp) -); +if ((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) +{ + rhoU_rhs = -fvc::div(phiUp); -U.ref() = rhoU()/rho(); -U.correctBoundaryConditions(); -rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); + rhoU = rkcoe1[nrk]*rhoU_save + + rkcoe2[nrk]*rhoU + + rkcoe3[nrk]*rhoU_rhs*runTime.deltaT(); + + U.ref() = rhoU()/rho(); + U.correctBoundaryConditions(); + rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); -if (!inviscid) + if (!inviscid) + { + rhoU_rhs = fvc::laplacian(turbulence->muEff(),U) + fvc::div(tauMC); + rhoU += rkcoe3[nrk]*rhoU_rhs*runTime.deltaT(); + + U.ref() = rhoU()/rho(); + U.correctBoundaryConditions(); + rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); + } + + Info <<"in rk"<muEff(), U) - - fvc::div(tauMC) + fvm::ddt(rhoU) + + fvc::div(phiUp) ); - rhoU = rho*U; + U.ref() = rhoU()/rho(); + U.correctBoundaryConditions(); + rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); + + if (!inviscid) + { + solve + ( + fvm::ddt(rho, U) - fvc::ddt(rho, U) + - fvm::laplacian(turbulence->muEff(), U) + - fvc::div(tauMC) + ); + rhoU = rho*U; + } + + Info << "\nmin / max mag(U) ; " << gMin(mag(U)()()) << " / " << gMax(mag(U)()()) << nl << endl; } -Info << "\nmin / max mag(U) ; " << gMin(mag(U)()()) << " / " << gMax(mag(U)()()) << nl << endl; + diff --git a/applications/solvers/dfHighSpeedFoam/rhoYEqn.H b/applications/solvers/dfHighSpeedFoam/rhoYEqn.H index 46b0c9a94..26beb31c8 100644 --- a/applications/solvers/dfHighSpeedFoam/rhoYEqn.H +++ b/applications/solvers/dfHighSpeedFoam/rhoYEqn.H @@ -22,13 +22,33 @@ tmp> mvConvection ); { - start = std::clock(); - combustion->correct(); - label flag_mpi_init; - MPI_Initialized(&flag_mpi_init); - if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM); - end = std::clock(); - time_monitor_chem += double(end - start) / double(CLOCKS_PER_SEC); + if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) + { + if(chemScheme == "wrate") + { + start = std::clock(); + + chemistry->calculateW(); + + label flag_mpi_init; + MPI_Initialized(&flag_mpi_init); + if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM); + end = std::clock(); + time_monitor_chem += double(end - start) / double(CLOCKS_PER_SEC); + } + } + else + { + start = std::clock(); + + combustion->correct(); + + label flag_mpi_init; + MPI_Initialized(&flag_mpi_init); + if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM); + end = std::clock(); + time_monitor_chem += double(end - start) / double(CLOCKS_PER_SEC); + } volScalarField Yt(0.0*Y[0]); @@ -44,23 +64,41 @@ tmp> mvConvection } if (i != inertIndex) - { - // original convection term - // fvScalarMatrix YiEqn - // ( - // fvm::ddt(rho, Yi) - // + mvConvection->fvmDiv(phi, Yi) - // == - // combustion->R(Yi) - // ); - - solve - ( - fvm::ddt(rhoYi[i]) - + fvc::div(phiYi[i]) - == - chemistry->RR(i) - ); + { + if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) + { + rhoYi_rhs = -fvc::div(phiYi[i]); + + if(chemScheme == "wrate") + { + rhoYi_rhs.ref() += chemistry->wrate(i); + Info <<"max reaction rate "<< Yi.name() << " is " << max(chemistry->wrate(i)).value() << endl; + } + + + rhoYi[i] = rkcoe1[nrk]*rhoYi_save[i] + + rkcoe2[nrk]*rhoYi[i] + + rkcoe3[nrk]*rhoYi_rhs*runTime.deltaT(); + } + else + { + // original convection term + // fvScalarMatrix YiEqn + // ( + // fvm::ddt(rho, Yi) + // + mvConvection->fvmDiv(phi, Yi) + // == + // combustion->R(Yi) + // ); + + solve + ( + fvm::ddt(rhoYi[i]) + + fvc::div(phiYi[i]) + == + chemistry->RR(i) + ); + } Yi=rhoYi[i]/rho; Yi.max(0.0); @@ -70,22 +108,34 @@ tmp> mvConvection const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf(); tmp DEff = chemistry->rhoD(i) + turbulence->mut()/Sct; - // original term - // YiEqn -= fvm::laplacian(DEff(), Yi) - mvConvection->fvmDiv(phiUc, Yi); + if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) + { + rhoYi_rhs = fvc::laplacian(DEff(), Yi) - fvc::div(phiUc,Yi,"div(phic,Yi)"); - fvScalarMatrix YiEqn - ( - fvm::ddt(rho, Yi) - fvc::ddt(rho, Yi) - - fvm::laplacian(DEff(), Yi) - + mvConvection->fvmDiv(phiUc, Yi) - ); + rhoYi[i] += rkcoe3[nrk]*rhoYi_rhs*runTime.deltaT(); + + Yi = rhoYi[i]/rho; + Yi.max(0.0); + } + else + { + // original term + // YiEqn -= fvm::laplacian(DEff(), Yi) - mvConvection->fvmDiv(phiUc, Yi); - YiEqn.relax(); + fvScalarMatrix YiEqn + ( + fvm::ddt(rho, Yi) - fvc::ddt(rho, Yi) + - fvm::laplacian(DEff(), Yi) + + mvConvection->fvmDiv(phiUc, Yi) + ); - YiEqn.solve("Yi"); + YiEqn.relax(); - Yi.max(0.0); + YiEqn.solve("Yi"); + Yi.max(0.0); + } + } // original term @@ -93,6 +143,7 @@ tmp> mvConvection // YiEqn.solve("Yi"); // Yi.max(0.0); + if((ddtSchemes == "RK2SSP") || (ddtSchemes == "RK3SSP")) Yi.correctBoundaryConditions(); rhoYi[i] = rho*Yi; Yt += Yi; } @@ -100,6 +151,7 @@ tmp> mvConvection Y[inertIndex] = scalar(1) - Yt; Y[inertIndex].max(0.0); + rhoYi[inertIndex] = rho*Y[inertIndex]; end = std::clock(); time_monitor_Y += double(end - start) / double(CLOCKS_PER_SEC); diff --git a/applications/solvers/dfHighSpeedFoam/updateFieldsSave.H b/applications/solvers/dfHighSpeedFoam/updateFieldsSave.H new file mode 100644 index 000000000..2f9e8c33a --- /dev/null +++ b/applications/solvers/dfHighSpeedFoam/updateFieldsSave.H @@ -0,0 +1,12 @@ +rho_save = rho; + +rhoU_save = rho*U; + +rhoE_save = rho*(ea + 0.5*magSqr(U)); + +forAll(rhoYi_save,i) +{ + rhoYi_save[i] = rho*Y[i]; +} + +Info <<"finish update saved fields"<< endl; \ No newline at end of file