🥭Python铅蓄放电热力学无量纲模拟分析

1. 非等压电池放电物理模型模拟高放电对流现象。 2. 配置电池和简化化学性质分析，分析液体热力学化学过程及数学计算方式。 3. 使用无量纲化系统确定系统中的主导效应，构建模型数值解析。 4. 推导出模型系统的三个解析近似解，并在接近小极限的情况下进行渐近分析。

Python铅蓄放电热力学无量纲模拟分析 | 亚图跨际亚图跨际 on Notion

🏈指点迷津 | Brief

🍇Python双曲偏微分方程

浅水方程是一组双曲偏微分方程，用于描述流体中压力面以下的流动（有时但不一定是自由表面）。浅水方程是从质量守恒定律和线性动量守恒定律（纳维-斯托克斯方程）推导出来的，即使浅水假设不成立（例如跨越水跃），该方程仍然成立。在水平床的情况下，科里奥利力、摩擦力和粘性力可以忽略不计，浅水方程为：

$$\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}$$

使用乘积法则展开上面的导数，得到浅水方程的非保守形式。由于速度不受基本守恒方程的影响，非保守形式在冲击或水跃中不成立。还包括科里奥利力、摩擦力和粘性力的适当项，以获得（对于恒定流体密度）：

$$\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}$$

通常情况下，u 和 v 的二次项（表示整体平流的影响）与其他项相比较小。这称为地转平衡，相当于说罗斯贝数很小。假设波高与平均高度相比非常小（h ≪ H），我们有（没有侧向粘性力）：

$$\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}$$

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)
 ​

PyPy