一种基于元胞自动机的非饱和地层入渗分析方法技术

本发明专利技术提供了一种基于元胞自动机的非饱和地层入渗分析方法，结合达西定律和元胞自动机模型，构建平稳型元胞空间进行迭代模拟计算，输出水力特征曲线等表征非饱和土中孔隙水的实时状态。本发明专利技术的技术方案构建的元胞自动机模型，可在自初始状态开始后，每次完成一个时间步长的计算，每一个元胞空间都处于固定状态，不再随时间变化，计算方法简单、有序，并计算过程有迹可循；通过平稳型元胞空间进行迭代模拟计算，输出水力特征曲线等表征非饱和土中孔隙水的实时状态，提高了EFGM离散渗流模型在非饱和土渗流分析领域的工程实用性；采用无网格的方式让模型建立更加简单高效，数值计算过程也能有效避免网格畸变、计算量偏大、精度不足等问题。足等问题。足等问题。

【技术实现步骤摘要】
一种基于元胞自动机的非饱和地层入渗分析方法


[0001]本专利技术涉及岩土工程、水动力学
，尤其涉及一种基于元胞自动机的非饱和地层入渗分析方法。

技术介绍

[0002]工程上引发边坡失稳的根本原因是土体内部某个滑动面上的剪应力达到其抗剪强度，而导致抗剪强度监督的因素有多种，如降雨入渗土体软化、孔隙水压力增大等。地下水运动形成的渗流场影响着一定区域内地下水的赋存规律，污染物质的迁移等诸多方面；因此若对于地下水运动有充分的认识，无论是对于地下水的开采使用，还是基坑降水方案设计等都具有十分重要的工程意义。其中，降雨入渗是一个十分复杂的过程，对此国内外学者对湿润锋的影响开展了理论研究，研究表明：水自上而下的渗透过程让土体的湿润锋呈现高度非线性变化，导致在使用传统有限元方法模拟渗流时，会出现网格畸变、计算量大、精度低甚至不收敛的现象。Phoon K，Tan T，Chong P.等人的研究中发现当网格单元密度过小时，会使数值结果在湿润锋附近产生振荡。而常用的无网格法中，自由单元伽辽金法(Element
‑
free Galerkin method,EFGM)精度高、形式简单，但应用于渗流分析中仍存在节点的密度分布问题，并且大多数局限于饱和渗流分析，对于非饱和土的研究开展较少。因此亟待提出面向非饱和土的精度高、计算方便，能有效避免因湿润锋的变化导致节点密度不够等问题的计算方法。

技术实现思路

[0003]针对以上技术问题，本专利技术公开了一种基于元胞自动机的非饱和地层入渗分析方法，有效避免了网格畸变、计算量偏大、精度不足、计算振荡等问题。
[0004]对此，本专利技术采用的技术方案为：
[0005]一种基于元胞自动机的非饱和地层入渗分析方法，其包括如下步骤：
[0006]步骤S1，确定初始条件，划分土层的元胞，设定元胞和节点
[0007]步骤S2，构造水力特征矩阵，包括位置水头矩阵h
z
、压力水头矩阵h
p
、元胞中心位置水头矩阵(h
z
)
ele
、元胞中心压力水头矩阵(h
p
)
ele
、元胞渗透系数矩阵k、元胞体积含水率系数矩阵、元胞孔隙水体积矩阵；
[0008]步骤S3，设定时间步长dt，步数N；
[0009]步骤S4，从节点1开始至节点n+1，依次计算每个节点位置水头h
z,i
和压力水头h
p,I
；
[0010]步骤S5，初始化每个元胞中心的水力特性：令θ
i
＝θ
r
，计算元胞中心的位置水头(h
z,i
)
ele
和压力水头(h
p,i
)
ele
，并根据式(1)的本构关系计算元胞渗透系数k
i
和元胞体积含水率θ
i
；
[0011][0012]上式中，k表示当前状态的渗透系数，ks表示饱和渗透系数，θ表示当前状态的体积
含水率，θ
s
表示饱和体积含水率，θ
r
表示残余体积含水率，
[0013]步骤S6，根据式(2)计算第一个元胞的水力梯度Δh0；
[0014]Δh0＝(h
p,1
+h
z,1
)
‑
[(h
p,1
)
ele
+(h
z,1
)
ele
]＝(0+h
z,1
)
‑
[(h
p,1
)
ele
+(h
z,1
)
ele
]ꢀꢀꢀ
(2)
[0015]步骤S7，根据达西定律，计算从地表入渗到第一个元胞的流量；
[0016]步骤S8，更新第一个元胞的状态；
[0017]步骤S9，依次计算第二个元胞的水力梯度、从第一个元胞中心入渗到第二个元胞中心的渗流量，更新第二个元胞的状态；重复此步骤，更新第三至第N个元胞；
[0018]步骤S10，依次根据时间步长，重复迭代计算N次后各元胞中心压力水头，以每个元胞的水力特征状态得到水力特征曲线。
[0019]其中，进行土层的元胞划分可以从上到下进行一维划分，也可以将土层划分二维和三维空间下渗作用的元胞空间。
[0020]此技术方案采用了元胞自动机理论，构建了平稳性动力学行为的元胞自动机模型，对非饱和土在降雨入渗下地下水水头分布和含水量的模拟计算方法，可在自初始状态开始后，每次完成一个时间步长dt的计算，自上而下更新土层的状态，直到开始下一个计算步，从节点1至节点n+1，再进行下一个dt，累计N次。并将达西定律应用于单个元胞的二维非饱和瞬态渗流分析，对非饱和土的渗流分析按深度离散为一定厚度的元胞单位，按单个时间步长进行每个元胞在入渗时的体积含水量、压力水头和渗透系数的计算和状态更新。通过地下水位以上沿深度变化的压力水头分布曲线表征此元胞空间内湿润锋的动态性，借此可进行剩余蓄水承载力的计算。
[0021]作为本专利技术的进一步改进，步骤S1中，确定初始条件包括：
[0022]设有n个元胞，n+1个节点，元胞厚度为l，截面积为A；
[0023]T＝0时，上部流量边界q0，底部为零压力水头h0＝0m，左右不透水；T＝0表示初始状态；
[0024]T>0时，上部流量q，土体饱和渗透系数ks，不饱和系数α，饱和体积含水率θ
s
，残余体积含水率θ
r
。
[0025]此技术方案中，元胞按照一维分析过程进行划分。
[0026]作为本专利技术的进一步改进，步骤S7中，从地表入渗到第一个元胞的流量采用如下公式(3)计算：
[0027][0028]v
e1
＝2k0·
Δh0/l
[0029][0030]其中，v
e1
表示渗透速度，k0表示当前状态的渗透系数，k
s
表示饱和渗透系数，α表示不饱和系数，l表示元胞厚度，A为截面积，θ
s
表示饱和体积含水率，表示第一个元胞最初未更新的体积含水率。
[0031]作为本专利技术的进一步改进，步骤S8中，所述元胞的状态包括如下参数：饱和度、水头、土体饱和渗透系数、残余体积含水率、不饱和系数等。
[0032]作为本专利技术的进一步改进，所述元胞的体积含水率为：
[0033][0034]其中：表示体积含水率，上标表示该元胞的第N
*
次更新(N
*
＜N)，0表示最初未更新，下标表示元胞序号。
[0035]作为本专利技术的进一步改进，步骤S9中，渗流量和压力水头在每个时间步长中采用如下方法计算：
[0036]设第j个元胞中心到第j+1个元胞中心的水头损失为Δh
本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种基于元胞自动机的非饱和地层入渗分析方法，其特征在于：其包括如下步骤：步骤S1，确定初始条件，划分土层的元胞，设定元胞和节点；步骤S2，构造水力特征矩阵，包括位置水头矩阵h
z
、压力水头矩阵h
p
、元胞中心位置水头矩阵(h
z
)
ele
、元胞中心压力水头矩阵(h
p
)
ele
、元胞渗透系数矩阵k、元胞体积含水率系数矩阵、元胞孔隙水体积矩阵；步骤S3，设定时间步长dt，步数N；步骤S4，从节点1开始至节点n+1，依次计算每个节点位置水头h
z,i
和压力水头h
p,i
；步骤S5，初始化每个元胞中心的水力特性：令θ
i
＝θ
r
，计算元胞中心的位置水头(h
z,i
)
ele
和压力水头(h
p,i
)
ele
，并根据式(1)的本构关系计算元胞渗透系数k
i
和元胞体积含水率θ
i
；上式中，k表示当前状态的渗透系数，ks表示饱和渗透系数，θ表示更新状态的体积含水率，θ
s
表示饱和体积含水率，θ
r
表示残余体积含水率，步骤S6，根据式(2)计算第一个元胞的水力梯度Δh0；Δh0＝(h
p,1
+h
z,1
)
‑
[(h
p,1
)
ele
+(h
z,1
)
ele
]＝(0+h
z,1
)
‑
[(h
p,1
)
ele
+(h
z,1
)
ele
]
ꢀꢀꢀ
(2)步骤S7，根据达西定律，计算从地表入渗到第一个元胞的流量；步骤S8，更新第一个元胞的状态；步骤S9，依次计算第二个元胞的水力梯度、从第一个元胞中心入渗到第二个元胞中心的渗流量，更新第二个元胞的状态；重复此步骤，更新第三至第N个元胞；步骤S10，依次根据时间步长，重复迭代计算N次后各元胞中心压力水头，以每个元胞的水力特征状态得到水力特征曲线。2.根据权利要求1所述的基于元胞自动机的非饱和地层入渗分析方法，其特征在于：步骤S1中，...

【专利技术属性】
技术研发人员：苏栋，黄茂隆，蒋鹏，
申请(专利权)人：深圳大学，
类型：发明
国别省市：