🎯要点
- 非等压电池放电物理模型模拟高放电对流现象。
- 配置电池和简化化学性质分析，分析液体热力学化学过程及数学计算方式。
- 使用无量纲化系统确定系统中的主导效应，构建模型数值解析。
- 推导出模型系统的三个解析近似解，并在接近小极限的情况下进行渐近分析。
🍁铅蓄电池分析方式
🍪语言内容分比
🍇Python双曲偏微分方程
浅水方程是一组双曲偏微分方程，用于描述流体中压力面以下的流动（有时但不一定是自由表面）。浅水方程是从质量守恒定律和线性动量守恒定律（纳维-斯托克斯方程）推导出来的，即使浅水假设不成立（例如跨越水跃），该方程仍然成立。在水平床的情况下，科里奥利力、摩擦力和粘性力可以忽略不计，浅水方程为：
∂
(
ρ
η
)
∂
t
+
∂
(
ρ
η
u
)
∂
x
+
∂
(
ρ
η
v
)
∂
y
=
0
∂
(
ρ
η
u
)
∂
t
+
∂
∂
x
(
ρ
η
u
2
+
1
2
ρ
g
η
2
)
+
∂
(
ρ
η
u
v
)
∂
y
=
0
∂
(
ρ
η
v
)
∂
t
+
∂
∂
y
(
ρ
η
v
2
+
1
2
ρ
g
η
2
)
+
∂
(
ρ
η
u
v
)
∂
x
=
0.
\begin{aligned} & \frac{\partial(\rho \eta)}{\partial t}+\frac{\partial(\rho \eta u)}{\partial x}+\frac{\partial(\rho \eta v)}{\partial y}=0 \\ & \frac{\partial(\rho \eta u)}{\partial t}+\frac{\partial}{\partial x}\left(\rho \eta u^2+\frac{1}{2} \rho g \eta^2\right)+\frac{\partial(\rho \eta u v)}{\partial y}=0 \\ & \frac{\partial(\rho \eta v)}{\partial t}+\frac{\partial}{\partial y}\left(\rho \eta v^2+\frac{1}{2} \rho g \eta^2\right)+\frac{\partial(\rho \eta u v)}{\partial x}=0 . \end{aligned}
∂t∂(ρη)+∂x∂(ρηu)+∂y∂(ρηv)=0∂t∂(ρηu)+∂x∂(ρηu2+21ρgη2)+∂y∂(ρηuv)=0∂t∂(ρηv)+∂y∂(ρηv2+21ρgη2)+∂x∂(ρηuv)=0.
使用乘积法则展开上面的导数，得到浅水方程的非保守形式。由于速度不受基本守恒方程的影响，非保守形式在冲击或水跃中不成立。还包括科里奥利力、摩擦力和粘性力的适当项，以获得（对于恒定流体密度）：
∂
h
∂
t
+
∂
∂
x
(
(
H
+
h
)
u
)
+
∂
∂
y
(
(
H
+
h
)
v
)
=
0
∂
u
∂
t
+
u
∂
u
∂
x
+
v
∂
u
∂
y
−
f
v
=
−
g
∂
h
∂
x
−
k
u
+
ν
(
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
)
∂
v
∂
t
+
u
∂
v
∂
x
+
v
∂
v
∂
y
+
f
u
=
−
g
∂
h
∂
y
−
k
v
+
ν
(
∂
2
v
∂
x
2
+
∂
2
v
∂
y
2
)
\begin{aligned} & \frac{\partial h}{\partial t}+\frac{\partial}{\partial x}((H+h) u)+\frac{\partial}{\partial y}((H+h) v)=0 \\ & \frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}-f v=-g \frac{\partial h}{\partial x}-k u+\nu\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right) \\ & \frac{\partial v}{\partial t}+u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+f u=-g \frac{\partial h}{\partial y}-k v+\nu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right) \end{aligned}
∂t∂h+∂x∂((H+h)u)+∂y∂((H+h)v)=0∂t∂u+u∂x∂u+v∂y∂u−fv=−g∂x∂h−ku+ν(∂x2∂2u+∂y2∂2u)∂t∂v+u∂x∂v+v∂y∂v+fu=−g∂y∂h−kv+ν(∂x2∂2v+∂y2∂2v)
通常情况下，u 和 v 的二次项（表示整体平流的影响）与其他项相比较小。这称为地转平衡，相当于说罗斯贝数很小。假设波高与平均高度相比非常小（h ≪ H），我们有（没有侧向粘性力）：
∂
h
∂
t
+
H
(
∂
u
∂
x
+
∂
v
∂
y
)
=
0
∂
u
∂
t
−
f
v
=
−
g
∂
h
∂
x
−
k
u
∂
v
∂
t
+
f
u
=
−
g
∂
h
∂
y
−
k
v
\begin{aligned} & \frac{\partial h}{\partial t}+H\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0 \\ & \frac{\partial u}{\partial t}-f v=-g \frac{\partial h}{\partial x}-k u \\ & \frac{\partial v}{\partial t}+f u=-g \frac{\partial h}{\partial y}-k v \end{aligned}
∂t∂h+H(∂x∂u+∂y∂v)=0∂t∂u−fv=−g∂x∂h−ku∂t∂v+fu=−g∂y∂h−kv
Python示例计算
import time
import numpy as np
import matplotlib.pyplot as plt
物理参数
L_x = 1E+6
L_y = 1E+6
g = 9.81
H = 100
f_0 = 1E-4
beta = 2E-11
rho_0 = 1024.0
tau_0 = 0.1
use_coriolis = True
use_friction = False
use_wind = False
use_beta = True
use_source = False
use_sink = False
param_string = "\n================================================================"
param_string += "\nuse_coriolis = {}\nuse_beta = {}".format(use_coriolis, use_beta)
param_string += "\nuse_friction = {}\nuse_wind = {}".format(use_friction, use_wind)
param_string += "\nuse_source = {}\nuse_sink = {}".format(use_source, use_sink)
param_string += "\ng = {:g}\nH = {:g}".format(g, H)
计算参数
N_x = 150
N_y = 150
dx = L_x/(N_x - 1)
dy = L_y/(N_y - 1)
dt = 0.1*min(dx, dy)/np.sqrt(g*H)
time_step = 1
max_time_step = 5000
x = np.linspace(-L_x/2, L_x/2, N_x)
y = np.linspace(-L_y/2, L_y/2, N_y)
X, Y = np.meshgrid(x, y)
X = np.transpose(X)
Y = np.transpose(Y)
param_string += "\ndx = {:.2f} km\ndy = {:.2f} km\ndt = {:.2f} s".format(dx, dy, dt)
定义摩擦数组值
if (use_friction is True):
kappa_0 = 1/(5*24*3600)
kappa = np.ones((N_x, N_y))*kappa_0
#kappa[0, :] = kappa_0
#kappa[-1, :] = kappa_0
#kappa[:, 0] = kappa_0
#kappa[:, -1] = kappa_0
#kappa[:int(N_x/15), :] = 0
#kappa[int(14*N_x/15)+1:, :] = 0
#kappa[:, :int(N_y/15)] = 0
#kappa[:, int(14*N_y/15)+1:] = 0
#kappa[int(N_x/15):int(2*N_x/15), int(N_y/15):int(14*N_y/15)+1] = 0
#kappa[int(N_x/15):int(14*N_x/15)+1, int(N_y/15):int(2*N_y/15)] = 0
#kappa[int(13*N_x/15)+1:int(14*N_x/15)+1, int(N_y/15):int(14*N_y/15)+1] = 0
#kappa[int(N_x/15):int(14*N_x/15)+1, int(13*N_y/15)+1:int(14*N_y/15)+1] = 0
param_string += "\nkappa = {:g}\nkappa/beta = {:g} km".format(kappa_0, kappa_0/(beta*1000))
定义风阻
if (use_wind is True):
tau_x = -tau_0*np.cos(np.pi*y/L_y)*0
tau_y = np.zeros((1, len(x)))
param_string += "\ntau_0 = {:g}\nrho_0 = {:g} km".format(tau_0, rho_0)
定义科里奥利数组
if (use_coriolis is True):
if (use_beta is True):
f = f_0 + beta*y
L_R = np.sqrt(g*H)/f_0
c_R = beta*g*H/f_0**2
else:
f = f_0*np.ones(len(y))
alpha = dt*f
beta_c = alpha**2/4
param_string += "\nf_0 = {:g}".format(f_0)
param_string += "\nMax alpha = {:g}\n".format(alpha.max())
param_string += "\nRossby radius: {:.1f} km".format(L_R/1000)
param_string += "\nRossby number: {:g}".format(np.sqrt(g*H)/(f_0*L_x))
param_string += "\nLong Rossby wave speed: {:.3f} m/s".format(c_R)
param_string += "\nLong Rossby transit time: {:.2f} days".format(L_x/(c_R*24*3600))
param_string += "\n================================================================\n"
其他处理
if (use_source):
sigma = np.zeros((N_x, N_y))
sigma = 0.0001*np.exp(-((X-L_x/2)**2/(2*(1E+5)**2) + (Y-L_y/2)**2/(2*(1E+5)**2)))
if (use_sink is True):
w = np.ones((N_x, N_y))*sigma.sum()/(N_x*N_y)
with open("param_output.txt", "w") as output_file:
output_file.write(param_string)
print(param_string)
u_n = np.zeros((N_x, N_y))
u_np1 = np.zeros((N_x, N_y))
v_n = np.zeros((N_x, N_y))
v_np1 = np.zeros((N_x, N_y))
eta_n = np.zeros((N_x, N_y))
eta_np1 = np.zeros((N_x, N_y))
h_e = np.zeros((N_x, N_y))
h_w = np.zeros((N_x, N_y))
h_n = np.zeros((N_x, N_y))
h_s = np.zeros((N_x, N_y))
uhwe = np.zeros((N_x, N_y))
vhns = np.zeros((N_x, N_y))
u_n[:, :] = 0.0
v_n[:, :] = 0.0
u_n[-1, :] = 0.0
v_n[:, -1] = 0.0
eta_n = np.exp(-((X-L_x/2.7)**2/(2*(0.05E+6)**2) + (Y-L_y/4)**2/(2*(0.05E+6)**2)))
eta_list = list(); u_list = list(); v_list = list()
hm_sample = list(); ts_sample = list(); t_sample = list()
hm_sample.append(eta_n[:, int(N_y/2)])
ts_sample.append(eta_n[int(N_x/2), int(N_y/2)])
t_sample.append(0.0)
anim_interval = 20
sample_interval = 1000
t_0 = time.perf_counter()
while (time_step < max_time_step):
u_np1[:-1, :] = u_n[:-1, :] - g*dt/dx*(eta_n[1:, :] - eta_n[:-1, :])
v_np1[:, :-1] = v_n[:, :-1] - g*dt/dy*(eta_n[:, 1:] - eta_n[:, :-1])
if (use_friction is True):
u_np1[:-1, :] -= dt*kappa[:-1, :]*u_n[:-1, :]
v_np1[:-1, :] -= dt*kappa[:-1, :]*v_n[:-1, :]
if (use_wind is True):
u_np1[:-1, :] += dt*tau_x[:]/(rho_0*H)
v_np1[:-1, :] += dt*tau_y[:]/(rho_0*H)
if (use_coriolis is True):
u_np1[:, :] = (u_np1[:, :] - beta_c*u_n[:, :] + alpha*v_n[:, :])/(1 + beta_c)
v_np1[:, :] = (v_np1[:, :] - beta_c*v_n[:, :] - alpha*u_n[:, :])/(1 + beta_c)
v_np1[:, -1] = 0.0
u_np1[-1, :] = 0.0
h_e[:-1, :] = np.where(u_np1[:-1, :] > 0, eta_n[:-1, :] + H, eta_n[1:, :] + H)
h_e[-1, :] = eta_n[-1, :] + H
h_w[0, :] = eta_n[0, :] + H
h_w[1:, :] = np.where(u_np1[:-1, :] > 0, eta_n[:-1, :] + H, eta_n[1:, :] + H)
h_n[:, :-1] = np.where(v_np1[:, :-1] > 0, eta_n[:, :-1] + H, eta_n[:, 1:] + H)
h_n[:, -1] = eta_n[:, -1] + H
h_s[:, 0] = eta_n[:, 0] + H
h_s[:, 1:] = np.where(v_np1[:, :-1] > 0, eta_n[:, :-1] + H, eta_n[:, 1:] + H)
uhwe[0, :] = u_np1[0, :]*h_e[0, :]
uhwe[1:, :] = u_np1[1:, :]*h_e[1:, :] - u_np1[:-1, :]*h_w[1:, :]
vhns[:, 0] = v_np1[:, 0]*h_n[:, 0]
vhns[:, 1:] = v_np1[:, 1:]*h_n[:, 1:] - v_np1[:, :-1]*h_s[:, 1:]
eta_np1[:, :] = eta_n[:, :] - dt*(uhwe[:, :]/dx + vhns[:, :]/dy)
if (use_source is True):
eta_np1[:, :] += dt*sigma
if (use_sink is True):
eta_np1[:, :] -= dt*w
u_n = np.copy(u_np1)
v_n = np.copy(v_np1)
eta_n = np.copy(eta_np1)
time_step += 1
if (time_step % sample_interval == 0):
hm_sample.append(eta_n[:, int(N_y/2)])
ts_sample.append(eta_n[int(N_x/2), int(N_y/2)])
t_sample.append(time_step*dt)
if (time_step % anim_interval == 0):
print("Time: \t{:.2f} hours".format(time_step*dt/3600))
print("Step: \t{} / {}".format(time_step, max_time_step))
print("Mass: \t{}\n".format(np.sum(eta_n)))
u_list.append(u_n)
v_list.append(v_n)
eta_list.append(eta_n)
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
