Fractional Step Method

2025년 6월 3일

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는 다음과 같은 단계로 진행됩니다.

  1. 식 \eqref{eq:intermediate-velocity}을 풀어 중간 속도 $\u^*$를 구합니다.
  2. 식 \eqref{eq:pressure-correction-poisson}을 풀어 압력 수정 $p’$을 구합니다.
  3. 식 \eqref{eq:pressure-correction}를 이용해 속도 $\u^{n+1}$을 구합니다.
  4. 식 \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이라고 부릅니다.

경계 조건

CFD