[系列]Python计算CFD问题<5>：一维Burgers方程
一维Burgers方程形式为：
[attach]69955[/attach]
从方程形式上来看，该方程是前面的非线性对流方程与扩散方程的混合体。采用前面提到的离散方法，可以将Burgers方程离散成：[attach]69956[/attach]
将其写成迭代的形式为：
[attach]69957[/attach]
除了迭代形式外，我们还必须知道初始条件和边界条件。
该方程有解析解：
[attach]69958[/attach]
因此本次计算设初始条件为：
[attach]69959[/attach]
设置边界条件为：
[attach]69960[/attach]
初始条件中含有偏导项，无法直接赋值，需要采用一些数学化简手段进行处理。
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([ufunc(t,x0,nu) for x0 in x]) #list 转化成 np.array
-
- plt.figure(figsize=(4,4),dpi=100)
- plt.plot(x,u,lw=2)
- plt.xlim([0,2*np.pi])
- plt.ylim([0,8])
-
- for n in range(nt):
- un = u.copy()
- for i in range(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])
- u[-1] = un[-1] - un[-1] * dt/dx * (un[-1] - un[-2]) + nu*dt/dx**2*(un[0]-2*un[-1]+un[-2])
-
- u_analytical = np.asarray([ufunc(nt*dt,xi,nu) for xi in x])
-
- 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([0,2*np.pi])
- plt.ylim([0,10])
- plt.legend()
- plt.show()
初始条件如图：
[attach]69961[/attach]
计算结果如图：
[attach]69962[/attach]
【来源：http://nbviewer.ipython.org/gith ... ing-Time-with-SymPy】
|
