计算力学OpenFOAM讲座-基础篇(陈丽萍)

OpenFOAM的基础篇陈丽萍南京工业大学城建学院暖通研究所[e-mail]第一部分OpenFOAM的简介OpenFOAM（OpenSourceFieldOperationandManipulation）第一部分OpenFOAM简介其官方网站为”库”（library）C++语言编写得到可执行的软件（application）求解器solver工具utility数值计算后处理流体流动、化学反应、换热、结构动力学、电磁场、金融评估等建模、网格、边界条件前处理计算结果显示和处理1.OpenFOAM框架2.OpenFOAM应用分类应用分类(1)直接利用标准求解器(2)自定义求解器(3)自定义离散方法替代商业软件免费喔这类需求较大，商业软件与之不可相比的更高级的应用下面再谈谈分类（2)和（3）的魅力商业软件中的所谓udf，usersubroutine和这是不可相比的。总之，将OpenFOAM作为类库来构建自己的求解器，这是其它软件无法实现的。分类（2）自定义求解器的魅力按照自己的求解流程来编写求解器，关注的重点是求解的流程。需要关心离散和求解的最底层的知识吗？开始初始值用CFD进行室内温度计算辐射传热计算Qri导热计算Qcdi由热平衡求新的对流传热量Qcvim+1收敛条件结束Qcvim赋新值对于研究离散格式、代数求解器用户来说，创建更高效更高精度的离散方法，这需要修改finiteVolume库和OpenFOAM库中对应的代码。分类（3）自定义离散方法的魅力尤其是对流项，尽管OpenFOAM已经提供了基于NVD和TVD的模板和40多种有名的高阶高精度格式，这仍然是不够的，对流项的离散仍然是目前CFD的重点研究方向。高级应用服务未来第二部分OpenFOAM的应用标准求解器solver的应用（以顶盖驱动流为例）算例Case的结构：FoamFile{version2.0;formatascii;classdictionary;objectblockMeshDict;}//********************//convertToMeters0.1;网格文件blockMeshDict文件头vertices((000)(100)(110)(010)(000.1)(100.1)(110.1)(010.1));blocks(hex(01234567)(20201)simpleGrading(111));edges();patches(wallmovingWall((3762))wallfixedWalls((0473)(2651)(1540))emptyfrontAndBack((0321)(4567)));算例所在的目录下，在运行终端窗口输入：blockMesh进行网格生成。网格生成过程中，如有错误，运行终端窗口将显示。查看网格在运行终端输入：paraFoamFoamFile{version2.0;formatascii;classdictionary;objectcontrolDict;}//******************//applicationicoFoam;startFromstartTime;startTime0;stopAtendTime;endTime0.5;deltaT0.005;文件为system/controlDict,典型的controlDict文件如下writeControltimeStep;writeInterval20;purgeWrite0;writeFormatascii;writePrecision6;uncompressed/compressed;writeCompressionuncompressed;timeFormatgeneral;timePrecision6;//时间文件夹精度runTimeModifiableyes;//在求解过程中是否允许修改以上参数文件为system/fvSolution，典型的文件为：（省去文件头）solvers//方程求解器{pPCG//压力采用预条件共轭梯度法（主要用于求解对称矩阵）{preconditionerDIC;tolerance1e-06;//残差relTol0;//迭代容差};UPBiCG//速度采用预条件双共轭梯度法（主要用于反对称矩阵）{preconditionerDILU;//预测器，对角不完全LUtolerance1e-05;//残差relTol0;//迭代容差};}PISO//piso控制参数{nCorrectors2;//修正次数nNonOrthogonalCorrectors0;//非正交修正次数pRefCell0;//压力参考cell的indexpRefValue0;//压力参考值}离散格式的选择：文件：system/fvSchemes，典型文件：（省去文件头）ddtSchemes{defaultEuler;//欧拉离散}gradSchemes//梯度离散{defaultGausslinear;//高斯方法，线性插值grad(p)Gausslinear;//压力的梯度离散}divSchemes//散度离散{defaultnone;//散度的离散（必须指定没有默认值）div(phi,U)Gausslinear;//对流项离散，高斯理论，用线性插值}laplacianSchemes//拉普拉斯项离散{defaultnone;//拉普拉斯项离散，必须指定laplacian(nu,U)Gausslinearcorrected;//扩散项离散采用高斯理论，采用线性插值laplacian((1|A(U)),p)Gausslinearcorrected;//压力方程离散采用高斯理论，线性插值}interpolationSchemes//插值格式{defaultlinear;//默认线性插值interpolate(HbyA)linear;//线性插值}snGradSchemes//梯度法向分量{defaultcorrected;//默认带有非正交修正}fluxRequired//是否计算流律{defaultno;//默认不计算p;//压力需要计算，因为需要利用压力流律修正速度}0目录中压力文件：FoamFile{version2.0;formatascii;classvolScalarField;objectp;}//***********************//dimensions[02-20000];internalFielduniform0;boundaryField{movingWall{typezeroGradient;}fixedWalls{typezeroGradient;}frontAndBack{typeempty;}}//***********************************//dimensions[01-10000];internalFielduniform(000);boundaryField{movingWall{typefixedValue;valueuniform(100);}fixedWalls{typefixedValue;valueuniform(000);}frontAndBack{typeempty;}}速度文件如下：transportProperties字典文件：FoamFile{version2.0;formatascii;classdictionary;locationconstant;objecttransportProperties;}//*********************//nunu[02-10000]0.01;在算例所在的目录下，在运行终端输入：icoFoam等值线或等值面第三部分求解器求解器结构：Info“ReadingtransportProperties\n”endl;//声明属性字典类对象。IOdictionarytransportProperties(IOobject(“transportProperties”,runTime.constant(),//位置mesh,IOobject::MUST_READ,IOobject::NO_WRITE));createFields.H文件dimensionedScalarnu(transportProperties.lookup(nu));//创建压力场InfoReadingfieldp\nendl;volScalarFieldp(IOobject(p,runTime.timeName(),mesh,IOobject::MUST_READ,IOobject::AUTO_WRITE),mesh);//提示并创建速度场InfoReadingfieldU\nendl;volVectorFieldU(IOobject(U,runTime.timeName(),mesh,IOobject::MUST_READ,IOobject::AUTO_WRITE),mesh);//创建界面流律#includecreatePhi.H//设定压力参考cell的index和参考值labelpRefCell=0;scalarpRefValue=0.0;//通过查询system/fvSolution下的PISO进行更新压力参考cell更新并设定参考值。setRefCell(p,mesh.solutionDict().subDict(PISO),pRefCell,pRefValue);intmain(intargc,char*argv[]){#includesetRootCase.H#includecreateTime.H#includecreateMesh.H#includecreateFields.H#includeinitContinuityErrs.H//**********************************//icoFoam.CInfo“\nStartingtimeloop\n”endl;//用runTime控制整个的循环for(runTime++;!runTime.end();runTime++){//显示运行时间Info“Time=”runTime.timeName()nlendl;#include“readPISOControls.H”//读入piso控制参数，写在runTime循环中#includeCourantNo.H“//构造速度方程fvVectorMatrixUEqn(fvm::ddt(U)+fvm::div(phi,U)-fvm::laplacian(nu,U));solve(UEqn==-fvc::grad(p));∂ρui∂τ∂ρujui∂xj=∂∂xjμt∂ui∂xj−∂p∂xi//---PISOloop----速度修正步for(intcorr=0;corrnCorr;corr++){//求解系数矩阵对角元素的倒数。volScalarFieldrUA=1.0/UEqn.A();U=rUA*UEqn.H();//更新速度//流通量修正phi=(fvc::interpolate(U)&mesh.Sf())+fvc::ddtPhiCorr(rUA,U,phi);//调整边值，以便满足连续性条件adjustPhi(phi,U,p);//非正交循环for(intnonOrth=0;nonOrth=nNonOrthCorr;nonOrth++){fvScalarMatrixpEqn//求解压力方程(fvm::laplacian(rUA,p)==fvc::div(phi));pEqn.setReference(pRefCell,pRefValue);pEqn.solve();//根据新求解的压力，更新面流通量if(nonOrth==nNonOrthCorr){phi-=pEqn.flux(