Navier-Stokes 方程式（1）：動脈循環之血壓-血流關係

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

         台大名譽教授

 

日期：2021/09/15

 

十九世紀中葉 Navier 與 Stokes 兩位學者基於動量守恆原則，建立了描述黏滯流體物質（viscous fluid substance）運動行為的一組偏微分方程，稱之為 Navier-Stokes 方程式。

 

廣義的 Navier-Stokes 方程式

假設流體為不可壓縮（密度 ρ = 常數）之非牛頓液（黏滯度 η ≠ 常數），那麼流體之本構方程式（constitutive equation）或廣義的 Navier-Stokes 方程式為：

 

ρ(Dv_i/Dt) = X_i－(∂P/∂x_i) + η∇²v_i + (∂η/∂x_j)[(∂v_j/∂x_i + ∂v_i/∂x_j)]

 

i,j = 1,2,3 ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（1）

 

其中

 

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

 

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

 

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

 

v_i = 流體流速（fluid velocity）

 

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

 

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

 

F_a = ρ(Dv_i/Dt) = 流體受力的合力，也就是慣性力（inertial force）

 

F_g = X_i = 重力（gravitational force）

 

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

 

F_visc = η∇²v_i + (∂η/∂x_j)[(∂v_j/∂x_i + ∂v_i/∂x_j)] = 黏滯力（viscous force）

 

方程式（1）非常複雜且具有非線性項的黏滯力，沒有解析解（analytic solution）。然而在若干條件的假設下，可將廣義的 Navier-Stokes 方程式線性化，以利血行力學（hemodynamics）的探討。

 

線性化的 Navier-Stokes 方程式

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

 

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

 

其中

 

F_visc = η∇²v_i

 

公式（2）即為線性化的 Navier-Stokes 方程式，也是 1955 年 Womersley 用來探討動脈血流（blood flows in arteries）之基礎。

 

（一）硬管的穩態血壓-血流關係：血管阻力（vascular resistance）

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

 

假設圓柱管為硬管（rigid tube），硬管裡的流體具有穩態層流。若將硬管放在同一水平面上，任何一處的流體元素（fluid element）所受的重力完全一樣，可以忽略不計，因此作用在流體元素的慣性力便等於驅動力與黏滯力之和。由於沒有擾流且流體流速不隨時間而改變，對發展完整的血流（fully developed flow）而言，u僅是半徑的函數而與時間無關，也就是：u(r)，v = 0，w = 0。此時任一流體元素所受的驅動力與黏滯力，大小相等、方向相反，因此流體元素所受的慣性力為零，也就是F_a = F_p + F_visc = 0。

 

流體之剪應力（shear stress，S）為單位面積（A = 2πrL）所受之黏滯摩擦力（F_visc），也就是S= F_visc/A。本質上黏滯性（viscosity，η）是阻礙流體變形的速率（rate of deformation），為流體剪應力對速度梯度（velocity of gradient，du/dr）之比值，也就是：

 

 η = S/(du/dr) ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（3）

 

想像圓柱管是由無窮多個同軸圓柱殼所構成的集合，因此在半徑為 r 之圓柱殼，流體元素所受的黏滯阻力為：

 

F_visc = A × S = 2πrLη(du/dr)

 

而推動流體元素前進的趨動力為：

 

F_p = πr²(P₁ – P₂)

 

P₁ = 硬管入口端的壓力，P₂ = 硬管末端的壓力，P₁ > P₂。

 

因為F_a = F_p + F_visc = 0，所以 F_p = －F_visc，也就是

 

du/dr = －r(P₁ – P₂)/ 2ηL ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（4）

 

解微分方程（4），配合邊界條件：u(R) = 0，u(0) < ∞，可得 u(r) 如下：

 

u(r) = (P1 – P2)(R2 – r2)/4ηL ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（5）

由公式（5）可知：u(r) 之速度輪廓（velocity profile）具拋物線特性；u(0) = u_max，u(R) = 0，u_mean = u_R/√2 = u_max/2。

 

對體積流量率（volume flow rate，Q）而言，

 

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

 

由於圓柱管具對稱性，因此從 0 到R 積分公式（6）可得：

 

Q = πR4(P1 – P2)/8ηL ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（7）

公式（7）即是所謂的 Poiseuille 定律。

 

Poiseuille 定律提供一個很重要的觀念，那便是週邊血管阻力（peripheral vascular resistance，R_p）：

 

R_p = (P₁ – P₂)/Q = 8ηL/πR⁴  ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（8）

 

注意：Poiseuille 定律僅適用於（1）具圓形截面積之圓柱硬管，其流體特性為穩態層流，（2）僅表達平均血壓與平均血流之關係（血管阻力），也就是穩態的血壓-血流關係。

 

根據流體黏滯性的定義：η = S/(du/dr)，稍作調整可得剪應力S = η(du/dr)，由此可知管壁阻礙血液流動之剪應力為：S|_r=R = η(du/dr)|_r=R = –4ηQ/πR³。由於管壁對血流所造成的剪應力與血流加諸於管壁的剪應力（τ_s）大小相等、方向相反，故

 

τ_s = –S|_r=R = 4ηQ/πR³ = R(P₁ – P₂)/2L  ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（9）

 

對高血壓之主動脈而言，高剪應力容易對管壁內皮細胞造成傷害而產生破口，這一破口在血流的沖襲下容易導致主動脈剝離（dissecting aneurysm），對生命造成威脅。

 

（二）彈性管的脈態血壓-血流關係：血管阻抗（vascular impedance）

假設管子夠長（R/L << 1）、管壁夠薄（h/R << 1，h = 管壁厚度）而且僅有些微的彈性（u_m/c₀<< 1，u_m = x 軸之平均流速，c₀ = 脈波傳播速度），那麼在沒有擾流（w = 0）、但流體流速隨時間而改變時，u 和 v 不僅是半徑和時間的函數，而且與距離有關：u(x,r,t)，v(x,r,t)。此時Navier-Stokes方程式可表為：

 

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

 

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

 

此外流體之連續性方程式（continuity equation）為：

 

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

 

由於彈性管最重要的波反射效應可用 x 軸之血壓與血流的變化表示，並不需要血流的全貌，因此一維方法（one-dimensional method）便足以表達流體流量之特性，所以公式（11）可以略去而不加以考慮。一維方法的理論即是傳輸線理論（transmission line theorem）。

 

公式（10）、公式（12）配合複雜的邊界條件可推導出下列公式：

 

∂Q/∂t + (A/ρ)(∂P/∂x) = 0 ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（13）

 

∂Q/∂x + (A/ρc₀)(∂P/∂t) = 0  ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（14）

 

其中A = 截面積，c₀ = (Eh/ρ2R)^1/2，E = 楊氏係數（Young’s modules）。將公式（13）對 x 作偏微、將公式（14）對 t 作偏微之後再經調整可得：

 

∂²P/∂t² = c₀²(∂²P/∂x²) ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（15）

 

∂²Q/∂t² = c₀²(∂²Q/∂x²) ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯ （16）

 

公式（15）、公式（16）即是所謂的一維波方程式（one-dimensional wave equations）。這意味著血壓和血流是由相同的波方程式所調控，且其傳播速度相同，但並不代表血壓波和血流波具有相同的相位（in phase），相反地，血壓波和血流波彼此互為反相（out of phase）。

 

波方程式之解可使用變數分離法為之，也就是 P(x,t) = P_x(x)P_t(t)。假設管子入口處（x = 0）之趨動壓為：P(0,t) = P₀e^iωt = P_x(0)P_t(t)，ω = 2πf = 角頻率。因為 P_x(0) = P₀；P_t(t) = e^iωt，所以P(x,t) = P_x(x) e^iωt。此時公式（15）便成為：

 

d²P_x/dx² + (ω²/c₀²) P_x = 0 ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯ （17）

 

方程式（17）之通解為：

 

P(x,t) = Be^iω(t-x/c0) + Ce^iω(t+x/c0) ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（18）

 

其中

 

Be^iω(t-x/c0) = P_f = 前進血壓波（forward pressure wave）

 

Ce^iω(t+x/c0) = P_b = 反射血壓波（backward pressure wave）

 

將入口處之血壓波帶入公式（18）可得 B = P₀，因此P(x,t) = P₀e^iω(t-x/c0)，稍作整理可得 p(x,t) = e^iω(t-x/c0)，其中 p = P/P₀。

 

令初始的前進波為：p_f = e^iω(t-x/c0)，而反射波為：p_b = Ce^iω(t+x/c0)。在管子末端（x = L），波反射係數（R_f）之定義為：

 

R_f = p_b(L,t)/p_f(L,t) ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（19）

 

由此可解得C = R_fe^-2iωL/c0，因此 p_b(x,t) = R_fe^{-iω(t+x/c0-2L/c0)}。由於管子任何一處之血壓波為前進波與反射波之和，因此

 

p(x,t) = p_f(x,t) + p_b(x,t) = e^iω(t-x/c0) + R_fe^{-iω(t+x/c0-2L/c0)}  ⋯⋯⋯⋯（20）

 

因為起始之前進血壓波為：P_f(x,t) = P₀e^iω(t-x/c0)，其所對應之前進血流波為：Q_f(x,t) = Be^iω(t-x/c0)，其中 B 為常數。由公式（13）、公式（14）可知：∂Q_f/∂t + (A/ρ)(∂P_f/∂x) = 0，∂Q_f/∂x +(A/ρc₀²)(∂P_f/∂t) = 0。由此可得 B = (A/ρc₀)P₀，所以

 

Q_f(x,t) = Be^iω(t-x/c0) = (A/ρc₀)P₀ e^iω(t-x/c0) = (A/ρc₀)P_f(x,t) ⋯⋯⋯（21）

 

同理，反射血壓波為：P_b(x,t) = R_fP₀e^{iω(t+x/c0-2L/c0)}，其所對應之反射血流波為：Q_b(x,t) = C e^{iω(t+x/c0-2L/c0)}，其中 C 為常數。將之帶入公式（13）、公式（14）可得 C = －(A/ρc₀)P₀，所以

 

Q_b(x,t) = －(A/ρc₀)P_b(x,t) ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（22）

 

【小小整理：

 

P(x,t) = P_f(x,t) + P_b(x,t)，Q(x,t) = Q_f(x,t) + Q_b(x,t)；

 

P_f(x,t) = P₀ e^iω(t-x/c0)，Q_f(x,t) = (A/ρc₀)P_f(x,t)；

 

P_b(x,t) = R_fP₀e^{iω(t+x/c0-2L/c0)}，Q_b(x,t) = －(A/ρc₀)P_b(x,t)；

 

其中 A = 截面積，c₀ = 波速。】

 

因此彈性管之脈態血壓-血流關係如下：

（1）特徵阻抗（characteristic impedance，Z_c）

根據定義：在沒有反射波存在的情況下，血壓-血流關係便是動脈管的特徵阻抗，也就是：

 

Z_c = P_f(x,t)/ Q_f(x,t) = ρc₀/A ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（23）

 

Z_c 受血管截面積與波速所調控。在截面積不變的情況下，特徵阻抗僅隨波速而改變，此時特徵阻抗可作為動脈管硬化的指標：Z_c愈大，血管的硬化程度就愈高，反之亦然。

 

（2）血管輸入阻抗（vascular input impedance，Z_i）

根據定義：血管輸入阻抗為血管入口處（x = 0）之血壓對血流的比值，也就是：

 

Z_i = [P(x,t)/ Q(x,t)]|_x=0  ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（24）

 

由於

 

P(x,t) = P_f(x,t) + P_b(x,t) = P₀ e^iω(t-x/c0) + R_fP₀e^{iω(t+x/c0-2L/c0)} ⋯⋯⋯⋯⋯⋯（25）

 

Q(x,t) = Q_f(x,t) + (Q_b x,t) = (1/Z_c){ P₀ e^iω(t-x/c0)－R_fP₀e^{iω(t+x/c0-2L/c0)}} ⋯（26）

 

將公式（25）、公式（26）帶入公式（24）可得：

 

Z_i = Z_c{[1+R_fe^-iω2L/c0]/ [1－R_fe^-iω2L/c0]} ⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯（27）

 

【小小整理：

 

R_f = P_b(L,t)/P_f(L,t)；

 

Z_c = P_f(x,t)/ Q_f(x,t) = ρc₀/A；

 

Z_i = [P(x,t)/ Q(x,t)]|_x=0 = Z_c{[1+R_fe^-iω2L/c0]/ [1－R_fe^-iω2L/c0]} 】

 

結語

體循環動脈系統之組成包含：（1）具緩衝、平滑心臟間竭性射血功能的主動脈和大動脈，也就是彈性血管（elastic arteries），（2）具灌流組織、器官功能的小動脈和細動脈，也就是阻性血管（resistance vessels）。主動脈和大動脈的物理性質可用特徵阻抗、波反射之時序與振幅加以描述；小動脈和細動脈的物理性質則表現在週邊血管阻力。無論是硬管的穩態血壓-血流關係或是彈性管的脈態血壓-血流關係都可由 Navier-Stokes 方程式推導之。

 

附錄

Navier-Stokes 方程式的其他血行力學應用

心臟瓣膜閉合的本質是流體的減速（deceleration）

在沒有擾流（turbulent flow）而且重力與黏滯力小到可以忽略不計時，公式（2）可寫成：

 

ρ(Dv_i/Dt) = －(∂P/∂x_i) 

 

若加速度為負（也就是減速）時，Dv_i/Dt < 0，那麼 ∂P/∂x_i > 0，也就是當 x_i2 趨近於 x_i1 時，P₂會大於 P₁，這意味著「流體減速時，在流速的方向上，遠端的壓力隨之增加。」

  

舒張初期，由於左心房壓（P₁）大於左心室壓（P₂），因此僧帽瓣（mitral valve）打開，使得血液由左心房流入左心室。但舒張末了時，血流的減速所造成的壓力梯度，使得左心室壓（P₂）大於左心房壓（P₁），導致瓣膜的關閉。對運作正常的心臟而言，瓣膜的關閉並無伴隨逆流（backward flow）或回流（regurgitation），因此瓣膜關閉的本質是血流的減速而非血液的逆流。