上海交通大学船舶海洋与建筑工程学院谢彬Numerical TESTs for PDEs解答5.2.2
/*---------------------------------------------------------------------------*\
License
Changed from OpenFoam 7 icoFoam to nse521Foam
Application
nse522Foam
Description
Designed for solving Numerical TESTs for PDEs 5.2.2.
\*---------------------------------------------------------------------------*/
#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());
sqr3 = 1.73205080757;
Point(1) = {-sqr3, 0, 0, 1e22};
Point(2) = {sqr3, 0, 0, 1e22};
Point(3) = {0, -3.0, 0, 1e22};
Point(4) = {-sqr3*0.8, -0.2, 0, 1e22};
Point(5) = {sqr3*0.8, -0.2, 0, 1e22};
Point(6) = {0, -2.6, 0, 1e22};
//+
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 1};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 4};
//+
Line Loop(1) = {1, 2, 3};
Line Loop(2) = {4, 5, 6};
Plane Surface(3) = {1,2};
Plane Surface(4) = {2};
Transfinite Line {1} = 251 Using Progression 1;
Transfinite Line {2} = 251 Using Progression 1;
Transfinite Line {3} = 251 Using Progression 1;
Transfinite Line {4} = 101 Using Progression 1;
Transfinite Line {5} = 101 Using Progression 1;
Transfinite Line {6} = 101 Using Progression 1;
// Transfinite Surface {1};
// Recombine Surface {1};
Extrude {0, 0, 0.1} {
Surface{3, 4}; Layers{1}; Recombine;
}
Physical Surface("frontAndBack") = {3, 4, 38, 55};
Physical Surface("movingWall") = {17};
Physical Surface("fixedWalls") = {21, 25};
Physical Volume("box") = {1, 2};
/*--------------------------------*- 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 30;
deltaT 5e-4; // grade(20*20) 0.0025; // 20*20 0.005;
writeControl adjustableRunTime;
writeInterval 2; // #calc "$endTime/20";
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable yes;
adjustTimeStep yes; // 调整时间步长
maxCo 0.2; // 最大 Co 数
maxDeltaT 5e-4; // 最大时间步长
// ************************************************************************* //
/*--------------------------------*- 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 "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.01;
// ************************************************************************* //
sqr3 = 1.73205080757;
Point(1) = {-sqr3, 0, 0, 1e22};
Point(2) = {sqr3, 0, 0, 1e22};
Point(3) = {0, -3.0, 0, 1e22};
Point(4) = {-sqr3*0.8, -0.2, 0, 1e22};
Point(5) = {sqr3*0.8, -0.2, 0, 1e22};
Point(6) = {0, -2.6, 0, 1e22};
//+
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 1};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 4};
//+
Line Loop(1) = {1, 2, 3};
Line Loop(2) = {4, 5, 6};
Plane Surface(3) = {1,2};
Plane Surface(4) = {2};
Transfinite Line {1} = 251 Using Progression 1;
Transfinite Line {2} = 251 Using Progression 1;
Transfinite Line {3} = 251 Using Progression 1;
Transfinite Line {4} = 101 Using Progression 1;
Transfinite Line {5} = 101 Using Progression 1;
Transfinite Line {6} = 101 Using Progression 1;
// Transfinite Surface {1};
// Recombine Surface {1};
Extrude {0, 0, 0.1} {
Surface{3, 4}; Layers{1}; Recombine;
}
Physical Surface("frontAndBack") = {3, 4, 38, 55};
Physical Surface("movingWall") = {17};
Physical Surface("fixedWalls") = {21, 25};
Physical Volume("box") = {1, 2};
/*--------------------------------*- 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 30;
deltaT 5e-4; // grade(20*20) 0.0025; // 20*20 0.005;
writeControl adjustableRunTime;
writeInterval 2; // #calc "$endTime/20";
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable yes;
adjustTimeStep yes; // 调整时间步长
maxCo 0.2; // 最大 Co 数
maxDeltaT 5e-4; // 最大时间步长
// ************************************************************************* //
clearvars; clc; % nse522Re1hplot.m % edit
yToPlot = '[0.0000],[-0.0571],[-0.0714],[-0.0857],[-0.1000],[-0.4429],[-0.8714],[-1.1857],[-1.3143],[-1.5000],[-1.5857],[-1.6429],[-1.6714],[-1.6857],[-1.9429],[-1.9571],[-2.1143],[-2.1857],[-2.2000],[-2.3000],[-3.0000]';
yToPlot = str2num(yToPlot); yToPlot = yToPlot';
dataPath = 'C:\myPC\myFile\SJTU_senior_1st_total\Sundry\grad\hw\hw1\data\5.2.2\';
dataName = '1hUY.txt';
dataName = [dataPath, dataName];
dataMat = importdata(dataName);
uToPlot = dataMat;
figure(1);
xllim = -1; xrlim = -xllim; yllim = xllim; yrlim = xrlim;
set(gca, 'xlim', [xllim, xrlim]); hold on;
set(gca, 'ylim', [yllim, yrlim]); hold on;
yToPlot = yToPlot .*(2/3) + ones(size(yToPlot, 1), 1);
plot(uToPlot, yToPlot, 'Black', 'LineWidth', 1); hold on;
dataName = '1hUYN.xy';
dataName = [dataPath, dataName];
dataMat = importdata(dataName);
yToPlot = dataMat(:,1); uToPlot = dataMat(:,2);
yToPlot = yToPlot .*(2/3) + ones(size(yToPlot, 1), 1);
plot(uToPlot, yToPlot, 'Black--', 'LineWidth', 1); hold on;
title('5.2.2. Re = 100 at Centre', 'FontSize', 15); hold on; % edit
xlabel('Scaled U or X', 'FontSize', 15); hold on;
ylabel('Scaled V or Y', 'FontSize', 15); hold on;
xToPlot = '[1.1547],[1.0854],[1.0623],[1.0277],[0.9815],[0.9122],[0.8545],[0.0000],[-0.4734],[-0.4965],[-0.6813],[-0.7852],[-0.8545],[-0.8891],[-0.9122],[-1.1547]';
xToPlot = str2num(xToPlot); xToPlot = xToPlot';
xToPlot = xToPlot .* (3^0.5/2);
dataName = '1hVXE.txt';
dataName = [dataPath, dataName];
dataMat = importdata(dataName);
vToPlot = dataMat;
plot(xToPlot, vToPlot, 'Black', 'LineWidth', 1); hold on;
dataName = '1hVXN.xy';
dataName = [dataPath, dataName];
dataMat = importdata(dataName);
xToPlot = dataMat(:,1); vToPlot = dataMat(:,3);
xToPlot = xToPlot .* (3^0.5/2);
plot(xToPlot, vToPlot, 'Black--', 'LineWidth', 1); hold on;
targetPath = 'C:\myPC\myFile\SJTU_senior_1st_total\Sundry\grad\hw\hw1\figs\5.2.2\';
targetName = 'nse522Re1h.png'; % edit
targetName = [targetPath, targetName];
legend('REF', 'NUM', 'Location','NW', 'FontSize', 15); hold on;
saveas(gcf, targetName); close all;
clearvars; clc; % nse522Re5hplot.m % edit
yToPlot = '[0.0000],[-0.0571],[-0.0714],[-0.0857],[-0.1000],[-0.4429],[-0.8714],[-1.1857],[-1.3143],[-1.5000],[-1.5857],[-1.6429],[-1.6714],[-1.6857],[-1.9429],[-1.9571],[-2.1143],[-2.1857],[-2.2000],[-2.3000],[-3.0000]';
yToPlot = str2num(yToPlot); yToPlot = yToPlot';
dataPath = 'C:\myPC\myFile\SJTU_senior_1st_total\Sundry\grad\hw\hw1\data\5.2.2\';
dataName = '5hUYE.txt'; % edit
dataName = [dataPath, dataName];
dataMat = importdata(dataName);
uToPlot = dataMat;
figure(1);
xllim = -1; xrlim = -xllim; yllim = xllim; yrlim = xrlim;
set(gca, 'xlim', [xllim, xrlim]); hold on;
set(gca, 'ylim', [yllim, yrlim]); hold on;
yToPlot = yToPlot .*(2/3) + ones(size(yToPlot, 1), 1);
% plot(uToPlot, yToPlot, 'Black', 'LineWidth', 1); hold on;
dataName = '2hUYN.xy'; % edit
dataName = [dataPath, dataName];
dataMat = importdata(dataName);
yToPlot = dataMat(:,1); uToPlot = dataMat(:,2);
yToPlot = yToPlot .*(2/3) + ones(size(yToPlot, 1), 1);
plot(uToPlot, yToPlot, 'Black:', 'LineWidth', 1); hold on;
dataName = '5hUYN.xy'; % edit
dataName = [dataPath, dataName];
dataMat = importdata(dataName);
yToPlot = dataMat(:,1); uToPlot = dataMat(:,2);
yToPlot = yToPlot .*(2/3) + ones(size(yToPlot, 1), 1);
plot(uToPlot, yToPlot, 'Black--', 'LineWidth', 1); hold on;
title('5.2.2. Re = 200 and 500 at C', 'FontSize', 15); hold on; % edit
xlabel('Scaled U or X', 'FontSize', 15); hold on;
ylabel('Scaled V or Y', 'FontSize', 15); hold on;
xToPlot = '[1.1547],[1.0854],[1.0623],[1.0277],[0.9815],[0.9122],[0.8545],[0.0000],[-0.4734],[-0.4965],[-0.6813],[-0.7852],[-0.8545],[-0.8891],[-0.9122],[-1.1547]';
xToPlot = str2num(xToPlot); xToPlot = xToPlot';
xToPlot = xToPlot .* (3^0.5/2);
dataName = '5hVXE.txt';
dataName = [dataPath, dataName];
dataMat = importdata(dataName);
vToPlot = dataMat;
% plot(xToPlot, vToPlot, 'Black', 'LineWidth', 1); hold on;
dataName = '2hVXN.xy';
dataName = [dataPath, dataName];
dataMat = importdata(dataName);
xToPlot = dataMat(:,1); vToPlot = dataMat(:,3);
xToPlot = xToPlot .* (3^0.5/2);
plot(xToPlot, vToPlot, 'Black:', 'LineWidth', 1); hold on;
dataName = '5hVXN.xy';
dataName = [dataPath, dataName];
dataMat = importdata(dataName);
xToPlot = dataMat(:,1); vToPlot = dataMat(:,3);
xToPlot = xToPlot .* (3^0.5/2);
plot(xToPlot, vToPlot, 'Black--', 'LineWidth', 1); hold on;
legend('Re200', 'Re500', 'Location','NW', 'FontSize', 15); hold on;
targetPath = 'C:\myPC\myFile\SJTU_senior_1st_total\Sundry\grad\hw\hw1\figs\5.2.2\';
targetName = 'nse522Re2&5h.png'; % edit
targetName = [targetPath, targetName];
saveas(gcf, targetName); close all;