1 模型结构

1.1 Richards方程

ZZZ方向Richards方程
∂θ∂t=∂∂z[D(θ)∂θ∂z−K(θ)]\frac{\partial \theta}{\partial t} = \frac{\partial}{\partial z}[D(\theta) \frac{\partial \theta}{\partial z} - K(\theta)] ∂t∂θ=∂z∂[D(θ)∂z∂θ−K(θ)]

式中，D(θ)D(\theta)D(θ)为扩散率，K(θ)K(\theta)K(θ)为导水率

1.2 差分格式

∂θ∂z=θi+12j+1−θi−12j+1Δz\frac{\partial \theta}{\partial z} = \frac{\theta_{i+\frac{1}{2}}^{j+1}-\theta_{i-\frac{1}{2}}^{j+1}}{\Delta z} ∂z∂θ=Δzθi+21j+1−θi−21j+1

∂θ∂t=θij+1−θijΔt\frac{\partial \theta}{\partial t} = \frac{\theta_{i}^{j+1}-\theta_{i}^{j}}{\Delta t} ∂t∂θ=Δtθij+1−θij

1.3 方程离散

θij+1−θijΔt=[D(θ)∂θ∂z−K(θ)]i+12j+1−[D(θ)∂θ∂z−K(θ)]i−12j+1Δz\frac{\theta_{i}^{j+1}-\theta_{i}^{j}}{\Delta t} = \frac{[D(\theta) \frac{\partial \theta}{\partial z} - K(\theta)]_{i+\frac{1}{2}}^{j+1} - [D(\theta) \frac{\partial \theta}{\partial z} - K(\theta)]_{i-\frac{1}{2}}^{j+1}}{\Delta z} Δtθij+1−θij=Δz[D(θ)∂z∂θ−K(θ)]i+21j+1−[D(θ)∂z∂θ−K(θ)]i−21j+1

[D(θ)∂θ∂z−K(θ)]i+12j+1=D(θ)i+12j+1(θi+12j+1−θi−12j+1Δz)i+12j+1−K(θ)i+12j+1=D(θ)i+12j+1θi+1j+1−θij+1Δz−K(θ)i+12j+1[D(\theta) \frac{\partial \theta}{\partial z} - K(\theta)]_{i+\frac{1}{2}}^{j+1} = D(\theta)_{i+\frac{1}{2}}^{j+1} (\frac{\theta_{i+\frac{1}{2}}^{j+1}-\theta_{i-\frac{1}{2}}^{j+1}}{\Delta z})_{i+\frac{1}{2}}^{j+1}- K(\theta)_{i+\frac{1}{2}}^{j+1} = D(\theta)_{i+\frac{1}{2}}^{j+1} \frac{\theta_{i+1}^{j+1}-\theta_{i}^{j+1}}{\Delta z}- K(\theta)_{i+\frac{1}{2}}^{j+1} [D(θ)∂z∂θ−K(θ)]i+21j+1=D(θ)i+21j+1(Δzθi+21j+1−θi−21j+1)i+21j+1−K(θ)i+21j+1=D(θ)i+21j+1Δzθi+1j+1−θij+1−K(θ)i+21j+1

[D(θ)∂θ∂z−K(θ)]i−12j+1=D(θ)i−12j+1(θi+12j+1−θi−12j+1Δz)i−12j+1−K(θ)i−12j+1=D(θ)i−12j+1θij+1−θi−1j+1Δz−K(θ)i−12j+1[D(\theta) \frac{\partial \theta}{\partial z} - K(\theta)]_{i-\frac{1}{2}}^{j+1} = D(\theta)_{i-\frac{1}{2}}^{j+1} (\frac{\theta_{i+\frac{1}{2}}^{j+1}-\theta_{i-\frac{1}{2}}^{j+1}}{\Delta z})_{i-\frac{1}{2}}^{j+1}- K(\theta)_{i-\frac{1}{2}}^{j+1} = D(\theta)_{i-\frac{1}{2}}^{j+1} \frac{\theta_{i}^{j+1}-\theta_{i-1}^{j+1}}{\Delta z}- K(\theta)_{i-\frac{1}{2}}^{j+1} [D(θ)∂z∂θ−K(θ)]i−21j+1=D(θ)i−21j+1(Δzθi+21j+1−θi−21j+1)i−21j+1−K(θ)i−21j+1=D(θ)i−21j+1Δzθij+1−θi−1j+1−K(θ)i−21j+1

则
θij+1−θijΔt=1Δz[D(θ)i+12j+1θi+1j+1−θij+1Δz−K(θ)i+12j+1−D(θ)i−12j+1θij+1−θi−1j+1Δz+K(θ)i−12j+1]\frac{\theta_{i}^{j+1}-\theta_{i}^{j}}{\Delta t} = \frac{1}{\Delta z}[D(\theta)_{i+\frac{1}{2}}^{j+1} \frac{\theta_{i+1}^{j+1}-\theta_{i}^{j+1}}{\Delta z}- K(\theta)_{i+\frac{1}{2}}^{j+1}-D(\theta)_{i-\frac{1}{2}}^{j+1} \frac{\theta_{i}^{j+1}-\theta_{i-1}^{j+1}}{\Delta z}+ K(\theta)_{i-\frac{1}{2}}^{j+1}] Δtθij+1−θij=Δz1[D(θ)i+21j+1Δzθi+1j+1−θij+1−K(θ)i+21j+1−D(θ)i−21j+1Δzθij+1−θi−1j+1+K(θ)i−21j+1]

ΔtD(θ)i−12j+1Δz2θi−1j+1+{1+ΔtΔz2[D(θ)i+12j+1−D(θ)i−12j+1]}θij+1−ΔtD(θ)i+12j+1Δz2θi+1j+1=θij−ΔtΔzK(θ)i+12j+1+ΔtΔzK(θ)i−12j+1\frac{\Delta t D(\theta)_{i-\frac{1}{2}}^{j+1}}{{\Delta z}^2}\theta_{i-1}^{j+1}+\{1+\frac{\Delta t}{{\Delta z}^2}[D(\theta)_{i+\frac{1}{2}}^{j+1}-D(\theta)_{i-\frac{1}{2}}^{j+1}]\}\theta_i^{j+1}-\frac{\Delta t D(\theta)_{i+\frac{1}{2}}^{j+1}}{{\Delta z}^2}\theta_{i+1}^{j+1} = \theta_i^j-\frac{\Delta t}{\Delta z}K(\theta)_{i+\frac{1}{2}}^{j+1}+\frac{\Delta t}{\Delta z}K(\theta)_{i-\frac{1}{2}}^{j+1} Δz2ΔtD(θ)i−21j+1θi−1j+1+{1+Δz2Δt[D(θ)i+21j+1−D(θ)i−21j+1]}θij+1−Δz2ΔtD(θ)i+21j+1θi+1j+1=θij−ΔzΔtK(θ)i+21j+1+ΔzΔtK(θ)i−21j+1

即
Aiθi−1j+1+Biθij+1+Ciθi+1j+1=DiA_i\theta_{i-1}^{j+1} + B_i\theta_{i}^{j+1}+C_i\theta_{i+1}^{j+1}=D_i Aiθi−1j+1+Biθij+1+Ciθi+1j+1=Di
式中
Ai=ΔtΔz2D(θ)i−12j+1A_i = \frac{\Delta t }{{\Delta z}^2}D(\theta)_{i-\frac{1}{2}}^{j+1} Ai=Δz2ΔtD(θ)i−21j+1

Bi=1+ΔtΔz2[D(θ)i+12j+1−D(θ)i−12j+1]B_i = 1+\frac{\Delta t}{{\Delta z}^2}[D(\theta)_{i+\frac{1}{2}}^{j+1}-D(\theta)_{i-\frac{1}{2}}^{j+1}] Bi=1+Δz2Δt[D(θ)i+21j+1−D(θ)i−21j+1]

Ci=−ΔtΔz2D(θ)i+12j+1C_i = -\frac{\Delta t }{{\Delta z}^2}D(\theta)_{i+\frac{1}{2}}^{j+1} Ci=−Δz2ΔtD(θ)i+21j+1

Di=θij−ΔtΔz[K(θ)i+12j+1−K(θ)i−12j+1]D_i = \theta_i^j-\frac{\Delta t}{\Delta z}[K(\theta)_{i+\frac{1}{2}}^{j+1}-K(\theta)_{i-\frac{1}{2}}^{j+1}] Di=θij−ΔzΔt[K(θ)i+21j+1−K(θ)i−21j+1]

采用半隐式解法，用jjj时刻的DDD和KKK代替j+1j+1j+1时刻的值，采用上游迁移法或
Ki±12=Ki+Ki±12K_{i \pm \frac{1}{2}}=\frac{K_i+K_{i\pm1}}{2} Ki±21=2Ki+Ki±1

Di±12=Di+Di±12D_{i \pm \frac{1}{2}}=\frac{D_i+D_{i\pm1}}{2} Di±21=2Di+Di±1

2 模型条件

2.1 边界条件

上下边界可以是土壤含水率固定值或时间序列θ(t)\theta(t)θ(t)，在三对角矩阵的第一行和最后一行添加。

2.2 初始条件

θx0=θ0\theta_x^0 = \theta_0 θx0=θ0

2.3 K(θ)K(\theta)K(θ)

非饱和土壤导水率KKK与土壤水吸力sss、或与土壤含水率θ\thetaθ的关系，常根据实测结果拟合为经验公式。常用的经验公式形式有
K(θ)=Kt(θθs)m或K(θ)=Kt(θ−θ0θs−θ0)mK(\theta) = K_t(\frac{\theta}{\theta_s})^m 或 K(\theta)=K_t(\frac{\theta - \theta_0}{\theta_s - \theta_0})^m K(θ)=Kt(θsθ)m或K(θ)=Kt(θs−θ0θ−θ0)m

K(θ)=Kte−β(θs−θ)K(\theta) = K_t e^{-\beta(\theta_s-\theta)} K(θ)=Kte−β(θs−θ)

式中，KtK_tKt为饱和导水率，θs\theta_sθs为饱和含水率，θ0\theta_0θ0为某一特征含水率，mmm为经验常数。

2.4 D(θ)D(\theta)D(θ)

定义非饱和土壤水的扩散率D(θ)D(\theta)D(θ)为导水率K(θ)K(\theta)K(θ)和比水容量C(θ)C(\theta)C(θ)的比值，即
D(θ)=K(θ)C(θ)=K(θ)/dθdψmD(\theta)=\frac{K(\theta)}{C(\theta)}=K(\theta)/\frac{d\theta}{d\psi_m} D(θ)=C(θ)K(θ)=K(θ)/dψmdθ
显然，非饱和土壤水的扩散率DDD同样是土壤含水率θ\thetaθ的或基质势ψm\psi_mψm的函数。其函数关系必须通过试验测定，常用的经验公式形式有
D(θ)=aebθD(\theta) = ae^{b\theta} D(θ)=aebθ

D(θ)=D0(θθs)mD(\theta) = D_0(\frac{\theta}{\theta_s})^m D(θ)=D0(θsθ)m

D(θ)=D0e−β(θ0−θ)D(\theta) = D_0e^{-\beta(\theta_0-\theta)} D(θ)=D0e−β(θ0−θ)

式中，θ0\theta_0θ0为某一特征含水率，aaa，bbb，D0D_0D0、mmm、β\betaβ为经验常数，取决于土壤质地和结构。

3 模型求解

3.1 三对角矩阵

三对角矩阵Ax=BAx=BAx=B形式如下：
(1000………………0A2B2C2A3B3C3A4B4C4An−2Bn−2Cn−2An−1Bn−1Cn−10………………0001)⏟A(θ1θ2θ3θ4θ5θn−2θn−1θn)⏟x=(θ1(t)D2D3D4D5Dn−2Dn−1θn(t))⏟B\underbrace{\begin{pmatrix} 1&0&0&0&\dots&\dots&\dots&\dots&\dots&\dots&0\\[4pt] A_{2}&B_{2}&C_{2}\\[4pt] & A_{3}&B_{3}&C_{3}\\[4pt] & & A_{4}&B_{4}&C_{4} \\[4pt] \\[4pt] & & &&&&&&\\[4pt] & & &&&&&&\\[4pt] & & &&&&&&\\[4pt] & & &&&&&&\\[4pt] & &&&&&& A_{n-2}&B_{n-2}&C_{n-2} \\[4pt] & & &&&&&&A_{n-1}&B_{n-1}&C_{n-1} \\[4pt] 0&\dots &\dots &\dots&\dots&\dots&\dots&0 & 0 &0 &1 \\[4pt] \end{pmatrix}}_{A} \underbrace{\begin{pmatrix} \theta_{1} \\[4pt] \theta_{2} \\[4pt] \theta_{3} \\[4pt] \theta_{4} \\[4pt] \theta_{5} \\[4pt] \\[4pt] \\[4pt] \\[4pt] \\[4pt] \theta_{n-2} \\[4pt] \theta_{n-1} \\[4pt] \theta_{n} \\[4pt] \end{pmatrix}}_{x}=\underbrace{\begin{pmatrix} \theta_{1}(t) \\[4pt] D_{2} \\[4pt] D_{3} \\[4pt] D_{4} \\[4pt] D_{5} \\[4pt] \\[4pt] \\[4pt] \\[4pt] \\[4pt] D_{n-2} \\[4pt] D_{n-1} \\[4pt] \theta_{n}(t) \\[4pt] \end{pmatrix}}_{B} A⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛1A200B2A3…0C2B3A4…0C3B4……C4………………An−20…Bn−2An−10…Cn−2Bn−100Cn−11⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞x⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛θ1θ2θ3θ4θ5θn−2θn−1θn⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=B⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛θ1(t)D2D3D4D5Dn−2Dn−1θn(t)⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞

3.2 各参数取值

参数	含义	取值
Z	土柱深度	1m
T	模拟时间	1h
dz	空间步长	1cm
dt	时间步长	0.0001h
θ1\theta_1θ1	上边界土壤含水量	0.30
θs\theta_sθs	饱和土壤含水量	0.45
θ0\theta_0θ0	初始土壤含水量	0.17
KtK_tKt	饱和导水率¹	6.72m/h
mkm_kmk	导水率经验常数¹	3
aaa	扩散率经验常数²	0.019
bbb	扩散率经验常数²	12.3804

3 求解结果

3.1 计算结果

3.2 MATLAB模型代码

%========================================================%程序说明
%创建时间：2020.12.22
%最后修改时间：2020.12.28
%程序目标：求Richards方程数值解
%离散方法：半隐式法
%边界条件：土壤含水率固定值或时间序列$\theta(t)$
%上边界：上边界取第一类边界条件，$\theta = \theta_1$
%下边界：底部土壤含水量恒定，等于饱和土壤含水量,$\theta = \theta_s$
%初始条件：
%$\theta_x^0 = \theta_0$
%上游迁移法计算K(i+-1/2)、D(i+-1/2)
%$K(\theta) = K_t(\frac{\theta}{\theta_s})^m$
%$D(\theta) = a*e^{b\theta}$
%方程及结果：https://lujiabo98.github.io/2020/12/22/Course03/%=========================================================%变量定义及初始化
clc,clear;%Z：土柱深度,m，T：时间,s，dz：空间步长,m，dt：时间步长,s，theta_1：上边界土壤含水量，theta_s：饱和土壤含水量，theta_0：初始土壤含水量,theta_r：土壤干旱体积含水量
Z=1; T=1*3600; dz=0.01; dt=0.0001*3600; theta_1 = 0.3; theta_s = 0.45; theta_0 = 0.17; theta_r = 0.132;
%K_t:饱和导水率,m/s; m_k：导水率经验常数; a：扩散率经验常数; b：扩散率经验常数;
K_t = 6.72/3600; m_k = 3; a = 0.019; b = 12.3804;%n：土柱分段数，N：断面数即theta_i变量数，lambda:网格比，t：时间分段数
n=Z/dz; N=n+1; lambda=dt/dz^2; t = round(T/dt);%x：theta_i解向量，B：三对角矩阵方程的常数项，I,J,M:三对角矩阵A的三列，K,D：非饱和土壤导水率和扩散率
x=zeros(1,N); B=zeros(1,N); I=zeros(1,N); J=zeros(1,N-1); M=zeros(1,N-1); K=zeros(1,N); D=zeros(1,N);
%K,D：非饱和土壤导水率和扩散率K(i+-1/2)、D(i+-1/2),K_p=K(i+1/2),K_m=K(i-1/2),D_p=D(i+1/2),D_m=D(i-1/2)
K_p=zeros(1,N); K_m=zeros(1,N); D_p=zeros(1,N); D_m=zeros(1,N);%theta为模拟时段内的土壤含水量过程，X：土壤含水量组合矩阵
theta = zeros(t+1,N); X=zeros(N,t+1);%初始化theta_i解向量x
for i=1:1:Nx(i) = theta_0;
end%上下边界土壤含水量初始化
x(1) = theta_1;
x(N) = theta_s;%将初始时刻的解向量放入组合矩阵中存储
X(:,1)=x';%初始化B,I,J,K向量
B(1) = x(1); B(N) = x(N);
I(1) = 1; I(N) = 1;
J(1) = 0; M(N-1) = 0;%初始化非饱和土壤导水率K(thtea)和扩散率D(theta)
K = K_t * (x / theta_s).^m_k;
D = a * exp(b * x) * 10^(-4) /60;%=========================================================
%程序主体部分
for k=1:1:t  %按时间层循环%上游迁移法计算K(i+-1/2)、D(i+-1/2)K_p = K; K_m(1) = K(1); K_m(2:end) = K(1:end-1);D_p = D; D_m(1) = D(1); D_m(2:end) = D(1:end-1);%B,I,J向量赋值%B(i)为D_i = \theta_i^j-\frac{\Delta t}{\Delta z}[K(\theta)_{i+\frac{1}{2}}^{j+1}-K(\theta)_{i-\frac{1}{2}}^{j+1}]B = x - dz * lambda * (K_p - K_m);B(1) = x(1); B(N) = x(N);%三对角矩阵顶角向量I赋值，I(i)为方程组系数B(i)I = 1 + lambda * (D_p - D_m);I(1) = 1; I(N) = 1;%三对角矩阵上一向量J赋值，J(i)为方程组系数C(i)J = - lambda * D_p(1:end-1);J(1) = 0;%三对角矩阵下一向量M赋值，M(i)为方程组系数A(i)M = lambda * D_m(2:end);M(N-1) = 0;    %组合I,J,M向量得到三对角矩阵AA = diag(I) + diag(J,1) + diag(M,-1);%解五对角矩阵差分方程组，并以列形式存储在组合矩阵X中X(:,k+1) = A\B';%这一次的解作为下一次循环的初始值x=X(:,k+1)';%更新非饱和土壤导水率K(thtea)和扩散率D(theta)K = K_t * (x / theta_s).^m_k;D = a * exp(b * x) * 10^(-4) /60;
end  %=========================================================
%结果输出
%从组合矩阵X中提取出土壤含水量
%X的每一列表示一个时段的计算结果
%theta的每一行表示一个时段的计算结果theta = X'; %绘制出土壤含水量三维图
figure(1);
mesh(theta);
xlabel('断面n/1cm');
ylabel('时间t/0.0001h');
zlabel('土壤含水量');
%=========================================================

4 参考文献

[1]土壤水动力学，1987.3

[2]李映强.赤红壤非饱和土壤水扩散率及其影响因素[J].华南农业大学学报,1998(02):3-5.

【课程】03 Richards方程数值解相关推荐

泊松方程与拉普拉斯方程数值解
泊松方程与拉普拉斯方程数值解 posted on 2019-01-24 22:01 luoganttcc 阅读(...) 评论(...) 编辑收藏
[CUPT]国一博主, 教你求解95%以上的方程(数值解)
大家好呀~ 博主在2021年cupt中取得了国家级一等奖的成绩, 那么也是在这里和大家分享自己的求解各种方程(数值解)的技巧, 具体讲解在B站:[cupt]国家级一等奖:教你求解95%以上的方程_哔哩 ...
matlab数值计算pdf_Gnuplot科学绘图（九）——栅格以及方程数值解估算
Gnuplot科学绘图系列内容 Gnuplot科学绘图(一)--从安装到简单函数绘图(文末有彩蛋) Gnuplot科学绘图(二)--坐标取值范围及刻度(文末有彩蛋) Gnuplot科学绘图(三)--点 ...
不动点迭代求解方程数值解
求解方程 2x−x3=0 解:可以采用不动点迭代的方式求出数值解. 其不动点是 x=3log2x ,故可以采用下式进行不动点迭代: xn+1=3log2xn 初始值为 x0=2 ,迭代终止条件是 |x ...
Java课程03总结
一:继承条件下的构造方法调用 package test;class Grandparent {public Grandparent(){System.out.println("GrandPa ...
Java项目课程03：涉及知识点
文章目录一.Java基本语法 (一)数据类型 (二)变量与常量 (三)运算符与表达式二.Java流程控制 (一)选择结构 (二)循环结构三.Java面向对象编程 (一)封装 (二)继承 (三)多 ...
TensorFlow(keras)入门课程--03 卷积介绍
目录 1 简介 2 使用卷积 3 开始编码 4 创建卷积 5 查看结果 6 了解池化 7 编写池化代码 1 简介在本节中,你将了解卷积以及它们在计算机视觉场景中如此强大的原因. 在上一节中,我们了解 ...
运筹学状态转移方程例子_强化学习第4期：H-J-B方程
在上一篇文章中,我们介绍了一种最简单的MDP--s与a都是有限的MDP的求解方法.其中,我们用到了动态规划的思想,并且推出了"策略迭代"."值迭代"这样的方法. ...
2022年Java项目课程目录
一.2022Java项目课程目录 Java项目课程01:课程概述 Java项目课程02:系统概述 Java项目课程03:涉及知识点 Java项目课程04:需求分析 Java项目课程05:系统设计 Ja ...

【课程】03 Richards方程数值解