Rでコールブルック・ホワイトの式

使用例

円管(\(D_H = 0.1\) m)、商業用鋼管(\(\varepsilon = 0.045\) mm)、\(Re = 10^5\) のときの管摩擦係数を求めます。

# 管摩擦係数 lambda を求める
lambda_value <- colebrook_white(epsilon = 0.045e-3, D_H = 0.1, Re = 1e5)
lambda_value

[1] 0.02011523

求めた \(\lambda\) を使って、残りの3変数を逆算し、元の値に戻ることを確認します。

# 逆算して元の値に戻るかを確認する
epsilon <- colebrook_white(lambda = lambda_value, D_H = 0.1, Re = 1e5) # epsilon
D_H <- colebrook_white(lambda = lambda_value, epsilon = 0.045e-3, Re = 1e5) # D_H
Re <- colebrook_white(lambda = lambda_value, epsilon = 0.045e-3, D_H = 0.1) # Re
list(epsilon = epsilon, D_H = D_H, Re = Re)

$epsilon
[1] 4.5e-05

$D_H
[1] 0.1

$Re
[1] 1e+05

複数の値を同時に与えることもできます。レイノルズ数を変えたときの \(\lambda\) の変化を見てみます。

# Re を変化させて lambda を求める
re_seq <- 10^(4:7)
data.frame(
  Re = re_seq,
  lambda = colebrook_white(epsilon = 0.045e-3, D_H = 0.1, Re = re_seq)
)

     Re     lambda
1 1e+04 0.03156739
2 1e+05 0.02011523
3 1e+06 0.01684944
4 1e+07 0.01635936

方程式の形

コールブルック・ホワイトの式に \(x = 1/\sqrt{\lambda}\) を代入し、全てを左辺に移すと次のようになります。

\[
\begin{eqnarray}
g(x) &=& x + 2\log_{10}\left(a + b\,x\right) \;=\; 0
\end{eqnarray}
\tag{5}\]

\[
\begin{eqnarray}
a &=& \frac{\varepsilon}{3.71\,D_H} \;\ge\; 0, \qquad b \;=\; \frac{2.51}{Re} \;>\; 0
\end{eqnarray}
\tag{6}\]

\(a\) は相対粗度に由来する定数、\(b\) は粘性に由来する定数で、どちらも \(x\) には依存しません。

求める \(\lambda\) は、\(x>0\) の範囲にある \(g(x)=0\) の解から \(\lambda = 1/x^2\) で得られます。

単調増加であること

\(g(x)\) を \(x\) で微分します。

\(\log_{10} u = \ln u / \ln 10\) を使うと、次のようになります。

\[
\begin{eqnarray}
g'(x) &=& 1 + \frac{2}{\ln 10}\cdot\frac{b}{a + b\,x}
\end{eqnarray}
\tag{7}\]

\(x>0\) では \(a + b\,x > 0\) であり、\(b>0\) です。したがって第2項は正になり、

\[
\begin{eqnarray}
g'(x) &>& 1 \;>\; 0
\end{eqnarray}
\tag{8}\]

となります。つまり \(g(x)\) は \(x>0\) の全域で厳密に単調増加です。

解が少なくとも1つ存在すること

\(g(x)\) は \(x>0\) で連続です。区間の両端での符号を調べます。

• \(x \to 0^+\) のとき、\(g(x) \to 2\log_{10} a\) です。\(a<1\)(\(\varepsilon/D_H < 3.71\))なら負になります。\(a=0\) の滑面管では \(g(x) \to -\infty\) です。
• \(x \to \infty\) のとき、\(x\) の項が支配的になり、\(g(x) \to +\infty\) です。

両端で符号が異なるため、中間値の定理により \(g(x)=0\) となる \(x\) が少なくとも1つ存在します。

解がただ1つであること

厳密に単調増加な関数は、同じ値を2回とりません。

仮に \(x_1 < x_2\) がともに解なら \(g(x_1) < g(x_2)\) となり、\(g(x_1) = g(x_2) = 0\) と矛盾します。

よって解はちょうど1つです。

さらに、\(\lambda = 1/x^2\) は \(x>0\) で一対一の対応ですので、\(\lambda\) の解もただ1つに決まります。

コードとの対応

コードでは、探索区間を interval = c(1e-6, 1e3) としています。

この区間で符号が変わること(\(g(10^{-6})<0\) かつ \(g(10^3)>0\))が、uniroot() を使うための条件です。

• \(g(10^3)>0\) は、\(Re\) が極端に大きくない限り、常に成り立ちます。
• \(g(10^{-6})<0\) は、ほぼ \(a<1\) と同じ条件です。コードの if (a >= 1) は、この条件を確認するためのものです。

実用上は、相対粗度 \(\varepsilon/D_H\) が 0.05 程度以下であれば、この条件には余裕があります。

項ごとの微分

微分の線形性より、次のようになります。

\[
\begin{eqnarray}
g'(x) &=& \frac{d}{dx}\,x \;+\; 2\,\frac{d}{dx}\log_{10}\left(a + b\,x\right)
\end{eqnarray}
\tag{10}\]

第1項は \[\dfrac{d}{dx}\,x = 1 \tag{11}\] です。

常用対数を自然対数に直す

底の変換公式 \(\log_{10} u = \dfrac{\ln u}{\ln 10}\) を使います。\(\ln 10\) は定数です。

\[
\begin{eqnarray}
\frac{d}{dx}\log_{10}\left(a + b\,x\right) &=& \frac{1}{\ln 10}\cdot\frac{d}{dx}\ln\left(a + b\,x\right)
\end{eqnarray}
\tag{12}\]

合成関数の微分

\(u = a + b\,x\) とおくと、\(\ln u\) の微分は合成関数の微分(連鎖律)から、

\[
\begin{eqnarray}
\frac{d}{dx}\ln u &=& \frac{1}{u}\cdot\frac{du}{dx}
\end{eqnarray}
\tag{13}\]

となります。\(a\) と \(b\) は定数ですので、

\[
\begin{eqnarray}
\frac{du}{dx} &=& \frac{d}{dx}\left(a + b\,x\right) \;=\; b
\end{eqnarray}
\tag{14}\]

です。したがって、

\[
\begin{eqnarray}
\frac{d}{dx}\ln\left(a + b\,x\right) &=& \frac{b}{a + b\,x}
\end{eqnarray}
\tag{15}\]

となります。

最終導出

式 10 に 式 11 、式 12、式 15 を代入します。

\[
\begin{eqnarray}
g'(x) &=& 1 + 2\cdot\frac{1}{\ln 10}\cdot\frac{b}{a + b\,x} \\
&=& 1 + \frac{2}{\ln 10}\cdot\frac{b}{a + b\,x}
\end{eqnarray}
\tag{16}\]

式 7 が導出されます。

補足(符号の確認)

\(x>0\) では \(a + b\,x > 0\) です(\(a \ge 0\)、\(b>0\) のため)。また \(\ln 10 \approx 2.303 > 0\) です。そのため第2項は常に正になり、

\[
\begin{eqnarray}
g'(x) &>& 1
\end{eqnarray}
\tag{17}\]

が成り立ちます。これが「\(g(x)\) は単調増加」の根拠です。

コードでの扱い

solve_lambda() には if (a >= 1) の確認を入れてあり、その場合は警告を出して NA を返します。

以上です。