[공기역학] 나비에-스토크스 방정식 (Navier-Stokes Equation) – 2

응력과 속도의 관계

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

두 평판 사이에 유체를 채워넣고 두 평판을 반대 방향으로 잡아당기면 유체에 전단응력 $\tau$가 걸리게 됩니다. 그러면 전단변형률(shear strain) $\gamma$가 발생하는데, 고체의 경우엔 $\tau$가 정해지면 $\gamma$가 정해지지만, 유체는 흐를 수 있기 때문에 똑같이 $\tau$를 줘도 처음에는 $\gamma$가 작다가 시간이 지날수록 점점 커져서 무한대로 변형됩니다. 그래서 고체에서 $\tau=G\gamma$라고 표현하는 것과 달리 전단응력을 전단변형률의 시간 변화율 $\frac{\mathrm D \gamma}{\mathrm D t}$로 나타내야 합니다.

뉴턴 유체

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

  1. $\tau$가 $\frac{\mathrm D \gamma}{\mathrm D t}$에 선형입니다.
  2. $\frac{\mathrm D \gamma}{\mathrm D t}=0$이면 $\tau=0$
  3. 1의 비례상수는 방향에 관계없이 동일합니다(유체의 등방성).

여기서 비례상수를 $\mu$를 유체의 점성(viscosity)이라고 합니다. 따라서 뉴턴 유체는 아래 식을 만족합니다.

\[ \tau = \mu \frac{\mathrm D \gamma}{\mathrm D t} \]

여기서부터는 유체가 뉴턴 유체라 가정하고 전단응력을 수식으로 나타내보겠습니다.

속도 차이에 의한 전단변형

먼저 간단하게 2차원 유동으로 살펴봅시다. 한 변의 길이가 $h$인 정사각형 모양 유체 요소 $ABCD$가 시간 $\Delta t$동안 이동하여 사각형 $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\Delta t$만큼 더 이동합니다. $\overline{AB}$의 길이 변화를 무시하고 $\Delta t$가 아주 작아 $\alpha$ 또한 매우 작다고 하면

\[ \alpha=\frac{\partial v}{\partial x}\Delta t \]

같은 방법으로

\[ \beta=\frac{\partial u}{\partial y}\Delta t \]

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

\[ \Delta \gamma_{xy}=\frac{\partial u}{\partial y}\Delta t+\frac{\partial v}{\partial x}\Delta t \]

따라서

\[ \frac{\mathrm D\gamma_{xy}}{\mathrm Dt}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \]

같은 방법으로 $yz$ 평면과 $zx$ 평면도 계산하면

\[ \frac{\mathrm D\gamma_{ij}}{\mathrm Dt}=\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}, \qquad i \neq j \]

따라서 뉴턴 유체에 대해

\[ \tau_{ij}=\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right), \qquad i\neq j \]

와도(Vorticity)

한편 $\Delta t$동안 $z$축을 기준으로 유체 요소가 회전한 각도 $\Delta\theta_z$는 $(\alpha-\beta)/2$로 표현할 수 있습니다. 따라서 각속도는

\[ \frac{\mathrm D\theta_z}{\mathrm Dt}=\frac{1}{2}\left(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}\right) \]

각속도의 $x$와 $y$성분도 구하면 다음과 같은 결론이 나옵니다.

\[ (\textrm{angular velocity of fluid element})=\frac{1}{2} \nabla\times\mathrm V \]

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

속도 차이에 의한 수직변형

전단변형을 설명할 때 썼던 그림이랑 비슷한 그림입니다. 이번엔 직각 이등변 삼각형 모양의 유체 요소 $ABC$이고, 유체 요소가 시간 $\Delta t$동안 이동하여 $A’B’C’$이 되었습니다. 그림에서는 길이 관계를 계산하기 쉽게 $A’$이 $A$와 겹치도록 평행이동시켰습니다. 그러면 역시 $B$와 $C$의 속도 차이 때문에 $\overline{BC}$의 길이가 달라집니다. 즉,

속도 차이는 수직변형도 발생시킨다.

자, 이제 $\overline{BC}$의 수직변형률(normal strain) 변화량을 구하기 위해 나중 길이 $\overline{B’C’}$을 계산합시다.

\[ \begin{align}
\overline{B’C’}&=\left[ \left( h+\frac{\partial u}{\partial x}h\Delta t-\frac{\partial u}{\partial y}h\Delta t \right)^2+
\left( h+\frac{\partial v}{\partial y}h\Delta t-\frac{\partial v}{\partial x}h\Delta t \right)^2 \right]^{1/2} \\
&\approx\sqrt 2h \left(1+\frac{1}{2}\frac{\partial u}{\partial x}\Delta t-\frac{1}{2}\frac{\partial u}{\partial y}\Delta t
+\frac{1}{2}\frac{\partial v}{\partial y}\Delta t-\frac{1}{2}\frac{\partial v}{\partial x}\Delta t \right) \\
&=\sqrt 2h \left( 1+\frac{\Delta\epsilon_x}{2}+\frac{\Delta\epsilon_y}{2}-\frac{\Delta\gamma_{xy}}{2} \right)
\end{align} \]

마지막에 나오는 $\Delta\epsilon_x$와 $\Delta\epsilon_y$는 위 그림에서 $\overline{AB}$와 $\overline{AC}$의 변형률을 계산하면 구할 수 있습니다.

\[ \Delta\epsilon_x=\frac{\partial u}{\partial x}\Delta t, \qquad
\Delta\epsilon_y=\frac{\partial v}{\partial y}\Delta t \]

여기서 $xy$ 좌표계를 반시계방향으로 45도 회전한 것을 $x’y’$ 좌표계라고 합시다. $\overline{BC}$는 $y’$ 방향이므로 이것의 수직변형률은 $\Delta\epsilon_{y’}$이 됩니다. 따라서,

\[ \Delta\epsilon_{y’} = \frac{\overline{B’C’}-\overline{BC}}{\overline{BC}}
= \frac{1}{2} (\Delta\epsilon_x+\Delta\epsilon_y-\Delta\gamma_{xy}) \]

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

\[ \Delta\epsilon_{x’}=\frac{1}{2}\left(\Delta\epsilon_x+\Delta\epsilon_y+\Delta\gamma_{xy}\right) \]

두 식을 빼면

\[ \Delta\epsilon_{x’}-\Delta\epsilon_{y’}=\Delta\gamma_{xy} \]

다음은 $xy$ 좌표계와 $x’y’$ 좌표계 사이 관계식을 구해서 위 식의 좌표계를 하나로 통일합니다. $\tau_{xx}$, $\tau_{yy}$, $\tau_{xy}$가 주어지면 응력 텐서를 써서 $\tau_{x’x’}$을 계산할 수 있습니다.

\[ \boldsymbol\sigma\left(\frac{1}{\sqrt 2}, \frac{1}{\sqrt 2}, 0\right)
=\left(\frac{\tau_{xx}+\tau_{xy}}{\sqrt 2},\frac{\tau_{xy}+\tau_{yy}}{\sqrt 2},0\right) \]

\[ \therefore \tau_{x’x’}=\boldsymbol\sigma\left(\frac{1}{\sqrt 2}, \frac{1}{\sqrt 2}, 0\right)\cdot\left(\frac{1}{\sqrt 2}, \frac{1}{\sqrt 2}, 0\right)=\frac{\tau_{xx}+\tau_{yy}}{2}+\tau_{xy} \]

같은 방법으로,

\[ \tau_{y’y’}=\frac{\tau_{xx}+\tau_{yy}}{2}-\tau_{xy} \]

두 식을 빼면

\[ \tau_{x’x’}-\tau_{y’y’}=2\tau_{xy}=2\mu\frac{\Delta\gamma_{xy}}{\Delta t} \]

가 됩니다. 마지막은 뉴턴 유체의 성질을 썼습니다. 이제 $\Delta\gamma_{xy}$를 $x’y’$ 좌표계에 대해 나타내었으므로 위 식에 대입합시다.

\[ \frac{\Delta\epsilon_{x’}}{\Delta t}-\frac{\Delta\epsilon_{y’}}{\Delta t}=\frac{\Delta\gamma_{xy}}{\Delta t}=\frac{\tau_{x’x’}-\tau_{y’y’}}{2\mu} \]

\[ \tau_{x’x’}-\tau_{y’y’}=2\mu\left(\frac{\Delta\epsilon_{x’}}{\Delta t}-\frac{\Delta\epsilon_{y’}}{\Delta t}\right) \]

위첨자 $’$을 다 떼버리면

\[ \tau_{xx}-\tau_{yy}=2\mu\left(\frac{\partial u}{\partial x}-\frac{\partial v}{\partial y}\right) \]

마찬가지로 $yz$ 평면과 $xz$ 평면에서 이걸 또 계산하면…

\[ \tau_{yy}-\tau_{zz} = 2\mu\left(\frac{\partial v}{\partial y}-\frac{\partial w}{\partial z}\right) \]

\[ \tau_{zz}-\tau_{xx}=2\mu\left(\frac{\partial w}{\partial z}-\frac{\partial u}{\partial x}\right) \]

이걸 연립하면 수직응력을 구할 수 있을 거 같은데, 순환하는 식이라 아쉽게도 안 됩니다. 여기서 영국의 수학자 조지 스토크스(George Stokes)는 뉴턴 유체의 수직응력에 관한 두 가지 가정을 세웠습니다(정확히 말하자면 이 가정을 만족하는 유체를 뉴턴 유체라 부르는 것에 가깝지만).

  1. 수직응력은 수직변형률의 시간변화율에 선형적이다.
  2. 정지한 유체에서 수직응력은 압력 뿐이다.

생각해보면 합리적인 가정이죠. 이 가정에 따라 수직응력을 다음과 같이 표현할 수 있습니다.

\[ \tau_{xx} = -p + C_1\frac{\partial u}{\partial x} + C_2\frac{\partial v}{\partial y} + C_3\frac{\partial w}{\partial z} \]

잠시 $\frac{\partial w}{\partial z}=0$인 경우를 생각해봅시다. 이때 $x$축을 중심으로 $y$축과 $z$축을 반시계 방향으로 90도 회전하면 새 $z$축은 원래 $y$축이므로 회전 후에는 $\frac{\partial v}{\partial x}=0$이고 원래 $\frac{\partial v}{\partial x}$가 회전 후 $\frac{\partial w}{\partial z}$가 됩니다. 그런데 이렇게 한다고 $\tau_{xx}$가 달라지진 않을 것이고, 뉴턴 유체는 등방성을 가지기 때문에 $C_1$, $C_2$, $C_3$도 그대로입니다. 즉,

\[ \tau_{xx} = -p + C_1\frac{\partial u}{\partial x} + C_2\left(\frac{\partial v}{\partial y}\right)_\text{before}
= -p + C_1\frac{\partial u}{\partial x} + C_3\left(\frac{\partial w}{\partial z}\right)_\text{after}
= -p + C_1\frac{\partial u}{\partial x} + C_3\left(\frac{\partial v}{\partial y}\right)_\text{before} \]

따라서 $C_2=C_3$이고, $\tau_{xx}$를 아래와 같이 조금 다르게 쓸 수 있습니다.

\[ \tau_{xx} = -p + A\frac{\partial u}{\partial x} + \lambda\left(\frac{\partial u}{\partial x}
+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right) = -p + A\frac{\partial u}{\partial x}
+ \lambda\nabla\cdot\mathbf V \]

역시 등방성에 의해 다음도 성립합니다.

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

이제 위에서 구한 식 중 첫 번째 식에 대입합니다.

\[ \begin{align}
\tau_{xx} – \tau_{yy} &= \left(-p + A\frac{\partial u}{\partial x} + \lambda\nabla\cdot\mathbf V\right)
-\left(-p + A\frac{\partial v}{\partial y} + \lambda\nabla\cdot\mathbf V\right) \\
&= A\frac{\partial u}{\partial x} – A\frac{\partial v}{\partial y} \\
&= 2\mu \left(\frac{\partial u}{\partial x} – \frac{\partial v}{\partial y}\right)
\end{align} \]

\[ \therefore A = 2\mu \]

최종적으로

\[ \tau_{ii} = -p + 2\mu\frac{\partial u_i}{\partial x_i} + \lambda\nabla\cdot\mathbf V \]

가 됩니다. 여기서 $\lambda$는 $\mu$랑은 별개로 존재하는 점성과 관련된 물리량으로, 체적점성(bulk viscosity 또는 volume viscosity)이라고 합니다. (어째 책마다 체적점성이 가리키는 대상이 다른데, 일단 여기서는 White의 Viscous Fluid Flow를 따르겠습니다.)

그리고 고체역학을 배우면 알겠지만, 이 결과는 (압력을 제외하면) 라메 상수(Lamé parameters)를 써서 응력을 변형률로 나타낸 식과 완전히 똑같습니다. 뉴턴 유체가 선형성을 가진 물질이니 필연적으로 같을 수 밖에 없죠. 또 고체에서 수직변형에 영률 $E$, 전단변형에 전단 탄성 계수 $G$를 쓰듯이 유체에서도 수직변형에 체적점성 $\lambda$, 전단변형에 점성 $\mu$를 써서 운동을 설명해야 합니다. 그러니, 갑자기 $\lambda$가 등장하는 게 전혀 이상할 것 없습니다.

스토크스 가설

세 방향의 수직응력을 모두 더한 뒤 3으로 나눠봅시다.

\[ -\frac{\tau_{xx}+\tau_{yy}+\tau_{zz}}{3}
=p-\left(\frac{2}{3}\mu+\lambda\right)\nabla\cdot\mathbf V \]

좌변은 수직응력의 평균값에 마이너스를 붙인 것으로, 기계적 압력(mechanical pressure) $p_m$이라고 합니다. 반면에 우변의 $p$는 분자 운동에 의한 열역학적 압력(thermodynamic pressure)입니다. 여기서 풍선을 예로 들어보면, 고무의 탄성력이 풍선 내부 기체에 기계적 압력을 주고 풍선 내부 기체는 운동하면서 생기는 열역학적 압력으로 기계적 압력을 상쇄한다고 볼 수 있습니다. 여기에서 착안하여 스토크스는 기계적 압력과 기체분자 운동론에 의한 열역학적 압력을 같은 것으로 보자는 스토크스 가설(Stokes’ hypothesis)을 제시합니다.

\[ \text{Stokes’ Hypothesis:} \qquad p_m = p \]

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

나비에-스토크스 방정식

전단응력 식과 수직응력 식을 합하면

\[ \tau_{ij}=(-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) \]

여기서

\[ \delta_{ij}=\begin{cases}
1 & \textrm{if }i=j \\
0 & \textrm{if }i\neq j
\end{cases} \]

는 크로네커 델타입니다. 이제 $\tau_{ij}$를 모두 계산했으므로 준비가 끝났습니다.

\[ \begin{align}
\rho\frac{\mathrm D u_i}{\mathrm Dt}
&=\rho g_i+\frac{\partial}{\partial x_j} \left[\left(-p+\lambda\nabla\cdot\mathbf V\right)\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_j}\left(-p+\lambda\nabla\cdot\mathbf V\right)\delta_{ij}
+\mu\frac{\partial}{\partial x_j}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right) \\
&=\rho g_i+\frac{\partial}{\partial x_i}\left(-p+\lambda\nabla\cdot\mathbf V\right)
+\mu\left(\nabla^2u_i+\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{align} \]

따라서

\[ \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) \]

이를 나비에-스토크스 방정식(Navier-Stokes equation) 또는 간단히 운동량 방정식(momentum equation)이라고 합니다. 좌변은 유체 요소의 밀도와 가속도의 곱, 우변 첫 항은 외력, 우변 둘째 항은 압력차, 셋째와 넷째 항은 점성에 의한 마찰력을 의미합니다.

비압축성 유동의 경우 우변 넷째 항이 0이 되어

\[ \rho\frac{\mathrm D\mathbf V}{\mathrm Dt}=\rho\mathbf g-\nabla p+\mu\nabla^2\mathbf V \]

이고 비점성 유동은 우변 셋째와 넷째 항이 0이 되어

\[ \rho\frac{\mathrm D\mathbf V}{\mathrm Dt}=\rho\mathbf g-\nabla p \]

인데 이를 특별히 오일러 방정식(Euler equation)이라고 부릅니다. 정지 상태의 유체인 경우 속도가 없으므로

\[ \nabla p=\rho\mathbf g \]

라는 익히 알려진 결과가 나옵니다.