代 宇， 沈 明， 陈 晖

(1. 福州大学数学与计算机科学学院， 福建 福州 350116； 2. 福州大学机械工程及自动化学院， 福建 福州 350116)

0 引言

近年来， 流体在多孔介质中的边界层流动问题备受关注. 多孔介质中流体流动的模型也普遍存在于工业、 医学等各个领域中， 比如海水入侵、 分离加工以及人工透析等. Ali等[1]通过同伦分析方法， 研究了在多孔介质可渗透拉伸板附近的非定常边界层流动； Staff等[2]分析了可压缩流体在多孔介质介质通道的滑移边界层流动； Zhu等[3]研究了在多孔介质中震荡的达西流动； Fang等[4]研究了具有抽吸的收缩表面上的非定常粘性流动问题； Khalid等[5]应用有限元方法探讨了理想流体在多孔介质中的二维稳态流动情况； Rasidi等[6]研究了磁场强度的变化对流体在多孔介质中驻点流动的影响. 对多孔介质中流体流动的研究有着巨大的前景.

对非线性方程求解的方法有很多种， 其中同伦分析方法是由Liao等[7-9]根据拓扑学中同伦的概念所提出的一种求解强非线性方程的有效方法. 目前同伦分析方法已经在流体力学、 海洋工程和金融等诸多领域得到了广泛的应用， 如Xu等[10]应用同伦分析方法研究了在水平通道中的充分发展生物流问题； Danarvand[11]研究了在旋转圆盘上的非中心驻点流动. 同伦分析方法对于非线性的偏微分方程的求解也同样适用， Fan等[12]分析了在伸缩板上的非稳态驻点流动和传热问题. 尤翔程[13]讨论了同伦分析方法在流体力学和时滞动力系统中的应用.

本研究分析不可压缩流体在多孔介质非线性延伸壁面上的驻点边界层流动. 引入无量纲变换， 将控制方程化简为非线性的偏微分方程， 再利用同伦分析方法求得该偏微分方程的非相似解. 最后讨论级数解的收敛性， 并分析多孔介质渗透率、 孔隙率等重要参数对流体速度分布和表面摩擦力的影响.

1 控制方程

考虑到孔隙率以及渗透率等因素的影响， 不可压缩流体在多孔介质中流动时动量方程的模型[14]为

(1)

其中：ρ表示流体的密度； 向量q表示流体速度； d/dt表示随体导数；p为压力；S表示应力张量； 多孔固体基质导致的达西阻力r满足：

(2)

其中：μ表示流体的运动粘性系数；ε是多孔介质的孔隙率；k代表多孔介质的渗透率.

引入达西流速V=(u,v)=εq， 并将式(2)代入方程(1)， 并对方程进行二维边界层， 近似可得

(3)

其中:Sxy是应力张量S的分量， 且满足：

(4)

对于不可压缩流体在多孔介质非线性延伸壁面上的驻点边界层流动， 假设边界层外部的主流速度满足ue=axm， 其中a是正常数，m是非线性参数， 当m=1时退化为线性的情况. 壁面的延伸速度为uw=bxm，b是正常数. 将外部主流流速以及式(4)代入方程(3)， 得：

(5)

将式(5)再代入式(3)中可得

(6)

考虑稳态条件下的该驻点边界层流动问题， 则控制方程为：

(7)

相应的边界条件为：

u(x, 0)=uw,v(x, 0)=0,u(x, ∞)=ue

(8)

引入下面的无量纲变换：

(9)

其中:ψ表示流函数， 则连续性方程自动满足， 且有

u=axmfη(x,η)

(10)

(11)

将式(10)和(11)代入方程(7)得：

(12)

为了简化方程(12)， 在非线性的情况下引入一个新的变换：

(13)

在新的坐标(ξ,η)下， 控制方程(12)变为：

(14)

相应的边界条件转化为：

f(ξ, 0)=0，fη(ξ, 0)=α，fη(ξ, ∞)=1

(15)

除速度场， 表面摩擦系数也是一个重要的物理参数， 其表达式为

2 同伦分析方法的应用

根据同伦分析方法的原则， 找到一个合适的初始猜测解对于解该非线性的偏微分方程至关重要， 其初始猜测解必须满足方程的边界条件. 本文选择的初始猜测解为：

f0(ξ,η)=(α-1)ηe-η+η

(16)

根据初始猜测值的选择， 相应的线性算子可以定义为：

(17)

而Lf[C1+C2η+C3e-η]=0. 其中C1=C1(ξ),C2=C2(ξ),C3=C3(ξ).

根据同伦分析方法的理论基础， 零阶形变方程构造为

(18)

由方程(12)， 非线性算子应取为：

(19)

(20)

(21)

通过泰勒级数展开可知：

(22)

(23)

n阶形变方程与边界条件为：

(24)

fη(ξ, 0)=0,fnη(ξ, 0)=0,fnη(ξ, ∞)=0

(25)

其中

(27)

所以方程(24)的通解为：

(28)

最后应用数学软件Mathematics对方程(24)进行求解， 即可得到方程(14)和(15)的级数解.

3 结果分析

根据同伦分析方法的理论， 同伦级数解(22)的敛散性和收敛速度依赖于辅助参数h的选取， 所以选择合适的辅助参数h至关重要. 图1、 2分别给出了α=0.5和α=2.0时在不同参数下的辅助参数h曲线图. 当α=0.5时， 在各个参数的变化下，h取值的合理范围为： -1.5≤h≤-0.5； 当α=2.0时， 在各个参数的变化下，h取值的合理范围为-1.2≤h≤-0.5. 结合图1和图2的结果，h的取值范围为： -1.2≤h≤-0.5. 由于同伦分析方法给予了很大的自由空间来选择辅助参数h， 所以在-1.2≤h≤-0.5范围内所选取的辅助参数的值均是合理的. 本研究选取h=-1. 表1列出了在35阶下的fηη(ξ, 0)值, 可以验证所求结果的准确性.

图1 α=0.5时辅助参数h曲线图Fig.1 h curves of fηη(ξ, 0) for α=0.5

图2 α=2.0时辅助参数h曲线图Fig.2 h curves of fηη(ξ, 0) for α=2.0

nα=0.5α=2.0nα=0.5α=2.011.11863 -2.9173151.23721 -3.22526101.24154 -3.23526151.24155 -3.23574201.24155 -3.23568251.24155 -3.23568301.24155 -3.23568351.24155 -3.23568

图3 ε=0.5, ξ=0.5, m=1.5时， 在不同P值下的速度分布Fig.3 Variations of velocity profile for different values of P when ε=0.5, ξ=0.5, m=1.5

讨论多孔介质孔隙率、 渗透率以及非线性参数等对边界层流动的影响. 当α=0.5时， 外部主流速度大于壁面延伸速度； 当α=2.0时， 外部主流速度小于壁面延伸速度. 图3描述了多孔介质渗透率对速度分布的影响， 在α=0.5和α=2.0两种情况下， 速度边界层厚度随着多孔介质渗透参数P的增大而减小. 图4分析了多孔介质孔隙率对速度分布的影响， 当α=0.5和α=2.0时， 速度边界层厚度均随着孔隙率的增大而增大. 图5研究了流体非线性参数m对速度分布的影响， 在壁面延伸速度大于外部主流速度和小于外部主流速度两种情况下， 随着m的增大， 边界层的厚度均减小， 流体在边界层内流动增强. 在ξ=0.5时， 图3～5均表明, 当α=0.5时， 边界层内的流体速度随着η的增大逐渐增大； 当α=2.0时， 边界层内的流体速度随着η的增大逐渐变小. 由于边界层内流体的速度是从壁面速度过渡到外部主流速度， 所以这一结果也与事实相符， 同时也证明结果的正确性. 此外， 在多孔介质中, 孔隙率参数和非线性参数对流体达西速度剖面的影响作用大于渗透参数的影响.

图4 P=2, ξ=0.5, m=1.5时， 在不同ε值下的速度分布Fig.4 Variations of velocity profile for different values of ε when P=2, ξ=0.5, m=1.5

图5 P=2, ε=0.5, ξ=0.5时， 在不同 m值下的速度分布Fig.5 Variations of velocity profile for different values of m when P=2, ε=0.5, ξ=0.5

图6～7分析了多孔介质渗透率和孔隙率对表面摩擦系数的影响. 当外部主流速度大于壁面延伸的速度时， 表面摩擦系数为正数； 当外部主流速度小于壁面延伸的速度时， 表面摩擦系数为负数. 在α=0.5和α=2.0两种情况下， 当多孔介质渗透参数P增大时， 表面摩擦系数的绝对值均增大; 当多孔介质孔隙率ε增大时， 表面摩擦系数的绝对值均减小， 这与图3、 4所得的结果一致.

图6 ε=0.5, ξ=0.5, m=1.5时， 不同P值下的表面摩擦力Fig.6 Variations of skin friction coefficient for different values P of when ε=0.5, ξ=0.5, m=1.5

图7 P=2, ξ=0.5, m=1.5时， 在不同ε值下的表面摩擦力Fig.7 Variations of skin friction coefficient for different values of ε when P=2, ξ=0.5, m=1.5

4 结语

本研究利用同伦分析方法求得多孔介质中非线性延伸壁面上的驻点边界层流动问题的非相似解. 引入达西流速， 分析多孔介质的孔隙率、 渗透率和非线性参数等对边界层内流体速度分布和表面摩擦系数的影响. 同伦分析方法对非线性的偏微分方程求解得到了应用， 并根据最后的收敛结果验证了结果的正确性.

在此基础上， 根据所得级数解对该边界层问题进行讨论， 研究结果表明： 在速度比例参数α=0.5和α=2.0两种情况下， 边界层内的流动随着多孔介质渗透参数P和非线性参数m的增大而增强， 但随着孔隙率的增大而减弱；α=0.5时表面摩擦力方向与α=2.0时的方向相反， 其绝对值随着多孔介质渗透参数P的增大而增大， 随孔隙率ε的增大而减小.

参考文献：

[1] ALI A, MEHMOOD A. Homotopy analysis of unsteady boundary layer flow adjacent to permeable stretching surface in a porous medium[J]. Communications in Nonlinear Science &Numerial Simulation, 2008, 13(2): 340-349.

[2] STAFF T P O. Correction: modeling and analysis of unsteady axisymmetric squeezing fluid flow through porous medium channel with slip boundary[J]. Plos One, 2015, 10(4): e0124851.

[3] ZHU T, MANHART M. Oscillatory darcy flow in porous media[J]. Transport in Porous Media, 2016, 111(2): 521-539.

[4] FANG T, ZHANG J, Yao S. Viscous flow over an unsteady shrinking sheet with mass transfer[J]. Chinese Physics Letters, 2009, 26(1): 172-175.

[5] KHALID A, KHAN I, KHAN A,etal. Unsteady MHD free convection flow of Casson fluid past over an oscillating vertical plate embedded in a porous medium[J]. Engineering Science & Technology An International Journal, 2015, 8(3): 309-317.

[6] KECHILS A, HASHIM I. Approximate analytical solution for MHD stagnation-point flow in porous media[J]. Communications in Nonlinear Science & Numerical Simulation, 2009, 14(4):1346-1354.

[7] LIAO S J. The proposed homotopy analysis technique for the solution of nonlinear problems[D]. Shanghai: Shanghai Jiao Tong University,1992.

[8 ] LIAO S J. An explicit, totally analytic approximate solution for Blasius’ viscous flow problems[J]. International Journal of Non-Linear Mechanics, 1999, 34(4): 759-778.

[9] LIAO S J, CHWANG A T. Application of homotopy analysis method in nonlinear oscillations[J]. Journal of Applied Mechanics, 1998, 65(4): 914-922.

[10] XU H, Pop I. Fully developed mixed convection flow in a horizontal channel filled by a nanofluid containing both nanoparticles and gyrotactic microorganisms[J]. European Journal of Mechanics, 2014, 46(12): 37-45.

[11] DINARVAND S. On explicit, purely analytic solutions of off-centered stagnation flow towards a rotating disc by means of HAM[J]. Nonlinear Analysis Real World Applications, 2010, 11(5): 3389-3398

[12] FAN T, XU H, Pop I. Unsteady stagnation flow and heat transfer towards a shrinking sheet[J]. International Communications in Heat & Mass Transfer, 2010, 37(10): 1440-1446.

[13] 尤翔程. 同伦分析方法在流体力学和时滞动力系统中的应用[D]. 上海： 上海交通大学, 2011.

[14] ZHAO J, ZHENG L, ZHANG X,etal. Unsteady natural convection heat transfer past a vertical flat plate embedded in a porous medium saturated with fractional Oldroyd-B fluid[J]. Journal of Heat Transfer, 2017, 139(1): 012501.