sprayFoam计算等压环境下点火延迟时间
-
目前在尝试使用sprayFoam等压环境下的C12H26点火延迟时间，由于sprayFoam中通过泊松方程求解压力（见原始pEqn.H），网格内监测压力显示不能保持等压状态。因此参考chemFoam的源代码，想把泊松方程改为完全气体状态方程P=rhoRT(见修改的pEqn.H)，同时在sprayFoam.C中注释掉压力修正的过程（见修改的sprayFoam.C）。另外能量方程也改为求解Enthalpy的形式。在OpenFOAM6环境下编译没问题，但是计算两个时间步后即浮点溢出了，计算设置见下楼的图1，错误提示下楼的图2。请问各位前辈有没有解决方法啊？
原始pEqn.Hrho = thermo.rho(); rho = max(rho, rhoMin); rho = min(rho, rhoMax); rho.relax(); volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); if (pimple.nCorrPISO() <= 1) { tUEqn.clear(); } if (pimple.transonic()) { surfaceScalarField phid ( "phid", fvc::interpolate(psi) *( fvc::flux(HbyA) + MRF.zeroFilter ( rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho) ) ) ); MRF.makeRelative(fvc::interpolate(psi), phid); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvm::div(phid, p) - fvm::laplacian(rhorAUf, p) == parcels.Srho() + fvOptions(psi, p, rho.name()) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { phi == pEqn.flux(); } } } else { surfaceScalarField phiHbyA ( "phiHbyA", ( fvc::flux(rho*HbyA) + rhorAUf*fvc::ddtCorr(rho, U, phi) ) ); MRF.makeRelative(fvc::interpolate(rho), phiHbyA); // Update the pressure BCs to ensure flux consistency constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); while (pimple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvc::div(phiHbyA) - fvm::laplacian(rhorAUf, p) == parcels.Srho() + fvOptions(psi, p, rho.name()) ); pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { phi = phiHbyA + pEqn.flux(); } } } #include "rhoEqn.H" #include "compressibleContinuityErrs.H" // Explicitly relax pressure for momentum corrector p.relax(); // Recalculate density from the relaxed pressure rho = thermo.rho(); rho = max(rho, rhoMin); rho = min(rho, rhoMax); rho.relax(); Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; U = HbyA - rAU*fvc::grad(p); U.correctBoundaryConditions(); fvOptions.correct(U); K = 0.5*magSqr(U); if (thermo.dpdt()) { dpdt = fvc::ddt(p); }修改的pEqn.H
{ rho = thermo.rho(); p = rho*Rspecific*thermo.T(); }修改的sprayFoam.C
/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. Application sprayFoam Description Transient solver for compressible, turbulent flow with a spray particle cloud. \*---------------------------------------------------------------------------*/ #include "fvCFD.H" #include "turbulentFluidThermoModel.H" #include "basicSprayCloud.H" #include "psiReactionThermo.H" #include "CombustionModel.H" #include "radiationModel.H" #include "SLGThermo.H" #include "pimpleControl.H" #include "fvOptions.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "postProcess.H" #include "setRootCaseLists.H" #include "createTime.H" #include "createMesh.H" #include "createControl.H" #include "createTimeControls.H" #include "createFields.H" #include "createFieldRefs.H" #include "compressibleCourantNo.H" #include "setInitialDeltaT.H" turbulence->validate(); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (runTime.run()) { #include "readTimeControls.H" #include "compressibleCourantNo.H" #include "setDeltaT.H" runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; parcels.evolve(); #include "rhoEqn.H" #include "YEqn.H" #include "EEqn.H" #include "pEqn.H" /* if (!pimple.frozenFlow()) { #include "rhoEqn.H" // --- Pressure-velocity PIMPLE corrector loop while (pimple.loop()) { // #include "UEqn.H" #include "YEqn.H" #include "EEqn.H" // --- Pressure corrector loop while (pimple.correct()) { #include "pEqn.H" } if (pimple.turbCorr()) { turbulence->correct(); } } rho = thermo.rho(); if (runTime.write()) { combustion->Qdot()().write(); } } else { if (runTime.writeTime()) { parcels.write(); } } */ Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return 0; } // ************************************************************************* // -
-
游荡
