Fractional Step Method를 이용한 시간 이산화
비압축성 운동량 방정식
\begin{equation}
\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u}
\end{equation}
을 대류항 $\N = \mathbf{u} \cdot \nabla \mathbf{u}$는 Adam-Bashforth 방법으로, 점성항 $\nu \nabla^2 \mathbf{u}$는는 Crank-Nicolson 방법으로 시간에 대해 이산화해봅시다.
\begin{equation}
\frac{\u^{n+1} - \u^n}{\Delta t} + \left( \frac{3}{2} \N^n - \frac{1}{2} \N^{n-1} \right) = -\frac{1}{\rho} \nabla p^{n+1/2} + \frac{\nu}{2} \left( \nabla^2 \u^{n+1} + \nabla^2 \u^n \right)
\label{eq:discretized-momentum}
\end{equation}
위와 같은 이산화는 시각 $n+\frac{1}{2}$에서 시간에 대해 2차 정확도를 가집니다.
\begin{align*}
\frac{\partial\mathbf{u}}{\partial t}^{n+1/2} &= \frac{\u^{n+1} - \u^n}{\Delta t} + O(\Delta t^2) \\
\N^{n+1/2} &= \frac{3}{2} \N^n - \frac{1}{2} \N^{n-1} + O(\Delta t^2) \\
\nabla^2 \mathbf{u}^{n+1/2} &= \frac{1}{2} \left( \nabla^2 \u^{n+1} + \nabla^2 \u^n \right) + O(\Delta t^2)
\end{align*}
이제 여느 CFD 알고리즘이 그러하듯이 식 \eqref{eq:discretized-momentum}에서 압력을 분리(segregate)해봅시다. 식 \eqref{eq:discretized-momentum}의 문제는 모르는 변수가 $\u^{n+1}$과 $p^{n+1/2}$ 두 개라서 이 식만 가지고는 해를 구할 수 없다는 것인데, fractional step method(FSM)는 모르는 $p^{n+1/2}$ 대신에 이미 알려져 있는 스칼라장 $q$를 대입하여 압력을 분리해냅니다. 하지만 이렇게 하면 $\u^{n+1}$이 이산화된 운동량 방정식을 만족하지 않겠죠. 새롭게 바뀐 방정식의 해를 중간 속도(intermediate velocity) $\u^*$라고 합시다. 즉, 중간 속도는 다음 방정식을 만족합니다.
\begin{equation}
\frac{\u^* - \u^n}{\Delta t} + \left( \frac{3}{2} \N^n - \frac{1}{2} \N^{n-1} \right) = -\frac{1}{\rho} \nabla q + \frac{\nu}{2} \left( \nabla^2 \u^* + \nabla^2 \u^n \right)
\label{eq:intermediate-velocity}
\end{equation}
식 \eqref{eq:intermediate-velocity}을 식 \eqref{eq:discretized-momentum}에서 빼면
\begin{equation*}
\frac{\u^{n+1} - \u^*}{\Delta t} = -\frac{1}{\rho} \nabla (p^{n+1/2} - q) + \frac{\nu}{2} \left( \nabla^2 \u^{n+1} - \nabla^2 \u^* \right)
\end{equation*}
또는
\begin{equation}
\frac{\u^{n+1} - \u^*}{\Delta t} - \frac{\nu}{2} \nabla^2 \left( \u^{n+1} - \u^* \right) = -\frac{1}{\rho} \nabla (p^{n+1/2} - q)
\label{eq:velocity-correction}
\end{equation}
이로부터 $\u^{n+1} - \u^*$는 어떤 스칼라장의 기울기로 표현됨을 유추할 수 있습니다.
\begin{equation}
\u^{n+1} - \u^* = -\frac{\Delta t}{\rho} \nabla p'
\label{eq:pressure-correction}
\end{equation}
앞의 계수 $-\Delta t/\rho$는 식을 간단하게 하기 위해 적당히 넣은 것입니다. 이 스칼라장 $p’$은 압력 수정(pressure correction)이라고 부릅니다. 식 \eqref{eq:pressure-correction}에 발산(divergence)을 취하면 비압축성 연속 방정식에 의해 $\u^{n+1}$의 발산이 0이므로 압력 수정에 대한 푸아송 방정식을 얻습니다.
\begin{equation}
\nabla^2 p' = \frac{\rho}{\Delta t} \nabla \cdot \u^*
\label{eq:pressure-correction-poisson}
\end{equation}
한편 식 \eqref{eq:pressure-correction}을 식 \eqref{eq:velocity-correction}에 대입하면
\begin{equation*}
\nabla p^{n+1/2} = \nabla q + \nabla p' - \frac{\nu \Delta t}{2} \nabla^2 (\nabla p')
\end{equation*}
일반적으로 라플라시안과 기울기는 교환할 수 있으므로 양변에서 기울기를 없애면 압력에 관한 식을 얻을 수 있습니다.
\begin{equation}
p^{n+1/2} = q + p' - \frac{\nu \Delta t}{2} \nabla^2 p'
\label{eq:pressure-update}
\end{equation}
정리하자면 fractional step method는 다음과 같은 단계로 진행됩니다.
- 식 \eqref{eq:intermediate-velocity}을 풀어 중간 속도 $\u^*$를 구합니다.
- 식 \eqref{eq:pressure-correction-poisson}을 풀어 압력 수정 $p’$을 구합니다.
- 식 \eqref{eq:pressure-correction}를 이용해 속도 $\u^{n+1}$을 구합니다.
- 식 \eqref{eq:pressure-update}을 이용해 압력 $p^{n+1/2}$을 구합니다.
남은 건 제일 첫 단계에 필요한 $q$를 결정하는 것 뿐입니다.
P1 Method
가장 간단한 방법은 $q$를 0으로 두는 것입니다. 이때 식 \eqref{eq:intermediate-velocity}과 식 \eqref{eq:pressure-update}은 각각 다음과 같이 됩니다.
\begin{gather*}
\frac{\u^* - \u^n}{\Delta t} + \left( \frac{3}{2} \N^n - \frac{1}{2} \N^{n-1} \right) = \frac{\nu}{2} \left( \nabla^2 \u^* + \nabla^2 \u^n \right) \\
p^{n+1/2} = p' - \frac{\nu \Delta t}{2} \nabla^2 p'
\end{gather*}
P2 Method
다음으로 $q$를 이전 시각의 압력 $p^{n-1/2}$로 두는 방법이 있습니다. 이때 식 \eqref{eq:intermediate-velocity}과 식 \eqref{eq:pressure-update}은 각각 다음과 같이 됩니다.
\begin{gather}
\frac{\u^* - \u^n}{\Delta t} + \left( \frac{3}{2} \N^n - \frac{1}{2} \N^{n-1} \right) = -\frac{1}{\rho} \nabla p^{n-1/2} + \frac{\nu}{2} \left( \nabla^2 \u^* + \nabla^2 \u^n \right) \nonumber \\
p^{n+1/2} = p^{n-1/2} + p' - \frac{\nu \Delta t}{2} \nabla^2 p'
\label{eq:p2-pressure-update-full-pressure}
\end{gather}
그런데 사실
\begin{equation*}
p^{n+1/2} = p^{n-1/2} + O(\Delta t)
\end{equation*}
이므로 식 \eqref{eq:p2-pressure-update-full-pressure}과 비교하면 $p’$ 자체가 $O(\Delta t)$입니다. 그리고 우변 마지막 항은 $O(\Delta t^2)$이고요. 따라서 압력이 시간에 대한 2차 정확도를 가지는 데에는 사실 이 마지막 항이 필요하지 않습니다. 즉, 다음만으로도 충분합니다.
\begin{equation}
p^{n+1/2} = p^{n-1/2} + p'
\label{eq:p2-pressure-update-not-full-pressure}
\end{equation}
물론 식 \eqref{eq:p2-pressure-update-full-pressure}은 식 \eqref{eq:p2-pressure-update-not-full-pressure}보다 더 높은 정확도를 보여줍니다. 식 \eqref{eq:p2-pressure-update-full-pressure}을 full pressure scheme이라고 부릅니다.
경계 조건