2.3 极坐标下的分数阶气、液迁移方程的数值解法

在实际工程中，受外界因素的影响，土壤中的水分迁移是不稳定的。为减少外界的影响，本节采用密闭的PVC圆柱管进行室内模拟，由于边界条件被限定，因此数值分析中不考虑土体体积的变化。水分迁移方程采用极坐标形式表示，以图2.3所示的土体单元为例。

假设 D _{w r} =D _{w η} =D _{w x} 、 D _{v r} =D _{v η} =D _{v x} ，柱坐标中任一点由3个坐标给出，相对于3个方向的通量为

在Δ t 时段内，在 r 方向的流入量为 q _r rΔηΔxΔt ，流出量为 Δr ） ΔηΔxΔt ，忽略高阶无穷小量，则流出量和流入量之差为

同理，在 η 方向的流出量和流入量之差为

在 x 方向的流出量和流入量之差为

上述3个方向的水量差总计为

单元体体积应为 ，略去高阶无穷小量后为 rΔrΔηΔx ，根据式（2.10），在Δ t 时间内单元体内水分变化量为

式（2.33）和式（2.34）相等，即

将式（2.29）代入式（2.35），即

式（2.36）为极坐标系下基于时间分数阶的土体气、液迁移方程。

本书所进行的水分迁移试验不考虑重力的影响，并且假定土体各向同性，水分只沿水平方向迁移，因此可以用一维情况来模拟计算。方程定解条件包括初始条件和边界条件。

初始条件表示为

其中， ω _i 、 ω _j 均为初始时刻的含水量； 为土柱的湿段部分； 为土柱的干段部分。

边界条件表示为

式中， Y 为迁移总时长。

在设计气态水迁移的试验时，为了避免液态水的参与，土柱的湿段和干段之间间隔了10mm的空气段，这种情况保证了水分只能以气态水形式迁移，但在湿段或干段内部却是气、液混合的。如果按照全程只考虑气态水迁移的含水量方程进行计算，那么势必与试验结果存在较大误差。在此分段进行含水量变化的数值计算，以保证每个阶段参数取值的准确性。

对迁移方程的数值计算，可以分为以下3个阶段。

（1）空气段。空气段整个截面只有水蒸气通过，时间分数阶土壤水分运动质量守恒方程为

作为单根均匀直径的孔隙，扩散系数取

式中， T 为绝对温度（K）。在恒温条件下，土壤孔隙中的液态水和水蒸气有足够的时间凝结或蒸发，达到或接近平衡状态，故水蒸气始终处于饱和状态 ^[27] 。因此，在恒温条件下，水平土柱中土壤孔隙中的空气流动可忽略不计，可以认为土壤中的气体是停滞的。水蒸气密度由下式确定：

式中， W _v 为水的分子量（0.018kg/mol）； R 为通用气体常数（8.314J/mol·K）； P _v 为孔隙中的实际水蒸气压力。因此

Edlefsen等 ^[98] 提出以下关系式来描述：

式中， P _v，sat 为纯水饱和水蒸气压力 ^[90] ，取值见式（2.6）； h _r 为相对湿度。

式中， u _a - u _w 为基质吸力。

那么

其中

则

当温度在 x 方向的变化量为0时，式（2.48）可简化为

（2）湿段。这部分土体处于脱湿过程中。对于液态水部分，扩散系数计算采用式（2.16）或式（2.17），在此将式（2.16）代入迁移方程式（2.15），液态水部分可得

水蒸气在土壤中的扩散系数为

那么水蒸气的迁移可以写为

因此

其中

式中， A 为水蒸气通过的土体截面积，是含水量的函数。在此对 A 的求导表示为

土壤基质吸力在 x 方向上的变化主要是由含水量 ω 的变化引起的，如果已知土壤的土水特征曲线（Soil Water Characteristic Curve，SWCC），那么

式中， k 为SWCC切线对应的斜率。在湿段部分， k 应根据土水特征曲线脱湿曲线计算得出。将式（2.50）和式（2.53）联合，即为湿段液态水和气态水的总迁移量。

（3）干段。这部分土体主要处于吸湿阶段，为混合水迁移。迁移公式与湿段相同，但此阶段 k 应根据土水特征曲线吸湿曲线计算得出。

上述3个阶段水分迁移方程可以总结为

在混合水迁移试验中，在基质势的推动下，水分从湿段向干段迁移，水分迁移形式包括液态水和气态水。在进行数值计算时，与气态水迁移方程类似，但是其中没有空气段，参数 M 、 N 按照式（2.58）和式（2.59）的土体段取值。需要注意的是，混合水迁移的土体有湿段和干段，两段处于脱湿和吸湿不同的过程，所以其参数 k 的取值有所不同。 Egfl4kqenynfkWlfglci3Dqk1tFwYKK2ALHPBw7HGT5Cz+L2Df2GkeH0DWaTAtGe