基于ANSYS模拟非稳定饱和_非饱和渗流场_岳庆河

Ss为储水率 , 其值较小 , 一般忽略不计 ;α在饱和区取 1, 非饱和 区取 0;H0 为已知水头 ;q为 边界流 量 ;n为土 体空隙 比 ;cos(n,
x)、cos(n, y)、cos(n, z)为 边 界 面 外 法 线 方 向 余 弦 ;t为 时 间
变量 。
其中式 (1)为连 续性方 程 , 式 (2)～ 式 (4)为边界 条件 ,
图 3 非稳定渗流模拟流程
(下转第 80页 )
· 80·
人 民 黄 河 2008年
tem[ J] .OperRes, 1995(3):187 -197. [ 2] LoucksD P.Somecommentsonlineardecisionrulesand
2.2 渗透参数调整与逸出点位置的确定
渗透 参数是 未知的 , 需 要根据 计算结果 不断调 整 , 以往进 行调整时只使用上一 时步末的水头场分布 , 会产生过 大的累积 误差 , 因此考虑时步初 (即为上一时步末 )与时步末水头场分布 来确定时步内的渗透 参数更为合理 [ 8] 。
逸出 按如 下步骤确定逸出点位 置 [ 9 -10] :①第 1次计算 , 对于初始时步 , 将 可能的逸出面节点 全部作 为未 知点 , 对于 其余 时步 , 将时 步初 的逸出面节点作为 已知 点 , 经 过有 限元 计算求 出渗 流场分 布 。 ②将水面以上的正压 边界点全部按逸出点处 理 , 求出 渗流场分 布后 , 再对水面以上零 压边界 点进 行判断 , 若 流通 量大于 等于 0, 即
第 30卷第 4期 人 民 黄 河 Vol.30, No.4 2008 年 4月 YELLOW RIVER Apr., 2008
【水利水电工程 】
渗流计算是堤 坝设 计及 除险加 固中 的一项 重要 内容 。 以 往的渗流计算 , 一般忽 略非饱 和区 渗流的 影响 , 但 是土的 基质 势和重力势的作用 , 使土 体中孔 隙水 产生非 饱和 流动 [ 1] , 只考
虑饱和区渗流的方法 不能 有效地 解决 水位变 化和 降雨等 引起 的非稳定渗流问题 。
该方法在判断可能逸出面 节点流通量后 , 从逸出 面节点中 成批排除流通量小于 0的节点 , 没有考虑 施加在假定 逸出面上 的固定水 头对节点 流通量的 影响 , 现节合图 2 进行说明 。 图 2 为取水位降落 后部 分坝 坡面 , 曲 线 A为上 一 时步 求得 的 浸润 线 , 曲线 B为此时步待求的浸润线 , 节点 2、6则分别为上一时步 求 得的逸出点及此时步待求的逸出点 。求解时 , 对节点 2 ～ 9施 加 固定水头 (H =z), 节点 2 ～ 6所加固定水头将影响与其相邻 单元的水头值 , 水头与渗透参数的高度非 线性关系可 能导致节 点 2 ～ 7在此时步流通量均流向域内 , 而排除在逸出面之外 , 从 而使逸出点错误地收 敛于节点 8, 浸 润线相 应地 由实际 曲线 B 变为曲线 C。
2.3 非线性迭代的程序结构
非饱和土的渗流 问题是非线性问题 , 饱和度 、渗透系 数 、逸 出面范围等只能在计 算过程中迭代确定 。 迭代分 3层进 行 (见 图 3):最内层迭代确定渗流参数 , 迭代中 边界条 件及时 步长不 变 , 渗流参数采用单元孔隙水压力值所对 应的参数不 断进行调 整 , 直到前后两次计算 各单元 孔压 变化均 小于 某一 定值 ;次外 层迭代确定渗流逸 出面边 界 , 迭代 时步长 不变 , 根 据内层 的计 算结果调整逸出面 边界范 围 , 直至 迭代收 敛 ;最外 层迭代 调整 时间步长 , 只有相邻时 步间的 水头 差在允 许范 围之 内时 , 才进 行下一时步计算 , 否则 此步计 算结 果作废 , 更 改时 间步长 重新 计算 。
1 饱和 -非饱和渗流数学模型及有限元解法
Neuman将非饱和区 渗透 系数 看做 孔隙 水负 压的 函数 , 使
饱和区与非饱和区 渗流控 制方 程得到 统一 [ 5] 。 不 考虑土 骨架
变形及水压缩性时 , 非稳定饱和 -非饱和渗流数学模型如下 :
x[ kx(h) Hx] + y[ ky(h) Hy] + z[ kz(h) Hz] =
逐步积分 , 这些积分只要设定时间步长及 相应的边界 条件与初
始条件后 , 即可由 ANSYS自动实现 。
2 非稳定饱和 -非饱和渗流场数值模拟
2.1 时间步长的自动选择
非稳 定渗流 场数值 模拟中 , 选取 合适的 时间步 长 , 既能保
收稿日期 :2007-10-09 作者简介 :岳庆河 (1981— ), 男 , 山东济南人 , 硕士研究 生 , 主 要研究方向为水工结 构 。 通讯作者 :刘福胜 (1964— ), 教授 。 E-mail:liufsh@
Г 1
为已知水头边界 , Г2 为已知流量边界 , Г3 为逸出段边界 (渗流场
边界见图 1), 式 (5)为初始状态 , H0 即为渗流场的初始水头 。
图 1 渗流场边界条件示意
热传导基本微分 方程为
kx
2xT2 +ky
2yT2 +kz 2zT2 =Cρ
T t
(6)
比较渗流微分方 程与热传导微分方程 , 可见二者 控制方程
已有不少学者 应用 ANSYS进 行了 渗流 场数 值模 拟研 究 。 许玉景等 [ 2] 首次使用 ANSYS热分析模块 结合单元生 死技术进 行了稳定渗流模拟 , 蒙富强 [ 3] 对此 进行了 进一 步充 实 , 罗 启北 等 [ 4] 使用 ANSYS进行了导流洞有压渗 流的动态 模拟 。 为解决 存在浸润面的 非稳定 渗流 场模 拟问题 , 笔 者对 ANSYS进 行二 次开发 , 编制了非稳定饱 和 -非饱 和渗流 程序 , 并 通过实 例验 证了程序的正确性 。
(4)
H(x, y, z, t0) =H0(x, y, z, t0 )
(5)
式中 :H =h+z, H为总水头 ;h为压力水头 , 在饱和区为 正值 ,
非 饱和区为负值 ;kx、ky、kz分别为沿坐标 x、y、z的孔隙水主渗透 系 数 , 是压力水头 h的函数 ;C= θh为容水度 , θ为介质含水率 ;
余时步时 , 将上一时步 确定的 逸出 面节点 作为 已知 点 , 计 算出 渗流场分布 。 ②不 改变 边界条 件 , 对渗 透参数 进行 循环迭 代 , 直到达到足够精度 后 , 再判断 逸出 点流通 量 , 若流 通量大 于等 于零 , 则该点为 逸出 点 ;反 之 , 在 可能 逸出 面 节点 中删 除 此节 点 。 重复这 2个步骤 , 直到找出流通量大 于等于 0的 可能逸出 点为止 。
基于 ANSYS模拟非稳定饱和 -非饱和渗流场
岳庆河 1 , 刘福胜1 , 刘佩玺2
(1.山东农业大学 水利土木工程学院 , 山东 泰安 271018;2.济宁市水利局 , 山东 济宁 272019)
摘 要 :堤坝非稳定渗流数值模拟中 , 浸润面 、逸出点事先未知 , 需要迭代计算 , 渗透参数需根据 计算结果不断调整 , 都增 加了数值模拟的难度 。 选用 Neuman饱和 -非饱和渗 流数学 模型进行 非稳定 渗流场 数值模 拟 , 分析了数 值模拟 中逸出 面一般求解方法存在的 问题并加以改进 , 使用 ANSYS开发出渗流数值模拟程序 , 经实例验证表明该法是实用有效的 。 关 键 词 :饱和 -非饱和 ;渗流 ;逸出面 ;ANSYS 中图分类号 :TV139.14 文献标识码 :A 文章编号 :1000-1379(2008)04-0076-02
-kx Hxcos(n, x)-ky Hycos(n, y)-kz Hzcos(n, z)≥ 0 说明该节点流向域外 , 则该点为逸出面节 点 , 否 则 , 该 点不是逸 出面节点 。 ③将前 一次 确 定的 逸 出 面节 点 作为 已 知点 处 理 。 重复这 3个步骤 , 直到计算收敛 , 再进行下一时步 计算 。
的形式完全相 同 , 这就为 应用 ANSYS热分 析模 块进 行渗 流场 分析奠定了基础 [ 2 -4] 。 对二维渗流微分 方程应用 变分法 , 并对
时间取隐式差分 , 可得相应的有限元格式 。 对于存在 时 间变 量 的瞬 态 热 分析 问 题的 求 解 [ 6] , ANSYS
在空间域上按有限元 方法求逐步积分 , 在 时间域上用 差分方法
[ C(h)+αSs]
H t
(1)
H(x, y, z, t) =H1 (x, y, z, t), 在 Г1 上
(2)
kx Hxcos(n, x)+ky Hycos(n, y)+kz Hzcos(n, z)=q, 在 Г2 上
(3)
H(x, y, z, t) =z(x, y, z, t), t> 0, 在 Г3 上
图 2 逸出面调整示意 笔者对此加以改 近 , 只判 断可能 逸出 点 (也就 是逸出 面上 最高位置节点 )的流通量 , 如果流向域外则 计算收敛 , 否则从可 能逸出面节点中删除此节点 , 进行下一 次循环 。 这样 在每次循 环中只排除一 个可 能 逸出 面节 点 , 避 免了 上 述错 误情 况 的发 生 。 具体 做法为 :①第 1次 计算 , 对于初 始时步 将可能 的逸出 面节点全部作为已 知水头 节点 (H =z), 求 出渗流 场分布 , 其
chanceconstraints[ J] .WaterResour, 1970, 6(2):668 671. [ 3] 赵鸣雁 , 程春田 , 李刚 .水 库群系 统优化 调度新进 展 [ J] . 水文 , 2005, 25(6):18 -23. [ 4] AhmedI.Onthedeterminationofmultireservoiroperation policyunderuncertainty[ D] .TucsonArizona:TheuniversityofArizona, 2001. [ 5] 谭维炎 , 黄守信 .应用随机动态规 划进行水电站水库的优 化调度 [ J] .水利学报 , 1982(7):1 -7. [ 6] 李爱玲 .梯级水电 站水库 群兴利随 机优化 调度数 学模型 与方法研究 [ J] .水利学报 , 1998(5):71 -74. [ 7] DorfmanR.Themulti-structureapproachindesignofwaterresourcesystems[ M] .HarvardUniversityPress, 1962. [ 8] WindsorJS.Optimizationmodelforreservoirfloodcontrol [ J] .WaterResour, 1973, 9(5):1103 -1114. [ 9] BarrosM, TsaiF, YangS, etal.OptimizationofLargescalehydropowersystem operations[ J] .JournalofWater Resour, 2003, 129(3):178 -188. [ 10] PengCS, BurasN.Practicalestimationofinflowintomultireservoirsystem[ J] .JournalofWaterResour, 2000, 126 (5):35 -40. [ 11] 黄强 .黄河干流水库联合实施调 度及智能决策支持系统 研究 [ D] .西安 :西安理 工大学 , 1995.
第 4期 岳庆河等 :基于 ANSYS模拟非稳定饱和 -非饱和渗流场
· 77·
证迭代收敛 , 又可避免 步长太 小而 引起累 积误 差过 大 , 且 节省 机时 。 当水位变化剧烈时 , 应选用较小的 时间步长以 提高计算 精度 ;水位变化缓慢时 , 则可用 较大的时间步长 节省机时 , 相邻 时步的步长变更不应 太大 [ 7] 。 选取参数 Δh记录相 邻时步的最 大单元水头差 , 设定水头差值的 允许范 围为 [ hmin, hmax] 。 Δh 在允许范围内 时继续 下一 时步 计算 , Δh>hmax时则 将时 间步 长减半 ;Δh<hmin时则将时间步长 增大 1.5倍 。 这样 , 在计算 过程中对时间步长就 进行了自动调整 。