Navier-Stokes 方程式（2）：穩態流速 vs 脈態流速（一）

文： 張國柱（Chang，Kuo-Chu）

         台大名譽教授

 

日期：2021/12/06

 

動脈血流是心臟-血管交互作用的結果，這牽涉到心臟所作的功率與動脈所消耗的功率彼此之間的平衡關係。動脈血流在循環系統中扮演非常重要的角色，而 Navier-Stokes 方程式則是探索動脈血流的唯一物理數學基礎。以物理的觀點言，Navier-Stokes 方程式是由一組特定的偏微分方程所組成，用來描述黏滯流體物質（viscous fluid substance）的運動行為。

 

Navier-Stokes 方程式

假設流體為不可壓縮（密度 ρ = 常數）之牛頓液（黏滯度 η = 常數），此時流體之 Navier-Stokes 方程式為：

 

ρ(Dv_i/Dt) = X_i－(∂p/∂x_i) + η∇²v_i    i,j = 1,2,3 ⋯⋯⋯⋯⋯⋯⋯⋯⋯（1）

 

其中

 

D = (∂/∂t+v_j∂/∂x_j)

 

∂/∂x_j = v₁∂/∂x₁ + v₂∂/∂x₂ + v₃∂/∂x₃

 

 

∇² = ∂²/∂x₁² + ∂²/∂x₂² +∂²/∂x₃²

 

v_i = 流體流速

 

∂v_i/∂t = 瞬時加速度（transient acceleration）

 

v_j∂v_i/∂x_j = 位變加速度（convective acceleration）

 

ρ(Dv_i/Dt) = 流體慣性力（inertial force）

 

X_i = 重力（gravitational force）

 

-∂p/∂x_i = 驅動力（driving pressure）

 

η∇²v_i = 黏滯力（viscous force）

 

假設圓形截面之圓柱管的長度為 L、半徑為 a，由於是圓柱管，此後之方程式就以圓柱坐標（x,r,θ）表示，所對應之流體流速為（u,v,w）。

 

（一）硬管的穩態和脈態流速（steady and pulsatile flow in a rigid tube）

想像充滿液體的圓柱管是由無窮多個同軸圓柱殼所構成的集合。假設圓柱管為硬管，將它放在同一水平面上，那麼任何一處的流體元素（fluid element）所受的重力完全一樣，因此重力可以忽略不計。由於是硬管，所以 v = 0。假設沒有擾流（w = 0），那麼對發展完整的流速（fully developed flow rate）而言，u僅是半徑與時間的函數，也就是：u(r,t)，v = 0，w = 0。此時之Navier-Stokes 方程式如下：

 

ρ(∂u/∂t)+ ∂p/∂x = η[(∂²u/∂r²)+ (∂u/r∂r)] ⋯⋯⋯⋯⋯⋯⋯⋯⋯（3）

 

公式（3）顯示：沿著管子而變的參數為壓力而非流速，也就是 p = p(x,t)；流速僅是半徑的函數，也就是 u(r,t)。方程式（3）可用於硬管的穩態流速與脈態流速分析。

 

【1】公式（3）之穩態解（steady-state solution）：穩態流速

假設硬管裡的流體流速不隨時間而改變，也就是說p = p_s(x)，u = u_s(r)。那麼公式（3）便成為：

 

dp_s/dx = η[(d²u_s/dr²)+ (du_s/rdr)] ⋯⋯⋯⋯⋯⋯⋯⋯⋯（4）

 

公式（4）的重要特性是此方程式與流體密度（ρ）無關。令 k_s = 壓力梯度 = dp_s/dx < 0。配合兩個邊界條件：u_s(a) = 0；|u_s(0)| < ∞，解方程式（4）可得具拋物線之穩態層流流速：

 

u(r) = (k_s/4η) (r² – a²)⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（5）

 

對體積流速（volumetric flow rate，q_s）而言，

 

q_s = ∫u(r) × 2πrdr  ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（6）

 

由於圓柱管具對稱性，因此從 0 到 a 作積分可得：

 

q_s = –k_sπa⁴/8η ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（7）

 

公式（7）即是所謂的 Poiseuille 體積流速。

 

【2】公式（3）之脈態解（pulsatile-state solution）：脈態流速

假設硬管裡的流體流速隨時間而改變，也就是 p(x,t) = p_s(x) + p_φ(x,t)；u(r,t) = u_s(r) + u_φ(r,t)。將之帶入公式（3）可得脈態分量如下：

 

ρ(∂u_φ/∂t)+ ∂p_φ/∂x = η[(∂²u_φ/∂r²)+ (∂u_φ/r∂r)] ⋯⋯⋯⋯⋯⋯⋯⋯⋯（8）

 

由於 p(x,t) = p_s(x) + p_φ(x,t)，因此 ∂p/∂x = dp_s/dx + ∂p_φ/∂x。令 k(t) = ∂p/∂x；k_s = dp_s/dx；k_φ(t) =∂p_φ/∂x。將 k_φ(t) = ∂p_φ/∂x 帶入公式（8）可得：

 

η[(∂²u_φ/∂r²)+ (∂u_φ/r∂r)] – ρ(∂u_φ/∂t) = k_φ(t) ⋯⋯⋯⋯⋯⋯⋯⋯⋯（9）

 

令 k_φ(t) = k_se^iωt；u_φ(r,t) = U_φ(r)e^iωt，並將之帶入公式（9）可得：

 

[(d²U_φ/dr²)+ (dU_φ/rdr)] – (iΩ²/a²)U_φ = k_s/η ⋯⋯⋯⋯⋯⋯⋯⋯⋯（10）

 

其中 Ω = a(ρω/η)^1/2，Ω 又稱為 Womersley 參數。因為 U_φ(r) 是複變函數，故 k_φ(t) 與 u_φ(r,t) 之相位並不相同。

 

公式（10）為貝索方程式（Bessel equation）的型態，其邊界條件為：U_φ(a) = 0；| U_φ(0)| < ∞。因此貝索方程式之解為

 

U_φ(r) = (ik_sa²/ηΩ²){1– [J₀(ζ)/J₀(Λ)]} ⋯⋯⋯⋯⋯⋯⋯⋯⋯（11）

 

其中 ζ(r) = Λ(r/a)，Λ = Ω[(i–1)/√2]，J₀ 為零階貝索函數。

 

因為 u_φ(r,t) = U_φ(r)e^iωt，所以

 

u_φ(r,t) =(ik_sa²/ηΩ²){1– [J₀(ζ)/J₀(Λ)]}e^iωt ⋯⋯⋯⋯⋯⋯⋯（12）

 

對脈態體積流速（pulsatile volumetric flow rate，q_p）而言，

 

q_p(t) = ∫u_φ(r,t) × 2πrdr  ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（13）

 

由於圓柱管具對稱性，因此從 0 到 a 作積分可得：

 

q_p(t) = (ik_sπa⁴/ηΩ²){1– [2J₁(Λ)/ΛJ₀(Λ)]}e^iωt ⋯⋯⋯⋯⋯⋯⋯⋯（14）

 

J₁ 為一階貝索函數。將公式（14）對公式（7）作常規化可得

 

q_p(t)/q_s = (–8/Λ²) {1– [2J₁(Λ)/ΛJ₀(Λ)]}e^iωt ⋯⋯⋯⋯⋯⋯⋯⋯⋯（15）

 

令 g = 2J₁(Λ)/ΛJ₀(Λ)，可將公式（15）改寫如下，

 

q_p(t)/q_s = (–8/Λ²) (1 – g) e^iωt ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（16）

 

注意：將 q_p(t) 作一個週期（T = 2π/ω）的積分，所得之結果為零。因為 q(t) = q_s + q_p(t)，所以一個週期的積分 (1/T) ∫q(t)dt = q_s。這個結果證實了硬管的脈態流速僅是流體的前、後移動（back and forth movements），並沒有實質的流體流向任何方向的淨流（net flow）。

 

（二）彈性管的脈態流速（pulsatile flow in an elastic tube）

假設管子夠長（a/L << 1）、管壁夠薄（h/a << 1）而且僅有些微的彈性（u_m/c₀ << 1），其中 h = 管壁厚度；u_m = 流速在 x 軸之平均值；c₀  = 無黏滯流速（inviscid flow）之波速 = (Eh/2ρr)^1/2，E = 楊氏係數（Young’s modulus）。由於是彈性管，因此 v ≠ 0。在沒有擾流（w = 0），且流體流速隨時間而改變時，u 和 v 不僅是半徑和時間的函數，而且與距離有關：u(x,r,t)，v(x,r,t)。此時之 Navier-Stokes方程式為：

 

ρ(∂u/∂t)+ ∂p/∂x = η[(∂²u/∂x²)+ (∂u/r∂r)] ⋯⋯⋯⋯⋯⋯⋯⋯（17）

 

ρ(∂v/∂t)+ ∂p/∂r = η[(∂²v/∂r²)+ (∂v/r∂r)－(v/r²)] ⋯⋯⋯⋯⋯（18）

 

公式（18）為彈性管之徑向運動方程式。此外流體之連續性方程式（continuity equation）為：

 

∂u/∂x + ∂v/∂r + v/r = 0 ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（19）

 

令 p(x,r,t) = P(r)e^iω(t-x/c)、u(x,r,t) = U(r)e^iω(t-x/c)、v(x,r,t) = v(r)e^iω(t-x/c)，並將之帶入公式（17～19），可得型態為貝索方程式之數學式：

 

(d²U/dr²) + (dU/rdr) – (iωρ/η)U = – iωp/ηc ⋯⋯⋯⋯⋯⋯⋯⋯⋯（20）

 

(d²V/dr²) + (dV/rdr) – [(1/r²) + (iωρ/η)]V = (1/η)(dP/dr) ⋯⋯⋯（21）

 

(dV/dr) + (V/r) – (iω/c)U = 0  ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（22）

 

配合邊界條件可得方程式之解：

 

U(r) = AJ₀(ζ) + B(1/ρc) ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（23）

 

V(r) = A(iωa/cΛ)J₁(ζ) + B(iωr/2ρc²) ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（24）

 

常數 A 和 B 由流體流速和管壁速度在 r = a 處相匹配時所決定。經由相匹配的流體壓力（fluid pressure）與剪應力（shear stress）共同作用於管壁的結果，可得：

 

A = {1/[ρcJ₀(Λ)]}GB

 

B = P(r)

 

G = [2 + z(2σ – 1)]/[z(2σ – g)]

 

z = E_σh/ρac²

 

E_σ = E/(1 – σ²)

 

E = Young’s modulus

 

σ = Poisson’s ratio

 

c = [2/(1 – σ²)z]^1/2c₀ = 黏滯流（viscous flow）之波速

 

g = 2J₁(Λ)/ΛJ₀(Λ)

 

 

將 A 帶入公式（23）可得：

 

U(r) = (B/ρc) {1 – G[J₀(ζ)/J₀(Λ)]} ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（25）

 

因為 u(x,r,t) = U(r)e^iω(t-x/c)，所以

 

u(x,r,t) = (B/ρc) {1 – G[J₀(ζ)/J₀(Λ)]}e^iω(t-x/c) ⋯⋯⋯⋯⋯⋯⋯（26）

 

由於 p(x,r,t) = P(r)e^iω(t-x/c) = Be^iω(t-x/c)，導致 ∂p/∂x = – (iω/c)Be^iω(t-x/c)。令 k_s = –(iω/c)B，可得 B= (ic/ω)k_s。將 B = (ic/ω)k_s 帶入公式（26）

 

[u(x,r,t)/u_sa] = (–4/Λ²) {1 – G[J₀(ζ)/J₀(Λ)]}e^iω(t-x/c) ⋯⋯⋯⋯⋯⋯（27）

 

其中 u_sa = – (k_sa²/4η)。

 

將 A 帶入公式（24）可得：

 

[v(x,r,t)/u_sa] = (2aω/iΛ²c) {(r/a) – G[2J₁(ζ)/ΛJ₀(Λ)]}e^iω(t-x/c) ⋯⋯（28）

 

因為 ζ(r) = Λ(r/a)，因此在管壁 r = a 處，ζ = Λ。

 

[v(x,a,t)/u_sa] = (2aω/iΛ²c) (1 – Gg) e^iω(t-x/c) ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（29）

 

假設管壁的徑向運動很小，對公式（27）從 0 到 a 作積分可得：

 

q_p(x,t)/q_s = (–8/Λ²) (1 – Gg) e^iω(t-x/c) ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（30）

 

注意：彈性因子 G 為複變函數，由頻率所決定。彈性管脈態流速與硬管脈態流速之差異在於彈性因子 G 和波速 c。由於彈性管的 v(x,r,t) 並不為零，也就是彈性管具有管壁的徑向位移，這將使得流體比較容易在彈性管內流動。

 

有關穩態流速與脈態流速進一步的討論，請參考【Navier-Stokes 方程式（3）：穩態流速 vs 脈態流速（二）】