简介
Time Accuracy算例是suntans源代码包中的自带算例(存储在/examples/accuracy目录下),目的是验证SUNTANS模型的时间精确性。根据Armfield&Street(2000)的描述,不可压缩流动数值模拟方法的时间精度很大程度上取决于离散格式中压力项的处理及计算网格的类型。Time Accuracy算例计算了不同时间步长内,SUNTANS模型对流体密度界面不稳定演化过程的模拟,以证明其数值方法在时间上有二阶精度。
计算域与网格
该算例计算区域为垂向的二维立面。区域的水平长度为100m,深度为100m。
水平方向上,该计算域由50个等边三角形网格组成;水平向网格示意图如下:
在垂直方向上,由50个相等的垂直网格组成;垂向网格示意图如下:
定解条件
- 边界条件:所有边壁(左、右、底面)均为壁面零梯度边界(即模型中的类型1边界);
- 初始条件:全域的初始水深为常数100m,即自由水面高程为0m;初始密度场ρ由下式确定,
ρ ρ 0 = − 1 2 Δ ρ ρ 0 t a n h [ 2 t a n h − 1 ( α ) δ ( z + D 2 − ζ ) ] \frac{\rho}{\rho_0 } = -\frac{1}{2} \frac{\Delta \rho}{\rho_0} tanh[\frac{2tanh^{-1}(\alpha)}{\delta} (z+\frac{D}{2}- \zeta)] ρ0ρ=−21ρ0Δρtanh[δ2tanh−1(α)(z+2D−ζ)]
式中,Δρ/ρ0=0.06:上下层流体的密度差;α = 0.99:确定密度跃层范围的参数;δ = 20 m:度跃层厚度;z:垂向空间坐标;D = 100 m:初始水深;ζ = cos(πx/L):表示了密度跃层的垂向密度分布, 其中L表示计算域长度L=100m。
由上述初始化过程得到的计算域密度分布如下图所示:
(上图来自SUNTANS用户手册Page34)
运行算例
打开终端,进入文件夹examples/accuracy下,输入指令
make test
即可成功运行。结果如下图所示。
在运行时,终端中会出现如下内容:
On run 1 of 6 (nsteps = 10, dt = 0.1)
On run 2 of 6 (nsteps = 20, dt = 0.05)
On run 3 of 6 (nsteps = 40, dt = 0.025)
On run 4 of 6 (nsteps = 80, dt = 0.0125)
On run 5 of 6 (nsteps = 160, dt = 0.00625)
On run 6 of 6 (nsteps = 320, dt = 0.003125)
上述信息表示,Time Accuracy算例运行了6次,且在各次重复运行时,时间步长逐步减半(时间步长从第1次的0.1s,到最后的0.003125s)。
注意: 若在运行中出现如下所示的错误:
可在运行make test前,修改Makefile文件。打开Makefile文件,找到第48·49行的内容,如下所示:
$(EXEC): accuracy.o
$(LD) -o $@ $(MATHLIB) accuracy.o
修改为
$(EXEC): accuracy.o
$(LD) -o $@ accuracy.o $(MATHLIB)
之后再输入make test 即可成功运行。
结果查看
成功运行后,终端中显示了如下结果:
Error results (Error(n)/Error(0)):
Your results:
--------------------------------------------------
dt0/dt U W S Q Q0 h
--------------------------------------------------
1 1.0 1.0 1.0 1.0 1.0 1.0
2 3.9 3.9 4.0 4.0 1.1 4.0
4 15.5 15.6 16.2 16.4 2.0 16.4
8 64.9 65.2 68.4 68.9 4.3 69.1
16 323.9 325.0 342.7 344.6 12.5 345.9
Reference results:
1 1.0 1.0 1.0 1.0 1.0 1.0
2 3.9 3.9 4.0 4.1 1.2 4.0
4 15.6 15.5 16.2 16.4 2.1 16.4
8 65.1 64.8 68.1 68.9 4.4 69.1
16 324.7 323.5 340.6 344.8 12.8 345.9
Difference (relative):
1 0.00 0.00 0.00 0.00 0.00 0.00
2 0.00 -0.01 -0.00 0.00 0.02 0.00
4 0.00 -0.01 -0.00 -0.00 0.02 0.00
8 0.00 -0.01 -0.01 -0.00 0.02 -0.00
16 0.00 -0.00 -0.01 0.00 0.02 0.00
运行结果中共含三个数据表。
对于各个数据表,其各行对应不同时间步长case下的运行结果;各表共有7列数据,第一列表示相对时间步长,之后依次为水平速度U、垂向速度W、(相对)密度S、非静压力Q、无外推法下的非静压力Q0 和自由水面高程h。
根据全局的计算结果,由下式得到不同case下的误差,
式中,上标ref表示参考值,其对应结果见上表中Reference results;所得的Error2 见上表中Your results。之后以第一组dt=0.1s的数据为基准,得到各个case的相对误差值,见上表中的Difference (relative)。