3 模型求解
模型的求解也分为两部分。运动波方程是一个非线性的对流型方程,求解采用了简单的一阶显式迎风格式。经过比较,一阶显式迎风格式对此问题的计算结果在幅值、相位、守恒性几个方面的综合效果较其他一些格式为优。我们曾使用过preissmann四点偏心格式,它在计算此类坡面薄层水流运动波方程时耗散过大,表现为达到平衡产流的时间偏快。我们也曾试过中心差分格式和高阶迎风格式,结果发现标准的中心差分格式仍然和在普通的线性对流型方程中的表现一样是绝对不稳定的。一些变形的中心差分格式能得到较好的结果,但在守恒性和稳定性上不如一阶迎风格式。高阶迎风格式的结果也与此类似,其中二阶迎风格式实际计算时稳定的currant数比一阶迎风小一半。除此之外,一阶迎风格式的计算是最为简单的,边界条件也易于处理。
将式(1)第一式写为
|
|
其中f代表通量,f=f(h)=q=1/nh5/3s01/2,sr为源项,sr=pcosθ-i。
计算格式形式按有限体积法写为
|
hin+1-hni/δt+fi+1/2-fi-1/2/δx=srni |
一阶迎风格式的通量写为
|
fi+1/2(h)=1/hhi5/3s01/2 (10) |
[tn,tn+1]时段入渗方程的求解首先根据两个因子判断tn+1时刻有无积水,若无积水可直接由r(tn+1)=r(tn)得到sr=0。若有积水,按照方程(5)给出的关系直接用newton法求解代数超越方程。若是均匀降雨、坡度不变且土壤的物理参数也不变,对所有的网格可以采用同样的积水时间,可以直接根据(5)式求解。若由于各种参数的变化造成各网格的产流时间不同,则还须考虑产流网格向相邻未产流网格汇入造成该网格的既定入渗曲线的改变。这时cu,cp的公式中还须加入相应的汇入量。
4 模型验证
morgali and linsley(1965)曾在长为72 feet(合22m),坡度s0=0.04的坡面上进行了无入渗降雨产流实验,降雨强度为3.66 inches/hour(1.55mm/min)。实验中分别采用了光滑和粗糙两种表面条件。运用本文模型对无渗透坡面上的降雨产流过程进行模拟计算。几种条件下单宽流量随时间变化过程的数值模拟结果与实验结果的对比如图1,图2,图3所示。所取阻力系数值均按原作者进行数值计算时所取值。结果表明,当阻力系数取值较准确时,数值解与实测数据符合较好。需要说明的是,实验中粗糙表面所设粗糙颗粒远大于在一般坡面上的尺度,因此其阻力系数也较一般无植被覆盖的平坦坡面大得多。退水过程中由于已无雨滴打击的影响,其阻力系数有所减小。
包括渗透过程的降雨产流过程计算我们采用了lima(1992)的实验数据来验证。该实验在长1m,宽0.5m,坡度s0=0.1的土质坡面上进行。降雨强度为0.03741mm/s。本文计算中涉及到的实测土壤参数和根据实测土壤参数与降雨总量率定的参数为:k=1.67×10-6m/s,θs=0.506,θi=0.0107,s=0.02m。计算结果和实验结果见图4。整个产流退水过程均符合较好。与原作者用土壤水分运动微分方程求得的结果也很接近。
|
图1 渗透光滑坡面涨水过程 |
图2 无渗透粗糙坡面涨水过程 |
|
图3 无渗透粗糙坡面退水过程 |
图4 有渗透坡面产流过程 |