Python氮氧甲烷乙烷乙烯丙烯气体和固体热力学模型计算

🎯要点

🎯固体和粒子：计算二态系统、简谐振子和爱因斯坦固体的内能和比热，比较爱因斯坦固体和德拜固体。模拟多个粒子的一维和二维随机游走，在数值上确认方差的线性趋势，模拟多个粒子的梯度下降，模拟双井的梯度下降。
🎯混合物：🖊建立理想气体及其混合物和多相纯物质的属性模型。
🎯计算给定条件下的气体数值：🖊氮给定条件下的焦耳-汤姆逊系数 | 🖊气态氧流等熵效率的比功和出口气体温度 | 🖊液态水流量泵的冷却负荷和轴功率 | 🖊甲烷与十二烷合并流温度 | 🖊制冷剂经饱和蒸气压缩后的出口压力 | 🖊液氨储存罐安全体积 | 🖊氮气流过阀门哪一部分变成液体 | 🖊二氧化碳和氧气合流绝热压缩和等温压缩的压缩功率 | 🖊乙烯经两阀门膨胀后的温度 | 🖊蒸馏塔中塔压下降时泄漏率会如何变化 | 🖊丙烯储存容器加热后初始压力和最终温度 | 🖊氧气流被压缩机压缩后，压缩机功率和氧气出口温度 | 🖊热力循环水经过恒压加热、等熵膨胀、恒压冷凝和等熵压缩后的热效率 | 🖊计算甲烷焦耳-汤姆逊系数 | 🖊氨被压缩后，每摩尔压缩机的负荷和出口温度 | 🖊蒸汽压缩循环模拟。
🎯Python和C++物理计算热力学 | 🎯Python物理差分方程解

🍇Python水波纹偏微分方程

有一个不可否认的事实：波一直在我们周围。无论是通过电磁辐射、空气中传播的声音，还是由一滴水引起的涟漪，波绝对无处不在。尽管其形式多种多样，但在所有波浪中都适用的一件事是控制它们的严格物理规则。在数学上，所有这些物理规则都浓缩为一个简单的微分方程之一：波动方程。
$\frac{\partial^2 u}{\partial t^2}=c^2\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)$
波动方程告诉我们任何波如何在空间中传播并随时间演化，它为我们提供了函数 u(t, x, y)，该函数给出了任意时间内任意点 (x, y) 处的波的高度 t。本质上，如果你告诉波动方程，你把一个球扔进池塘，它会告诉你球产生的涟漪如何随着时间的推移而发展，并允许你预测由此引起的水的变化。就像你告诉波动方程你尽可能大声尖叫一样，它会告诉你你创建的声波将如何传播并与环境相互作用。总而言之，通过求解波动方程，我们将能够轻松模拟我们的纹波。

想象一下，两个人最初紧紧地握着一根绳子（红色），然后一个人快速上下移动手，在绳子上释放出波浪。在波穿过这条绳子后拍一张快照，我们想弄清楚它将如何随着时间的推移继续演化。为了举例说明，我们假设一根无重力的无摩擦绳索。为了开始解决这个问题，我们首先定义一些变量。

绳索在位置 x 和时间 t 处的高度将被称为函数 u(t, x) 的值。考虑到这一点，我们要考虑什么可能会导致绳索随着时间的推移而发生变化。由于我们确定不存在重力或摩擦力，因此作用在绳子上的唯一力就是张力，我们可以想象，张力在绳子的最小和最大高度处最大。

本质上，这是对凹度的陈述，因为张力在沿着绳索的高弯曲点处达到最高，并且随着远离这些点而减弱。这里最终可以确定的是，在绳索上的任何给定点，它所承受的张力与其凹度成正比。由于微积分告诉我们函数的凹性是由其在空间中的二阶导数定义的，因此我们可以将其表述为：
$F_{\text {tension }}=k \frac{\partial^2 u}{\partial x^2}$
其中 k 是一个简单的比例常数。现在我们必须找到一种方法，将力用 u 表示，并将其转化为我们可以解决和使用的东西。研究牛顿第二定律，它告诉我们 F = ma，允许快速替换，结果为：
$\frac{\partial^2 u}{\partial x^2}$
其中 m 是该点的质量，a 是加速度，u 是高度，k 是比例常数。正如物理学告诉我们的那样，任何处理位置的函数的加速度只是其相对于时间的二阶导数，这意味着我们可以将方程重新定义为：
$\frac{\partial^2 u}{\partial t^2}=\frac{k}{m} \frac{\partial^2 u}{\partial x^2}$
离散化波动方程
$\frac{\partial^2 u}{\partial t^2}=c^2\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)$
在代码中，网格域将如下所示：

Lx = 10 
Nx = 80 
Ly = 10 
Ny = 80 

x_vec = numpy.linspace(0, Lx, Nx) 
dx = x_vec[2] - x_vec[1] 

y_vec = numpy.linspace(0, Ly, Ny)
dy = y_vec[2] - y_vec[1]

dt 可以选择，这就是对我有用的数字，c 保持为 1 以保持方程简单。 dt 越小越好，因为它将导致更准确但更慢的模拟。

dt = .025 
Nt = 4000 
c = 1 
u = numpy.zeros([Nt, len(x), len(y)])

如前所述，我们希望在池的中心出现扰动（中心为 (nx/2, ny/2)），因此 t-dt 和 t 的初始条件如下所示：

u[0, Nx // 2, Ny // 2] = numpy.sin(0) 
u[1, Nx // 2, Ny // 2] = numpy.sin(1/10)

最后，要解决所有问题，我们要做的就是迭代时间并将所有值代入离散方程，从而得到：

for t in range(1, Nt-1):
    for x in range(1, Nx-1):
        for y in range(1, Ny-1):
            if (t < 100):
                u[t, Nx // 2, Ny // 2] = numpy.sin(t / 10)
            u[t+1, x, y] = c**2 * dt**2 * ( ((u[t, x+1, y] - 2*u[t, x, y] + u[t, x-1, y])/(dx**2)) + ((u[t, x, y+1] - 2*u[t, x, y] + u[t, x, y-1])/(dy**2)) ) + 2*u[t, x, y] - u[t-1, x, y]

现在您已经通过求解波动方程正式模拟了水波纹！只需迭代计算值并在代码中为每个步骤及时绘制一个曲面：

fig = pyplot.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y = numpy.meshgrid(x_vec, y_vec)
for t in range(0, Nt):
    surf = ax.plot_surface(X, Y, u[t], color='b', shade=True,
                           linewidth=0, antialiased=False)

    ax.view_init(elev=45)
    ax.set_zlim(-.0001, 2.4)
    pyplot.axis('off')

    pyplot.pause(.0001)
    pyplot.cla()