圧損計算の基本について

バルブ、配管、蒸留塔の基本的な圧損計算についてまとめていきます。

・バルブ、オリフィスの圧損計算(亜音速、超音速)

・安全弁の圧損計算(亜音速、超音速)

・安全弁のサイジング(ガス、液、2相流)

・配管の圧損計算(亜音速、超音速、2相流)

・偏流解析

・拡大管/縮小管、曲がり配管の圧損計算

・蒸留塔の圧損計算(トレイ塔、充填塔)

バルブの圧損計算

バルブCV、バルブ特性について

CV値の定義

Cv値はバルブ容量を表す係数で、60°Fの水をCv[gal/min]流した時の圧損が1[lbg/in2]となります。

Cv値は、K値と以下の関係があります。(K値については、配管の圧損計算->代表的なFitting類のK値をご参照ください)

$$C_v=\frac{29.9 \cdot d^2}{\sqrt K}$$

なお、dはバルブポート径[inch]になります。

非圧縮性流体のCVサイジング(乱流域)

・亜音速域のサイジング : \(\Delta P \geq F_L^2(P_1-F_FP_V)\)

$$C=\frac{Q}{N_1} \cdot \sqrt{\frac{\rho_{in} / \rho_0}{\Delta P}}$$

・超音速域のサイジング

$$C=\frac{Q}{N_1 F_L} \cdot \sqrt{\frac{\rho_{in} / \rho_0}{P_{in} - F_F P_v}}$$

※\(N_1 : 8.65E^{-2}[-] , Q[m^3/hr], P[kPa]\)

※\(\rho_0 \): 15C, 1atmにおける水の密度

※\(F_L \): Liquid Pressure Recovery Factor (グローブ弁で0.8~0.9程度)

※\(F_F \): Liquid Critical Pressure ratio factor (\(0.96-0.28 \sqrt {\frac{P_v}{P_c}}\)で近似可能)

圧縮性流体のCVサイジング(乱流域)

・亜音速域のサイジング

・超音速域のサイジング

参照:IEC(International Electrotechnical Commission) 60534-2-1

バルブ特性

代表的なバルブタイプに応じた特性線図。

オリフィスの圧損計算

非圧縮性流体のサイジング式

$$q= C \cdot A \sqrt{\frac{2 \Delta p}{\rho}}$$ $$C= \frac{C_d}{\sqrt{1-\beta^4}} $$

Reynolds数>2e5の領域では、C≒0.61

※Reynolds数は、入口径、入口物性より求める

安全弁の圧損計算

亜音速域の基本式

API RP 520参照

$$F=\frac{44.75159 \cdot Position \cdot Area \cdot K_d \cdot F_2 \cdot \sqrt{DP \cdot R \cdot MW}}{MW}$$

ただし、液の場合\(F_2\)は1、ガスの場合は以下式より求めます

$$F_2=\sqrt{\left(\frac{C_p/C_v}{C_p/C_v-1}\right) \cdot \left(\frac{P_x}{P_i}\right)^{2C_v/C_p} \cdot \left(\frac{1-\left(\frac{P_x}{P_i}\right)^{\frac{C_p/C_v-1}{C_p/C_v}}}{1-\frac{P_x}{P_i}}\right)}$$

F:モル流量、Kd:Discharge Coefficient、Cp/Cv:比熱比、R:モル密度、Pi:入口圧力、Px:出口圧力

超音速域の基本式

$$F=\frac{0.02119 \cdot Position \cdot C \cdot Area \cdot K_d \cdot \sqrt{DP \cdot UGC \cdot R \cdot MW}}{MW}$$

ただし、Cは以下式より求めます

$$C=520 \cdot \sqrt{C_p/C_v \cdot \left(\frac{2}{C_p/C_v+1}\right) \cdot \left(\frac{C_p/C_v+1}{C_p/C_v-1}\right)}$$

安全弁のサイジング

ガス、液サービスのサイジング


必要面積の算出

ASME Codeで規定されている必要面積は以下の式より求めます。

$$ A=\frac{277.8W}{K_d K_b K_c K_v G}$$

各パラメーターは以下の通り。

Kd : discharge coefficient/p>

ガスサービスの場合は、0.975。

液サービスの場合は、subcooled:0.85, 飽和液:0.65を使用。

Kb : backpressure correction factor ()

KbはAPI記載のチャートより求めます。

Kc : combination crrection factor

PRVの上流にrapture diskが有る場合:0.9、無い場合:1.0。

Kv : viscosity correction factor

API記載のチャートから求めるか、以下の式により算出。

$$K_v=\left(0.9935+\frac{2.878}{Re^{0.5}}+\frac{342.75}{Re^{1.5}}\right)^{-1.9}$$

2相流


HEM(Homogeneous Equilibrium Method)

HEMの前提①:流体を仮想的な単相流体として取り扱うことができる。

HEMの前提②:Nozzleが4in以上あれば、熱的、機械的に平衡と考えることができる。

上記の前提がそぐわない場合、HNE(Homogeneous Nonequilibrium Method)で検討します。

・OMEGA method(cf.Leung Omega Method)の計算の流れ(API 520 Annex C参照)

①OMEGAの算出

$$ \omega=9\left(\frac{v_9}{v_0} -1 \right)$$

ここで、\(v_9\)はPRV 入口圧力の90%におけるspecific volume [ft3/lb] 。

また、\(v_0\)はPRV 入口圧力におけるspecific volume [ft3/lb] 。

②臨界条件の判定

上式で求めた\(\omega\)を使って、臨界圧力比\(\eta_c\)を求めます。

$$ \eta_{c}^2 + \left(\omega^2-2\omega\right)\left(1-\eta_c\right)^2+2\omega^2ln\eta_c+2\omega^2\left(1-\eta_c\right)=0 $$ $$ P_c = \eta_{c} P_{0} $$ $$ P_c \geq P_a \Rightarrow critical flow $$ $$ P_c < P_a \Rightarrow subcritical flow $$

③質量流束G[kg/s-m2]の算出

Critical Flowの場合

$$ G = \eta_c \sqrt{\frac{P_0}{v_0\omega}}$$

SubCritical Flowの場合

$$ G = \frac{\sqrt{-2\left(\omega ln \eta_a+(\omega-1)(1-\eta_a)\right)}}{\omega (1/\eta_a-1)+1}\sqrt{\frac{P_a}{v_a}}$$

ここで、\(eta_a\)はback pressure ratio[-] 。

$$ \eta_a=\frac{P_a}{P_0}$$
HDI(Homogeneous Direct Integration)

安全弁入口から出口まで圧力で直接積分していきます。

$$ G^2=(\rho_1^2) \cdot (-2\int_{P_0}^{P_1}\frac{dP}{\rho})$$

・積分部分を数値計算で離散的に積分していく場合

$$ \int_{P_0}^{P_1}\frac{dP}{\rho} \approx \sum 2\frac{P_i-P_{i-1}}{\rho_{i}+\rho_{i-1}}$$

配管の圧損計算

亜音速域の基本式

層流域の摩擦係数式(Hagen-Poiseuille lawからの導出)

$$f=\frac{64}{Re}$$

平滑管の摩擦係数式

適用範囲:Re>4000

$$\frac{1}{\sqrt{f}}=2log_{10}(Re\sqrt{f})-0.8$$

遷移域の摩擦係数式(Colebrook-White式)

$$\frac{1}{\sqrt{f}}=1.14-2log_{10}(\epsilon/D+9.35/Re\sqrt{f})$$

乱流域の摩擦係数式(Colebrook式)

Colebrook式適用範囲:Re>4000

$$\frac{1}{ \sqrt {f} }=-2log\left(\frac{\epsilon/D}{3.7}+\frac{2.51}{Re\sqrt{f}}\right)$$

Colebrook式は収束計算が必要となりますので、計算例として、ExcelのGoalSeakを用いる方法、VBAを用いた計算方法をご紹介します。

なお、本サイトでは以下Darcy式を圧損計算の基本式として考えています。

$$DP=f \cdot \frac{L}{D} \frac{v^2}{2g}$$

機械系ではFanningの式を用いる場合が多いかもしれませんが、その場合は、前述の摩擦係数は1/4になります。

まず計算しやすいように、以下のような前提、数値で圧損計算を行います。

流量W:100[t/hr]、配管径D:1.0[m]、配管長L:100[m]、配管粗さε:0.05[mm]

密度ρ:1[kg/m3]、粘度μ:1[cP]、フィッティング無し。

流速を計算します。

$$V=\frac{W}{\rho} \cdot \frac{1}{D^2/4 \cdot \pi}=\frac{100000[kg/hr]}{1[kg/m3]} \cdot \frac{1}{1^2/4 \cdot \pi [m^2]}=127324[m/hr]$$

続いてReynolds数を求めます。

$$Re=\frac{V \cdot \rho \cdot D}{\mu}=\frac{127324/3600[m/s] \cdot 1[kg/m3] \cdot 1[m]}{1/1000[Pa \cdot s]}=35368[-]$$

摩擦係数fは、Colebrook式で求める場合、上述のように収束計算が必要となります。

まずはGoalSeekで求める方法です。以下のように、適当なfを初期値として与え、Colebrook式の左辺、右辺の計算を行い、左辺と右辺の差が0となるようにGoalSeekで設定します。

GoaslSeek

GoalSeekだとエラー値の評価が甘いですが、Solverを使えば、Error値を小さくすることができます。

VBAで求める場合は、例えば逐次代入法を用いて解くことができます。(Do Loopで、右辺の値の2乗の逆数を更新していき、前回値との差分がtolerance以下になるまで計算を繰り返します。)

以下のように関数として登録しておき、ExcelからVBA関数を呼び出すことが可能です。

fcalVBA

また、MillerのFlow Measurement(Flow Measurement Engineering Handbook 3rd (third) Edition by Miller Richard published by McGraw-Hill Professional (1996))ではNewton法による解法が紹介されています。

具体的には、以下のように陽関数の形に変換し、微分値を求め、初期値f0から計算を行います。

$$F= - f^{-1/2} -2log\left(\frac{\epsilon/D}{3.7}+\frac{2.51}{Re\sqrt{f}}\right)$$ $$F'= \frac{1}{2} f^{-3/2} -\frac{4.034}{(\epsilon/D)Ref^{3/2}+9.287f}$$ $$f0= 0.25 \left[log \left( \frac{0.2703 \epsilon}{D} + \frac{5.74}{Re^{0.9}} \right) \right]^{-2} $$

諸条件から圧損を求める簡易ツールをこちらに作成しました。簡易版ですが、確認用にご使用ください。

DPTool

超音速域の基本式

代表的なFitting類のK値(Resistance Coefficients)

エルボーやティー等、フィッテイングの圧損を考える場合、前述のDP計算において\(K=f \cdot L/D\)に以下の補正をかけます。

$$K=f \cdot \frac{L}{D}$$

各フィッティングのK値については、以下代表的なものを記載いたします※Crane Engineering Handbook "K" Factor Tableより一部抜粋

Fitting K valueK
Gate valves (Wedge Disc)8f
Swing Check valves 100f
Stop Check valves 400f
Glove valves 340f
Ball valves 3f
Butterfly valves (2 to 8")45f
Butterfly valves (10 to 14")35f
Butterfly valves (16 to 24")25f
Standard elbows (90°)30f
Standard elbows (45°)16f
Standard tees (Flow through branch) 60f
Standard tees (Flow through run tees) 20f
Pipe Engrance(Inward Projecting) 0.78
Pipe Exit 1.0

ザックリ相当長ベースで入力したい場合は、以下式にて相当長\(L_{eq}\)を入力します。

$$K=f \cdot \frac{L+L_eq}{D}$$

急拡大、急縮小による損失

急拡大、急縮小等、K値の補正が必要な場合は、以下式を使用します。(徐々に管径が変わる場合は、こちらをご参照ください)

$$K=f \cdot \frac{L+L_eq}{D} + K_{additional}$$

急拡大

$$K_{additional}=\left(1-\frac{d_{1}^2}{d_{2}^2}\right)^2$$

急縮小

$$K_{additional}=0.5\left(1-\frac{d_{1}^2}{d_{2}^2}\right)$$

※1,2の添え字はそれぞれ小口径、大口径配管を示しています。

2相流配管の圧損(Beggs-Brill式)

2相流配管の圧損計算式としてBeggs-Brill式を紹介します。

出典

Beggs, H.D. and Brill, J.P.; "A study of two-phase flow in inclined pipe;" Journal of Petroleum Technology, 1973, p.607–617

Beggs-Brill式の精度について

Baker, A., Nielsen, K., and Gabb, A; "Pressure loss, liquid-holdup calculations developed" Oil and Gas Journal, 1988, p.55–59

2相流配管の圧損(Lockhart-Martineli式)

考え方はBBMHよりも単純で気相、液相が単相で流れた場合の圧損比をLockhart-MartineliのパラメーターXの相関式で表現し、以下のように計算します。

1.液側、ガス側のReynolds数を求める。⇒各相がTurbulentかLaminarか判断する。⇒使用する係数式が決まる

2.液側の圧力損失を計算する。(液側の重量流量、密度からColebrook式等で通常通り計算する)

3.1で求めた係数と2で求めた液側の圧力損失を掛け合わせて、2相流の圧力損失を求める。

Lockhart-Martineli式の精度について

Perryのhandbookに以下の記載があります。

出典

偏流解析

簡単のため、まずは水のみ、密度一定とした場合の偏流計算を考えてみます。

下図のように、圧力P1のタンクから、各ユーザー(User1、User2)に、流れる流量を考えます。

FlowSimple1

上図で、送出元と各ユーザーの到達圧、密度が分かっている場合、配管径が決まれば、ヘッダー圧、各ユーザーへの流量は簡単に求めることができます。

例えば、以下のような問題設定で考えてみます。

FlowSimple2

まず未知数を見てみます。赤字で示していますが、各配管での流量、圧損とヘッダー圧の計7変数が未知数の状態ですので、7つの式が必要となります。

今回は密度一定、つまり非圧縮性の流体を考えていますので、物質収支、圧力バランスで蓄積項を考慮する必要は無く、以下のように考えられます。今回はNewton法で求めるため、予めf(x)=0の形に変形しておきます。

FlowSimpleEq1

これで7変数に対し、7つの式を立てることができ、自由度は満足していますので、後は解くだけです。

とはいえ、流量-圧損の関係式が非線形となっており、式を見てパッと解ける連立方程式にはなっていないため、今回は①Pythonを使ってNewton-Raphson法で解く方法と、②Excelのsolverを使う方法をご紹介します。

Pythonで求める場合

まずは境界条件と適当な初期値を決めます。初期値は適当ですが、一応マスバランスは満たした形で決めておきます。

FlowSimplePy1

f(x)=の誤差がtolerance内に収まるまで繰り返し計算を行います。ここで、逆流が発生した場合に、ルート内が負とならないように場合分けしておきます。

FlowSimplePy2

あとはJacobianを手計算で求め、その逆行列と関数をかけ合わせれば良いのですが、ここでnumpyでは非常に便利な関数があります。

linalg.solve(Jacobian、関数)で、Jacobianと関数の積を一気に求めることができます。

FlowSimplePy3

刻み0.5、tolerance:e-5で、収束回数:14回。最終的な解は以下のように求まります。

FlowSimplePy4

なお、前述のように、逆流が発生する場合の式処理も行っておけば、以下のようにユーザー側の圧力が上昇した場合の逆流も計算することが可能です。

FlowSimplePy5

なお、以下のように収束計算の初めに摩擦係数を計算させてしまえば(収束計算になります)、配管径、配管長、フィッティング情報からconductance値を求めることも可能です。

FlowSimplePy6

Cf)Csharpで求める場合

上述のPythonを使ったケース同様、逆行列計算さえできれば簡単に解くことができます。

C++では、数値計算ライブラリとしてEigenが有名ですが、Csharpでも同様のものを探していたところ、Math.Net.Numericsというものがありました。

詳細は、Csharpのページでゆくゆくまとめていきたいと思いますが、ひとまず暫定版として、上述のPythonと同様の簡易サンプルを以下に作成しました。

Newton法による簡易偏流計算

簡単な使い方を以下に記載します。

DriftCalcsample

Excelで求める場合

Excelを使う場合、いくつか方法が考えられます。①VBAを使って上述のようにNewton-Raphson法を用いる方法、②Excelのシート上でNewton-Raphson法を用いる方法(Excelであれば逆行列計算、行列積も関数で行うことができます。)、③①と②の複合(②だけで行うと収束計算分セルを拡張する必要があるので、VBAと複合して使った方が見やすい結果になります)、④Solverを使う方法等です。

①~③はPythonを使う方法と考え方は同様なので省略し、④の方法を紹介します。

問題設定はPythonの場合と同様とします。

拡大管/縮小管の圧損計算

Craneの計算式に基づく計算(AppendixA A-26参照)

$$DP= K \cdot \rho \cdot \frac{v^2}{2} $$

※DP[Pa]

拡大管の場合

$$ \cdot \theta <=45 $$ $$ \ K = 2.6 * sin ( \frac{\theta}{2} ) * ( 1- \beta^2 )^2$$ $$ \cdot \theta >45 $$ $$ \ K = ( 1- \beta^2)^2$$

縮小管の場合

$$ \cdot \theta <=45 $$ $$ \ K = 0.8* sin ( \frac{\theta}{2} ) * ( 1- \beta^2) $$ $$ \cdot \theta >45 $$ $$ \ K = 0.5 * \sqrt{sin ( \frac{\theta}{2} ) } * ( 1- \beta^2) $$

拡大管/縮小管の圧損計算式の出典

・Crane Valve ; "Flow of Fluids Through Valves, Fittings & Pipe"

こちらのリンクから計算可能です

蒸留塔の圧損計算

トレイ塔の渇き圧損

  Fair式

充填塔の渇き圧損

充填塔の渇き圧損計算式としてSticlmair式をご紹介します。

粒子径の相当長[m]:\(Dia_P\)は、ベッドのVoid fraction、充填物の比表面積[1/m]:PackAreaを用いて下式となります。

$$Dia_P=\frac{6 \cdot \left(1-Void\right)}{PackArea}$$

充填塔の空塔速度[m/sec]:\(U_v\)は、ガスの通過流量[kg-mol/sec]:\left(F_v\rght)と、ガス密度[kg-mol/m3]:\(R_v\)、カラムの断面積[m2]:Aを用いて下式となります。

$$U_v=\frac{F_v}{R_v A}$$

抵抗係数[-]:\(f_0\)は、Stichlmair係数[-]:C1~C3とガスのレイノルズ数[-]:\(Re_{vap})を用いて下式となります。

$$f_0=\frac{C_1}{Re_{vap}}+\frac{C_2}{2 \cdot \sqrt{Re_{vap}}}+C_3$$ $$c=\frac{\left(- \frac{C_1}{Re_{vap}} - \frac{C_2}{2 \cdot \sqrt{Re_{vap}}}\right)}{f_0}$$ $$\Delta P_d = \frac{1}{1000 \cdot KJ^2}\frac{3}{4}f_0 \frac{1-Void}{Void^{4.65}}R_{gas} \frac{PackHeight}{Dia_P}U_v^2$$

最終的な渇き圧損[kPa]:\( \Delta P\)は、下式となります。

$$\Delta P = \Delta P_d \left(\frac{1-Void \cdot \left(1-\frac{Beta}{Void}\right)}{1-Void}^{\frac{2+c}{3}}\right) \cdot \left(1-\frac{Beta}{Void}\right)^{-4.65}$$

Stichlmair係数C1~C3、比表面積、Void Fractionは各充填物タイプに応じて求められています。

以下に一部抜粋して示します。詳細は出展をご参照ください。

規則充填物

PackingType/sizearea[1/m]void[-]C1[-]C2[-]C3[-]
MontZB1/3003000.97230.9
B1/2002000.98241.0
B1/1001000.99271.0
SulzerMellapak(plastic)/250Y2500.85110.32

出典

Stichlmair, J.J., Bravo, J.L., Fair, J.R.,; “General model for prediction of pressure drop and capacity of counter current gas/liquid Packed Columns” Gas Separation & Purification,. Vol. 3, March 1989, p. 19-28

充填塔のフラッディング計算

Stichlmairの式   

充填塔の渇き圧損計算式としてSticlmair式をご紹介します。

  

以下の流れで収束計算を行います。

  

a.ガス流速を仮定します。

  

b.乾き圧損を計算します。

  

c.濡れ圧損を計算します。

  

d.下式の残差を求め、収束誤差を超えている場合は再度ガス流速をアップデートし、aに戻ります。

こちらのリンクから計算可能です

出典

Stichlmair, J.J., Bravo, J.L., Fair, J.R.,; “General model for prediction of pressure drop and capacity of counter current gas/liquid Packed Columns” Gas Separation & Purification,. Vol. 3, March 1989, p. 19-28

トレイ塔のフラッディング計算

Souders-Brownの係数 $$\sqrt{\frac{\rho_{vap}}{\rho_{liq}-\rho_{vap}}}$$

トレイ塔のウィーピング計算

当サイトに不具合、ご意見等ございましたらCEsolutionにお知らせください。

Page Top