SUNTANS模型学习（3）——学习cylinder算例

简介

Cylinder Flow（圆柱绕流）是流体数值模型的经典算例。对于平面对称布置的圆柱绕流工况，在理想状态下，圆柱两侧绕流流线呈对称状。但随着雷诺数的增大，圆柱上下游的流线会逐渐失去对称性，而在圆柱后侧的流场中产生卡门涡街。卡门涡街是流体力学中重要的现象，对于卡门涡街的模拟也常被用于模型检验。
在suntans模型代码包中，内置有一圆柱绕流算例cylinder。模拟的配置文件位于 /examples/cylinder 文件夹中。在本算例中，模型预设了两套网格配置（三角形网格、三角形和矩形混合的网格），以及垂向平均流动与三维流动的参数配置。模拟工况汇总如下：

通过模拟本算例，主要学习参数文件suntans.dat中各主要参数的设定及其物理含义。

网格配置

两套网格配置如下图所示（上图为三角形网格，下图为三角形、矩形混合网格）：

其中的红点为各个网格的Voronoi点。
如果模拟三维流动，则需要设定垂向网格配置。在本算例中，我们设定垂向网格总数为Nkmax=20；垂向网格上疏下密（垂向网格由suntans.dat文件中的rstretch参数控制，这个下一节中细讲）。

参数配置

对于本算例，模型参数由suntans.dat文件给出；文件中的参数共包含以下5个部分：

1. Input file for SUNTANS：主要是网格与物理参数
2. Grid Files：网格文件名称
3. Output Data Files：输出文件的名称
4. Input Data Files：输入文件名称
5. Use output to calculate Strouhal number in mfiles/strouhal.m：用于计算Strouhal数的输出结果

Input file for SUNTANS部分

这是这个参数文件的主体部分，这里我们以三角形网格的三维流动模拟为例，看看各个控制参数。以下就是各参数，输入方式为【参数名】+【空格】+【参数值】；#号后侧的为注释；即

variablename variablevalue # Description of variable

以下是suntans.dat中第一部分的内容，在此附上中文说明：

Nkmax			20	# 最大的垂向网格数字，以此为例，垂向网格最多的地方有20层网格。
stairstep		0	# 该值可取1或0；1表示每层垂向网格间距在水平方向上不变，不过这也有可能导致某处最低网格底面低于实际深度；
                    #0表示在底面附近网格使用partial-stepping，使单元的底面与深度的实际值相一致。
rstretch		-1.05	# 垂向网格拉伸系数，取值范围(-1.1, 1)或(1, 1.1)；当 r<0，底部网格得到了加密，反之，则加密顶上的网格
CorrectVoronoi		-1	# 关于是否启用三角网格修正，包含两种方法，即该值有两个取值：1或-1。
                        # 当取1时，对Voronoi点距离过近的网格进行校正；当取-1时，将钝角三角形的Voronoi点移动到它的形心。
VoronoiRatio		0	# Voronoi点的校正因子。当不启用CorrectVoronoi时，该值取0；
                        # 当CorrectVoronoi=1，该值取[0,1]；当CorrectVoronoi=-1，该值取角度[0,90]；（详见手册）
vertgridcorrect 	0 	#  对最小水深D_min<Δz的情况进行校正；当该值为1，遇此情况时有Δz=D_min，从而使Nkmax = D_max /D_min
IntDepth 		0	# 如果取1，Voronoi点的水深由depth文件给出（空间插值）；否则，水深在initialization.c中设定。
dzsmall			0.1	# 校正前网格最小间距比【这个参数不再使用，不必设置】
scaledepth 		0 	# 对水深的缩放（一般用于测试）；当取值为1，则该功能启动
scaledepthfactor 	0 	# 水深缩放比例因子，在scaledepth=1的时候生效
thetaramptime	        1  	# ramp time（详见手册），给初始化时一个缓慢的物理量变化过程，来抑制瞬态的正压效应；在本例子中，该值取1s。
beta			0	# 盐度扩散系数（本例中取0）
theta			0.55	# 控制时间数值积分的显式/隐式/半隐式的系数（对于水动力控制方程）
# θ = 0: 显式(不稳定)；θ = 0.5: Crank-Nicolson ；θ = 1.0: 隐式；θ = 0.55: 常用设置。
thetaS			0.55	# 控制时间数值积分格式的系数（对于对流-扩散方程）
thetaB			0.55	# 【不再使用的参数】
kappa_s                 0       # 盐度的垂直扩散系数
kappa_sH                0       # 盐度的水平扩散系数（详见手册）
gamma 			1 	# 温度的扩散系数
kappa_T                 0       # 温度的垂直扩散系数
kappa_TH                0       #  温度的水平扩散系数（详见手册）
nu 			5e-4   	# 水的垂直方向粘度 (m^2 s^-1) （对于层流）
nu_H 			5e-4 	# 水的水平方向粘度 (m^2 s^-1) （对于层流）
tau_T			0 	# 风应力
z0T	                0.0 	# 水面边界的粗糙度 （单位：m），用于计算CdT
z0B	                0.0	# 底面边界的粗糙度 （单位：m），用于计算CdB
CdT	                0.0 	# 如果z0T=0，则顶部的阻力系数由这个常数给出；否则，这个参数不生效。如果用free-slip条件，则取-1。
CdB	                -1	# 如果z0B=0，则顶部的阻力系数由这个常数给出；否则，这个参数不生效。如果用free-slip条件，则取-1。
CdW			0.0	# 侧向边界阻力系数；只有当壁面没有延伸到自由表面时，才计算侧壁上的阻力（例如漫滩情况下，不计算侧向阻力）
turbmodel		0	# 紊流模型 (0 - 不使用, 1 - MY25)
AB          3 # Adam-Bashforth法的阶数（2或3，对应2阶或3阶A-B方法）
dt 			2.5e-3	# 时间步长
grav			9.81e10 # 重力加速度（默认为9.81，此算例中采用了一个较大值以控制雷诺数）
Cmax 			1 	# Courant数的最大允许值
nsteps 			1600    # 时间步总数
ntout   		80  	# 输出文件的步频，此处表示每80个dt输出一次
ntprog   		5 	# 报告运算进度的频率 (%)
ntconserve 		1	# 输出 conserved data 的步频（例如质量、体积、势能等），详见phys.c中函数OutputData
nonhydrostatic		1	# 0 - 应用静水压假定, 1 = 非静压模拟
cgsolver		1	# 线性方程组的求解方法 0 = GS, 1 = CG 【这个参数已不再使用】
maxiters		1000	# CG 法（共轭梯度法）迭代的最大次数
qmaxiters		2000	# CG 法迭代的最大次数 （泊松方程）
qprecond		2	# 是否对泊松方程进行预处理（0 - 不处理；1 - Diagonal预处理；2 - Block-diagonal预处理）
epsilon			1e-10 	# CG法中收敛性判别所允许的最小值
qepsilon		1e-5	# CG法中收敛性判别所允许的最小值（压力泊松方程）
resnorm			0	# 在计算CG求解器的公差准则时是否对残差进行归一化（0和1对应了不同的误差形式）
relax			1	# CG求解器的松弛系数
amp 			0.8	# 用于指定boundaries.c中的参数prop->amp （振幅）
omega 			0	# 用于指定boundaries.c中的参数prop->omega （圆频率）
flux 			0	# 用于指定boundaries.c中的参数prop->flux （通量或流量）
timescale		0	#  用于指定boundaries.c中的参数prop->timescale（应用于开边界的时间尺度）
volcheck		0	# 是否检查体积守恒（0-不检查；1-检查）
masscheck		0	# 是否检查质量守恒（0-不检查；1-检查）
interp			1	# 插值方法（0-经典Perot方法；1-二次插值）
nonlinear		2	# 动量方程对流项格式（0-无对流；1-迎风格式；2-中心差分格式）
laxWendroff		0       # Lax Wendroff扩散系数（0-不应用Lax Wendroff法；1-当nonlinear=2时，涡粘系数通过Lax Wendroff法给出）
newcells		0	# 1 if adjust momentum in surface cells as the volume changes, 0 otherwise
wetdry			0       # 是否采用干-湿网格设置 （1-使用, 0-不使用）
Coriolis_f              0	# 科氏力参数 f=2*Omega*sin(phi)
sponge_distance	        0 	# 海绵边界层衰减距离（详见phys.c）
sponge_decay	        0	# 海绵边界层衰减时间系数 
readSalinity		0	# 0-垂直盐度分布在initilization.c中设定；1-从InitSalinityFile中读取垂直盐度分布，文件必须为二进制格式，且必须按照从上到下的顺序包含N_kmax个双精度值。
readTemperature		0	# 与readSalinity类似，当取1时，从InitTemperatureFile中读取垂直温度分布

Grid Files部分

网格文件均是ASCII的二进制文件。suntans.dat文件中需要设定它们的文件名。

pslg   full_circle.dat 	# 平面直线图形文件；详见手册第3章。 (input) 【本例中不采用】
points points.dat	# 网格节点文件 (input)
edges edges.dat		# 网格边配置文件 (input)
cells cells.dat		# 网格连接关系文件 (input)
depth depth.dat		# 水深文件 (若INTERPDEPTH=1，则应用该文件) (input)
nodes nodes.dat   # 节点的拓扑关系 (output) 
celldata celldata.dat	# 每个线程的网格中心数据 (output)
edgedata edgedata.dat	# 每个线程的网格边界数据t (output)
vertspace vertspace.dat	# 垂向网格间距 (output)
topology topology.dat	# 每个线程的网格拓扑关系

各个文件的格式详见手册第7章。

Output Data Files和Input Data Files部分

########################################################################
#  Output Data Files
########################################################################
FreeSurfaceFile   	fs.dat   #自由水面高程
HorizontalVelocityFile 	u.dat   #水平流速
VerticalVelocityFile 	w.dat   #垂向流速
SalinityFile 		s.dat  #盐度场数据
BGSalinityFile 		s0.dat  #
TemperatureFile		T.dat  #温度
PressureFile		q.dat  #动水压力
VerticalGridFile 	g.dat  #垂向网格
ConserveFile		e.dat  #守恒性检查
ProgressFile	        step.dat  #模拟进度
StoreFile		store.dat  #保存数据，用于热启动的文件
StartFile		start.dat  #用于热启动的文件
EddyViscosityFile	nut.dat  #涡粘系数
ScalarDiffusivityFile	kappat.dat  #标量扩散系数
########################################################################
# Input Data Files 【本例中不采用】
########################################################################
InitSalinityFile	sinit.dat  #初始盐度场
InitTemperatureFile	Tinit.dat  #初始温度场

Use output to calculate Strouhal number部分

ProfileVariables	u	         # 输出的变量值
DataLocations		dataxy.dat       # 变量所在位置的x，y坐标
ProfileDataFile		profdata.dat     # 垂向剖面信息
ntoutProfs		10               # 每次输出配置文件数据所间隔的时间步；此处含义为每10个时间步输出配置文件数据
NkmaxProfs		1                # 输出文件的最顶层
numInterpPoints		1                # 是否”仅输出最邻近的点“

注：这个部分在手册中也并未详细描述；在此先略过。

其余参数

在目录/examples/cylinder下输入以下指令以运行三角形网格浅水流动：

make test

指令框中将会显示一些参数未被设定，但模型运行时自动采用了其默认值。详细信息如下：

Warning in MPI_GetValue...retrieved default value of maxFaces=3.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of ntoutStore=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of TVDsalt=4.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of TVDtemp=4.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of TVDturb=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of TVDmomentum=3.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of conserveMomentum=1.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of thetaM=-1.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of mergeArrays=1.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of computeSediments=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of calcage=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of agemethod=1.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of calcaverage=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of latitude=29.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of gmtoffset=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of metmodel=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of varmodel=1.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of nugget=0.10 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of sill=0.90 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of range=100000.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of outputNetcdf=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of netcdfBdy=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of readinitialnc=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of Lsw=2.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of Cda=1100.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of Ce=1400.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of Ch=1400.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of hprecond=1.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of prettyplot=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of linearFS=0.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of fixdzz=1.00 which was not specified in data/suntans.dat.
Warning in MPI_GetValue...retrieved default value of smoothbot=0.00 which was not specified in data/suntans.dat.

这里面的参数包括：

maxFaces=3.00 # 网格的最大边界数，默认应用三角形网格，故取3
ntoutStore=0.00 # 
TVDsalt=4.00 # 求解盐度控制方程所选用的TVD格式
TVDtemp=4.00 # 求解温度控制方程所选用的TVD格式
TVDturb=0.00 # 求解紊流模型控制方程所选用的TVD格式
TVDmomentum=3.00 # 求解动量方程所选用的TVD格式
conserveMomentum=1.00 # 是否进行动量的守恒性检查
thetaM=-1.00 # 求解动量方程所选用的TVD格式
mergeArrays=1.00 # 求解垂向对流项的格式；当大于0.5表示隐式，当等于-1表示在UPredictor()中使用了原来的保守方案
computeSediments=0.00 # 是否求解泥沙模块
calcage=0.00 # 是否模拟水龄
agemethod=1.00 # 水龄模拟的方式 1-使用河流边界，2-使用内部源项
calcaverage=0.00 # 是否计算平均水龄

##以下几个参数用于热量和太阳辐射模块##
latitude=29.00 # 纬度，在计算太阳辐射时候使用
gmtoffset=0.00 # 用于修正太阳辐射的参数
metmodel=0.00 # 关于太阳辐射热量模型 0-无气象数据输入，1-COARE3.0,计算辐射短波和长波
varmodel=1.00 # 太阳辐射模型中所用插值方法 0-反距离比例插值，2-克里金法
nugget=0.10 # 克里金插值中所需给定的参数
sill=0.90 # 克里金插值中所需给定的参数
range=100000.00 # 变差函数参数范围
outputNetcdf=0.00 # 输出文件格式 0-二进制格式，1-netcdf格式
netcdfBdy=0.00 # netcdf文件的输出步频
readinitialnc=0.00 # 是否采用netcdf格式的初始条件
Lsw=2.00 # Light extinction depth (m)
Cda=1100.00 # 阻力系数
Ce=1400.00 # 传热系数
Ch=1400.00 # 传热系数
##########################

hprecond=1.00 # 对于自由水面方程求解食肉采用预处理求解器 0-不使用，1-Jacobi预处理求解器
prettyplot=0.00 # 是否使用二次插值得出的输出值，若为1则使用；否则，请指定
linearFS=0.00 # 对自由水面求解是否采用线性方法
fixdzz=1.00 # 对底网格高度的处理 0-不调整，1-假设底部单元高度必须大于dzsmall*dz[Nkmax-1]，2-假设底部单元高度必须大于dzsmall
smoothbot=0.00 # 在SetFluxHeight和ComputeVelocityVector中是否进行平滑处理

模型运行

考虑到有四个模拟工况，模型代码中预载了下面四个指令来分别运行。

1. 三角形网格二维浅水流动模拟

make test

2. 三角形网格三维流动模拟

make test3d

3. 混合网格二维浅水流动模拟

make test-quad

4. 混合网格三维流动模拟

make test-quad-3d

注意：在每次模拟前，记得先使用make clobber取消编译；之后再使用上述指令进行编译+运行；运行结果生成在 /examples/cylinder/data中。
还需要注意的是使用make clobber会自动删除 /examples/cylinder/data 这个目录的文件夹（删除掉这个data文件夹），所以每次运行的结果记得先存好。

模拟结果

由上一节的内容可知，模型的时间步长dt=2.5e-3s，共运行1600步；ntout=80，故每个文件包含了 1+floor(Nsteps/Ntout) = 21个时间步的流场信息。
利用代码包/mfiles文件夹下的getcdata.m读取二进制文件，可以获得每一时间步及每个网格的流场数据，使用getcdata函数的示例参考SUNTANS模型学习（2）——用MATLAB工具提取Time Accuracy算例的计算数据。
以下展示二维浅水模拟结果中形成的卡门涡街（水平流速云图，T = 1600*dt）：

1. 三角形网格
   
2. 混合网格