【来源：http://nbviewer.ipython.org/github/barbagroup/CFDPython/blob/master/lessons/03_CFL_Condition.ipynb】

对于前面的关于1D对流问题，计算过程中，采用的网格数为41，时间步数为25，时间步长为0.025，对于给定的条件，此参数配置能够获得较好的计算结果。那么这些计算参数能否影响最终计算结果？比如说增大或减小网格数量？增大或减小时间步长？下面来测试一下。

我们采用案例1的条件进行测试。修改程序如下：

import numpy as np

import matplotlib.pyplot as plt

 

def linearconv(nx):

dx=2.0/(nx-1)

nt=20

dt=0.025

c=1

 

u=np.ones(nx)

u[.5/dx:1/dx+1]=2

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,lw=2,color='red')

plt.title('nx=$%i$'%nx)

plt.xlabel('$x$')

plt.ylabel('$u(x)$')

这里函数的形式，以方便后面的比较。这里分别比较nx=41、61、81、101、121时候，ux分布情况。

 

看出问题了没？

随着NX增大，即网格的加密，计算结果似乎越来越离谱了？是不是有点儿颠覆了俺们的常规认识了呢？这就是计算稳定性所导致的问题。为了找到导致突变的网格数，我们从80开始增大nx。从图中可以看出，在nx=80的时候，u(x)的分布是可以接受的，我们从80开始慢慢增加nx。

当nx增大到82时，情况发生突变，说明此处是一个坎儿。

由于计算中采用了固定的时间步长，因此随着网格的加密，空间步长与时间步长的比值将会越来越小。这里定义库朗数：

式中u为波速，在本例中u=1；σ为库朗数；σ_max为所采用的离散算法确保稳定性的值。对于本例采用的离散算法，σ_max=1，可以计算得nx最大值为81。

据此，可以修改计算程序，通过库朗数自动调整时间步长，以保证计算过程稳定。

修改后的计算程序为：

#-*- coding=utf-8 -*-

import numpy as np

import matplotlib.pyplot as plt

 

def linearconv(nx):

dx=2.0/(nx-1)

nt=20

c=1

sigma = 0.5 #库朗数

dt=sigma*dx/c #计算时间步长

 

u=np.ones(nx)

u[.5/dx:1/dx+1]=2

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,lw=2,color='red')

plt.title('nx=$%i$'%nx)

plt.xlabel('$x$')

plt.ylabel('$u(x)$')

 

linearconv(80)

计算nx超过80后的ux曲线，如图所示。

 

 

从图中可以看出，通过库朗数定义计算时间步长，可以有效的避免计算不稳定情况。上面各图中因为采用的NX不同，所以导致时间步长有差异，在统一的时间步数情况下，最终计算时长存在差异，有感兴趣的童鞋可以通过调整时间步数以使它们的时间点一致。

总结：对于显式非稳态问题，库朗数是控制其计算稳定性的重要参数。减小时间步长或增大网格都有助于计算稳定。