나비에-스토크스 방정식

개요

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

\begin{equation} \frac{\uD(m\V)_\mathrm{sys}}{\uD t} = \frac{\partial}{\partial t} \CV \rho\V\ud\mathcal{V} + \CS \V(\rho\V\cdot\mathbf{n})\ud S \end{equation}

여기서 좌변은 검사 체적 내 유체의 운동량 변화율이므로 뉴턴의 운동법칙에 의해 유체가 받는 외력의 총합이 됩니다. 이걸 $\mathbf{F}_\mathrm{ext}$라고 합시다. 이 방정식의 $x$ 방향 성분만 고려하면

\begin{equation} F_{\mathrm{ext},x} = \frac{\partial}{\partial t} \CV \rho u \ud\mathcal{V} + \CS \rho u \V\cdot\mathbf{n}\ud S \end{equation}

발산 정리를 써서 면적분을 선적분으로 바꿔줍니다.

\begin{equation} F_{\mathrm{ext},x} = \CV \left[ \frac{\partial(\rho u)}{\partial t} + \nabla\cdot(\rho u\V) \right] \ud \mathcal{V} \end{equation}

그런데

\begin{align} \frac{\partial(\rho u)}{\partial t} + \nabla\cdot(\rho u\V) &= \rho\frac{\partial u}{\partial t} + \frac{\partial\rho}{\partial t}u + u\nabla\cdot(\rho\V) + \rho\V\cdot\nabla u \nonumber \\ &= \rho\frac{\partial u}{\partial t} + \rho\V\cdot\nabla u + u\cancelto{0}{\left[ \frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\V) \right]} \nonumber \\ &= \rho\left[ \frac{\partial u}{\partial t} + \V\cdot\nabla u \right] \nonumber \\ &= \rho \frac{\uD u}{\uD t} \label{eq:exp1} \end{align}

이고 식 \eqref{eq:exp1}에 대입하면

\begin{equation} F_{\mathrm{ext},x} = \CV \rho\frac{\uD u}{\uD t}\ud\mathcal{V} \end{equation}

$y$와 $z$ 방향에 대해서도 똑같이 해서 각 방향에 대한 식 세 개를 얻고, 셋을 합쳐 벡터 방정식 하나로 만들 수 있습니다.

\begin{equation} \mathbf{F}_\mathrm{ext} = \CV\rho\frac{\uD\V}{\uD t}\ud\mathcal{V} \label{eq:fext} \end{equation}

$\uD\V/\uD t$가 유체 요소의 가속도라는 것을 생각해보면 사실 이 식은 운동 방정식 $F=ma$를 유체에 적요한 것과 같습니다.

유체 요소에 작용하는 힘

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

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

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

체적력

단위 질량당 체적력(즉, 체적력에 의한 가속도)를 $\mathbf{g}$라고 하면 체적력은 자명하게도 아래와 같습니다.

\begin{equation} \mathbf{F}_\mathrm{body} = \CV \rho\mathbf{g}\ud\mathcal{V} \end{equation}

표면력

유체역학이든 고체역학이든 둘 다 연속적인 물질을 다루는 연속체 역학(continuum mechanics)이기 때문에 대부분의 기초 개념을 공유합니다. 잠시 고체역학으로 돌아가서 표면력에 대해 생각해보면, 물체에 작용하는 표면력을 응력(stress)을 사용해서 나타내죠.

단위 면적당 작용하는 표면력을 응력이라고 한다.

그리고 이 응력은 같은 점에서도 면의 방향(단위 법선 벡터 $\mathbf{n}$)에 따라 달라지기 때문에 응력 텐서(stress tensor)를 이용해 방향의 함수로 나타낸다는 점도 고체역학에서 기초적으로 다루는 내용입니다. 단위 법선 벡터가 $\mathbf{n}$인 면에 작용하는 응력을 $\boldsymbol{\tau}(\mathbf{n})$이라 하면

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

여기서 $\tau_{xy}$는 법선 벡터가 $x$ 방향인 면에 작용하는 응력의 $y$ 방향 성분입니다. 나머지도 마찬가지고요. 때로는 간편함을 위해 $\tau_{xy}$ 같이 쓰는 대신 $\tau_{12}$처럼 텐서 표기법을 쓰기 쉽게 표기하기도 합니다.

그 다음으로 나오는 내용은 응력 텐서의 대칭성인데, 과연 유체역학에서도 이게 성립할까요? 각 변이 축에 평행한 직육면체 모양의 유체 요소에 $\tau_{xy}$와 $\tau_{yx}$만 작용하고 나머지 응력 성분은 모두 0이라고 해봅시다.

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

\begin{equation} \tau_{xy}\delta x\delta y\delta z - \tau_{yx}\delta x\delta y\delta z \end{equation}

이고 회전관성은

\begin{equation} \frac{\rho\delta x\delta y\delta z}{12}(\delta x^2+\delta y^2) \end{equation}

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

\begin{equation} (\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 \end{equation}

입니다. 이제 극한 $\delta x\rightarrow0$, $\delta y\rightarrow0$, $\delta z\rightarrow0$을 취하면 $\tau_{xy}=\tau_{yx}$이고, 같은 방법으로 $\tau_{yz}=\tau_{zy}$, $\tau_{zx}=\tau_{xz}$가 됩니다. 즉, 응력 텐서는 대칭 텐서입니다. 단, 외력에 의한 체적토크가 발생할 경우(전자기력이 작용하는 경우가 대표적입니다.) 회전 운동 방정식에 체적토크 항이 들어가게 되어 $\tau_{xy}$와 $\tau_{yx}$가 달라집니다. 보통 공기역학에선 그런 일 없으니 무시해도 괜찮습니다.

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

\begin{equation} \ud \mathbf{F}_\mathrm{surf} = \begin{bmatrix} \tau_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & \tau_{yy} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \tau_{zz} \end{bmatrix} \begin{bmatrix} n_x \\ n_y \\ n_z \end{bmatrix} \ud S \end{equation}

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

\begin{align} F_{\textrm{surface}, x} &= \CS (\tau_{xx}, \tau_{xy}, \tau_{xz}) \cdot \mathbf{n} \ud S \nonumber \\ &= \CV \nabla \cdot (\tau_{xx}, \tau_{xy}, \tau_{xz}) \ud \mathcal V \nonumber \\ &= \CV \left( \frac{\partial\tau_{xx}}{\partial x} +\frac{\partial\tau_{xy}}{\partial y} +\frac{\partial\tau_{xz}}{\partial z} \right) \ud \mathcal V \end{align}

$y$, $z$ 방향도 같은 방법으로 하면 일반적으로 아래 식이 성립함을 알 수 있습니다.

\begin{equation} F_{\mathrm{surf}, i} = \CV \frac{\partial\tau_{ij}}{\partial x_j} \ud\mathcal V \end{equation}

코시 운동량 방정식

최종적으로 식 \eqref{eq:fext}은

\begin{equation} \CV \rho \frac{\uD\V}{\uD t}\ud\mathcal{V} = \mathbf{F}_\mathrm{ext} = \mathbf{F}_\mathrm{body} + \mathbf{F}_\mathrm{surf} \end{equation}

또는 체적력과 표면력 식을 대입하고 텐서 표기법을 사용하면

\begin{gather} \CV \rho \frac{\uD u_i}{\uD t}\ud\mathcal{V} = \CV \rho g_i \ud \mathcal{V} + \CV \frac{\partial\tau_{ij}}{\partial x_j} \ud \mathcal{V} \\ \therefore \rho \frac{\uD u_i}{\uD t}\ud = \rho g_i + \frac{\partial\tau_{ij}}{\partial x_j} \label{eq:cauchy} \end{gather}

식 \eqref{eq:cauchy}는 가장 기본적인 유체 요소의 운동 방정식으로, 특별히 코시 운동량 방정식(Cauchy momentum equation)이라고 부릅니다.

그렇다면, 응력 $\tau_{ij}$는 어떻게 구할까요?

응력과 속도의 관계

먼저 유체에 전단응력(shear stress)이 가해지는 상황을 가정해 봅시다.

두 평판 사이에 유체를 채워넣고 평판을 서로 반대 방향으로 잡아당기면 유체에 전단응력 $\tau$가 걸리게 됩니다. 그럼 전단변형률(shear strain) $\gamma$가 발생하는데, 고체의 경우엔 $\tau$와 $\gamma$ 사이에 관계식을 세우는 것과 달리 유체는 흐를 수 있기 때문에 $\tau$를 주면 $\gamma$가 시간에 따라 점점 커져서 무한대로 변형됩니다. 대신, 유체에서는 $\tau$와 전단변형률의 시간 변화율 $\uD\gamma/\uD t$ 사이에 관계식을 세워야 합니다.

뉴턴 유체

하지만 전단응력과 전단변형률의 시간 변화율 사이의 관계식은 유체마다 다릅니다! 다행히 유체는 대부분 아래와 같은 간단한 성질을 만족하는데, 이들을 뉴턴 유체(Newtonian fluid)라고 합니다.

뉴턴 유체는 다음 성질을 만족하는 유체이다.
  1. 응력은 변형률의 시간변화율에 선형적이다.
  2. 변형률이 변하지 않으면 응력은 수직응력으로 작용하는 압력 뿐이다.
  3. 이 관계는 방향에 관계없이 동일하다(등방성).

따라서, 뉴턴 유체에서 전단응력과 전단변형률의 시간변화율은 정비례 관계이며, 이 비례상수를 점성(viscosity)라 부르고 흔히 $\mu$라고 표기합니다. 즉,

\begin{equation} \tau = \mu \frac{\uD\gamma}{\uD t} \end{equation}

이제 남은 건 $\uD\gamma/\uD t$를 속도로 표현하는 것 뿐입니다.

속도 차이에 의한 전단변형

먼저 간단하게 2차원 유동으로 살펴봅시다. 한 변의 길이가 $h$인 정사각형 모양 유체 요소 $ABCD$가 시간 $\Dt$ 동안 이동하여 사각형 $A’B’C’D’$이 되었습니다. 네 꼭짓점의 속도가 서로 다르다면 $A’B’C’D’$의 모양은 당연히 정사각형이 아닙니다. 위 그림처럼 약간 찌그러진 사각형이 될 것이고, 이는 전단변형률이 존재한다는 얘기죠. 다시 말해, 속도 차이가 전단변형률을 만듭니다.

이때 $A$에서 $y$ 방향 속도가 $v$이면 $B$에서 $y$ 방향 속도는 $v+\frac{\partial v}{\partial x}h$이므로 $B$는 $A$보다 $y$방향으로 $\frac{\partial v}{\partial x}h\Dt$만큼 더 이동합니다. $\overline{AB}$의 길이 변화가 매우 작고 $\Dt$가 매우 작아 $\alpha$ 또한 매우 작다고 가정하면 다음과 같은 근사식이 성립함니다.

\begin{equation} \alpha = \frac{\partial v}{\partial x}\Dt \end{equation}

같은 방법으로

\begin{equation} \beta = \frac{\partial u}{\partial y}\Dt \end{equation}

$\Dt$ 동안 전단변형률의 변화량은 $\alpha + \beta$이므로

\begin{gather} \Delta \gamma_{xy} = \frac{\partial v}{\partial x}\Dt + \frac{\partial u}{\partial y}\Dt \\ \therefore \frac{\uD\gamma_{xy}}{\uD t} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \end{gather}

$yz$ 평면과 $zx$ 평면도 계산하면 일반적으로 다음이 성립함니다.

\begin{equation} \frac{\uD\gamma_{ij}}{\uD t} = \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}, \quad\quad i \neq j \label{eq:strain-rate} \end{equation}

한편 사각형 $ABCD$를 $A’B’C’D’$으로 보내는 변환은 전단변환 말고도 회전변환 또한 포함합니다. $\Dt$ 동안 $z$축을 기준으로 유체 요소가 회전한 각도 $\Delta\theta_z$는 $(\alpha-\beta)/2$로 표현할 수 있습니다. 따라서 $z$ 방향 각속도는

\begin{equation} \frac{\uD\theta_z}{\uD t} = \frac{1}{2}\left( \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \right) \label{eq:z-angvel} \end{equation}

각속도의 나머지 성분도 구하면 다음과 같은 결론이 나옵니다.

\begin{equation} (\mathrm{각속도}) = \frac{1}{2} \nabla \times \V \end{equation}

즉, 유체 요소의 각속도는 속도의 회전의 정확히 절반입니다. 여기서 속도의 회전 $\nabla\times\V$를 와도(vorticity)라 하고 흔히 기호 $\boldsymbol\omega$ (또는 $\boldsymbol\zeta$, $\boldsymbol\xi$)로 표기합니다.

식 \eqref{eq:strain-rate}과 \eqref{eq:z-angvel}를 보면 전단변형률의 시간변화율과 각속도 모두 속도의 기울기에 관련있다는 점을 알 수 있습니다. 수학적으로는 속도의 기울기 $\nabla\V$를 대칭 텐서와 비대칭 텐서로 분해한 것이 각각 전단변형률의 시간변화율과 각속도와 관련있으며 변형률의 시간변화율 텐서(strain rate tensor) $\mathbf{S}$와 와도 텐서(vorticity tensor) $\boldsymbol\Omega$라고 부릅니다.

\begin{gather} (\nabla\V)_{ij} = \frac{\partial u_i}{\partial x_j} \\ \mathbf{S} = \mathrm{sym}(\nabla\V) = \frac{1}{2}(\nabla\V + \nabla\V^T) \\ \boldsymbol\Omega = \mathrm{skw}(\nabla\V) = \frac{1}{2}(\nabla\V - \nabla\V^T) \end{gather}

자명하게 $\mathbf{S}+\boldsymbol\Omega=\nabla\V$입니다. 두 텐서와 전단변형률의 시간변화율 및 각속도 사이에는 아래와 같은 관계가 성립합니다.

\begin{gather} S_{ij} = \frac{1}{2} \frac{\uD\gamma_{ij}}{\uD t}, \quad\quad i \neq j \\[1em] \boldsymbol\Omega = \begin{bmatrix} 0 & \uD\theta_z/\uD t & \uD\theta_y/\uD t \\ -\uD\theta_z/\uD t & 0 & \uD\theta_x/\uD t \\ -\uD\theta_y/\uD t & -\uD\theta_x/\uD t & 0 \end{bmatrix} = \frac{1}{2} \begin{bmatrix} 0 & \omega_z & \omega_y \\ -\omega_z & 0 & \omega_x \\ -\omega_y & -\omega_x & 0 \end{bmatrix} \end{gather}

속도 차이에 의한 수직변형

뉴턴 유체의 정의에 의해 수직응력은 수직변형률의 시간변화율에 선형적입니다. 그리고 변형률이 일정하면 수직응력은 오직 압력 $p$ 뿐이므로, $x$방향 수직응력 $\tau_{xx}$는 다음과 같이 표현할 수 있습니다.

\begin{equation} \tau_{xx} = -p + C_1 \frac{\uD\epsilon_x}{\uD t} + C_2 \frac{\uD\epsilon_y}{\uD t} + C_3 \frac{\uD\epsilon_z}{\uD t} \label{eq:tauxx} \end{equation}

좌표계를 $x$축을 기준으로 90도 회전시켜 새로운 좌표계 $x’y’z’$을 만들어 봅시다. $x’$이 $x$와, $y’$이 $z$와, $z’$이 $y$와 나란합니다. 따라서,

\begin{align} \tau_{xx} &= \tau_{x'x'} \\ \epsilon_x &= \epsilon_{x'} \\ \epsilon_y &= \epsilon_{z'} \\ \epsilon_z &= \epsilon_{y'} \end{align}

식 \eqref{eq:tauxx}에 대입합시다.

\begin{equation} \tau_{x'x'} = -p + C_1 \frac{\uD\epsilon_{x'}}{\uD t} + C_3 \frac{\uD\epsilon_{y'}}{\uD t} + C_2 \frac{\uD\epsilon_{z'}}{\uD t} \label{eq:prime1} \end{equation}

그런데 등방성에 의해 식 \eqref{eq:tauxx}과 똑같은 형태의 관계가 $x’y’z’$ 좌표계에서도 성립해야 합니다.

\begin{equation} \tau_{x'x'} = -p + C_1 \frac{\uD\epsilon_{x'}}{\uD t} + C_2 \frac{\uD\epsilon_{y'}}{\uD t} + C_3 \frac{\uD\epsilon_{z'}}{\uD t} \label{eq:prime2} \end{equation}

식 \eqref{eq:prime1}과 \eqref{eq:prime2}을 비교하면 $C_2=C_3$라는 결론을 얻습니다. 이제 $C_2=C_3=\lambda$, $C_1=A+\lambda$라고 두면 식 \eqref{eq:tauxx}는

\begin{equation} \tau_{xx} = -p + (A+\lambda) \frac{\uD\epsilon_x}{\uD t} + \lambda \frac{\uD\epsilon_y}{\uD t} + \lambda \frac{\uD\epsilon_z}{\uD t} \label{eq:tauxx2} \end{equation}

일단 여기까지 하고 다음 그림을 보죠.

앞과 비슷하지만, 이번엔 직각이등변 삼각형 형태의 유체 요소 $ABC$가 시간 $\Dt$ 동안 이동하여 $A’B’C’$이 되었습니다. 그림에서는 길이 관계를 계산하기 쉽게 $A’$이 $A$와 겹치도록 평행이동했습니다. 여기서 $x$ 방향 수직변형률은 변 $\overline{AB}$의 길이 변형률에 해당하므로, $\Dt$ 동안 $x$ 방향 수직변형률의 변화량은

\begin{equation} \Delta\epsilon_x = \frac{\overline{A'B'} - \overline{AB}}{\overline{AB}} = \frac{\overline{A'B'} - h}{h} \end{equation}

$\overline{A’B’}$의 길이는 피타고라스 정리에 $\Dt$가 매우 작다는 근사를 써서 계산합니다.

\begin{gather} \begin{aligned} \overline{A'B'} &= \sqrt{ \left(h + \frac{\partial u}{\partial x}h\Dt\right)^2 + \left(\frac{\partial v}{\partial x}h\Dt\right)^2 } \\ &\approx h \left( 1 + 2\frac{\partial u}{\partial x}\Dt \right)^{1/2} \\ &\approx h \left( 1 + \frac{\partial u}{\partial x}\Dt \right) \end{aligned} \\[1em] \therefore \Delta\epsilon_x = \frac{\partial u}{\partial x} \Dt \end{gather}

같은 방법으로,

\begin{align} \Delta\epsilon_y &= \frac{\partial v}{\partial y} \Dt \\[5pt] \Delta\epsilon_z &= \frac{\partial w}{\partial z} \Dt \end{align}

따라서,

\begin{align} \frac{\uD\epsilon_x}{\uD t} &= \frac{\partial u}{\partial x} = 2S_{11} \\[5pt] \frac{\uD\epsilon_y}{\uD t} &= \frac{\partial v}{\partial y} = 2S_{22} \\[5pt] \frac{\uD\epsilon_z}{\uD t} &= \frac{\partial w}{\partial z} = 2S_{33} \end{align}

식 \eqref{eq:tauxx2}에 대입합시다.

\begin{equation} \tau_{xx} = -p + (A+\lambda)\frac{\partial u}{\partial x} + \lambda\frac{\partial v}{\partial y} + \lambda\frac{\partial w}{\partial z} = -p + A\frac{\partial u}{\partial x} + \lambda\nabla\cdot\V \end{equation}

등방성을 생각하면 $\tau_{yy}$와 $\tau_{zz}$는 다음과 같아야 합니다.

\begin{align} \tau_{yy} &= -p + A\frac{\partial v}{\partial y} + \lambda\nabla\cdot\V \\ \tau_{zz} &= -p + A\frac{\partial w}{\partial z} + \lambda\nabla\cdot\V \end{align}

남은 건 $A$와 $\lambda$를 결정하는 겁니다. 위 그림에서 좌표계를 $z$축에 대해 45도 돌려서 새 좌표계 $x’y’z’$를 만들어보죠. 그럼 변 $\overline{BC}$는 $y’$ 방향입니다. 그러므로

\begin{gather} \Delta\epsilon_{y'} = \frac{\overline{B'C'} - \overline{BC}}{\overline{BC}} = \frac{\overline{B'C'} - \sqrt{2}h}{\sqrt{2}h} \\[1em] \begin{aligned} \overline{B'C'} &= \sqrt{\left( h+\frac{\partial u}{\partial x}h\Dt-\frac{\partial u}{\partial y}h\Dt \right)^2+ \left( h+\frac{\partial v}{\partial y}h\Dt-\frac{\partial v}{\partial x}h\Dt \right)^2} \\ &\approx h\left( 2 + 2\frac{\partial u}{\partial x}\Dt - 2\frac{\partial u}{\partial y}\Dt + 2\frac{\partial v}{\partial y}\Dt - 2\frac{\partial v}{\partial x}\Dt \right)^{1/2} \\ &\approx \sqrt{2}h\left( 1 + \half\frac{\partial u}{\partial x}\Dt - \half\frac{\partial u}{\partial y}\Dt + \half\frac{\partial v}{\partial y}\Dt - \half\frac{\partial v}{\partial x}\Dt \right) \\ &= \sqrt{2}h\left( 1 + \frac{\Delta\epsilon_x}{2} + \frac{\Delta\epsilon_y}{2} - \frac{\Delta\gamma_{xy}}{2} \right) \end{aligned} \\[1em] \therefore\Delta\epsilon_{y'} = \half(\Delta\epsilon_x + \Delta\epsilon_y - \Delta\gamma_{xy}) \end{gather}

같은 방법으로 계산하면 $\Delta\epsilon_{x’}$은 다음과 같습니다.

\begin{equation} \Delta\epsilon_{x'}=\half(\Delta\epsilon_x + \Delta\epsilon_y + \Delta\gamma_{xy}) \end{equation}

두 식을 빼면

\begin{gather} \Delta\epsilon_{x'} - \Delta\epsilon_{y'} = \Delta\gamma_{xy} \\[5pt] \frac{\uD\gamma_{xy}}{\uD t} = \frac{\partial u'}{\partial x'} - \frac{\partial v'}{\partial y'} \\[5pt] \therefore\tau_{xy} = \mu\frac{\uD\gamma_{xy}}{\uD t} = \mu \left( \frac{\partial u'}{\partial x'} - \frac{\partial v'}{\partial y'} \right) \label{eq:tauxy} \end{gather}

$\tau_{xy}$를 $x’y’z’$ 좌표계로 변환해서 좌표계를 통일합시다. $xyz$ 좌표계에서 $x’y’z’$으로 변환하는 회전변환 행렬은

\begin{equation} \mathbf{Q} = \begin{bmatrix} 1/\sqrt2 & -1/\sqrt2 & 0 \\ 1/\sqrt2 & 1/\sqrt2 & 0 \\ 0 & 0 & 1 \end{bmatrix} \end{equation}

이므로 $\boldsymbol\tau’={\tau_{i’j’}}$를 $\boldsymbol\tau$로 변환하면

\begin{gather} \boldsymbol\tau = \mathbf{Q}\boldsymbol\tau'\mathbf{Q}^T = \begin{bmatrix} (\tau_{x'x'}+\tau_{y'y'}-2\tau_{x'y'})/2 & (\tau_{x'x'}-\tau_{y'y'})/2 & (\tau_{x'z'}-\tau_{y'z'})/\sqrt2 \\ (\tau_{x'x'}-\tau_{y'y'})/2 & (\tau_{x'x'}+\tau_{y'y'}+2\tau_{x'y'})/2 & (\tau_{x'z'}+\tau_{y'z'})/\sqrt2 \\ (\tau_{z'x'}-\tau_{z'y'})/\sqrt2 & (\tau_{z'x'}+\tau_{z'y'})/\sqrt2 & \tau_{z'z'} \end{bmatrix} \\[1em] \therefore\tau_{xy} = \half(\tau_{x'x'}-\tau_{y'y'}) \end{gather}

식 \eqref{eq:tauxy}에 대입합시다.

\begin{gather} \begin{aligned} \mu \left( \frac{\partial u'}{\partial x'} - \frac{\partial v'}{\partial y'} \right) &= \half(\tau_{x'x'}-\tau_{y'y'}) \\ &= \half\left[ \left( -p + A\frac{\partial u'}{\partial x'} + \lambda\nabla\cdot\V' \right) - \left( -p + A\frac{\partial v'}{\partial y'} + \lambda\nabla\cdot\V' \right) \right] \\ &= \half A\left( \frac{\partial u'}{\partial x'} - \frac{\partial v'}{\partial y'} \right) \end{aligned} \\[1em] \therefore A = 2\mu \end{gather}

최종적으로 수직응력은 다음과 같습니다.

\begin{align} \tau_{xx} &= -p + 2\mu\frac{\partial u}{\partial x} + \lambda\nabla\cdot\V \\ \tau_{yy} &= -p + 2\mu\frac{\partial v}{\partial y} + \lambda\nabla\cdot\V \\ \tau_{zz} &= -p + 2\mu\frac{\partial w}{\partial z} + \lambda\nabla\cdot\V \end{align}

$A$만 구하고 $\lambda$는 따로 구하지 않았는데, 사실 $\lambda$는 $\mu$와 별개로 존재하는 물리량입니다. 고체에서 영률 $E$가 수직변형, 전단 탄성 계수 $G$가 전단변형에 연관돼있듯이 뉴턴 유체에서는 $\mu$가 전단변형, $\lambda$가 수직변형에 연관됩니다. 그리고 사실, 압력을 제외하면 위 식은 고체역학의 라메 상수(Lamé parameters)를 사용한 식과 완전히 똑같습니다. 선형적인 물질이라는 점에서 닮은 것이죠.

이 새로운 물리량 $\lambda$는 체적점성(bulk viscosity 또는 volume viscosity)라고 부릅니다. 만약 비압축성 유동이라면 $\nabla\cdot\V=0$이므로 체적점성이 의미가 없고, 압축성 유동이라도 체적점성의 영향이 크지는 않습니다. 굉장히 빠른 압축과 팽창이 일어나는 경우(예를 들어 충격파나 변화가 심한 소리 파동)에는 체적점성의 영향이 커집니다.

체적점성의 값을 추정하는 방법이 있습니다. 세 방향의 수직응력을 모두 더한 뒤 3으로 나눠봅시다.

\begin{equation} -\frac{\tau_{ii}}{3} = p - \left(\frac{2}{3}\mu+\lambda\right)\nabla\cdot\V \end{equation}

좌변은 수직응력의 평균의 음수로 기계적 압력(mechanical pressure)이라고 부릅니다. 우변의 $p$는 우리가 자주 보는, 분자 운동에 의한 열역학적 압력(thermodynamic pressure)입니다. 여기서 풍선을 예로 들어보면, 고무의 탄성력이 풍선 내부 기체에 기계적 압력을 주고 풍선 내부 기체는 운동하면서 생기는 열역학적 압력으로 기계적 압력을 상쇄한다고 볼 수 있습니다. 여기에서 착안하여 조지 스토크스(George Stokes)는 기계적 압력과 기체분자 운동론에 의한 열역학적 압력을 같은 것으로 보자는 스토크스 가설(Stokes’ hypothesis)을 제시합니다. 압축과 팽창이 빠르지 않다면 압축, 팽창이 일어나더라도 기체 분자가 빠르게 이동하여 열역학적 압력과 기계적 압력 사이에 평형이 이루어진다고 가정하는 겁니다.

\begin{equation} \text{Stokes' hypothesis}: \qquad p_m = p \end{equation}

이를 적용하면 $\lambda=-\frac{2}{3}\mu$가 됩니다.

나비에-스토크스 방정식

전단응력과 수직응력 식을 다음과 같이 하나로 합칠 수 있습니다.

\begin{equation} \tau_{ij} = (-p + \lambda\nabla\cdot\V)\delta_{ij} + \mu\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) = (-p + \lambda\nabla\cdot\V)\delta_{ij} + 2\mu S_{ij} \label{eq:stokes_law} \end{equation}

여기서

\begin{equation} \delta_{ij} = \begin{cases} 1, & i = j \\ 0, & i \neq j \end{cases} \end{equation}

는 크로네커 델타입니다. 이제 식 \eqref{eq:stokes_law}을 식 \eqref{eq:cauchy}에 대입하면 유체에 대한 운동량 방정식, 즉 나비에-스토크스 방정식(Navier-Stokes equation)을 얻습니다.

\begin{equation} \begin{aligned} \rho\frac{\uD u_i}{\uD t} &= \rho g_i + \frac{\partial}{\partial x_j} \left[ (-p+\lambda\nabla\cdot\V) \delta_{ij} +\mu \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) \right] \\ &=\rho g_i + \frac{\partial}{\partial x_i} (-p+\lambda\nabla\cdot\V) + \frac{\partial}{\partial x_j}\left[ \mu \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) \right] \end{aligned} \end{equation}

보통은 $\mu$와 $\lambda$의 변화를 무시하여 다음과 같이 벡터 형식으로 나타냅니다.

\begin{gather} \begin{aligned} \rho\frac{\mathrm D u_i}{\mathrm Dt} &=\rho g_i+\frac{\partial}{\partial x_i}\left(-p+\lambda\nabla\cdot\mathbf V\right) +\mu\left(\frac{\partial u_i}{\partial x_j\partial x_j}+\frac{\partial}{\partial x_i}\frac{\partial u_j}{\partial x_j}\right) \\ &=\rho g_i-\frac{\partial p}{\partial x_i}+\lambda\frac{\partial}{\partial x_i} (\nabla\cdot\mathbf V) +\mu\left[\nabla^2 u_i+\frac{\partial}{\partial x_i}(\nabla\cdot\mathbf V)\right] \\ &=\rho g_i-\frac{\partial p}{\partial x_i}+\mu\nabla^2u_i+(\mu+\lambda)\frac{\partial}{\partial x_i}(\nabla\cdot\mathbf V) \end{aligned} \\[1em] \therefore \rho\frac{\mathrm D\mathbf V}{\mathrm Dt} =\rho\mathbf g-\nabla p+\mu\nabla^2\mathbf V+(\mu+\lambda)\nabla(\nabla\cdot\mathbf V) \label{eq:ns_vec} \end{gather}

특수한 경우

비압축성 유동은 $\nabla\cdot\V=0$이므로 식 \eqref{eq:ns_vec}의 우변 제일 마지막 항이 0이 됩니다.

\begin{equation} \rho\frac{\uD\V}{\uD t}=\rho\mathbf{g}-\nabla p+\mu\nabla^2\V \end{equation}

양변을 $\rho$로 나누면

\begin{equation} \frac{\uD\V}{\uD t}=\mathbf{g}-\frac{\nabla p}{\rho}+\frac{\mu}{\rho}\nabla^2\V \end{equation}

여기서 점성을 밀도로 나눈 값을 동점성(kinematic viscosity)이라 하고 $\nu$로 표기합니다. 동점성의 물리적 의미는 뒤에서 설명하겠습니다.

한편 비점성 유동은 점성이 없는 유동입니다. 이 경우 응력은 압력 뿐이므로 $\mu=\lambda=0$을 대입하면

\begin{equation} \rho\frac{\uD\V}{\uD t} = \rho\mathbf{g}-\nabla p \end{equation}

이 방정식은 특별히 오일러 방정식(Euler equation)이라고 부릅니다.