[공기역학] 운동량 방정식 / 나비에-스토크스 방정식(Navier-Stokes Equation) 유도 – 1

레이놀즈 수송정리에 운동량을 대입해봅시다. ($B=m\mathbf V$, $\beta=\mathbf V$)

\[ \frac{\mathrm D (m\mathbf V)_\textrm{sys}}{\mathrm D t} = \frac{\partial}{\partial t} \iiint \limits_\textrm{CV} \rho \mathbf V \mathrm d \mathcal V + \iint \limits_\textrm{CS} \mathbf V (\rho \mathbf V \cdot \mathbf n) \mathrm d S \]

여기서 좌변은 유체 요소의 운동량 변화율이므로 유체 요소에 작용하는 외력 $\mathbf F_\textrm{ext}$와 동일합니다. 이 방정식의 $x$ 방향 성분만 고려하면

\[ F_{\textrm{ext}, x} = \frac{\partial}{\partial t} \iiint \limits_\textrm{CV} \rho u \mathrm d \mathcal V + \iint \limits_\textrm{CS} \rho u \mathbf V \cdot \mathbf n \mathrm d S \]

여기서 $\mathbf V = (u, v, w)$입니다. 이제 발산 정리를 써서 면적분을 선적분으로 바꿔줍시다.

\[ F_{\textrm{ext}, x} = \iiint \limits_\textrm{CV} \left[ \frac{\partial (\rho u)}{\partial t} + \nabla \cdot (\rho u \mathbf V) \right] \mathrm d \mathcal V \]

그런데

\[ \begin{align} \frac{\partial (\rho u)}{\partial t} + \nabla \cdot (\rho u \mathbf V) &= \rho \frac{\partial u}{\partial t} + \frac{\partial \rho}{\partial t} u + u \nabla \cdot (\rho \mathbf V) + \rho \mathbf V \cdot \nabla u \\ &= \rho \frac{\partial u}{\partial t} + \rho \mathbf V \cdot \nabla u + u \underbrace{\left[ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf V) \right]}_{= 0 \ (\textrm{continuity equation})} \\ &= \rho \left[ \frac{\partial u}{\partial t} + \mathbf V \cdot \nabla u \right] \\ &= \rho \frac{\mathrm D u}{\mathrm D t} \end{align} \]

이고 위 식에 대입하면 

\[ F_{\textrm{ext}, x} = \iiint \limits_\textrm{CV} \rho \frac{\mathrm D u}{\mathrm D t} \mathrm d \mathcal V \]

$y$와 $z$ 방향에 대해서도 똑같이 해주면 각 방향에 대한 스칼라 식이 나오고, 셋을 합쳐 다음과 같은 벡터 방정식 하나로 만들 수 있습니다.

\[ \mathbf F_\textrm{ext} = \iiint \limits_\textrm{CV} \rho \frac{\mathrm D \mathbf V}{\mathrm D t} \mathrm d \mathcal V \]

$\frac{\mathrm D \mathbf V}{\mathrm D t}$가 유체 요소의 가속도라는 것을 생각해보면 사실 이 식은 운동 방정식 $F=ma$을 유체역학에 적용한 것이 됩니다.

유체 요소에 작용하는 힘

이제 좌변의 값을 구해야 합니다. 유체 요소에 작용하는 힘은 크게 두 가지로 나눌 수 있습니다.

  • 체적력(body force) 부피에 작용하는 힘입니다. 물체 외부에서 물체 내부로 파고 들어 원격으로 작용하는 힘 정도로 생각할 수 있겠습니다. 중력이나 전자기력이 이에 해당합니다.
  • 표면력(surface force) 표면에 작용하는 힘으로, 압력 등이 이에 해당합니다.

체적력과 표면력의 가장 큰 차이는 물체의 변형에서 나타납니다. 길쭉한 원통이 있을 때 길이의 한쪽 방향으로는 중력을 가하고, 반대 방향으로 정확히 똑같은 크기의 전자기력을 가하면 두 힘이 상쇄되어 원통은 움직이거나 변형되지 않습니다. 하지만 원통의 양 밑면에 같은 크기의 인장력을 가하면 두 힘이 크기가 같고 방향이 반대일지라도 원통의 길이가 늘어나게 됩니다.

단위 질량당 체적력(즉, 체적력에 의한 가속도)을 $\mathbf g$라고 하면(보통 공기역학에서 체적력은 중력 뿐이기 때문에 보통은 중력 가속도의 기호를 그대로 씁니다) 체적력은

\[ \mathbf F_\textrm{body} = \iiint \limits_\textrm{CV} \rho \mathbf g \mathrm d \mathcal V \]

문제는 표면력인데, 심히 복잡합니다.

응력 텐서

유체 속에 법선 벡터가 $\mathbf n$, 넓이가 $\delta S$인 아주 작은 면이 놓여 있고, 이 면에 작용하는 표면력이 $\delta \mathbf F$일 때 면에 작용하는 응력(stress)를 다음과 같이 정의합니다. (여기서 면은 양 쪽이 아니라 한 쪽만 생각합니다.)

\[ \boldsymbol \tau(\mathbf n) = \lim_{\delta S \rightarrow 0} \frac{\delta \mathbf F}{\delta S} \]

면의 방향이 반대가 되면 응력의 방향도 반대가 됩니다.

증명) 기둥 모양의 아주 작은 유체 요소를 생각해봅시다.

양 밑면의 법선 벡터는 각각 $\mathbf n$과 $-\mathbf n$이고 밑면의 넓이는 $\delta S$, 높이는 $\delta h$입니다. 운동방정식을 세우면

\[ \boldsymbol \tau(\mathbf n) \delta S + \boldsymbol \tau(-\mathbf n) \delta S + \mathbf F_\text{side} + \rho \delta S \delta h \mathbf g = \rho \delta S \delta h \frac{\mathrm D \mathbf V}{\mathrm D t} \]

$\delta h \rightarrow 0$을 취하면 $\mathbf F_\text{side}$(옆면에 작용하는 힘), $\rho\delta S\delta h\mathbf g$, 우변이 0이 되고 양 밑면이 서로 맞붙습니다. 따라서 다음이 성립합니다. ■

\[ \boldsymbol \tau(-\mathbf n) = – \boldsymbol \tau(\mathbf n) \]

이제 응력의 좌표 변환을 구해봅시다. 아래와 같은 삼각뿔 모양 유체 요소에 작용하는 알짜 표면력은

\[ \boldsymbol\tau(-\mathbf i)\delta A_x+\boldsymbol\tau(-\mathbf j)\delta A_y+\boldsymbol\tau(-\mathbf k)\delta A_z +\boldsymbol\tau(\mathbf n)\delta A = \boldsymbol\tau(\mathbf n)\delta A-\boldsymbol\tau(\mathbf i)\delta A_x -\boldsymbol\tau(\mathbf j)\delta A_y-\boldsymbol\tau(\mathbf k)\delta A_z \]

여기서 $\delta A_x$와 $\delta A$의 이면각은 $\mathbf i$와 $\mathbf n=(n_x, n_y, n_z)$ 사이의 각과 같으므로 그 코사인 값은 $n_x$이고(방향 코사인) 따라서

\[ \delta A_x = n_x \delta A \]

마찬가지로

\[ \begin{align} \delta A_y &= n_y \delta A \\ \delta A_z &= n_z \delta A \end{align} \]

위 식에 대입하면 알짜 표면력은

\[ \left[ \boldsymbol\tau(\mathbf n)-\boldsymbol\tau(\mathbf i)n_x-\boldsymbol\tau(\mathbf j)n_y-\boldsymbol\tau(\mathbf k)n_z \right] \delta A \]

한편 유체 요소의 부피는

\[ \sqrt{\frac{2}{9} \delta A_x \delta A_y \delta A_z} = \sqrt{\frac{2}{9} n_x n_y n_z \delta A^3} \]

이므로(직접 계산해봅시다) 운동방정식은 아래와 같습니다.

\[ \left[ \boldsymbol\tau(\mathbf n)-\boldsymbol\tau(\mathbf i)n_x-\boldsymbol\tau(\mathbf j)n_y-\boldsymbol\tau(\mathbf k)n_z \right] \delta A = \rho \sqrt{\frac{2}{9} n_x n_y n_z \delta A^3} \left( \frac{\mathrm D \mathbf V}{\mathrm D t} – \mathbf g \right) \]

$\delta A \rightarrow 0$을 취하면

\[ \boldsymbol\sigma(\mathbf n) = \boldsymbol\sigma(\mathbf i)n_x+\boldsymbol\tau(\mathbf j)n_y+\boldsymbol\tau(\mathbf k)n_z \]

$\boldsymbol\sigma(\mathbf i)$, $\boldsymbol\tau(\mathbf j)$, $\boldsymbol\tau(\mathbf k)$를 $x$, $y$, $z$ 성분으로 다음과 같이 분해합시다.

\[ \begin{align}
\boldsymbol\tau(\mathbf i) &= \tau_{xx} \mathbf i + \tau_{xy} \mathbf j + \tau_{xz} \mathbf k \\
\boldsymbol\tau(\mathbf j) &= \tau_{yx} \mathbf i + \tau_{yy} \mathbf j + \tau_{yz} \mathbf k \\
\boldsymbol\tau(\mathbf k) &= \tau_{zx} \mathbf i + \tau_{zy} \mathbf j + \tau_{zz} \mathbf k
\end{align} \]

$\tau_{xy}$의 의미는 법선 벡터가 $x$ 방향인 면에 작용하는 응력의 $y$방향 성분입니다(나머지도 마찬가지). 때로는 간편함을 위해 $\tau_{xy}$ 대신 $\tau_{12}$를 쓰기도 합니다(이렇게 표기하면 $\tau_{ij} \text{ for } i=1,2,3, \ j=1,2,3$같이 쓸 수 있습니다). 이걸 대입하면

\[ \boldsymbol\tau(\mathbf n) = \begin{bmatrix}
\tau_{xx}n_x + \tau_{yx}n_y + \tau_{zx}n_z \\
\tau_{xy}n_x + \tau_{yy}n_y + \tau_{zy}n_z \\
\tau_{xz}n_x + \tau_{yz}n_y + \tau_{zz}n_z
\end{bmatrix} = \begin{bmatrix}
\tau_{xx} & \tau_{yx} & \tau_{zx} \\
\tau_{xy} & \tau_{yy} & \tau_{zy} \\
\tau_{xz} & \tau_{yz} & \tau_{zz}
\end{bmatrix} \begin{bmatrix} n_x \\ n_y \\ n_z \end{bmatrix} \]

이때 $3\times 3$ 행렬

\[ \begin{bmatrix} \tau_{xx} & \tau_{yx} & \tau_{zx} \\ \tau_{xy} & \tau_{yy} & \tau_{zy} \\ \tau_{xz} & \tau_{yz} & \tau_{zz} \end{bmatrix} \]

응력 텐서(stress tensor)라고 합니다.

응력 텐서의 대칭성

각 변이 축에 평행한 직육면체 모양의 유체 요소에 $\tau_{xy}$와 $\tau_{yx}$가 작용한다고 해봅시다.

직육면체 중심을 지나고 $z$축에 평행한 직선에 대한 알짜토크는

\[ \tau_{xy}\delta x\delta y\delta z-\tau_{yx}\delta x\delta y\delta z \]

이고 회전관성은

\[ \frac{\rho\delta x\delta y\delta z}{12} (\delta x^2+\delta y^2) \]

이므로 회전 운동 방정식을 세우면

\[ (\tau_{xy}-\tau_{yx})\delta x\delta y\delta z=\frac{\rho\delta x\delta y\delta z}{12}(\delta x^2+\delta y^2)\alpha \]

입니다. 역시나 $\delta x\rightarrow 0$, $\delta y\rightarrow 0$, $\delta z\rightarrow 0$을 취하면 $\tau_{xy}=\tau_{yx}$이고, 같은 방법으로 $\tau_{yz}=\tau_{zy}$, $\tau_{zx}=\tau_{xz}$가 됩니다. 따라서 응력 텐서는 대칭 행렬입니다.

단, 외력에 의한 체적토크가 발생할 경우(전자기력이 작용하는 경우가 대표적입니다) 회전 운동 방정식에 체적토크 항이 들어가게 되어 $\tau_{xy}$와 $\tau_{yx}$가 달라집니다. 공기역학에선 그런 일 없으니 무시해도 괜찮습니다.

이제 본격적으로 표면력을 구해봅시다. 검사체적의 아주 작은 표면 $\mathrm d S$에 작용하는 표면력은

\[ \mathrm d \mathbf F_\textrm{surface}=\begin{bmatrix} \tau_{xx} & \tau_{yx} & \tau_{zx} \\ \tau_{xy} & \tau_{yy} & \tau_{zy} \\ \tau_{xz} & \tau_{yz} & \tau_{zz} \end{bmatrix} \begin{bmatrix} n_x \\ n_y \\ n_z \end{bmatrix} \mathrm d S \]

$x$방향 성분만 떼어다 면적분하면

\[ F_{\textrm{surface}, x}=\iint \limits_\textrm{CS} (\tau_{xx}, \tau_{yx}, \tau_{zx}) \cdot \mathbf n \mathrm d S =\iiint\limits_\textrm{CV}\nabla\cdot(\tau_{xx}, \tau_{yx}, \tau_{zx})\mathrm d \mathcal V =\iiint\limits_\textrm{CV}\left(\frac{\partial\tau_{xx}}{\partial x}+\frac{\partial\tau_{yx}}{\partial y}+\frac{\partial\tau_{zx}}{\partial z}\right) \mathrm d \mathcal V \]

$y$, $z$방향도 같은 방법으로 하면 일반적으로 아래 식이 성립합니다.

\[ F_{\textrm{surface}, i}=\iiint\limits_\textrm{CV}\frac{\partial\tau_{ji}}{\partial x_j}\mathrm d\mathcal V \]

가 됩니다. (텐서 표기법을 사용했습니다.) 최종적으로

\[ \iiint \limits_\textrm{CV} \rho \frac{\mathrm D u_i}{\mathrm D t} \mathrm d\mathcal V =\iiint \limits_\textrm{CV} \rho g_i \mathrm d\mathcal V +\iiint \limits_\textrm{CV}\frac{\partial\tau_{ji}}{\partial x_j}\mathrm d\mathcal V \]

또는

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

이게 가장 기본적인 유체 요소의 운동 방정식이 되고, 특별히 코시 운동량 방정식(Cauchy momentum equation)이라고 부릅니다. 이제 $\tau_{ji}$의 값을 구해봅시다.