12步解N-S方程之第三步

在前2步中我们已经学会如何使用python编写程序解决简单的一维线性/非线性对流问题。也通过调节参数，观察到了一些有意思的现象。在这些参数中，有一个组合形成的参数非常重要，它就是CFL数。之所以重要，是因为CFL数的取值会影响数值求解。

在第一步中，如果我们将程序稍作改动，使之变为参数输入的形式。代码如下：

# -*- coding: utf-8 -*-
#This is step 3 of python cfd courses.
#Text and code provided under a Creative Commons Attribution license, CC-BY
#. (c) Lorena A. Barba, 2013. Thanks: Gilbert Forsyth for help writing the
#notebooks. NSF for support via CAREER award #1149784.
#
# date author description
# 2014/08/31 wubin 1-d linear convection
#
# pu/pt+c*pu/px=0
# a little equation can teach us so much!!!!
#
import numpy as np
import matplotlib.pyplot as plt
import time,sys #load some system utilities
from IPython.core.display import clear_output #used for inline animation
#
def linearconv(nx):
dx=2.0/(nx-1)
nt=20
dt=0.025
c=1.0
#
u=np.ones(nx)
u[0.5/dx:1+1.0/dx]=2.0
un=np.ones(nx)
u0=np.empty_like(u)
u0=u.copy()
#
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.figure()
xc=np.linspace(0,2,nx)
plt.plot(xc,u,linewidth=2.0,label="converged")
plt.plot(xc,u0,linewidth=2.0,label="initial")
plt.legend(loc='upper right')
plt.savefig("F:\PYTHON\python_cfd\step3.png",dpi=150)
plt.show()

在ipython交互中运行linearconv(nx)函数，nx分别取41、61、81和85。得到如下结果:

nx=41

nx=61

nx=81

nx=85

从图可以看出，随着网格点间距逐渐减小，数值求解的精度在降低，最终发散掉（最后一图）。

这是为什么呢？

我们求解的过程如下：从初始的一个预估值开始，通过上一步的值来估计下一步的结果，循环下去直至收敛。每一个时间间隔dt上“波”传播的距离是dx，而dx与nx相关，保持数值求解稳定的关系式为:

    CFL=0.5
dt=CFL*dx

再次运行程序，得到如下结果：

CFL=0.5,nx=41

CFL=0.5,nx=61

CFL=0.5,nx=81

CFL=0.5,nx=101

CFL=0.5,nx=121

改进后，不管nx怎么增大，都能保证数值求解不发散。