[问题讨论]使用Python学习CFD初级理论系列一Courant数（4/10）

前面的文章中，并没有完全解决CFL的问题。有眼力好的童鞋们可能发现了，利用Courant数进行控制之后，虽然计算不会崩溃，但计算精度却下降了。 

如没有采用Courant数进行控制之前，网格数为81时，计算结果如图所示。 

采用Courant = 0.5进行控制后，相同网格数81时，计算结果如下图所示。 

注：为了便于比较，这里将时间步数增加了两倍，使他们的时间保持一致。

可以看到，采用Courant数为0.5进行控制之后，计算出现了一定程度的数值误差。

那么有人可能会问了，既然Courant数的引入会导致数值扩散误差，为什么还要引入这个呢？其实很好理解，引入了虽然会有一定的误差，但好歹能计算出个结果，不用就直接发散崩溃了。工程应用中，能定性分析总好过完全不能分析吧。

那么问题来了，这个Courant数到底取多少才合适呢？为了说明这个问题，将代码修改一下。修改后的代码如下：

import numpy as np
import 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/。感兴趣的同学可以去官方主页了解更多信息。

本文转载自微信公众号“CFD之道”，有删减，感谢源作者。