不可壓縮庫艾特流的隱式求解（附完整程式碼）

入門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)