粘土心墙坝渗流场分析

第3期YUNNANWATERPOWER

37

粘土心墙坝渗流场分析

刘永豪,涂兴怀,李飞燕

(西华大学能源与环境学院,四川成都　610039)

摘　要:通过渗流场和温度场的基本微分方程及边界条件的比较分析,将ANSYS软件的热分析模块应用于渗流场的分析,并采用生死单元技术,通过迭代算法计算自由水面位置(浸润线),解决了某实际工程粘土心墙土石坝渗流稳定问题的求解。该方法可以解决复杂边界、多种介质的渗流问题,同时为其他实际工程设计应用提供了强有力的手段。关键词:ANSYS软件;生死单元;渗流场;浸润线

中图分类号:TV222.2　　　文献标识码:A　　　文章编号:1006-3951(2008)03-0037-04

TheAnalysisofClayCoreEarth-RockFillDamSeepageFlowField

LIUYong-hao,TUXing-Huai,LIFei-yan

(SchoolofEnergyandEnvironment,XihuaUniversity,Chengdu610039,SichuanProvince,China)

Abstract:Accordingtothecomparisonandanalysisoftheseepageflowfieldandtemperaturefield'sfundamentaldiffer-entialequationsandboundaryconditionstheANSYSsoftware'sthermalanalysismoduleisappliedtotheanalysisofseep-ageflowfieldandtheelementbirth-deathtechniqueisadoptedtosolvetheseepageflowstabilityproblemoftheclaycoreearth-rockfilldamforaspecificprojectbycalculatingthefreewatersurfaceposition(saturationline)usingtheite-rativecomputingmethod.Thismethodcansolvetheseepageflowproblemwithcomplexboundariesandmediums.

Keywords:ANSYSsoftware;elementbirth-death;seepageflowfield;saturationline

1　前言

从20世纪初开始,渗流对工程的影响已为工程界广泛重视,许多工程技术人员及学者从工程实践和理论两方面进行了大量研究,并取得了许多有价值的成果,既解决了工程中的实际问题,又丰富和发展了渗流理论。随着相关学科的不断发展和完善,各类工程实践提出的渗流问题日益广泛和复杂。伴随着实验方法的逐渐改进,以及计算机在渗流计算中的普遍应用,渗流逐步形成具有自己的理论、研究方法和应用范围的独立学科。目前,主要的研究方法有流网法、差分法、有限元法、水力学法和模型试验。随着有限元技术的成熟,有限元法成为渗流分析的主要数值方法,对渗流场可以达到比较精确的模拟。ANSYS软件是美国ANSYS公司的产品,是基于有限元技术的大型通用分析仿真平台,广泛应用到各个行业领域,具有强大的前、后处理功能和求解器。在ANSYS软件的热分析计算模块中,定义有“渗流材料”,具有专门的渗流场模拟计算功能,能得

出令人满意的结果。

2　计算原理

2.1　微分方程的相似性

对于不可压缩各项异性非均质无源稳定渗流微分方程为:

H H H(kx)+(ky)+(kz)=0, x x y y z z式中:kx,ky,kz—x,、y、z方向的渗透系数

对于无热源的各项异性非均质稳定热传导微分方程为:

T)T)T)(kx + (ky + (kz =0,

x x y y z z式中:kx,ky,kz—x,、y、z方向的热传导系数2.2　初始条件和边界条件的相似性

渗流场的初始条件:H t=0=f0(x,y,z,0)在Ψ内;

热传导温度场的初始条件:T T0(x,y,z,t=0=0)在Ψ内。

收稿日期:2008-04-03

　作者简介:刘永豪(1984-),男,山东菏泽人,硕士研究生,研究方向:水工结构分析及优化技术。　　　　　　　　　　　　　　　　　　　　　　　　　　38云南水力发电　　　　　　　　　　　　　　　　　　　　　2008年第3期

第一类边界条件:

渗流场:H x,y,z,t)在Γ温度Γ=f1(1上;1

场:T x,y,z,t)在ΓΓ=T1(1在上。1

第二类边界条件:

H

渗流场:kn 2=f2(x,y,z,t)在Γ温2上; nΓ

T =f(度场:kn 在Γ2x,y,z,t)2上。2 nΓ

式中:Ψ—渗流和温度场区域,即为Γ1和Γ2所围成的研究区域;H=H(x,y,z,t)、T=T(x,y,z,t)为域内各点水头值和温度值,为待求函数;f0(x,y,z,0)、T0(x,y,z,0)—域内各点上的初始水头值和温度值;f1(x,y,z,t)、T1(x,y,z,t)—Γ1上的已知水头值和温度值;f2(x,y,z,t)、T2(x,y,z,t)—Γ2上的已知流量;n—Γ2的法线方向;Γ1—已知水头值和温度值的边界曲面;Γ2—已知流量值和热流量值的边界曲面。

通过以上分析可知,土体渗流问题是温度场问题的一种特殊形式。只需将温度场介质换成土体介质,热传导率换成土体渗透系数,温度换成渗流总水头,同时取域内热源密度为零,边界条件相应地变为已知水头分布和渗透流速分布,就可以用ANSYS软件进行渗流场的计算。

死功能,“杀死”该单元,找出所有位于自由面以上的单元,将其“杀死”;

3)重新求出各节点的水头值H,检验是否满足 H-H ≤ε,ε趋向于0,若满足,结束计算,如不满足,则重新施加边界条件,重新分析渗流,转入第2)步,直到满足精度要求为止。

逸出点的求解和渗流自由面上的结点一样,作为未知结点,通过计算来确定其位置,但是修改调整逸出点的位置是沿着溢出面(即坝坡面)进行的。

由计算步骤来看,分析将需要反复进行大量重复操作,包括验证单元是否高于浸润线、激活或杀死单元、施加相应边界条件等,若全部采用人机交互进行,工作量很大。本文采用ANSYS提供的APDL语言将模型参数化,使问题得以简化。

在迭代计算中,有可能需要激活已杀死的单元。对于自由面穿过的单元即过渡单元,该方法并没有给予解决。由于过渡单元的复合缝面处渗透系数发生突变,解的稳定性受到影响,这个缺陷可以通过或者采用ANSYS的网格自适应功能在浸润线附近以及水头变化较大区域进行局部网格加密来解决。这样就可以得到足够精度和足够光滑的浸润线。

n+1

n

n

4　应用算例

按照以上的思路和方法,运用ANSYS软件中的温度场,对如图1所示的均质矩形土坝稳定渗流问题进行计算,上游水头为6m,下游水头1m,坝高6m,坝宽4m,渗透系数1m d。单元采用平面三角形单元形式划分,如图1。计算后的渗流场如图2～4。图中均清晰地显示出了求出的渗流自由面位置。表1中列举了渗流自由面上典型点位置和通过甘油模型试验得出的渗流自由面典型位置两者间的比较。从表中看,两者的计算结果基本一致,最大误差只有1.4%,说明了该方法的准确和可靠性。

表1　典型点位置水头对比单位表

典型点ABCD

试验值 m5.625.084.363.23

计算值 m5.6315.1124.3683.278

误差 %0.190.620.181.4

3　应用分析方法

土石坝渗流问题实际上是存在自由面的无压渗

流问题。自由面的位置需要反复的迭代来确定,故属于边界非线性问题。采用有限单元法求解无压渗流场时,通常有变网格法和固定网格法两种。而文中利用的ANSYS中的生死单元技术就是固定网格法的一种。采用固定网格法求解渗流自由面,目的是通过消除土坝非饱和区单元对总体刚度矩阵的贡献,从而达到精确求解。ANSYS提供单元的生死,就可以消除某个单元对总体刚度矩阵的贡献。“杀死”一个单元,并不是从模型中删除该单元,而是在单元刚度矩阵前乘以一个很小的因子,使该单元对总体刚度矩阵没有贡献。这和Bathe提出的单元渗透矩阵调整法在原理上是一致的,因此利用ANSYS的生死单元计算渗流自由面是可行的。其具体步骤是:

1)首先假设整个渗流区域的渗透系数都是饱和渗透系数,求解出各节点的水头H;

2)比较各单元节点的Y坐标值与水头H,如果n

H大于Y,说明该点位于自由面以上,利用单元生n

5　工程实例

5.1　工程概况

某粘土心墙土石坝,坝顶高程为327.50m,设计　　　　　　　　　　　　　　　　　　　刘永豪,涂兴怀,李飞燕　粘土心墙坝渗流场分析39

最大坝高48.0m,坝顶宽6.0m,坝顶长311.0m。大坝桩号0+292断面是主断面,其坝体形状见图5。

图1　网格划分图　　　　　　　　　　　　　　　　　图2　自由面位置图

图3　坝体内总水头等势线图　　　　　　　　　　　　　图4　坝体内流速矢量图

设计坝坡根据地质提供的材料参数确定:上游坝坡高程327.5～312.5m之间为1∶1.8,312.5～297.5m之间为1∶1.9,297.5m以下为1∶2.0;下游坝坡高程327.5～312.5m之间为1∶1.75,312.5～297.5m之间为1∶1.8,297.5m以下为1∶1.9,上下游坝坡在高程312.5m、297.5m上均设置2.0m宽的马道。粘土心墙顶高程325.7m,坡度为1∶0.3。反滤层分为两层,第一层为粒径小于5mm的砂反滤层,第二层为粒径小于80mm的混合砂砾料反滤层。下游的砂反滤层水平宽度为2m,混合反滤层水平宽度为3m;上游的砂反滤层水平宽度为1m,混合反滤层水平宽度为2m。5.2　计算参数和条件

计算工况为上游水位321.5m,下游水位289.5

-3

m;上下游坝壳为砂质砾,渗透系数为2.5×10m

-5

s;上下游反渗层为粉砂,渗透系数为1.0×10m s;心墙为粉质粘土,渗透系数为1.0×10m s。坝基

-3

渗透系数是2.0×10m s。计算的范围为:以上游坝踵为起点,向上游取20m,以下游坝脚为起点,向-7

下游取20m,坝基取至强风化花岗岩层,深度为4m。5.3　有限元分析

建立模型,划分网格。如图5;计算结果如图6、7、8。

如图6、7、8所示,分别给出了浸润线、等势线和渗透流速分布图,采用ANSYS提供的彩色云图及矢量方式表示,可以更加直观地观察渗流场情况,也可以通过列表给出确切计算结果。在图7和图8中,利用生死单元技术将非饱和区域的单元,也即是压力小于0的单元杀死,仅留压力大于0的饱和区域,这样就在二维渗流面上出现一条压力等于0的线,也即是水头值等于坐标值的线———浸润线。通过图7可以看出,坝体内的水头主要在心墙内降幅比较大,说明渗透系数较小的心墙起到了主要防渗作用。通过图8可以看出,上游坝面和心墙内的流速都比下游坝面小,最大流速出现在心墙底端和下游坝面与下游水位交界点处,这与实际是相符的。从图6、7、8中可以看出,用ANSYS温度场模块计算出渗流　　　　　　　　　　　　　　　　　　　　　　　　　　40云南水力发电　　　　　　　　　　　　　　　　　　　　　2008年第3期

场的浸润线、等势线、流速分布,总体上是符合一般心墙土石坝渗流场规律的。通过具体的列表数值,

如坡降、水头、流速、渗流流量等,用于为工程设计服务,是令人满意的。

图5　某粘土心墙土石坝剖面网格图　　　　　　　　　　　　　　图6　渗流自由面位置图

图7　坝体内总水头等势线图　　　　　　　　　　　　　　图8　坝体内流速矢量图

6　结论

参考文献:

通过以上事例说明,采用ANSYS热模块分析土坝渗流场渗流计算是切实可行的,生死单元技术能够解决浸润线位置计算问题,通过浸润线的确定可以用于土坝结构计算,计算分析土石坝的防渗设计和稳定安全情况,从而为工程应用提供了有效的解决方案。利用ANSYS强大的前后处理功能能够为工程应用提供极大的便利,能够在解决复杂边界、多种介质的复杂渗流问题中取得令人满意的效果。强大的后处理功能可以处理复杂的数据。ANSYS为土石坝的渗流分析提供了强有力的手段。

〔1〕　蒋玉川,张建海.弹性力学与有限单元法〔M〕.北京:科学出版

社,2006.

〔2〕　顾慰慈.渗流控制力理论及应用〔M〕.北京:中国建材工业出版

社,2000.

〔3〕　郝文化.ANSYS土木工程应用实例〔M〕.北京:中国水利水电出

版社,2005.

〔4〕　博弈创作室.APDL参数化有限元分析技术及其应用实例〔M〕.

北京:中国水利水电出版社,2004.

〔5〕　贺晓明.基于ANSYS的大坝渗流分析研究〔D〕.2006.

〔6〕　蒙富强.基于ANSYS的土石坝稳定渗流场的数值模拟〔D〕.

2005.

一方有难　八方支援　危难时刻见真情

　　四川汶川发生8.0级特大地震后,云南全体水电员工情系灾区,守望相助,踊跃捐款捐物,奉献爱心。据不完全统计,至5月底云南电网公司员工捐

款1315.10万元;中水顾问集团昆明勘测设计研究院捐款141.13万元;水电十四局首次捐款100万

元;云南华能澜沧江水电有限公司捐款300万元,……。目前云南全体水电员工的爱心行动仍在继续,爱心有如插上翅膀,飞向地震灾区。大家纷纷表

示,汶川地震,感同身受,华夏同心,我们有责任为抗震救灾尽己所能,为灾区群众敬献绵薄之力。