OpenFOAM｜算例 06 建筑风场

本案例利用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;

}

// 这里指定算例几何文件

surfaces ("buildings.obj");

// 添加面特征提取信息

#includeEtc "caseDicts/surface/surfaceFeaturesDict.cfg"

注：在实际应用过程中，可以使用命令foamGet surfaceFeatures直接添加字典文件，之后修改字典文件中的算例几何信息即可。

”

字典文件中包含的surfaceFeaturesDict.cfg文件位于$WM_PROJECT_DIR/etc/caseDicts/surface文件夹中。文件内容如下所示，此文件一般情况下不需要修改。

includedAngle 150;

subsetFeatures

{

nonManifoldEdges yes;

openEdges yes;

}

trimFeatures

{

minElem 0;

minLen 0;

}

writeObj yes;

运行surfaceFeatures命令后会在constant/triSurface文件夹下创建一个名为buildings.eMesh的文件，该文件中包含了所有释放出的几何面信息，在后续的snappyHexMesh生成网格时需要使用到此文件。

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

(

($:backgroundMesh.xMin $:backgroundMesh.yMin $:backgroundMesh.zMin)

($:backgroundMesh.xMax $:backgroundMesh.yMin $:backgroundMesh.zMin)

($:backgroundMesh.xMax $:backgroundMesh.yMax $:backgroundMesh.zMin)

($:backgroundMesh.xMin $:backgroundMesh.yMax $:backgroundMesh.zMin)

($:backgroundMesh.xMin $:backgroundMesh.yMin $:backgroundMesh.zMax)

($:backgroundMesh.xMax $:backgroundMesh.yMin $:backgroundMesh.zMax)

($:backgroundMesh.xMax $:backgroundMesh.yMax $:backgroundMesh.zMax)

($:backgroundMesh.xMin $:backgroundMesh.yMax $:backgroundMesh.zMax)

);

// 创建一个block，并指定沿三个方向的网格数量

blocks

(

hex (0 1 2 3 4 5 6 7)

(

$:backgroundMesh.xCells

$:backgroundMesh.yCells

$:backgroundMesh.zCells

)

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;

}

// * * * * * * * * * * * * * * * * * * * //

// 配置文件，一般不改

#includeEtc "caseDicts/mesh/generation/snappyHexMeshDict.cfg"

castellatedMesh on;

snap on;

addLayers off;

geometry

{

buildings

{

type triSurfaceMesh;

file "buildings.obj"; //指定几何文件

}

refinementBox

{

type searchableBox;

min ( 0 0 0);

max (250 180 90);

}

};

castellatedMeshControls

{

features

(

{ file "buildings.eMesh"; level 1; } //指定几何面信息文件

);

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;

}

// * * * * * * * * * * * * * * * //

5

(

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

{

p

{

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

计算完毕后利用命令reconstruct收集计算结果。

reconstructPar

然后就可以进行后处理了。

7 计算结果

压力分布

速度分布

离地10m横截面上速度分布