【课程】03 Richards方程数值解

卢家波

已于 2024-04-27 17:41:00 修改

阅读量6.4k

点赞数 5

分类专栏： # 课程文章标签： Richards 数值解 MATLAB 土壤水

于 2021-01-07 22:22:23 首次发布

本文链接：https://blog.csdn.net/weixin_43012724/article/details/112341344

版权

课程专栏收录该内容

9 篇文章 11 订阅

订阅专栏

1 模型结构

1.1 Richards方程

$Z$ 方向Richards方程
$\frac{\partial \theta}{\partial t} = \frac{\partial}{\partial z}[D(\theta) \frac{\partial \theta}{\partial z} - K(\theta)]$

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

1.2 差分格式

$\frac{\partial \theta}{\partial z} = \frac{\theta_{i+\frac{1}{2}}^{j+1}-\theta_{i-\frac{1}{2}}^{j+1}}{\Delta z}$

$\frac{\partial \theta}{\partial t} = \frac{\theta_{i}^{j+1}-\theta_{i}^{j}}{\Delta t}$

1.3 方程离散

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

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

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

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

即
$A_i\theta_{i-1}^{j+1} + B_i\theta_{i}^{j+1}+C_i\theta_{i+1}^{j+1}=D_i$
式中
$A_i = \frac{\Delta t }{{\Delta z}^2}D(\theta)_{i-\frac{1}{2}}^{j+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}]$

$C_i = -\frac{\Delta t }{{\Delta z}^2}D(\theta)_{i+\frac{1}{2}}^{j+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}]$

采用半隐式解法，用 $j$ 时刻的 $D$ 和 $K$ 代替 $j + 1$ 时刻的值，采用上游迁移法或
$K_{i \pm \frac{1}{2}}=\frac{K_i+K_{i\pm1}}{2}$

$D_{i \pm \frac{1}{2}}=\frac{D_i+D_{i\pm1}}{2}$

2 模型条件

2.1 边界条件

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

2.2 初始条件

$\theta_x^0 = \theta_0$

2.3 $K(\theta)$

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

$K(\theta) = K_t e^{-\beta(\theta_s-\theta)}$

式中， $K_t$ 为饱和导水率， $\theta_s$ 为饱和含水率， $\theta_0$ 为某一特征含水率， $m$ 为经验常数。

2.4 $D(\theta)$

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

$D(\theta) = D_0(\frac{\theta}{\theta_s})^m$

$D(\theta) = D_0e^{-\beta(\theta_0-\theta)}$

式中， $\theta_0$ 为某一特征含水率， $a$ ， $b$ ， $D_0$ 、 $m$ 、 $\beta$ 为经验常数，取决于土壤质地和结构。

3 模型求解

3.1 三对角矩阵

三对角矩阵 $A x = 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}$

3.2 各参数取值

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

3 求解结果

3.1 计算结果

result

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:N
    x(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向量得到三对角矩阵A
    A = 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.

作者简介

很高兴认识您！
我叫卢家波，河海大学水文学及水资源博士研究生，研究兴趣为高效洪水淹没预测、洪水灾害预警、机器学习、替代模型和降阶模型。
变化环境下，极端洪水事件多发，我希望能通过研究为水灾害防御做出贡献，为人民服务。
欢迎交流讨论和研究合作，vx Jiabo_Lu。
主页 https://lujiabo98.github.io
简历 https://lujiabo98.github.io/file/CV_JiaboLu_zh.pdf
博客 https://blog.csdn.net/weixin_43012724?type=blog
来信请说明博客标题及链接，谢谢。

卢家波

关注

5
点赞
踩
34

收藏

觉得还不错? 一键收藏
打赏
3
评论
【课程】03 Richards方程数值解

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(θ)为.
复制链接

扫一扫