IAPWS95/thermal_conductivity.h マニュアル

(The documentation of IAPWS95/thermal_conductivity.h)
Last Update: 2023/9/25


IAPWS95/thermal_conductivity.hでは 熱伝導率を計算する関数が定義されている。
Functions to compute the thermal conductivity are defined in IAPWS95/thermal_conductivity.h.
熱伝導率の計算式はWagner and Pruss (2002)には含まれておらず、 Huber et al. (2012)において別途定義されている。 2023年9月21日現在、 NIST Chemistry WebBook において同論文の式が用いられている。
The equation of the thermal conductivity is not included in Wagner and Pruss (2002) but is defined separately in Huber et al. (2012). The equation in that paper is used in NIST Chemistry WebBook as of 21 September 2023.
Huber et al. (2012)において熱伝導率\(\lambda\)は \[\begin{equation} \frac{\lambda}{\lambda_{ref}} =\bar{\lambda}_0(\tau) \bar{\lambda}_1(\delta,\tau) +\bar{\lambda}_2(\delta,\tau) \label{eq.lambda} \end{equation}\] の形で与えられる(同論文(2)式)。 ここで\(\lambda_{ref}\equiv 0.001\) [J s\(^{-1}\) m\(^{-1}\) K\(^{-1}\)]は 基準となる熱伝導率を表す。 また\(\delta\equiv\rho/\rho_c\)は無次元化した密度、 \(\tau\equiv T_c/T\)は無次元化した温度の逆数であり、 これらはIAPWS-95全体の表記に合わせてある (そのためHuber et al. (2012)の表記とは異なる)。
Huber et al. (2012) gives the thermal conductivity \(\lambda\) in the form of Eq. (\ref{eq.lambda}), which is from Eq. (2) of the paper, where \(\lambda_{ref}\equiv 0.001\) [J s\(^{-1}\) m\(^{-1}\) K\(^{-1}\)] is a reference thermal conductivity; \(\delta\equiv\rho/\rho_c\) is a nondimensional density and \(\tau\equiv T_c/T\) is the inverse of a nondimensional temperature, defined consistently with the notation of IAPWS-95. Note that the notation is different from that of Huber et al. (2012).
\(\bar{\lambda}_0\)は低密度極限での熱伝導率を表し、 \[\begin{equation} \bar{\lambda}_0(\tau)= \frac{1}{\sqrt{\tau}\sum_{k=0}^4 L_k \tau^k} \label{eq.lambda0} \end{equation}\] で与えられる(論文(4)式)。ここで\(L_k\)は表1のように与えられる。
The factor \(\bar{\lambda}_0\) represents the thermal conductivity in the limit of zero density and is given as Eq. (\ref{eq.lambda0}), which is from Eq. (4) of the paper. The values of \(L_k\) are given in Table 1.
表1. (\ref{eq.lambda0})式の係数\(L_k\)。 Huber et al. (2012)の表2に基づく。
Table 1. The coefficients \(L_k\) of Eq. (\ref{eq.lambda0}) from Table 2 of Huber et al. (2012).
\(k\) \(L_k\)
0 2.443221e-03
1 1.323095e-02
2 6.770357e-03
3 -3.454586e-03
4 4.096266e-04

\(\bar{\lambda}_1\)は密度増加の熱伝導率への寄与を表し、 \[\begin{equation} \bar{\lambda}_1(\delta,\tau)= \exp\left[ \delta \sum_{i=0}^4 (\tau-1)^i \sum_{j=0}^5 L_{ij} (\delta-1)^j \right] \label{eq.lambda1} \end{equation}\] で与えられる(同論文(5)式)。 ここで\(L_{ij}\)は表2のように与えられる。
The factor \(\bar{\lambda}_1\) represents the contribution to the thermal conductivity due to increasing density and is given as Eq. (\ref{eq.lambda1}) which is from Eq. (5) of the paper. The values of \(L_{ij}\) are given in Table 2.
表2. (\ref{eq.lambda1})式の係数\(L_{ij}\)。 Huber et al. (2012)の表4に基づく。
Table 2. The coefficients \(L_{ij}\) of Eq. (\ref{eq.lambda1}) from Table 4 of Huber et al. (2012).
\(j=0\) \(j=1\) \(j=2\) \(j=3\) \(j=4\) \(j=5\)
\(i=0\) 1.60397357 -0.646013523 0.111443906 0.102997357 -0.0504123634 0.00609859258
\(i=1\) 2.33771842 -2.78843778 1.53616167 -0.463045512 0.0832827019 -0.00719201245
\(i=2\) 2.19650529 -4.54580785 3.55777244 -1.40944978 0.275418278 -0.0205938816
\(i=3\) -1.21051378 1.60812989 -0.621178141 0.0716373224 0 0
\(i=4\) -2.7203370 4.57586331 -3.18369245 1.1168348 -0.19268305 0.012913842

\(\bar{\lambda}_2\)は臨界点近傍での振る舞いを表し、 論文(13)式, (10)式, (9)式, (16)式, (17)式により \[\begin{equation} \bar{\lambda}_2(\delta,\tau)= \Lambda\frac{\delta}{\tau} \frac{c_p(\delta,\tau)}{R} \frac{\mu_{ref}}{\mu(\delta,\tau)} Z(\delta,\tau) \label{eq.lambda2} \end{equation}\] \[\begin{eqnarray} Z(\delta,\tau) &=& \frac{2}{\pi y(\delta,\tau)}\left\{ \left[ \frac{c_p(\delta,\tau)-c_v(\delta,\tau)}{c_p(\delta,\tau)} \tan^{-1}y(\delta,\tau) +\frac{c_v(\delta,\tau)}{c_p(\delta,\tau)}y(\delta,\tau) \right] \right. \nonumber \\ & & \left. -\left[ 1-\exp\frac{-3\delta^2 y(\delta,\tau)}{y(\delta,\tau)^3+3\delta^2} \right] \right\} \label{eq.Z} \end{eqnarray}\] \[\begin{equation} y(\delta,\tau)= \bar{q}_D \xi(\delta,\tau) \label{eq.y} \end{equation}\] \[\begin{equation} \xi(\delta,\tau)= \xi_0\left[\frac{\Delta\bar{\chi}(\delta,\tau)}{\Gamma_0}\right]^{\nu/\gamma} \label{eq.xi} \end{equation}\] \[\begin{equation} \Delta\bar{\chi}(\delta,\tau)= \max\left\{0, \delta\frac{P_c}{\rho_c}\left[ \frac{1}{\left(\PartialDiff{P}{\rho}\right)_T} -\frac{1}{\left.\left(\PartialDiff{P}{\rho}\right)_T\right|_{\tau=\tau_R}} \frac{\tau}{\tau_R} \right] \right\} \label{eq.Delta_chi} \end{equation}\] と書ける。 ここで\(c_p\)は定圧比熱、\(c_v\)は定積比熱、\(\mu\)は粘性係数、\(P\)は圧力、 \(P_c\)は臨界点での圧力、\(R\)は気体定数、 \(\mu_{ref}\equiv 1\times 10^{-6}\) [Pa s]は基準となる粘性係数であり、 他の定数は表3のように与えられる。
The factor \(\bar{\lambda}_2\) represents the behaviour of the thermal conductivity near the critical point, given as Eqs. (\ref{eq.lambda2})-(\ref{eq.Delta_chi}) according to Eqs. (13), (10), (9), (16), and (17) of the paper. In these equations, \(c_p\) is the isobaric heat capacity, \(c_v\) is the isochoric heat capacity, \(\mu\) is the viscosity, \(P\) is a pressure, \(P_c\) is the pressure at the critical point, \(R\) is the gas constant, and \(\mu_{ref}\equiv 1\times 10^{-6}\) [Pa s] is a reference viscosity. The other constants in these equations are given in Table 3.
(\ref{eq.Z})式の分母に\(y(\delta,\tau)\)が存在するため、 このままでは計算できないケースが存在する。 この問題について検討する。 以下では簡単のため、引数の\(\delta\), \(\tau\)は省略する。 まず(\ref{eq.Delta_chi})式より\(\Delta\bar{\chi}\geq 0\)であり、 この不等式を(\ref{eq.xi})(\ref{eq.y})に代入すると \(\xi\geq 0\), \(y\geq 0\)であることが分かる。 特に\(\Delta\bar{\chi}=0\)のとき、\(\xi=0\), \(y=0\)となる。 \(y=0\)での\(Z\)の振る舞いを調べるため、(\ref{eq.Z})式を \[\begin{equation} Z=\frac{2}{\pi y}f(y), \hspace{1em} f(y)=\left[\frac{c_p-c_v}{c_p}\tan^{-1}y+\frac{c_v}{c_p}y\right] -\left[1-\exp\frac{-3\delta^2 y}{y^3+3\delta^2}\right] \label{eq.Z.divide} \end{equation}\] と分割する。 \[\begin{eqnarray} f’(y) &=& \left[\frac{c_p-c_v}{c_p}\frac{d}{dy}\left(\tan^{-1}y\right) +\frac{c_v}{c_p}\right] -\left[\exp\left(\frac{-3\delta^2 y}{y^3+3\delta^2}\right) \frac{d}{dy}\left(\frac{3\delta^2 y}{y^3+3\delta^2}\right)\right] \nonumber \\ &=& \left[\frac{c_p-c_v}{c_p}\frac{1}{y^2+1}+\frac{c_v}{c_p}\right] -\exp\left(\frac{-3\delta^2 y}{y^3+3\delta^2}\right) \frac{3\delta^2(y^3+3\delta^2)-3\delta^2 y\cdot 3y^2}{(y^3+3\delta^2)^2} \nonumber \\ &=& \frac{c_p-c_v}{c_p}\frac{1}{y^2+1}+\frac{c_v}{c_p} -\exp\left(\frac{-3\delta^2 y}{y^3+3\delta^2}\right) \frac{3\delta^2y^3+9\delta^4-9\delta^2 y^3}{(y^3+3\delta^2)^2} \nonumber \\ &=& \frac{c_p-c_v}{c_p}\frac{1}{y^2+1}+\frac{c_v}{c_p} -\exp\left(\frac{-3\delta^2 y}{y^3+3\delta^2}\right) \frac{-6\delta^2y^3+9\delta^4}{(y^3+3\delta^2)^2} \label{eq.fdash} \end{eqnarray}\] であるので \[\begin{eqnarray} f(0) &=& \left[\frac{c_p-c_v}{c_p}0+\frac{c_v}{c_p}0\right] -\left[1-\exp\frac{0}{0+3\delta^2}\right] \nonumber \\ &=& [0+0]-[1-\exp(0)] \nonumber \\ &=& 0 \label{eq.f0} \end{eqnarray}\] \[\begin{eqnarray} f’(0) &=& \frac{c_p-c_v}{c_p}\frac{1}{0+1}+\frac{c_v}{c_p} -\exp\left(\frac{0}{0+3\delta^2}\right) \frac{-6\delta^2 0+9\delta^4}{(0+3\delta^2)^2} \nonumber \\ &=& \frac{c_p-c_v}{c_p}+\frac{c_v}{c_p} -\exp(0)\frac{9\delta^4}{9\delta^4} \nonumber \\ &=& 0 \label{eq.fdash0} \end{eqnarray}\] となり、\(f(y)\)の\(y=0\)のまわりでのテイラー展開は \(y^2\)以降の項しか残らないことが分かる。 したがって(\ref{eq.Z.divide})式より \(y=0\)の極限で\(Z=0\)であることが分かる。 コンピュータプログラム内ではこのことを用いて \(y=0\)の場合と\(y>0\)の場合で場合分けして\(Z\)を計算する。 \(y=0\)の場合には\(Z=0\)になるという上記の結論を直接用いる。
Because Eq. (\ref{eq.Z}) has \(y(\delta,\tau)\) in the denominator, the computation fails in some cases. Let us consider this issue. The arguments \(\delta\) and \(\tau\) are omitted for simplicity. Eq. (\ref{eq.Delta_chi}) shows that \(\Delta\bar{\chi}\geq 0\), and inserting this relation into Eqs. (\ref{eq.xi}) and (\ref{eq.y}) results in \(\xi\geq 0\) and \(y\geq 0\). Especially, \(\xi=0\) and \(y=0\) occurs in case of \(\Delta\bar{\chi}=0\). To survey the behaviour of \(Z\) around \(y=0\), divide Eq. (\ref{eq.Z}) into two equations as (\ref{eq.Z.divide}). The derivative \(f’(y)\) is calculated as (\ref{eq.fdash}), and thus \(f(0)=0\) and \(f’(0)=0\) hold (Eqs. \ref{eq.f0} and \ref{eq.fdash0}). This means that the Taylor expansion of \(f(y)\) around \(y=0\) has only the \(y^2\) and higher-order terms. Therefore, from Eq. (\ref{eq.Z.divide}), we obtain \(Z=0\) at the limit of \(y=0\). The computer program computes \(Z\) using branches for \(y=0\) and \(y>0\); for \(y=0\), the conclusion \(Z=0\) obtained above is used directly.
表3. (\ref{eq.lambda2})-(\ref{eq.Delta_chi})式の定数。 Huber et al. (2012)の表5に基づく。
Table 5. The constants in Eqs. (\ref{eq.lambda2})-(\ref{eq.Delta_chi}) from Table 5 of Huber et al. (2012).
定数
Constant 		値
Value
\(\Lambda\) 177.8514
\(\bar{q}_D\) \(\frac{1}{0.4\times 10^{-9}}\) [m\(^{-1}\)]
\(\nu\) 0.630
\(\gamma\) 1.239
\(\xi_0\) \(0.13\times 10^{-9}\) [m]
\(\Gamma_0\) 0.06
\(\tau_R\) \(\frac{1}{1.5}\)

このヘッダファイル内で定義されている関数を以下に示す。 各関数の詳細は関数名をクリックしてリンク先を参照のこと。
Functions defined in this header file are listed below. For details of individual functions, click the links.
関数名
Function name 		機能・用途
Purpose
IAPWS95_calculate_thermal_conductivity_lambda0 \(\bar{\lambda}_0(\tau)\)を計算する。
Compute \(\bar{\lambda}_0(\tau)\).
IAPWS95_calculate_thermal_conductivity_lambda1 \(\bar{\lambda}_1(\delta,\tau)\)を計算する。
Compute \(\bar{\lambda}_1(\delta,\tau)\).
IAPWS95_calculate_thermal_conductivity_lambda2 \(\bar{\lambda}_2(\delta,\tau)\)を計算する。
Compute \(\bar{\lambda}_2(\delta,\tau)\).
IAPWS95_calculate_thermal_conductivity Huber et al. (2012)の計算式を用いて 指定された密度・温度における熱伝導率を計算する。
Compute the thermal conductivity at given density and temperature using the equations of Huber et al. (2012).


◆検証 (Validation)

関数IAPWS95_calculate_thermal_conductivity を用いて熱伝導率を計算し、 Huber et al. (2012)の表7, 表8と厳密に一致する結果が得られることを確認した。 但し、ymaeda_opentoolsでは密度\(\rho=0\)での計算はできないので \(\rho=1\times 10^{-10}\) [kg m\(^{-3}\)]での計算で置き換えた。
The thermal conductivities computed by function IAPWS95_calculate_thermal_conductivity were exactly identical to the numerical values shown in Tables 7 and 8 of Huber et al. (2012); the density \(\rho=0\), which is unable to treat by ymaeda_opentools, was replaced by \(\rho=1\times 10^{-10}\) [kg m\(^{-3}\)].

◆引用文献 (References)