songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge commit 'OpenCFD/master' into olesenm

Browse Source
This commit is contained in:
Mark Olesen 2009-02-06 10:40:19 +01:00
commit d3ec38f0aa
55 changed files with 1479 additions and 630 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
View File

@ -2,7 +2,7 @@
volScalarField kappaEff
(
"kappaEff",
turbulence->nu() + turbulence->nut()/Pr
turbulence->nu() + turbulence->nut()/Prt
);
fvScalarMatrix TEqn

11
applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
View File

@ -1,4 +1,4 @@
// Solve the Momentum equation
// Solve the momentum equation
tmp<fvVectorMatrix> UEqn
(
@ -13,8 +13,13 @@
(
UEqn()
==
-fvc::grad(pd)
+ beta*gh*fvc::grad(T)
-fvc::reconstruct
(
(
fvc::snGrad(pd)
- betaghf*fvc::snGrad(T)
) * mesh.magSf()
)
).initialResidual();
maxResidual = max(eqnResidual, maxResidual);

4
applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
View File

@ -53,8 +53,8 @@
incompressible::RASModel::New(U, phi, laminarTransport)
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
Info<< "Calculating field beta*(g.h)\n" << endl;
surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf()));
label pdRefCell = 0;
scalar pdRefValue = 0.0;

76
applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pdEqn.H
View File

@ -1,41 +1,49 @@
volScalarField rUA = 1.0/UEqn().A();
U = rUA*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, pd);
phi += fvc::interpolate(beta*gh*rUA)*fvc::snGrad(T)*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
(
fvm::laplacian(rUA, pd) == fvc::div(phi)
);
volScalarField rUA("rUA", 1.0/UEqn().A());
surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
pdEqn.setReference(pdRefCell, pdRefValue);
U = rUA*UEqn().H();
UEqn.clear();
// retain the residual from the first iteration
if (nonOrth == 0)
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, pd);
surfaceScalarField buoyancyPhi = -betaghf*fvc::snGrad(T)*rUAf*mesh.magSf();
phi -= buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
eqnResidual = pdEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pdEqn.solve();
fvScalarMatrix pdEqn
(
fvm::laplacian(rUAf, pd) == fvc::div(phi)
);
pdEqn.setReference(pdRefCell, pdRefValue);
// retain the residual from the first iteration
if (nonOrth == 0)
{
eqnResidual = pdEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pdEqn.solve();
}
if (nonOrth == nNonOrthCorr)
{
// Calculate the conservative fluxes
phi -= pdEqn.flux();
// Explicitly relax pressure for momentum corrector
pd.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rUAf);
U.correctBoundaryConditions();
}
}
if (nonOrth == nNonOrthCorr)
{
phi -= pdEqn.flux();
}
#include "continuityErrs.H"
}
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
pd.relax();
U -= rUA*(fvc::grad(pd) - beta*gh*fvc::grad(T));
U.correctBoundaryConditions();

4
applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/readTransportProperties.H
View File

@ -9,5 +9,5 @@
// reference kinematic pressure [m2/s2]
dimensionedScalar pRef(laminarTransport.lookup("pRef"));
// Prandtl number
dimensionedScalar Pr(laminarTransport.lookup("Pr"));
// turbulent Prandtl number
dimensionedScalar Prt(laminarTransport.lookup("Prt"));

2
applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H
View File

@ -23,7 +23,7 @@
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pd + rhoEff*gh + pRef
pd + rhoEff*(g & mesh.C()) + pRef
);
p.write();
}

16
applications/solvers/heatTransfer/buoyantFoam/UEqn.H
View File

@ -9,4 +9,18 @@
UEqn.relax();
solve(UEqn == -fvc::grad(pd) - fvc::grad(rho)*gh);
if (momentumPredictor)
{
solve
(
UEqn
==
-fvc::reconstruct
(
(
fvc::snGrad(pd)
+ ghf*fvc::snGrad(rho)
) * mesh.magSf()
)
);
}

37
applications/solvers/heatTransfer/buoyantFoam/buoyantFoam.C
View File

@ -41,42 +41,41 @@ Description
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readEnvironmentalProperties.H"
#include "createFields.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "readEnvironmentalProperties.H"
# include "createFields.H"
# include "initContinuityErrs.H"
# include "readTimeControls.H"
# include "compressibleCourantNo.H"
# include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
# include "readTimeControls.H"
# include "readPISOControls.H"
# include "compressibleCourantNo.H"
# include "setDeltaT.H"
#include "readTimeControls.H"
#include "readPISOControls.H"
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "rhoEqn.H"
#include "rhoEqn.H"
# include "UEqn.H"
#include "UEqn.H"
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
# include "hEqn.H"
# include "pEqn.H"
#include "hEqn.H"
#include "pEqn.H"
}
turbulence->correct();

1
applications/solvers/heatTransfer/buoyantFoam/createFields.H
View File

@ -58,6 +58,7 @@
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());
dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));

103
applications/solvers/heatTransfer/buoyantFoam/pEqn.H
View File

@ -1,60 +1,67 @@
bool closedVolume = pd.needReference();
rho = thermo->rho();
volScalarField rUA = 1.0/UEqn.A();
U = rUA*UEqn.H();
phi =
fvc::interpolate(rho)
*(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
)
- fvc::interpolate(rho*rUA*gh)*fvc::snGrad(rho)*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
bool closedVolume = pd.needReference();
rho = thermo->rho();
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
U = rUA*UEqn.H();
surfaceScalarField phiU
(
fvm::ddt(psi, pd)
+ fvc::ddt(psi)*pRef
+ fvc::ddt(psi, rho)*gh
+ fvc::div(phi)
- fvm::laplacian(rho*rUA, pd)
fvc::interpolate(rho)
*(
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
)
);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
pdEqn.solve(mesh.solver(pd.name() + "Final"));
}
else
{
pdEqn.solve(mesh.solver(pd.name()));
fvScalarMatrix pdEqn
(
fvm::ddt(psi, pd)
+ fvc::ddt(psi)*pRef
+ fvc::ddt(psi, rho)*gh
+ fvc::div(phi)
- fvm::laplacian(rhorUAf, pd)
);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pdEqn.solve(mesh.solver(pd.name() + "Final"));
}
else
{
pdEqn.solve(mesh.solver(pd.name()));
}
if (nonOrth == nNonOrthCorr)
{
phi += pdEqn.flux();
}
}
if (nonOrth == nNonOrthCorr)
U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
U.correctBoundaryConditions();
p == pd + rho*gh + pRef;
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
phi += pdEqn.flux();
p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
/fvc::domainIntegrate(thermo->psi());
rho = thermo->rho();
}
}
p == pd + rho*gh + pRef;
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
U.correctBoundaryConditions();
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
/fvc::domainIntegrate(thermo->psi());
pd == p - (rho*gh + pRef);
rho = thermo->rho();
}

10
applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
View File

@ -11,7 +11,15 @@
eqnResidual = solve
(
UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh
UEqn()
==
-fvc::reconstruct
(
(
fvc::snGrad(pd)
+ ghf*fvc::snGrad(rho)
) * mesh.magSf()
)
).initialResidual();
maxResidual = max(eqnResidual, maxResidual);

1
applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
View File

@ -53,6 +53,7 @@
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());
dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));

98
applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
View File

@ -1,55 +1,65 @@
volScalarField rUA = 1.0/UEqn().A();
U = rUA*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p);
phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
(
fvm::laplacian(rho*rUA, pd) == fvc::div(phi)
);
volScalarField rUA = 1.0/UEqn().A();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
pdEqn.setReference(pdRefCell, pdRefValue);
// retain the residual from the first iteration
if (nonOrth == 0)
U = rUA*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p);
surfaceScalarField buoyancyPhi = ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
phi -= buoyancyPhi;
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
eqnResidual = pdEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pdEqn.solve();
fvScalarMatrix pdEqn
(
fvm::laplacian(rhorUAf, pd) == fvc::div(phi)
);
pdEqn.setReference(pdRefCell, pdRefValue);
// retain the residual from the first iteration
if (nonOrth == 0)
{
eqnResidual = pdEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pdEqn.solve();
}
if (nonOrth == nNonOrthCorr)
{
// Calculate the conservative fluxes
phi -= pdEqn.flux();
// Explicitly relax pressure for momentum corrector
pd.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rhorUAf);
U.correctBoundaryConditions();
}
}
if (nonOrth == nNonOrthCorr)
#include "continuityErrs.H"
p == pd + rho*gh + pRef;
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
phi -= pdEqn.flux();
p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
/fvc::domainIntegrate(thermo->psi());
}
}
#include "continuityErrs.H"
rho = thermo->rho();
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
// Explicitly relax pressure for momentum corrector
pd.relax();
p = pd + rho*gh + pRef;
U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
U.correctBoundaryConditions();
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
/fvc::domainIntegrate(thermo->psi());
pd == p - (rho*gh + pRef);
}
rho = thermo->rho();
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

1
applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/Make/options
View File

@ -1,4 +1,5 @@
EXE_INC = \
-I../buoyantSimpleFoam \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-I$(LIB_SRC)/turbulenceModels \

18
applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/UEqn.H
View File

@ -1,18 +0,0 @@
// Solve the Momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
- fvm::Sp(fvc::div(phi), U)
+ turbulence->divDevRhoReff(U)
);
UEqn().relax();
eqnResidual = solve
(
UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh
).initialResidual();
maxResidual = max(eqnResidual, maxResidual);

1
applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/createFields.H
View File

@ -54,6 +54,7 @@
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());
dimensionedScalar pRef("pRef", p.dimensions(), thermo->lookup("pRef"));

54
applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/pEqn.H
View File

@ -1,54 +0,0 @@
volScalarField rUA = 1.0/UEqn().A();
U = rUA*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
bool closedVolume = adjustPhi(phi, U, p);
phi -= fvc::interpolate(rho*gh*rUA)*fvc::snGrad(rho)*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
(
fvm::laplacian(rho*rUA, pd) == fvc::div(phi)
);
pdEqn.setReference(pdRefCell, pdRefValue);
// retain the residual from the first iteration
if (nonOrth == 0)
{
eqnResidual = pdEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pdEqn.solve();
}
if (nonOrth == nNonOrthCorr)
{
phi -= pdEqn.flux();
}
}
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
pd.relax();
p = pd + rho*gh + pRef;
U -= rUA*(fvc::grad(pd) + fvc::grad(rho)*gh);
U.correctBoundaryConditions();
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
/fvc::domainIntegrate(thermo->psi());
pd == p - (rho*gh + pRef);
}
rho = thermo->rho();
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

2
applications/solvers/heatTransfer/chtMultiRegionFoam/Make/files
View File

@ -7,6 +7,8 @@ derivedFvPatchFields/solidWallHeatFluxTemperatureCoupled/solidWallHeatFluxTemper
derivedFvPatchFields/solidWallTemperatureCoupled/solidWallTemperatureCoupledFvPatchScalarField.C
derivedFvPatchFields/solidWallMixedTemperatureCoupled/solidWallMixedTemperatureCoupledFvPatchScalarField.C
fluid/compressibleCourantNo.C
chtMultiRegionFoam.C
EXE = $(FOAM_APPBIN)/chtMultiRegionFoam

11
applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
View File

@ -36,16 +36,10 @@ Description
#include "turbulenceModel.H"
#include "fixedGradientFvPatchFields.H"
#include "regionProperties.H"
#include "compressibleCourantNo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "solveContinuityEquation.C"
#include "solveMomentumEquation.C"
#include "compressibleContinuityErrors.C"
#include "solvePressureDifferenceEquation.C"
#include "solveEnthalpyEquation.C"
#include "compressibleCourantNo.C"
int main(int argc, char *argv[])
{
@ -58,7 +52,6 @@ int main(int argc, char *argv[])
# include "createSolidMeshes.H"
# include "createFluidFields.H"
# include "createSolidFields.H"
# include "initContinuityErrs.H"
@ -89,6 +82,7 @@ int main(int argc, char *argv[])
{
Info<< "\nSolving for fluid region "
<< fluidRegions[i].name() << endl;
# include "setRegionFluidFields.H"
# include "readFluidMultiRegionPISOControls.H"
# include "solveFluid.H"
}
@ -97,6 +91,7 @@ int main(int argc, char *argv[])
{
Info<< "\nSolving for solid region "
<< solidRegions[i].name() << endl;
# include "setRegionSolidFields.H"
# include "readSolidMultiRegionPISOControls.H"
# include "solveSolid.H"
}

31
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
View File

@ -1,10 +1,25 @@
tmp<fvVectorMatrix> UEqn = solveMomentumEquation
// Solve the Momentum equation
tmp<fvVectorMatrix> UEqn
(
momentumPredictor,
Uf[i],
rhof[i],
phif[i],
pdf[i],
ghf[i],
turb[i]
fvm::ddt(rho, U)
+ fvm::div(phi, U)
+ turb.divDevRhoReff(U)
);
UEqn().relax();
if (momentumPredictor)
{
solve
(
UEqn()
==
-fvc::reconstruct
(
(
fvc::snGrad(pd)
+ ghf*fvc::snGrad(rho)
) * mesh.magSf()
)
);
}

58
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.C
View File

@ -1,58 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Continuity errors for fluid meshes
\*---------------------------------------------------------------------------*/
void compressibleContinuityErrors
(
scalar& cumulativeContErr,
const volScalarField& rho,
const basicThermo& thermo
)
{
dimensionedScalar totalMass = fvc::domainIntegrate(rho);
scalar sumLocalContErr =
(
fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass
).value();
scalar globalContErr =
(
fvc::domainIntegrate(rho - thermo.rho())/totalMass
).value();
cumulativeContErr += globalContErr;
const word& regionName = rho.mesh().name();
Info<< "time step continuity errors (" << regionName << ")"
<< ": sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr
<< endl;
}

21
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleContinuityErrors.H Normal file
View File

@ -0,0 +1,21 @@
{
dimensionedScalar totalMass = fvc::domainIntegrate(rho);
scalar sumLocalContErr =
(
fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass
).value();
scalar globalContErr =
(
fvc::domainIntegrate(rho - thermo.rho())/totalMass
).value();
cumulativeContErr[i] += globalContErr;
Info<< "time step continuity errors (" << mesh.name() << ")"
<< ": sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr[i]
<< endl;
}

9
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C
View File

@ -22,13 +22,12 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Calculates and outputs the mean and maximum Courant Numbers for the fluid
regions
\*---------------------------------------------------------------------------*/
scalar compressibleCourantNo
#include "compressibleCourantNo.H"
#include "fvc.H"
Foam::scalar Foam::compressibleCourantNo
(
const fvMesh& mesh,
const Time& runTime,

26
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveContinuityEquation.C → applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.H
View File

@ -23,15 +23,27 @@ License
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Solve continuity equation
Calculates and outputs the mean and maximum Courant Numbers for the fluid
regions
\*---------------------------------------------------------------------------*/
void solveContinuityEquation
(
volScalarField& rho,
const surfaceScalarField& phi
)
#ifndef compressibleCourantNo_H
#define compressibleCourantNo_H
#include "fvMesh.H"
namespace Foam
{
solve(fvm::ddt(rho) + fvc::div(phi));
scalar compressibleCourantNo
(
const fvMesh& mesh,
const Time& runTime,
const volScalarField& rho,
const surfaceScalarField& phi
);
}
#endif
// ************************************************************************* //

4
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H
View File

@ -7,8 +7,8 @@
(
fluidRegions[regionI],
runTime,
rhof[regionI],
phif[regionI]
rhoFluid[regionI],
phiFluid[regionI]
),
CoNum
);

91
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
View File

@ -1,15 +1,16 @@
// Initialise fluid field pointer lists
PtrList<basicThermo> thermof(fluidRegions.size());
PtrList<volScalarField> rhof(fluidRegions.size());
PtrList<volScalarField> Kf(fluidRegions.size());
PtrList<volVectorField> Uf(fluidRegions.size());
PtrList<surfaceScalarField> phif(fluidRegions.size());
PtrList<compressible::turbulenceModel> turb(fluidRegions.size());
PtrList<volScalarField> DpDtf(fluidRegions.size());
PtrList<volScalarField> ghf(fluidRegions.size());
PtrList<volScalarField> pdf(fluidRegions.size());
PtrList<basicThermo> thermoFluid(fluidRegions.size());
PtrList<volScalarField> rhoFluid(fluidRegions.size());
PtrList<volScalarField> KFluid(fluidRegions.size());
PtrList<volVectorField> UFluid(fluidRegions.size());
PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
PtrList<volScalarField> DpDtFluid(fluidRegions.size());
PtrList<volScalarField> ghFluid(fluidRegions.size());
PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
PtrList<volScalarField> pdFluid(fluidRegions.size());
List<scalar> initialMassf(fluidRegions.size());
List<scalar> initialMassFluid(fluidRegions.size());
dimensionedScalar pRef
(
@ -24,8 +25,8 @@
Info<< "*** Reading fluid mesh thermophysical properties for region "
<< fluidRegions[i].name() << nl << endl;
Info<< " Adding to pdf\n" << endl;
pdf.set
Info<< " Adding to pdFluid\n" << endl;
pdFluid.set
(
i,
new volScalarField
@ -42,16 +43,15 @@
)
);
Info<< " Adding to thermof\n" << endl;
thermof.set
Info<< " Adding to thermoFluid\n" << endl;
thermoFluid.set
(
i,
basicThermo::New(fluidRegions[i]).ptr()
);
Info<< " Adding to rhof\n" << endl;
rhof.set
Info<< " Adding to rhoFluid\n" << endl;
rhoFluid.set
(
i,
new volScalarField
@ -64,12 +64,12 @@
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
thermof[i].rho()
thermoFluid[i].rho()
)
);
Info<< " Adding to Kf\n" << endl;
Kf.set
Info<< " Adding to KFluid\n" << endl;
KFluid.set
(
i,
new volScalarField
@ -82,12 +82,12 @@
IOobject::NO_READ,
IOobject::NO_WRITE
),
thermof[i].Cp()*thermof[i].alpha()
thermoFluid[i].Cp()*thermoFluid[i].alpha()
)
);
Info<< " Adding to Uf\n" << endl;
Uf.set
Info<< " Adding to UFluid\n" << endl;
UFluid.set
(
i,
new volVectorField
@ -104,8 +104,8 @@
)
);
Info<< " Adding to phif\n" << endl;
phif.set
Info<< " Adding to phiFluid\n" << endl;
phiFluid.set
(
i,
new surfaceScalarField
@ -118,29 +118,29 @@
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(rhof[i]*Uf[i])
linearInterpolate(rhoFluid[i]*UFluid[i])
& fluidRegions[i].Sf()
)
);
Info<< " Adding to turb\n" << endl;
turb.set
Info<< " Adding to turbulence\n" << endl;
turbulence.set
(
i,
autoPtr<compressible::turbulenceModel>
(
compressible::turbulenceModel::New
(
rhof[i],
Uf[i],
phif[i],
thermof[i]
rhoFluid[i],
UFluid[i],
phiFluid[i],
thermoFluid[i]
)
).ptr()
);
Info<< " Adding to DpDtf\n" << endl;
DpDtf.set
Info<< " Adding to DpDtFluid\n" << endl;
DpDtFluid.set
(
i,
new volScalarField
@ -150,9 +150,9 @@
surfaceScalarField
(
"phiU",
phif[i]/fvc::interpolate(rhof[i])
phiFluid[i]/fvc::interpolate(rhoFluid[i])
),
thermof[i].p()
thermoFluid[i].p()
)
)
);
@ -162,8 +162,8 @@
("environmentalProperties");
dimensionedVector g(environmentalProperties.lookup("g"));
Info<< " Adding to ghf\n" << endl;
ghf.set
Info<< " Adding to ghFluid\n" << endl;
ghFluid.set
(
i,
new volScalarField
@ -172,12 +172,21 @@
g & fluidRegions[i].C()
)
);
ghfFluid.set
(
i,
new surfaceScalarField
(
"ghf",
g & fluidRegions[i].Cf()
)
);
Info<< " Updating p from pd\n" << endl;
thermof[i].p() == pdf[i] + rhof[i]*ghf[i] + pRef;
thermof[i].correct();
thermoFluid[i].p() == pdFluid[i] + rhoFluid[i]*ghFluid[i] + pRef;
thermoFluid[i].correct();
initialMassf[i] = fvc::domainIntegrate(rhof[i]).value();
initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
}

20
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
View File

@ -1,9 +1,17 @@
solveEnthalpyEquation
{
tmp<fvScalarMatrix> hEqn
(
rhof[i],
DpDtf[i],
phif[i],
turb[i],
thermof[i]
fvm::ddt(rho, h)
+ fvm::div(phi, h)
- fvm::laplacian(turb.alphaEff(), h)
==
DpDt
);
hEqn().relax();
hEqn().solve();
thermo.correct();
Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T())
<< endl;
}

1
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/initContinuityErrs.H Normal file
View File

@ -0,0 +1 @@
List<scalar> cumulativeContErr(fluidRegions.size(), 0.0);

91
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
View File

@ -1,60 +1,75 @@
{
bool closedVolume = false;
bool closedVolume = pd.needReference();
rhof[i] = thermof[i].rho();
rho = thermo.rho();
volScalarField rUA = 1.0/UEqn().A();
Uf[i] = rUA*UEqn().H();
surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
phif[i] =
fvc::interpolate(rhof[i])
U = rUA*UEqn().H();
surfaceScalarField phiU
(
fvc::interpolate(rho)
*(
(fvc::interpolate(Uf[i]) & fluidRegions[i].Sf())
+ fvc::ddtPhiCorr(rUA, rhof[i], Uf[i], phif[i])
(fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
)
- fvc::interpolate(rhof[i]*rUA*ghf[i])
*fvc::snGrad(rhof[i])
*fluidRegions[i].magSf();
);
// Solve pressure difference
# include "pdEqn.H"
phi = phiU - ghf*fvc::snGrad(rho)*rhorUAf*mesh.magSf();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
(
fvm::ddt(psi, pd)
+ fvc::ddt(psi)*pRef
+ fvc::ddt(psi, rho)*gh
+ fvc::div(phi)
- fvm::laplacian(rho*rUA, pd)
);
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pdEqn.solve(mesh.solver(pd.name() + "Final"));
}
else
{
pdEqn.solve(mesh.solver(pd.name()));
}
if (nonOrth == nNonOrthCorr)
{
phi += pdEqn.flux();
}
}
// Correct velocity field
U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
U.correctBoundaryConditions();
// Update pressure field (including bc)
p == pd + rho*gh + pRef;
DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
// Solve continuity
# include "rhoEqn.H"
// Update pressure field (including bc)
thermof[i].p() == pdf[i] + rhof[i]*ghf[i] + pRef;
DpDtf[i] = fvc::DDt
(
surfaceScalarField("phiU", phif[i]/fvc::interpolate(rhof[i])),
thermof[i].p()
);
// Update continuity errors
compressibleContinuityErrors(cumulativeContErr, rhof[i], thermof[i]);
// Correct velocity field
Uf[i] -= rUA*(fvc::grad(pdf[i]) + fvc::grad(rhof[i])*ghf[i]);
Uf[i].correctBoundaryConditions();
# include "compressibleContinuityErrors.H"
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
thermof[i].p() +=
(
dimensionedScalar
(
"massIni",
dimMass,
initialMassf[i]
)
- fvc::domainIntegrate(thermof[i].psi()*thermof[i].p())
)/fvc::domainIntegrate(thermof[i].psi());
pdf[i] == thermof[i].p() - (rhof[i]*ghf[i] + pRef);
rhof[i] = thermof[i].rho();
p += (massIni - fvc::domainIntegrate(psi*p))/fvc::domainIntegrate(psi);
rho = thermo.rho();
}
// Update thermal conductivity
Kf[i] = thermof[i].Cp()*turb[i].alphaEff();
K = thermoFluid[i].Cp()*turb.alphaEff();
// Update pd (including bc)
pd == p - (rho*gh + pRef);
}

14
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pdEqn.H
View File

@ -1,14 +0,0 @@
solvePressureDifferenceEquation
(
corr,
nCorr,
nNonOrthCorr,
closedVolume,
pdf[i],
pRef,
rhof[i],
thermof[i].psi(),
rUA,
ghf[i],
phif[i]
);

2
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/rhoEqn.H
View File

@ -1 +1 @@
solveContinuityEquation(rhof[i], phif[i]);
solve(fvm::ddt(rho) + fvc::div(phi));

15
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setInitialDeltaT.H
View File

@ -1,15 +0,0 @@
if (adjustTimeStep)
{
if (CoNum > SMALL)
{
runTime.setDeltaT
(
min
(
maxCo*runTime.deltaT().value()/CoNum,
maxDeltaT
)
);
}
}

18
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H Normal file
View File

@ -0,0 +1,18 @@
const fvMesh& mesh = fluidRegions[i];
basicThermo& thermo = thermoFluid[i];
volScalarField& rho = rhoFluid[i];
volScalarField& K = KFluid[i];
volVectorField& U = UFluid[i];
surfaceScalarField phi = phiFluid[i];
compressible::turbulenceModel& turb = turbulence[i];
volScalarField& DpDt = DpDtFluid[i];
const volScalarField& gh = ghFluid[i];
const surfaceScalarField& ghf = ghfFluid[i];
volScalarField& pd = pdFluid[i];
volScalarField& p = thermo.p();
const volScalarField& psi = thermo.psi();
volScalarField& h = thermo.h();
const dimensionedScalar massIni("massIni", dimMass, initialMassFluid[i]);

56
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveEnthalpyEquation.C
View File

@ -1,56 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Solve enthalpy equation
\*---------------------------------------------------------------------------*/
void solveEnthalpyEquation
(
const volScalarField& rho,
const volScalarField& DpDt,
const surfaceScalarField& phi,
const compressible::turbulenceModel& turb,
basicThermo& thermo
)
{
volScalarField& h = thermo.h();
tmp<fvScalarMatrix> hEqn
(
fvm::ddt(rho, h)
+ fvm::div(phi, h)
- fvm::laplacian(turb.alphaEff(), h)
==
DpDt
);
hEqn().relax();
hEqn().solve();
thermo.correct();
Info<< "Min/max T:" << min(thermo.T()) << ' ' << max(thermo.T())
<< endl;
}

2
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
View File

@ -12,4 +12,4 @@
# include "pEqn.H"
}
}
turb[i].correct();
turb.correct();

73
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solvePressureDifferenceEquation.C
View File

@ -1,73 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Solve pressure difference equation
\*---------------------------------------------------------------------------*/
void solvePressureDifferenceEquation
(
const label corr,
const label nCorr,
const label nNonOrthCorr,
bool& closedVolume,
volScalarField& pd,
const dimensionedScalar& pRef,
const volScalarField& rho,
const volScalarField& psi,
const volScalarField& rUA,
const volScalarField& gh,
surfaceScalarField& phi
)
{
closedVolume = pd.needReference();
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pdEqn
(
fvm::ddt(psi, pd)
+ fvc::ddt(psi)*pRef
+ fvc::ddt(psi, rho)*gh
+ fvc::div(phi)
- fvm::laplacian(rho*rUA, pd)
);
//pdEqn.solve();
if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
{
pdEqn.solve(pd.mesh().solver(pd.name() + "Final"));
}
else
{
pdEqn.solve(pd.mesh().solver(pd.name()));
}
if (nonOrth == nNonOrthCorr)
{
phi += pdEqn.flux();
}
}
}

6
applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H Normal file
View File

@ -0,0 +1,6 @@
// fvMesh& mesh = solidRegions[i];
volScalarField& rho = rhos[i];
volScalarField& cp = cps[i];
volScalarField& K = Ks[i];
volScalarField& T = Ts[i];

5
applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
View File

@ -3,10 +3,9 @@
{
solve
(
fvm::ddt(rhosCps[i], Ts[i]) - fvm::laplacian(Ks[i], Ts[i])
fvm::ddt(rho*cp, T) - fvm::laplacian(K, T)
);
}
Info<< "Min/max T:" << min(Ts[i]) << ' ' << max(Ts[i])
<< endl;
Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;
}

2
applications/solvers/multiphase/interFoam/alphaEqn.H
View File

@ -6,7 +6,7 @@
phic = min(interface.cAlpha()*phic, max(phic));
surfaceScalarField phir = phic*interface.nHatf();
for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
surfaceScalarField phiAlpha =
fvc::flux

44
applications/solvers/multiphase/twoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C
View File

@ -81,17 +81,20 @@ Foam::phaseModel::phaseModel
{
Info<< "Reading face flux field " << phiName << endl;
phiPtr_ = new surfaceScalarField
phiPtr_.reset
(
IOobject
new surfaceScalarField
(
phiName,
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
IOobject
(
phiName,
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
)
);
}
else
@ -112,18 +115,21 @@ Foam::phaseModel::phaseModel
}
}
phiPtr_ = new surfaceScalarField
phiPtr_.reset
(
IOobject
new surfaceScalarField
(
phiName,
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fvc::interpolate(U_) & mesh.Sf(),
phiTypes
IOobject
(
phiName,
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fvc::interpolate(U_) & mesh.Sf(),
phiTypes
)
);
}
}

6
applications/solvers/multiphase/twoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H
View File

@ -69,7 +69,7 @@ class phaseModel
volVectorField U_;
//- Fluxes
surfaceScalarField* phiPtr_;
autoPtr<surfaceScalarField> phiPtr_;
public:
@ -133,12 +133,12 @@ public:
const surfaceScalarField& phi() const
{
return *phiPtr_;
return phiPtr_();
}
surfaceScalarField& phi()
{
return *phiPtr_;
return phiPtr_();
}
};

2
applications/utilities/surface/surfaceMeshConvertTesting/Make/options
View File

@ -2,4 +2,4 @@ EXE_INC = \
-I$(LIB_SRC)/triSurface/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude
EXE_LIBS = -lsurfMesh
EXE_LIBS = -ltriSurface -lsurfMesh

2
etc/cshrc
View File

@ -158,7 +158,7 @@ case Linux:
setenv WM_CXXFLAGS '-mabi=64 -fPIC'
setenv WM_LDFLAGS '-mabi=64 -G0'
setenv WM_MPLIB MPI
;;
breaksw
default:
echo Unknown processor type `uname -m` for Linux
breaksw

18
src/Pstream/Allwmake
View File

@ -6,17 +6,17 @@ wmake libso dummy
case "$WM_MPLIB" in
GAMMA)
wmake libso gamma
;;
wmake libso gamma
;;
LAM | *MPI* )
export WM_OPTIONS=${WM_OPTIONS}$WM_MPLIB
set +x
echo
echo "Note: ignore spurious warnings about missing mpicxx.h headers"
set -x
wmake libso mpi
;;
WM_OPTIONS=${WM_OPTIONS}$WM_MPLIB
set +x
echo
echo "Note: ignore spurious warnings about missing mpicxx.h headers"
set -x
wmake libso mpi
;;
esac
# ----------------------------------------------------------------- end-of-file

3
src/meshTools/Make/files
View File

@ -55,14 +55,15 @@ indexedOctree/treeDataTriSurface.C
searchableSurface = searchableSurface
$(searchableSurface)/distributedTriSurfaceMesh.C
$(searchableSurface)/searchableBox.C
$(searchableSurface)/searchableCylinder.C
$(searchableSurface)/searchablePlane.C
$(searchableSurface)/searchablePlate.C
$(searchableSurface)/searchableSphere.C
$(searchableSurface)/searchableSurface.C
$(searchableSurface)/searchableSurfaces.C
$(searchableSurface)/searchableSurfacesQueries.C
$(searchableSurface)/triSurfaceMesh.C
$(searchableSurface)/searchableSurfaceWithGaps.C
$(searchableSurface)/triSurfaceMesh.C
topoSets = sets/topoSets
$(topoSets)/cellSet.C

468
src/meshTools/searchableSurface/searchableCylinder.C Normal file
View File

@ -0,0 +1,468 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "searchableCylinder.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(searchableCylinder, 0);
addToRunTimeSelectionTable(searchableSurface, searchableCylinder, dict);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::pointIndexHit Foam::searchableCylinder::findNearest
(
const point& sample,
const scalar nearestDistSqr
) const
{
pointIndexHit info(false, sample, -1);
vector v(sample - point1_);
// Decompose sample-point1 into normal and parallel component
scalar parallel = (v & unitDir_);
// Remove the parallel component
v -= parallel*unitDir_;
scalar magV = mag(v);
if (parallel <= 0)
{
// nearest is at point1 end cap. Clip to radius.
if (magV < ROOTVSMALL)
{
info.setPoint(point1_);
}
else
{
//info.setPoint(point1_ + min(magV, radius_)*v/magV);
info.setPoint(point1_ + radius_*(v/magV));
}
}
else if (parallel >= magDir_)
{
// nearest is at point2
if (magV < ROOTVSMALL)
{
info.setPoint(point2_);
}
else
{
info.setPoint(point2_ + min(magV, radius_)*v/magV);
}
}
else
{
if (magV < ROOTVSMALL)
{
info.setPoint(point1_ + parallel*unitDir_);
}
else
{
// Project onto radius
info.setPoint(point1_ + parallel*unitDir_ + radius_*v/magV);
}
}
if (magSqr(sample - info.rawPoint()) < nearestDistSqr)
{
info.setHit();
info.setIndex(0);
}
return info;
}
// From http://www.gamedev.net/community/forums/topic.asp?topic_id=467789 -
// intersection of cylinder with ray
void Foam::searchableCylinder::findLineAll
(
const point& start,
const point& end,
pointIndexHit& near,
pointIndexHit& far
) const
{
near.setMiss();
far.setMiss();
// Line as P = start + t*V
const vector V(end-start);
//Pout<< "point1:" << point1_ << " point2:" << point2_
// << " start:" << start << " end:" << end << endl;
const vector x = (start-point1_) ^ unitDir_;
const vector y = V ^ unitDir_;
const scalar d = sqr(radius_);
// Second order equation of the form a*t^2 + b*t + c
const scalar a = (y&y);
const scalar b = 2*(x&y);
const scalar c = (x&x)-d;
const scalar disc = b*b-4*a*c;
//Pout<< "a:" << a << " b:" << b << " c:" << c << " disc:" << disc
// << endl;
if (disc < 0)
{
return;
}
else if (disc < ROOTVSMALL)
{
// Single solution
if (mag(a) > ROOTVSMALL)
{
scalar t = -b/(2*a);
if (t >= 0 && t <= 1)
{
near.setPoint(start + t*V);
near.setHit();
near.setIndex(0);
}
}
}
else
{
if (mag(a) > ROOTVSMALL)
{
scalar sqrtDisc = sqrt(disc);
scalar t1 = (-b + sqrtDisc)/2*a;
scalar t2 = (-b - sqrtDisc)/2*a;
if (t1 < t2)
{
if (t1 >= 0 && t1 <= 1)
{
near.setPoint(start + t1*V);
near.setHit();
near.setIndex(0);
}
if (t2 >= 0 && t2 <= 1)
{
far.setPoint(start + t2*V);
far.setHit();
far.setIndex(0);
}
}
else
{
if (t2 >= 0 && t2 <= 1)
{
near.setPoint(start + t2*V);
near.setHit();
near.setIndex(0);
}
if (t1 >= 0 && t1 <= 1)
{
far.setPoint(start + t1*V);
far.setHit();
far.setIndex(0);
}
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchableCylinder::searchableCylinder
(
const IOobject& io,
const point& point1,
const point& point2,
const scalar radius
)
:
searchableSurface(io),
point1_(point1),
point2_(point2),
magDir_(mag(point2_-point1_)),
unitDir_((point2_-point1_)/magDir_),
radius_(radius)
{
Pout<< "point1_:" << point1_ << endl;
Pout<< "point2_:" << point2_ << endl;
Pout<< "magDir_:" << magDir_ << endl;
Pout<< "unitDir_:" << unitDir_ << endl;
Pout<< "radius_:" << radius_ << endl;
}
Foam::searchableCylinder::searchableCylinder
(
const IOobject& io,
const dictionary& dict
)
:
searchableSurface(io),
point1_(dict.lookup("point1")),
point2_(dict.lookup("point2")),
magDir_(mag(point2_-point1_)),
unitDir_((point2_-point1_)/magDir_),
radius_(readScalar(dict.lookup("radius")))
{
Pout<< "point1_:" << point1_ << endl;
Pout<< "point2_:" << point2_ << endl;
Pout<< "magDir_:" << magDir_ << endl;
Pout<< "unitDir_:" << unitDir_ << endl;
Pout<< "radius_:" << radius_ << endl;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::searchableCylinder::~searchableCylinder()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::wordList& Foam::searchableCylinder::regions() const
{
if (regions_.empty())
{
regions_.setSize(1);
regions_[0] = "region0";
}
return regions_;
}
void Foam::searchableCylinder::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
List<pointIndexHit>& info
) const
{
info.setSize(samples.size());
forAll(samples, i)
{
info[i] = findNearest(samples[i], nearestDistSqr[i]);
}
}
void Foam::searchableCylinder::findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
info.setSize(start.size());
pointIndexHit b;
forAll(start, i)
{
// Pick nearest intersection. If none intersected take second one.
findLineAll(start[i], end[i], info[i], b);
if (!info[i].hit() && b.hit())
{
info[i] = b;
}
}
}
void Foam::searchableCylinder::findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
info.setSize(start.size());
pointIndexHit b;
forAll(start, i)
{
// Discard far intersection
findLineAll(start[i], end[i], info[i], b);
if (!info[i].hit() && b.hit())
{
info[i] = b;
}
}
}
void Foam::searchableCylinder::findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >& info
) const
{
info.setSize(start.size());
forAll(start, i)
{
pointIndexHit near, far;
findLineAll(start[i], end[i], near, far);
if (near.hit())
{
if (far.hit())
{
info[i].setSize(2);
info[i][0] = near;
info[i][1] = far;
}
else
{
info[i].setSize(1);
info[i][0] = near;
}
}
else
{
if (far.hit())
{
info[i].setSize(1);
info[i][0] = far;
}
else
{
info[i].clear();
}
}
}
}
void Foam::searchableCylinder::getRegion
(
const List<pointIndexHit>& info,
labelList& region
) const
{
region.setSize(info.size());
region = 0;
}
void Foam::searchableCylinder::getNormal
(
const List<pointIndexHit>& info,
vectorField& normal
) const
{
normal.setSize(info.size());
normal = vector::zero;
forAll(info, i)
{
if (info[i].hit())
{
vector v(info[i].hitPoint() - point1_);
// Decompose sample-point1 into normal and parallel component
scalar parallel = v & unitDir_;
if (parallel < 0)
{
normal[i] = -unitDir_;
}
else if (parallel > magDir_)
{
normal[i] = -unitDir_;
}
else
{
// Remove the parallel component
v -= parallel*unitDir_;
normal[i] = v/mag(v);
}
}
}
}
void Foam::searchableCylinder::getVolumeType
(
const pointField& points,
List<volumeType>& volType
) const
{
volType.setSize(points.size());
volType = INSIDE;
forAll(points, pointI)
{
const point& pt = points[pointI];
vector v(pt - point1_);
// Decompose sample-point1 into normal and parallel component
scalar parallel = v & unitDir_;
if (parallel < 0)
{
// left of point1 endcap
volType[pointI] = OUTSIDE;
}
else if (parallel > magDir_)
{
// right of point2 endcap
volType[pointI] = OUTSIDE;
}
else
{
// Remove the parallel component
v -= parallel*unitDir_;
if (mag(v) > radius_)
{
volType[pointI] = OUTSIDE;
}
else
{
volType[pointI] = INSIDE;
}
}
}
}
// ************************************************************************* //

225
src/meshTools/searchableSurface/searchableCylinder.H Normal file
View File

@ -0,0 +1,225 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::searchableCylinder
Description
Searching on cylinder
SourceFiles
searchableCylinder.C
\*---------------------------------------------------------------------------*/
#ifndef searchableCylinder_H
#define searchableCylinder_H
#include "treeBoundBox.H"
#include "searchableSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
/*---------------------------------------------------------------------------*\
Class searchableCylinder Declaration
\*---------------------------------------------------------------------------*/
class searchableCylinder
:
public searchableSurface
{
private:
// Private Member Data
//- 'left' point
const point point1_;
//- 'right' point
const point point2_;
//- length of vector point2-point1
const scalar magDir_;
//- normalised vector point2-point1
const vector unitDir_;
//- Radius squared
const scalar radius_;
//- Names of regions
mutable wordList regions_;
// Private Member Functions
//- Find nearest point on cylinder.
pointIndexHit findNearest
(
const point& sample,
const scalar nearestDistSqr
) const;
//- Find intersection with cylinder
void findLineAll
(
const point& start,
const point& end,
pointIndexHit& near,
pointIndexHit& far
) const;
//- Disallow default bitwise copy construct
searchableCylinder(const searchableCylinder&);
//- Disallow default bitwise assignment
void operator=(const searchableCylinder&);
public:
//- Runtime type information
TypeName("searchableCylinder");
// Constructors
//- Construct from components
searchableCylinder
(
const IOobject& io,
const point&,
const point&,
const scalar radius
);
//- Construct from dictionary (used by searchableSurface)
searchableCylinder
(
const IOobject& io,
const dictionary& dict
);
// Destructor
virtual ~searchableCylinder();
// Member Functions
virtual const wordList& regions() const;
//- Whether supports volume type below
virtual bool hasVolumeType() const
{
return true;
}
//- Range of local indices that can be returned.
virtual label size() const
{
return 1;
}
// Multiple point queries.
virtual void findNearest
(
const pointField& sample,
const scalarField& nearestDistSqr,
List<pointIndexHit>&
) const;
virtual void findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
virtual void findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
//- Get all intersections in order from start to end.
virtual void findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >&
) const;
//- From a set of points and indices get the region
virtual void getRegion
(
const List<pointIndexHit>&,
labelList& region
) const;
//- From a set of points and indices get the normal
virtual void getNormal
(
const List<pointIndexHit>&,
vectorField& normal
) const;
//- Determine type (inside/outside/mixed) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
virtual void getVolumeType
(
const pointField&,
List<volumeType>&
) const;
// regIOobject implementation
bool writeData(Ostream&) const
{
notImplemented("searchableCylinder::writeData(Ostream&) const");
return false;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

1
src/postProcessing/foamCalcFunctions/Make/files
View File

@ -7,6 +7,7 @@ field/magSqr/magSqr.C
field/magGrad/magGrad.C
field/div/div.C
field/randomise/randomise.C
field/interpolate/interpolate.C
basic/add/add.C

136
src/postProcessing/foamCalcFunctions/field/interpolate/interpolate.C Normal file
View File

@ -0,0 +1,136 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "interpolate.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace calcTypes
{
defineTypeNameAndDebug(interpolate, 0);
addToRunTimeSelectionTable(calcType, interpolate, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::calcTypes::interpolate::interpolate()
:
calcType()
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::calcTypes::interpolate::~interpolate()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::calcTypes::interpolate::init()
{
Foam::argList::validArgs.append("interpolate");
argList::validArgs.append("fieldName1 .. fieldNameN");
}
void Foam::calcTypes::interpolate::preCalc
(
const argList& args,
const Time& runTime,
const fvMesh& mesh
)
{
if (args.additionalArgs().size() < 2)
{
Info<< nl << "must specify one or more fields" << nl;
args.printUsage();
FatalError.exit();
}
}
void Foam::calcTypes::interpolate::calc
(
const argList& args,
const Time& runTime,
const fvMesh& mesh
)
{
const stringList& params = args.additionalArgs();
for (label fieldi=1; fieldi<params.size(); fieldi++)
{
const word fieldName(params[fieldi]);
IOobject fieldHeader
(
fieldName,
runTime.timeName(),
mesh,
IOobject::MUST_READ
);
// Check field exists
if (fieldHeader.headerOk())
{
bool processed = false;
writeInterpolateField<scalar>(fieldHeader, mesh, processed);
writeInterpolateField<vector>(fieldHeader, mesh, processed);
writeInterpolateField<sphericalTensor>
(
fieldHeader,
mesh,
processed
);
writeInterpolateField<symmTensor>(fieldHeader, mesh, processed);
writeInterpolateField<tensor>(fieldHeader, mesh, processed);
if (!processed)
{
FatalError
<< "Unable to process " << fieldName << nl
<< "No call to interpolate for fields of type "
<< fieldHeader.headerClassName() << nl << nl
<< exit(FatalError);
}
}
else
{
Info<< " No " << fieldName << endl;
}
}
}
// ************************************************************************* //

138
src/postProcessing/foamCalcFunctions/field/interpolate/interpolate.H Normal file
View File

@ -0,0 +1,138 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::calcTypes::interpolate
Description
Interpolates volume fields to surface fields for each time.
SourceFiles
interpolate.C
\*---------------------------------------------------------------------------*/
#ifndef interpolate_H
#define interpolate_H
#include "calcType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace calcTypes
{
/*---------------------------------------------------------------------------*\
Class interpolate Declaration
\*---------------------------------------------------------------------------*/
class interpolate
:
public calcType
{
// Private Member Functions
//- Disallow default bitwise copy construct
interpolate(const interpolate&);
//- Disallow default bitwise assignment
void operator=(const interpolate&);
protected:
// Member Functions
// Calculation routines
//- Initialise - typically setting static variables,
// e.g. command line arguments
virtual void init();
//- Pre-time loop calculations
virtual void preCalc
(
const argList& args,
const Time& runTime,
const fvMesh& mesh
);
//- Time loop calculations
virtual void calc
(
const argList& args,
const Time& runTime,
const fvMesh& mesh
);
// I-O
//- Write interpolate fields
template<class Type>
void writeInterpolateField
(
const IOobject& header,
const fvMesh& mesh,
bool& processed
);
public:
//- Runtime type information
TypeName("interpolate");
// Constructors
//- Construct null
interpolate();
// Destructor
virtual ~interpolate();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace calcTypes
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "writeInterpolateField.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

53
applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveMomentumEquation.C → src/postProcessing/foamCalcFunctions/field/interpolate/writeInterpolateField.C
View File

@ -22,36 +22,41 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Description
Solve momentum equation and return matrix for use in pressure equation
\*---------------------------------------------------------------------------*/
tmp<fvVectorMatrix> solveMomentumEquation
template<class Type>
void Foam::calcTypes::interpolate::writeInterpolateField
(
const bool momentumPredictor,
volVectorField& U,
const volScalarField& rho,
const surfaceScalarField& phi,
const volScalarField& pd,
const volScalarField& gh,
const compressible::turbulenceModel& turb
const IOobject& header,
const fvMesh& mesh,
bool& processed
)
{
// Solve the Momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
+ turb.divDevRhoReff(U)
);
typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfaceFieldType;
UEqn().relax();
if (momentumPredictor)
if (header.headerClassName() == fieldType::typeName)
{
solve(UEqn() == -fvc::grad(pd) - fvc::grad(rho)*gh);
}
Info<< " Reading " << header.name() << endl;
fieldType field(header, mesh);
return UEqn;
Info<< " Calculating interpolate" << header.name() << endl;
surfaceFieldType interpolateField
(
IOobject
(
"interpolate" + header.name(),
mesh.time().timeName(),
mesh,
IOobject::NO_READ
),
fvc::interpolate(field)
);
interpolateField.write();
processed = true;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

4
src/postProcessing/functionObjects/field/fieldAverage/fieldAverage/fieldAverage.C
View File

@ -82,7 +82,7 @@ void Foam::fieldAverage::initialise()
// Add mean fields to the field lists
forAll(faItems_, i)
{
const word fieldName = faItems_[i].fieldName();
const word& fieldName = faItems_[i].fieldName();
if (obr_.foundObject<volScalarField>(fieldName))
{
addMeanField<scalar>(i, meanScalarFields_);
@ -117,7 +117,7 @@ void Foam::fieldAverage::initialise()
{
if (faItems_[i].prime2Mean())
{
const word fieldName = faItems_[i].fieldName();
const word& fieldName = faItems_[i].fieldName();
if (!faItems_[i].mean())
{
FatalErrorIn("Foam::fieldAverage::initialise()")

2
src/thermophysicalModels/radiation/radiationModel/radiationModel/newRadiationModel.C
View File

@ -86,7 +86,7 @@ autoPtr<radiationModel> radiationModel::New
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End radiation
} // End namespace radiation
} // End namespace Foam
// ************************************************************************* //

4
src/turbulenceModels/compressible/RAS/RASModel/RASModel.C
View File

@ -88,7 +88,7 @@ RASModel::RASModel
dimensioned<scalar>::lookupOrAddToDict
(
"kappa",
subDict("wallFunctionCoeffs"),
wallFunctionDict_,
0.4187
)
),
@ -97,7 +97,7 @@ RASModel::RASModel
dimensioned<scalar>::lookupOrAddToDict
(
"E",
subDict("wallFunctionCoeffs"),
wallFunctionDict_,
9.0
)
),

6
wmake/wmake
View File

@ -163,8 +163,10 @@ then
exit $?
elif [ ! -d $MakeDir ]
then
$make -k -f $WM_DIR/MakefileApps \
FOAM_APPS="`find . -maxdepth 1 \( -type d -a ! -name "." -a ! -name "Optional" -a ! -name "Make" \) -printf "%f "`"
# FOAM_APPS=$(find . -maxdepth 1 \( -type d -a ! -name "." -a ! -name Optional -a ! -name Make \) -printf "%f ")
# avoid 'find' with '-printf' ... not entirely portable
FOAM_APPS=$(for d in *; do [ -d "$d" -a "$d" != Optional -a "$d" != Make ] && echo "$d"; done | xargs)
$make -k -f $WM_DIR/MakefileApps FOAM_APPS="$FOAM_APPS"
exit $?
fi