An improved numerical manifold method for solving heterogeneous seepage problem
-
摘要: 精确的流速场计算是模拟非均质多孔介质渗流、溶质运移以及污染物迁移等过程的基础。常规方法计算得到流速场精度低、连续性差。高阶流形方法可以显著提高流速场计算精度,但对于高水力梯度问题,流速场依然是不连续的,且常常伴随着线性相关问题,导致求解精度降低,甚至求解失败。针对这一问题,采用两种改进的权函数构造流形单元上的总体近似函数,编制相应的程序,消除高阶覆盖函数带来的线性相关问题,且分别得到节点连续和全局连续的流速/梯度解。通过若干算例对改进权函数程序的收敛性和精度进行分析。Abstract: An accurate velocity filed is important for simulating Darcy flow, solution transport and contaminant migration in heterogeneous porous media. Generally, the velocity results calculated by the conventional numerical methods have low precision and poor continuity. To overcome this defect, two improved weight functions are used to construct the global approximation on manifold elements and two improved NMM codes are developed to eliminate the linear dependence problems caused by the higher-order overburden functions, and the velocity/gradient solutions of node continuity and global continuity are obtained, respectively. Furthermore, the convergence and accuracy of the improved weight function NMM codes are analyzed through several examples.
-
Keywords:
- numerical manifold method /
- heterogeneity /
- Darcy velocity /
- weight function /
- high order
-
0. 引言
静载试验为目前公认的最可靠的桩基承载力测试方法,但由于自身的构造特点及技术局限,堆载法、锚桩法和自平衡法均存有不同程度的缺陷[1]。随着岩土锚固技术的逐步完善,基桩自锚试验成为一种新型检测基桩承载力的方法,利用桩底锚固技术为试验提供反力,基桩、地基土和桩底锚杆构成相对封闭的内力传递系统,由于桩底土(岩)层在基桩和锚杆荷载作用下发生位移及变形,导致基桩承载性能发生变化。自锚试验所得Q-S曲线与静载试验存有多大的差异以及相互之间如何实现转化,需进一步分析研究。
目前,学者们针对自锚试验桩-土-锚系统中的锚杆荷载传递机理研究成果较为丰硕,主要有现场试验、数值模拟以及理论等方面的研究[2-6],然而,对于三者耦合作用的研究较为少见。自锚试验中桩-土(岩)界面荷载传递模型对基桩承载性状有重要影响,目前针对基桩承载性状的计算方法主要有：弹性理论法、荷载传递法、剪切位移法、数值计算法等[7-21]。基于上述研究背景,本文利用室内模型试验和数值模拟结果提出一种非线性的桩-土荷载传递模型,进而基于桩土界面荷载传递法实现基桩Q-S曲线的转化。
1. 模型试验及数值模拟
1.1 室内模型试验
本次试验共分为锚固深度L0依次为12.5, 25, 50 cm的3组自锚试验以及一组静载试验(对比试验)。其中,模型桩嵌入混合料层的深度H0均为20 cm,桩长90 cm,锚杆采用直径为φ8的钢绞线模拟,试验装置全貌如图1所示,桩周土和混合料其他具体参数详见文献[22]。
自锚试验与静载试验中基桩的Q-S曲线变化如图2所示。
![]()
图 2 自锚试验与静载试验Q-S曲线图Figure 2. Q-S curves of self-anchoring tests and static load tests从图2可知,相同竖向荷载下,静载试验中基桩的沉降大于自锚试验,本次静载试验基桩的极限承载力为7.2 kN,自锚试验中基桩的极限承载力分别为9.9 kN(12.5 cm)、9.3 kN(25 cm)、8.6 kN(50 cm),自锚试验中,桩底锚杆起到桩锚强化效应。当锚杆受到向上的拉拔荷载,桩端锚固体有向上的运动趋势,一方面类似于给基桩施加了向上的预应力,另一方面桩端土体的压力大于静止土压力,土体位移不为零,致使桩端土挤密压实,锚固深度越小,桩端挤压效果愈加明显。相同荷载下,锚固荷载引起土层位移的竖向影响深度以及每一深度的位移值都在增大,导致基桩位移总变化量较小。随着锚杆锚固深度的减小,基桩的极限承载力逐渐增大,但并非线性增加。在桩-土-锚共同作用下,当锚固深度逐渐增大时,桩端受到的影响逐渐远离锚固荷载的应力扩散范围,其极限承载力值愈发接近桩底无锚时的情况。
1.2 数值模拟
桩底有无锚杆对Q-S曲线的刚性系数和极限承载力均存有影响,为了进一步探讨锚固深度对基桩承载性状的影响,利用有限元进行数值模拟,结合室内模型试验。桩周土体采用二维轴对称可变形的壳体,计算范围以桩轴线为对称轴,左右岩土体宽各50 cm,桩、填料参数与室内模型试验均一致。在试验基础上增加6组基桩试验,锚固深度依次为45, 40, 35, 30, 20, 15 cm,模拟所用的材料参数如表1所示,室内实测与数值模拟Q-S曲线如图3所示。
表 1 材料参数Table 1. Material parameters单元 弹性模量/GPa 泊松比v 密度/(g·cm-3) 黏聚力/kPa 内摩擦角/(°) 黄土 0.00282 0.26 1.560 14.91 23.7 混合料 0.35 0.30 1.736 146.3 38.0 桩 71 0.20 2.1 — — 锚杆 195 0.25 2.6 — — 锚固体 210 0.28 2.8 — — 从图3可知,实测值与模拟值较为接近,曲线吻合度较高,因此,可用该模型探究其他锚固深度对基桩性状的影响。其中,当荷载等级为2.1 kN时,锚固深度分别为12.5 cm和25 cm的竖向位移云图如图4所示,对于桩-土-锚共同体系,在多种力耦合作用情况下,桩端出现了塑变区和成拱区两个性质不同的区域,受拉拔荷载的影响,塑变区和成拱区延伸到桩端附近以上几倍的桩径范围之外,同时,由于上覆土压力的限制作用使其形成一稳定区域。桩端土体通常强度较高,在高应力作用下,桩端附近桩土相对位移增大,且使桩端附近范围内土体应力增大。自锚试验中,在加载末期,桩侧摩阻力逐渐达到极限值,对应于基桩顶部荷载的增量,桩端阻力增量所占比例逐渐增大,使得桩端出现塑限区并且逐渐扩展,基桩急剧下沉,导致桩端发生刺入破坏。
![]()
图 4 不同锚固深度下的竖向位移云图(2.1 kN)Figure 4. Cloud diagram of vertical displacement at different anchoring depths(2.1 kN)试验结果得静载试验与自锚试验的弹性斜率Kp与极限承载力Qu的值如表2所示,其中,数值解极限承载力采用切线交汇法确定。
表 2 弹性斜率与极限承载力Table 2. Elastic slopes and ultimate bearing capacities参数 静载 L0=12.5 L0=15 L0=20 L0=25 Qu/kN 7.2 9.9 9.76 9.54 9.3 Kp/(kN·mm-1) 0.83 3.73 3.35 2.66 2.09 L0=30 L0=35 L0=40 L0=45 L0=50 Qu/kN 9.09 8.80 8.57 8.42 8.6 Kp(kN·mm-1) 1.78 1.55 1.31 1.22 1.16 注： 表中自锚试验的L0(锚固深度)单位为cm。利用模拟结果得Quz(自锚试验极限承载力)和Quc(静载试验极限承载力)的比值以及Kpz(自锚试验Q-S曲线弹性斜率)和Kpc(静载试验Q-S曲线弹性斜率)的比值与基桩锚固深度的关系如图5所示。当锚固深度取值趋近于无穷时,Quz/Quc=1.078, Kpz/Kpc=1.046,其变化规律与图2揭示的机理类似,即当锚固深度足够大时,桩端处于锚固荷载的应力扩散范围之外,锚固荷载引起的土层位移达不到桩端或小到可以忽略。然而在实际工程中,软土地区基桩埋深较大,桩端一般达至基岩顶面或桩端位于较密的持力层,若锚固深度过大,施工难度不但大幅度提升,而且也不经济。因而基桩锚固深度值不可能过大,其端部通常位于锚固区的影响范围之内,锚杆的存在必定引起基桩承载特性的变化,因此,研究自锚试验与静载试验结果之间的差异,将自锚试验结果向静载试验转化对实际工程具有重要意义。
![]()
图 5 锚固深度与极限承载力及弹性斜率比值关系图Figure 5. Relationship among anchoring depth, ultimate bearing capacity and elastic slope从图5可知,上述关系满足式(1)的指数关系：
(1) 2. 自锚试验数据的转化
2.1 基于实测数据的地基土参数分析
为准确获得基桩的桩土荷载关系,需要在实测资料的基础上加以分析处理,基于试验得出基桩Q-S曲线及每级荷载下的桩身轴力,进而获得桩土荷载传递资料。假设在荷载Qi作用下,桩身各测点的轴力为Ni,桩顶位移s1为轴力N1所对应的位移,如图6所示。
根据静力平衡方程,可得两测点间的平均剪应力为
(2) 式中,r为桩的半径,hi为两测点间的距离,ΔNi为第i个测点与i+1测点的轴力差值。两测点间平均应变可表示为
(3) 式中,E为弹性模量,A为横截面积。两测点间桩身变形量为
(4) 各测点桩身位移为
(5) 桩土平均位移计算式可表示为
(6) 式中,
基于能量原理可求解出锚固荷载作用下的桩底土体位移[6],模型桩(锚固深度25 cm)底土层位移分布如图7所示。
针对锚固深度25 cm的基桩桩土界面荷载传递关系进行非线性拟合分析,如图8所示,桩土界面荷载传递模型用对数曲线y=a-bln(x+c)拟合精度高,由于桩底有无锚杆对桩土相对位移与基桩摩阻力的荷载传递关系无关[6],因而基桩的桩土荷载传递关系可用于计算桩底无锚时的Q-S曲线,对数荷载传递模型表达式为
(7) 式中,Δs为桩土相对位移,a, b, c为待求参数。
特别指出,在实际工程中利用实测数据估算桩土相对位移时,一般假定桩土之间无相对滑移,通常所说的桩土相对位移只是桩相对较远处不发生位移的土的位移。因此,现场实测得到的某深度处的桩土相对位移,其实就是指该深度处的桩身位移[23],故第i段的桩土相对位移值Δsi如下所示：
(8) 式中,st为桩顶沉降,Lj为第j段桩的桩长,εj和εj+1分别表示第j和第j+1的断面钢筋应变。
2.2 荷载沉降计算方法
采用Matlab软件编制程序,结合下述迭代算法可求解基桩承载性状[24],其主要计算步骤如下：
(1)沿桩身长度方向将桩长分为n段；
(2)假设作用在桩顶的荷载Pt=0,因此桩段1的荷载Pt1=0,假定桩顶有向下的小位移St1,则桩段1有向下的位移zt1,此时St1=zt1,起初桩土相对位移sc1=zt1;
(3)利用式(7)可求解出桩段1的剪应力τ1;
(4)此时可求解出桩段1的总侧摩阻力T1, T1=2πrl1τs1;
(5)根据桩段力的平衡可得：桩段1底部的力Pb1=Pt1-T1;
(6)桩段1中点处的弹性压缩量可表示为：zc1=
(7)桩段1的修正相对位移sc1′=sc1+zc1;
(8)比较桩段1内的桩身修正相对位移sc1′与初始值sc1的差,如果|sc1′-sc1|≥1×10-6m,重复步骤(2)~(7),直至|sc1′-sc1|≤1×10-6m为止；
(9)利用式zb1=
(10)更新Pt1和zt1: Pt1=Pb1, zt1=zb1,重复步骤(2)~(9)直至桩段n,然后Pt1和zt1分别储存为P1和s1;
(11)假定一系列桩顶位移St…Stn,重复步骤(2)~(10),直至得到一系列P1…Pn和s1…sn。
根据上述荷载沉降计算,可实现自锚试验中基桩Q-S曲线的转化,转化结果如图9所示。从图可知,采用对数荷载传递模型模拟桩侧阻力和桩土相对位移以及桩端阻力与桩端位移间的关系,运用迭代法转化后的基桩Q-S曲线与静载试验结果较为吻合,与基于能量原理的桩身能量差分方程求解结果较为接近,可见此方法较为合理,计算过程简便可行,可作为自锚试验测试结果向静载试验转化的一种有效方法。
![]()
图 9 自锚试验基桩Q-S曲线转化图Figure 9. Transformation of Q-S of foundation piles from self-anchoring tests3. 结论
本文结合室内模型试验以及数值模拟研究,基于桩土界面的荷载传递实现自锚试验Q-S曲线的转化,其主要结论为
(1)相同竖向荷载下,静载试验中基桩的沉降量大于自锚试验,其极限承载力小于自锚试验；随着锚固深度的减小,基桩的极限承载力逐渐增大,但并非线性增加。
(2)自锚试验中,随着锚固深度逐渐增大,桩端受到的影响逐渐远离锚固荷载的应力扩散范围,Quz/Quc和Kpz/Kpc的值分别与基桩锚固深度呈指数关系。
(3)采用对数荷载传递模型模拟桩侧阻力和桩土相对位移以及桩端阻力与桩端位移间的关系,运用迭代法可较好地实现自锚试验基桩Q-S曲线向静载试验的转化。
-
![]()
图 11 y=100处流速计算结果与解析解的对比图
Figure 11. Comparison between numerical and analytical solutions
![]()
图 16 x方向流速等值线图局部放大示意图
Figure 16. Partial enlargement view of Darcy velocity in x direction
![]()
图 19 不同均匀程度对水头等值线的影响
Figure 19. Contour maps of hydraulic head under different heterogeneities
表 1 不同NMM程序总体系数矩阵的亏秩数与总计算自由度
Table 1 Number of nullity and total degrees of freedom
kh P0-W0 P1-W0 P1-WA P1-WB 4 0/25 10/75 0/75 0/75 6 0/49 14/147 0/147 0/147 8 0/81 18/243 0/243 0/243 10 0/121 22/363 0/363 0/363 12 0/169 26/507 0/507 0/507 注： P为覆盖函数的阶次,其中P0为常数,P1为一阶泰勒展开多项式。
下载: 导出CSV
-
[1] 薛禹群, 谢春红. 地下水数值模拟[M]. 北京: 科学出版社, 2007: 175-178. XUE Yu-qun, XIE Chun-hong. Numerical Simulation for Groundwater[M]. Beijing: Science Press, 2007: 175-178. (in Chinese)
[2] BATU V. A finite element dual mesh method to calculate Nodal Darcy velocities in nonhomogeneous and anisotropic aquifers[J]. Water Resources Research, 1984, 20(11): 1705-1717. doi: 10.1029/WR020i011p01705
[3] YEH G T. On the computation of Darcian velocity and mass balance in the finite element modeling of groundwater flow[J]. Water Resources Research, 1981, 17(5): 1529-1534. doi: 10.1029/WR017i005p01529
[4] 谢春红, 赵文良, 张天岭, 等. 地下水不稳定渗流达西速度计算新方法[J]. 岩土工程学报, 1996, 18(1): 68-74. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC601.009.htm XIE Chun-hong, ZHAO Wen-liang, ZHANG Tian-ling, et al. A new method for calculating the Darcy velocity of unstable groundwater seepage[J]. Chinese Journal of Geotechnical Engineering, 1996, 18(1): 68-74. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC601.009.htm
[5] 张湘伟, 章争荣, 吕文阁, 等. 数值流形方法研究及应用进展[J]. 力学进展, 2010, 40(1): 1-12. https://www.cnki.com.cn/Article/CJFDTOTAL-LXJZ201001003.htm ZHANG Xiang-wei, ZHANG Zheng-rong, LÜ Wen-ge, et al. Advances and perspectives in numerical manifold method and its applications[J]. Advances in Mechanics, 2010, 40(1): 1-12. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-LXJZ201001003.htm
[6] MA G W, AN X M, HE L. The numerical manifold method: a review[J]. International Journal of Computational Methods, 2010, 7(1): 1-32. doi: 10.1142/S0219876210002040
[7] ZHANG Z R, ZHANG X W, YAN J H. Manifold method coupled velocity and pressure for Navier-Stokes equations and direct numerical solution of unsteady incompressible viscous flow[J]. Computers and Fluids, 2010, 39(8): 1353-1365. doi: 10.1016/j.compfluid.2010.04.005
[8] 陈远强, 杨永涛, 郑宏, 等. 饱和-非饱和渗流的数值流形法研究与应用[J]. 岩土工程学报, 2019, 41(2): 338-347. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201902014.htm CHEN Yuan-qiang, YANG Yong-tao, ZHENG Hong, et al. Saturated-unsaturated seepage by numerical manifold method[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(2): 338-347. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201902014.htm
[9] 李伟, 郑宏. 基于数值流形法的渗流问题边界处理新方法[J]. 岩土工程学报, 2017, 39(10): 1867-1873. doi: 10.11779/CJGE201710015 LI Wei, ZHENG Hong. New boundary treatment for seepage flow problem based on numerical manifold method[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(10): 1867-1873. (in Chinese) doi: 10.11779/CJGE201710015
[10] WANG Y, HU M S, ZHOU Q L, et al. Energy-work-based numerical manifold seepage analysis with an efficient scheme to locate the phreatic surface[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2014, 38(15): 1633-1650. doi: 10.1002/nag.2280
[11] WANG Y, HU M S, ZHOU Q L, et al. A new second-order numerical manifold method model with an efficient scheme for analyzing free surface flow with inner drains[J]. Applied Mathematical Modelling, 2016, 40(2): 1427-1445. doi: 10.1016/j.apm.2015.08.002
[12] SHI G H. Manifold methods of material analysis[C]//Transactions of the 9th Army Conference on Applied Mathematics and Computing. Minneapolis, 1991.
[13] 唐旭海, 郑超, 吴圣川, 等. 节点应力连续的四边形单元[J]. 应用数学和力学, 2009, 30(12): 1427-1439. doi: 10.3879/j.issn.1000-0887.2009.12.004 TANG Xu-hai, ZHENG Chao, WU Sheng-chuan, et al. A novel four -node quadrilateral element with continuous nodal stress[J]. Applied Mathematics and Mechanics, 2009, 30(12): 1427-1439. (in Chinese) doi: 10.3879/j.issn.1000-0887.2009.12.004
[14] YANG Y T, SUN G H, ZHENG H, et al. A four-node quadrilateral element fitted to numerical manifold method with continuous nodal stress for crack analysis[J]. Computers & Structures, 2016, 177: 69-82.
[15] 王勖成. 有限单元法[M]. 北京: 清华大学出版社, 2003: 112-113. WANG Xu-cheng. Finite Element Method[M]. Beijing: Tsinghua University Press, 2003: 112-113. (in Chinese)
[16] YANG Y T, SUN G H, ZHENG H. A high-order numerical manifold method with continuous stress/strain field[J]. Applied Mathematical Modelling, 2020, 78: 576-600. doi: 10.1016/j.apm.2019.09.034
[17] 谢一凡. 改进多尺度有限单元法求解二维地下水流问题[D]. 南京: 南京大学, 2015: 45-46. XIE Yi-fan. Modified Multiscale Finite Element Method for 2-D Groundwater Flow Problems[D]. Nanjing: Nanjing University, 2015: 45-46. (in Chinese)
[18] AN X M, LI L X, MA G W, et al. Prediction of rank deficiency in partition of unity-based methods with plane triangular or quadrilateral meshes[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(5/6/7/8): 665-674.
[19] 赵文凤, 谢一凡, 吴吉春. 一种模拟节点达西渗透流速的双重网格多尺度有限单元法[J]. 岩土工程学报, 2020, 42(8): 1474-1481. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202008018.htm ZHAO Wen-feng, XIE Yi-fan, WU Ji-chun. A dual-mesh multiscale finite element method for simulating nodal Darcy velocities in aquifers[J]. Chinese Journal of Geotechnical Engineering, 2020, 42(8): 1474-1481. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC202008018.htm
[20] GRIFFITHS D V, FENTONT G A. The Random Finite Element Method (RFEM) in Steady Seepage Analysis[M]//CISM Courses and Lectures. Vienna: Springer, 2007: 225-241.
[21] LACY S J, PREVOST J H. Flow through porous media: a procedure for locating the free surface[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1987, 11(6): 585-601.
计量
- 文章访问数:
- HTML全文浏览量: 0
- PDF下载量:
