入門CFD，主要參考書目《計算流體力學基礎及其應用》（John D.Anderson?著，吳頌平等?譯）

實現了?第 7.3?節?亞聲速-超聲速等熵噴管流動的CFD解法?的代碼，采用的是?MacCormack 方法。代碼中增加了動態顯示無量綱溫度和無量綱密度的功能，參考的是?python中plot實現即時數據動態顯示方法，最終結果如下：圖中紅點代表喉口所在位置。

不足之處，歡迎指正 ！

???

?? ? ? ?

? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ?

完整代碼如下：

# -*- coding: utf-8 -*- """ Created on Wed Aug 5 19:11:15 2020@author: PANSONG """#### Laval 噴管 擬一維計算 ###############import numpy as np import matplotlib.pyplot as plt#%% 參數，采用無量綱量 #噴管截面形狀 def section(x):return 1+2.2*(x-1.5)**2# 初始條件 def initial_condition(x):rho = 1-0.314*xT = 1-0.2314*xV = (0.1+1.09*x)*np.sqrt(T)return rho,T,V# 理想氣體比熱比 gamma = 1.4## 離散化 max_x = 3 dx = 0.1 num_point = int(max_x/dx + 1) x = np.linspace(0,3,num_point)max_Iteration = 10000#%% 動態繪圖 plt.ion() #開啟interactive mode 成功的關鍵函數 plt.figure() plt.title('Dimensionless Parameters Evolution') plt.xlabel('Iteration')iterations = []#%% CFD 計算 #### 初始化，無量綱量 A = section(x)## 這里給兩個初值，一個 0，一個線性，避免迭代過程中 殘差判斷時 重復判斷 i array0 = np.ones(shape=x.shape)*0.01 rho,T,V = initial_condition(x) rho_history = np.vstack([array0,rho]) T_history = np.vstack([array0,T]) V_history = np.vstack([array0,V])####### MacCormack 方法 ################################### # 向前，向后差分計算，輸入 N 個數據，輸出 N-2 個偏導數 def forward_partial(y,dx):return (y[2:]-y[1:-1])/dxdef backward_partial(y,dx):return (y[1:-1]-y[:-2])/dxdef new_rho_T(original,partial,dt):new = original + np.hstack([0,partial,0])*dt # P_rho_t 首尾補 0 new[-1] = 2*new[-2]-new[-3] # 出口處又外插值確定，入口處保持不變return newdef new_V(V,P_V_t,dt):new_V = V + np.hstack([0,P_V_t,0])*dtnew_V[0] = 2*new_V[1] - new_V[2] #入口的速度可變，外插值確定new_V[-1] = 2*new_V[-2]-new_V[-3]return new_VC = 1 # 柯朗數,小于等于 1;從精度考慮，盡可能接近 1 # 關于柯朗數，通過對線性雙曲型偏微分方程的穩定性分析，要求 C<=1，對應數值解的依賴區域包含精確解依賴區域 # 從解的精度考慮，應使數值解的依賴區域盡可能等于精確解依賴區域，即 C 盡可能接近 1 # 對于非線性偏微分方程，起指導性作用，參見參考書目 P221for i in range(max_Iteration):# 繪制 喉口處 無量綱溫度 和 無量綱密度iterations.append(i)plt.plot(iterations,rho_history[1:,15],'r-',linewidth=1.0)plt.plot(iterations,T_history[1:,15],'b-',linewidth=1.0)plt.legend(['Density','Temperature'])plt.pause(1e-3)########################################## rho_error = max(np.abs((rho_history[i]-rho_history[i-1])/rho_history[i-1]))V_error = max(np.abs((V_history[i]-V_history[i-1])/V_history[i-1]))T_error = max(np.abs((T_history[i]-T_history[i-1])/T_history[i-1]))error = max(rho_error,V_error,T_error) ## 取最大相對誤差if error < 1e-4:breaka = np.sqrt(T) #無量綱當地聲速 dt = min(C*dx/(V+a)) # 計算最大允許時間推進步長 ###### 預估步 ########P_rho_t = -V[1:-1]*forward_partial(rho, dx) - rho[1:-1]*forward_partial(V, dx) - rho[1:-1]*V[1:-1]*forward_partial(np.log(A), dx)rho_pred = new_rho_T(rho,P_rho_t,dt)P_V_t = -V[1:-1]*forward_partial(V, dx) - (1.0/gamma)*(forward_partial(T, dx)+(T[1:-1]/rho[1:-1])*forward_partial(rho, dx))V_pred = new_V(V,P_V_t,dt)P_T_t = -V[1:-1]*forward_partial(T, dx) - (gamma-1.0)*T[1:-1]*(forward_partial(V, dx)+V[1:-1]*forward_partial(np.log(A), dx))T_pred = new_rho_T(T,P_T_t,dt)###### 校正步 ########P_rho_t_pred = ( -V_pred[1:-1]*backward_partial(rho_pred, dx) -rho_pred[1:-1]*backward_partial(V_pred, dx) -rho_pred[1:-1]*V_pred[1:-1]*backward_partial(np.log(A), dx) )P_rho_t_av = 0.5*( P_rho_t_pred + P_rho_t )P_V_t_pred = ( -V_pred[1:-1]*backward_partial(V_pred, dx)-1.0/gamma*backward_partial(T_pred, dx)-1.0/gamma*T_pred[1:-1]/rho_pred[1:-1]*backward_partial(rho_pred, dx) )P_V_t_av = 0.5*( P_V_t_pred + P_V_t )P_T_t_pred = ( -V_pred[1:-1]*backward_partial(T_pred, dx)-(gamma-1.0)*T_pred[1:-1]*backward_partial(V_pred, dx)-(gamma-1.0)*T_pred[1:-1]*V_pred[1:-1]*backward_partial(np.log(A), dx) )P_T_t_av = 0.5*( P_T_t_pred + P_T_t )rho = new_rho_T(rho,P_rho_t_av,dt)rho_history = np.vstack([rho_history,rho])V = new_V(V,P_V_t_av,dt)V_history = np.vstack([V_history,V])T = new_rho_T(T,P_T_t_av,dt)T_history = np.vstack([T_history,T])plt.ioff() #關閉動態繪圖#%% 計算結果可視化 def plot_results(x,y,x0,y0,title='value'):plt.figure()plt.plot(x,y)plt.scatter(x0,y0,color='red')plt.xlabel('Dimensionless distance')plt.title('Steady '+title+' in different position')plot_results(x,rho_history[-1,:],1.5,rho_history[-1,15],title='density') plot_results(x,T_history[-1,:],1.5,T_history[-1,15],title='temperature')p = rho_history[-1,:]*T_history[-1,:] plot_results(x,p,1.5,p[15],title='pression')Ma = V_history[-1,:]/np.sqrt(T_history[-1,:]) plot_results(x,Ma,1.5,Ma[15],title='Mach number')plt.show()

?

總結

以上是生活随笔為你收集整理的亚声速-超声速等熵喷管拟一维流动的CFD解法（附完整代码）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。