0 引言

常规油气探测中一般视孔隙流体为牛顿流体。但是，在油气勘探和开发过程中普遍存在非牛顿流体流动，这涉及到许多重要的工程，如用加热方法提高重油采收率；用聚合物溶液、油和泡沫溶液的乳状液充当驱替液；用含有悬浮颗粒的化学溶剂充当页岩气压裂液等。当这类流体与固体之间发生相对运动时会产生剪切应力，应力大小跟流体黏性和流体剪切应变率有关。非牛顿流体剪切应力与剪切应变之间一般不满足线性关系，黏性系数是剪切变化率的复杂函数。在化学工程、流变学、石油工程等领域，对非牛顿流体流变学进行了大量的研究。Tan等^[1]研究了Maxwell黏性流体在平板之间的流动现象，卢明辉^[2]在流体的本构模型中引入了流体抵抗剪切变形的作用，Tong等^[3]给出了圆管中非牛顿流体流动的解析解。但是，孔隙介质中非牛顿流体对弹性波场传播的影响机制仍不清楚^[4]。

孔隙介质波动力学的研究多采用理想流体或具有黏性的牛顿流体，这使公式推导较为简便。但是，致密油储层地质特征复杂，岩石包含矿物颗粒和固态有机质等物质，流体非均质性强^[5]，毛细管压力大^[6]。这类岩石中主要发育微纳米级孔喉连通体系，孔隙结构复杂^[7]、连通性差，孔径范围在纳米、微米量级，且孔隙度一般小于10%、渗透率一般小于0.1mD。研究表明流体在宏观与微观孔隙结构中的流动行为并不一致，当孔隙尺度逐渐变小时流体的连续介质流动特征发生改变^[8]。Thomas等^[9]发现在微纳米孔中流体呈现接近固体的规则分子排布，其黏性特征已发生改变，表明致密孔隙中流体是介于固体与理想流体之间的复杂状态物质。同时，在微纳米孔隙通道条件下，非牛顿流体表现出“非达西”流动特征，因此常规牛顿流体假设已不适用于致密储层，采用非牛顿本构模型更为合理。

Sochi^[10]综述了孔隙介质非牛顿流体的流动，重点叙述含屈服应力的流体、黏弹性流体的性质及流动模型，并讨论了测试流体屈服应力的不同方法。Pearson等^[11]建立了复杂非牛顿流体在孔隙介质中流动的数学模型。Dong等^[12]研究了重油非牛顿流体在孔隙介质中的流动规律，实验结果表明流体流速与压力梯度呈非线性关系，且部分实验中流体流动具有启动压力梯度。Xiong等^[4]通过对三维孔隙网络中的Maxwell流体渗透率建模计算，发现渗透率变化受频率、流体流变学属性和孔隙连通性的影响。Talon等^[13]在格子玻尔兹曼框架下，针对含屈服应力非牛顿流体在孔隙介质中的流动问题，推导了广义达西方程以描述流体的流动，并分析各参数的计算方法。Markov^[14]研究了饱和Maxwell非牛顿流体孔隙介质瑞利面波沿边界层的传播问题，建立了波动方程并进行频散和衰减分析，指出有可能通过分析波速频散和衰减估计Deborah数。

Teeuw等^[15]通过实验手段研究了Bentheim砂岩岩心中高分子溶液的流动特征，实验结果显示，对于较大高分子的溶液，孔隙壁对分子的排斥作用引起了明显的滑移效应，导致孔隙介质中出现较大流动速率；Greaves等^[16]进行了Elginshire砂岩的高分子溶液注入实验，发现高分子溶液在孔隙中的黏性大于在实验室容器中，认为这是由溶液与孔隙壁间的相互作用引起的；Cannella等^[17]在Berea岩心实验中也发现相同的现象；Hejri等^[18]通过实验研究了Ottawa砂岩样本中高分子溶液的流变学行为，发现在溶液流速高于阈值时，非固结砂岩中高分子溶液流动规律不满足线性达西定律，而是满足幂率流模型。

孔隙流体属性是重要的岩石物理参数^[19]，尽管对孔隙尺度非牛顿流体流动现象的研究较为深入，但基于非牛顿流体—固体骨架耦合作用的孔隙介质波动方程研究却尚未引起足够的重视。传统孔隙介质波场传播模型假设孔隙中填充了经典牛顿流体，然而大量理论和实验研究表明，现有牛顿流体孔隙介质波动方程已不适用于低孔、低渗致密油储层的复杂流动条件。为了更好地利用地震技术发现非常规油藏和提高剩余油采收率，人们试图研究非牛顿流体对多孔介质动力学的影响。一种观点认为孔隙流体的非线性流变性可能是低频范围的主要影响机制^[20]，这需要对非牛顿流体有更深入的了解。幂律应变—应力关系是一种表征多孔介质(如Bentheim砂岩^[15])中聚合物水溶液非线性模型的常见本构关系；同时，Maxwell流体是分析黏弹性特性的另一个常用模型，能够描述特定频率下振荡流速提高现象。Tsiklauri等^[21]将波场的共振现象归因于孔隙中非牛顿流体，实验证实了在振荡压力梯度下黏弹性流体的动力响应^[22]。现有非牛顿流体模型研究仍处于不断探索阶段，对实际流体流变行为的机理仍不甚了解，需要更灵活的本构关系解释实验现象。

非牛顿流体模型类型众多，主要可分为黏弹性、非时变黏性和时变黏性流体三类。黏弹性非牛顿流体模型主要包括Maxwell和开尔文流体模型，通过弹簧和阻尼的串联、并联实现弹性和黏性的结合。Maxwell流体的流变特性一般采用串联的弹簧阻尼模型，然而该模型通常过于理想化，无法描述真实的非牛顿流体。近年来，分数阶微分模型受到了相当大的关注，在描述物质黏弹性属性方面具有很大优势。分数阶导数模型使用的是非整数阶导数描述应变—应力关系，其优点在于它的非局部结构，适用于对黏弹性流体的记忆特性进行建模。最近的一些研究分析了Maxwell流体的流动^{[1, 23-24]}，实验结果证实了聚合物流变学与分数阶导数Maxwell(fractional derivative Maxwell，fdMaxwell)模型一致^[25]。目前多孔介质中fdMaxwell流体对弹性波传播影响的研究尚未展开，主要是因为非牛顿流体与固体之间的黏性和弹性耦合机制尚不清楚。

本文从非牛顿流体黏性系数的非线性本构关系出发，基于Biot理论中非泊肃叶流动的框架，考虑振动压力梯度作用下的孔隙微管黏性流体流动，用fdMaxwell模型描述多孔介质中弹性波传播的非牛顿流体效应，通过建立孔隙介质流体动量方程和非牛顿流体本构关系，得到含非牛顿流体的孔隙介质动力学方程组，通过平面波分析方法，得到非牛顿流体效应下波场频散和衰减的计算方法。

1 fdMaxwell流体的黏弹性耗散

考虑在多孔介质中存在一个封闭表面∂Ω，内部包含单元体Ω，单元体Ω由饱和Maxwell黏弹性流体的孔隙介质组成，多孔介质骨架是各向同性的、均匀的线性弹性体。单元体内部物质满足连续性法则。其质量连续性方程为

$ \frac{{\partial {\rho _{\rm{f}}}}}{{\partial t}} + \nabla \cdot \left( {{\rho _{\rm{i}}}{\mathit{\boldsymbol{v}}_{\rm{f}}}} \right) = 0 $

(1)

动量守恒方程为

$ {\rho _{\rm{f}}}\frac{{{\rm{D}}{\mathit{\boldsymbol{v}}_{\rm{f}}}}}{{{\rm{D}}t}} = - \nabla p + \nabla \cdot \mathit{\boldsymbol{\tau }} $

(2)

式中：$\frac{{{\rm{D}}{{\mathit{\boldsymbol{v}}}_{\rm{f}}}}}{{{\rm{D}}t}}$表示流体速度v_f的物质导数；p为流体压力；ρ_f为流体密度；τ是应力张量。对于fdMaxwell模型，本构方程^[26-27]为

$ \mathit{\boldsymbol{\tau }} + {\lambda ^\alpha }{\rm{d}}_t^\alpha \mathit{\boldsymbol{\tau }} = \eta {\lambda ^{\beta - 1}}{\rm{d}}_t^{\beta - 1}\mathit{\boldsymbol{\dot \gamma }} $

(3)

式中：η为流体黏性系数；λ=η/μ_f为弛豫时间，其中μ_f为流体剪切模量；α和β为分数阶导数的阶次；变形速率张量被定义为

$ \mathit{\boldsymbol{\dot \gamma }} = \nabla {\mathit{\boldsymbol{v}}_{\rm{f}}} + {\left( {\nabla {\mathit{\boldsymbol{v}}_{\rm{f}}}} \right)^{\rm{T}}} - 2/3 \cdot \nabla \cdot {\mathit{\boldsymbol{v}}_{\rm{f}}}\mathit{\boldsymbol{I}} $

(4)

表示应变对时间的一阶导数，其中I表示二阶单位张量。

分数阶导数定义^[28]为

$ {\rm{d}}_t^\alpha f(t) = \frac{{{{\rm{d}}^\alpha }f(t)}}{{{\rm{d}}{t^\alpha }}} = \frac{1}{{\mathit{\Gamma }(m - \alpha )}}\int_0^t {\frac{{{f^{(m)}}(\tau )}}{{{{(t - \tau )}^{\alpha + m - 1}}}}} {\rm{d}}\tau $

(5)

式中：m为正整数；m-1≤α≤m；Γ为伽马函数。为了获得流体相的动态方程，对本构方程两边取散度，带入动量守恒方程，可得

$ \begin{array}{*{20}{c}} {\left( {1 + {\lambda ^\alpha }{\rm{D}}_t^\alpha } \right){\rho _{\rm{f}}}\frac{{{\rm{D}}{\mathit{\boldsymbol{v}}_{\rm{f}}}}}{{{\rm{D}}t}} = - \left( {1 + {\lambda ^\alpha }{\rm{D}}_t^\alpha } \right)\nabla p + }\\ {\eta {\lambda ^{\beta - 1}}{\rm{D}}_t^{\beta - 1}\left[ {{\nabla ^2}{\mathit{\boldsymbol{v}}_{\rm{f}}} + \frac{1}{3}\nabla \left( {\nabla \cdot {\mathit{\boldsymbol{v}}_{\rm{f}}}} \right)} \right]} \end{array} $

(6)

上式是流体满足的动力学方程。

对于圆柱形孔隙微管道，流体运动方程存在解析解。在微纳米孔隙介质中，管道半径比孔隙界面的波长要小得多(“长波长近似”假设)，因此固壁局部位移可以看作是常数。在低速条件下，孔隙流体的压缩性可以忽略不计，则质量连续性条件可写为

$ \operatorname{div} \boldsymbol{v}_{\mathrm{f}}=0 $

(7)

在微米量级的孔隙管道里流动的黏弹性流体，其雷诺数比失稳临界值要小得多，层流假设依然合理。因此，在柱坐标系(r，θ，z)下，可以假设唯一的非零速度分量是沿z轴(孔隙管道轴向)方向，而且速度只取决于径向坐标，因此有

$ \left\{\begin{array}{l} \boldsymbol{v}_{\mathrm{f}}=\left(0,0, v_{\mathrm{fz}}\right) \\ \mathrm{D} / \mathrm{D} t=\partial / \partial t \end{array}\right. $

(8)

在孔隙介质中，流体速度可以表示为

$ \boldsymbol{v}_{\mathrm{f}}=\boldsymbol{w}+\dot{\boldsymbol{u}} $

(9)

式中：u是固体位移；$ {\mathit{\boldsymbol{\dot u}}}$为u对时间的一阶导数；w(r)为流体相对速度。如果w(r)与z坐标无关且轴对称，那么流体速度的控制方程可以表示为

$ \begin{array}{*{20}{c}} {{\rho _{\rm{f}}}\left( {1 + {\lambda ^\alpha }{\rm{d}}_t^\alpha } \right)\mathit{\boldsymbol{\dot w}} = - \left( {1 + {\lambda ^\alpha }{\rm{d}}_t^\alpha } \right)\left( {\nabla p + {\rho _{\rm{i}}}\mathit{\boldsymbol{\ddot u}}} \right) + }\\ {\eta {\lambda ^{\beta - 1}}{\rm{d}}_t^{\beta - 1}\left( {{\nabla ^2}\mathit{\boldsymbol{w}} + {\nabla ^2}\mathit{\boldsymbol{\dot u}}} \right)} \end{array} $

(10)

式中ü表示固体位移对时间的二阶导数。在孔隙微管表面上，$ {\mathit{\boldsymbol{\dot u}}}$(z)与轴向坐标有关，在圆柱坐标下的速度拉普拉斯算子为$ {\nabla ^2}{\mathit{\boldsymbol{w}}} = 1/r \cdot \partial \left( {r\partial w/\partial r} \right)/\partial r$和$ {\nabla ^2}{\mathit{\boldsymbol{\dot u}}} = {\partial ^2}{\mathit{\boldsymbol{\dot u}}}/\partial {z^2}$。如果使用$ {\partial ^2}{\mathit{\boldsymbol{w}}}/\partial {r^2} = {\mathit{\boldsymbol{O}}}\left( {{\mathit{\boldsymbol{w}}}/{r^2}} \right)$和$ {\partial ^2}{\mathit{\boldsymbol{\dot u}}}/\partial {r^2} = {\mathit{\boldsymbol{O}}}\left( {{\mathit{\boldsymbol{\dot u}}}/{\mathit{L}^2}} \right)$作为$ {\nabla ^2}{\mathit{\boldsymbol{w}}}、{\nabla ^2}{\mathit{\boldsymbol{\dot u}}}$的近似估计(其中L是发生重大变化的距离(波长))，根据“长波长假设”(r/L«1)，则$ {\partial ^2}{\mathit{\boldsymbol{w}}}/\partial {r^2} \gg {\partial ^2}{\mathit{\boldsymbol{\dot u}}}/\partial {r^2}$。在波场作用下，速度、位移和压强随时间变化满足波函数形式(exp(-iωt))，则控制方程可写为

$ \frac{{{\partial ^2}\mathit{\boldsymbol{w}}}}{{\partial {r^2}}} + \frac{1}{r}\frac{{\partial \mathit{\boldsymbol{w}}}}{{\partial r}} + \frac{{{\rm{i}}\omega A}}{\psi } \cdot \mathit{\boldsymbol{w}} = - \frac{A}{\psi }\mathit{\boldsymbol{f}} $

(11)

式中：ψ=η/ρ_f；f=- ($ \nabla p + {\rho _{\rm{f}}}{\mathit{\boldsymbol{\ddot u}}}$)/ρ_f；ω为波场圆频率；A=(-iD_ec²)^1-β+(-iD_ec²)^α-β+1，其中c= $ a\sqrt {\omega /\psi } $，D_e=ψλ/a²为Deborah数，a为孔隙管道半径。在流—固边界上应用无滑移边界条件w|_r=a=0，可得方程的解为

$ \mathit{\boldsymbol{w}} = - A\left[ {1 - {J_0}\left( {lr} \right)/{J_1}(la)} \right]\mathit{\boldsymbol{f}}/\left( {{l^2}\psi } \right) $

(12)

式中：J₀为零阶贝塞尔函数；J₁为一阶贝塞尔函数；l= $ \sqrt {{\rm{i}}\omega \mathit{A}{\rm{/}}\psi } $。在流体和固体交界面上的剪切应力为速度梯度和黏性系数的函数

$ {\left. \mathit{\boldsymbol{\tau }} \right|_{r = a}} = {\left. {\frac{\eta }{A} \cdot \frac{{\partial \mathit{\boldsymbol{w}}}}{{\partial r}}} \right|_{r = a}} = \frac{\eta }{{l\psi }} \cdot \frac{{{J_1}(la)}}{{{J_0}(la)}}\mathit{\boldsymbol{f}} $

(13)

将平均速度〈w〉=2/a²∫₀^arwdr和总应力〈τ〉=2πaτ|_r=a结合起来，可得剪切作用力为〈τ〉=8πηF(c)〈w〉，其中黏性耗散函数可表示为

$ F(c)=-\frac{\mathrm{i} c}{4 \varDelta} \cdot \frac{J_{1}(c \varDelta)}{J_{0}(c \Delta)}\left[1-\frac{2 J_{1}(c \varDelta)}{c \Delta J_{0}(c \varDelta)}\right] $

(14)

式中

$ \varDelta=\sqrt{\mathrm{i}} \cdot \sqrt{\left(D_{\mathrm{e}} c^{2}\right)^{1-\beta}(-\mathrm{i})^{1-\beta}+\left(D_{\mathrm{e}} c^{2}\right)^{\alpha-\beta+1}(-\mathrm{i})^{\alpha-\beta+1}} $

(15)

Deborah数定义为D_e=λ/λ_ψ，表示了黏性效应的弛豫时间和特征时间的比值，刻画了流体—固相系统是否处于耗散状态或弹性状态，其中λ_ψ=a²/ψ。具有较低Deborah数的材料表现为牛顿黏性流体；当Deborah数增加时，弹性行为占主导地位。

2 饱和fdMaxwell流体的孔隙介质波动方程

在固体壁面上的流—固耦合剪切作用力为〈τ〉=8πηF(c)〈w〉。在Biot理论中，流—固耦合引起的黏性耗散项是ηϕ²($ {\mathit{\boldsymbol{\dot U}}}$-$ {\mathit{\boldsymbol{\dot u}}}$)/κ，其中κ是渗透率，ϕ是孔隙度，$ {\mathit{\boldsymbol{\dot U}}}$、$ {\mathit{\boldsymbol{\dot u}}}$是流体和固体速度。考虑到非牛顿流体效应的修正函数F(c)，黏性耗散项可写为ηϕ²×F(c)($ {\mathit{\boldsymbol{\dot U}}}$-$ {\mathit{\boldsymbol{\dot u}}}$)/κ，其中ηϕ²F(c)κ为一般形式流—固耦合黏性项。借鉴Biot理论高频拓展方法^[29]，可得包含非牛顿流体(分数阶Maxwell流体模型)的孔隙介质波动方程

$ \left\{\begin{array}{l} \mu_{\mathrm{S}} \nabla^{2} \boldsymbol{u}+\nabla\left[\left(A+\mu_{\mathrm{S}}\right) \nabla \cdot \boldsymbol{u}+G \nabla \cdot \boldsymbol{U}\right]= \\ \rho_{11} \frac{\partial^{2} \boldsymbol{u}}{\partial t^{2}}+\rho_{12} \frac{\partial^{2} \boldsymbol{U}}{\partial t^{2}}-\frac{\eta \phi^{2} F(c)}{\kappa}(\dot{\boldsymbol{U}}-\dot{\boldsymbol{u}}) \\ \nabla(G \nabla \cdot \boldsymbol{u}+R \nabla \cdot \boldsymbol{U})= \\ \rho_{12} \frac{\partial^{2} \boldsymbol{u}}{\partial t^{2}}+\rho_{22} \frac{\partial^{2} \boldsymbol{U}}{\partial t^{2}}+\frac{\eta \phi^{2} F(c)}{\kappa}(\dot{\boldsymbol{U}}-\dot{\boldsymbol{u}}) \end{array}\right. $

(16)

式中：μ_S是固体相剪切模量；ρ₁₁、ρ₁₂、ρ₂₂分别为Biot理论中固体、流体和流固相互作用密度系数；G和R是Biot理论的弹性系数。将散度和旋度算子应用于波动方程固体和流体位移，令ζ₁= $ \nabla $·u、ξ₁= $ \nabla $·U、s= $ \nabla $×u、S= $ \nabla $×U，可得纵波波动方程

$ \left\{\begin{array}{l} \left(A+2 \mu_{s}\right) \nabla^{2} \zeta_{1}+G \nabla^{2} \xi_{1}= \\ \rho_{11} \frac{\partial^{2} \zeta_{1}}{\partial t^{2}}+\rho_{12} \frac{\partial^{2} \xi_{1}}{\partial t^{2}}-\frac{\eta \phi^{2} F(c)}{\kappa} \frac{\partial}{\partial t}\left(\xi_{1}-\zeta_{1}\right) \\ G \nabla^{2} \zeta_{1}+R \nabla^{2} \xi_{1}= \\ \rho_{12} \frac{\partial^{2} \zeta_{1}}{\partial t^{2}}+\rho_{22} \frac{\partial^{2} \xi_{1}}{\partial t^{2}}+\frac{\eta \phi^{2} F(c)}{\kappa} \frac{\partial}{\partial t}\left(\xi_{1}-\zeta_{1}\right) \end{array}\right. $

(17)

横波波动方程为

$ \left\{\begin{array}{l} \mu_{\mathrm{S}} \nabla^{2} \boldsymbol{s}=\rho_{11} \frac{\partial^{2}}{\partial t^{2}} \boldsymbol{s}+\rho_{12} \frac{\partial^{2}}{\partial t^{2}} \boldsymbol{S}-\frac{\eta \phi^{2} F(c)}{\kappa} \frac{\partial}{\partial t}(\boldsymbol{S}-\boldsymbol{s}) \\ \rho_{12} \frac{\partial^{2} \boldsymbol{s}}{\partial t^{2}}+\rho_{22} \frac{\partial^{2} \boldsymbol{S}}{\partial t^{2}}+\frac{\eta \phi^{2} F(c)}{\kappa} \frac{\partial}{\partial t}(\boldsymbol{S}-\boldsymbol{s})=0 \end{array}\right. $

(18)

根据平面波假设，四种波场可表示为ζ₁=C₁×exp[i(ωt-k·x)]、ξ₁=C₂exp[i(ωt-k·x)]、s=T₁exp[i(ωt-k·x)]和S=T₂exp[i(ωt-k·x)]，其中k是波数，C₁、C₂、T₁、T₂是波(矢量)常数。则纵波方程变为

$ \left[\begin{array}{lll} a_{11}\left(\frac{k}{\omega}\right)^{2}+b_{11} & a_{12}\left(\frac{k}{\omega}\right)^{2}+b_{12} \\ a_{21}\left(\frac{k}{\omega}\right)^{2}+b_{21} & a_{22}\left(\frac{k}{\omega}\right)^{2}+b_{22} \end{array}\right]\left[\begin{array}{l} C_{1} \\ C_{2} \end{array}\right]=0 $

(19)

式中：a₁₁=A+2μ_S；b₁₁=-ρ₁₁+iηϕ²F(c)/(κω)；a₁₂= a₂₁=G；b₁₂=b₂₁=-ρ₁₂-iηϕ²F(c)/(κω)；a₂₂=R；b₂₂=-ρ₂₂-iηϕ²F(c)/(κω)。

将波场表达式代入横波方程，经化简可得

$ \left[\begin{array}{l} \left(a_{11}^{\mathrm{S}} \frac{k^{2}}{\omega^{2}}+b_{11}^{\mathrm{S}}\right) \boldsymbol{I} \qquad b_{12}^{\mathrm{S}} \boldsymbol{I} \\ b_{21}^{\mathrm{S}} \boldsymbol{I} \qquad \left(a_{22}^{\mathrm{S}} \frac{k^{2}}{\omega^{2}}+b_{22}^{\mathrm{S}}\right) \boldsymbol{I} \end{array}\right]\left[\begin{array}{l} \boldsymbol{T}_{1} \\ \boldsymbol{T}_{2} \end{array}\right]=0 $

(20)

式中：a₁₁^S=μS；b₁₁^S=-ρ₁₁+iηϕ²F(c)/(κω)；b₁₂^S=-ρ₁₂-iηϕ²F(c)/(κω)；b₂₁^S=-ρ₁₂-iηϕ²F(c)/(κω)；b₂₂^S=-ρ₂₂+iηϕ²F(c)/(κω)。

波动方程存在非零解的条件是系数行列式为零，利用平面波分析方法得到速度频散和衰减。因此，可得纵波方程的解

$ (k / \omega)^{2}=\frac{-B_{\mathrm{P}} \pm \sqrt{B_{\mathrm{P}}^{2}-4 M_{\mathrm{P}} E_{\mathrm{P}}}}{2 M_{\mathrm{P}}} $

(21)

式中

$ M_{\mathrm{P}}=\left(A+2 \mu_{\mathrm{S}}\right) R-G^{2} $

(22)

$ \begin{aligned} B_{\mathrm{P}}=&-\rho_{11} R-\left(A+2 \mu_{\mathrm{S}}\right) \rho_{22}+2 \rho_{12} G+\\ & \frac{\mathrm{i} \eta \phi^{2}}{\kappa \omega} F(c)\left(A+2 \mu_{\mathrm{s}}+2 G+R\right) \end{aligned} $

(23)

$ E_{\mathrm{P}}= \rho_{11} \rho_{22}-\rho_{12}^{2}-\frac{\mathrm{i} \eta \phi^{2}}{\kappa \omega} F(c)\left(\rho_{11}+\rho_{22}+2 \rho_{12}\right) $

(24)

求解上面关于k/ω的方程可以得到纵波复速度 V_P^*=ω/k，进而可得相速度V_P=Re(V_P^*) 和衰减因子Q_P^-1=Im[(V_P^*)²]/Re[(V_P^*)²]。

横波的k/ω满足

$ B_{\mathrm{S}}(k / \omega)^{2}+\varPsi_{\mathrm{S}}=0 $

(25)

式中系数为：B_S=-ρ₂₂μ_S+iF(c)μ_Sηϕ²/(κω)；Ψ_S=ρ₁₁ρ₂₂-ρ₁₂²μ_S-iF(c)(ρ₁₁+ρ₂₂+2ρ₁₂)ηϕ²/(κω)。求解上式可得横波复速度V_S^*=ω/k，则相速度V_S=Re(V_S^*)、衰减因子Q_S^-1=Im[(V_S^*)²]÷Re[(V_S^*)²]。

3 数值算例 3.1 流变参数和分数阶导数对流体黏性耗散系数的影响

溶液的流动特性和黏弹性取决于临界剪切率，在临界剪切率之上可以观察到剪切稀化(Shear Thinning)和剪切稠化(Shear Thickening)现象。这种非牛顿流体行为受Deborah数的控制，该参数通过可测量参数计算获得，包括流体密度ρ_f、黏度η、弛豫时间λ和孔隙通道半径a。应用甘油和CPyCL/NaSal(氯化钠和水杨酸钠)溶液定量计算流—固耦合黏性耗散函数F(c)。两种流体的参数^[22]为：甘油(牛顿流体)的密度为1250kg/m³，黏性为1 Pa·s；CPyCL/NaSal溶液(非牛顿流体)的密度为1050 kg/m³，黏度为60 Pa·s。

甘油的弛豫时间设定为一个很小的值(10^-30s)，代表牛顿流体力学行为，其分数阶导数为(0，1.00)和(0，1.05)时的F(c)函数计算结果如图 1所示。从图 1可以看出，当频率接近0时，甘油的黏性耗散修正函数值等于1，表明低频极限下黏性耗散系数退化到Biot理论的低频情况；当频率逐渐增大时，黏性耗散修正函数趋近于线性增长，实部和虚部逐渐平行于一条渐近线(Asymptotic line)。该算例给出的甘油黏性耗散函数随频率变化规律与Biot理论给出的牛顿流体行为完全一致^[22]，这一结论与Tsiklauri等^[30]的研究吻合。计算结果表明，当甘油弛豫时间趋近于零时，牛顿流体可由fdMaxwell模型获得。

值得指出的是，对于孔隙含有牛顿流体的情况，Biot^[29]给出了高频下孔隙介质波动方程的修正方法，在低频情况下牛顿流体黏性耗散修正函数趋近于1，高频情况下黏性耗散系数随频率呈现线性增长，其渐近线为$ F\left( c \right) = \frac{c}{3}\frac{{1 + {\rm{i}}}}{{\sqrt 2 }}$。对于fdMaxwell流体，根据前面给出的黏性耗散修正函数F(c)定义式可知，当频率接近无穷大时，渐近线满足F(c)=- $ \frac{{{\rm{i}}c}}{{4\mathit{\Delta }}}$。可以看出，非牛顿流体对孔隙介质流固耦合的影响体现在参数Δ中，该参数是流体Deborah数D_e的函数，反映了非牛顿流体处于流体和固体之间的特殊流变学属性；同时，分数阶本构关系的阶数也出现在流固耦合黏性耗散机制中，反映了真实流体本构关系的复杂性。

分数阶导数(与应变延迟有关)对F(c)有显著影响。从图 1可以看出，阶数β增加5%导致F(c)大幅下降，这将显著降低流体与固体之间的黏性作用。因此，分数阶导数的阶数β刻画了应变与应力变化率的关系。

CPyCL/NaSal溶液的弛豫时间设置为1.9s，用于解释其非牛顿流体行为^[22]，其分数阶导数为(1，1.0)和(1，1.5)的F(c)函数计算结果如图 2所示。

首先，当α=1和β=1.0时，模型退化为传统的Maxwell模型。CPyCL/NaSal溶液的F(c)在低频率范围内快速衰减，然后呈现多个峰值，这些峰值反映出CPyCL/NaSal溶液具有近似固体的弹性共振现象。采用α=1和β=1.5说明分数阶导数对流固耦合黏性耗散机制的影响。在低频率状态下仍然可以观察到耗散函数的快速衰减，然而并未出现多个峰值。计算结果表明，分数阶导数应变—应力关系控制了非牛顿流体的共振行为，对流体从耗散态到弹性系统的过渡具有重要影响。

3.2 fdMaxwell流体分数阶导数对纵波速度频散和衰减的影响

采用平均孔隙度为21%的French Vosgian砂岩^[31](表 1)分析纵波速度频散和衰减对黏弹性流体的分数阶导数的依赖关系。

盐水的弛豫时间被设定为一个微小的值(10ns)模拟牛顿流体的行为。

图 3为由Biot理论和Maxwell模型(分数阶导数α=1和β=1)计算的砂岩饱含盐水的波速频散和衰减曲线，可见，无论是地震频带还是超声波频带(小于1MHz)，非牛顿波动方程的纵波波速度与Biot的理论匹配得很好，微小的差异发生在超高频带，这是因为在波动方程中摩擦损耗被重新定义了。该算例表明，在饱含牛顿流体(盐水)的多孔介质中，Biot理论在弹性波传播方面很有效，在此情况下，黏弹性效应可以忽略不计，两种理论表现一致。

由Biot理论和Maxwell模型(分数阶导数α=1和β=1.0)和fdMaxwell模型(分数阶导数α=1和β=1.5)计算饱含CPyCL/NaSal溶液的French Vosgian砂岩的纵波速度频散和衰减曲线如图 4所示。由图可见，当β从1.0增至1.5时，频散和衰减曲线会变得更加平滑，较大的β值压制了黏弹性共振现象。此外，fdMaxwell模型的速度跃升和衰减峰值会向低频方向移动(出现在大约1MHz)，但不像Maxwell模型那样低(大约10kHz)；fdMaxwell模型的纵波速度与Biot理论相同，与Maxwell模型有很大不同。fdMaxwell模型的导数阶数表征了Maxwell模型和Biot理论之间的中间流变动态。分数阶导数的值对fdMaxwell流体饱和的多孔介质弹性波速具有明显的影响。

虽然分数阶导数对表征纯弹性和纯黏性之间的黏弹性行为起重要作用，但通过实验确定分数阶导数的阶数仍是一个挑战。许多研究试图从物理上解释分数导数，例如分形流变模型将聚合物链的变形机理与分数阶微分方程联系起来。但是，宏观应力和微结构变形之间的关系只能部分解决真实的聚合物液体，如何提供合理的分数导数参数拟合真实的黏弹性流体仍然十分困难。虽然目前的理想化模型过分简化了分子水平的细节，但含有fdMaxwell流体的多孔介质波速预测模型仍然提供了对非牛顿流体效应机理的深入认识。

3.3 不同温度下含重油砂岩纵波速度预测

利用含非牛顿流体波动方程对不同温度(T)下含重油(API=9.2)砂岩的波速进行了预测。温度变化范围为20~80℃，骨架密度为2630kg/m³，孔隙度为0.326，渗透率为9D^[32]。研究表明，在10~150℃范围内，骨架体积模量和剪切模量随温度线性降低^[33]，如图 5所示。

重油黏性、密度、剪切模量随温度变化曲线^{[32, 34]}如图 6所示，随着温度的上升，重油黏性和剪切模量下降，黏性系数下降趋势逐渐变缓。随着温度升高，重油作为黏滞流体的动态渗透率明显升高(图 7)，在20~80℃范围内，给出了三个不同频率对应的动态渗透率变化情况，其中25Hz处于地震波频率范围，100kHz和1MHz处于超声波频率范围。从图中可以看出，不同频率的渗透率变化规律基本相同。

在1MHz、100kHz和25Hz频率下模拟了纵波速度随温度变化(随温度的升高呈现下降趋势)，并与实验数据^[32]做了对比(图 8~图 10)。不同频率下，纵波速度随着温度升高逐渐下降，随着频率降低，这种速度下降更剧烈。在实验频率(1MHz)附近，非牛顿流体模型预测速度值更加接近实测数据(图 8)。当频率降到100kHz(图 9)和25Hz(图 10)时，牛顿流体模型和非牛顿流体模型预测速度均低于实验测量值，这是由于速度的频散引起的。计算结果显示，在不同频率下，非牛顿流体波动方程与常规波动方程预测结果均存在明显差异，非牛顿流体本构关系对流固耦合带来新的变化，对孔隙介质波场频散和衰减影响不可忽视。非牛顿流体模型的预测结果在实验频率下与观测数据吻合较好。

4 结论

本文提出了基于分数阶Maxwell流体模型的非牛顿流体孔隙介质波动方程方法。根据重新定义的黏性耗散函数，Biot理论被推广到非牛顿流体饱和情况，能够解释更复杂的流—固耦合效应。利用波动方程定量研究流变学参数对多孔弹性介质波场传播的影响，研究表明弹性波速度频散和衰减不仅取决于密度、弹性模量和孔隙度，还取决于孔隙流体流变参数，如Deborah数和分数阶导数的阶数等参数。理论结果与实验数据对比分析表明，在牛顿和非牛顿流体饱和的砂岩中，波速存在显著差异。同时，分数阶导数本构关系控制着从黏性耗散状态到弹性机制的转变，对储层岩石波场传播特征具有重要影响。但是，不同流体如何选取分数阶参数是一个没有完全解决的问题，仍然需要根据实验或者新理论模型进行深入的研究。