/*---------------------------------------------------------------------------*\
License
Changed from OpenFoam 7 icoFoam to nse521Foam
Application
nse521Foam
Description
Designed for solving Numerical TESTs for PDEs 5.2.1.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "pisoControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
pisoControl piso(mesh);
#include "createFields.H"
#include "initContinuityErrs.H"
scalar my_Re = 0;
scalar my_U = 1;
scalar my_a = 1;
scalar my_nu = mag(nu).value();
my_Re = my_U*my_a/my_nu;
scalar my_tmp = 0;
scalar minLength = 100;
forAll(U,celli)
{
my_tmp = mesh.V()[celli];
if (minLength>=my_tmp)
{
minLength = my_tmp;
}
else {}
}
// using std::pow;
// minLength = pow(minLength,1/3);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "CourantNo.H"
// Momentum predictor
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);
volScalarField p_ex(p); // measure previous p
if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
}
// --- PISO loop
while (piso.correct())
{
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);
adjustPhi(phiHbyA, U, p);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU);
// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (piso.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}
}
#include "continuityErrs.H"
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
Info<< "Re = " << my_Re << endl;
Info<< "min volume = " << minLength << endl;
scalar pL1=0;
scalar up = 0;
scalar low = 0.00001;
forAll(p,celli)
{
up += mag(p_ex[celli]-p[celli]) * mesh.V()[celli];
low += mag(p_ex[celli])*mesh.V()[celli];
}
pL1 = up/low;
Info << "L1 p changing = " << pL1 << endl;
scalar pL2 = 0;
up = 0;
low = 0.00001;
forAll(p,celli)
{
up += mag(p_ex[celli]-p[celli])*mag(p_ex[celli]-p[celli])*mesh.V()[celli];
low += p_ex[celli]*p_ex[celli]*mesh.V()[celli];
}
using std::sqrt;
pL2 = sqrt(up/low);
Info << "L2 p changing = " << pL2 << endl;
scalar pL3 = 0;
scalar up_max = 0;
scalar low_max = 0.00001;
forAll(p,celli)
{
my_tmp = mag(p_ex[celli]-p[celli]);
if (up_max<=my_tmp)
{
up_max = my_tmp;
}
else {}
my_tmp = mag(p_ex[celli]);
if (low_max<=my_tmp)
{
low_max = my_tmp;
}
else {}
}
pL3 = up_max / low_max;
Info << "L3 p changing = " << pL3 << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
dimensionedScalar nu
(
"nu",
dimViscosity,
transportProperties.lookup("nu")
);
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field X\n" << endl;
volVectorField X
(
IOobject
(
"X",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh
);
forAll(X,celli)
{
scalar xx = mesh.C()[celli].x();
scalar yy = mesh.C()[celli].y();
X[celli].x() = xx;
X[celli].y() = yy;
X[celli].z() = 0;
}
#include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
mesh.setFluxRequired(p.name());
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0.5 0 0)
(1 0 0)
(0 0.5 0)
(0.5 0.5 0)
(1 0.5 0)
(0 1 0)
(0.5 1 0)
(1 1 0)
(0 0 0.1)
(0.5 0 0.1)
(1 0 0.1)
(0 0.5 0.1)
(0.5 0.5 0.1)
(1 0.5 0.1)
(0 1 0.1)
(0.5 1 0.1)
(1 1 0.1)
);
blocks
(
hex (0 1 4 3 9 10 13 12) (100 100 1) simpleGrading (2.5 2.5 1)
hex (1 2 5 4 10 11 14 13) (100 100 1) simpleGrading (0.4 2.5 1)
hex (3 4 7 6 12 13 16 15) (100 100 1) simpleGrading (2.5 0.4 1)
hex (4 5 8 7 13 14 17 16) (100 100 1) simpleGrading (0.4 0.4 1)
);
edges
(
);
boundary
(
movingWall
{
type wall;
faces
(
(6 15 16 7)
(7 16 17 8)
);
}
fixedWalls
{
type wall;
faces
(
(3 12 15 6)
(0 9 12 3)
(0 1 10 9)
(1 2 11 10)
(2 5 14 11)
(5 8 17 14)
);
}
frontAndBack
{
type empty;
faces
(
(0 3 4 1)
(1 4 5 2)
(3 6 7 4)
(4 7 8 5)
(9 10 13 12)
(10 11 14 13)
(12 13 16 15)
(13 14 17 16)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0; //lastTime;
stopAt endTime;
endTime 15;
deltaT 2e-3; // grade(20*20) 0.0025; // 20*20 0.005;
writeControl adjustableRunTime;
writeInterval 0.5; // #calc "$endTime/20";
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable yes;
adjustTimeStep yes; // 调整时间步长
maxCo 0.2; // 最大 Co 数
maxDeltaT 2e-3; // 最大时间步长
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class polyBoundaryMesh;
location "constant/polyMesh";
object boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
3
(
movingWall
{
type wall;
inGroups List<word> 1(wall);
nFaces 200;
startFace 79600;
}
fixedWalls
{
type wall;
inGroups List<word> 1(wall);
nFaces 600;
startFace 79800;
}
frontAndBack
{
type empty;
inGroups List<word> 1(empty);
nFaces 80000;
startFace 80400;
}
)
// ************************************************************************* //
Point(1) = {0, 0, 0, 1e22};
Point(2) = {1.0, 0, 0, 1e22};
Point(3) = {1.0, 1.0, 0, 1e22};
Point(4) = {0, 1.0, 0, 1e22};
//+
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
//+
Line Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};
Transfinite Line {3,1} = 201 Using Bump 0.2;
Transfinite Line {2,4} = 201 Using Bump 0.2;
Transfinite Surface {1};
// Recombine Surface {1};
Extrude {0, 0, 0.1} {
Surface{1}; Layers{1}; Recombine;
}
Physical Surface("frontAndBack") = {26, 1};
Physical Surface("movingWall") = {21};
Physical Surface("fixedWalls") = {13,25,17};
Physical Volume("box") = {1};,