多孔介质中单相气体局部流动的均质化建模

李树光, 曲 凯

(大连海事大学 理学院, 辽宁 大连 116026)

0 引 言

多孔介质中的气体流动广泛存在于自然现象和工业技术中,因而对其开展数学建模研究具有重要的科学意义和应用价值[1],如地质层中天然气的输运[2],飞机油箱内多孔阀的油气过滤,飞行器外壳热防护材料的气动烧蚀等．然而,由于测量技术的局限性和多孔介质空间结构的复杂性,使得这种流动的研究非常困难,大多数工作仅能获得多孔介质中气体流动的宏观结果．许多研究是根据Darcy定律及其修正,从宏观现象分析孔隙中气体的流动．这种情况下,单孔中所发生的局部效应仅作为整体特征的简单估计(如Darcy定律中平均流速和实验中渗透系数的确定),未充分考虑孔隙中流动的真实过程．孔隙中的流动对渗流的整体影响非常重要[3]:孔隙的几何特征和局部流动决定了渗流的渗透性．因此,有必要准确描述孔隙中的局部流动．

许多学者应用多种方法研究了多孔介质孔隙中的局部流动,如均质化和体积平均模型[4],孔隙网络(“管束”)模型[5],分数阶模型[6]等．每种方法都有其优缺点,对于渗流问题的研究,描述孔隙中发生的实际物理过程以及微宏观过程之间的理论联系至关重要．因此,本文采用了文献[7]中提出的针对周期性复合材料建模的渐近均质法,该方法能够有效模拟多孔介质孔隙中的物理过程．目前,渐近均质法已经被广泛应用于解决多孔材料的热力学问题[8]、复合材料有效性能预测[9-10],以及材料渗透性能估计[11]等问题．

应用渐近均质法获得的描述局部流动的数学模型具有一些特殊性[12],主要是:通常具有积分-微分类型的均匀化限制和周期性边界条件所构成的约束条件．这使得局部问题的理论意义和求解方法有别于经典的数学模型．一些研究数值求解了周期单元上的局部问题,包括罚函数有限元法[13]、光滑粒子流体动力学(SPH)方法[14]等．但是少有讨论局部问题的特殊性和复杂约束条件的求解方法．

本文的目的是研究周期单元上气体的局部流动问题,提出了一种基于对称性和反对称性扩展的简化方法,开发求解局部问题的最小二乘有限元法,通过微管中Poiseuille流动的解析解验证了相关模型和算法．

1 多孔介质中气体流动的连续介质力学模型

1.1 多孔介质的几何模型

考虑如图1所示的正交各向异性多孔介质．建立如下假设: ① 多孔介质的微观结构具有周期性特征; ② 多孔介质无死孔; ③ 周期单元是孔隙中充满气体的纤维和纤维材料所组成的区域; ④ 周期单元关于局部直角坐标系中的坐标平面具有几何和物理的对称性,这使得能够将计算区域简化为1/8周期单元．

(a) 多孔介质的微观结构 (b) 多孔介质的周期单元 (c) 多孔介质的1/8周期单元 (a) The microstructure of the porous medium (b) Periodic cell of the porous medium (c) 1/8 periodic cell of the porous medium

图1中,V表示整个孔隙所占的区域,其边界(固-气边界)用Σ表示,Vξ表示周期单元,Σξ表示孔隙中气体与固体的边界．

1.2 气体流动的力学方程组

多孔介质中气体的流动由可压缩Navier-Stokes方程组描述:

(1)

其中,ρ是密度,v是速度,T是Cauchy应力张量,⊗表示张量积．

对气体介质设置如下假设: ① 孔隙中的气体介质为理想气体; ② 气体是各向同性的Newton黏性介质,黏性系数非常小且不为零; ③ 气体的质量力密度为零; ④ 气体的运动是等温等压的; ⑤ 气体的体积黏度系数[12]为零,即λ+2μ/3=0,其中λ和μ是黏度系数．

根据各向同性线黏性介质的本构关系,Cauchy应力张量可写为

T=-pE+σ,

(2)

其中,p是压力,σ=λ(∇·v)E+2με是黏性应力张量,ε=(∇⊗v+∇⊗vT)/2是应变率张量．

考虑理想气体的状态方程:

p=ρRθ0,

(3)

其中,R是气体常数;θ0是恒定温度．

在固气界面,设置无滑移边界条件:

v|Σ=0．

(4)

初始时刻t=t0的给定压力为

p|t=t0=p0．

(5)

2 多孔介质中气体流动的局部建模

2.1 局部建模中渐近均质法的尺度放大

根据多孔介质的周期性假设,速度v、压力p和密度ρ由小参数κ的渐近展开式表示:

(6)

根据气体黏性的假设,可以得到μ=μ0κ2．

将式(6)和黏性假设代入系统(1)—(5)并以小参数κ的幂次进行分组后,得到周期单元上的局部问题如下:

(7)

进一步考虑三维多孔结构,气体沿着局部坐标轴运动时,由于局部问题(7)的线性性质,其解可表示为“输入数据”∇xp(0)的线性函数:

(8)

(9)

(10)

2.2 局部问题的物理意义

通过分析系统(9)可得结论:在局部层面上,当系统(9)固定一个α的值,也就确定一个伪不可压缩线性流体的稳态流动问题,并且系统(9)的解仅由孔的几何形状决定．因此,系统(9)适用于在上述假设框架下计算多孔介质中气体的流动．

2.3 局部问题在1/8周期单元上的简化

进一步,考虑将局部问题简化为1/8周期单元上的问题,从而进行数值求解．因此,根据解的对称和反对称原理得到以下关于解的延续性定理[11-12]．

(11)

2,3时,可采用类似的置换方式进行．图2中的符号“+”表示从第一象限扩展到其他象限时函数不更改符号(即对称扩展),符号“-”表示符号更改(即反对称扩展)．

利用系统(11)的解进行对称和反对称扩展,在1/8周期单元的边界平面上设置如下边界条件:

(12)

2.4 渗透性的确定

(13)

3 局部问题的最小二乘有限元数值解法

(14)

(15)

(16)

接下来,考虑有限元公式(14)—(16)的近似函数．由于最小二乘公式(14)中近似函数的最小要求是线性的,因此可采用C0基函数[17]．尽管在所用的近似空间之间没有兼容性限制(即inf-sup条件),在本文中仍然使用满足inf-sup条件的近似函数．

4 数 值 结 果

本节中,首先验证了所提出的模型和算法的可靠性和准确性,然后模拟了在图1所示结构中的流动．

4.1 微管中Poiseuille流动的数学模型和数值模拟

为验证模型和算法,考虑单通道结构中α=1时模型(9)的解析解．微管的1/8周期单元如图3所示．

图3 圆柱形单通道结构的1/8周期单元Fig. 3 The 1/8 periodic cell of the cylindrical single channel structure

(17)

(18)

圆柱形单通道结构的1/8周期单元中,系统(18)在柱坐标系下可表述为

(19)

其中,R是圆柱形单通道的半径．根据式(19)可得Cauchy问题:

(20)

在试验中,微管的半径为0.4．因此,孔隙度的精确结果为0.502 654 82．计算微管中Poiseuille流动采用的是2 674个单元、4 698个节点的四面体网格．图4中给出了计算结果,孔隙度的数值结果为0.502 654 65．由此,数值结果验证了模型和数值算法的可靠性和准确性．

(a) Poiseuille流动的数值解 (b) 与精确解的比较 (a) The numerical solution of the Poiseuille flow (b) Comparison of numerical and exact solutions

4.2 正交各向异性多孔介质中局部流动的数值结果

考虑图1所示的几何模型,孔隙中球体的无量纲半径为0.35,圆柱(连接通道)的无量纲半径为0.1．模型采用的参数为:压力量纲p0=106Pa,速度量纲v0=1 m/s,多孔介质的量纲长度L=1 m,周期单元的量纲长度l=10-4m,气体(氮气)的黏度μ=0.018 448 Pa·s．

表1 多孔介质孔隙率和渗透系数的计算结果

(a) 速度分量 速度分量

5 结 论

本文应用渐近均质法建立了描述三维周期性多孔结构孔隙尺度下单相气体局部流动的数学模型,获得了周期单元上的局部问题,明确了局部问题的数学特殊性和物理意义,提出了一种基于对称性和反对称性扩展的简化方法．通过对局部问题的分析,能够准确地获得多孔介质的渗透率．

本文在分析局部问题的基础上,提出了求解局部问题的最小二乘有限元方法,克服了由平均算子和周期性边界条件引起的数值困难．最后,通过理论分析获得了微管中Poiseuille流动的解析解,并验证了所提出的数学模型和数值算法的可靠性和准确性．研究结果表明,应用渐近均质法能够解决多孔介质中的单相气体流动问题,获得多孔介质的气体渗透性．