Advances in Porous Flow 渗流力学进展, 2013, 3, 24-28 doi:10.12677/apf.2013.31004 Published Online March 2013 (http://www.hanspub.org/journal/apf.html) Copyright © 2013 Hanspub 24 Numerical Method and Figures Based on the Network Structure Model of Two-Phase Seepage Flow* Lijiang Zhou, Qingyu Li, Qian Wang, Xianzhi Xu# University of Science and Technology of China, Hefei Email: #xuxz@ustc.edu.cn Received: Dec. 6th, 2012; revised: Dec. 14th, 2012; accepted: Dec. 26th, 2012 Abstract: This paper investigates a method to describing the porous medium, which is simplifying the pore part of po- rous medium into network structure model consist of node and pipeline. We simulate seepage flow of water and air in porous medium through this model, obtain the figures of gas-liquid distribution when steady, and cipher out the prob- ability of that water flows through the whole medium. Keywords: Two-Phase Seepage Flow; Node; Pipeline; Network Structure; Hydrophilic; Hydrophobic; Figures of Gas-Liquid Distribution 两相渗流的拓扑网络模型算法和图形* 周李姜,李清宇,王 乾,徐献芝# 中国科学技术大学,合肥 Email: #xuxz@ustc.edu.cn 收稿日期:2012 年12 月6日;修回日期:2012年12月14 日;录用日期:2012 年12月26 日 摘 要:本文探究一种描述多孔介质的方法,将多孔介质的孔隙部分简化为由节点和管道组成的拓扑网络模型, 模拟水和空气两相在介质中的渗流过程,得到流动稳定时介质中两相的分布图形,并计算出水能够透过整个介 质的概率。 关键词:两相渗流;节点;管道;拓扑网络;亲水;憎水;气液分布图形 1. 引言 关于多孔介质中渗流过程的研究,前人已做过许 多工作[1-6]。本文试图提出一种简明的描述多孔介质的 模型,利用这个模型模拟二维的两相渗流过程。 多孔介质包含孔隙部分和固体骨架部分,由于流 体只会在孔隙部分流动,所以模型只研究孔隙部分, 而将固体骨架部分排除在模型之外。 多孔介质可能由多种材料构成,不同材料的润湿 性也会不同。模型只研究整个多孔介质中的两相分 布,而不涉及介质某个孔隙中的两相分布,故对于润 湿性的不同本文做了合理的简化。 本文在多孔介质的简化模型中,对两相渗流进行 计算,得到了不同材料构成的多孔介质的两相渗流图 形,并对图形进行了分析。 2. 一种简单的两相渗流物理模型 如图 1所示,多孔介质是均匀、各向同性的,介 质上方是空气,下方是水。水和空气的压力相同。初 始状态下介质的孔隙中全部是空气,介质下方的水由 于表面张力作用通过介质向上渗透,同时水也可以在 介质中向各个方向流动(包括向下流动)。 *基金项目:国家自然科学基金,项目批准号 10872193,2009.1- 2011.1。 #通讯作者。 两相渗流的拓扑网络模型算法和图形 Copyright © 2013 Hanspub 25 空气 水 Figure 1. The two phase flow of physical model 图1. 两相渗流的物理模型 液气两相渗流的基本方程有(不考虑空气在水中 的溶解): w ww K vgradP (液相运动方程) o oo K vgradP (气相运动方程) ww ww div vs t (液相连续方程) 1 ggg w div vs t ( 气相连续方程) 对于两相渗流问题,可以利用基本方程求解,但 求解过程比较复杂,并且本模型中不存在压力,运动 方程不完全适用。因此我们试图使用另一种方法进行 计算,得到渗流稳定时介质中的气液分布。 3. 简化模型 如图 2所示,我们忽略掉介质中的固体骨架部分, 用节点堆砌成介质,每个节点通过管道与它周围节点 相连,那么介质的孔隙就简化为了节点和管道,流体 沿各个方向的管道流动。 如果介质由多种材料构成,体现在网络模型中就 是水在不同管道中的润湿性可能不同。在保证计算结 果准确性的基础上,为使计算过程更简单,我们将所 有管道分为两类:极端亲水管道和极端憎水管道。具 体来说,就是将水能完全通过的管道一律视为极端亲 水的,管道中水的流动速度没有差别;将水不能完全 通过的管道一律视为极端憎水的,管道中没有水的存 在。两种管道分别均匀分布在整个模型中。 4. 基于拓扑网络模型的渗流过程算法 我们设置了不同的亲水和憎水管道所占比例,建 立不同比例下的拓扑网络模型(网络尺寸为 60 × 60)。 如图 3所示,初始状态下的网络图形如下(亲水、憎水 管道各占 50%,蓝色管道表示亲水的,红色管道表示 憎水的): 算法一:欧拉算法,基于坐标进行计算。底部节 点全部赋值为有水,从坐标底部向上一层一层扫描, 再从坐标顶部向下扫描,反复进行该过程,直到扫描 的结果不再变化为止,最后得到气液分布图。 算法二:拉格朗日算法,基于质点进行计算。选 取网格最下层第一个节点作为出发的水头,水头每到 一个节点,都可能流向它周围所有管道,直到水头所 到节点周围所有管道都有水或憎水为止。然后选取最 下层第二个节点,重复上述计算,直到底部所有节点 计算完成为止,最后得到气液分布图。 我们分别使用上述两种算法进行计算。虽然得出 的结果完全相同,但与欧拉算法相比,拉格朗日算法 略为简洁一些,而且计算过程更接近实际物理情形。 在以后的计算中,我们全部使用拉格朗日算法。该算 法是由 Matlab 软件编程实现的,程序代码见附录。 Figure 2. Porous medium topology network model 图2. 多孔介质的拓扑网络模型 010 2030 4050 60 0 10 20 30 40 50 60 x y Figure 3. The initial state of the topology network model 图3. 初始状态下的拓扑网络模型 两相渗流的拓扑网络模型算法和图形 Copyright © 2013 Hanspub 26 5. 不同管道比例下的气液分布图 我们基于拓扑网络模型,使用拉格朗日算法对两 相渗流过程进行计算,得到不同管道比例下的气液分 布图。一些典型的图形(如图 4~9)可见,(蓝色管道中 为水,红色管道中为空气): 我们计算了水渗出介质的概率,取的亲水管道比 例区间为 40%~60%,概率曲线图如如图 10: 6. 分析讨论 1) 在绘制拓扑网络气液分布图形的过程中我们 发现,虽然各次绘制的图形细节上有所不同,但整个 气液分布图仍然有一定规律可循。当亲水管道比例小 于40%时,水基本只分布在网格的底层;当亲水管到 比例大于 60%时,水基本能够均匀分布在整个网格 010 2030 40 50 60 0 10 20 30 40 50 60 x y Figure 4. Close water pipeline accounted for 20%, hydrophobic pipeline 80% of the gas liquid profile 图4. 亲水管道占20%、憎水管道占80%的气液分布图 010 20 304050 60 0 10 20 30 40 50 60 x y Figure 5. Close water pipeline accounted for 40%, hydrophobic pipeline 60% of the gas liquid profile 图5. 亲水管道占40%、憎水管道占60%的气液分布图 010 20 30 40 5060 0 10 20 30 40 50 60 x y Figure 6. Close water pipeline accounted for 50%, hydrophobic pipeline 50% of the gas liquid profile (water can seep medium) 图6. 亲水管道占50%、憎水管道占50%的气液分布图(水能够渗出 介质) 010 20 30 40 5060 0 10 20 30 40 50 60 x y Figure 7. Close water pipeline accounted for 50%, hydrophobic pipeline 50% of the gas liquid profile (water can’t seepage me- dium) 图7. 亲水管道占50%、憎水管道占50%的气液分布图(水不能渗出 介质) 010 2030 40 5060 0 10 20 30 40 50 60 x y Figure 8. Close water pipeline accounted for 60%, hydrophobic pipeline 40% of the gas liquid profile 图8. 亲水管道占60%、憎水管道占40%的气液分布图 两相渗流的拓扑网络模型算法和图形 Copyright © 2013 Hanspub 27 010 2030 40 5060 0 10 20 30 40 50 60 x y Figure 9. Close water pipeline accounted for 80%, hydrophobic pipeline 20% of the gas liquid profile 图9. 亲水管道占80%、憎水管道占20%的气液分布图 0.40.42 0.440.46 0.480.50.52 0.540.56 0.580. 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 10. Chart water seepage probability curve 图10. 水渗出概率曲线 中。当亲水、憎水管到约各占一半时,会出现多种多 样的图形,有些图形中水较均匀地分布在网格的下半 部分,有些图形中水则较零散地分布在整个网格中。 2) 在亲水管道约占 50%的图形中,有些图形出现 了大面积的被水环绕的无水区域(如图 8),通常称之为 圈闭的死气区;还有一些图形的有水区域会沿着很细 的路径向上延伸(如图 9),通常称之为指进。 3) 从水渗出概率曲线(图10)中可以看出,当亲水 管道比例小于约 44%时,水基本无法渗出介质,当亲 水管道比例大于约 56%时,水一定能渗出介质,当亲 水管道比例处于 44%到56%之间时,水渗出介质是有 一定概率的,且概率随比例增大而增大。 7. 结论 本文将多孔介质简化为由节点和管道组成的拓 扑网络模型,并采用拉格朗日算法编程计算,绘制出 两相渗流情形下介质中气液分布的网络图,图形能够 表达更加细观的圈闭和指进等渗流图形,还可以由此 模型得出水渗出整个介质的概率。 参考文献 (References) [1] J. Bockris, S. Srinivasan. Fuel cell: Their electrochemistry. New York: McGraw-Hill Book Company, 1969: 230-253. [2] F. Huang. Three-dimensional reconstruction and simulation of porous media. University of Science and Technology of China, 2007: 5-122. [3] S. P. Guo, D. S. Yu and W. D. Wu. The physical characteristics of the porous media concerning flow in viscera. Acta Mechanica Sinica, 1982, (1): 29-36, 116-118. [4] X. Z. Xu, M. Zhu, F. Li and H. Song. Research on microstruc- ture of the gas diffusion electrode by the small pore diffusion law. Chinese Battery Industry, 2008, 13(1): 1-3. [5] B. M. Yu, P. Xu and J. Chen. Study on the seepage characteris- tics in fractured network porous media. Journal of Xi’an Shiyou University (Natural Science Edition), 2007, (2): 7-8, 37-38, 41. [6] P. Chen, C. F. Wei, J. L. Wang, P. P. Yi and H. F. Cao. Numerical analysis of seepage processes in unsaturated porous media under nearly saturated conditions. Rock and Soil Mechanics, 2012, (1): 299-304. 两相渗流的拓扑网络模型算法和图形 Copyright © 2013 Hanspub 28 附录:主要计算程序代码 M = 61;N = 60; A = rand(M,N);B = rand(N,M) ;P = zeros(M); Z = zeros(1,2*M*M);top = −1; p = input(''); for j = 1:M P(1,j) = 1; end for j = 1:M*N if A(j) < p A(j) = 1; else A(j) = 0; end if B(j) < p B(j) = 1; else B(j) = 0; end end for k = 1:M top = top + 2; Z(top) = 1;Z(top + 1) = k; for j = 1:M*M if Z(top) < M&&B(Z(top),Z(top + 1)) = 1&&P(Z(top) + 1,Z(top + 1)) = 0 Z(top + 2) = Z(top) + 1; Z(top + 3) = Z(top + 1); top = top + 2; P(Z(top),Z(top + 1)) = 1; elseif Z(top + 1) < M&&A(Z(top),Z(top + 1)) = 1&&P (Z(top),Z(top + 1) + 1) = 0 Z(top + 2) = Z(top); Z(top + 3) = Z(top + 1) + 1; top = top + 2; P(Z(top),Z(top + 1)) = 1; elseif Z(top) > 1&&B(Z(top)-1,Z(top + 1)) = 1&&P(Z(top)-1,Z(top + 1)) = = 0 Z(top + 2) = Z(top) − 1; Z(top + 3) = Z(top + 1); top = top + 2; P(Z(top),Z(top + 1)) = 1; elseif Z(top + 1) > 1&&A(Z(top),Z(top + 1)−1) = 1&&P(Z(top),Z(top + 1) − 1) = 0 Z(top + 2) = Z(top); Z(top + 3) = Z(top + 1) − 1; top = top + 2; P(Z(top),Z(top + 1)) = 1; else top = top − 2; end end end end if top = −1 break; end end end for m = 1:M for n = 1:N if P(m,n) = 0&&P(m,n + 1) = 0 A(m,n) = 0; end if P(n,m) = 0&&P(n + 1,m) = 0 B(n,m) = 0; end end end for m = 1:M for n = 1:N x = n−1:n; if A(m,n) = 1 plot(x,x*0 + m − 1); else plot(x,x*0 + m − 1,'r'); end holdon; if B(n,m) = 1 plot(x*0 + m − 1,x); else plot(x*0 + m − 1,x,'r'); end holdon; end end |