[系列]Python计算CFD问题
[系列]Python计算CFD问题<1>：一维线性对流问题#-*-coding=utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import time,sys
nx=41 #空间节点数
dx=2.0/(nx-1) #网格间距
nt=25 #时间步数
dt=0.025 #时间步长
c=1
#初始化节点
u=np.ones(nx)
u=2
plt.plot(np.linspace(0,2,nx),u,lw=2)
#离散求解
un=np.ones(nx)
for n in range(nt):
un=u.copy()
for i in range(1,nx):
u=un-c*dt/dx*(un-un)
#显示计算结果
plt.plot(np.linspace(0,2,nx),u,lw=2,color='red')
plt.xlabel("$x$")
plt.ylabel("$u(x)$")
显示的计算结果如图所示：
图中蓝色线条为初始速度分布，红色线条为0.625s（25个时间步，时间步长0.025s）后的速度分布。
利用matplotlib库也很容易制作动画，有兴趣的童鞋可以试着将各时间步上的速度分布以动画的形式输出。 [系列]Python计算CFD问题<2>：一维非线性对流问题
#-*-coding=utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import time,sys
nx=41 #空间节点数
dx=2.0/(nx-1) #网格间距
nt=25 #时间步数
dt=0.025 #时间步长
c=1
#初始化节点
u=np.ones(nx)
u=2
plt.plot(np.linspace(0,2,nx),u,lw=2)
#离散求解
un=np.ones(nx)
for n in range(nt):
un=u.copy()
for i in range(1,nx):
u=un-un*dt/dx*(un-un)
#显示计算结果
plt.plot(np.linspace(0,2,nx),u,lw=2,color='red')
plt.xlabel("$x$")
plt.ylabel("$u(x)$")
显示的计算结果如图所示：
图中蓝色线条为初始速度分布，红色线条为0.625s（25个时间步，时间步长0.025s）后的速度分布。
【来源：http://nbviewer.ipython.org/gith ... ons/02_Step_2.ipynb】 本帖最后由 funi 于 2015-10-29 21:09 编辑
[系列]Python计算CFD问题<3>：关于CFL
对于前面的关于1D对流问题，计算过程中，采用的网格数为41，时间步数为25，时间步长为0.025，对于给定的条件，此参数配置能够获得较好的计算结果。那么这些计算参数能否影响最终计算结果？比如说增大或减小网格数量？增大或减小时间步长？下面来测试一下。
我们采用案例1的条件进行测试。修改程序如下：
import numpy as np
import matplotlib.pyplot as plt
def linearconv(nx):
dx=2.0/(nx-1)
nt=20
dt=0.025
c=1
u=np.ones(nx)
u[.5/dx:1/dx+1]=2
un=np.ones(nx)
for n in range(nt):
un=u.copy()
for i in range(1,nx):
u=un-c*dt/dx*(un-un)
plt.plot(np.linspace(0,2,nx),u,lw=2,color='red')
这里函数的形式，以方便后面的比较。这里分别比较nx=41、61、81、101、121时候，ux分布情况。
看出问题了没？
随着NX增大，即网格的加密，计算结果似乎越来越离谱了？是不是有点儿颠覆了俺们的常规认识了呢？这就是计算稳定性所导致的问题。为了找到导致突变的网格数，我们从80开始增大nx。从图中可以看出，在nx=80的时候，u(x)的分布是可以接受的，我们从80开始慢慢增加nx。
当nx增大到82时，情况发生突变，说明此处是一个坎儿。
由于计算中采用了固定的时间步长，因此随着网格的加密，空间步长与时间步长的比值将会越来越小。这里定义库朗数：
式中u为波速，在本例中u=1；σ为库朗数；σmax为所采用的离散算法确保稳定性的值。对于本例采用的离散算法，σmax=1，可以计算得nx最大值为81。
据此，可以修改计算程序，通过库朗数自动调整时间步长，以保证计算过程稳定。
修改后的计算程序为：
#-*- coding=utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
def linearconv(nx):
dx=2.0/(nx-1)
nt=20
c=1
sigma = 0.5 #库朗数
dt=sigma*dx/c #计算时间步长
u=np.ones(nx)
u[.5/dx:1/dx+1]=2
un=np.ones(nx)
for n in range(nt):
un=u.copy()
for i in range(1,nx):
u=un-c*dt/dx*(un-un)
plt.plot(np.linspace(0,2,nx),u,lw=2,color='red')
linearconv(80)
计算nx超过80后的ux曲线，如图所示。
从图中可以看出，通过库朗数定义计算时间步长，可以有效的避免计算不稳定情况。上面各图中因为采用的NX不同，所以导致时间步长有差异，在统一的时间步数情况下，最终计算时长存在差异，有感兴趣的童鞋可以通过调整时间步数以使它们的时间点一致。
总结：对于显式非稳态问题，库朗数是控制其计算稳定性的重要参数。减小时间步长或增大网格都有助于计算稳定。
【来源：http://nbviewer.ipython.org/gith ... CFL_Condition.ipynb】 [系列]Python计算CFD问题<4>：一维扩散方程
与对流方程不同，扩散方程是2阶的偏微分方程。一维扩散方程可用下面的方程进行描述：
该方程可以离散为：
写成迭代的形式为：
采用与前例相同的初始条件，及0.5<=x<=1时，u=1，其他位置u=2。扩散系数mu=0.3。
程序代码为：
利用函数diffision可查看各时刻速度分布（时间步数分别为0、50、100、150、200、250）。
放到一起看起来可能更直观一些：
【来源：hhttp://nbviewer.ipython.org/github/barbagroup/CFDPython/blob/master/lessons/04_Step_3.ipynb】 [系列]Python计算CFD问题<5>：一维Burgers方程
一维Burgers方程形式为：
从方程形式上来看，该方程是前面的非线性对流方程与扩散方程的混合体。采用前面提到的离散方法，可以将Burgers方程离散成：
将其写成迭代的形式为：
除了迭代形式外，我们还必须知道初始条件和边界条件。
该方程有解析解：
因此本次计算设初始条件为：
设置边界条件为：
初始条件中含有偏导项，无法直接赋值，需要采用一些数学化简手段进行处理。
Python中提供了sympy库，该库是一个符号计算库，可以帮助进行化简。
完整的程序如下：
# -*-coding:utf-8 -*-
import numpy as np
import sympy
from sympy.utilities.lambdify import lambdify
import matplotlib.pylab as plt
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*np.pi)**2/(4*nu*(t+1)))
phiprime = phi.diff(x)
u = -2*nu*(phiprime/phi)+4
ufunc = lambdify((t,x,nu),u)
nx = 101
nt=100
dx = 2*np.pi/(nx-1)
nu=0.07
dt=dx*nu
x= np.linspace(0,2*np.pi,nx)
un = np.empty(nx)
t=0
u=np.asarray() #list 转化成 np.array
plt.figure(figsize=(4,4),dpi=100)
plt.plot(x,u,lw=2)
plt.xlim()
plt.ylim()
for n in range(nt):
un = u.copy()
for i in range(nx-1):
u = un - un * dt/dx *(un - un) + nu*dt/dx**2*(un-2*un+un)
u[-1] = un[-1] - un[-1] * dt/dx * (un[-1] - un[-2]) + nu*dt/dx**2*(un-2*un[-1]+un[-2])
u_analytical = np.asarray()
plt.figure(figsize=(6,6),dpi=100)
plt.plot(x,u,marker="o",color="blue",lw=2,label='Computational')
plt.plot(x,u_analytical,color="yellow",label='analytical')
plt.xlim()
plt.ylim()
plt.legend()
plt.show()
初始条件如图：
计算结果如图：
【来源：http://nbviewer.ipython.org/github/barbagroup/CFDPython/blob/master/lessons/05_Step_4.ipynb#Saving-Time-with-SymPy】 [系列]Python计算CFD问题<6>：2D线性对流问题
在第一个例子中，我们提到了解一维线性对流问题。在本例中，我们将利用python数值求解二维线性对流问题。
2D线性对流问题可以用此方程进行表达：
时间项采用向前差分，空间项采用先后差分。离散后的方程为：
换成迭代格式为：
初始条件为：
边界条件：
编制计算程序为：
from mpl_toolkits.mplot3d import Axes3D ##New Library required for projected 3d plots
import numpy as np
import matplotlib.pyplot as plt
###variable declarations
nx = 81
ny = 81
nt = 100
c = 1
dx = 2.0/(nx-1)
dy = 2.0/(ny-1)
sigma = .2
dt = sigma*dx
x = np.linspace(0,2,nx)
y = np.linspace(0,2,ny)
u = np.ones((ny,nx)) ##create a 1xn vector of 1's
un = np.ones((ny,nx)) ##
###Assign initial conditions
u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2
###Plot Initial Condition
fig = plt.figure(figsize=(11,7), dpi=100)
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(x,y)
surf = ax.plot_surface(X,Y,u[:])
u = np.ones((ny,nx))
u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2
for n in range(nt+1): ##loop across number of time steps
un = u.copy()
for i in range(1, len(u)):
for j in range(1, len(u)):
u = un - (c*dt/dx*(un - un))-(c*dt/dy*(un-un))
u = 1
u[-1,:] = 1
u[:,0] = 1
u[:,-1] = 1
fig = plt.figure(figsize=(11,7), dpi=100)
ax = fig.gca(projection='3d')
surf2 = ax.plot_surface(X,Y,u[:])
输出结果：
初始条件：
最终结果：
做CFD计算的童鞋可能更容易理解云图，这里以云图的方式显示：
只需修改图形显示代码：
plt.figure()
im=plt.imshow(u,interpolation='bilinear',origin='lower',cmap=cm.rainbow)
plt.contour(X,Y,u)
plt.colorbar(im,orientation='vertical')
【来源：http://nbviewer.ipython.org/github/barbagroup/CFDPython/blob/master/lessons/07_Step_5.ipynb】 [系列]Python计算CFD问题<7>：2D对流问题
在前面的例子中，我们处理了2D线性对流问题，虽然问题是二维的，但实际上还是一维的，因为只有计算域是2D，物理变量只有X方向速度是变量。在本例中，我们来处理真正的2D对流问题。
2D对流问题可以用此方程组进行表达：
时间项采用向前差分，空间项采用先后差分。离散后的方程为：
换成迭代格式为：
初始条件为：
边界条件：
编制计算程序为：
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.cm as cm
nx=101
ny=101
nt=80
dx=2.0/(nx-1)
dy=2.0/(ny-1)
sigma=0.2
dt=sigma*dx
x=np.linspace(0,2,nx)
y=np.linspace(0,2,ny)
u=np.ones((nx,ny))
v=np.ones((nx,ny))
un=np.ones((nx,ny))
vn=np.ones((nx,ny))
u[.5/dx:1/dx+1,.5/dy:1/dy+1]=2
v[.5/dx:1/dx+1,.5/dy:1/dy+1]=2
for n in range(nt+1):
un=u.copy()
vn=v.copy()
u=un-(un*dt/dx*(un-un))-vn*dt/dy*(un-un)
v=vn-(un*dt/dx*(vn-vn))-vn*dt/dy*(vn-vn)
u=1
u[-1,:]=1
u[:,0]=1
u[:,-1]=1
v=1
v[-1,:]=1
v[:,0]=1
v[:,-1]=1
X,Y=np.meshgrid(x,y)
plt.figure()
im=plt.imshow(u,interpolation='bilinear',origin='lower',cmap=cm.rainbow)
plt.contour(X,Y,v)
plt.colorbar(im,orientation='vertical')
输出结果：
初始条件（u，v分布）
最终结果（u,v分布）：
【来源：http://nbviewer.ipython.org/github/barbagroup/CFDPython/blob/master/lessons/07_Step_5.ipynb】 大拇指，厉害。先收藏了。
页:
[1]