채널 유동 수치해석 (Channel Flow DNS)

3차원 비압축성 채널 유동(channel flow) 을 직접 수치해석(DNS)하기 위한 솔버를 정리한 글이다. 큰 줄기는 다음 한 줄로 요약된다 — 지배방정식 무차원화 → Projection method → (Adams–Bashforth + Crank–Nicolson) 시간적분 → (Spectral + FDM) 공간 이산화 → 비균일·엇갈림(staggered) 격자 → TDMA → 압력 Poisson.

🌊

핵심 설계 — 벽에 평행한 \(x\)(주류)·\(z\)(스팬) 방향은 주기적이므로 스펙트럴(Fourier) 로 정확히 미분하고, 벽에 수직한 \(y\) 방향만 유한차분(FDM) + 비균일 격자로 다룬다. 비선형 항은 explicit(Adams–Bashforth), 점성 항은 implicit(Crank–Nicolson)으로 처리해 안정성과 정확도를 함께 잡는다.

채널 유동: 위·아래 벽 사이의 유동. \(x,z\) 는 주기적, \(y\) 는 벽으로 막혀 있다

1. 지배방정식과 무차원화

채널 유동(Channel flow)의 Navier–Stokes 운동량 방정식부터 시작한다. Momentum equation 은 벡터형과 성분형으로 다음과 같다.

\[ \rho\left(\frac{\partial \mathbf u}{\partial t} + \mathbf u\cdot(\nabla \mathbf u)\right) = \rho \mathbf g + \nabla p + \mu\nabla^2 \mathbf u \;\cdots\,① \]

\[ \rho\left(\frac{\partial \mathbf u}{\partial t} + u\frac{\partial \mathbf u}{\partial x} + v\frac{\partial \mathbf u}{\partial y} + w\frac{\partial \mathbf u}{\partial z}\right) = \rho \mathbf g + \frac{\partial p}{\partial x_j} + \mu\nabla^2 \mathbf u \]

무차원화(nondimensionalization) 를 한다. 대표 길이 \(L\), 대표 속도 \(U\) 를 기준으로 \(\dfrac{u}{U}=u^{*},\ \dfrac{x}{L}=x^{*}\) 로 정규화하여 대입하면

\[ \rho\left(U\frac{\partial u^*}{\partial t} + \frac{U^2}{L}u^*\frac{\partial u^*}{\partial x^*} + \frac{U^2}{L}v^*\frac{\partial u^*}{\partial y^*} + \frac{U^2}{L}w^*\frac{\partial u^*}{\partial z^*}\right) = \rho \mathbf g + \frac{\partial p}{\partial x_j} + \frac{U}{L^2}\mu\frac{\partial^2 u^*/U \cdot U}{\partial x_j^{*2}} \]

\[ \rho\left(U\frac{\partial u^*}{\partial t} + \frac{U^2}{L}\,u^*\cdot(\nabla u^*)\right) = \rho \mathbf g + \frac{1}{L}\frac{\partial p}{\partial x_j} + \frac{U}{L^2}\mu\,\nabla^2 u^* \;\cdots\,② \]

②식의 양변에 \(\dfrac{L}{U^2}\) 를 곱해 정리하면 — Result(nondimensionalization: Velocity, length) —

\[ \frac{L}{U}\frac{\partial u^*}{\partial t} + u^*\cdot(\nabla u^*) = \frac{L}{U^2}\mathbf g + \frac{1}{\rho U^2}\frac{\partial p}{\partial x_j} + \frac{1}{Re}\nabla^2 u^* \]

\[ \frac{\partial u^*}{\partial(tU/L)} + u^*\cdot(\nabla u^*) = g^* + \frac{\partial p^*}{\partial x_j} + \frac{1}{Re}\nabla^2 u^* \;\cdots\,③ \]

즉 단 하나의 무차원 수 레이놀즈 수 만 남는다.

\[ t^* = \frac{tU}{L}, \qquad g^* = \frac{g}{U^2/L}, \qquad p^* = \frac{p}{\rho U^2}, \qquad Re = \frac{UL}{\nu} \]

정리된 무차원 운동량 방정식 ③ 을 이제 Numerical method 로 풀어 Solution 을 구해야 한다(이후 \(^*\) 는 생략).

\[ \frac{\partial \mathbf u}{\partial t} + \mathbf u\cdot(\nabla\mathbf u) = \mathbf g + \nabla p + \frac{1}{Re}\nabla^2\mathbf u, \qquad \nabla\cdot\mathbf u = 0 \]

2. Projection method

비압축성 유동의 난점은 압력–속도 결합(closure problem) 이다. 압력에 대한 독립 방정식이 없기 때문이다. Projection method 는 이를 2단계로 분리한다.

Projection method: 압력 없이 중간속도를 구하고, Poisson으로 의사압력을 푼 뒤, 발산이 0이 되도록 보정

먼저 압력을 빼고 중간속도(intermediate velocity) \(\hat{\mathbf u}\) 를 구한 뒤, 의사압력(pseudo pressure) \(\phi\) 에 대한 Poisson 방정식을 풀고, 마지막에 보정해 발산이 0인 속도장을 얻는다.

\[ \nabla^2 \phi = \frac{1}{\Delta t}\,\nabla\cdot\hat{\mathbf u}, \qquad \mathbf u^{\,n+1} = \hat{\mathbf u} - \Delta t\,\nabla\phi \]

3. 시간 적분 — Adams–Bashforth + Crank–Nicolson

Intermediate velocity 를 구하기 위해 ④식을 우선 풀어야 한다. 텐서 표기(tensor notation)로 쓰면 — 정상유체라 증명은 되었으나 편의상 위와 같이 표현한다(참고) —

\[ \frac{\partial \hat u_i}{\partial t} + \hat u_j\frac{\partial}{\partial x_j}(\hat u_i \hat u_j) = \frac{1}{Re}\frac{\partial^2}{\partial x_j^2}\hat u_i \]

여기서 항의 성질에 따라 시간 적분 방법을 다르게 쓴다 — 비선형(대류) 항은 Adams–Bashforth, 점성(viscous) 항은 Crank–Nicolson 으로 차분한다. 시간 항은 Euler 로 이산화한다.

\[ \frac{\hat u_i - u_i^{\,t}}{\Delta t} + \underbrace{\left(\tfrac{3}{2}\frac{\partial}{\partial x_j}(u_i^{t}u_j^{t}) - \tfrac{1}{2}\frac{\partial}{\partial x_j}(u_i^{t-1}u_j^{t-1})\right)}_{\text{Adams–Bashforth (대류)}} = \underbrace{\frac{1}{Re}\frac{\partial^2}{\partial x_j^2}\,\tfrac{1}{2}\big(\hat u_i + u_i^{t}\big)}_{\text{Crank–Nicolson (점성)}} \]

증분 \(\delta u_i = \hat u_i - u_i^{\,t}\) 로 치환하여 \(\delta u_i\) 에 대해 정리하면,

\[ \frac{\delta u_i}{\Delta t} = -\left(\tfrac{3}{2}\frac{\partial}{\partial x_j}(u_i^{t}u_j^{t}) - \tfrac{1}{2}\frac{\partial}{\partial x_j}(u_i^{t-1}u_j^{t-1})\right) + \frac{1}{2Re}\frac{\partial^2}{\partial x_j^2}\delta u_i + \frac{1}{Re}\frac{\partial^2}{\partial x_j^2}u_i^{t} \;\cdots\,⑤ \]

\(\delta u_i\) 에 대해 한 번 더 정리하면 LHS(미지수 \(\delta u_i\) 의 implicit 연산자)와 RHS(직전 시점으로 계산되는 \(f_{(x,y,z)}\))로 깔끔히 나뉜다.

\[ \underbrace{\left(1 - \frac{\Delta t}{2Re}\frac{\partial^2}{\partial x_j^2}\right)\delta u_i}_{\text{LHS}} = \underbrace{-\frac{\Delta t}{2}\left(3\frac{\partial}{\partial x_j}(u_i^{t}u_j^{t}) - \frac{\partial}{\partial x_j}(u_i^{t-1}u_j^{t-1})\right) + \frac{\Delta t}{Re}\frac{\partial^2}{\partial x_j^2}u_i^{t}}_{\text{RHS} = f_{(x,y,z)}} \;\cdots\,⑥ \]

4. 공간 이산화 — Spectral(x, z) + FDM(y)

방향별로 가장 알맞은 이산화를 섞어 쓴다.

FDM(Finite difference method) — Factorization 과 TDMA

⑥식의 좌변(LHS)을 \(y\) 뿐 아니라 모든 방향을 FDM 으로 푼다고 하면, 라플라시안 연산자를 방향별로 factorization 한다.

\[ \left(1 - \frac{\Delta t}{2Re}\frac{\partial^2}{\partial x_j^2}\right)\delta u_i \;\approx\; \left(1 - \frac{\Delta t}{2Re}\frac{\partial^2}{\partial x_1^2}\right)\left(1 - \frac{\Delta t}{2Re}\frac{\partial^2}{\partial x_2^2}\right)\left(1 - \frac{\Delta t}{2Re}\frac{\partial^2}{\partial x_3^2}\right)\delta u_i \quad(\text{Factorization}) \]

Factorization 에 의해 생기는 Error 는 어쩔 수 없다. 2D case 로 곱을 풀어 보면, 원래 식에 없던 교차항이 Error term 으로 더해진다.

\[ \left(1 - \tfrac{\Delta t}{2Re}\tfrac{\partial^2}{\partial x_1^2}\right)\left(1 - \tfrac{\Delta t}{2Re}\tfrac{\partial^2}{\partial x_2^2}\right) = 1 - \tfrac{\Delta t}{2Re}\tfrac{\partial^2}{\partial x_1^2} - \tfrac{\Delta t}{2Re}\tfrac{\partial^2}{\partial x_2^2} + \underbrace{\frac{\Delta t^2}{4Re^2}\frac{\partial^2}{\partial x_1^2}\frac{\partial^2}{\partial x_2^2}}_{\text{Error term}} \]

TDMA 를 풀기 위한 꼴로 factorization 한 식을 재정리하면 다음과 같다. 각 인수를 \(g_{(x,y,z)}\), \(h_{(x,y,z)}\), \(\delta u_i\) 로 두고 순서대로 푼다.

\[ \underbrace{\left(1 - \tfrac{\Delta t}{2Re}\tfrac{\partial^2}{\partial x_1^2}\right)}_{g_{(x,y,z)}}\underbrace{\left(1 - \tfrac{\Delta t}{2Re}\tfrac{\partial^2}{\partial x_2^2}\right)}_{h_{(x,y,z)}}\left(1 - \tfrac{\Delta t}{2Re}\tfrac{\partial^2}{\partial x_3^2}\right)\delta u_i = -\frac{\Delta t}{2}\left(3\frac{\partial}{\partial x_j}(u_i^{t}u_j^{t}) - \frac{\partial}{\partial x_j}(u_i^{t-1}u_j^{t-1})\right) + \frac{\Delta t}{Re}\frac{\partial^2}{\partial x_j^2}u_i^{t} \;\cdots\,⑦ \]

우변 \(\text{RHS}=f_{(x,y,z)}\). Protocol 은 다음과 같다 — 1ˢᵗ TDMA 로 \(g_{(x,y,z)}\) → 2ⁿᵈ TDMA 로 \(h_{(x,y,z)}\) → 3ʳᵈ TDMA 로 \(\delta u_i{}_{(x,y,z)}\) 를 구한다. 즉 위의 순서대로 풀어 \(\delta u\) 를 구하면 되고, 이후 아래의 Poisson equation 을 풀어 다음 step 을 구한다.

\[ \text{Divergence}\left(\frac{\partial u^{\,t+0.5}}{\partial t^*} = \nabla\phi^{*}\right) \;\Longrightarrow\; \frac{\partial}{\partial t}\frac{\partial u_j^{\,t+1}}{\partial x_j} = \nabla^2\phi_{(x,y,z)} \quad(\text{Poisson equation}) \]

Divergence 를 취하는 이유는, Navier–Stokes 를 풀 때 우리는 momentum eqn 을 풀지만 divergence-free condition(continuity equation) 으로 맞춰주지 않으면 안 되기 때문이다. 이에 따라 위 momentum eqn 에 Divergence 를 취해 풀어낸다. ⑦식에 대해 \(x,y,z\) 각 방향으로 Central difference method 등을 사용해 차분하고 풀면 그것이 곧 Finite difference method 가 된다(앞서 Cavity 에서 풀어보았으니 참조).

Spectral method 를 간단히 보자. Fourier transform 은 모든 물리량을 파동으로 이루어져 있다고 본다.

\[ f(x) = \int_{-\infty}^{\infty} F(\alpha)\,e^{i\alpha x}\,d\alpha, \qquad F(\alpha) = \int_{-\infty}^{\infty} f(x)\,e^{-i\alpha x}\,dx \qquad(\alpha:\ \text{wave number}) \]

이를 Discretization 시키면 \(f(x)\) 는 각 basis 파동 \(\hat f_n(x)\) 들의 합으로 표현된다.

\[ f(x) = \sum_{n} F(\alpha_n)e^{i\alpha_n x} = \hat f_1(x) + \hat f_2(x) + \hat f_3(x) + \cdots \]

미분하면 FDM 과 달리 Spectral 은 exact 하게 대수적으로 떨어진다(미분 → 곱셈).

\[ \frac{df(x)}{dx} = \sum_{n} i\alpha_n F(\alpha_n)e^{i\alpha_n x} \;\Longrightarrow\; \frac{d\hat f_n}{dx} = i\alpha_n \hat f_n \;\Longrightarrow\; \frac{d^2\hat f_n}{dx^2} = -\alpha_n^2\,\hat f_n \]

앞의 Navier–Stokes 를 구성하는 하나의 Fourier mode \(\hat u_n\) 으로 해석하면, \(x,z\) 방향(파수 \(\alpha,\gamma\))의 2계 미분이 다음처럼 떨어진다.

\[ \frac{\partial^2(\delta\hat u_i)}{\partial x_1^2} = -\alpha^2\,\delta\hat u_i, \qquad \frac{\partial^2(\delta\hat u_i)}{\partial x_3^2} = -\gamma^2\,\delta\hat u_i \]

이를 ⑥ 식에 넣으면 \(x,z\) 의 라플라시안이 상수로 떨어져, 실제로 차분해야 할 방향은 \(y\) 하나만 남는다.

\[ \left(1 + \frac{\Delta t}{2Re}\alpha^2 + \frac{\Delta t}{2Re}\gamma^2 - \frac{\Delta t}{2Re}\frac{\partial^2}{\partial x_2^2}\right)\delta\hat u_i = f_{(x,y,z)} \;\cdots\,⑧ \]

덕분에 \(x,z\) 방향에는 차분 오차나 factorization 이 필요 없다(FDM 으로 풀 경우엔 \(\big(1-\tfrac{\Delta t}{2Re}\partial_{x_1}^2\big)\big(1-\tfrac{\Delta t}{2Re}\partial_{x_2}^2\big)\big(1-\tfrac{\Delta t}{2Re}\partial_{x_3}^2\big)\) 로 factorization 한 뒤 방향마다 TDMA 를 풀어야 하고, factorization 에서 \(\tfrac{\Delta t^2}{4Re^2}\) 차수의 Error term 이 생긴다). 다만 비선형 항을 스펙트럴로 계산하면 aliasing 오차가 생기므로 곱(convolution)을 구할 때 zero-padding(2/3 law) 으로 제거한다. 반면 \(y\) 방향은 벽 때문에 주기적이지 않아 FDM 으로 차분한다.

5. 비균일 격자 (벽 근처 클러스터링)

벽 근처에서는 경계층의 급격한 변화를 잡기 위해 해상도를 높여야 한다. 그래서 \(y\) 방향은 Chebyshev 형 의 비균일 격자(벽 쪽으로 촘촘)를 쓴다.

비균일 격자: 벽(위·아래)으로 갈수록 격자가 촘촘해진다 (Chebyshev 분포)

간격이 일정하지 않으므로 미분도 Taylor expansion 으로 다시 유도해야 한다. 인접 간격을 \(dy_1\)(아래), \(dy_2\)(위)라 하고 점 \(y\) 주변을 전개하면

\[ f(y+dy_2) = f(y) + dy_2 f'(y) + \frac{dy_2^{\,2}}{2!}f''(y) + \frac{dy_2^{\,3}}{3!}f'''(y) + \cdots \;\cdots\,⑨ \]

\[ f(y-dy_1) = f(y) - dy_1 f'(y) + \frac{dy_1^{\,2}}{2!}f''(y) - \frac{dy_1^{\,3}}{3!}f'''(y) + \cdots \;\cdots\,⑩ \]

2계 미분 \(f''(y)\)

\(f'(y)\) 항을 소거하기 위해 \(\dfrac{dy_1}{dy_2}\times⑨ + ⑩\) 을 만든다.

\[ \frac{dy_1}{dy_2}f(y+dy_2) + f(y-dy_1) = \left(\frac{dy_1}{dy_2}+1\right)f(y) + \left(\frac{dy_1 dy_2}{2!}+\frac{dy_1^{\,2}}{2!}\right)f''(y) + \left(\frac{dy_1 dy_2^{\,2}}{3!}-\frac{dy_1^{\,3}}{3!}\right)f'''(y) + \cdots \]

\(f''(y)\) 의 계수는 \(\dfrac{dy_1}{2}(dy_1+dy_2)\) 이므로, 이에 대해 정리하면 2차 오차(second order error) 의 비균일 중심차분 2계 미분을 얻는다.

\[ f''(y) = \frac{2}{dy_2(dy_1+dy_2)}f(y+dy_2) - \frac{2}{dy_1 dy_2}f(y) + \frac{2}{dy_1(dy_1+dy_2)}f(y-dy_1) + O(dy) \;\cdots\,⑪ \]

1계 미분 \(f'(y)\)

반대로 \(⑨-⑩\) 으로 \(f''(y)\) 를 일부 소거하면 1계 미분이 나온다.

\[ f(y+dy_2) - f(y-dy_1) = (dy_1+dy_2)f'(y) + \frac{dy_2^{\,2}-dy_1^{\,2}}{2!}f''(y) + \cdots \]

\[ f'(y) = \frac{f(y+dy_2) - f(y-dy_1)}{dy_1+dy_2} - \frac{dy_2^{\,2}-dy_1^{\,2}}{2(dy_1+dy_2)}f''(y) + \cdots \quad(\text{first order error}) \]

격자가 균일(\(dy_1=dy_2=dy\)) 해지면 \(dy_2^{\,2}-dy_1^{\,2}=0\) 이 되어 1차 오차항이 사라지고, 익숙한 2차 정확도 중심차분으로 돌아간다.

\[ f'(y) = \frac{f(y+dy) - f(y-dy)}{2\,dy} + O(dy^2), \qquad f''(y) = \frac{f(y+dy) - 2f(y) + f(y-dy)}{dy^2} + O(dy^2) \]

이렇게 어떤 grid 를 쓰든 points 사이의 관계를 Taylor expansion 으로 대수적으로 나타내는 작업을 마무리했다. 이 차분식이 곧 \(y\) 방향 FDM 의 계수가 된다.

6. 격자 배치 — Collocated vs Staggered

속도와 압력을 격자 위 어디에 둘지 가 또 하나의 핵심이다. 흔히 빨간 dashed-line 을 staggered grid, 파란 선을 collocated grid 라고 오해 하기 쉬운데, 우리가 쓰는 격자는 파란 선 하나다. 둘의 진짜 차이는 “속도 성분과 압력이 한 지점에 몰려 있느냐(Coupling) 아니면 따로 떨어져 있느냐(Decoupling)” 일 뿐이다.

속도·압력이 한 점 → 결합은 쉽지만 압력 checkerboard 디커플링 위험

u 는 세로 면(→), v 는 가로 면(↑), p 는 셀 중심 → 디커플링 회피

3D 로 보면, 한 점에 세 속도 성분이 모이면 Coupling, 면마다 흩어지면 Decoupling 이다.

한 점에 세 속도가 모이면 Coupling, 면마다 흩어지면 Decoupling

이 솔버는 두 방식을 섞는다. 앞서 staggered grid 는 FDM 을 쓰기 때문에 \(u_1,u_3\) 가 decoupling 되어 있었는데, \(x,z\) 방향은 스펙트럴이라 exact 한 미분이 가능 하므로 \(u_1,u_3\) 를 압력 위치(cell center)에 직접 coupling 시킨다. 결과적으로 \(u_1,u_3\) 는 cell center 에 coupling 되고, FDM 으로 풀어야 하는 \(u_2(=v)\) 방향만 엇갈림(staggered) 으로 둔다 — 즉 “staggered grid with spectral” 구성이다. 벽에서는 no-slip 을 위해 ghost cell 을 두고(\(u_1=-u_0\) 등), 액추에이터(blower/suction)가 있을 때는 벽면 \(v\) 에 그 값을 준다.

7. 선형 시스템 — TDMA

스펙트럴 방향은 미분이 대수적으로 떨어지므로, 실제로 풀어야 하는 것은 \(y\) 방향 한 줄 뿐이다. 앞서 구한 중앙차분(3개의 점이 필요)을 grid 에 적용하면, FDM 의 인접 3점이 만드는 삼중대각(tridiagonal) 계를 TDMA(Thomas algorithm) 로 한 번에 푼다.

\[ a\,\delta\hat u_{j-1} + b\,\delta\hat u_{j} + c\,\delta\hat u_{j+1} = f_{(x,y,z)} \]

\[ a = -\frac{\Delta t}{2Re}\frac{2}{dy_j(dy_j+dy_{j-1})}, \quad b = 1 + \frac{\Delta t}{2Re}\alpha^2 + \frac{\Delta t}{2Re}\gamma^2 - \frac{\Delta t}{2Re}\Big(\frac{1}{dy_j}+\frac{1}{dy_{j-1}}\Big)\frac{2}{dy_j+dy_{j-1}}, \quad c = -\frac{\Delta t}{2Re}\frac{2}{dy_j(dy_j+dy_{j-1})} \]

[Figure 01 — staggered & collocated grid] \(y\) 축은 총 129 개 점(128 구간)으로 둔다. 변수 개수는 \(u,w\): 130 개, \(v\): 129 개다. collocated grid 기준 위·아래 2칸의 ghost cell 을 위해 다음 ghost cell interval 조건과 no-slip 경계조건, 그리고 staggered–collocated 위치 관계식을 둔다.

\[ dy_0 = dy_1,\quad dy_{129} = dy_{128} \qquad(\text{ghost interval}) \]

\[ u_0 = -u_1,\quad u_{128} = -u_{129} \qquad(\text{no-slip}) ,\qquad dy_k^{*} = \frac{dy_k + dy_{k+1}}{2} \]

TDMA solver — Matrix form (Thomas algorithm). \(u,w\)(staggered grid 기준)와 \(v\)(collocated grid 기준)에 대해 각각 삼중대각 행렬을 만들어 푼다. 경계 행에는 경계조건이 반영된다.

\[ \begin{pmatrix} b_1 & c_1 & & \\ a_2 & b_2 & c_2 & \\ & \ddots & \ddots & \ddots \\ & & a_{n} & b_{n} \end{pmatrix} \begin{pmatrix} \delta\hat u_1 \\ \delta\hat u_2 \\ \vdots \\ \delta\hat u_{n} \end{pmatrix} = \begin{pmatrix} f_1 \\ f_2 \\ \vdots \\ f_{n} \end{pmatrix} \]

\(u,w\) 는 ghost cell 을 고려한 no-slip(\(u_{y'=0}=-u_{y'=-0}\) 등)을, \(v\) 는 액추에이터 유무(No actuator: 양 벽 \(v=0\) / With actuator: \(v_0=\) bottom act, \(v_{128}=\) top act)를 첫·마지막 행에 반영한다. \(u,w\) point 는 staggered grid 위에, \(v\) point 는 collocated grid 위에 있으므로 위치에 맞게 보간(interpolation)하여 적용한다.

이산 RHS \(f_{(x,y,z)}\) — Fourier · 차분 명시형

⑥·⑦에서 기호로만 두었던 우변 \(\text{RHS}=f_{(x,y,z)}\) 를 실제로 푸는 형태로 펼치면 다음과 같다. \(x,z\) 는 스펙트럴이라 \(\partial_{x_1}\!\to i\alpha\), \(\partial_{x_3}\!\to i\gamma\) 로 떨어지고, \(y\)(=\(x_2\)) 방향만 차분으로 남는다. 대류항은 AB(3·\(t\) − \(t{-}1\)), 점성항은 CN 으로 계산한다.

\[ \begin{aligned} f_{(x,y,z)} = &-\frac{\Delta t}{2}\Bigg[\,3\Big( i\alpha\,\widehat{(u_1u_j)}^{\,t} + \tfrac{\partial}{\partial x_2}\widehat{(u_2u_j)}^{\,t} + i\gamma\,\widehat{(u_3u_j)}^{\,t}\Big) - \Big( i\alpha\,\widehat{(u_1u_j)}^{\,t-1} + \tfrac{\partial}{\partial x_2}\widehat{(u_2u_j)}^{\,t-1} + i\gamma\,\widehat{(u_3u_j)}^{\,t-1}\Big)\Bigg] \\[4pt] &+\frac{\Delta t}{Re}\Big(\,-\alpha^2\,\hat u_i^{\,t} + \big(A\,\hat u_{(y'-dy_k)} + B\,\hat u_{(y')} + C\,\hat u_{(y'+dy_{k+1})}\big)^{t} - \gamma^2\,\hat u_i^{\,t}\Big) \end{aligned} \]

여기서 \(y\) 방향 2계 미분 계수 \(A,B,C\) 는 비균일 격자의 중심차분 계수(앞 5절 ⑪)와 같다.

\[ A = \frac{2}{dy_k(dy_{k+1}+dy_k)}, \qquad B = -\Big(\frac{1}{dy_{k+1}}+\frac{1}{dy_k}\Big)\frac{2}{dy_{k+1}+dy_k}, \qquad C = \frac{2}{dy_{k+1}(dy_{k+1}+dy_k)} \]

비선형 항을 스펙트럴 공간에서 곱하면 aliasing 이 생기므로, 곱(convolution) \( \widehat{(u_iu_j)} \) 를 구할 때 zero-padding(2/3 law) — 고파수 1/3 을 0 으로 채워 — aliasing 을 제거한다.

Staggered ↔ Collocated 보간 (2차 정확도 유지)

\(u,w\)(staggered)와 \(v\)(collocated)는 격자 위치가 어긋나 있으므로, 대류항의 \(y\)-미분 \(\partial_{x_2}\widehat{(u_2u_j)}\) 를 “우리가 원하는 point”에서 평가하려면 보간이 필요하다. 한 칸 앞·뒤 차분을 각각 구하면 둘 다 1차이지만,

\[ \frac{\partial}{\partial x_2}\widehat{(u_2u_j)}\Big|_{y'} = \frac{1}{dy_{k+1}}\Big(\widehat{(u_2u_j)}_{(y'+dy_{k+1})} - \widehat{(u_2u_j)}_{(y')}\Big) + O(dy^2) \]

\[ \frac{\partial}{\partial x_2}\widehat{(u_2u_j)}\Big|_{y'} = \frac{1}{dy_{k}}\Big(\widehat{(u_2u_j)}_{(y')} - \widehat{(u_2u_j)}_{(y'-dy_{k})}\Big) + O(dy^2) \]

두 식을 각 \(0.5\) 씩 더해(평균) 합치면, 1차 오차가 상쇄되어 2차 정확도(Second order error) 가 유지된다.

\[ \frac{\partial}{\partial x_2}\widehat{(u_2u_j)}\Big|_{y'} = \frac{0.5}{dy_{k+1}}\Big(\widehat{(u_2u_j)}_{(y'+dy_{k+1})} - \widehat{(u_2u_j)}_{(y')}\Big) + \frac{0.5}{dy_{k}}\Big(\widehat{(u_2u_j)}_{(y')} - \widehat{(u_2u_j)}_{(y'-dy_{k})}\Big) + O(dy^2) \]

참고로 위 Interpolation form 대신 Taylor expansion 으로 직접 구한 form 을 쓰면, 비균일 격자에서 난류가 잘 안 잡히는 경우가 있어 보간형을 쓴다.

8. 압력 — Poisson equation

Intermediate velocity 를 구했으니, 다음 step 의 속도로 update 하도록 의사압력 \(\phi\) 를 푼다. 발산(Divergence)을 취한 ⑤식에서 Poisson equation 이 나온다.

\[ \frac{\partial}{\partial t}\frac{\partial u_j^{\,t+1}}{\partial x_j} = \nabla^2\phi_{(x,y,z)} \;\Longrightarrow\; \nabla^2\phi = \frac{1}{\Delta t}\left(\frac{\partial \hat u_1}{\partial x_1}+\frac{\partial \hat u_2}{\partial x_2}+\frac{\partial \hat u_3}{\partial x_3}\right) \]

여기서도 \(x,z\) 는 스펙트럴, \(y\) 는 FDM 이다. Fourier mode 로 바꾸면 \(x,z\) 의 라플라시안이 \(-\alpha^2,-\gamma^2\) 로 떨어진다. 압력은 셀 중심(cell center) 에 위치한다(staggered 기준).

\[ -\alpha^2\tilde\phi - \gamma^2\tilde\phi + \frac{\partial^2 \tilde\phi}{\partial y^2} = \frac{1}{\Delta t}\left(i\alpha\,\tilde{\hat u}_1 + \frac{\partial \tilde{\hat u}_2}{\partial y} + i\gamma\,\tilde{\hat u}_3\right) \]

\(y\) 방향 2계 미분을 앞의 비균일 중심차분으로 차분하면 다시 \(y\) 방향 삼중대각 계가 되고, TDMA 로 푼다.

\[ A\,\tilde\phi_{j-1} + B\,\tilde\phi_{j} + C\,\tilde\phi_{j+1} = D \]

마지막으로 \(\mathbf u^{\,n+1} = \hat{\mathbf u} - \Delta t\,\nabla\phi\) 로 보정하면 발산이 0 인 속도장을 얻는다 — 이렇게 한 time step 이 끝난다.

\[ \frac{\partial^2 \tilde\phi}{\partial y^2} - (a^2+\gamma^2)\,\tilde\phi = \frac{1}{\Delta t}\,\widetilde{(\nabla\cdot\hat{\mathbf u})} \]