前面的文章中,并没有完全解决CFL的问题。有眼力好的童鞋们可能发现了,利用Courant数进行控制之后,虽然计算不会崩溃,但计算精度却下降了。
如没有采用Courant数进行控制之前,网格数为81时,计算结果如图所示。
采用Courant = 0.5进行控制后,相同网格数81时,计算结果如下图所示。
注:为了便于比较,这里将时间步数增加了两倍,使他们的时间保持一致。
可以看到,采用Courant数为0.5进行控制之后,计算出现了一定程度的数值误差。
那么有人可能会问了,既然Courant数的引入会导致数值扩散误差,为什么还要引入这个呢?其实很好理解,引入了虽然会有一定的误差,但好歹能计算出个结果,不用就直接发散崩溃了。工程应用中,能定性分析总好过完全不能分析吧。
那么问题来了,这个Courant数到底取多少才合适呢?为了说明这个问题,将代码修改一下。修改后的代码如下:
import numpy as npimport matplotlib.pyplot as plt# 将nt作为参数是为了控制时间步数,后面观察是方便控制时间步数一致def linearconv(nt,sigma): nx = 401 dx = 2 / ( nx - 1 ) c = 1 dt = sigma * dx u = np.ones(nx) u[int(0.5 / dx):int(1 / dx + 1)] = 2 plt.figure(figsize=(8,2)) plt.plot(np.linspace(0,2,nx),u,'b',lw=3,label='t=0') un = np.ones(nx) for n in range(nt): un = u.copy() for i in range(1,nx): u[i] = un[i]- c* dt / dx * (un[i] - un[i-1]) plt.plot(np.linspace(0,2,nx),u,'r',lw=3,label='t=%.1fs'%(nt*dt)) plt.title('$\sigma$=%.1f'%(sigma)) plt.xlabel('x(m)') plt.ylabel('u(m/s)') plt.legend() plt.show()
这里以Courant数作为参数,来观察库朗数不同条件下的图形情况。在代码中,网格数量及网格尺寸均为固定。
以下是库朗数sigma分别为0.2,0.8,1.0,2.0情况下,0.4s时刻波形情况。
linearconv(401,0.2)linearconv(101,0.8)linearconv(81,1)linearconv(41,2)
波形如下图所示。
从图中可以看出,sigma在0.2~1 s时刻,计算表现良好,随着sigma增大,数值扩散越小,波形保真度越高。然而sigma取2时,计算发散。
那sigma从1.0~2.0之间发生了什么呢?缩小范围看看。
linearconv(121,0.8)linearconv(101,1.0)linearconv(81,1.2)linearconv(71,1.4)
图形如下图所示。
从上面的图形分析,在上面的条件中,sigma取值为1为临界值,大于1计算会发散,小于1计算稳定,等于1计算精度最高。
linearconv(101,1.01)linearconv(101,1.001)
得到的图形如图所示。
证实了前面的猜测。
下面来分析Courant数sigma到底是个什么东西。sigma的表达式为:
式中,c为速度。从量纲上来分析,sigma是一个无量纲数,是一个倍数或者比率,可以简单的认为sigma是波在一个时间步长内穿越的网格数量。通常情况下,我们控制在一个时间步长内波形最多穿越1个网格。即控制:
在前面的程序中,
所以得到:
如上面的程序中,c=1,所以sigma必须小于等于1。下面改变程序中的c为2,如果c=2,那么sigma的临界值应该为0.5,下面来验证一下。
与猜测的情况完全相同,代码在sigma大于0.5时计算发散。
小结:上面分析了Courant数在计算稳定性和控制数值扩散上的作用,不过需要注意的是,Courant数除了与流动参数有关外,还与选用的离散格式有关。关于离散格式,留在后面再讨论。关于Courant的深层讨论,还需要参阅专业的计算流体力学理论,涉及到稳定性分析的内容,这里就不展开讨论了。
注:本系列教程来自国外一个使用Python进行CFD初级理论学习的项目,源项目网址为:http://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/。感兴趣的同学可以去官方主页了解更多信息。
免责声明:本文系网络转载或改编,未找到原创作者,版权归原作者所有。如涉及版权,请联系删