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

CF入门资料推荐CFD参考书推荐CFD代码勘误

syphomn

5691人浏览 · 2022-04-21 19:12:15

syphomn · 2022-04-21 19:12:15 发布

前言

最近在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 plt


def 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, dy


def 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**-5
    min_parameter = 10**-16
    Alpha_computing = lambda a1, a2, a3, a4, a6: (a1 - 2 * a2 + a3) / dx**2 + (
        a4 - 2 * a2 + a6) / dy**2
    r[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_temp1
        beta = 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:
            break
        else:
            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) / 3
    rhs_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_value


def boundary_condition_updating(dx, dy, omega, phi):
    omega[0] = (-4 * phi[1] + 0.5 * phi[2]) / dx**2
    omega[-1] = (-4 * phi[-2] + 0.5 * phi[-3]) / dx**2
    omega[:, 0] = (-4 * phi[:, 1] + 0.5 * phi[:, 2]) / dy**2
    omega[:, -1] = (-4 * phi[:, -2] + 0.5 * phi[:, -3]) / dy**2 - 3 / dy
    return omega


def 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, phi


def 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, phi


def 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_time


def 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_time


if __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