计算流体力学CFD入门教程介绍

前言

最近在github上看到一个Oklahoma State University，Suraj Pawar and Omer San等人写的
CFD入门教程CFD: 'Julia A Learning Module Structuriing an Introductory Course on Computational Fluid Dynamics'(https://github.com/surajp92/CFD_Julia), 并附有PDF教程（https://www.mdpi.com/2311-5521/4/3/159）。该教程质量比较高，在这里分享给大家。

教程主要内容

| --- | --- |
| 01 | 1D heat equation: Forward time central space (FTCS) scheme |
| 02 | 1D heat equation: Third-order Runge-Kutta (RK3) scheme |
| 03 | 1D heat equation: Crank-Nicolson (CN) scheme |
| 04 | 1D heat equation: Implicit compact Pade (ICP) scheme |
| 05 | 1D inviscid Burgers equation: WENO-5 with Dirichlet and periodic boundary condition |
| 06 | 1D inviscid Burgers equation: CRWENO-5 with Dirichlet and periodic boundary conditions |
| 07 | 1D inviscid Burgers equation: Flux-splitting approach with WENO-5|
| 08 | 1D inviscid Burgers equation: Riemann solver approach with WENO-5 using Rusanov solver |
| 09 | 1D Euler equations: Roe solver, WENO-5, RK3 for time integration |
| 10 | 1D Euler equations: HLLC solver, WENO-5, RK3 for time integration |
| 11 | 1D Euler equations: Rusanov solver, WENO-5, RK3 for time integration |
| 12 | 2D Poisson equation: Finite difference fast Fourier transform (FFT) based direct solver |
| 13 | 2D Poisson equation: Spectral fast Fourier transform (FFT) based direct solver |
| 14 | 2D Poisson equation: Fast sine transform (FST) based direct solver for Dirichlet boundary |
| 15 | 2D Poisson equation: Gauss-Seidel iterative method |
| 16 | 2D Poisson equation: Conjugate gradient iterative method |
| 17 | 2D Poisson equation: V-cycle multigrid iterative method |
| 18 | 2D incompressible Navier-Stokes equations (cavity flow): Arakawa, FST, RK3 schemes |
| 19 | 2D incompressible Navier-Stokes equations (vortex merging): Arakawa, FFT, RK3 schemes |
| 20 | 2D incompressible Navier-Stokes equations (vortex merging): Hybrid RK3/CN approach |
| 21 | 2D incompressible Navier-Stokes equations (vortex merging): Pseudospectral solver, 3/2 dealiasing, Hybrid RK3/CN approach |
| 22 | 2D incompressible Navier-Stokes equations (vortex merging): Pseudospectral solver, 2/3 dealiasing, Hybrid RK3/CN approach |

说明

上述教程虽说是面向高年级本科生以及研究生的CFD入门教程，但也只给出了求解算法及其Julia实现，缺少一定的推导过程，数值计算基础较差的人比较难以理解。因此我收集了一些参考资料，以便读者能更容易地理解教程内容。

另外，github上的代码几乎是用Julia完成的，我也用Python（少部分使用了Matlab）重写了一遍。Python代码及参考资料都可在https://download.csdn.net/download/liuqihang11/85195758下载，不过需要付费9.9元。无法使用github，对Julia代码不熟悉以及需要一些数值计算理论辅助的可考虑下载我整理的文档，否则直接点击开头的两个链接下载PDF以及Julia代码即可。

另外需要注意，Suraj Pawar等人给的PDF教程有些地方有错，这是我按照PDF所给算法编写Python代码时发现的（事实上他们给的Julia代码没有错误，PDF文档中的应该是编辑错误）。下面我直接给大家指出来：

(1) Page 15/77, 3rd term of Equation(40) ,
beta2=13/12(ui − 2ui+1 + ui+2)2 + 14(3ui − 4ui+1 + ui+2)2
instead of
13/12(ui − 2ui+1 + ui+2)2 + 14(3ui − 4ui+1 + 3ui+2)2

(2) Page 25/77,Equation(59)，
fi+1/2 = 1/2 fLi+1/2 + fRi+1/2 − (ci+1/2)/2 (uRi+1/2 − uLi+1/2)
instead of
fi+1/2 = 1/2 fLi+1/2 + fRi+1/2 − (ci+1/2)/2 (uLi+1/2 − uRi+1/2)

(3) Page 32/77,Equation(79)，
(SL − uL)-rhoR*uR instead of (SL − uL)*rhoR*uR

(4) Page 32/77,Equation(81)，
SL*uL*FL instead of SL*uL-FL
SR*uR*FR instead of SR*uR-FR

(5) Page 32/77,Equation(105)，
(2/delta_x**2+2/delta_y**2)**(-1) instead of 2/delta_x**2+2/delta_y**2

(6) Page 51/77,Equation(120)，
e=(...)/4 instead of e=(...)/2

上述内容看起来或许有点晦涩，不过比照PDF就很清楚了。

文档说明

最后的最后介绍一下https://download.csdn.net/download/liuqihang11/85195758里面的内容吧。

（1）CFD_Julia-master

这个文件是在https://github.com/surajp92/CFD_Julia'上直接下载的，里面是CFD的julia代码，这并非本人创作内容，放在这里只是为了给大家提供方便。

（2）CFD Julia A Learning Module Structuriing an Introductory Course on Computational Fluid Dynamics

这个文件是在Fluids | Free Full-Text | CFD Julia: A Learning Module Structuring an Introductory Course on Computational Fluid Dynamics上下载的，它是CFD教程的核心内容，介绍了一些基本的CFD算法，(1)中的CFD_Julia-master文件不过是此PDF的一个Julia实现，同样，这也并非本人创作内容，放在这里只是为了给大家提供方便；

（3）CFD study

这个文件是我根据上述PDF教程自己编写的Python代码，与Suraj Pawar等人给的Julia代码功能基本一致。写python时，我虽然尽可能使用了数组运算来提高速度，但计算效率仍旧会比Julia低一些，没办法，Python老毛病了。这里也建议大家用自己擅长的语言重新实现PDF中的算法，以便加深自己对CFD计算过程的理解。

这里给出了一个对顶盖方腔驱动流(lid driven cavity)进行CFD模拟的Python代码作为示例：

#!/usr/bin/env python3
# -*- coding: utf-8 -*-import time
import copy
import numpy as np
from matplotlib import cm
from matplotlib import pyplot as pltdef discrete_computing(x_range, y_range, nx, ny):x = np.linspace(0, x_range, nx + 1)y = np.linspace(0, y_range, ny + 1)dx = x[1] - x[0]dy = y[1] - y[0]return dx, dydef Conjugate_Gradient_Algorithm(dx, dy, nx, ny, f):phi = np.zeros((nx + 1, ny + 1))r = np.zeros_like(phi)q = np.zeros_like(phi)phi_p1 = copy.deepcopy(phi)epsilon = 10**-5min_parameter = 10**-16Alpha_computing = lambda a1, a2, a3, a4, a6: (a1 - 2 * a2 + a3) / dx**2 + (a4 - 2 * a2 + a6) / dy**2r[1:-1, 1:-1] = f[1:-1, 1:-1] - Alpha_computing(phi[2:, 1:-1], phi[1:-1, 1:-1], phi[:-2, 1:-1], phi[1:-1, 2:],phi[1:-1, :-2])p = copy.deepcopy(r)while 1:q[1:-1, 1:-1] = Alpha_computing(p[2:, 1:-1], p[1:-1, 1:-1],p[:-2, 1:-1], p[1:-1,2:], p[1:-1, :-2])alpha_temp1 = (r[1:-1, 1:-1]**2).sum()alpha_temp2 = (q[1:-1, 1:-1] * p[1:-1, 1:-1]).sum()alpha = alpha_temp1 / (alpha_temp2 + min_parameter)phi_p1[1:-1, 1:-1] = phi[1:-1, 1:-1] + alpha * p[1:-1, 1:-1]r[1:-1, 1:-1] -= alpha * q[1:-1, 1:-1]beta_temp1 = (r[1:-1, 1:-1]**2).sum()beta_temp2 = alpha_temp1beta = beta_temp1 / (beta_temp2 + min_parameter)p[1:-1, 1:-1] = r[1:-1, 1:-1] + beta * p[1:-1, 1:-1]if (abs(phi_p1 - phi) > epsilon).sum() == 0:breakelse:phi = copy.deepcopy(phi_p1)return phi_p1[1:-1, 1:-1]def rhs_computing(dx, dy, re_number, omega, phi):rhs_value = np.ones_like(omega)jacobi1 = ((omega[2:, 1:-1] - omega[:-2, 1:-1]) *(phi[1:-1, 2:] - phi[1:-1, :-2]) -(omega[1:-1, 2:] - omega[1:-1, :-2]) *(phi[2:, 1:-1] - phi[:-2, 1:-1])) / (4 * dx * dy)jacobi2 = (omega[2:, 1:-1] *(phi[2:, 2:] - phi[2:, :-2]) - omega[:-2, 1:-1] *(phi[:-2, 2:] - phi[:-2, :-2]) - omega[1:-1, 2:] *(phi[2:, 2:] - phi[:-2, 2:]) + omega[1:-1, :-2] *(phi[2:, :-2] - phi[:-2, :-2])) / (4 * dx * dy)jacobi3 = (omega[2:, 2:] *(phi[1:-1, 2:] - phi[2:, 1:-1]) - omega[:-2, :-2] *(phi[:-2, 1:-1] - phi[1:-1, :-2]) - omega[:-2, 2:] *(phi[1:-1, 2:] - phi[:-2, 1:-1]) + omega[2:, :-2] *(phi[2:, 1:-1] - phi[1:-1, :-2])) / (4 * dx * dy)jacobi = (jacobi1 + jacobi2 + jacobi3) / 3rhs_value[1:-1, 1:-1] = -jacobi + (omega[2:, 1:-1] - 2.0 * omega[1:-1, 1:-1] + omega[:-2, 1:-1]) / (re_number * dx * dx) + (omega[1:-1, 2:] - 2.0 * omega[1:-1, 1:-1] +omega[1:-1, :-2]) / (re_number * dy * dy)return rhs_valuedef boundary_condition_updating(dx, dy, omega, phi):omega[0] = (-4 * phi[1] + 0.5 * phi[2]) / dx**2omega[-1] = (-4 * phi[-2] + 0.5 * phi[-3]) / dx**2omega[:, 0] = (-4 * phi[:, 1] + 0.5 * phi[:, 2]) / dy**2omega[:, -1] = (-4 * phi[:, -2] + 0.5 * phi[:, -3]) / dy**2 - 3 / dyreturn omegadef RK3(nx, ny, dx, dy, dt, re_number, omega, phi):omega1 = np.zeros_like(omega)rhs_value = rhs_computing(dx, dy, re_number, omega, phi)omega1[1:-1, 1:-1] = omega[1:-1, 1:-1] + dt * rhs_value[1:-1, 1:-1]omega1 = boundary_condition_updating(dx, dy, omega1, phi)phi[1:-1, 1:-1] = Conjugate_Gradient_Algorithm(dx, dy, nx, ny, -omega1)rhs_value = rhs_computing(dx, dy, re_number, omega1, phi)omega1[1:-1, 1:-1] = 0.75 * omega[1:-1, 1:-1] + 0.25 * omega1[1:-1, 1:-1] + 0.25 * dt * rhs_value[1:-1, 1:-1]omega1 = boundary_condition_updating(dx, dy, omega1, phi)phi[1:-1, 1:-1] = Conjugate_Gradient_Algorithm(dx, dy, nx, ny, -omega1)rhs_value = rhs_computing(dx, dy, re_number, omega1, phi)omega[1:-1, 1:-1] = (1.0 / 3.0) * omega[1:-1, 1:-1] + (2.0 / 3.0) * omega1[1:-1, 1:-1] + (2.0 / 3.0) * dt * rhs_value[1:-1, 1:-1]omega = boundary_condition_updating(dx, dy, omega, phi)phi[1:-1, 1:-1] = Conjugate_Gradient_Algorithm(dx, dy, nx, ny, -omega)return omega, phidef lid_driven_cavity_problem_solver(nx, ny, nt, dx, dy, dt, re_number):omega = np.zeros((nx + 1, ny + 1))phi = np.zeros_like(omega)for _ in range(nt):omega, phi = RK3(nx, ny, dx, dy, dt, re_number, omega, phi)return omega, phidef drow_contour(data1, data2):plt.rc('font', family='Times New Roman')fig = plt.figure("lid_driven_cavity_problem",figsize=(14, 6),constrained_layout=True)level1 = np.arange(-195, 65 + 5, 5)level2 = np.arange(-0.104, 0.002 + 0.002, 0.002)ax1 = plt.subplot(1, 2, 1)cset1 = ax1.contourf(data1.T, level1, cmap=cm.jet, vmin=-180, vmax=60)ax1.set_title("Vorticity field")ax1.set_xlabel("X")ax1.set_ylabel("Y")cbar1 = fig.colorbar(cset1, ax=ax1)cbar1.set_ticks([-180, -150, -120, -90, -60, -30, 0, 30, 60])ax2 = plt.subplot(1, 2, 2)cset2 = ax2.contourf(data2.T, level2, cmap=cm.jet)ax2.set_title("Streamfunction")ax2.set_xlabel("X")ax2.set_ylabel("Y")cbar2 = fig.colorbar(cset2, ax=ax2)cbar2.set_ticks([-0.096, -0.084, -0.072, -0.060, -0.048, -0.036, -0.024, -0.012, 0.000])# plt.savefig('lid_driven_cavity_problem.svg', dpi=500)end_time = time.time()plt.show()return end_timedef result_showing(x_range, y_range, t_range, nx, ny, dt, re_number):dx, dy = discrete_computing(x_range, y_range, nx, ny)nt = int(t_range / dt)omega, phi = lid_driven_cavity_problem_solver(nx, ny, nt, dx, dy, dt,re_number)end_time = drow_contour(omega, phi)return end_timeif __name__ == "__main__":start_time = time.time()end_time = result_showing(1.0, 1.0, 10, 32, 32, 0.01, 100.0)print('Time Consuming: {:.2f}s'.format(end_time - start_time))

计算结果如下：

Time Consuming: 17.91s

上述代码基本没有注释以及变量解释，不过结合PDF内容就很容易理解了，基本不需要额外的注释。

（4）Reference

这个里面是一些数值计算算法的详细介绍。因为PDF教程里面只给了求解算法，没有具体介绍算法推导原理，这里搜集了一些基本数值算法，以供数值计算基础薄弱的读者使用。主要包含以下内容：

1）FFT的详细推导

求解泊松方程时会用到该算法。其实这部分内容可在我发布的https://blog.csdn.net/liuqihang11/article/details/124026272上查看。

2）场论基础

这部分介绍了一些场论基础知识，比如梯度，散度、旋度以及Nabla算子的性质及基本运算等，化简二维不可压缩流体的NS方程时会涉及一些矢量运算。

3）多网格技术

求解泊松方程时会使用到多网格技术。

4）共轭梯度法

求解泊松方程时会用到共轭梯度法。

5）循环三对角矩阵求解

求解周期性边界条件下的一维博格斯方程时会用到循环TDMA法。

6）CFD专业词汇

由于PDF是英文的，对英文不太熟悉的读者对照该文档能更容易理解PDF内容。

7）NumericalRecipesin

详细介绍了大量数值计算算法，包罗万象。

8）euler_1d

介绍了一维欧拉方程的推导过程。