🍍Python湍流隐式模型耗散粘性方程和大涡流模拟

Python | 物理 | 数学 | 偏微分方程 | 波动方程 | 二维变速 | 算子 | 双曲波 | 耗散系统 | 粘性伯格斯方程 | 算法 | 傅里叶 | 线性和非线性 | 对流扩散 | 湍流 | 模型

Python湍流隐式模型耗散粘性方程和大涡流模拟 | 亚图跨际亚图跨际 on Notion

🏈指点迷津 | Brief

🎯要点

🎯达朗贝尔一维波动通解，二维变速模拟 | 🎯达朗贝尔算子解双曲波形微分方程 | 🎯耗散系统粘性伯格斯方程快速傅里叶变换算法 | 🎯二维线性和非线性对流扩散解和湍流隐式建模

📜偏微分方程用例：Python自动造波器椭圆曲线波孤子解

📜有限差分用例：Python微磁学磁倾斜和西塔规则算法

🍇Python一维粘性伯格斯方程

一维空间中粘性伯格斯方程的一般形式是耗散系统：

$$\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}=\nu \frac{\partial^2 u}{\partial x^2}$$

此项 $u \partial u / \partial x$ 也可以重写为$\partial\left(u^2 / 2\right) / \partial x$。当扩散项不存在时（即 $\nu=0$​ ），粘性伯格斯程变为无粘伯格斯方程：

$$\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}=0$$

这是可以产生不连续性（冲击波）的守恒方程的原型。

当人们检查等式的左侧时，$\nu$ 的小值形成锐梯度的原因就变得直观清楚了。此项 $\partial / \partial t+u \partial / \partial x$ 显然是一个波算子，描述以 $u$ 速度沿正 $x$ 方向传播的波。由于波速为$u$，表现出较大$u$值的区域将比表现出较小u值的区域更快地向右传播；换句话说，如果 $u$ 最初沿$x$方向减小，则位于背面的较大 $u$ 将赶上位于正面的较小 $u$ 。右侧扩散项的作用本质上是阻止梯度变得无穷大。

在此，我们将使用非线性对流和扩散，仅创建一维伯格方程。

$$\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}=\nu \frac{\partial^2 u}{\partial x^2}$$

我们可以离散化这个微分，采用以下形式：

$$\frac{u_i^{n+1}-u_i^n}{\Delta t}+u_i^n \frac{u_i^n-u_{i-1}^n}{\Delta x}=v \frac{u_{i+1}^n-2 u_i^n+u_{i-1}^n}{\Delta x^2}$$

在解决了未知数之后，我们得到了用 Python 编码的算法。

$$u_i^{n+1}=u_i^n-u_i^n \frac{\Delta t}{\Delta x}\left(u_i^n-u_{i-1}^n\right)+v \frac{\Delta t}{\Delta x^2}\left(u_{i+1}^n-2 u_i^n+u_{i-1}^n\right)$$

速度的初始条件是使用以下函数创建的：

$$\begin{aligned} u & =-\frac{2 v}{\phi} \frac{\partial \phi}{\partial x }+4 \\ \phi & =\exp \left(\frac{- x ^2}{4 v}\right)+\exp \left(\frac{-( x -2 \pi)^2}{4 v}\right) \end{aligned}$$

边界条件意味着周期性，由下式给出：

$$u(0)=u(2 \pi)$$

有一个解析解：

$$\begin{aligned} & u =-\frac{2 v}{\phi} \frac{\partial \phi}{\partial x }+4 \\ & \phi=\exp \left(\frac{-( x -4 t )^2}{4 v( t +1)}\right)+\exp \left(\frac{-( x -4 t -2 \pi)^2}{4 v( t +1)}\right) \end{aligned}$$

解该问题的Python代码如下：

 import numpy as np
 import sympy as sp
 import pylab as pl
 pl.ion()
 ​
 x, nu, t = sp.symbols('x nu t')
 phi = sp.exp(-(x-4*t)**2/(4*nu*(t+1))) + sp.exp(-(x-4*t-2*np.pi)**2/(4*nu*(t+1)))
 ​
 phiprime = phi.diff(x)
 u = -2*nu*(phiprime/phi)+4
 ​
 from sympy.utilities.lambdify import lambdify
 ufunc = lambdify ((t, x, nu), u)
 ​
 nx = 101 
 nt = 100 
 dx = 2*np.pi/(nx-1) 
 nu = 0.07 
 dt = dx*nu 
 T = nt*dt 
 ​
 grid = np.linspace(0, 2*np.pi, nx) 
 un = np.empty(nx) 
 t = 0
 u = np.asarray([ufunc(t, x, nu) for x in grid])
 ​
 pl.figure(figsize=(11,7), dpi=100)
 pl.plot(grid,u, marker='o', lw=2)
 pl.xlim([0,2*np.pi])
 pl.ylim([0,10])
 pl.xlabel('X')
 pl.ylabel('Velocity') 
 pl.title('1D Burgers Equation - Initial condition')
 ​
 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(T, xi, nu) for xi in grid])
 ​
 pl.figure(figsize=(11,7), dpi=100)
 pl.plot(grid, u, marker='o', lw=2, label='Computational')
 pl.plot(grid, u_analytical, label='Analytical')
 pl.xlim([0, 2*np.pi])
 pl.ylim([0,10])
 pl.legend()
 pl.xlabel('X')
 pl.ylabel('Velocity') 
 pl.title('1D Burgers Equation - Solutions')