Numerical analysis of fully coupled thermoporoelastic behavior of two-dimensional saturated porous media flat plate under LTNE condition
-
摘要:
针对目前多孔介质热流固耦合数学模型多为不完全耦合模型,且物理模型多为一维或一维轴对称模型的现状,对二维饱和多孔介质平板热流固完全耦合问题进行数学建模和数值分析. 采用强耦合方式实现热流固3个物理场的完全耦合. 在平板左侧施加30 ℃流体和固体温度边界条件及固定位移边界条件,平板其他边界温度场和位移场假设为自然边界条件,平板四周加载0 Pa的孔隙压力边界条件. 利用COMSOL Multiphysics有限元分析软件的偏微分方程(PDEs)模式实现上述完全耦合模型的求解,得到渗流场、应变场以及双温度场的数值解. 数值结果表明,随时间增加,流体和固体温度、应变及孔隙压力沿x轴方向传递,同时发现x轴方向应变远大于y轴方向应变. 此外,随时间增加,y轴方向应变和孔隙压力最大值逐渐减小. 研究建立的数学模型和数值解有助于深入认识二维饱和多孔介质热流固完全耦合力学行为.
Abstract:Most of current thermoporoelastic models of porous media are incomplete coupled ones, and the physical models mostly adopt one-dimensional or one-dimensional axisymmetric models. In view of this point, a fully coupled model of a two-dimensional saturated porous media flat plate was presented and the corresponding numerical simulation was carried out. The fully coupled of three physical fields (heat, fluid flow, and solid stress/deformation) were realized by means of strong coupling. The boundary conditiorns were assumed as follows: As for the temperature field and the displacement field, 30 ℃ fluid and solid temperature boundary conditions and fixed displacement boundary conditions were exerted on the left side of the plate, and naturual boundary conditions were prescribed on the other sides of the plate. Meanwhile, 0 Pa pore pressure boundary conditions were loaded around the plate. The numerical solutions of flow field, strain field and two temperature fields were obtained by the PDEs pattern of COMSOL Multiphysics FEA software. Numerical results showed that fluid and solid temperature, strain and pore pressure transfer along the x-axis direction with increasing time, and the strain in the x-direction is much larger than that in the y-direction. In addition, it was found that the maximum values of y-direction strain and pore pressure gradually decrease with time. The presented mathematical model and numerical solutions are of benefit to provide deep insights into the fully thermoporoelastic coupled behavior of two-dimensional saturated porous media.
-
无论在生产还是生活中,多孔材料都非常常见,如岩土、海绵、泡沫隔音材料、锂离子电池隔膜等[1-4]. 研究发现,如果多孔介质温度场方程采用单一温度场模型描述[5-7],得到的温度场数值解或解析解和实验结果往往有所差别. 为更精准地描述多孔介质中温度场分布情况,学者们提出基于局部热非平衡(Local Thermal Non-Equilibrium)条件的双温度模型,即LTNE传热理论. 该理论将多孔介质骨架温度和孔隙流体温度分开考虑,建立固相和液相各自的温度场模型. 该双温度模型结果更为可靠[8]. 多孔介质工作时通常处于多物理场共同作用的环境,因此除了温度场,对多孔介质的研究还会考虑应力应变场[9]、渗流场[10]以及传质[11]等.
多孔介质的复杂结构会导致解析求解或数值求解时一定的困难. 为方便求解,物理模型多采用较为简单的模型,如多孔介质圆管[12-13]、圆环[14]、平板[15]、无限大平面[16]等,且多涉及到一维、一维轴对称、二维. 涉及三维的案例较少,且采用的数学模型也较为简单. 此外,考虑到物理模型通常为部分填充多孔介质[17-18]、复合金属泡沫[19]、多相流[20-21]、裂隙流[22-23]等.
多孔介质数学模型发展较早[24]. 已有不少工作对多孔介质的温度场−渗流场−应力场(以下简称热流固)耦合数学模型进行了研究. 如孔祥言等[25]建立弱耦合的控制方程组,该控制方程组基于线弹性理论,并考虑孔隙率和渗透率随压力和温度的变化. Abousleiman等[7]对横观各向同性储层中倾斜井眼情形进行三场耦合研究,得到解析解. 结果表明钻井深度对整体稳定性存在潜在影响. Zhang等[26]基于局部动力方法建立热流固耦合模型,通过引入颗粒熵的概念得到非弹性力学本构关系,文中模型能够更加准确地描述饱和土壤的热固结问题.
除上述耦合方式外,还可以通过强耦合的方式实现热流固的三场耦合. Yang等[27]基于多孔热弹性理论,对中间有孔的无限大饱和多孔介质平板受冲击热和机械载荷工况时的情形进行研究,采用拉普拉斯变换得到温度场、孔隙压力场及应力位移场. Gao等[16]对类似的物理模型进行研究,只是将孔的边界条件考虑为可渗透与不可渗透两种情形. 这两种耦合虽为强耦合方式,但均为非完全耦合,且物理模型均为一维轴对称. 鉴于此,本研究采用完全强耦合模型,研究基于LTNE的二维多孔介质平板热流固耦合问题.
1. 物理问题和数学模型
1.1 控制方程组
饱和多孔介质平板示意图如图1所示. 其中灰蓝色区域为多孔介质平面,x和y坐标方向如图所示,设平板长和宽分别为0.30 和0.15 m. 假设多孔介质为均质、各向同性和线弹性介质;且被单一流体充满,流体流动符合达西定律;骨架变形符合小变形假设;多孔介质为低渗透,忽略流体对流对传热的影响.
在上述假设之下,基于Li等[12-13, 27-28]研究成果,推导得到LTNE条件下二维饱和多孔介质热流固完全耦合控制方程组为
G(∂2u∂x2+∂2u∂y2)+(λ+G)∂2u∂x2−αB∂P∂x−βTK∂Ts∂x=0 (1) G(∂2v∂x2+∂2v∂y2)+(λ+G)∂2v∂y2−αB∂P∂y−βTK∂Ts∂y=0 (2) −c(∂2P∂x2+∂2P∂y2)+aPP∂P∂t+αB(∂2u∂x∂t+∂2v∂y∂t)+aPTs∂Ts∂t+aPTf∂Tf∂t=0 (3) −(1−ϕ)ks(∂2Ts∂x2+∂2Ts∂y2)+aTsTs∂Ts∂t+T0sβTK(∂2u∂x∂t+∂2v∂y∂t)+T0saPTs∂P∂t+h(Ts−Tf)=0 (4) −ϕkf(∂2Tf∂x2+∂2Tf∂y2)+aTfTf∂Tf∂t+T0faPTf∂P∂t+h(Tf−Ts)=0 (5) 式中:
u 、v 、P 、Ts 和Tf 为控制方程组的因变量,分别表示多孔介质x和y方向的位移、孔隙流体压力以及孔隙流体和固体骨架温度;λ 和G 为Lame常数;αB 为Biot系数;βT 为固体体积热膨胀系数;K 为多孔介质体积弹性模量.k 和μ 分别为渗透率和流体黏度,aPP 为孔隙压力随时间变化量的贡献系数;aPTS 和aPTf 分别为固体和流体温度对孔隙压力的贡献系数;ϕ 为多孔介质的孔隙度;ks 和kf 分别为固体热导率和流体热导率;ρs 和ρf 为固体和流体的密度;C(v)s 和C(p)f 分别为固体的恒容热容和流体的恒压热容;h 为固液界面传热系数,其余参数为aTsTs=(1−f)ρsC(v)s ,aTfTf=ϕρfC(P)f ,c=kμ .方程(1)和(2)分别为x和y方向应力场平衡方程. 由于固体温度相比于加权温度对热膨胀贡献更精确[27],因此方程只使用了固体温度,没有使用加权温度. 观察方程(1)和(2)不难发现,方程中除主因变量位移场外,还包含因变量孔隙压力和温度对位移场的贡献. 换言之,在应力场平衡方程中,位移和压力及温度是完全强耦合的. 方程(3)为渗流场方程,而方程(4)和(5)为LTNE条件下的双温度方程. 同理,容易看出方程(3)~(5)也是完全强耦合形式. 总之,上述控制方程组属多孔介质热流固完全强耦合方程组,较弱耦合/解耦模型而言,更能准确反映多孔介质多场耦合的物理实质.
1.2 边界条件
为求解控制方程组式(1)~(5),给定如下位移场、渗流场和温度场边界条件.
位移场边界条件为
u(x,y,t)|x=0=v(x,y,t)|x=0=0 (6) ∇u(x,y,t)|x=0.3,y=0,y=0.15=∇v(x,y,t)|x=0.3,y=0,y=0.15=0 (7) 孔隙压力场边界条件为
P(x,y,t)|x=0,x=0.3,y=0,y=0.15=0 (8) 流体和固体温度场边界条件为
Ts(x,y,t)|x=0=Tf(x,y,t)|x=0=30℃ (9) ∇Ts(x,y,t)|x=0.3,y=0,y=0.15=∇Tf(x,y,t)|x=0.3,y=0,y=0.15=0 (10) 另假设初始时刻(t = 0),位移和孔隙压力均为0,流体和固体初始温度为0 ℃.
控制方程组式(1)~(5)和边界条件式(6)~(10)及初始条件构成本研究问题的数学模型和初边值问题.
2. 数值方法
鉴于上述数学模型解析解求解困难,本研究使用COMSOL Multiphysics有限元分析软件对其进行数值求解. COMSOL Multiphysics软件开发PDEs(偏微分方程组)求解模式,即用户可自定义所研究问题的复杂控制方程组,并按照软件约定将该复杂偏微分方程组导入软件的PDEs模式,从而实现数值求解.
如前文所述,该数学模型控制方程组式(1)~(5)属自行推导得到的热流固三场完全强耦合方程组,区别于软件模块组合内置的弱耦合方程组. 因此,可利用PDEs模式对方程组式(1)~(5)进行求解. 为保证结果精度,分别对不同网格数(
180×80 ,190×90 和200×100 )下不同坐标点((0.1,0.075),(0.2,0.075)和(0.3,0.075))处500 s时的y向应变值进行对比分析,见表1. 可以发现其相对误差绝对值在1‱以内. 为提高仿真结果精确性,本研究采用200×100 的网格进行计算.表 1 不同网格数下孔隙压力的相对误差Table 1. Relative error of pore pressure under different grid numbers网格数量
坐标180×80 190×90 200×100 (0.1,0.075) 7.96373×10−8 7.96366×10−8 7.96311×10−8 y向应变相对误差/% −0.00088 −0.00691 (0.2,0.075) 1.45365×10−8 1.45372×10−8 1.45339×10−8 y向应变相对误差/% 0.00482 −0.02271 (0.3,0.075) 2.41816×10−9 2.41843×10−9 2.41734×10−9 y向应变相对误差/% 0.01116 −0.04509 此外,本研究通过PDEs模块两种数学应用模式(一般型和系数型)进行计算,验证结果统一性. 所研究问题相关物性参数取值见表2. 部分参数取自文献[27-29],值得说明的是,多孔介质为砂岩,流体为水.
表 2 流体和固体物性参数Table 2. Fluid and solid physical parameters参数/单位 取值 Lame常数G/GPa 6.8 Lame常数λ/GPa 3.8 Biot系数αB 0.74 流度c/(m2·s−1) 1.4×10−3 排空体积模量K/GPa 8.4 参考温度T0/℃ 30 固体密度ρs/(kg⋅m−3) 2600 流体密度ρf/(kg⋅m−3) 1000 孔隙度ϕ 0.4 固体体积热膨胀系数βT/K−1 3.3×10−5 固体温度对孔隙压力的贡献系数aPTS/Pa−1 −7.17×10−5 流体温度对孔隙压力的贡献系数aPTf/Pa−1 −4.78×10−5 固体热导率ks/(W·(m·K)−1) 2.4 流体热导率kf/(W·(m·K)−1) 0.6 恒容热容C(v)s/(J·(kg·K)−1) 920 恒压热容C(p)f/(J·(kg·K)−1) 4200 固液界面传热系数h/(W·(m3·K)−1) 50 孔隙压力随时间变化量的贡献系数aPP/Pa−1 1.4×10−3 3. 结果与讨论
3.1 验 证
考虑到本研究所考察的基于LTNE的二维饱和多孔介质热流固完全耦合模型的解析解难以获取,其数值解也鲜见报道. 为方便起见,此处将二维模型退化为一维轴对称情形,即物理模型与文献[29]相同,为一中间钻孔的无限大多孔介质平面. 对比二者在同一时刻t = 2817 s时的温度场分布特征不难发现:二者结果吻合很好,在一定程度上验证了本研究数学模型和数值结果的正确性,如图2所示.
接下来给出本研究所考察二维问题的有限元解,并对双温度场、孔隙压力场和应变场的分布特征进行分析.
3.2 温度分布
图3为第1000 s时固体温度和流体温度的分布情况. 从图中可以看出,固体温度比流体温度传递的更远,这是由于固体温度热导率大于流体温度的热导率. 图4为500 s和1000 s时,
y=0.075m 处(如无特殊说明,本研究均取该处值)流体和固体温度的分布情况. 从图中可以看出,随着时间增加,温度沿x轴方向逐渐向远处传递.3.3 孔隙压力分布
图5为0.1、1、100、1000 s时孔隙压力分布情况. 由图可知,初始时孔隙压力在
x=0 附近最大,之后沿x轴方向孔隙压力较大值的范围在逐渐扩大,同时靠近边界附近孔隙压力逐渐减小并趋于0. 这是由于随着时间增加,边界条件的设定以及孔隙压力传递和释放导致的. 图6为孔隙压力随x的变化情况. 图中反映出相同特征,同时也可以发现压力最大值在减小,这同样是由于随着时间增加压力释放导致的.3.4 应变分布
图7为1000 s时x方向应变及0.1、1、100、1000 s时y方向应变分布图. 由图7(a)可以看出,x方向应变在
x=0 处最大,沿x轴方向逐渐变小. 由图7(b)至7(e)可以看出,y方向应变初始时刻在x=0 两个角处最大,然后扩张连在一起,之后两个角处应变减小且中心部分应变沿x轴移向远处. 这是因为初期由于位移边界条件的限制,热膨胀导致两角处出现大的y向应变. 随着时间增加,内部孔隙压力变化,导致出现上述应变变化.图8为应变随x轴变化情况. 由图8(a)可以看出,x方向应变随时间增加而增加,并且向x轴正方向移动. 由图8(b)可以看出,y向应变的最值随时间增加先增大后减小,并沿x轴正向移动. 同时,y向应变相比于x向应变要小很多.
4. 结 语
利用COMSOL Multiphysics有限元分析软件对LTNE条件下饱和多孔介质平板热流固三场完全耦合数学模型进行数值求解,分析应变场、孔隙压力场及双温度场的分布和演化特征, 主要得出如下结论.
1)y向应变要远小于x向应变. x向应变随时间增加逐渐增大,而y向应变在初始时在约束边两脚处最大,之后最大值向着中心移动.
2)孔隙压力在初始一段时间内,温度最高时,压力值最大,之后压力最大值沿x轴向远处移动,同时压力最大值在降低.
3)随着时间增加,流体和固体温度沿x轴向远处移动.
-
表 1 不同网格数下孔隙压力的相对误差
Table 1. Relative error of pore pressure under different grid numbers
网格数量
坐标180×80 190×90 200×100 (0.1,0.075) 7.96373×10−8 7.96366×10−8 7.96311×10−8 y向应变相对误差/% −0.00088 −0.00691 (0.2,0.075) 1.45365×10−8 1.45372×10−8 1.45339×10−8 y向应变相对误差/% 0.00482 −0.02271 (0.3,0.075) 2.41816×10−9 2.41843×10−9 2.41734×10−9 y向应变相对误差/% 0.01116 −0.04509 表 2 流体和固体物性参数
Table 2. Fluid and solid physical parameters
参数/单位 取值 Lame常数G/GPa 6.8 Lame常数λ/GPa 3.8 Biot系数αB 0.74 流度c/(m2·s−1) 1.4×10−3 排空体积模量K/GPa 8.4 参考温度T0/℃ 30 固体密度ρs/(kg⋅m−3) 2600 流体密度ρf/(kg⋅m−3) 1000 孔隙度ϕ 0.4 固体体积热膨胀系数βT/K−1 3.3×10−5 固体温度对孔隙压力的贡献系数aPTS/Pa−1 −7.17×10−5 流体温度对孔隙压力的贡献系数aPTf/Pa−1 −4.78×10−5 固体热导率ks/(W·(m·K)−1) 2.4 流体热导率kf/(W·(m·K)−1) 0.6 恒容热容C(v)s/(J·(kg·K)−1) 920 恒压热容C(p)f/(J·(kg·K)−1) 4200 固液界面传热系数h/(W·(m3·K)−1) 50 孔隙压力随时间变化量的贡献系数aPP/Pa−1 1.4×10−3 -
[1] JIANG P X, REN Z P. Numerical investigation of forced convection heat transfer in porous media using a thermal non-equilibrium model[J] . International Journal of Heat and Fluid Flow,2001,22:102 − 110. doi: 10.1016/S0142-727X(00)00066-7 [2] VAFAI K, KIM S J. Forced convection in a channel filled with a porous medium: An exact solution[J] . Journal of Heat Transfer,1989,111:1103 − 1106. doi: 10.1115/1.3250779 [3] NIELD D A, BEJAN A. Convection in porous media [M]. fifth edition. Berlin: Springer, 2017. [4] 马德正, 李培超, 张恒运. 锂离子电池隔膜在压缩过程中的流固耦合效应[J] . 储能科学与技术,2021,10(2):483 − 490. [5] NAIR R, ABOUSLEIMAN Y, ZAMAN M. A finite element thermoporoelastic model for dualporosity media[J] . International Journal for Numerical and Analytical Methods in Geomechanics,2004,28:875 − 898. doi: 10.1002/nag.336 [6] WU B S, ZHANG X, JEFFREY R G, et al. A semi-analytic solution of a wellbore in a nonisothermal low-permeability porous medium under non-hydrostatic stresses[J] . International Journal of Solids and Structures,2012,49(13):1472 − 1484. doi: 10.1016/j.ijsolstr.2012.02.035 [7] ABOUSLEIMAN Y, EKBOTE S. Solutions for the inclined borehole in a thermoporoelastic transversely isotropic medium[J] . Journal of Applied Mechanics,2005,72(1):102 − 114. doi: 10.1115/1.1825433 [8] 路朗, 辛成运, 刘忠鑫. 多孔介质局部非热平衡模型研究综述[J] . 热能动力工程,2019,34(7):1 − 8. [9] 陈帅, 李波波, 张尧, 等. 页岩气储层微观渗流机理研究[J] . 中国科学: 技术科学,2021,51(5):580 − 590. [10] LI P C, ZHANG J L, WANG K Y, et al. Heat transfer characteristics of thermally developing forced convection in a porous circular tube with asymmetric entrance temperature under LTNE condition[J] . Applied Thermal Engineering,2019,154:326 − 331. doi: 10.1016/j.applthermaleng.2019.03.109 [11] 周峰, 刘琪英, 王晨光, 等. 基于COMSOL的果糖脱水传热传质数值模拟[J] . 太阳能学报,2019,40(6):1677 − 1683. [12] LI P C, ZHANG J L, WANG K Y, et al. Analysis of thermally developing forced convection heat transfer in a porous medium under local thermal non-equilibrium condition: A circular tube with asymmetric entrance temperature[J] . International Journal of Heat and Mass Transfer,2018,127:880 − 889. doi: 10.1016/j.ijheatmasstransfer.2018.08.081 [13] YUE F L, LI P C, Zhao C Y. Numerical investigation of thermally developing non-Darcy forced convection in a porous circular duct with asymmetric entrance temperature under LTNE condition[J] . Transport in Porous Media,2021,136(2):639 − 655. doi: 10.1007/s11242-020-01533-7 [14] WANG K Y, TAVAKKOLI F, VAFAI K. Analysis of gaseous slip flow in a porous micro-annulus under local thermal non-equilibrium condition: An exact solution[J] . International Journal of Heat and Mass Transfer,2015,89:1331 − 1341. doi: 10.1016/j.ijheatmasstransfer.2015.06.001 [15] KUZNETSOV A V, NIELD D A. Forced convection in a channel partly occupied by a bidisperse porous medium: Asymmetric case [J]. International Journal of Heat and Mass Transfer, 2010, 53(23/24): 5167−5175. [16] GAO J J, DENG J G, LAN K, et al. Porothermoelastic effect on wellbore stability in transversely isotropic medium subjected to local thermal non-equilibrium[J] . International Journal of Rock Mechanics and Mining Sciences,2017,96:66 − 84. [17] YASSER M, MEHDI M. Analytical investigation of heat transfer enhancement in a channel partially filled with a porous material under local thermal non-equilibrium condition[J] . International Journal of Thermal Sciences,2011,50(12):2386 − 2401. doi: 10.1016/j.ijthermalsci.2011.07.008 [18] CHIKH S, BOUMEDIEN A, BOUHADEF K, et al. Analytical solution of non-Darcian forced convection in an annular duct partially filled with a porous medium[J] . International Journal of Heat and Mass Transfer,1995,38:1543 − 1551. doi: 10.1016/0017-9310(94)00295-7 [19] XU Z G, GONG Q. Numerical investigation on forced convection of tubes partially filled with composite metal foams under local thermal non-equilibrium condition[J] . International Journal of Thermal Sciences,2018,133:1 − 12. [20] TONG F G, JING L R, ZIMMERMAN R W. A fully coupled thermo-hydro-mechanical model for simulating multiphase flow, deformation and heat transfer in buffer material and rock masses[J] . International Journal of Rock Mechanics & Mining Sciences,2010,47(2):205 − 217. [21] LEWIS R W, ROBERTS P J, SCHREFLFLER B A. Finite element modeling of two phase heat and fluid flow in deforming porous media[J] . Transport in Porous Media,1989,4:319 − 334. doi: 10.1007/BF00165778 [22] GELET R, LORET B, KHALILI N. Borehole stability analysis in a thermoporoelastic dual porosity medium[J] . International Journal of Rock Mechanics and Mining Sciences,2012,50:65 − 76. doi: 10.1016/j.ijrmms.2011.12.003 [23] GHASSEMI A, NYGREN A, CHENG A. Effects of heat extraction on fracture aperture: a poro- thermoelastic analysis[J] . Geothermics,2008,37:525 − 539. doi: 10.1016/j.geothermics.2008.06.001 [24] BOER R D. 多孔介质理论发展史上的重要成果 [M]. 刘占芳, 严波, 译. 重庆: 重庆大学出版社, 1995. [25] 孔祥言, 李道伦, 徐献芝, 等. 热−流−固耦合渗流的数学模型研究[J] . 水动力学研究与进展,2005,20(2):269 − 275. [26] ZHANG Z, CHENG X. A fully coupled THM model based on a non-equilibrium thermodynamic approach and its application[J] . International Journal for Numerical and Analytical Methods in Geomechanics,2017,41(4):527 − 554. doi: 10.1002/nag.2569 [27] YANG Y, KLAUS G, TOM S. Thermo-osmosis effect in saturated porous medium[J] . Transport in Porous Media,2014,104(2):253 − 271. doi: 10.1007/s11242-014-0332-5 [28] LI W D, CHEN M, JIN Y, et al. Effect of local thermal non-equilibrium on thermoporoelastic response of a borehole in dual-porosity media[J] . Applied Thermal Engineering,2018,142:166 − 183. doi: 10.1016/j.applthermaleng.2018.06.055 [29] HE L W, JIN Z H. A local thermal nonequilibrium poroelastic theory for fluid saturated porous media[J] . Journal of Thermal Stresses,2010,33:799 − 813. doi: 10.1080/01495739.2010.482358 -