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)$이고요. 따라서 이 마지막 항이 없어도 식 \eqref{eq:p2-pressure-update-full-pressure}은 여전히 시간에 대한 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이라고 부릅니다.
Fractional Step Method의 행렬 표현
식 \eqref{eq:discretized-momentum}의 양변에 $\frac{1}{\rho}\nabla q$를 더하고 $\u^{n+1}$과 $p^{n+1/2}-q$에 대한 방정식임이 잘 드러나도록 정리해봅시다.
\begin{equation*}
\left( \frac{1}{\Delta t} - \frac{\nu}{2} \nabla^2 \right) \u^{n+1} + \frac{1}{\rho} \nabla \left( p^{n+1/2} - q \right) = \left( \frac{1}{\Delta t} + \frac{\nu}{2} \nabla^2 \right) \u^n - \left( \frac{3}{2} \N^n - \frac{1}{2} \N^{n-1} \right) - \frac{1}{\rho} \nabla q
\end{equation*}
이 식과 연속방정식을 합쳐 $\u^{n+1}$과 $p^{n+1/2} - q$에 대한 행렬 방정식을 쓸 수 있습니다.
\begin{equation}
\begin{bmatrix}
\A & \frac{1}{\rho} \G \\
\D & 0
\end{bmatrix}
\begin{bmatrix}
\u^{n+1} \\
p^{n+1/2} - q
\end{bmatrix}
=
\begin{bmatrix}
\mathbf{r} \\
0
\end{bmatrix}
\label{eq:matrix-form-fractional-step}
\end{equation}
여기서
\begin{align*}
\A &= \frac{1}{\Delta t} \I - \frac{\nu}{2} \L \\
\mathbf{r} &= \left( \frac{1}{\Delta t} \I + \frac{\nu}{2} \L \right) \u^n - \left( \frac{3}{2} \N^n - \frac{1}{2} \N^{n-1} \right) - \frac{1}{\rho} \G q \\
\end{align*}
이고 $\I$, $\G$, $\D$, $\L$은 각각 항등, 기울기, 발산, 라플라시안 연산자입니다.
핵심은 식 \eqref{eq:matrix-form-fractional-step}의 행렬을 다음과 같이 LU 분해할 수 있다는 것입니다. [4]
\begin{equation*}
\begin{bmatrix}
\A & \frac{1}{\rho} \G \\
\D & 0
\end{bmatrix}
=
\begin{bmatrix}
\A & 0 \\
\D & -\frac{1}{\rho} \D \A^{-1} \G
\end{bmatrix}
\begin{bmatrix}
\I & \frac{1}{\rho} \A^{-1} \G \\
0 & \I
\end{bmatrix}
\end{equation*}
따라서 식 \eqref{eq:matrix-form-fractional-step}은 LU 분해를 통해 두 단계로 나누어 풀 수 있습니다.
\begin{align*}
\begin{bmatrix}
\A & 0 \\
\D & -\frac{1}{\rho} \D \A^{-1} \G
\end{bmatrix}
\begin{bmatrix}
\u^* \\
p^*
\end{bmatrix}
&=
\begin{bmatrix}
\mathbf{r} \\
0
\end{bmatrix} \\
\begin{bmatrix}
\I & \frac{1}{\rho} \A^{-1} \G \\
0 & \I
\end{bmatrix}
\begin{bmatrix}
\u^{n+1} \\
p^{n+1/2} - q
\end{bmatrix}
&=
\begin{bmatrix}
\u^* \\
p^*
\end{bmatrix}
\end{align*}
전부 정리하면 다음과 같습니다.
\begin{gather}
\A \u^* = \mathbf{r} \label{eq:matrix-form-intermediate-velocity} \\
\frac{1}{\rho} \D \A^{-1} \G \left( p^{n+1/2} - q \right) = \D \u^* \label{eq:matrix-form-pressure-difference} \\
\u^{n+1} = \u^* - \frac{1}{\rho} \A^{-1} \G \left( p^{n+1/2} - q \right) \label{eq:matrix-form-velocity-update}
\end{gather}
이제 슬슬 LU 분해가 위에서 fractional step method에서 중간 속도를 도입하는 것과 동일하다는 걸 느끼실 겁니다. 식 \eqref{eq:matrix-form-intermediate-velocity}은 중간 속도를 구하는 식 \eqref{eq:intermediate-velocity}과 완전히 같습니다. 그리고 식 \eqref{eq:matrix-form-pressure-difference}로 압력 $p^{n+1/2}$을 구하고, 식 \eqref{eq:matrix-form-velocity-update}으로 속도 $\u^{n+1}$을 구하는 것까지 비슷합니다.
위에서 언급했듯이 일반적으로 $\L$과 $\G$는 서로 교환할 수 있으므로 $\A$와 $\G$는 교환 가능합니다. 따라서 $\A^{-1}$도 $\G$와 교환할 수 있고, 식 \eqref{eq:matrix-form-pressure-difference}와 식 \eqref{eq:matrix-form-velocity-update}에 적용하면
\begin{gather*}
\frac{1}{\rho} \D \G \A^{-1} \left( p^{n+1/2} - q \right) = \D \u^* \\
\u^{n+1} = \u^* - \frac{1}{\rho} \G \A^{-1} \left( p^{n+1/2} - q \right)
\end{gather*}
또한 $(1 / \Delta t) \A^{-1} (p^{n+1/2} - q)$를 압력 수정 $p’$으로 정의하면
\begin{gather*}
\frac{\Delta t}{\rho} \D \G p' = \D \u^* \\
\u^{n+1} = \u^* - \frac{\Delta t}{\rho} \G p' \\
p^{n+1/2} = q + \Delta t \A p'
\end{gather*}
가 되어 각각 식 \eqref{eq:pressure-correction-poisson}, 식 \eqref{eq:pressure-correction}, 식 \eqref{eq:pressure-update}과 완벽히 동일합니다.
라플라시안과 기울기의 교환
앞서 압력 수정과 관련한 부분에서 “라플라시안 $\L$과 기울기 $\G$는 교환 가능하다”라는 성질을 사용했는데, 사실 공간에 대한 이산화를 적용하면 이 성질이 성립하지 않을 수 있습니다. 즉, $\A^{-1}$와 $\G$가 교환 가능하지 않아 다음과 같이 오차가 발생할 수 있습니다.
\begin{equation}
\A^{-1} \G = \G \A^{-1} + \mathbf{Q}
\label{eq:commutative-error}
\end{equation}
이 경우 식 \eqref{eq:matrix-form-pressure-difference}와 식 \eqref{eq:matrix-form-velocity-update}에서 $\A^{-1}$와 $\G$를 교환하면 오차항이 붙게 됩니다.
\begin{gather*}
\frac{\Delta t}{\rho} \D \G p' = \D \u^* - \frac{1}{\rho} \D \mathbf{Q} \left( p^{n+1/2} - q \right) \\
\u^{n+1} = \u^* - \frac{\Delta t}{\rho} \G p' - \frac{1}{\rho} \mathbf{Q} \left( p^{n+1/2} - q \right)
\end{gather*}
여기서 시간에 대한 정확도를 분석해봅시다. $\A = O(1/\Delta t)$이므로 $\A^{-1} = O(\Delta t)$이고, 식 \eqref{eq:commutative-error}에 의해 $\mathbf{Q}$ 역시 $O(\Delta t)$입니다. 그러니 추가된 오차항은 시간에 대해 1차 정확도를 갖고, 지금까지 방정식을 2차 정확도로 이산화해왔던 것과 비교하면 정확도에 악영향을 줍니다. 다행히, 대부분의 유한 차분법을 이용한 공간 이산화는 라플라시안과 기울기가 교환 가능하지만, 경계 부근에서는 경계 조건과 결합되면서 이 성질이 깨질 수 있으므로 공간 이산화 시 주의해야 합니다. [5]
참고 문헌
- J. Kim and P. Moin, Application of a fractional-step method to incompressible Navier-Stokes equations, J. Comput. Phys., 59, 308-323 (1985).
- S. Armfield and R. Street, An analysis and comparison of the time accuracy of fractional-step methods for the Navier-Stokes equations on staggered grids, Int. J. Numer. Meth. Fluids, 38, 255-282 (2002).
- S. Armfield and R. Street, The pressure accuracy of fractional-step methods for the Navier-Stokes equations on Staggered Grids, ANZIAM J., 44, C20-C39 (2003).
- J. Perot, An analysis of the fractional step method, J. Comput. Phys., 108, 51-58 (1993).
- J. Perot, Comments on the fractional step method, J. Comput. Phys., 121, 190-191 (1995).
- D. Brown, Accuracy of projection methods for the incompressible Navier-Stokes equations, Numerical Simulations of Incompressible Flows, 101-107 (2003).