本案例演示利用SketchUp得到山地地形，利用SCDM生成流体计算区域，利用Fluent Meshing生成计算网格，利用OpenFOAM计算风环境的基本流程。
1 三维地形生成
三维地形的生成方法有很多，比如利用实测等高线数据、卫星地图遥测数据、DEM数据库等点云数据逆向生成三维几何，不过这些数据都不太容易获取。本文选择利用SketchUp Pro 2020进行地形获取，此方法虽然精度不高，但用于大范围环境计算也还勉强凑合。
-
首先启动SketchUp，利用菜单文件→ 地理位置 → 添加位置…打开地区选择对话框
-
找到想要计算的地形区域（图中白色框内的区域），右侧面板选择提供商Digital Globe（这个是免费的，另外一个提供商Nearmap精度更高，不过需要花银子），点击导入按钮将地形图导入到模型中
注：如果打不开卫星地图，显示网络无法连接的话，则可能需要爬墙。这里的卫星地图是基于谷歌地图的）
”
此时在SketchUp中显示的是平面图，需要将其显示为三维几何。
-
点击菜单项文件 → 地理位置 → 显示地形，以三维方式显示几何
几何模型如图所示。
-
选中几何模型，选择菜单文件 → 导出 → 三维模型输出几何文件
导出几何模型zone.obj。
注：这里可以选择导出为obj或stl，建议导出obj格式。
”
2 计算区域创建
利用SCDM导入并处理几何模型，同时生成计算区域。
-
启动SCDM，导入几何文件 zone.obj -
点击按钮分离即将几何模型中的面片分离
-
点击按钮Auto Skin，并选择图形窗口中的几何模型，如下图所示操作，生成三维曲面
-
删除模型树中多余的面片
-
地形面如下图所示
-
创建一个草图面
-
创建一个矩形草图
-
拉伸矩形草图形成三维几何
-
利用曲面分割实体，将下半部的几何删除掉
-
删除曲面，并调整计算区域尺寸
发现几何尺寸单位没有导入进来，这里对几何三个方向各放大1000倍，将尺寸单位转化为米。
完成后的几何模型尺寸如下图所示量级。
创建计算边界：inlet、outlet、ground、side1、side2及top。
利用Fluent Meshing生成计算网格，如下图所示。
-
输出cas文件。注意在输出文件对话框中取消选项Write Binary File，如下图所示
3 求解计算
3.1 准备文件
从OpenFOAM算例库中拷贝一个与想要计算的问题比较类似的案例。
-
创建工程路径，如文件夹 i:/OpenFOAM/case。利用下面的命令从simpleFoam中导入算例文件windAroundBuildings：
cp -r /opt/openfoam8/tutorials/incompressible/simpleFoam/windAroundBuildings/ .
mv windAroundBuildings/ wind
cd wind
-
将前面生成的网格文件 zone.cas拷贝到文件夹wind中 -
删除文件夹中多余的文件，剩下的wind文件夹中的文件结构如下图所示
3.2 转换网格
-
利用命令转换计算网格
fluent3DMeshToFoam zone.cas
哦豁，出错了。
看错误提示是说在zone.cas文件的第1040318行有无法识别的字符[。
用文本编辑器打开zone.cas文件，找到第1040318行，这一行的内容是关于cad几何的，可以直接删除掉。
注：我不甚清楚为何Fluent Meshing导出的cas文件中会存储这么多的几何模型相关信息，事实上该文件的后面一半的数据全部都可以删除，不会影响到网格。
”
删除此行信息后保存cas文件。重新利用命令fluent3DMeshToFoam zone.cas，此时网格顺利转化，如下图所示。
-
打开文件 constant/polyMesh/boundary查看边界信息，文件内如下所示
FoamFile
{
version 2.0;
format ascii;
class polyBoundaryMesh;
location "constant/polyMesh";
object boundary;
}
// * * * * * * * * * * * * * * * * * * * //
6
(
ground
{
type wall;
inGroups List 1(wall);
nFaces 5173;
startFace 570669;
}
inlet
{
type patch;
nFaces 1288;
startFace 575842;
}
outlet
{
type patch;
nFaces 1094;
startFace 577130;
}
side1
{
type symmetry;
inGroups List 1(symmetry);
nFaces 1069;
startFace 578224;
}
side2
{
type symmetry;
inGroups List 1(symmetry);
nFaces 999;
startFace 579293;
}
top
{
type symmetry;
inGroups List 1(symmetry);
nFaces 4987;
startFace 580292;
}
)
可以看到边界类型完整地转化过来了。
3.3 设置边界条件
这里考虑到高度对风速、湍动能及湍流耗散率的影响，计算方式如下所示。
式中，为摩擦速度；为卡门常数，；为湍流粘度系数，常取；为竖直方向高度；为表面粗糙度；为z方向最小坐标值。
其中：
式中，为高度位置的参考速度；为参考高度。
在计算边界设置时，速度入口采用atmBoundaryLayerInletVelocity边界类型；入口湍动能采用atmBoundaryLayerInletK，入口湍流耗散率采用atmBoundaryLayerInletEpsilon。
1、U文件
U文件内容如下所属。
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type atmBoundaryLayerInletVelocity;
Uref 10.0;
Zref 20;
zDir (0 1 0);
flowDir (1 0 0);
z0 uniform 0.1;
zGround uniform -1937.019;
}
outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
value uniform (0 0 0);
}
ground
{
type slip;
}
"(side1|side2|top)"
{
type symmetry;
}
}
2、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
{
type uniformFixedValue;
uniformValue 0;
}
ground
{
type zeroGradient;
}
"(side1|side2|top)"
{
type symmetry;
}
}
3、k文件
k文件内容如下所示。
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object k;
}
// * * * * * * * * * * * * * * * * * * //
kInlet 1.5;
dimensions [0 2 -2 0 0 0 0];
internalField uniform $kInlet;
boundaryField
{
inlet
{
type atmBoundaryLayerInletK;
Uref 10.0;
Zref 20;
zDir (0 1 0);
flowDir (1 0 0);
z0 uniform 0.1;
zGround uniform -1937.019;
}
outlet
{
type inletOutlet;
inletValue uniform $kInlet;
value uniform $kInlet;
}
ground
{
type zeroGradient;
}
"(side1|side2|top)"
{
type symmetry;
}
#includeEtc "caseDicts/setConstraintTypes"
}
4、epsilon文件
epsilon文件如下所示。
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object epsilon;
}
// * * * * * * * * * * * * * * * //
epsilonInlet 0.03;
dimensions [0 2 -3 0 0 0 0];
internalField uniform $epsilonInlet;
boundaryField
{
inlet
{
type atmBoundaryLayerInletEpsilon;
Uref 10.0;
Zref 20;
zDir (0 1 0);
flowDir (1 0 0);
z0 uniform 0.1;
zGround uniform -1937.019;
}
outlet
{
type inletOutlet;
inletValue uniform $epsilonInlet;
value uniform $epsilonInlet;
}
ground
{
type zeroGradient;
}
"(side1|side2|top)"
{
type symmetry;
}
#includeEtc "caseDicts/setConstraintTypes"
}
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;
}
ground
{
type calculated;
value uniform 0;
}
"(side1|side2|top)"
{
type symmetry;
}
#includeEtc "caseDicts/setConstraintTypes"
}
3.4 求解控制
1、指定controlDict
controlDict字典文件内容如下所示。
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object controlDict;
}
// * * * * * * * * * * * * * * * * * * //
application simpleFoam;
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 5000;
deltaT 1;
writeControl timeStep;
writeInterval 50;
purgeWrite 0;
writeFormat ascii;
writePrecision 12;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
functions
{
#includeFunc residuals
}
// ******************************** //
libs ("libatmosphericModels.so");
2、fvSchemes文件
文件内容如下所示。
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,U) bounded Gauss upwind;
div(phi,epsilon) bounded Gauss upwind;
div(phi,k) bounded Gauss upwind;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
default Gauss linear limited corrected 0.33;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default limited corrected 0.33;
}
3、fvSolution文件
文件内容如下所示：
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver GAMG;
tolerance 1e-7;
relTol 0.1;
smoother GaussSeidel;
}
"(U|k|epsilon)"
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-8;
relTol 0.1;
nSweeps 1;
}
}
SIMPLE
{
nNonOrthogonalCorrectors 0;
residualControl
{
p 1e-3;
U 1e-4;
"(k|epsilon)" 1e-4;
}
}
relaxationFactors
{
fields
{
p 0.3;
}
equations
{
U 0.7;
k 0.7;
epsilon 0.7;
}
}
cache
{
grad(U);
}
3.5 并行计算
算例模型较大，采用并行计算可以有效地提高计算效率。
利用命令添加decomposeParDict：
foamGet decompseParDict
修改字典文件：
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object decomposeParDict;
}
// * * * * * * * * * * * * //
numberOfSubdomains 48;
method scotch;
运行命令进行计算：
decompsePar
mpirun -np 48 simpleFoam -parallel
计算完毕后可以采用两种方式：
-
利用命令 reconstruct汇总结果数据后进行后处理 -
利用 paraFoam -builtin直接进入后处理，在paraview中选择Case Type为Decomposed Case
4 计算结果
计算过程中可以使用下面的命令查看残差。
foamMonitor -l postProcessing/residuals/0/residuals.dat
计算过程收敛残差如下图所示。
-
地面上速度分布
-
地面上压力分布
-
切面上速度分布
算例相关文件下载：
本篇文章来源于微信公众号: CFD之道
评论前必须登录！
注册