悬挂式止水帷幕深基坑降水方案模拟优化方法与流程

1.本发明属于地下工程基坑降水技术领域，具体涉及一种悬挂式止水帷幕深基坑降水方案模拟优化方法。


背景技术：

2.随着城市的建筑不断立体化发展，地下开挖深度不断加深，基坑底部距离含水层顶板越来越近。在进行深基坑工程施工时，了解地下含水层的分布及渗透性分析愈发重要，因此施工前一般要先进行基坑降水。由于深基坑降水会破坏原有的水土应力平衡、造成地下水资源严重浪费，因此，深基坑降水中广泛使用封闭式止水帷幕。而针对含水层厚度较大的地区时，由于封闭式止水帷幕技术难度大、造价成本高、实际工程应用效果不明确，因此采用悬挂式止水帷幕来达到止水的效果。
3.现有技术中，由于地下水渗流场变化及其复杂，具有明显的三维非稳定渗流特性，传统的解析模型难以满足深基坑降水的模拟预测需求，无法有效控制地面变形，减少地下水资源浪费。公开号为cn110263366b的专利公开了一种确定悬挂式止水帷幕插入降水含水层中深度的方法，所述方法根据工程概况、地质信息，结合抽水试验建立三维有限元模型，在给定降水井滤管长度的情况下，通过改变止水帷幕插入降水含水层中的深度，分析止水帷幕内外两侧的水力梯度及坑外地表沉降的变化规律，从而确定止水帷幕插入降水含水层中的合理深度，使之满足工程安全、环境效应及经济性的要求。但其只可计算悬挂式止水帷幕插入降水含水层中的深度，无法模拟完整的降水方案。


技术实现要素：

4.本发明的目的是提供一种悬挂式止水帷幕深基坑降水方案模拟优化方法，为基坑降水方案设计提供可靠的参考意见和计算依据。
5.本发明提供了如下的技术方案：
6.本技术提出一种悬挂式止水帷幕深基坑降水方案模拟优化方法，包括如下步骤：
7.s1.根据研究域的地质、水文地质特征和含水层间的水力联系，建立深基坑降水水文地质概念模型；
8.s2.根据水文地质概念模型，取坐标x、y和z轴的取向与含水层介质各向异性的主方向一致，建立三维非稳定渗流数学模型和含水层土体压缩量的地面沉降数学模型，将三维非稳定渗流数学模型与地面沉降数学模型进行耦合，并通过有限元程序研制、求解，用于通过渗流场、应力场的耦合作用，完整模拟研究域地下水水位和地面沉降的动态趋势；
9.s3.根据非稳定渗流与地面沉降三维耦合数值模型，初始化模型参数，结合已知的地下水位观测值和地面沉降量观测值，通过抽水试验数据进行参数反演，获取满足预测精度和置信度且具有良好收敛性和稳定性的非稳定渗流与地面沉降三维耦合数值模型；
10.s4.进行悬挂式止水帷幕和降水井过滤器不同组合深度的基坑降水与地面沉降模拟计算，结合施工技术难度、成本、环境影响，获取最佳深基坑高效降水技术方案。
11.优先地，步骤s2具体包括如下步骤：
12.s21.以研究域中心海平面为坐标原点，根据研究域水文地质特性、研究域维护连续墙埋藏深度及抽水井滤水管设置位置，将研究域进行垂直分层；
13.s22.根据水均衡原理和达西定律，建立三维非稳定渗流数学模型，具体包括用于描述研究域地下水运动规律的偏微分方程和定解条件，并通过有限单元法求解；
14.s23.建立研究域地下水位下降引起的含水层土体压缩量的地面沉降数学模型，具体包括潜水含水层的弹性变形量、非弹性变形量和承压含水层的弹性变形量、非弹性变形量，并由各含水层的压缩量计算总地面沉降量；
15.s24.综合孔隙度、渗透系数和储水率的动态变化因素，将三维非稳定渗流数学模型和地面沉降数学模型通过水头项进行耦合，建立非稳定渗流与地面沉降三维耦合数值模型，并通过有限元程序研制求解。
16.优先地，步骤s22中，所述三维非稳定渗流数学模型包括控制方程和定解条件，控制方程如下：
[0017][0018]
其中，(x，y，z)为空间位置坐标；k
xx
、k
yy
、k
zz
为各向异性主方向渗透系数；h为点(x，y，z)在t时刻的水头值；w为源汇项；t为时间；ω为研究域；ss为点(x，y，z)处的贮水率，具体为：s
s
＝ρ
w
g(α nβ
w
)，n为土体的孔隙度，α为多孔介质的体积压缩系数，β
w
为水的体积压缩系数，即为水的密度，g为重力加速度；
[0019]
定解条件包括初始条件和至少一种边界条件，将研究域由中心向外进行水平剖分，确定边界条件，其中，研究域内部边界为降水井点，外部侧向边界为第一类边界条件，底边界为第二类边界条件，顶部为自由面边界，
[0020]
所述初始条件为：h(x,y,z,t)|
t＝0
＝h0(x,y,z,t0)，(x,y,z)∈ω，t＝0；
[0021]
所述边界条件包括第一类边界条件、第二类边界条件和自由面边界条件，其中：
[0022]
第一类边界条件为：
[0023]
第二类边界条件为：
[0024][0025]
自由面边界条件为：
[0026]
其中，h0(x，y，z，t0)为点(x，y，z)处的初始水头值，h1(x，y，z，t)为点(x，y，z)在t时刻第一类边界上的已知水头值，q(x，y，z，t)为t时刻第二类边界上单位面积的补给量，cos(n，x)、cos(n，y)、cos(n，z)分别为第二类边界外法线方向与坐标轴方向夹角的余弦，μ为饱和差或给水度，表示自由面改变单位高度下从含水层单位截面积上吸收或排出的水量，n
z
为自由面外法线方向向量n＝{n
x
,n
y
,n
z
}的第3个分量，γ1、γ2、γ3分别为第一类边界、第二类边界和自由面边界。
[0027]
优先地，步骤s22中，具体包括如下步骤：
[0028]
s221.将研究域离散化，即将研究域剖分成多个单元体，单元体与单元体之间通过单元体顶点一节点相联系；
[0029]
s222.在标准的单元体内建立任意一点的坐标与节点之间相互关系的形函数n
i
，既然满足在节点i，ni＝1，在其他节点ni＝0，在单元体内任意点上的形函数之和等于1；
[0030]
s223.对于不规则的单元体利用等参变换，将单元体坐标的变换与单元体内水头插值函数采用相同数目的节点参数和相同的形函数进行变换，在空间上，使(ε，η，ζ)局部坐标系中形状简单、规则的标准单元体，在(x，y，z)整体坐标系中变换为具有曲线或曲面边界的形状复杂的单元；
[0031]
s224.对三维非稳定渗流数学模型的边值问题的求解转化为对泛函i(h)的极值问题的求解，即泛函i(h)的偏导数为零，其中，泛函i(h)的表达式如下：
[0032][0033]
单元泛函i
e
(h)的表达式如下：
[0034][0035]
s225.对单元泛函i
e
(h)偏导数为零，得到分别对求偏导，并整理得出单元支配表达式为：集成后可得整体有限元支配表达式为：并转化为时间差分格式，
[0036]
其中，[k]为总体渗透矩阵，由单元渗透矩阵[k]
e
集合而成；[s]为储水矩阵，由各单元储水矩阵[s]
e
迭加而得；[g]为流量补给矩阵，由自由面穿过的各单元流量补给矩阵[g]
e
迭加而得；{f}为节点流量矩阵，由各单元水量矩阵{f}
e
迭加而得；{h}为待求节点水头向量。
[0037]
优先地，步骤s225中，若是内部单元，则
[0038]
若为自由面边界单元，则若为有压渗流场，则针对自由面边界，
[0039]
优先地，步骤s22中，边界条件与控制方程联立前需对选取的边界调节进行处理，其中，对降水井点采用置大数法进行处理，对自由面边界采用改进复合单元渗透矩阵调整法进行处理以稳定三维非稳定渗流数学模型，其中，将自由面边界的几何区域分为虚单元、
实单元和过渡单元并通过迭代求解，其中对于过渡单元采用连续的折现型罚函数法将计算区域拓展至整个单元，使待定自由面边界固化。
[0040]
优先地，整体有限元支配表达式的时间差分格式包括显示差分格式、中心差分格式和隐式差分格式，其中，隐式差分格式是无条件稳定的。
[0041]
优先地，步骤s23中，
[0042]
潜水含水层的弹性变形量为：δb＝
‑
δh(1
‑
n n
w
)μ
ske
b0＝
‑
δhμ
fe
；
[0043]
潜水含水层的非弹性变形量为：δb
*
＝
‑
δh(1
‑
n n
w
)μ
skv
b0＝
‑
δhμ
fv
；
[0044]
承压含水层的弹性变形量为：δb＝
‑
δhμ
ske
b0＝
‑
δhμ
fe
；
[0045]
承压含水层的非弹性变形量为：δb
*
＝
‑
δhμ
skv
b0＝
‑
δhμ
fv
；
[0046]
其中，μ
ske
为土体骨架成分的弹性储水率，μ
skv
为土体骨架成分的非弹性储水率，μ
fe
为土体骨架成分的弹性储水因子，μ
fv
为土体骨架成分的非弹性储水因子；μ
fe
、μ
fv
为将含水层及包含于其中的弱含水层视为一个整体后，在整个含水层规模上的等效参数；
[0047]
而μ
ske
与μ
skv
的模型为：
[0048]
其中，g为剪切模量，v为泊松比，g为重力常数，c
c
为土体的压缩系数，σ0′
为初始有效应力，e0为初始孔隙比，ρ
w
为水的密度，σ
′
max
为土层的前期固结应力。
[0049]
优先地，步骤s24的具体步骤如下：
[0050]
s241.结合孔隙流体在多孔介质中的流动规律及其对多孔介质本身的变形或者强度造成的影响进行土体非线性固结模型改进，即考虑多孔介质内应力场与渗流场之间的相互耦合作用来处理土体的非线性渗流问题，其中：
[0051]
固结过程中的孔隙度n为：其中，土层厚度为m
i
，由于降水引起的垂直沉降量为s
i
；
[0052]
渗透系数k的变化公式为：其中，n0为初始孔隙度；k0为初始渗透系数；
[0053]
地下水渗流连续性方程的右端项中的储水率μ
s
为：μ
s
＝γ
w
(α nβ
w
)，其中，γ
w
为水的重度；α为土的体积压缩系数；n为孔隙度；β
w
为水的体积压缩系数；
[0054]
根据固结试验曲线e
‑
lgσ
′
的斜率c
c
代入储水率计算公式中，可得储水率μ
s
随e变化的关系为：
[0055]
其中，e0为初始孔隙比；σ
′
为有效应力；
[0056]
s242.将整体有限元支配表达式的时间差分格式中的隐式差分格式中的系数矩阵进行处理，获取对称正定方程组，并结合边界条件，将对称正定方程组的系数矩阵进行预处理，用于减少等价问题的条件数，之后采用共轭梯度法初始化参数后进行迭代求解，获取初
始条件下的渗流场水头，并计算任意时刻渗流场水头，用于提高收敛速度和稳定性；
[0057]
s243.采用fortran95语言将上述整个计算过程编制成计算机程序，形成计算机模型，所述计算机模型包括主程序和子程序，所述主程序用于读取文件、输入参数和程序结束时输出计算值，子程序包括形成主元素序号指示矩阵模块、形成整体劲度矩阵模块、形成常数项模块、置大数法修正已知水头模块、方程求解模块、沉降计算模块和参数变化模块。
[0058]
优先地，步骤s242中，具体步骤如下：
[0059]
s2421.整体有限元支配表达式的时间差分格式中的隐式差分格式为：
[0060]
对其系数矩阵进行处理，令则隐式差分格式可简化为：[a]{h}
t δt
＝{b}，其中，[a]称为总刚度矩阵；{b}称为常数项；
[0061]
s2422.假设a为对称正定矩阵，则对称正定方程组为ax＝b；
[0062]
s2423.采用对称逐步超松驰迭代法的分裂矩阵作为预处理矩阵m，具体公式为：m＝(2
‑
ω)
‑1(d/ω l)(d/ω)
‑1(d/ω l)
t
，
[0063]
其中，d为a的对角阵；l为a的严格下三角矩阵；ω为松驰因子且0<ω<2；
[0064]
s2424.采用预处理共轭梯度法进行迭代求解，具体迭代格式为：
[0065][0066]
其中，预处理共轭梯度法每步迭代的主要计算量为(2r
a
6)n次乘法运算，r
a
为对称正定矩阵a的各行非零元素个数的平均值。
[0067]
本发明的有益效果是：在充分分析、掌握研究域工程地质、水文地质和环境地质条件的基础上，耦合地下水渗流场与地面沉降土体变形量，建立超深基坑降水地下水三维非稳定渗流与地面沉降耦合模型，并采用有限单元法进行求解，通过抽水试验反演，使模型的识别满足一定精度和置信度，并提高模型的计算速度和稳定性，为基坑降水方案设计提供可靠的计算依据。
附图说明
[0068]
附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：
[0069]
图1是本发明的方法流程图。
具体实施方式
[0070]
本技术提出一种悬挂式止水帷幕深基坑降水方案模拟优化方法，包括如下步骤：
[0071]
s1.根据研究域的地质、水文地质特征和含水层间的水力联系，建立深基坑降水水文地质概念模型。
[0072]
s2.根据水文地质概念模型，取坐标x、y和z轴的取向与含水层介质各向异性的主方向一致，建立三维非稳定渗流数学模型和含水层土体压缩量的地面沉降数学模型，将三维非稳定渗流数学模型与地面沉降数学模型进行耦合，并通过有限元程序研制、求解，用于通过渗流场、应力场的耦合作用，完整模拟研究域地下水水位和地面沉降的动态趋势。
[0073]
步骤s2具体包括如下步骤：
[0074]
s21.以研究域中心海平面为坐标原点，根据研究域水文地质特性、研究域维护连续墙埋藏深度及抽水井滤水管设置位置，将研究域进行垂直分层。
[0075]
s22.根据水均衡原理和达西定律，建立三维非稳定渗流数学模型，具体包括用于描述研究域地下水运动规律的偏微分方程和定解条件，并通过有限单元法求解。
[0076]
步骤s22中，三维非稳定渗流数学模型包括控制方程和定解条件，控制方程如下：
[0077][0078]
其中，(x，y，z)为空间位置坐标；k
xx
、k
yy
、k
zz
为各向异性主方向渗透系数；h为点(x，y，z)在t时刻的水头值；w为源汇项；t为时间；ω为研究域；ss为点(x，y，z)处的贮水率，具体为：s
s
＝ρ
w
g(α nβ
w
)，n为土体的孔隙度，α为多孔介质的体积压缩系数，β
w
为水的体积压缩系数，即为水的密度，g为重力加速度；
[0079]
定解条件包括初始条件和至少一种边界条件，将研究域由中心向外进行水平剖分，确定边界条件，其中，研究域内部边界为降水井点，外部侧向边界为第一类边界条件，底边界为第二类边界条件，顶部为自由面边界，
[0080]
初始条件为：h(x,y,z,t)|
t＝0
＝h0(x,y,z,t0)，(x,y,z)∈ω，t＝0；
[0081]
边界条件包括第一类边界条件、第二类边界条件和自由面边界条件，其中：
[0082]
第一类边界条件为：
[0083]
第二类边界条件为：
[0084][0085]
自由面边界条件为：
[0086]
其中，h0(x，y，z，t0)为点(x，y，z)处的初始水头值，h1(x，y，z，t)为点(x，y，z)在t时刻第一类边界上的已知水头值，q(x，y，z，t)为t时刻第二类边界上单位面积的补给量，cos(n，x)、cos(n，y)、cos(n，z)分别为第二类边界外法线方向与坐标轴方向夹角的余弦，μ为饱和差或给水度，表示自由面改变单位高度下从含水层单位截面积上吸收或排出的水量，n
z
为自由面外法线方向向量n＝{n
x
,n
y
,n
z
}的第3个分量，γ1、γ2、γ3分别为第一类边界、第二类边界和自由面边界。
[0087]
步骤s22中，边界条件与控制方程联立前需对选取的边界调节进行处理，其中，对降水井点采用置大数法进行处理，对自由面边界采用改进复合单元渗透矩阵调整法进行处理以稳定三维非稳定渗流数学模型，其中，将自由面边界的几何区域分为虚单元、实单元和过渡单元并通过迭代求解，其中对于过渡单元采用连续的折现型罚函数法将计算区域拓展至整个单元，使待定自由面边界固化。
[0088]
步骤s22中，具体包括如下步骤：
[0089]
s221.将研究域离散化，即将研究域剖分成多个单元体，单元体与单元体之间通过单元体顶点一节点相联系；
[0090]
s222.在标准的单元体内建立任意一点的坐标与节点之间相互关系的形函数n
i
，既然满足在节点i，ni＝1，在其他节点ni＝0，在单元体内任意点上的形函数之和等于1；
[0091]
s223.对于不规则的单元体利用等参变换，将单元体坐标的变换与单元体内水头插值函数采用相同数目的节点参数和相同的形函数进行变换，在空间上，使(ε，η，ζ)局部坐标系中形状简单、规则的标准单元体，在(x，y，z)整体坐标系中变换为具有曲线或曲面边界的形状复杂的单元；
[0092]
s224.对三维非稳定渗流数学模型的边值问题的求解转化为对泛函i(h)的极值问题的求解，即泛函i(h)的偏导数为零，其中，泛函i(h)的表达式如下：
[0093][0094]
单元泛函i
e
(h)的表达式如下：
[0095][0096]
s225.对单元泛函i
e
(h)偏导数为零，得到分别对求偏导，并整理得出单元支配表达式为：
[0097]
集成后可得整体有限元支配表达式为：并转化为时间差分格式，
[0098]
其中，[k]为总体渗透矩阵，由单元渗透矩阵[k]
e
集合而成；[s]为储水矩阵，由各单元储水矩阵[s]
e
迭加而得；[g]为流量补给矩阵，由自由面穿过的各单元流量补给矩阵[g]
e
迭加而得；{f}为节点流量矩阵，由各单元水量矩阵{f}
e
迭加而得；{h}为待求节点水头向量。整体有限元支配表达式的时间差分格式包括显示差分格式、中心差分格式和隐式差分格式，其中，隐式差分格式是无条件稳定的。
[0099]
步骤s225中，若是内部单元，则
[0100]
若为自由面边界单元，则若为有压渗流场，则针对自由面边界，
[0101]
s23.建立研究域地下水位下降引起的含水层土体压缩量的地面沉降数学模型，具体包括潜水含水层的弹性变形量、非弹性变形量和承压含水层的弹性变形量、非弹性变形量，并由各含水层的压缩量计算总地面沉降量，其中，
[0102]
潜水含水层的弹性变形量为：δb＝
‑
δh(1
‑
n n
w
)μ
ske
b0＝
‑
δhμ
fe
；
[0103]
潜水含水层的非弹性变形量为：δb
*
＝
‑
δh(1
‑
n n
w
)μ
skv
b0＝
‑
δhμ
fv
；
[0104]
承压含水层的弹性变形量为：δb＝
‑
δhμ
ske
b0＝
‑
δhμ
fe
；
[0105]
承压含水层的非弹性变形量为：δb
*
＝
‑
δhμ
skv
b0＝
‑
δhμ
fv
；
[0106]
其中，μ
ske
为土体骨架成分的弹性储水率，μ
skv
为土体骨架成分的非弹性储水率，μ
fe
为土体骨架成分的弹性储水因子，μ
fv
为土体骨架成分的非弹性储水因子；μ
fe
、μ
fv
为将含水层及包含于其中的弱含水层视为一个整体后，在整个含水层规模上的等效参数；
[0107]
而μ
ske
与μ
skv
的模型为：
[0108]
其中，g为剪切模量，v为泊松比，g为重力常数，c
c
为土体的压缩系数，σ0′
为初始有效应力，e0为初始孔隙比，ρ
w
为水的密度，σ
′
max
为土层的前期固结应力。
[0109]
s24.综合孔隙度、渗透系数和储水率的动态变化因素，将三维非稳定渗流数学模型和地面沉降数学模型通过水头项进行耦合，建立非稳定渗流与地面沉降三维耦合数值模型，并通过有限元程序研制求解。
[0110]
步骤s24的具体步骤如下：
[0111]
s241.结合孔隙流体在多孔介质中的流动规律及其对多孔介质本身的变形或者强度造成的影响进行土体非线性固结模型改进，即考虑多孔介质内应力场与渗流场之间的相互耦合作用来处理土体的非线性渗流问题，其中：
[0112]
固结过程中的孔隙度n为：其中，土层厚度为m
i
，由于降水引起的垂直沉降量为s
i
；
[0113]
渗透系数k的变化公式为：其中，n0为初始孔隙度；k0为初始渗透系数；
[0114]
地下水渗流连续性方程的右端项中的储水率μ
s
为：μ
s
＝γ
w
(α nβ
w
)，其中，γ
w
为水的重度；α为土的体积压缩系数；n为孔隙度；β
w
为水的体积压缩系数；
[0115]
根据固结试验曲线e
‑
lgσ
′
的斜率c
c
代入储水率计算公式中，可得储水率μ
s
随e变化的关系为：
[0116]
其中，e0为初始孔隙比；σ
′
为有效应力。
[0117]
s242.将整体有限元支配表达式的时间差分格式中的隐式差分格式中的系数矩阵进行处理，获取对称正定方程组，并结合边界条件，将对称正定方程组的系数矩阵进行预处理，用于减少等价问题的条件数，之后采用共轭梯度法初始化参数后进行迭代求解，获取初始条件下的渗流场水头，并计算任意时刻渗流场水头，用于提高收敛速度和稳定性。
[0118]
步骤s242中，具体步骤如下：
[0119]
s2421.整体有限元支配表达式的时间差分格式中的隐式差分格式为：
[0120][0121]
对其系数矩阵进行处理，令则隐式差分格式可简化为：[a]{h}
t δt
＝{b}，其中，[a]称为总刚度矩阵；{b}称为常数项；
[0122]
s2422.假设a为对称正定矩阵，则对称正定方程组为ax＝b。
[0123]
s2423.采用对称逐步超松驰迭代法的分裂矩阵作为预处理矩阵m，具体公式为：m＝(2
‑
ω)
‑1(d/ω l)(d/ω)
‑1(d/ω l)
t
，
[0124]
其中，d为a的对角阵；l为a的严格下三角矩阵；ω为松驰因子且0<ω<2。
[0125]
s2424.采用预处理共轭梯度法进行迭代求解，具体迭代格式为：
[0126][0127]
其中，预处理共轭梯度法每步迭代的主要计算量为(2r
a
6)n次乘法运算，r
a
为对称正定矩阵a的各行非零元素个数的平均值。
[0128]
s243.采用fortran95语言将上述整个计算过程编制成计算机程序，形成计算机模型，计算机模型包括主程序和子程序，主程序用于读取文件、输入参数和程序结束时输出计算值，子程序包括形成主元素序号指示矩阵模块、形成整体劲度矩阵模块、形成常数项模块、置大数法修正已知水头模块、方程求解模块、沉降计算模块和参数变化模块。
[0129]
s3.根据非稳定渗流与地面沉降三维耦合数值模型，初始化模型参数，结合已知的地下水位观测值和地面沉降量观测值，通过抽水试验数据进行参数反演，获取满足预测精度和置信度且具有良好收敛性和稳定性的非稳定渗流与地面沉降三维耦合数值模型。
[0130]
s4.进行悬挂式止水帷幕和降水井过滤器不同组合深度的基坑降水与地面沉降模拟计算，结合施工技术难度、成本、环境影响，获取最佳深基坑高效降水技术方案。
[0131]
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。