ダイアライザ血液層の圧力・流量分布の解析解を求めてみる

臨床工学の分野の中でも，透析ってそれなりに数理的に遊んでみる余地のある分野かなと思います．物質移動や透析効率など，教科書にも数式のよく出てくる分野です．今回は，中空糸型ダイアライザの長手方向の圧力分布，流量分布について解析解を求めてみようと思います．

ダイアライザとは

腎臓の機能を代行する治療を一般に「透析」と呼び，その治療に必要な人工臓器をダイアライザと言います．ダイアライザにも幾つかの種類がありますが，現在臨床使用されているものの殆どは中空糸型と呼ばれるものです．今回はこの中空糸型ダイアライザ内の流れについてモデル化し，その理論解析を試みてみようと思います．下にその代表的な製品を挙げておきます．

トレライトNV

東レ株式会社が製造販売しているポリスルホン膜の中空糸型ダイアライザ．（画像は東レのWEBサイトから引用）

中空糸型ダイアライザのモデル化

寸法パラメータを規定する

中空糸型ダイアライザの中身はどうなっているでしょうか．

右のカットイメージを見てもらえば分かるように，中空糸型ダイアライザの中身はパイプ状の透析膜が無数に束ねられた構造をしています．その両端は中空糸の開口部を残してポッティングされ膜の内側と外側が隔てられています．この寸法形状を規定するパラメータとしては，中空糸直径 $d$ ，中空糸有効長 $L$ ，中空糸本数 $n$ ，有効膜面積 $A$ が挙げられ，その関係は次の通りです．

(1) $\begin{equation*} \pi dLn = A \end{equation*}$

変数は4つありますが，陽に与えられるパラメータとしては $A, \, L, \, d$ を採用します．

出典：Fresenius Medical Care , FX high-flux and FX low-flux dialyzers

流れる流体の量と圧のパラメータを規定する

ダイアライザの模式図を右に示します．左の赤い部分が血液の入り口，青い部分が出口です．下に付いている2つの突起は透析液の出入り口を表していて，通常は透析液が血液と対向流となるように流します．

今回は中空糸の内側に流れる血液のダイアライザ長手方向に関する圧力分布と流量分布を求めることを目的としています．この分布をダイアライザの入口からの距離 $x$ を変数とした関数として表し，それぞれ圧力 $p(x)$ と流量 $q(x)$ とおきます．この分布はもちろん既に調べられていて，およそ右のグラフのようになると考えられています．しかし，具体的な関数系も含めて記述されているものは少ないのでこれを明らかにしていきます．

各ポート出入り口における圧力と流量については，今後次の文字を使用します．

入口圧： $P_I$
静脈圧（出口圧）： $P_V$
透析液圧： $P_D$

脱血流量： $Q_B$
透析液流量： $Q_D$
除水速度： $Q_S$

この時，関数 $p(x)$ と $q(x)$ の境界条件は次の通りになります．

(2) $\begin{align*} p(0) &= P_I & p(L) &= P_V \\ q(0) &= Q_B & q(L) &= Q_B -Q_S \end{align*}$

中空糸内流れの関係式

中空糸の内側の流れを考えると，中空糸の内側と外側とでの交通はあるものの，局所的にはハーゲン・ポアズイユ流れを仮定できます．つまり，中空糸1本当たりの流量は圧力勾配と中空糸直径の4乗に比例すると考えられます．

(3) $\begin{equation*} \frac{q(x)}{n} \propto \frac{\mathrm{d}p}{\mathrm{d}x} d^4 \end{equation*}$

式(1)を用いて $n$ を用いない形に変形したうえで比例定数を $K_1$ として次の式を立てることができます．

(4) $\begin{equation*} \frac{\mathrm{d}p}{\mathrm{d}x} = - \frac{L K_1}{Ad^3} q(x) \end{equation*}$

層間流れの関係式

中空糸の内側から外側への流体の移動量について考えます．中空糸は円筒状のケースの中を半ば埋め尽くす程に多数固定され，血液と透析液はその内側と外側にそれぞれ接しています．ここではその構造を右の図に示すように簡略化し，血液層と透析液層がひとつの透析膜を介して接してると考えます．入口ポートからの距離が $x$ の微小な領域を切り出して流量のバランス（連続の式）を考えることで微分方程式を立式します．

血液層での流量変化は血液層と透析液層との間の限外濾過によってもたらされると考え，その濾過流量は圧力較差と膜面積に比例すると仮定します．透析液層の圧力は全ての場所で一定の $P_D$ とし，血液層の圧力は $p(x)+\frac{1}{2}\mathrm{d}p$ で代表されるとすると，ダイアライザ全体の膜面積 $A$ と比例定数 $K_2$ を用いて次の関係が成り立ちます．ただし，浸透圧などは無視します．

(5) $\begin{equation*} q(x) = q(x) + \mathrm{d}q + K_2 \left( p(x) + \frac{1}{2}\mathrm{d}p - P_D \right) \frac{A}{L} \mathrm{d}x \end{equation*}$

微小量の積は無視し，整理すると次の式となります．

(6) $\begin{equation*} \frac{\mathrm{d}q}{\mathrm{d}x} = - \frac{ A K_2}{L} \left( p(x) - P_D \right) \end{equation*}$

微分方程式を解く

解くべき連立微分方程式は先に立てた次の式です．

(7) $\begin{align*} \begin{cases} \frac{\mathrm{d}p}{\mathrm{d}x} &= - \frac{L K_1}{Ad^3} q(x) \\ \frac{\mathrm{d}q}{\mathrm{d}x} &= - \frac{ A K_2}{L} \left( p(x) - P_D \right) \end{cases} \end{align*}$

今回はちょっと手抜きします．

WolframAlpha で一般解を求める

ちょっとした式変形や，微分方程式を解いてもらったり，グラフを描いてもらうのに便利なサービスとして WolframAlpha が有名です．これでサクッと一般解を見つけてもらいましょう．解きたい式は独立変数を $x$ とした2つの関数 $p(x)$ および $q(x)$ についての連立一階微分方程式です．WolframAlpha で連立微分方程式を解かせる方法はちょっとググった感じでは無さげだったので，ひとつの二階微分方程式に書き直します．

(8) $\begin{equation*} \frac{\mathrm{d}^2 p}{\mathrm{d}x^2} = \frac{K_1 K_2}{d^3} \left( p(x) - P_D \right) \end{equation*}$

これをこんな感じで解いてもらいます．定数係数は全て簡素化して入力しました．すると次のような結果が得られます．

DSolve[{p''[x] == a (-b + p[x])}, p[x], x]

{p[x] == b + E^(Sqrt[a] x) Subscript[c, 1] + Subscript[c, 2]/E^(Sqrt[a] x)}

定数係数を戻して分かりやすく書くと次のようになります．ただし， $C_1$ および $C_2$ は積分定数で， $K = \sqrt{\frac{K_1 K_2}{d^3}}$ とおいて表記を簡潔にしています．

(9) $\begin{equation*} p(x) = P_D + C_1 e^{Kx} + C_2 e^{-Kx} \end{equation*}$

(10) $\begin{equation*} q(x) = -\frac{Ad^3K}{LK_1} \left( C_1 e^{Kx} - C_2 e^{-Kx} \right) \end{equation*}$

特殊解を求める

積分定数2つに対して境界条件を4つも用意しているので，そのうち2つを利用して特殊解を求めます．

圧力分布の境界条件を利用する場合

境界条件として $p(0) = P_I , p(L) = P_V$ を用いると，積分定数は次のようになり，特殊解が得られます．

(11) $\begin{equation*} C_1 = \frac{1}{e^{2KL}-1} \left\{ (P_V - P_D)e^{KL} - (P_I -P_D) \right\} \end{equation*}$

(12) $\begin{equation*} C_2 = \frac{e^{KL}}{e^{2KL}-1} \left\{ (P_I-P_D)e^{KL} - (P_V-P_D) \right\} \end{equation*}$

(13) $\begin{equation*} p(x) = P_D + \frac{1}{e^{2KL}-1} \left[ \left\{ (P_V - P_D)e^{KL} - (P_I -P_D) \right\}e^{Kx} + e^{KL} \left\{ (P_I-P_D)e^{KL} - (P_V-P_D) \right\}e^{-Kx} \right] \end{equation*}$

(14) $\begin{equation*} q(x) = -\frac{Ad^3K}{LK_1}\frac{1}{e^{2KL}-1} \left[ \left\{ (P_V - P_D)e^{KL} - (P_I -P_D) \right\}e^{Kx} - e^{KL} \left\{ (P_I-P_D)e^{KL} - (P_V-P_D) \right\}e^{-Kx} \right] \end{equation*}$

ちなみに，このようにして WolframAlpha で一気に特殊解を得ることもできます．

DSolve[{p''[x] == a (-b + p[x]), p[0] == c, p[L] == d}, p[x], x]

{p[x] == (b (-1 + E^(Sqrt[a] L)) (E^(Sqrt[a] L) - E^(Sqrt[a] x)) (-1 + E^(Sqrt[a] x)) + c (E^(2 Sqrt[a] L) - E^(2 Sqrt[a] x)) + d E^(Sqrt[a] L) (-1 + E^(2 Sqrt[a] x)))/(E^(Sqrt[a] x) (-1 + E^(2 Sqrt[a] L)))}

(15) $\begin{equation*} p(x) = \frac{e^{-Kx}}{e^{2KL}-1} \left\{ P_D \left( e^{KL} - 1 \right) \left( e^{Kx} - 1 \right) \left( e^{KL} - e^{Kx} \right) + P_I \left( e^{2KL} - e^{2Kx} \right) + P_V e^{KL} \left( e^{2Kx} - 1 \right) \right\} \end{equation*}$

一見異なる式に見えますが，式変形により同値であることが確認できます．

流量分布の境界条件を利用する場合

境界条件として $q(0) = Q_B , q(L) = Q_B -Q_S$ を用いると，積分定数は次のようになり，特殊解が得られます．

(16) $\begin{equation*} C_1 = \frac{LK_1}{Ad^3K} \frac{1}{e^{2KL}-1} \left\{ Q_B - (Q_B - Q_S)e^{KL} \right\} \end{equation*}$

(17) $\begin{equation*} C_2 = - \frac{LK_1}{Ad^3K} \frac{e^{KL}}{e^{2KL}-1} \left\{ (Q_B - Q_S) - Q_B e^{KL} \right\} \end{equation*}$

(18) $\begin{equation*} p(x) = P_D + \frac{LK_1}{Ad^3K} \frac{1}{e^{2KL}-1} \left[ \left\{ Q_B - (Q_B - Q_S)e^{KL} \right\}e^{Kx} - e^{KL} \left\{ (Q_B - Q_S) - Q_B e^{KL} \right\}e^{-Kx} \right] \end{equation*}$

(19) $\begin{equation*} q(x) = - \frac{1}{e^{2KL}-1} \left[ \left\{ Q_B - (Q_B - Q_S)e^{KL} \right\}e^{Kx} +e^{KL} \left\{ (Q_B - Q_S) - Q_B e^{KL} \right\}e^{-Kx} \right] \end{equation*}$

$P_I$ や $P_V$ に対して $Q_B$ や $Q_S$ は治療条件として陽に与えられるものなので，こちらの方法の方がより自然な表現となります．

圧力境界条件解と流量境界条件解の関係性

圧力分布の境界条件を利用した場合でも，流量分布の境界条件を利用した場合でも，それにより求まる積分定数 $C_1 , C_2$ は式 (11)と式 (16)，式 (12)と式 (17) でそれぞれ同じ値となるはずです．

(20) $\begin{equation*} \frac{1}{e^{2KL}-1} \left\{ (P_V - P_D)e^{KL} - (P_I -P_D) \right\} = \frac{LK_1}{Ad^3K} \frac{1}{e^{2KL}-1} \left\{ Q_B - (Q_B - Q_S)e^{KL} \right\} \end{equation*}$

(21) $\begin{equation*} \frac{e^{KL}}{e^{2KL}-1} \left\{ (P_I-P_D)e^{KL} - (P_V-P_D) \right\} = - \frac{LK_1}{Ad^3K} \frac{e^{KL}}{e^{2KL}-1} \left\{ (Q_B - Q_S) - Q_B e^{KL} \right\} \end{equation*}$

この2式を整理すると次の2式となります．

(22) $\begin{equation*} P_I = P_D + \frac{LK_1}{Ad^3K} \frac{1}{e^{2KL}-1} \left\{ Q_B e^{2KL} - 2(Q_B - Q_S)e^{KL} + Q_B \right\} \end{equation*}$

(23) $\begin{equation*} P_V = P_D + \frac{LK_1}{Ad^3K} \frac{1}{e^{2KL}-1} \left\{ -(Q_B - Q_S) e^{2KL} + 2Q_B e^{KL} - (Q_B - Q_S) \right\} \end{equation*}$

$K, \, K_1, \, L, \, A, \, d$ はダイアライザの材料定数として， $Q_B, \, Q_S$ は治療条件として， $P_V$ は運転条件として given なパラメータであると考えるとそれらによって $P_D, \, P_I$ が計算可能であることが分かります．

陽に与えられるパラメータのみで表現する

式(23)を利用して式(18)を given なパラメータのみで表すと次のようになります．ただし， $Q_S = \alpha Q_B$ としています．

(24) $\begin{equation*} p(x) = P_V + \frac{LK_1}{Ad^3K} \frac{Q_B}{e^{2KL}-1} \Bigl[ \left\{ 1 - (1 - \alpha)e^{KL} \right\} e^{Kx} + \left\{ e^{KL} - (1 - \alpha) \right\} e^{K(L-x)} + (1 - \alpha)(e^{2KL}+1) -2e^{KL} \Bigr] \end{equation*}$

また，式(19)も $Q_S$ の代わりに $\alpha$ を用いると次のようになります．

(25) $\begin{equation*} q(x) = \frac{Q_B}{e^{2KL}-1} \left[ \left\{ (1 - \alpha)e^{KL} -1 \right\}e^{Kx} +e^{KL} \left\{ e^{KL} - (1 - \alpha) \right\}e^{-Kx} \right] \end{equation*}$

除水なしの場合

除水なし，すなわち $\alpha = 0$ の場合は次のようになります．

(26) $\begin{equation*} \left. p(x) \right|_{\alpha =0} = P_V + \frac{LK_1}{Ad^3K} \frac{Q_B}{e^{KL}+1} \left( - e^{Kx} + e^{K(L-x)} + e^{KL} -1 \right) \end{equation*}$

(27) $\begin{equation*} \left. q(x) \right|_{\alpha =0} = \frac{Q_B}{e^{KL}+1} \left( e^{Kx} + e^{K(L-x)} \right) \end{equation*}$

長くなってきたので，今回はここまで！