长江科学院院报

长江科学院院报 >

2024 , Vol. 41 >Issue 12: 133 - 137

DOI: https://doi.org/10.11988/ckyyb.20230742

岩土工程

数值流形法渗流分析中处理弱不连续界面的新方法

  • 贾真 , 1 ,
  • 杨冬梅 2 ,
  • 郑宏 1
展开
  • 1 北京工业大学 城市建设学部,北京 100124
  • 2 北京市政建设集团有限责任公司,北京 100048

贾 真(1996-),男,山东菏泽人,博士研究生,主要从事计算岩土力学研究。E-mail:

Copy editor: 罗娟

收稿日期: 2023-07-10

  修回日期: 2023-11-02

  网络出版日期: 2024-01-10

基金资助

国家自然科学基金项目(52130905)

国家自然科学基金项目(52079002)

收起

A New Method to Deal with Weak Discontinuous Interface in Seepage Analysis Based on Numerical Manifold Method

  • JIA Zhen , 1 ,
  • YANG Dong-mei 2 ,
  • ZHENG Hong 1
Expand
  • 1 Faculty of Architecture, Civil and Transportation Engineering, Beijing University of Technology, Beijing100124, China
  • 2 Beijing Municipal Construction Co., Ltd., Beijing 100048, China

Received date: 2023-07-10

  Revised date: 2023-11-02

  Online published: 2024-01-10

Fold

摘要

数值流形法采用数学覆盖和物理覆盖的双覆盖系统,具有灵活处理边界问题、高效进行网格划分以及便捷提高近似精度等优点,是一种很有前景的数值方法。不同于传统数值流形法根据界面来切割数学覆盖形成物理覆盖的做法,数值流形法基于弱不连续物理覆盖,在流行单元中利用折射定律构造出一种新的权函数,以此建立局部近似,并将其应用在稳定渗流问题中。通过对典型算例的计算分析,结果表明该方法在解决不连续界面问题中具有准确性和便利性。

关键词： 数值流形法; 渗流分析; 非均质; 材料界面; 折射定律; 弱不连续

本文引用格式

贾真 , 杨冬梅 , 郑宏 . 数值流形法渗流分析中处理弱不连续界面的新方法[J]. 长江科学院院报, 2024 , 41(12) : 133 -137 . DOI: 10.11988/ckyyb.20230742

Abstract

The numerical manifold method is a promising numerical method which uses a dual coverage system consisting of both mathematical and physical coverages. It has the advantages of flexible handling of boundary problems, efficient meshing and convenient improvement of approximation accuracy. Different from the traditional numerical manifold method which cuts the mathematical coverage to form the physical coverage according to the interface, a new weight function is constructed by using the refraction law in the manifold element to establish the local approximation based on the weak discontinuous physical coverage, and it is applied to the steady seepage problem. The accuracy and convenience of the method in solving discontinuous interface problems is demonstrated through the analysis of the calculation of typical seepage flow problems.

Key words： numerical manifold method; seepage analysis; heterogeneity; material interface; refraction law; weak discontinuous

开放科学(资源服务)标识码(OSID):

0 引言

非均质性是一种普遍存在于各种岩土中的基本性质,通常与裂隙、孔洞、材料界面等特征相关。岩土工程中的水流流动受非均质性的影响很大,因此非均质渗流分析具有重要的研究意义。
传统的实验方法在研究非均质渗流问题时具有较大的局限,因此该问题的研究广泛采用数值方法。其中材料界面具有弱不连续性,经典的有限元方法[1⇓-3]要求网格节点与界面重合以及局部网格细化,这将会使得网格划分复杂化以及全局自由度增加,进而导致计算工作量的大幅增加。为了减少经典有限元网格划分的不便,扩展有限元法(Extended Finite Element Method,XFEM)[4⇓-6]和广义有限元法(Generalized Finite Element Method,GFEM)[7⇓-9]结合单位分解法(Partition of Nuity Method,PUM)[10]对有限元进行推广,网格划分独立于问题域,通过附加函数来扩展近似解,能够高效、高精度地解决复杂实际问题,但这种灵活性的代价是增加了对附加函数数值积分的复杂性,并使方程组的规模和求解难度加大。
数值流形法(Numerical Manifold Method,NMM)[11]在求解不连续问题上具有独特优势,该方法最具创新性的特点之一是它采用了双重覆盖系统,即数学覆盖和物理覆盖,数学网格不受界面所影响,这一优点意味着NMM可以始终使用规则网格来离散问题域。物理覆盖通过界面从数学覆盖中切割得到,同时,NMM的自由度在物理覆盖上定义。通过这种双覆盖网格系统,数值流形法在界面处理上十分灵活高效。An等[12]通过构造特殊的跳跃函数提出了模拟材料不连续的新方法。An等[13]构造出特殊的物理覆盖及局部近似,对任意不连续问题进行了模拟。李伟等[14]使用MLS-NMM技术模拟了稳定渗流中的边界问题。Hu等[15]结合新的拉格朗日乘子法,建立了模拟非均质介质中流体流动的不连续方法。
目前对于非均质渗流模拟中权函数的研究较少,本文采用三角形网格覆盖问题域,引入了一种特殊类型的物理覆盖,并利用界面上水头和法向流量的连续性来构建全局近似值,即折射定律[16]:
h - = h + 　 , q n - = q n + 　 。
式中: h为水头;qn为界面的外法向流量;“-”和“+”分别为界面法向量的负方向和正方向。通过2个典型的渗流算例的计算,验证了此方法的准确性。

1 处理材料界面的新方法

1.1 物理覆盖的生成

图1(a)为两种材料的问题域Ω,由2个子域Ω1和Ω2组成,Γ1-2是Ω1和Ω2的界面。
图1 基于NMM对含有2种材料的问题域建模

Fig.1 Modelling of the problem domain with two materials based on NMM

图1(b)为该问题的传统NMM模型,每个与材料界面相交的数学片(Mathematical Patch, MP) 被分割成两个独立的部分,形成两个物理片 (Physic Patch, PP)。图1(b)中黄色标记的区域为一个数学片,由类似于有限元法中共有同一结点的单元所组成。
图1(c)为相同问题的新NMM模型,采用与图1(b)相同的网格,每个与材料界面相交的MP只形成一个PP,这种PP具有弱不连续性,为区别常规PP,将其称为弱不连续PP。图1(c)中绿色标记的三角形单元是由3个弱不连续PP覆盖形成的流形单元(Manifold Element, ME)。该ME同样具有弱不连续性,将其称为弱不连续ME,在其中采用传统的整体逼近无法满足材料界面处的弱不连续性,为此,需要为弱不连续ME构造特殊的权函数。

1.2 弱不连续ME的权函数

弱不连续ME内的全局近似是与其重叠弱不连续PP的整体逼近,即
h ( x , y ) = ∑ i = 1 3 N i ( x , y ) h i 　 。
式中: h ( x , y )为弱不连续ME内的近似水头; N i ( x , y )为弱不连续物理片 i的权函数,该权函数构建的条件是弱不连续ME中水头连续并且材料界面的流量也是连续的; h i为弱不连续物理片 i的水头。
如图2所示,将一个弱不连续ME分为两部分,分别为三角形145对应的第1部分和四边形2354对应的第2部分,图中 k 1和 k 2为不同的渗透系数。不同于有限元法在同一结点上有一个未知量,本方法在弱不连续ME中同一结点上有两个未知量。在每一部分中,近似值 h j ( x , y )用三角形的形状函数(面积坐标) L 1 、 L 2和 L 3来表示,即
h j ( x , y ) = ∑ i = 1 3 L i ( x , y ) u i j , 　 j ∈ 1,2 　 。
式中: L 1 ( x , y ) 、 L 2 ( x , y ) 、 L 3 ( x , y )分别为弱不连续ME内一点所对应1、2、3点的面积坐标;未知量 u i j与 h i有关。
图2 弱不连续ME

Fig.2 Weak-discontinuous ME

弱不连续ME两部分各自的整体逼近(式(3))共同构成弱不连续ME的整体逼近(式(2))。
根据形状函数的性质,弱不连续流形单元3个节点的水头为
h 1 ( x 1 , y 1 ) = u 1 1 = h 1 　 , h 2 ( x 2 , y 2 ) = u 2 2 = h 2 　 , h 2 ( x 3 , y 3 ) = u 3 2 = h 3 　 。
式中: h 1 ( x 1 , y 1 )为第一部分中1点处的水头; u 1 1为第一部分中1点处的未知量; h 1为弱不连续物理片1的水头;其余符号类同。
基于折射定律,两部分在界面处4、5点的水头相同:
( 1 - L 1 4 ) u 2 1 - L 1 4 u 1 2 = ( 1 - L 1 4 ) h 2 - L 1 4 h 1 　 , ( 1 - L 1 5 ) u 3 1 - L 1 5 u 1 2 = ( 1 - L 1 5 ) h 3 - L 1 5 h 1 　 。 　
式中 L 1 4为节点1的形函数在4点处的值,其余符号类同。
同时,两部分在界面处的流量相同。
k 1 ∂ h 1 ∂ n = k 2 ∂ h 2 ∂ n 　 。
式中: h 1为第一部分中的水头; h 2为第二部分中的水头;n为界面的外法向量。
将式(5)和式(6)联立得
A x = B h 　 。
其中,
x = u 2 1 　 u 3 1 　 u 1 2 T , 　 h = h 1 　 h 2 　 h 3 T 　 , A = 1 - L 1 4 0 - L 1 4 0 1 - L 1 5 - L 1 5 k 1 d 2 k 1 d 3 - k 2 d 1 　 , B = - L 1 4 1 - L 1 4 0 - L 1 5 0 1 - L 1 5 - k 1 d 1 k 2 d 2 k 2 d 3 　 。
矩阵A、B中,d1=c x 3 - x 2+s y 3 - y 2,d2=c x 1 - x 3+s y 1 - y 3,d3=c x 2 - x 1+s y 2 - y 1。
其中 c , s为从4点到5点的单位向量的方向数。
由于矩阵A可逆,式(7)可以表示为变换矩阵C,即
C = A - 1 B 　 。
因此,得到了弱不连续ME两部分所对应的权函数,在三角形145中
h 1 ( x , y ) = N 1 ( x , y ) h 　 。
其中,N1(x,y)= N 1 1 　 N 2 1 　 N 3 1, N 1 1=L1+C11L2+C21L3, N 2 1=C12L2+C22L3, N 3 1=C13L2+C23L3。
式中 C 11为变换矩阵C中第一行第一列所对应的值,其余符号类同。
在四边形2354中
h 2 ( x , y ) = N 2 ( x , y ) h 　 。
式中,N2(x,y)= N 1 2 　 N 2 2 　 N 3 2,其中 N 1 2=C31L1, N 2 2=C32L1+L2, N 3 2=C33L1+L3。
图3(a)、图3(b)分别为弱不连续ME中第1、2部分1点处的权函数示意图。各部分的权函数之和为1,满足单位分解的要求。
图3 弱不连续ME的权函数示意图

注:弱不连续ME被平铺,权函数 N 1 1 　 N 2 1 　 N 3 1中的3个值分别为1、2、3点处垂直虚线的高度,例如 N 1 1为1点处垂直虚线的高度。

Fig.3 Diagram of the weight function of weak-discontinuous ME

值得注意的是,两个相邻的弱不连续ME与材料界面的同一个交点处的水头可能存在不相等的情况,这样的点取其在2个弱不连续ME中的平均值,适当细化网格可以降低误差,从而得到较为精确的解。

2 渗流控制方程及离散

本文将使用该方法求解二维有压渗流问题,假设材料是各向同性的,问题域内不存在“源”或“汇”。
在渗流区域Ω内,总水头 h满足连续性方程
$ -\nabla \cdot(k \nabla h)=0$
水头边界条件Γh为
h = h - 　 。
流量边界条件Γq为
k ∂ h ∂ n = q - 　 。
材料界面条件Γm为
h - = h + 　 , k - ∂ h - ∂ n = k + ∂ h + ∂ n 　 。
式中:$ \nabla$为梯度算子; h -为已知水头边界; q -为已知流量边界; k为渗透系数;“-”和“+”分别为界面Γm外法向量的负方向和正方向。
根据式(11)和边界条件式(12)—式(14)得到泛函为
$ I(h)=\frac{1}{2} \int_{\Omega} k \nabla h \cdot \nabla h \mathrm{~d} \Omega-\int_{\Gamma} \bar{q} h \mathrm{~d} \Gamma$
本文采用罚函数法来确定水头边界条件Γh,因此,式(15)可以进一步表示为
I * ( h ) = I ( h ) + ∫ Γ h 1 2 p ‖ h - h - ‖ 2 d Γ 　 。 　
式中 p为水头边界条件的罚因子。
由 δ I * ( h ) = 0,得到离散控制方程
K h = q 　 。
式中K、h和q分别代表离散控制方程的总渗透系数矩阵、总水头列阵和总流量列阵。
通过式(17)即可得到总水头场(即势函数),若求渗流区域的流网,需求得相应的流函数,具体过程参见文献[17]。

3 积分方案

本文采用三角形背景网格进行数值积分,界面处弱不连续ME拆分出的部分可能为四边形或三角形。若是四边形,则通过Delaunay剖分为2个三角形。利用Hammer积分求得其单元积分。

4 算例

为验证该方法对于稳定渗流分析的正确性,将进行以下算例。

4.1 单一材料界面问题

采用文献[16]中的算例,计算模型的尺寸见图4。图4的左侧上游边界水头为10 m,右侧下游边界水头为2 m,顶部和底部边界为不透水边界,问题域包含2种不同材料,材料界面位于x=1 m处。2个材料域的渗透系数为k1=1×10-8 m/s,k2=1×10-7 m/s。
图4 单一材料界面尺寸

Fig.4 Sizes of single material interface

本例中的材料界面位于弱不连续PP内,无需切割PP,减少了工作量。本文方法计算的水头分布如图5所示,在 y = 1处的水头与解析解几乎一致(如图6)。证实了本文方法求解此问题的准确性。
图5 网格划分和水头分布

Fig.5 Meshing and hydraulic head distribution

图6 y=1 m处本文方法与解析解的对比

Fig.6 Comparison of our solution and analytical solution at profile y=1 m

4.2 多层材料介质问题

采用文献[18]中的算例,此矩形区域有压渗流模型尺寸如图7(a)所示,图中的底边、左边和右上边黑实线所示边界为不透水边界,虚线为材料界面。3层材料的渗透系数为:k1=0.1 cm/s、k2=0.2 cm/s,k3=0.3 cm/s。网格布置如图7(b)所示。
图7 多层材料介质尺寸和网格划分

Fig.7 Sizes of multilayer material medium and meshing

本文方法计算的流网如图8(a)所示,由理论分析,流网满足在材料界面处的变化趋势,即势线和流线发生转折。此外,由于假设每层介质都是各向同性的,所以势线和流线相互垂直,同时势线垂直于不透水边界。本文的计算结果(如图8(a))与文献[18]中的结果(如图8(b))基本一致,并且能够反映出流网的特征。
图8 多层材料介质流网

Fig.8 Flow net of multilayer material medium

图9为本文方法和有限元法在不同横截面处的流量对比,可以看出,本文方法计算出的截面流量几乎与有限元法结果相同。
图9 不同横截面处流量与有限元法结果的对比

Fig.9 Comparison of flow at different cross sections between our solution and FEM

5 结论

本文基于数值流形法的双覆盖网格系统,将材料界面看作内部界面,由物理覆盖所包含,建立了不同于传统NMM的物理片。通过构造出新的权函数,实现了材料界面处水头和法向流量的连续性,满足折射定律。不同于其他方法复杂的权函数形式,本文建立的弱不连续权函数仅利用了折射定律来构造,求解简单且有较高的精度。此外,网格无需与材料界面相匹配,减少网格划分的工作量,实现了高效便捷的计算。2个算例表明了该方法求解稳定渗流问题是可行的,计算结果与解析解或已有的数值结果几乎一致。下一步将对其他单元类型和覆盖开展研究,并利用界面处的切向力和法向力相等来构造新的权函数,扩展到力学分析中。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
崔皓东, 朱岳明, 张家发, 等. 深埋洞室群围岩渗流场分析及渗控效果初步评价[J]. 长江科学院院报, 2009, 26(10): 71-75.

(CUI Hao-dong, ZHU Yue-ming, ZHANG Jia-fa, et al. Analysis on Seepage Field and Preliminary Evaluation on Seepage Control Effect of Rockmass around Deep Buried Caverns[J]. Journal of Yangtze River Scientific Research Institute, 2009, 26(10): 71-75.) (in Chinese)

[2]
张俊霞, 李莉, 朱登峰. 西霞院水库电站厂房坝段三维渗流计算分析[J]. 长江科学院院报, 2009, 26(10):67-70.

(ZHANG Jun-xia, LI Li, ZHU Deng-feng. Three-dimensional Seepage Analysis of Powerhouse Dam Section of Xixia Institute Anti-regulation Reservoir[J]. Journal of Yangtze River Scientific Research Institute, 2009, 26(10):67-70.) (in Chinese)

[3]
张伟, 许季军, 陈劲松. 江口水电站坝址区三维渗流计算分析[J]. 长江科学院院报, 2002, 19(2):31-33.

(ZHANG Wei, XU Ji-jun, CHEN Jin-song. Analysis on Three-dimensional FEM Seepage of Jiangkou Hydroelectric Station’s Dam Area[J]. Journal of Yangtze River Scientific Research Institute, 2002, 19(2):31-33.) (in Chinese)

[4]
BELYTSCHKO T, GRACIE R, VENTURA G. AReview of Extended/Generalized Finite Element Methods for Material Modeling[J]. Modelling and Simulation in Materials Science and Engineering, 2009, 17(4): 043001.

[5]
FRIES T P, BELYTSCHKO T. The Extended/Generalized Finite Element Method: an Overview of the Method and Its Applications[J]. International Journal for Numerical Methods in Engineering, 2010, 84(3): 253-304.

[6]
董玉文, 任青文, 苏琴. 接触摩擦问题的扩展有限元数值模拟方法[J]. 长江科学院院报, 2009, 26(5): 45-49.

(DONG Yu-wen, REN Qing-wen, SU Qin. Extended Finite Element Method of Frictional Contact Problem[J]. Journal of Yangtze River Scientific Research Institute, 2009, 26(5): 45-49.) (in Chinese)

[7]
STROUBOULIS T, BABUŠKA I, COPPS K. The Design and Analysis of the Generalized Finite Element Method[J]. Computer Methods in Applied Mechanics and Engineering, 2000, 181(1/2/3):43-69.

[8]
STROUBOULIS T, COPPS K, BABUŠKA I. The Generalized Finite Element Method[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190(32/33): 4081-4193.

[9]
李录贤, 刘书静, 张慧华, 等. 广义有限元方法研究进展[J]. 应用力学学报, 2009, 26(1):96-108,214.

(LI Lu-xian, LIU Shu-jing, ZHANG Hui-hua, et al. Researching Progress of Generalized Finite Element Method[J]. Chinese Journal of Applied Mechanics, 2009, 26(1):96-108,214.) (in Chinese)

[10]
MELENK J M, BABUŠKA I. Approximation with Harmonic and Generalized Harmonic Polynomials in the Partition of Unity Method[J]. Computer Assisted Mechanics and Engineering Sciences, 1997, 4(3/4): 607-632.

[11]
郑宏. 数值流形法[M]. 北京: 科学出版社, 2022.

(ZHENG Hong. Numerical Manifold Method[M]. Beijing: Science Press, 2022.) (in Chinese)

[12]
AN X, MA G, CAI Y, et al. A New Way to Treat Material Discontinuities in the Numerical Manifold Method[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(47/48): 3296-3308.

[13]
AN X, MA G, CAI Y, et al. Arbitrary Discontinuities in the Numerical Manifold Method[J]. International Journal of Computational Methods, 2011, 8(2): 315-347.

[14]
李伟, 郑宏. 基于数值流形法的渗流问题边界处理新方法[J]. 岩土工程学报, 2017, 39(10): 1867-1873.

(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)

[15]
HU M, WANG Y, RUTQVIST J. Development of a Discontinuous Approach for Modeling Fluid Flow in Heterogeneous Media Using the Numerical Manifold Method[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2015, 39(17): 1932-1952.

[16]
HU M, WANG Y, RUTQVIST J. An Effective Approach for Modeling Fluid Flow in Heterogeneous Media Using Numerical Manifold Method[J]. International Journal for Numerical Methods in Fluids, 2015, 77(8): 459-476.

[17]
AALTO J. Finite Element Seepage Flow Nets[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1984, 8(3): 297-303.

[18]
TRACY F T, RADHAKRISHNAN N. Automatic Generation of Seepage Flow Nets by Finite Element Method[J]. Journal of Computing in Civil Engineering, 1989, 3(3): 268-284.

Options
文章导航

/

〈 〉

*严禁使用本系统存储、处理、传输涉密信息*
版权所有 © 《长江科学院院报》编辑部
技术支持：北京玛格泰克科技发展有限公司