本案例利用OpenFOAM中的simpleFoam求解器计算建筑物外流场。
注：案例位于\incompressible\simpleFoam\wind
”1 计算模型
建筑物模型如下图所示。
外流计算域在blockMeshDict中进行指定。外流场几何尺寸x轴方向[-20, 330]，y方向[-50,230]，z方向[0,140]。
本算例计算网格采用snappyHexMesh进行生成。其实也可以利用第三方网格工具生成后导入到OpenFOAM中。
2 计算网格
网格生成命令为：
surfaceFeatures
blockMesh
snappyHexMesh -overwrite
网格生成完毕后如下图所示。
2.1 构造面特征
在进行网格生成之前，需要先利用程序surfaceFeatures构造几何特征。该程序需要利用字典文件surfaceFeaturesDict，此字典文件位于system路径下，文件内容为：
FoamFile
version 2.0;
format ascii;
class dictionary;
object surfaceFeaturesDict;
// 这里指定算例几何文件
// 添加面特征提取信息
注：在实际应用过程中，可以使用命令foamGet surfaceFeatures直接添加字典文件，之后修改字典文件中的算例几何信息即可。
includedAngle 150;
subsetFeatures
nonManifoldEdges yes;
openEdges yes;
trimFeatures
minElem 0;
minLen 0;
writeObj yes;
2.2 生成背景网格
背景网格使用blockMesh生成。
blockMesh命令需要配合blockMeshDict字典文件一起运行。本算例的blockMeshDict字典文件内容如下所示：
FoamFile
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
// * * * * * * * * * * * * * * * * * * * * //
// 定义计算域尺寸及网格生成所需的数据
backgroundMesh
xMin -20; // L = 350
xMax 330;
yMin -50; // L = 280
yMax 230;
zMin 0;
zMax 140;
xCells 25;
yCells 20;
zCells 10;
convertToMeters 1;
// 本算例外流场是一个长方体
// 这里指定长方体的六个定点坐标
vertices
);
// 创建一个block，并指定沿三个方向的网格数量
blocks
hex (0 1 2 3 4 5 6 7)
simpleGrading (1 1 1)
);
edges
);
// 下面定义边界
// 创建了4个边界：inlet\outlet\ground\frontAndBack
boundary
inlet
type patch;
faces
(0 3 7 4)
);
outlet
type patch;
faces
(1 5 6 2)
);
ground
type wall;
faces
(0 1 2 3)
);
// 边界类型为symmetry
frontAndBack
type symmetry;
faces
(0 4 5 1)
(3 2 6 7)
(4 7 6 5)
);
);
mergePatchPairs
);
2.3 生成体网格
体网格利用snappyHexMesh生成，需要在system文件夹中定义字典文件snappyHexMeshDict。本算例中该文件定义如下所示。
FoamFile
version 2.0;
format ascii;
class dictionary;
object snappyHexMeshDict;
// * * * * * * * * * * * * * * * * * * * //
// 配置文件，一般不改
castellatedMesh on;
snap on;
addLayers off;
geometry
buildings
type triSurfaceMesh;
refinementBox
type searchableBox;
min ( 0 0 0);
max (250 180 90);
};
castellatedMeshControls
features
);
refinementSurfaces
buildings
level (3 3);
patchInfo { type wall; }
refinementRegions
refinementBox
mode inside;
levels ((1E15 2));
locationInMesh (1 1 1);
snapControls
explicitFeatureSnap true;
implicitFeatureSnap false;
addLayersControls
layers
"CAD.*"
nSurfaceLayers 2;
relativeSizes true;
expansionRatio 1.2;
finalLayerThickness 0.5;
minThickness 1e-3;
meshQualityControls
{}
writeFlags
// scalarLevels
// layerSets
// layerFields
);
mergeTolerance 1e-6;
snappyHexMeshDict文件较为复杂，该文件模板位于$WM_PROJECT_DIR/etc/caseDicts/mesh/generation文件中，应用过程中可以使用命令foamGet snappyHexMeshDict将该文件直接添加到system文件夹中，之后修改文件中与算例几何相关的文件名即可。
3 物理模型
本算例计算的是不可压缩湍流流动，未考虑传热。因此仅需在constant文件夹中定义属性文件transportProperties及momentumTransport即可。
注：在早期OpenFOAM版本中，momentumTransport文件名为turbulenceProperties。
在transportProperties文件中定义流体物性，如下所示：
FoamFile
version 2.0;
format ascii;
class dictionary;
object transportProperties;
// * * * * * * * * * * * * * * * * * * * //
transportModel Newtonian;
// 定义运动粘度为1.5e-5 m2/s
nu [0 2 -1 0 0 0 0] 1.5e-05;
在momentumTransport文件中定义湍流模型。本算例采用kEpsilon湍流模型。文件内容如下：
FoamFile
version 2.0;
format ascii;
class dictionary;
object RASProperties;
// * * * * * * * * * * * * * * * * * * * * * //
simulationType RAS;
RAS
RASModel kEpsilon;
turbulence on;
printCoeffs on;
4 边界条件
边界条件的指定位于0文件夹，其中需要指定p、U、k、espilon、nut等物理量。边界名称在文件constant/polyMesh/boundary中指定，该文件如下所示。
FoamFile
version 2.0;
format ascii;
class polyBoundaryMesh;
location "constant/polyMesh";
object boundary;
// * * * * * * * * * * * * * * * //
inlet
type patch;
nFaces 494;
startFace 558513;
outlet
type patch;
nFaces 200;
startFace 559007;
ground
type wall;
inGroups List 1(wall);
nFaces 6310;
startFace 559207;
frontAndBack
type symmetry;
inGroups List 1(symmetry);
nFaces 1000;
startFace 565517;
buildings
type wall;
inGroups List 1(wall);
nFaces 22723;
startFace 566517;
可以看到，一共定义了5个边界：inlet、outlet、ground、frontAndBack以及buildings。其中边界inlet与outlet类型为patch；frontAndBack的类型为symmetry；ground与buildings边界类型为wall。在边界文件中需要指定这些边界条件。OpenFOAM提供了一些便利条件，可以给相同类型的边界指定统一的边界条件，如wall、symmetry类型的边界。
4.1 p文件
p文件中指定各边界上的压力分布。
FoamFile
version 2.0;
format ascii;
class volScalarField;
object p;
// * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
inlet
type zeroGradient;
// 指定边界outlet类型为总压出口
// 不可压缩流体中仅需要指定p0即可
outlet
type totalPressure;
p0 uniform 0;
// 这里是简化用法
// 等效于"(ground|building)"
wall
type zeroGradient;
#includeEtc "caseDicts/setConstraintTypes"
这个边界条件指定文件中并没有特别复杂的地方，下面两点需要注意。
1、以边界类型指定边界条件
如文件中的：
wall
type zeroGradient;
其意义为指定所有wall类型的边界（本算例为ground与buildings）条件。其等效于：
"(ground|buildings)"
type zeroGradient;
2、包含文件setConstraintTypes
在边界文件的最后一行#includeEtc "caseDicts/setConstraintTypes"是包含了一个文件setContraintType，该文件位于文件路径/opt/openfoam8/etc/caseDicts/setConstraintTypes，在该文件中为一些常用的约束边界指定了类型。在本算例中，该语句等效于：
symmetry
type symmetry;
或写成更直观的形式：
frontAndBack
type symmetry;
因此p文件可以改写为：
FoamFile
version 2.0;
format ascii;
class volScalarField;
object p;
// * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
inlet
type zeroGradient;
outlet
type totalPressure;
p0 uniform 0;
"(ground|buildings)"
type zeroGradient;
frontAndBack
type symmetry;
4.3 U文件
U文件内容如下所示。
FoamFile
version 2.0;
format ascii;
class volVectorField;
object U;
// * * * * * * * * * * * * * * //
Uinlet (10 0 0);
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
// 指定入口速度为X方向10 m/s
inlet
type fixedValue;
value uniform $Uinlet;
// 指定出口速度类型为pressureInletOutletVelocity
/* 此边界条件适用于指定了压力的压力边界。
对于流出，应用零梯度条件；对于流入，速度由边界面的法向分量获得。
*/
outlet
type pressureInletOutletVelocity;
value uniform (0 0 0);
wall
type noSlip;
#includeEtc "caseDicts/setConstraintTypes"
4.3 k文件
k文件内容如下所示。
FoamFile
version 2.0;
format ascii;
class volScalarField;
object k;
// approx k = 1.5*(I*U)^2 ; I = 0.1
kInlet 1.5;
dimensions [0 2 -2 0 0 0 0];
internalField uniform $kInlet;
boundaryField
inlet
type fixedValue;
value uniform $kInlet;
outlet
type inletOutlet;
inletValue uniform $kInlet;
value uniform $kInlet;
wall
type kqRWallFunction;
value uniform $kInlet;
#includeEtc "caseDicts/setConstraintTypes"
4.4 epsilon文件
epsilon文件如下所示。
FoamFile
version 2.0;
format ascii;
class volScalarField;
object epsilon;
// Cmu^0.75 * k^1.5 / L ; L =10
epsilonInlet 0.03;
dimensions [0 2 -3 0 0 0 0];
internalField uniform $epsilonInlet;
boundaryField
inlet
type fixedValue;
value uniform $epsilonInlet;
outlet
type inletOutlet;
inletValue uniform $epsilonInlet;
value uniform $epsilonInlet;
wall
type epsilonWallFunction;
value uniform $epsilonInlet;
#includeEtc "caseDicts/setConstraintTypes"
4.5 nut文件
nut文件内容如下所示。
FoamFile
version 2.0;
format ascii;
class volScalarField;
object nut;
// * * * * * * * * * * * * * * * //
dimensions [0 2 -1 0 0 0 0];
internalField uniform 0;
boundaryField
inlet
type calculated;
value uniform 0;
outlet
type calculated;
value uniform 0;
wall
type nutkWallFunction;
value uniform 0;
#includeEtc "caseDicts/setConstraintTypes"
5 求解控制
5.1 controlDict文件
controlDict文件如下所示。
FoamFile
version 2.0;
format ascii;
class dictionary;
object controlDict;
// * * * * * * * * * * * * * //
application simpleFoam;
startFrom latestTime;
startTime 0;
stopAt endTime;
// 迭代计算400步
endTime 400;
deltaT 1;
writeControl timeStep;
writeInterval 50;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
functions
5.2 fvSchemes文件
fvSchemes文件指定求解算法。本算例的fvSchemes文件如下所示。
FoamFile
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
// 指定采用稳态计算
ddtSchemes
default steadyState;
gradSchemes
default Gauss linear;
limited cellLimited Gauss linear 1;
grad(U) $limited;
grad(k) $limited;
grad(epsilon) $limited;
divSchemes
default none;
div(phi,U) bounded Gauss linearUpwind limited;
turbulence bounded Gauss limitedLinear 1;
div(phi,k) $turbulence;
div(phi,epsilon) $turbulence;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
laplacianSchemes
default Gauss linear corrected;
interpolationSchemes
default linear;
snGradSchemes
default corrected;
wallDist
method meshWave;
5.3 fvSolution文件
fvSolution文件指定线性方程组的求解方法。本算例中fvSolution文件如下所示。
FoamFile
version 2.0;
format ascii;
class dictionary;
object fvSolution;
// * * * * * * * * * * * * * * //
solvers
solver GAMG;
smoother GaussSeidel;
tolerance 1e-6;
relTol 0.1;
"(U|k|omega|epsilon)"
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-6;
relTol 0.1;
SIMPLE
residualControl
p 1e-4;
U 1e-4;
"(k|omega|epsilon)" 1e-4;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
potentialFlow
nNonOrthogonalCorrectors 10;
relaxationFactors
fields
p 0.3;
equations
U 0.7;
"(k|omega|epsilon).*" 0.7;
6 并行计算
前面演示的是采用串行计算，实际上对于稍微复杂的3D模型，建议使用并行计算。若想要采用并行计算，可以先往system文件夹中添加字典文件decomposeParDict，
foamGet decomposeParDict
并修改文件内容为：
FoamFile
version 2.0;
format ascii;
class dictionary;
object decomposeParDict;
// 设置使用40核进行计算
numberOfSubdomains 40;
// 设置并行方法为scotch，该方法无需任何参数
method scotch;
然后在终端中输入命令：
decomposePar
这样就可以将计算区域分解。然后计算可以使用下面的命令：
mpirun -np 40 simpleFoam -parallel