上海交通大学船舶海洋与建筑工程学院谢彬Numerical TESTs for PDEs解答5.2.2

/*---------------------------------------------------------------------------*\LicenseChanged from OpenFoam 7 icoFoam to nse521FoamApplicationnse522FoamDescriptionDesigned 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 predictorfvVectorMatrix UEqn(fvm::ddt(U)+ fvm::div(phi, U)- fvm::laplacian(nu, U));volScalarField p_ex(p); // measure previous pif (piso.momentumPredictor()){solve(UEqn == -fvc::grad(p));}// --- PISO loopwhile (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 consistencyconstrainPressure(p, U, phiHbyA, rAU);// Non-orthogonal pressure corrector loopwhile (piso.correctNonOrthogonal()){// Pressure correctorfvScalarMatrix 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 % edityToPlot = '[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;

上海交通大学船舶海洋与建筑工程学院谢彬Numerical TESTs for PDEs解答5.2.2相关推荐

上海交通大学船舶海洋与建筑工程学院谢彬Numerical TESTs for PDEs解答1.2.1
求解器adv121Foam.C /*---------------------------------------------------------------------------*\Chang ...
上交船舶海洋与建筑工程学院考研经验
我是刚考上的18届的上交船舶海洋与建筑工程学院的研究生,今天给大家讲下我们学院的考研情况吧.19年的还没出,先给大家讲下18年上交船舶海洋与建筑工程学院考研报录比吧. 从上表可以看出,我们学院的力学. ...
上海交通大学安泰经济与管理学院荣获WRDS-SSRN创新奖
费城--(美国商业资讯)--沃顿研究数据服务中心(Wharton Research Data Services, WRDS)和SSRN欣然宣布,上海交通大学安泰经济与管理学院荣获2020年亚太区WRD ...
上海交通大学计算机系过敏意,上海交通大学-电子信息与电气工程学院-电子信息与电气工程学院...
2011年7月15-17日,由中国密码学会密码算法专业委员会与密码芯片专业委员会联合主办.上海交通大学计算机系承办.上海市计算机学会信息安全专业委员会协办的2011年密码算法与密码芯片学术会议在上海交 ...
上海交大计算机专业报录比,上海交通大学2020硕士研究生各学院各专业报录比汇总！含数据分析...
原标题:上海交通大学2020硕士研究生各学院各专业报录比汇总!含数据分析今天给各位#2022考研#准备报考上海交通大学的考生们整理了交大2020年,各学院各专业的实际报考人数和实际录取人数情况.且包 ...
【20保研】上海交通大学电子信息与电气工程学院2020级研究生优秀生源选拔通知...
点击文末的阅读原文或者公众号界面左下角的保研夏令营或者公众号回复"夏令营"是计算机/软件等专业的所有保研夏令营信息集合,会一直更新的. 补充说明: (1)在 "四.材料清 ...
上交大计算机推免生源学校,2020上海交通大学电子信息与电气工程学院推免夏令营通知...
为了2020考研的考生可以更方便的了解院校推免推免工作,及时获取院校推免信息,中公考研小编为大家分享"2020上海交通大学电子信息与电气工程学院推免夏令营通知"相关内容, 希望可以 ...
上海交通大学安泰经济与管理学院：大疫当前谈供应链思维
伴随新冠肺炎病例增加,包括医疗设施.物资人员等医疗供需备受关注,实际上,大到国家的宏观调控,小到每家每户的饮食安排,处理好供需矛盾是全社会面临的共性问题.上海交通大学安泰经济与管理学院院长.上海交通大 ...
上海交大计算机学院奖学金,上海交通大学-电子信息与电气工程学院-学生工作办公室...
各位同学: 由于中船奖学金申请人数不够,故申报截止时间推迟到3月23日17:00.公示时间推迟到3月26日-3月28日17:00. 请仍想申请本奖学金的同学于3月23日17:00前将材料交至指定地点, ...

上海交通大学船舶海洋与建筑工程学院谢彬Numerical TESTs for PDEs解答5.2.2

上海交通大学船舶海洋与建筑工程学院谢彬Numerical TESTs for PDEs解答5.2.2相关推荐

最新文章

热门文章