什么是储水模型？非饱和渗流方程怎么求解？

非饱和渗流微分方程的基本形式

在引入非饱和方程前，我们首先来看一下非常简单的饱和渗流方程（达西定律）：

（水头与压力的关系为，Z为高程）

不必多说，这个简单的微分方程就描述了：在重力存在的情况下，水会从水头高处流向水头低处这一自然现象，渗透系数K与多孔材料的性质相关。
此方程是线性方程，因为K的值与H无关，是一个常量。

但是在多孔介质中，由于水表面张力的存在，自然会产生水压力为负压的非饱和区域。负压区域的水并未充满整个介质空间，饱和度Se≤1，渗透系数会随饱和度的减小而减小，于是非饱和渗流方程表示非线性的以下形式：

方程因此变为非线性，因为K的值与H(p)有关，而非常量。

渗流场的求解

单位时间单位体积内流体质量守恒（j代表流体速度）

对于求解稳态流场，我们可以用有限元法，将边界条件和上述控制方程代入即可迭代求解（非饱和渗流方程因非线性存在需要多次迭代）。

但是对于求解瞬态问题，多了一个时间变量t，因此只有上面这个速度方程是不够的（多一个变量多一个方程），因为流体是存在一定的压缩性的，某一时刻某单位体积内质量的变化量，等于流进流出质量的差值，于是对于瞬态问题还有以下质量守恒方程（不考虑物理化学反应生成）。

其中θ表示水体积分数，ρ为水密度，u是速度场。方程左边两项分别表示：单位时间单位体积内水质量（密度与体积分数乘积）变化，等于流进流出质量之差。

然后我们对这个微分方程稍微做进行一定的处理，使其出现水密度ρ和压力p的关系项以及水体积分数θ与压力p的关系项，因为这两项的值相对更容易获得。

由此接下来我们将转换到对和的研究与获取。

水与多孔介质的可压缩性

在研究水的体积可压缩特性（单位体积内流体的体积对压力p的变化率）时，有以下定义，称之为压缩系数：

，单位Pa^-1

压缩系数的倒数，可以认为是水的“体积弹性模量”，单位为Pa。

同理对于多孔介质来说，其压缩系数等于固体颗粒的压缩系数k_s乘上固体部分占总体积的比率1-θ_s。所以对于同种固体材料，多孔介质的孔隙比越大则压缩系数越小，体积弹性模量越大，压缩性越大。

因此含水多孔介质的可压缩性由固体的可压缩性和液体的可压缩性组成，其中液体可视为近似不可压缩，多孔介质的储水来源于孔隙水和孔隙的压缩。

滞留模型

在研究渗透系数，含水量和孔隙水压力的关系时，Van Genuchten给出了一组描述K(p)-p以及θ-p关系的经验方程：

k_r称为折减系数，α、l、n、m为材料参数，渗透系数K可以用饱和渗透系数K_s乘上相对渗透系数k_r来表示，同时相对渗透系数k_r随饱和度减小而减小，饱和度又随孔隙水压力的减小而减小，因此渗透系数可以表示为水压力的而函数，在饱和时的值远大于非饱和。在孔隙水压由负转正的过程中，多孔介质中液体体积分数θ的分布从一个最低值θ_r增大到孔隙率大小θ_s，同时有效饱和度Se也从0增大到1。单位容水度C_m描述了θ随p变化的过程。下图描述了这种变化。（注意此方程组只是经验方程，非严格数学方法推导而来）

Van Genuchten 给出的水分滞留曲线

Richard方程的求解

接下来我们就可以将进一步分解：

于是质量守恒方程分解为

（其中）

方程中（此处做了θ约等于Seθ的近似），称之为储水系数，与多孔介质的物理力学性质有关；C_m为单位容水度，表示了单位压力引起的体积含水量变化值，与多孔介质的细观性质有关。(若饱和度Se=1，方程写为，即是饱和渗流方程）。

Richard方程中由于ρ，S，Se，C_m都随压力p发生变化，所以此微分方程存在强非线性。不过由于我们进行了一系列数学分解，使得Richard方程中给定初始条件后仅这一项为未知项其他均为已知项，此时便可使用对时间t的向后差分法，一步步求解整个过程中的速度场与孔压场。（例如给定初始时刻孔压场，首先根据各个非饱和参数与压力的关系求得K、Se、θ、C_m、S，再根据速度方程求出流场，最后代入质量守恒方程求得，给定一个很小的时间步长例如s，求得，则，重复上述步骤，直至求得给定时刻的孔压以及速度场。）