Fractional Step Method

개요

Fractional step method(FSM)는 비압축성 비정상(unsteady) 유동을 효과적으로 풀어내는 기법으로, 중간 변수를 도입하여 반복 없이 다음 time step의 속도와 압력을 구할 수 있다는 특징이 있습니다. 이 글에서는 Kim & Moin[1]의 설명을 바탕으로 Perot[2]가 제시한 연산자 분해의 관점으로 FSM을 소개합니다.

지배방정식의 이산화

비압축성 운동량 방정식

\begin{equation*} \frac{\partial \u}{\partial t} + \u \cdot \nabla \u = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \u \end{equation*}

을 time step $n+1/2$에서 이산화해봅시다. 이때 시간에 대한 미분 $\partial\u/\partial t$은 중앙차분법으로, 대류항 $\N = \u \cdot \nabla \u$는 Adam-Bashforth 방법으로, 점성항 $\nu \nabla^2 \u$는는 Crank-Nicolson 방법으로 시간에 대해 이산화해봅시다. [1] 즉, 아래처럼 시간에 대해 2차 정확도를 가지도록 이산화하겠다는 뜻입니다.

\begin{align*} \frac{\partial \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 \u^{n+1/2} &= \frac{1}{2} \left( \nabla^2 \u^{n+1} + \nabla^2 \u^n \right) + O(\Delta t^2) \end{align*}

결과는 다음과 같습니다.

\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}

양변에 $(1/\rho)\nabla q$를 더하고 적당히 정리해줍니다. 여기서 $q$는 값이 알려져 있는 적당한 스칼라장입니다.

\begin{equation*} \left( 1 - \frac{\nu \Delta t}{2} \nabla^2 \right) \u^{n+1} + \frac{\Delta t}{\rho} \nabla \left( p^{n+1/2} - q \right) = \left( 1 + \frac{\nu \Delta t}{2} \nabla^2 \right) \u^n - \Delta t \left( \frac{3}{2} \N^n - \frac{1}{2} \N^{n-1} \right) - \frac{\Delta t}{\rho} \nabla q \equiv \mathbf{r} \end{equation*}

이제 공간 이산화를 적용해서 위 식을 선형계로 바꾸면 다음과 같은 형태가 됩니다. 편의상 2차원을 가정하고, x 방향과 y 방향을 분리해서 썼습니다.

\begin{align} \left( \I - \frac{\nu \Delta t}{2} \L \right) u^{n+1} + \frac{\Delta t}{\rho} \G_x^\opp \left( p^{n+1/2} - q \right) &= r_x + b_{\mathrm{mom},x} \label{eq:discretized-u-momentum} \\ \left( \I - \frac{\nu \Delta t}{2} \L \right) v^{n+1} + \frac{\Delta t}{\rho} \G_y^\opp \left( p^{n+1/2} - q \right) &= r_y + b_{\mathrm{mom},y} \label{eq:discretized-v-momentum} \end{align}

여기서 $\I$는 항등 연산자, $\L$은 속도에 대한 이산화된 라플라시안 연산자, $\G_x^\opp$와 $\G_y^\opp$는 압력에 대한 각각 x 방향과 y 방향의 이산화된 기울기 연산자이고, $b_{\mathrm{mom},x}$와 $b_{\mathrm{mom},y}$는 $\u^{n+1}$과 $p^{n+1/2} - q$의 경계 조건 때문에 생기는 상수항입니다.

한편 time step $n+1$에서의 연속 방정식

\begin{equation*} \nabla \cdot \u^{n+1} = 0 \end{equation*}

도 똑같이 공간에 대해 이산화하면

\begin{equation} \G_x^\opv u^{n+1} + \G_y^\opv v^{n+1} = b_\mathrm{cont} \label{eq:discretized-continuity} \end{equation}

마찬가지로 $\G_x^\opv$와 $\G_y^\opv$는 속도에 대한 이산화된 기울기 연산자이고, $b_\mathrm{cont}$는 경계 조건 때문에 생기는 상수항입니다.

식 \eqref{eq:discretized-u-momentum}, \eqref{eq:discretized-v-momentum}, \eqref{eq:discretized-continuity}를 합쳐 행렬로 표현할 수 있습니다.

\begin{equation} \begin{bmatrix} \A & 0 & \frac{\Delta t}{\rho} \G_x^\opp \\ 0 & \A & \frac{\Delta t}{\rho} \G_y^\opp \\ \G_x^\opv & \G_y^\opv & 0 \end{bmatrix} \begin{bmatrix} u^{n+1} \\ v^{n+1} \\ p^{n+1/2} - q \end{bmatrix} = \begin{bmatrix} r_x + b_{\mathrm{mom},x} \\ r_y + b_{\mathrm{mom},y} \\ b_\mathrm{cont} \end{bmatrix} \label{eq:matrix-representation} \end{equation}

여기서 $\A = \I - (\nu \Delta t / 2) \L$입니다. 식 \eqref{eq:matrix-representation}를 직접 풀면 전형적인 coupled solver가 되겠지만, 이건 너무 비효율적이니 작은 선형계로 분리하여 풀어보도록 합시다.

연산자 분해로서의 Fractional Step Method

핵심은 식 \eqref{eq:matrix-representation}의 좌변에 등장하는 연산자 행렬을 다음과 같이 LU 분해할 수 있다는 점입니다.

\begin{equation} \begin{bmatrix} \A & 0 & \frac{\Delta t}{\rho} \G_x^\opp \\ 0 & \A & \frac{\Delta t}{\rho} \G_y^\opp \\ \G_x^\opv & \G_y^\opv & 0 \end{bmatrix} = \begin{bmatrix} \A & 0 & 0 \\ 0 & \A & 0 \\ \G_x^\opv & \G_y^\opv & -\frac{\Delta t}{\rho} ( \G_x^\opv \A^{-1} \G_x^\opp + \G_y^\opv \A^{-1} \G_y^\opp ) \end{bmatrix} \begin{bmatrix} \I & 0 & \frac{\Delta t}{\rho} \A^{-1} \G_x^\opp \\ 0 & \I & \frac{\Delta t}{\rho} \A^{-1} \G_y^\opp \\ 0 & 0 & \I \end{bmatrix} \label{eq:operator-decomposition} \end{equation}

이를 활용하면 중간 속도(intermediate velocity) $\u^*$를 도입하여 식 \eqref{eq:matrix-representation}를 두 단계로 나누어 풀 수 있습니다.

1단계

\begin{gather} \A u^* = r_x + b_{\mathrm{mom},x} \nonumber \\ \A v^* = r_y + b_{\mathrm{mom},y} \nonumber \\ \frac{\Delta t}{\rho} \left( \G_x^\opv \A^{-1} \G_x^\opp + \G_y^\opv \A^{-1} \G_y^\opp \right) \left( p^{n+1/2} - q \right) = \G_x^\opv u^* + \G_y^\opv v^* + b_\mathrm{cont} \label{eq:discretized-pressure-equation} \end{gather}

2단계

\begin{gather} u^{n+1} = u^* - \frac{\Delta t}{\rho} \A^{-1} \G_x^\opp \left( p^{n+1/2} - q \right) \label{eq:x-velocity-update} \\ v^{n+1} = v^* - \frac{\Delta t}{\rho} \A^{-1} \G_y^\opp \left( p^{n+1/2} - q \right) \label{eq:y-velocity-update} \\ \end{gather}

한 가지 문제가 있다면 $\A^{-1}$를 구하는 것이 여전히 매우 비효율적이라는 것인데, 이 문제는 두 가지 방법으로 해결 가능합니다. 첫 번째는 $\A$의 정의로부터 $\A^{-1}$를 시간에 대한 2차 정확도로 근사하는 것입니다.

\begin{equation*} \A^{-1} = \left( \I - \frac{\nu \Delta t}{2} \L \right)^{-1} = \I + \frac{\nu \Delta t}{2} \L + O(\Delta t^2) \end{equation*}

이걸 식 \eqref{eq:discretized-pressure-equation}에 대입하면 좌변의 행렬이 꽤 큰 스텐실을 필요로 합니다. 기울기 두 개와 라플라시안을 곱하니 공간에 대한 4차 미분이고, 공간에 대한 2차 정확도를 위해선 적어도 5개나 6개 점이 필요하죠. 이게 부담이라면 $\A^{-1}$를 그냥 $\I$로 근사하고 시간에 대한 정확도는 1차로 낮추는 방법도 있습니다. [2]

두 번째 방법은 $\A$와 $\G_x^\opp$, $\G_y^\opp$가 교환 가능하다고 가정하는 것입니다. 벡터 미적분학의 항등식에 따라 원래 라플라시안과 기울기는 교환할 수 있는데, 만약 이산화된 라플라시안 연산자와 압력에 대한 기울기 연산자도 그렇다면 $\A$와 $\G_x^\opp$, $\G_y^\opp$ 역시 교환할 수 있습니다. 그렇다면 $\A^{-1}$도 $\G_x^\opp$, $\G_y^\opp$와 교환할 수 있으므로, 다음과 같이 식 \eqref{eq:discretized-pressure-equation} – \eqref{eq:y-velocity-update}를 다시 쓸 수 있습니다.

\begin{gather*} \frac{\Delta t}{\rho} \left( \G_x^\opv \G_x^\opp + \G_y^\opv \G_y^\opp \right) \A^{-1} \left( p^{n+1/2} - q \right) = \G_x^\opv u^* + \G_y^\opv v^* + b_\mathrm{cont} \\ u^{n+1} = u^* - \frac{\Delta t}{\rho} \G_x^\opp \A^{-1} \left( p^{n+1/2} - q \right) \\ v^{n+1} = v^* - \frac{\Delta t}{\rho} \G_y^\opp \A^{-1} \left( p^{n+1/2} - q \right) \end{gather*}

여기서 압력 수정(pressure correction) $p’$을 $\A^{-1} (p^{n+1/2} - q)$로 정의하면 식이 훨씬 간단해집니다.

\begin{gather} \frac{\Delta t}{\rho} \left( \G_x^\opv \G_x^\opp + \G_y^\opv \G_y^\opp \right) p' = \G_x^\opv u^* + \G_y^\opv v^* + b_\mathrm{cont} \label{eq:pressure-correction} \\ u^{n+1} = u^* - \frac{\Delta t}{\rho} \G_x^\opp p' \nonumber \\ v^{n+1} = v^* - \frac{\Delta t}{\rho} \G_y^\opp p' \nonumber \\ p^{n+1/2} - q = \A p' \nonumber \end{gather}

두 번째 방법의 문제는 이산화된 라플라시안과 기울기 연산자가 교환 가능하다는 보장이 없다는 것입니다. 일반적으로 생각할 수 있는 중앙 차분법을 사용하는 공간 이산화 기법을 사용한다면 대개 교환 가능하지만, 경계 부근에서는 경계 조건과 결합되면서 이 성질이 깨질 수 있습니다. 특히나 이산화된 라플라시안은 속도에 대한 연산자이지만 이산화된 기울기는 압력에 대한 연산자이기 때문에 서로 다른 경계 조건이 적용되어 더욱 그럴 가능성이 큽니다. 이 경우에는 오차에 의하여 시간에 대한 정확도가 1차로 떨어지게 됩니다. [3]

이하에서는 압력 수정을 도입하는 두 번째 방법을 가정합니다.

비엇갈림 격자 상에서의 Fractional Step Method

식 \eqref{eq:pressure-correction}의 좌변에 등장하는 연산자는 기울기 연산자를 두 번 포함하는데, 이는 속도와 압력이 같은 위치에서 계산되는 비엇갈림 격자(non-staggered grid) 상에서는 문제가 될 수 있습니다. 균일한 2차원 직교 격자를 가정하고, 연산자가 어떤 스텐실을 갖는지 계산해봅시다. (간단히 중앙 차분법을 가정하고, 경계 조건은 고려하지 않습니다.)

\begin{align*} \left. \left( \G_x^\opv \G_x^\opp + \G_y^\opv \G_y^\opp \right) \phi \right|_{i,j} &= \frac{\G_x^\opp \phi |_{i+1,j} - \G_x^\opp \phi |_{i-1,j}}{2 \Delta x} + \frac{\G_y^\opp \phi |_{i,j+1} - \G_y^\opp \phi |_{i,j-1}}{2 \Delta y} \\ &= \frac{\left( \phi_{i+2,j} - \phi_{i,j} \right) / (2 \Delta x) - \left( \phi_{i,j} - \phi_{i-2,j} \right) / (2 \Delta x)}{2 \Delta x} \\ &\qquad+ \frac{\left( \phi_{i,j+2} - \phi_{i,j} \right) / (2 \Delta y) - \left( \phi_{i,j} - \phi_{i,j-2} \right) / (2 \Delta y)}{2 \Delta y} \\ &= \frac{\phi_{i+2,j} - 2 \phi_{i,j} + \phi_{i-2,j}}{4 \Delta x^2} + \frac{\phi_{i,j+2} - 2 \phi_{i,j} + \phi_{i,j-2}}{4 \Delta y^2} \end{align*}

보통의 조밀한 5점 스텐실(5-point stencil)을 사용한 라플라시안과 달리 한 점씩 더 멀리 떨어진, 상대적으로 성긴 스텐실의 라플라시안이 된 것을 볼 수 있습니다. 이는 비효율적일 뿐만 아니라 비엇갈림 격자에서 흔히 문제시되는 체커보드 압력 현상을 유발하는 원인이 됩니다. 스텐실 특성상 짝수 번째 격자점끼리, 홀수 번째 격자점끼리만 연관되고, 바로 인접한 격자점끼리는 아무런 상관 관계가 없어 해가 따로 놀기 때문이죠. 그림 1의 예시는 왜 체커보드라고 부르는지 그 이유를 잘 보여줍니다.

한편 비엇갈림 격자에 대비되는 엇갈림 격자(staggered grid)를 살펴봅시다. 이 경우 압력의 기울기는 속도 계산에 필요하므로 셀 면에서, 속도의 기울기는 셀 중심에서 계산합니다. 즉,

\begin{align*} \left. \G_x^{\opf\opp} p \right|_{i-1/2,j} &= \frac{p_{i,j} - p_{i-1,j}}{\Delta x} \\ \left. \G_y^{\opf\opp} p \right|_{i,j-1/2} &= \frac{p_{i,j} - p_{i,j-1}}{\Delta y} \\ \left. \G_x^{\opf\opv} U \right|_{i,j} &= \frac{U_{i+1/2,j} - U_{i-1/2,j}}{\Delta x} \\ \left. \G_x^{\opf\opv} V \right|_{i,j} &= \frac{V_{i,j+1/2} - V_{i,j-1/2}}{\Delta y} \end{align*}

여기서 $U$와 $V$는 엇갈림 속도로, 셀 면에서 계산되며 구분을 위해 대문자로 표기했습니다. 연산자들 또한 비엇갈림 격자와 구분할 수 있게 윗첨자 $\opf$를 붙였습니다. 엇갈림 속도의 자세한 위치 관계는 그림 2를 참고해주세요.

따라서 엇갈림 격자에서 위 연산자는

\begin{equation} \begin{aligned} \left. \left( \G_x^{\opf\opv} \G_x^{\opf\opp} + \G_y^{\opf\opv} \G_y^{\opf\opp} \right) \phi \right|_{i,j} &= \frac{\G_x^{\opf\opp} \phi |_{i+1/2,j} - \G_x^{\opf\opp} \phi |_{i-1/2,j}}{\Delta x} + \frac{\G_y^{\opf\opp} \phi |_{i,j+1/2} - \G_y^{\opf\opp} \phi |_{i,j-1/2}}{\Delta y} \\ &= \frac{\left( \phi_{i+1,j} - \phi_{i,j} \right) / \Delta x - \left( \phi_{i,j} - \phi_{i-1,j} \right) / \Delta x}{\Delta x} \\ &\qquad+ \frac{\left( \phi_{i,j+1} - \phi_{i,j} \right) / \Delta y - \left( \phi_{i,j} - \phi_{i,j-1} \right) / \Delta y}{\Delta y} \\ &= \frac{\phi_{i+1,j} - 2 \phi_{i,j} + \phi_{i-1,j}}{\Delta x^2} + \frac{\phi_{i,j+1} - 2 \phi_{i,j} + \phi_{i,j-1}}{\Delta y^2} \end{aligned} \label{eq:staggered-pressure-laplacian} \end{equation}

가 되어 조밀한 5점 스텐실을 사용하게 되고, 체커보드 압력 문제가 발생하지 않습니다.

그렇다면 비엇갈림 격자에서 이 장점을 흡수할 수는 없을까요? 한 가지 아이디어는 연속 방정식을 엇갈림 속도에 대한 식으로 쓰고, 엇갈림 속도 자체는 셀 중심 속도의 보간(interpolation)으로 정의하는 것입니다.

\begin{gather} \G_x^{\opf\opv} U^{n+1} + \G_y^{\opf\opv} V^{n+1} = b_\mathrm{cont} \label{eq:face-velocity-continuity} \\ U^{n+1} = \T_x u^{n+1} + b_{\mathrm{interp},x} \label{eq:u-interpolation} \\ V^{n+1} = \T_y v^{n+1} + b_{\mathrm{interp},y} \label{eq:v-interpolation} \\ \end{gather}

$\T_x$와 $\T_y$는 각각 $u$와 $v$에 대한 보간 연산자이고 $b_{\mathrm{interp},x}$와 $b_{\mathrm{interp},y}$는 속도의 경계 조건에 의한 상수항입니다. 위와 마찬가지로 식 \eqref{eq:face-velocity-continuity}, \eqref{eq:u-interpolation}, \eqref{eq:v-interpolation}를 식 \eqref{eq:discretized-u-momentum}, \eqref{eq:discretized-v-momentum}과 합쳐 행렬로 쓸 수 있습니다.

\begin{equation} \begin{bmatrix} \A & 0 & 0 & 0 & \frac{\Delta t}{\rho} \G_x^\opp \\ 0 & \A & 0 & 0 & \frac{\Delta t}{\rho} \G_y^\opp \\ \T_x & 0 & -\I & 0 & 0 \\ 0 & \T_y & 0 & -\I & 0 \\ 0 & 0 & \G_x^{\opf\opv} & \G_y^{\opf\opv} & 0 \end{bmatrix} \begin{bmatrix} u^{n+1} \\ v^{n+1} \\ U^{n+1} \\ V^{n+1} \\ p^{n+1/2} - q \end{bmatrix} = \begin{bmatrix} r_x + b_{\mathrm{mom},x} \\ r_y + b_{\mathrm{mom},y} \\ b_{\mathrm{interp},x} \\ b_{\mathrm{interp},y} \\ b_\mathrm{cont} \end{bmatrix} \label{eq:matrix-representation-interpolation} \end{equation}

좌변의 행렬은 아래와 같이 LU 분해됩니다.

\begin{equation*} \begin{bmatrix} \A & 0 & 0 & 0 & 0 \\ 0 & \A & 0 & 0 & 0 \\ \T_x & 0 & -\I & 0 & 0 \\ 0 & \T_y & 0 & -\I & 0 \\ 0 & 0 & \G_x^{\opf\opv} & \G_y^{\opf\opv} & -\frac{\Delta t}{\rho}\left( \G_x^{\opf\opv}\T_x\A^{-1}\G_x^\opp + \G_y^{\opf\opv}\T_y\A^{-1}\G_y^\opp \right) \end{bmatrix} \begin{bmatrix} \I & 0 & 0 & 0 & \frac{\Delta t}{\rho}\A^{-1}\G_x^\opp \\ 0 & \I & 0 & 0 & \frac{\Delta t}{\rho}\A^{-1}\G_y^\opp \\ 0 & 0 & \I & 0 & \frac{\Delta t}{\rho}\T_x\A^{-1}\G_x^\opp \\ 0 & 0 & 0 & \I & \frac{\Delta t}{\rho}\T_y\A^{-1}\G_y^\opp \\ 0 & 0 & 0 & 0 & \I \end{bmatrix} \end{equation*}

똑같이 기울기 연산자와 라플라시안 연산자가 교환된다고 가정하고 압력 수정을 도입하면 엇갈림 중간 속도를 포함하는 두 단계의 방정식을 얻습니다.

1단계

\begin{gather} \A u^* = r_x + b_{\mathrm{mom},x} \nonumber \\ \A v^* = r_y + b_{\mathrm{mom},y} \nonumber \\ U^* = \T_x u^* + b_{\mathrm{interp},x} \nonumber \\ V^* = \T_y v^* + b_{\mathrm{interp},y} \nonumber \\ \frac{\Delta t}{\rho}\left( \G_x^{\opf\opv}\T_x\G_x^\opp + \G_y^{\opf\opv}\T_y\G_y^\opp \right) p' = \G_x^{\opf\opv}U^* + \G_y^{\opf\opv}V^* + b_\mathrm{cont} \label{eq:pressure-correction-interp} \end{gather}

2단계

\begin{gather*} u^{n+1} = u^* - \frac{\Delta t}{\rho} \G_x^\opp \\ v^{n+1} = v^* - \frac{\Delta t}{\rho} \G_y^\opp \\ U^{n+1} = U^* - \frac{\Delta t}{\rho} \T_x \G_x^\opp \\ V^{n+1} = V^* - \frac{\Delta t}{\rho} \T_y \G_y^\opp \end{gather*}

그런데 식 \eqref{eq:pressure-correction-interp}의 연산자는 식 \eqref{eq:staggered-pressure-laplacian}과 다르죠. 선형 보간을 가정하고 스텐실을 계산해보면 원래의 성긴 라플라시안 연산자와 동일합니다.

\begin{align*} \left. \left( \G_x^{\opf\opv}\T_x\G_x^\opp + \G_y^{\opf\opv}\T_y\G_y^\opp \right) \phi \right|_{i,j} &= \frac{\T_x\G_x^\opp \phi |_{i+1/2,j} - \T_x\G_x^\opp \phi |_{i-1/2,j}}{\Delta x} + \frac{\T_y\G_y^\opp \phi |_{i,j+1/2} - \T_y\G_y^\opp \phi |_{i,j-1/2}}{\Delta y} \\ &= \frac{(\G_x^\opp \phi |_{i+1,j} + \G_x^\opp \phi |_{i,j}) / 2 - (\G_x^\opp \phi |_{i,j} + \G_x^\opp \phi |_{i-1,j}) / 2}{\Delta x} \\ &\qquad+ \frac{(\G_y^\opp \phi |_{i,j+1} + \G_y^\opp \phi |_{i,j}) / 2 - (\G_y^\opp \phi |_{i,j} + \G_y^\opp \phi |_{i,j-1}) / 2}{\Delta y} \\ &= \frac{\G_x^\opp \phi |_{i+1,j} - \G_x^\opp \phi |_{i-1,j}}{2\Delta x} + \frac{\G_y^\opp \phi |_{i,j+1} - \G_y^\opp \phi |_{i,j-1}}{2\Delta y} \\ &= \frac{\phi_{i+2,j} - 2\phi_{i,j} + \phi_{i-2,j}}{4\Delta x^2} + \frac{\phi_{i,j+2} - 2\phi_{i,j} + \phi_{i,j-2}}{4\Delta y^2} \end{align*}

사실, 엇갈림 속도를 선형 보간으로 계산하게 되면 본질적으로 달라지는 건 하나도 없습니다. 당장 위와 같이 압력 수정의 라플라시안도 같고, 엇갈림 속도의 발산 역시 원래대로 속도의 발산을 계산한 것과 다르지 않거든요. (물론 upwind나 QUICK scheme을 사용하면 얘기가 다르지만)

식 \eqref{eq:staggered-pressure-laplacian}의 연산자를 그대로 쓰려면 행렬 분해가 다음과 같은 꼴이 되어야 합니다. (연산자가 잘 드러나도록 일부러 $\A^{-1}$를 마지막에 썼습니다.)

\begin{equation*} \begin{bmatrix} \A & 0 & 0 & 0 & 0 \\ 0 & \A & 0 & 0 & 0 \\ \T_x & 0 & -\I & 0 & -\frac{\Delta t}{\rho}\left( \T_x\G_x^\opp - \G_x^{\opf\opp} \right)\A^{-1} \\ 0 & \T_y & 0 & -\I & -\frac{\Delta t}{\rho}\left( \T_y\G_y^\opp - \G_y^{\opf\opp} \right)\A^{-1} \\ 0 & 0 & \G_x^{\opf\opv} & \G_y^{\opf\opv} & -\frac{\Delta t}{\rho}\left( \G_x^{\opf\opv}\G_x^{\opf\opp} + \G_y^{\opf\opv}\G_y^{\opf\opp} \right)\A^{-1} \end{bmatrix} \begin{bmatrix} \I & 0 & 0 & 0 & \frac{\Delta t}{\rho}\G_x^\opp\A^{-1} \\ 0 & \I & 0 & 0 & \frac{\Delta t}{\rho}\G_y^\opp\A^{-1} \\ 0 & 0 & \I & 0 & \frac{\Delta t}{\rho}\G_x^{\opf\opp}\A^{-1} \\ 0 & 0 & 0 & \I & \frac{\Delta t}{\rho}\G_y^{\opf\opp}\A^{-1} \\ 0 & 0 & 0 & 0 & \I \end{bmatrix} \end{equation*}

비교해보면 엇갈림 중간 속도의 보간법이 달라졌습니다.

\begin{align*} U^* &= \T_x u^* - \frac{\Delta t}{\rho} \left( \T_x \G_x^\opp - \G_x^{\opf\opp} \right) p' + b_{\mathrm{interp},x} \\ V^* &= \T_y v^* - \frac{\Delta t}{\rho} \left( \T_y \G_y^\opp - \G_y^{\opf\opp} \right) p' + b_{\mathrm{interp},y} \end{align*}

이 식은 비엇갈림 격자에서 체커보드 압력 현상을 없애기 위해 사용되는 Rhie-Chow 보간과 유사한 형태를 하고 있습니다. 큰 차이점이라면 압력 대신 압력 수정을 사용한다는 점 정도네요.

한 가지 문제점은, 압력 수정 방정식의 우변은 엇갈림 중간 속도의 발산인데, 정작 엇갈림 중간 속도를 보간하는 데에 압력 수정이 필요하다는 것입니다. 이는 Rhie-Chow 보간을 사용한 SIMPLE 알고리즘처럼 반복법으로 압력 수정을 계산하는 방법으로 해결할 수 있지만 반복 없이 다음 단계를 계산할 수 있다는 FSM의 장점을 포기히게 됩니다. 또 다른 방법으로는 엇갈림 중간 속도 보간에서 압력 수정 항을 없애는 것이 있는데, 이렇게 하면 방정식이 살짝 바뀌게 됩니다.

\begin{align*} &\begin{bmatrix} \A & 0 & 0 & 0 & 0 \\ 0 & \A & 0 & 0 & 0 \\ \T_x & 0 & -\I & 0 & 0 \\ 0 & \T_y & 0 & -\I & 0 \\ 0 & 0 & \G_x^{\opf\opv} & \G_y^{\opf\opv} & -\frac{\Delta t}{\rho}\left( \G_x^{\opf\opv}\G_x^{\opf\opp} + \G_y^{\opf\opv}\G_y^{\opf\opp} \right)\A^{-1} \end{bmatrix} \begin{bmatrix} \I & 0 & 0 & 0 & \frac{\Delta t}{\rho}\G_x^\opp\A^{-1} \\ 0 & \I & 0 & 0 & \frac{\Delta t}{\rho}\G_y^\opp\A^{-1} \\ 0 & 0 & \I & 0 & \frac{\Delta t}{\rho}\G_x^{\opf\opp}\A^{-1} \\ 0 & 0 & 0 & \I & \frac{\Delta t}{\rho}\G_y^{\opf\opp}\A^{-1} \\ 0 & 0 & 0 & 0 & \I \end{bmatrix} \\ &\qquad= \begin{bmatrix} \A & 0 & 0 & 0 & \frac{\Delta t}{\rho}\G_x^\opp \\ 0 & \A & 0 & 0 & \frac{\Delta t}{\rho}\G_y^\opp \\ \T_x & 0 & -\I & 0 & \frac{\Delta t}{\rho}\left( \T_x\G_x^\opp - \G_x^{\opf\opp} \right)\A^{-1} \\ 0 & \T_y & 0 & -\I & \frac{\Delta t}{\rho}\left( \T_y\G_y^\opp - \G_y^{\opf\opp} \right)\A^{-1} \\ 0 & 0 & \G_x^{\opf\opv} & \G_y^{\opf\opv} & 0 \end{bmatrix} \end{align*}

식 \eqref{eq:matrix-representation-interpolation}와 비교하면 이번엔 엇갈림 속도에 Rhie-Chow 보간이 적용된 꼴이 되었습니다.

\begin{align*} U^{n+1} &= \T_x u^{n+1} + \frac{\Delta t}{\rho} \left( \T_x \G_x^\opp - \G_x^{\opf\opp} \right) p' + b_{\mathrm{interp},x} \\ V^{n+1} &= \T_y v^{n+1} + \frac{\Delta t}{\rho} \left( \T_y \G_y^\opp - \G_y^{\opf\opp} \right) p' + b_{\mathrm{interp},y} \end{align*}

이 방법의 단점은 엇갈림 속도를 단순히 속도의 선형 보간으로 계산할 수 없기 때문에 엇갈림 속도의 발산과 속도의 발산이 달라진다는 것입니다. 즉, 연속 방정식은 엇갈림 속도에 대해서만 성립하고, 셀 중심에서 속도(또는 속도를 포함하는 값)의 기울기를 계산할 때에도 엇갈림 속도를 사용하는 것이 더 정확합니다. [4, 5] 대표적으로 대류항은 경계 조건을 무시할 때 다음과 같이 게산할 수 있습니다. 여기서 $\odot$은 두 벡터의 원소별 곱(pointwise product)을 의미합니다.

\begin{align*} \N_x^n &= \G_x^{\opf\opv} \left[ U^n \odot \left( \T_x u^n \right) \right] + \G_y^{\opf\opv} \left[ V^n \odot \left( \T_y u^n \right) \right] \\ \N_y^n &= \G_x^{\opf\opv} \left[ U^n \odot \left( \T_x v^n \right) \right] + \G_y^{\opf\opv} \left[ V^n \odot \left( \T_y v^n \right) \right] \end{align*}

이렇게 (속도의 단순 보간으로 정의되지 않는) 엇갈림 속도를 별도로 계산하는 것은 Zang et al.[6]에서 소개된 바 있습니다.

참고 문헌

J. Kim and P. Moin, Application of a fractional-step method to incompressible Navier-Stokes equations, J. Comput. Phys., 59, 308-323 (1985).
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).
R. Mittal, H. Dong, M. Bozkurttas, F.M. Najjar, A. Vargas, and A. von Loebbecke, A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries, J. Comput. Phys., 227, 4825-4852 (2008).
D. Kim and H. Choi, A Second-Order Time-Accurate Finite Volume Method for Unsteady Incompressible Flow on Hybrid Unstructured Grids, J. Comput. Phys., 162, 411-428 (2000).
Y. Zang, R. L. Street, and J. R. Koseff, A non-staggered grid, fractional step method for time-dependent incompressible Navier–Stokes equations in curvilinear coordinates, J. Comput. Phys., 114, 18-33 (1994).

CFD