[系列]Python计算CFD问题

funi

[系列]Python计算CFD问题<1>：一维线性对流问题

1. #-*-coding=utf-8 -*-
   
2. import numpy as np
   
3. import matplotlib.pyplot as plt
   
4. import time,sys
   
5. 
   
6. nx=41 #空间节点数
   
7. dx=2.0/(nx-1) #网格间距
   
8. nt=25 #时间步数
   
9. dt=0.025 #时间步长
   
10. c=1
    
11. 
    
12. #初始化节点
    
13. u=np.ones(nx)
    
14. u[0.5/dx:1/dx+1]=2
    
15. plt.plot(np.linspace(0,2,nx),u,lw=2)
    
16. 
    
17. #离散求解
    
18. un=np.ones(nx)
    
19. for n in range(nt):
    
20.     un=u.copy()
    
21.     for i in range(1,nx):
    
22.         u[i]=un[i]-c*dt/dx*(un[i]-un[i-1])
    
23. 
    
24. #显示计算结果
    
25. plt.plot(np.linspace(0,2,nx),u,lw=2,color='red')
    
26. plt.xlabel("$x$")
    
27. plt.ylabel("$u(x)$")

复制代码

显示的计算结果如图所示：

图中蓝色线条为初始速度分布，红色线条为0.625s（25个时间步，时间步长0.025s）后的速度分布。
利用matplotlib库也很容易制作动画，有兴趣的童鞋可以试着将各时间步上的速度分布以动画的形式输出。

funi

[系列]Python计算CFD问题<2>：一维非线性对流问题

1. #-*-coding=utf-8 -*-
   
2. import numpy as np
   
3. import matplotlib.pyplot as plt
   
4. import time,sys
   
5. 
   
6. nx=41 #空间节点数
   
7. dx=2.0/(nx-1) #网格间距
   
8. nt=25 #时间步数
   
9. dt=0.025 #时间步长
   
10. c=1
    
11. 
    
12. #初始化节点
    
13. u=np.ones(nx)
    
14. u[0.5/dx:1/dx+1]=2
    
15. plt.plot(np.linspace(0,2,nx),u,lw=2)
    
16. 
    
17. #离散求解
    
18. un=np.ones(nx)
    
19. for n in range(nt):
    
20.     un=u.copy()
    
21.     for i in range(1,nx):
    
22.         u[i]=un[i]-un[i]*dt/dx*(un[i]-un[i-1])
    
23. 
    
24. #显示计算结果
    
25. plt.plot(np.linspace(0,2,nx),u,lw=2,color='red')
    
26. plt.xlabel("$x$")
    
27. plt.ylabel("$u(x)$")

复制代码

显示的计算结果如图所示：

图中蓝色线条为初始速度分布，红色线条为0.625s（25个时间步，时间步长0.025s）后的速度分布。
【来源：http://nbviewer.ipython.org/gith ... ons/02_Step_2.ipynb】

funi

[系列]Python计算CFD问题<3>：关于CFL
对于前面的关于1D对流问题，计算过程中，采用的网格数为41，时间步数为25，时间步长为0.025，对于给定的条件，此参数配置能够获得较好的计算结果。那么这些计算参数能否影响最终计算结果？比如说增大或减小网格数量？增大或减小时间步长？下面来测试一下。
我们采用案例1的条件进行测试。修改程序如下：

1. import numpy as np
   
2. import matplotlib.pyplot as plt
   
3. 
   
4. def linearconv(nx):
   
5. dx=2.0/(nx-1)
   
6. nt=20
   
7. dt=0.025
   
8. c=1
   
9. 
   
10. u=np.ones(nx)
    
11. u[.5/dx:1/dx+1]=2
    
12. un=np.ones(nx)
    
13. 
    
14. for n in range(nt):
    
15. un=u.copy()
    
16. for i in range(1,nx):
    
17. u[i]=un[i]-c*dt/dx*(un[i]-un[i-1])
    
18. 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。
据此，可以修改计算程序，通过库朗数自动调整时间步长，以保证计算过程稳定。
修改后的计算程序为：

1. #-*- coding=utf-8 -*-
   
2. import numpy as np
   
3. import matplotlib.pyplot as plt
   
4. 
   
5. def linearconv(nx):
   
6. dx=2.0/(nx-1)
   
7. nt=20
   
8. c=1
   
9. sigma = 0.5 #库朗数
   
10. dt=sigma*dx/c #计算时间步长
    
11. 
    
12. u=np.ones(nx)
    
13. u[.5/dx:1/dx+1]=2
    
14. un=np.ones(nx)
    
15. 
    
16. for n in range(nt):
    
17. un=u.copy()
    
18. for i in range(1,nx):
    
19. u[i]=un[i]-c*dt/dx*(un[i]-un[i-1])
    
20. plt.plot(np.linspace(0,2,nx),u,lw=2,color='red')
    
21. linearconv(80)

复制代码

计算nx超过80后的ux曲线，如图所示。






从图中可以看出，通过库朗数定义计算时间步长，可以有效的避免计算不稳定情况。上面各图中因为采用的NX不同，所以导致时间步长有差异，在统一的时间步数情况下，最终计算时长存在差异，有感兴趣的童鞋可以通过调整时间步数以使它们的时间点一致。
总结：对于显式非稳态问题，库朗数是控制其计算稳定性的重要参数。减小时间步长或增大网格都有助于计算稳定。
【来源：http://nbviewer.ipython.org/gith ... CFL_Condition.ipynb】

funi

[系列]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/gith ... ons/04_Step_3.ipynb】

funi

[系列]Python计算CFD问题<5>：一维Burgers方程
一维Burgers方程形式为：

从方程形式上来看，该方程是前面的非线性对流方程与扩散方程的混合体。采用前面提到的离散方法，可以将Burgers方程离散成：
将其写成迭代的形式为：

除了迭代形式外，我们还必须知道初始条件和边界条件。
该方程有解析解：

因此本次计算设初始条件为：

设置边界条件为：

初始条件中含有偏导项，无法直接赋值，需要采用一些数学化简手段进行处理。
Python中提供了sympy库，该库是一个符号计算库，可以帮助进行化简。
完整的程序如下：

1. # -*-coding:utf-8 -*-
   
2. import numpy as np
   
3. import sympy
   
4. from sympy.utilities.lambdify import lambdify
   
5. import matplotlib.pylab as plt
   
6. 
   
7. x,nu,t = sympy.symbols("x,nu,t")
   
8. phi = sympy.exp(-(x-4*t)**2/(4*nu*(t+1))) + sympy.exp(-(x-4*t-2*np.pi)**2/(4*nu*(t+1)))
   
9. phiprime = phi.diff(x)
   
10. u = -2*nu*(phiprime/phi)+4
    
11. ufunc = lambdify((t,x,nu),u)
    
12. 
    
13. nx = 101
    
14. nt=100
    
15. dx = 2*np.pi/(nx-1)
    
16. nu=0.07
    
17. dt=dx*nu
    
18. 
    
19. x= np.linspace(0,2*np.pi,nx)
    
20. un = np.empty(nx)
    
21. t=0
    
22. u=np.asarray([ufunc(t,x0,nu) for x0 in x]) #list 转化成 np.array
    
23. 
    
24. plt.figure(figsize=(4,4),dpi=100)
    
25. plt.plot(x,u,lw=2)
    
26. plt.xlim([0,2*np.pi])
    
27. plt.ylim([0,8])
    
28. 
    
29. for n in range(nt):
    
30. un = u.copy()
    
31. for i in range(nx-1):
    
32. 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])
    
33. u[-1] = un[-1] - un[-1] * dt/dx * (un[-1] - un[-2]) + nu*dt/dx**2*(un[0]-2*un[-1]+un[-2])
    
34. 
    
35. u_analytical = np.asarray([ufunc(nt*dt,xi,nu) for xi in x])
    
36. 
    
37. plt.figure(figsize=(6,6),dpi=100)
    
38. plt.plot(x,u,marker="o",color="blue",lw=2,label='Computational')
    
39. plt.plot(x,u_analytical,color="yellow",label='analytical')
    
40. plt.xlim([0,2*np.pi])
    
41. plt.ylim([0,10])
    
42. plt.legend()
    
43. plt.show()

复制代码

初始条件如图：

计算结果如图：

【来源：http://nbviewer.ipython.org/gith ... ing-Time-with-SymPy】

funi

[系列]Python计算CFD问题<6>：2D线性对流问题
在第一个例子中，我们提到了解一维线性对流问题。在本例中，我们将利用python数值求解二维线性对流问题。
2D线性对流问题可以用此方程进行表达：

时间项采用向前差分，空间项采用先后差分。离散后的方程为：

换成迭代格式为：

初始条件为：

边界条件：

编制计算程序为：

1. from mpl_toolkits.mplot3d import Axes3D ##New Library required for projected 3d plots
   
2. 
   
3. import numpy as np
   
4. import matplotlib.pyplot as plt
   
5. 
   
6. ###variable declarations
   
7. nx = 81
   
8. ny = 81
   
9. nt = 100
   
10. c = 1
    
11. dx = 2.0/(nx-1)
    
12. dy = 2.0/(ny-1)
    
13. sigma = .2
    
14. dt = sigma*dx
    
15. 
    
16. x = np.linspace(0,2,nx)
    
17. y = np.linspace(0,2,ny)
    
18. 
    
19. u = np.ones((ny,nx)) ##create a 1xn vector of 1's
    
20. un = np.ones((ny,nx)) ##
    
21. 
    
22. ###Assign initial conditions
    
23. u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2
    
24. ###Plot Initial Condition
    
25. fig = plt.figure(figsize=(11,7), dpi=100)
    
26. ax = fig.gca(projection='3d')
    
27. X, Y = np.meshgrid(x,y)
    
28. surf = ax.plot_surface(X,Y,u[:])
    
29. 
    
30. u = np.ones((ny,nx))
    
31. u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2
    
32. 
    
33. for n in range(nt+1): ##loop across number of time steps
    
34. un = u.copy()
    
35. for i in range(1, len(u)):
    
36. for j in range(1, len(u)):
    
37. u[i,j] = un[i, j] - (c*dt/dx*(un[i,j] - un[i-1,j]))-(c*dt/dy*(un[i,j]-un[i,j-1]))
    
38. u[0,:] = 1
    
39. u[-1,:] = 1
    
40. u[:,0] = 1
    
41. u[:,-1] = 1
    
42. fig = plt.figure(figsize=(11,7), dpi=100)
    
43. ax = fig.gca(projection='3d')
    
44. surf2 = ax.plot_surface(X,Y,u[:])

复制代码

输出结果：
初始条件：

最终结果：

做CFD计算的童鞋可能更容易理解云图，这里以云图的方式显示：


只需修改图形显示代码：

1. plt.figure()
   
2. im=plt.imshow(u,interpolation='bilinear',origin='lower',cmap=cm.rainbow)
   
3. plt.contour(X,Y,u)
   
4. plt.colorbar(im,orientation='vertical')

复制代码

【来源：http://nbviewer.ipython.org/gith ... ons/07_Step_5.ipynb】

funi

[系列]Python计算CFD问题<7>：2D对流问题
在前面的例子中，我们处理了2D线性对流问题，虽然问题是二维的，但实际上还是一维的，因为只有计算域是2D，物理变量只有X方向速度是变量。在本例中，我们来处理真正的2D对流问题。
2D对流问题可以用此方程组进行表达：

时间项采用向前差分，空间项采用先后差分。离散后的方程为：

换成迭代格式为：

初始条件为：
边界条件：

编制计算程序为：

1. import matplotlib.pyplot as plt
   
2. import numpy as np
   
3. import matplotlib.cm as cm
   
4. 
   
5. nx=101
   
6. ny=101
   
7. nt=80
   
8. 
   
9. dx=2.0/(nx-1)
   
10. dy=2.0/(ny-1)
    
11. sigma=0.2
    
12. dt=sigma*dx
    
13. 
    
14. x=np.linspace(0,2,nx)
    
15. y=np.linspace(0,2,ny)
    
16. 
    
17. u=np.ones((nx,ny))
    
18. v=np.ones((nx,ny))
    
19. un=np.ones((nx,ny))
    
20. vn=np.ones((nx,ny))
    
21. 
    
22. u[.5/dx:1/dx+1,.5/dy:1/dy+1]=2
    
23. v[.5/dx:1/dx+1,.5/dy:1/dy+1]=2
    
24. 
    
25. for n in range(nt+1):
    
26. un=u.copy()
    
27. vn=v.copy()
    
28. u[1:,1:]=un[1:,1:]-(un[1:,1:]*dt/dx*(un[1:,1:]-un[0:-1,1:]))-vn[1:,1:]*dt/dy*(un[1:,1:]-un[1:,0:-1])
    
29. v[1:,1:]=vn[1:,1:]-(un[1:,1:]*dt/dx*(vn[1:,1:]-vn[0:-1,1:]))-vn[1:,1:]*dt/dy*(vn[1:,1:]-vn[1:,0:-1])
    
30. 
    
31. u[0,:]=1
    
32. u[-1,:]=1
    
33. u[:,0]=1
    
34. u[:,-1]=1
    
35. 
    
36. v[0,:]=1
    
37. v[-1,:]=1
    
38. v[:,0]=1
    
39. v[:,-1]=1
    
40. 
    
41. X,Y=np.meshgrid(x,y)
    
42. plt.figure()
    
43. im=plt.imshow(u,interpolation='bilinear',origin='lower',cmap=cm.rainbow)
    
44. plt.contour(X,Y,v)
    
45. plt.colorbar(im,orientation='vertical')

复制代码

输出结果：
初始条件（u，v分布）


最终结果（u,v分布）：


【来源：http://nbviewer.ipython.org/gith ... ons/07_Step_5.ipynb】

fengshao2004

大拇指，厉害。先收藏了。