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

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

{% embed url="https://viadean.notion.site/Python-cb4fdb803fc7450fa939585a2f231012" %}

{% content-ref url="../zhi-dian-mi-jin-brief" %}
[zhi-dian-mi-jin-brief](../zhi-dian-mi-jin-brief)
{% endcontent-ref %}

### 🎯要点

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

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

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

{% @mermaid/diagram content="pie title 语言分比
 "Python":90
 "C++":40" %}

{% @mermaid/diagram content="pie title 内容分比
 "物理":90
 "数学":89
 "偏微分方程":80
 "流体力学":70
 "算法":40
 "对流扩散":20
 "波形传播":10
 "数值模拟":40" %}

### 🍇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代码如下：

```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')
```
