线性水库模型的解析解与数值解对比

1. 解析解

1.1 模型：

1.2 解析解

2. Code

import numpy as  np
import matplotlib.pyplot as plt
#全局字体样式
plt.rcParams['font.sans-serif'] = ['Times new roman']
#全局字体大小
plt.rcParams['font.size'] = 15
#负号正常显示
plt.rcParams['axes.unicode_minus'] = False
fig1=plt.figure(1)
h0=0
a=0.131
S=0.25
t1=np.arange(0,30,1)
t=np.array([0,5,7,10,12,30])
W=np.array([0,0,0.01,0,0.02,0])
plt.step(t,-W,lw=3,color='C9')
plt.xlabel(r'$t(day)$',fontsize=18)
plt.ylabel(r'$I(m/d)$',fontsize=18)
ax=plt.gca()
axloc=['bottom','top','left','right']
for i in axloc:
    ax.spines[i].set_linewidth(2)
ht=[]
for ti in t1:
    if ti<=t[0]:
        h=np.exp(-(a/S)*ti)*(h0)
        ht.append(h)
    else:
        inner_sum = 0
        tii=np.array([ti])
        i=np.sum(tii.reshape(-1,1)>t,axis=1)[0]
        for j in range(1,i):
            inner_sum+=W[j]/a*(np.exp((a/S)*t[j])-np.exp((a/S)*t[j-1]))
        h=np.exp(-(a/S)*ti)*(W[i]/a*(np.exp((a/S)*ti)-np.exp((a/S)*t[i-1]))+inner_sum+h0)
        ht.append(h)
        
        
fig2=plt.figure(2)
plt.plot(t1,ht,lw=3,color='C9')
df1=np.loadtxt(r'C:\Users\Administrator\Desktop\Gelhar.txt')
plt.scatter(df1[:,0],df1[:,1],color='red',facecolor='none',linewidths=2,s=80)
plt.legend(["Analytical Solution","Numerical Solution"])
plt.xlabel(r'$t(day)$',fontsize=18)
plt.ylabel(r'$z(m)$',fontsize=18)
ax=plt.gca()
axloc=['bottom','top','left','right']
for i in axloc:
    ax.spines[i].set_linewidth(2)
plt.show()

2.1 result

附上comsol数值模型，有需要可以下载：
Gelhar模型——comsol数值模型