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

🎯要点

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

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

📜有限差分用例：Python微磁学磁倾斜和西塔规则算法
在这里插入图片描述

🍇Python一维粘性伯格斯方程

一维空间中粘性伯格斯方程的一般形式是耗散系统：
$\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}=\nu \frac{\partial^2 u}{\partial x^2}$
此项 $\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}$
边界条件意味着周期性，由下式给出：
$\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')