不可压缩库艾特流的隐式求解（附完整代码）

入门CFD，主要参考书目《计算流体力学基础及其应用》（John D.Anderson 著，吴颂平等 译）

实现了 第 9.3 节 数值方法：隐式Crank-Nicolson  的代码，问题比较简单，主要是托马斯算法的实现，原书中没有提到对角占优的问题，该算法具体可参考10分钟理解托马斯算法（tridiagonal matrix algorithm，Thomas algorithm）。代码运行结果如下：

                                           

不足之处，欢迎指正！

完整代码如下：

# -*- coding: utf-8 -*-
"""
Created on Wed Aug 12 16:02:41 2020

@author: PANSONG
"""

###################################
######## 不可压库艾特流动 ##########
###################################

import numpy as np
import matplotlib.pyplot as plt;plt.close('all')


### 流动条件: 采用无量纲形式，只有一个参数雷诺数，最后的定常解固定: u=y
ReD = 5000 

### 计算参数
E_list = [1,5,10] # 表征时间间隔
N = 21 # 网格点数
y = np.linspace(0,1,N)
max_iteration = 500

### 计算：隐式 Crank-Nicolson 方法
### 托马斯算法（追赶法）解三对角方程
def Thomas_algo(b,d,a,c):
    ## 输入 b，d, a 为三条对角线上的系数，c 为方程右边的值;ndarray(),行向量
    ## 增加一个方程组有唯一解的判断 ？    
    b = np.insert(b,obj=0,values=0) ## 在开头插入 0
    a = np.append(a,values=0) ## 在末尾插入 0
    ## 可以使用托马斯算法判断，对角占优
    if (abs(d)-abs(b)-abs(a)).all()>0: 
        
        d_prime = np.zeros(shape=d.shape)
        c_prime = np.zeros(shape=c.shape)
        
        for i in range(len(d_prime)):
            if i == 0:
                d_prime[i] = d[i]
                c_prime[i] = c[i]
            else:        
                d_prime[i] = d[i] - b[i]*a[i-1]/d_prime[i-1]
                c_prime[i] = c[i] - c_prime[i-1]*b[i]/d_prime[i-1]
                
        u = np.zeros(shape=c.shape)
        u[-1] = c_prime[-1]/d_prime[-1]
        
        for i in range(len(u)-1):
            u[-(i+2)] = (c_prime[-(i+2)]-a[-(i+2)]*u[-(i+1)])/d_prime[-(i+2)]
            
        return u
    
    else:
        print('该方程组不适用托马斯算法 ！')
    
dic_list = []
for E in E_list:
    
    u = np.zeros(N)
    u[-1] = 1
    u_history = np.vstack([0.01*np.ones(u.shape),u]) ## 增加一行辅助值，由于计算相对误差
    
    A = -E/2.0*np.ones(shape=u.shape[0]-3)
    B = (E+1.0)*np.ones(shape=u.shape[0]-2)
    
    for i in range(max_iteration):
    
        if  (abs(u_history[-1]-u_history[-2])/(u_history[-2]+1e-10)<1e-04).all():
            
            break
        
        K = (1.0-E)*u[1:-1]+E/2.0*(u[2:]+u[:-2])
        K[-1] = K[-1]-A[-1]
        u[1:-1] = Thomas_algo(A,B,A,K)
        u_history = np.vstack([u_history,u])
    
    u_history = np.delete(u_history,0,axis=0) # 删除添加的辅助行
    
    dic = {'E':E,'u_history':u_history}
    dic_list.append(dic)

######## 结果可视化 #############################################
def find_E(E,dic_list):
    u_history_target = []
    for i in range(len(dic_list)):
        if dic_list[i]['E']==E:
            u_history_target = dic_list[i]['u_history']
            break
    if u_history_target.any():
        return u_history_target
    else:      
        print('Not Found !')

### 瞬态发展过程 ####
E_target = 1
u_history_target = find_E(E_target,dic_list)
time_target = [0,2,12,36,60,240]## 根据需要设置，注意不要超出 u_history 的范围
legend_list = ['0$\Delta$t','2$\Delta$t','12$\Delta$t','36$\Delta$t','60$\Delta$t','240$\Delta$t']

for i in time_target:
    plt.plot(u_history_target[i,:],y)

plt.xlabel('u/${u_e}$')
plt.ylabel('y/D')
plt.title('Evolution of u/${u_e}$')
plt.legend(legend_list)