ANSYS软件在平原水库土工膜气胀计算中的应用

　　国内外对于土工膜膜下气体传递和聚集的研究多集中在垃圾填埋场的防渗体系及其排气状态的稳态分析，对于平原水库土工膜的气胀探索尚不充分。ANSYS软件发展至今已相当完善，可以进行结构、热学、电磁场以及各种耦合场的稳态和瞬态有限元分析，其强大的前、后处理功能可以方便地表述模型各点状态沿时间变化的过程，并且利用物理场相似的原理已经在土石坝渗流中取得了应用。本文将从非饱和土气场和温度场的可比拟性出发，用ANSYS热分析模块解决膜下气场的瞬态数值模拟。

　　1. 非饱和土气场与温度场的相似

　　非饱和土中气场与热传导问题的温度场的相似性存在于基本理论、微分方程以及边界条件等三方面。

　　1.1 基本理论相似性

　　热传导问题的傅里叶定律为：

　　(1) 式中为单位时间内单位面积的热流量;为导热系数;为温度;为方向上的温度梯度。

　　非饱和土中的气体扩散遵从的Fick定律的表达形式为：

　　(2) 式中为通过单位面积土的空气流量;为土中空气流动的传导系数;为空气浓度;为方向上的浓度梯度。

　　由式(1)、(2)可知，非饱和土中气体传导和温度场中热传导的驱动势能都是浓度梯度，在数值计算中具有相同的表达形式。

　　1.2 微分方程相似性

　　三维热传导方程的表达形式为：

　　(3) 式中为介质比热容;为密度;为导热系数;为热源强度。

　　在介质各向同性和不考虑热源的情况下(3)式可以简化为：

　　(4) 非饱和土中非稳态气流的气象偏微分方程如下：

　　(5) 式中为气体孔隙压力;为孔隙水压力;为时间;为气象固结系数，与传导系数有关;为常数;为传导系数。在各向同性问题中，传导系数不沿程变化，则、和为零。只取地下水位之上的土体作为研究对象，则可假定此部分的孔隙水压力为零，亦为零，则(5)式化简为

　　(6) 与各向同性的热传导方程(4)一致。

　　由式(4)、(6)可知，非饱和土中气体传导和温度场中热传导可以用相同的微分方程描述，因此在数值计算中属于同一类问题。

　　1.3 边界条件相似性

　　对于温度场

　　①初始条件：，即初始时刻温度场各点的温度分布;

　　②第一类边界条件：，边界上点在时刻的温度;

　　③第二类边界条件：，边界上点在时刻热流强度函数，特别地，对于绝热边界。

　　对于非饱和土气场

　　①初始条件：，即初始时刻非饱和土气场各点的孔隙气压分布;

　　②第一类边界条件：，边界上点在时刻的孔隙气压;

　　③第二类边界条件：，边界上点在时刻气流函数，特别地，对于不透气边界。

　　由温度场和非饱和土气场的边界条件可知，两种物理场的边界条件在本质上也是相等的。因此在利用温度场分析非饱和土气体问题时，可以将各种参数对应输入进行求解。

　　2. 应用ANSYS计算非饱和气场

　　2.1 前处理

　　ANSYS中热分析模块专为解决热传导问题而设置，在计算非饱和土气场问题时还需做一定的转化。在前处理中，热分析模块要求输入材料的导热系数(Conductivity)、密度(Density)以及比热容(Specific Heat)。比热容的物理意义为单位质量物质升高单位温度所需热量，其表达式为

　　(7) 因此在输入材料的比热容时，需要考虑的是单位质量的土体在气体输入时孔隙气压升高的能力，由理想气体状态方程

　　(8) 可知，在体积不变的情况下土体的孔隙气压与体积内空气的物质的量成正比，而该方程也可以写成类似于比热容公式的表达形式

　　(9) 孔隙气压在计算初始状态下一般为标准大气压101.325kPa。在不考虑耦合场的情况下，为简便起见，初始状态下单位体积内空气的物质的量取亦取101.325，则土体“比气容“可取1，其余各项赋值1。在ANSYS定义材料性质、输入边界和荷载条件时，各项参数的对比和输入如下

　　表1. ANSYS分析非饱和土气场输入的对应参数

　　在建模和划分网格阶段，需利用对称原理将原本庞大的计算模型分割。在处理对称图形时可以将模型分割成1/2或1/4，在前处理中可使模型细节所占比例扩大和方便网格细化，在后处理中也大大节约了计算时间，同时对于一些关键节点的状态追踪也更为直观。

　　2.2 后处理

　　ANSYS软件包含了强大的后处理功能，为研究非饱和土气体运动的规律提供了强有力的技术支持。对于平原水库气胀现象的研究，在排气设施布置、地下水位上升速率和土体性质各不相同的情况下利用ANSYS后处理可以从不同的模型计算实例中找出各方面影响因素的重点，为之后的工程设计提供理论依据。对于非稳态分析，可以通过对某一特定节点或单元的追踪了解一个特定模型可能出现的最不利状况及其发生位置，并对提出的防治措施提出方案比较。

　　3. 计算实例

　　某水库地下水位为土工膜下5.0m处，经过2个月上升至距膜下1.3m处，土的渗气系数为0.2m/h，排气盲沟为矩形网格排列，相距20m，膜上覆土压重30kPa，计算简图如图1所示。现利用ANSYS软件分析此问题。

　　图1 计算模型的简图和有限元网格划分

　　①-排气盲沟 ②-非饱和土区域 ③-饱和土区域

　　图1为本实例的简图和其有限元计算网格。尽管排气盲沟在整个模型中所占的比例很小，但是在模型建立时仍不能简单地将其视为直线，否则造成的误差将相当大

　　图2 计算结果的孔隙气压力分布云图

　　图2为该模型1/4计算结果最终时刻孔隙气压力的分布云图。从本图中可以观察到，随着时间的推移，中心区域膜下非饱和土气体无法及时消散而发生累积。因此最易发生气胀现象的部位为膜下区域的中心。

　　图3 膜下气场中心点孔隙气压力变化

　　图3是计算模型中心点孔隙气压力随时间的变化图。本例中膜下中心区域的最大孔隙气压力为126kPa(标准大气压)，小于膜上覆土压重。

　　4. 结论

　　本文分析了非饱和土中气场和热传导问题温度场的基本理论、微分方程和边界条件，证明了应用ANSYS热分析模块进行非饱和土气胀问题的可行性，并为今后进一步研究非饱和土气体提供了理论基础和数值模拟方法。该方法建模简单，计算时间短，并且通过本文算例的计算结果与实验数据吻合良好。

　　在前处理中，改变材料参数可以研究不同因素对于非饱和土气体运行的影响;在后处理中，通过ANSYS的POST1和POST26模块可以方便地观察计算模型的整体孔隙气压力云图分布以及某一点的孔隙气压力随时间的变化和趋势，可以为工程设计和修改提供计算依据，具有普遍的应用意义。

　　参考文献

　　[1] 杜晖. 平原水库土工膜防渗技术研究与应用.[硕士论文].河海大学. 2006年2月。

　　[2] T. D. Stark & H. Choi. Methane gas migration through geomembranes. Geosynthetics International, 2005, 12, No. 1

　　[3] A. Bouazzaa, T. Vangpaisalb, H. Abuel-Nagaa, J. Kodikara. Analytical modelling of gas leakage rate through a geosynthetic clay liner–geomembrane composite liner due to a circular defect in the geomembrane. Geotextiles and Geomembranes 26 .2008.

　　[4] 贺晓明，李智录，王淑贤. 利用ANSYS热分析模块分析渗流场问题的探究. 工程地质计算机应用[J]. 2005年第4期.

　　[5] 刘玉鑫. 大学物理通用教程·热学[M]. 北京大学出版社. 2001.2

　　[6] Delwyn G. Fredlund & Harianto Rahardjo. Soil Mechanics for Unsaturated Soils [M]. 中国建筑工业出版社. 北京 1997.5

　　[7] 梁昆淼. 数学物理方法(第四版)[M]. 高等教育出版社. 2009.8