【Python4CFD】笔记step1-4

CFDpython - 12 steps to N-S equation

最近在等师兄数据，偶然发现github上有一个很有意思的项目，用python（其实是jupyter）学习CFD，借此机会来摸摸鱼，顺便记录一下做点笔记

00-前言

1. 课程链接CFD Python, a.k.a. the 12 steps to Navier-Stokes
2. 视频链接ME702-Youtube
3. 需要的库环境：

ipywidgets==7.2.1
jupyter==1.0.0
numpy==1.14.3
matplotlib==2.2.2
scipy==1.1.0
sympy==1.1.1

3. 前言部分主要介绍了python+numpy的简单应用，包括以下六个部分：

00-1 Libraries

1. 这一部分就是介绍怎么安装库，个人建议用pycharm+anaconda
2. 注意不要用from xxx import *
3. linspace是一个很有用的函数，用来建立等距的行向量

# np.linspace
myarray = numpy.linspace(0, 5, 10)
# 0-5之间产生10个点
# array([ 0.        ,  0.55555556,  1.11111111,  1.66666667,  2.22222222,
#        2.77777778,  3.33333333,  3.88888889,  4.44444444,  5.        ])

00-2 Variables

1. python无需像c一样声明变量的数据类型
2. 但是要注意，整数/整数=浮点数（int/int=float)

a = 5        #a is an integer 5
b = 'five'   #b is a string of the word 'five'
c = 5.0      #c is a floating point 5  
#%%
type(a) # int
#%%
type(b) # str
#%%
type(c) # float

00-3 Whitespace in Python

1. python中的语句关系由缩进来控制的，而不是像c一样的中括号
2. 这个的话用的例子是loop

00-4 Slicing Arrays

1. 切片，和c一样，从0开始，N-1结束
2. 冒号表示从start到end-1

myvals = numpy.array([1, 2, 3, 4, 5])

myvals[0], myvals[4] # output:1,5

myvals[5] # output: error

myvals[0:3] # output:1 2 3

00-5 Assigning Array Variables

1. 感觉像是python里面比较坑的机制了：引用机制，如果中途修改变量，会导致原变量也进行修改
2. 这个时候需采用copy来解决

a1 = numpy.linspace(1,5,5) # a1=1 2 3 4 5

b = a1

b[3] = 17 # b = a1 = 1 2 3 17 5

a2 = numpy.linspace(1,5,5) # a1=1 2 3 4 5

c = a2.copy()

c[3] = 17 # a2=1 2 3 4 5; c=1 2 3 17 5

00-6 learn more

这一节学到个新东西，原来jupyter里面还可以展示youtube，下次试试能不能展示bilibili

from IPython.display import YouTubeVideo
# a short video about using NumPy arrays, from Enthought
YouTubeVideo('...')

01-Step 1: 1-D Linear Convection 一维线性对流方程

1. 对于一维线性对流方程：
   $$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0$$
2. 其本质是一个双曲型偏微分方程，假设初值$u=u_0$根据特征线可以知道其解析解为:
   $$u=u_0(x\pm ct)$$
3. 该例子是对时间域采用向后一阶差分，对空间域采用向前一阶差分，通过有限差分法对方程进行离散：
   $$\frac{u_i^{n+1}-u_i^n}{\Delta t}+c\frac{u_i^n-u_{i-1}^n}{\Delta x}=0$$
   $$u_i^{n+1}=u_i^n-c\frac{\Delta t}{\Delta x}\left(u_i^n-u_{i-1}^n\right)$$
4. 所有代码如下：

# Remember: comments in python are denoted by the pound sign
import numpy                       #here we load numpy
from matplotlib import pyplot      #here we load matplotlib
import time, sys                   #and load some utilities

#this makes matplotlib plots appear in the notebook (instead of a separate window)
%matplotlib inline   

# 绘制网格：空间节点有41个，时间节点为25个，可以理解为解域是一个41*25的均匀网格
nx = 41  # 加大这个值，会让结果更为精确（即加密网格，但不是所有方程都会精确）
dx = 2 / (nx-1)
nt = 25    #nt is the number of timesteps we want to calculate
dt = .025  #dt is the amount of time each timestep covers (delta t)
c = 1      #assume wavespeed of c = 1

# 设置初始条件：在x属于0.5-1的时候是2，其余时刻是1
u = numpy.ones(nx)      #numpy function ones()
u[int(.5 / dx):int(1 / dx + 1)] = 2  #setting u = 2 between 0.5 and 1 as per our I.C.s
pyplot.plot(numpy.linspace(0, 2, nx), u);

# 设置temp变量，用来表示上一时刻的值，即n时刻，最终计算n+1时刻的值
un = numpy.ones(nx) #initialize a temporary array

# 开始求解方程:原作者说，这种方式比较浪费资源，之后课程有更好的方式
for n in range(nt):  #loop for values of n from 0 to nt, so it will run nt times
    un = u.copy() ##copy the existing values of u into un
    for i in range(1, nx): 
        u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])

pyplot.plot(numpy.linspace(0, 2, nx), u);

5. 网格可以通过如下代码进行观察：

xx=numpy.linspace(0,2,nx)
tt=numpy.linspace(0,dt*nt,nt)

x,y=numpy.meshgrid(xx,tt,indexing='xy')
pyplot.scatter(x,y)

02-Step 2: Nonlinear Convection非线性对流方程

1. 用u替换原来的c：
   $$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=0$$
2. 其有限差分形式为：
   $$u_i^{n+1}=u_i^n-u_i^n\frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n)$$
3. 此处仅需改变求解代码：

for n in range(nt):
    un = u.copy()
    for i in range(1, nx):
        u[i] = un[i] - un[i] * dt / dx * (un[i] - un[i-1]) 

4. 这里结果其实很差，出现了数值振荡，因为markdown的原因就不贴图了

03-Convergence and the CFL Condition

1. 这一节介绍收敛性和库朗数，作者说：“In my experience, CFD students love to make things blow up”还是很搞笑的
2. 这里强烈推荐去看原jupyter notebook，作者用不同的网格数x来显示一维对流方程的结果，发现在nx=85的时候出现了数值振荡
3. 库朗数：声波波速与单位网格移动速度之比，这个比值应该小于1，如果大于1，说明实际波速大于网格移动速度，那么网格无法完全计算声波各处的值，从而出现耗散
   $$\sigma =\frac{u\Delta t}{\Delta x}\le {\sigma }_{\max }$$
4. 可以将dt与dx绑定的形式来确保结果收敛：

import numpy
from matplotlib import pyplot

def linearconv(nx):
    dx = 2 / (nx - 1)
    nt = 20    #nt is the number of timesteps we want to calculate
    c = 1
    sigma = .5
    
    dt = sigma * dx

    u = numpy.ones(nx) 
    u[int(.5/dx):int(1 / dx + 1)] = 2

    un = numpy.ones(nx)

    for n in range(nt):  #iterate through time
        un = u.copy() ##copy the existing values of u into un
        for i in range(1, nx):
            u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])
        
    pyplot.plot(numpy.linspace(0, 2, nx), u)

5. 这里也正好解释了nx变大后，dx变小，从而导致库朗数变大，超过1之后就开始发散

04-Step 3: Diffusion Equation in 1-D一维扩散方程

1. 一维扩散方程形式如下：
   $$\frac{\partial u}{\partial t}=\nu \frac{{\partial }^{2}u}{\partial x^2}$$
2. 二阶中心差分：现在扩散项进行泰勒展开，然后两者相加即可。

$$u_{i+1}=u_i+\Delta x\frac{\partial u}{\partial x}{|}_i+\frac{\Delta x^2}{2}\frac{{\partial }^{2}u}{\partial x^2}{|}_i+\frac{\Delta x^3}{3!}\frac{{\partial }^{3}u}{\partial x^3}{|}_i+O(\Delta x^4)$$

$$u_{i-1}=u_i-\Delta x\frac{\partial u}{\partial x}{|}_i+\frac{\Delta x^2}{2}\frac{{\partial }^{2}u}{\partial x^2}{|}_i-\frac{\Delta x^3}{3!}\frac{{\partial }^{3}u}{\partial x^3}{|}_i+O(\Delta x^4)$$

$$u_{i+1}+u_{i-1}=2u_i+\Delta x^2\frac{{\partial }^{2}u}{\partial x^2}{|}_i+O(\Delta x^4)$$

$$\frac{{\partial }^{2}u}{\partial x^2}=\frac{u_{i+1}-2u_i+u_{i-1}}{\Delta x^2}+O(\Delta x^2)$$
3. 同理，非稳态项采用一阶向前差分，扩散项采用二阶中心差分，最终得到：
$$\frac{u_i^{n+1}-u_i^n}{\Delta t}=\nu \frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{\Delta x^2}$$

$$u_i^{n+1}=u_i^n+\frac{\nu \Delta t}{\Delta x^2}(u_{i+1}^n-2u_i^n+u_{i-1}^n)$$
4. 全部代码如下：

import numpy
from matplotlib import pyplot
%matplotlib inline

nx = 41
dx = 2 / (nx - 1)
nt = 20    #the number of timesteps we want to calculate
nu = 0.3   #the value of viscosity，粘性系数
sigma = .2 #sigma is a parameter, we'll learn more about it later
dt = sigma * dx**2 / nu # 这个时候库朗数变了，多了粘性项，具体看最终形式（v*dt/dx**2）

# 设置初始条件，同之前的step
u = numpy.ones(nx)
u[int(.5 / dx):int(1 / dx + 1)] = 2

# 存放前一时刻-n-时刻的速度值
un = numpy.ones(nx)

# 开始迭代
for n in range(nt):
    un = u.copy()
    for i in range(1, nx - 1):
        u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2 * un[i] + un[i-1])
        # 这里nu * dt / dx**2===sigma
        
pyplot.plot(numpy.linspace(0, 2, nx), u);

05-Step 4: Burgers’ Equation

051-Burgers方程介绍

1. 一维Burger’s Equation就是将以上两个简单对流扩散方程进行的叠加，方程形式如下：
   $$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\nu \frac{{\partial }^{2}u}{\partial x^2}$$
2. 同理，对非稳态项采用一阶向前差分，对对流项采用一阶向后差分，对扩散项采用二阶中心差分，离散后得到如下结果：
   $$\frac{u_i^{n+1}-u_i^n}{\Delta t}+u_i^n\frac{u_i^n-u_{i-1}^n}{\Delta x}=\nu \frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{\Delta x^2}$$
   $$u_i^{n+1}=u_i^n-u_i^n\frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n)+\nu \frac{\Delta t}{\Delta x^2}(u_{i+1}^n-2u_i^n+u_{i-1}^n)$$
3. 初始条件和边界条件这块儿我没懂：边界条件有点复杂，但是解析解是怎么确定的呢？
   • 初始条件长这样：
     $$\begin{array}{rl}u& =-\frac{2\nu }{\varphi }\frac{\partial \varphi }{\partial x}+4\\ \varphi & =\exp \left(\frac{-x^2}{4\nu }\right)+\exp \left(\frac{-(x-2\pi {)}^2}{4\nu }\right)\end{array}$$
   • 解析解长这样：
     $$\begin{array}{rl}u& =-\frac{2\nu }{\varphi }\frac{\partial \varphi }{\partial x}+4\\ \varphi & =\exp \left(\frac{-(x-4t{)}^2}{4\nu }\right)+\exp \left(\frac{-(x-4t-2\pi {)}^2}{4\nu }\right)\end{array}$$
   • 边界条件为：周期性边界条件
     $$u(0)=u(2\pi )$$

052- Sympy

1. sympy类似于Matlab中的symbolic，是一种符号表示的库
2. 以下代码是将初始条件转换为符号函数，再通过sympy进行求导$\frac{\partial \varphi }{\partial x}$:

import numpy 
import sympy

from sympy import init_printing # 用于输出打印
init_printing(use_latex=True)

# 用symbols来表示\phi
x, nu, t=sympy.symbols('x nu t')
phi = (sympy.exp(-(x - 4 * t)**2 / (4 * nu * (t + 1))) +
       sympy.exp(-(x - 4 * t - 2 * sympy.pi)**2 / (4 * nu * (t + 1))))

# 对x进行求导
phiprime = phi.diff(x)
print(phiprime)
# output
# -(-8*t + 2*x)*exp(-(-4*t + x)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)) - (-8*t + 2*x - 4*pi)*exp(-(-4*t + x - 2*pi)**2/(4*nu*(t + 1)))/(4*nu*(t + 1))     

3. Lambdify：将symbols函数转换为匿名函数，这样就可以通过输入输出值，来获取函数值

from sympy.utilities.lambdify import lambdify

u = -2 * nu * (phiprime / phi) + 4

ufunc = lambdify((t, x, nu), u)
print(ufunc(1, 4, 3))
# output: 3.49170664206

053-求解复杂初始条件+周期性边界条件的一维Burgers方程

1. 有限差分解代码如下：(在前文sympy的基础上)

from matplotlib import pyplot
%matplotlib inline

###variable declarations
nx = 101
nt = 100
dx = 2 * numpy.pi / (nx - 1)
nu = .07
dt = dx * nu

x = numpy.linspace(0, 2 * numpy.pi, nx)
un = numpy.empty(nx)
t = 0

# 计算u
u = numpy.asarray([ufunc(t, x0, nu) for x0 in x])
u

# 画图
pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.plot(x, u, marker='o', lw=2)
pyplot.xlim([0, 2 * numpy.pi])
pyplot.ylim([0, 10]);

2. 方程再展示一下：
   $$u_i^{n+1}=u_i^n-u_i^n\frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n)+\nu \frac{\Delta t}{\Delta x^2}(u_{i+1}^n-2u_i^n+u_{i-1}^n)$$
3. 周期性边界条件求解代码如下：

for n in range(nt):
    un = u.copy()

    # 计算1-N-1的结果
    for i in range(1, nx-1):
        u[i] = un[i] - un[i] * dt / dx *(un[i] - un[i-1]) + nu * dt / dx**2 *(un[i+1] - 2 * un[i] + un[i-1])

    # 计算边界结果
    # 对于i=0的计算结果，取u[-1]=u[0]，即最后一个点的结果=第一个点的结果
    u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + nu * dt / dx**2 *(un[1] - 2 * un[0] + un[-2])
    u[-1] = u[0]
        
u_analytical = numpy.asarray([ufunc(nt * dt, xi, nu) for xi in x])

    # 画图
    pyplot.figure(figsize=(11, 7), dpi=100)
    pyplot.plot(x,u, marker='o', lw=2, label='Computational')
    pyplot.plot(x, u_analytical, label='Analytical')
    pyplot.xlim([0, 2 * numpy.pi])
    pyplot.ylim([0, 10])
    pyplot.legend();