Go to the documentation of this file.00001 rho = thermo.rho();
00002 rho = max(rho, rhoMin);
00003 rho = min(rho, rhoMax);
00004 rho.relax();
00005
00006 volScalarField rUA = 1.0/UEqn().A();
00007 U = rUA*UEqn().H();
00008 UEqn.clear();
00009
00010 bool closedVolume = false;
00011
00012 if (transonic)
00013 {
00014 surfaceScalarField phid
00015 (
00016 "phid",
00017 fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
00018 );
00019
00020 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
00021 {
00022 fvScalarMatrix pEqn
00023 (
00024 fvm::div(phid, p)
00025 - fvm::laplacian(rho*rUA, p)
00026 );
00027
00028
00029 pEqn.relax(mesh.relaxationFactor("pEqn"));
00030
00031 pEqn.setReference(pRefCell, pRefValue);
00032
00033
00034 if (nonOrth == 0)
00035 {
00036 eqnResidual = pEqn.solve().initialResidual();
00037 maxResidual = max(eqnResidual, maxResidual);
00038 }
00039 else
00040 {
00041 pEqn.solve();
00042 }
00043
00044 if (nonOrth == nNonOrthCorr)
00045 {
00046 phi == pEqn.flux();
00047 }
00048 }
00049 }
00050 else
00051 {
00052 phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
00053 closedVolume = adjustPhi(phi, U, p);
00054
00055 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
00056 {
00057 fvScalarMatrix pEqn
00058 (
00059 fvm::laplacian(rho*rUA, p) == fvc::div(phi)
00060 );
00061
00062 pEqn.setReference(pRefCell, pRefValue);
00063
00064
00065 if (nonOrth == 0)
00066 {
00067 eqnResidual = pEqn.solve().initialResidual();
00068 maxResidual = max(eqnResidual, maxResidual);
00069 }
00070 else
00071 {
00072 pEqn.solve();
00073 }
00074
00075 if (nonOrth == nNonOrthCorr)
00076 {
00077 phi -= pEqn.flux();
00078 }
00079 }
00080 }
00081
00082
00083 #include <finiteVolume/continuityErrs.H>
00084
00085
00086 p.relax();
00087
00088 U -= rUA*fvc::grad(p);
00089 U.correctBoundaryConditions();
00090
00091
00092
00093 if (closedVolume)
00094 {
00095 p += (initialMass - fvc::domainIntegrate(psi*p))
00096 /fvc::domainIntegrate(psi);
00097 }
00098
00099 rho = thermo.rho();
00100 rho = max(rho, rhoMin);
00101 rho = min(rho, rhoMax);
00102 rho.relax();
00103 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;