[공기역학] 에너지 방정식 (Energy Equation) 유도

열역학 제1법칙은 에너지 보존법칙으로서, 열에너지와 일을 받은 만큼 계의 내부 에너지가 변한다는 법칙입니다. 이걸 운동하는 유체에서 말하면, 열에너지와 일을 받은 만큼 내부 에너지와 운동 에너지의 합이 변합니다. 즉,

\[ \frac{\mathrm{D}}{\mathrm{D}t}\left(E+\frac{1}{2}mV^2\right)=\dot{Q}+\dot{W} \]

여기서 $E$는 유체 요소의 총 내부 에너지입니다. 따라서 레이놀즈 수송정리에 $B=E+\frac{1}{2}mV^2$, $\beta=e+\frac{V^2}{2}$을 넣으면($e$는 단위 질량당 내부 에너지)

\[ \frac{\mathrm{D}}{\mathrm{D}t}\left(E+\frac{1}{2}mV^2\right)_\textrm{sys}
=\frac{\partial}{\partial t}\iiint\limits_\textrm{CV}\rho\left(e+\frac{V^2}{2}\right)\mathrm{d}\mathcal{V}
+\iint\limits_\textrm{CS}\left(e+\frac{V^2}{2}\right)\rho\mathbf{V}\cdot\mathbf{n}\mathrm{d}S
=\dot{Q}+\dot{W} \]

항상 그랬듯이 발산 정리를 써서 부피적분으로 통일합니다.

\[ \begin{align}
& \iiint\limits_\textrm{CV}\left\{\frac{\partial}{\partial t}\left[\rho\left(e+\frac{V^2}{2}\right)\right]
+\nabla\cdot\left[\rho\left(e+\frac{V^2}{2}\right)\mathbf{V}\right]\right\}\mathrm{d}\mathcal{V} \\
&= \iiint\limits_\textrm{CV}\left[\rho\frac{\partial}{\partial t}\left(e+\frac{V^2}{2}\right)
+\frac{\partial\rho}{\partial t}\left(e+\frac{V^2}{2}\right)+\left(e+\frac{V^2}{2}\right)\nabla\cdot(\rho\mathbf{V})
+\rho\mathbf{V}\cdot\nabla\left(e+\frac{V^2}{2}\right)\right]\mathrm{d}\mathcal{V} \\
&= \iiint\limits_\textrm{CV}\bigg\{\rho\bigg(\underbrace{\frac{\partial}{\partial t}+\mathbf{V}\cdot\nabla}
_{\mathrm{D}/\mathrm{D}t}\bigg)\left(e+\frac{V^2}{2}\right)
+\left(e+\frac{V^2}{2}\right)\bigg[\underbrace{\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{V})}_{=0\ \textrm{(continuity)}}\bigg]\bigg\}\mathrm{d}\mathcal{V} \\
&= \iiint\limits_\textrm{CV}\rho\frac{\mathrm{D}}{\mathrm{D}t}\left(e+\frac{V^2}{2}\right)\mathrm{d}\mathcal{V} \\
&= \dot{Q}+\dot{W}
\end{align} \]

유체가 받는 열에너지

푸리에 해석으로 유명한 과학자 조제프 푸리에(Jean-Baptiste Joseph Fourier)는 열에너지가 어떻게 이동하는지에 대한 푸리에의 법칙(Fourier’s law)을 발표했습니다.

  1. 열에너지는 온도 기울기의 반대 방향으로 이동한다.
  2. 온도 기울기에 수직으로 놓인 단위 면적을 통과해서 단위 시간동안 흐르는 열에너지의 양(열유속, heat flux)은 온도 기울기에 비례한다.

즉,

\[ \mathbf{q}=-k\nabla T \]

이때 비례상수 $k$를 열전도율(thermal conductivity)이라고 합니다. 만약 단위 면적이 온도 기울기에 수직이 아니라면 단위 법선 벡터를 내적하면 됩니다.

\[ q = -k\nabla T\cdot\mathbf{n} \]

따라서 단위 시간동안 검사 체적 내로 들어오는 열에너지는 이걸 면적분 하면 되겠죠.

\[ \dot{Q}=-\iint\limits_\textrm{CS}q\mathrm{d}S
=\iint\limits_\textrm{CS}k\nabla T\cdot\mathbf{n}\mathrm{d}S
=\iiint\limits_\textrm{CV}\nabla\cdot(k\nabla T)\mathrm{d}\mathcal{V} \]

앞에 음수가 붙는 이유는 나가는 열에너지가 아니라 들어오는 열에너지기 때문입니다.

유체가 받는 일

단위 시간당 유체가 받는 일, 즉 일률은 유명한 공식 $\mathbf{F}\cdot\mathbf{V}$를 써서 구할 수 있습니다. 먼저 체적력은

\[ \dot{W}_\textrm{body}=\iiint\limits_\textrm{CV}\rho\mathbf{g}\cdot\mathbf{V}\mathrm{d}\mathcal{V} \]

표면력은

\[ \begin{align}
\dot{W}_\textrm{surface}&=\iint\limits_\textrm{CS}\mathrm{d}\mathbf{F}_\textrm{surface}\cdot\mathbf{V} \\
&=\iint\limits_\textrm{CS}u_i\mathrm{d}\mathbf{F}_{\textrm{surface},\ i} \\
&=\iint\limits_\textrm{CS}u_i(\tau_{xx_i},\tau_{yx_i},\tau_{zx_i})\cdot\mathbf{n}\mathrm{d}S \\
&=\iiint\limits_\textrm{CV}\nabla\cdot\left[u_i(\tau_{xx_i},\tau_{yx_i},\tau_{zx_i})\right]\mathrm{d}\mathcal{V} \\
&=\iiint\limits_\textrm{CV}\frac{\partial(u_i\tau_{ji})}{\partial x_j}\mathrm{d}\mathcal{V} \\
&=\iiint\limits_\textrm{CV}\left[u_i\frac{\partial\tau_{ji}}{\partial x_j}+\frac{\partial u_i}{\partial x_j}\tau_{ji}\right]\mathrm{d}\mathcal{V}
\end{align} \]

이 식과 열에너지 식을 제일 처음 식에 대입하면

\[ \rho\frac{\mathrm{D}}{\mathrm{D}t}\left(e+\frac{V^2}{2}\right)
=\nabla\cdot(k\nabla T)+\rho\mathbf{g}\cdot\mathbf{V}
+u_i\frac{\partial\tau_{ji}}{\partial x_j}+\frac{\partial u_i}{\partial x_j}\tau_{ji} \]

그런데 나비에-스토크스 정리를 유도할 때

\[ \rho\frac{\mathrm{D}u_i}{\mathrm{D}t}=\rho g_i+\frac{\partial\tau_{ji}}{\partial x_j} \]

에서

\[ u_i\left(\rho g_i+\frac{\partial\tau_{ji}}{\partial x_j}\right)
=\rho\mathbf{g}\cdot\mathbf{V}+u_i\frac{\partial\tau_{ji}}{\partial x_j}
=\rho u_i\frac{\mathrm{D}u_i}{\mathrm{D}t}
=\rho\frac{\mathrm{D}}{\mathrm{D}t}\left(\frac{V^2}{2}\right) \]

이므로

\[ \rho\frac{\mathrm{D}e}{\mathrm{D}t}=\nabla\cdot(k\nabla T)+\frac{\partial u_i}{\partial x_j}\tau_{ji} \]

$\tau_{ji}$를 대입해서 우변 마지막 항을 계산합니다.

\[ \begin{align}
\frac{\partial u_i}{\partial x_j}\tau_{ji}
&=\frac{\partial u_i}{\partial x_j}\left[(-p+\lambda\nabla\cdot\mathbf{V})\delta_{ij}
+\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\right] \\
&=-p\nabla\cdot\mathbf{V}+\lambda(\nabla\cdot\mathbf{V})^2+\mu\frac{\partial u_i}{\partial x_j}
\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)
\end{align} \]

여기서 제일 오른쪽 두 항을 점성 소산(viscous dissipation)이라고 하며, 점성(=마찰)으로 인해 생기는 열에너지를 나타냅니다. 자세히 식을 풀어보자면,

\[ \begin{align}
\Phi&=\mu\left[2\left(\frac{\partial u}{\partial x}\right)^2+2\left(\frac{\partial v}{\partial y}\right)^2
+2\left(\frac{\partial w}{\partial z}\right)^2+\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2
+\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)^2
+\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)^2\right] \\
&\qquad+\lambda\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)^2 \\
&=\mu\left[\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2
+\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)^2
+\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)^2\right] \\
&\qquad+\frac{2}{3}\mu\left[\left(\frac{\partial u}{\partial x}-\frac{\partial v}{\partial y}\right)^2
+\left(\frac{\partial v}{\partial y}-\frac{\partial w}{\partial z}\right)^2
+\left(\frac{\partial w}{\partial z}-\frac{\partial u}{\partial x}\right)^2\right] \\
&\qquad+\left(\lambda+\frac{2}{3}\mu\right)\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)^2
\end{align} \]

$\mu\ge0$이므로 $\Phi\ge0$이기 위해선(마찰열이 음수일 순 없습니다) $\lambda\ge-\frac{2}{3}\mu$여야 한다는 계산이 나옵니다. 스토크스 가설은 가능한 $\lambda$의 최솟값인 셈이죠. 원래 식으로 넘어가서,

\[ \rho\frac{\mathrm{D}e}{\mathrm{D}t}=\nabla\cdot(k\nabla T)-p\nabla\cdot\mathbf{V}+\Phi \]

여기서 속도의 발산을 밀도의 물질 도함수로 나타내는 식을 대입합니다.

\[ p\nabla\cdot\mathbf{V}=p\left(-\frac{1}{\rho}\frac{\mathrm{D}\rho}{\mathrm{D}t}\right)
=\rho\frac{\mathrm{D}}{\mathrm{D}t}\left(\frac{p}{\rho}\right)-\frac{\mathrm{D}p}{\mathrm{D}t} \]

따라서

\[ \rho\frac{\mathrm{D}e}{\mathrm{D}t}=\nabla\cdot(k\nabla T)-\rho\frac{\mathrm{D}}{\mathrm{D}t}\left(\frac{p}{\rho}\right)
+\frac{\mathrm{D}p}{\mathrm{D}t}+\Phi \]

또는

\[ \rho\frac{\mathrm{D}}{\mathrm{D}t}\left(e+\frac{p}{\rho}\right)=\nabla\cdot(k\nabla T)+\frac{\mathrm{D}p}{\mathrm{D}t}+\Phi \]

여기서 $e+p/\rho$를 단위 질량당 엔탈피(enthalpy) $h$로 정의합니다.

\[ \rho\frac{\mathrm{D}h}{\mathrm{D}t}=\nabla\cdot(k\nabla T)+\frac{\mathrm{D}p}{\mathrm{D}t}+\Phi \]

이 식이 바로 에너지 방정식입니다.