首页/文章/ 详情

基于Python的CFD编程入门（5）二维线性对流问题

易木木响叮当

1年前浏览476

有了基于Python的CFD编程入门（1）— 一维线性对流问题的基础，很容易拓展至二维问题。

理论基础

根据不可压动量方程：

这次保留速度的x、y方向分量，再令速度的一阶导系数为常数，这样就得到了二维线性对流（2D Linear Convection）方程：

波速依旧为c。

对时间导数采用前向差分（上标n），对空间导数采用后向差分（由于空间为二维平面，因此下标有两个，分别为i、j）离散化对流方程：

调整格式：

初始条件：

在t = 0时：

x	i	y	j	u、v
0	0	0	0	1
0<x≤0.5	0<x≤5	0<x≤0.5	0<x≤5	1
0.5<x≤1	5<x≤10	0.5<x≤1	5<x≤10	2
1<x<2	10<x<20	1<x<2	10<x≤20	1
2	20	2	20	1

边界条件：

t	n	u、v
0<t≤0.5	0≤n≤50	1

代码实现

首先定义求解函数：

def convection(nt, nx, ny, tmax, xmax, ymax, c):
   
   # Increments
   dt = tmax/(nt-1)
   dx = xmax/(nx-1)
   dy = ymax/(ny-1)

   # Initialise data structures
   import numpy as np
   u = np.ones(((nx,ny,nt)))
   v = np.ones(((nx,ny,nt)))
   x = np.zeros(nx)
   y = np.zeros(ny)

   # Boundary conditions
   u[0,:,:] = u[nx-1,:,:] = u[:,0,:] = u[:,ny-1,:] = 1
   v[0,:,:] = v[nx-1,:,:] = v[:,0,:] = v[:,ny-1,:] = 1

   # Initial conditions
   u[int((nx-1)/4):int((nx-1)/2),int((ny-1)/4):int((ny-1)/2),0]=2
   v[int((nx-1)/4):int((nx-1)/2),int((ny-1)/4):int((ny-1)/2),0]=2

   # Loop
   for n in range(0,nt-1):
      for i in range(1,nx-1):
         for j in range(1,ny-1):
            u[i,j,n+1] = (u[i,j,n]-c*dt*((1/dx)*(u[i,j,n]-u[i-1,j,n])+
                                         (1/dy)*(u[i,j,n]-u[i,j-1,n])))
            v[i,j,n+1] = (v[i,j,n]-c*dt*((1/dx)*(v[i,j,n]-v[i-1,j,n])+
                                         (1/dy)*(v[i,j,n]-v[i,j-1,n])))

   # X Loop
   for i in range(0,nx):
      x[i] = i*dx

   # Y Loop
   for j in range(0,ny):
      y[j] = j*dy

   return u, v, x, y

接着定义结果绘图函数：

def plot_3D(u,x,y,time,title,label):
   """
   Plots the 2D velocity field
   """

   import matplotlib.pyplot as plt
   import numpy as np
   fig=plt.figure(figsize=(11,7),dpi=100)
   ax = fig.add_subplot(projection='3d')
   ax.set_xlabel('x (m)')
   ax.set_ylabel('y (m)')
   ax.set_zlabel(label)
   X,Y=np.meshgrid(x,y)
   surf = ax.plot_surface(X,Y,u[:,:,time],rstride=2,cstride=2,cmap='viridis')
   plt.title(title)
   plt.show()

最后调用函数进行运算：

u,v,x,y = convection(101, 81, 81, 0.5, 2.0, 2.0, 0.5)
plot_3D(u,x,y,0,'Figure 1: c=0.5m/s, nt=101, nx=81, ny=81, t=0sec','u (m/s)')
plot_3D(u,x,y,100,'Figure 2: c=0.5m/s, nt=101, nx=81, ny=81, t=0.5sec','u (m/s)')
plot_3D(v,x,y,0,'Figure 3: c=0.5m/s, nt=101, nx=81, ny=81, t=0sec','v (m/s)')
plot_3D(v,x,y,100,'Figure 4: c=0.5m/s, nt=101, nx=81, ny=81, t=0.5sec','v (m/s)')