계면을 어떻게 다룰 것인가#

압축성 다상유동에서 막히는 자리는 플럭스가 아니다. 두 유체가 맞닿은 면을 셀 안에 어떻게 적어 둘 것인가.

갈래는 세 개다.

1. Volume of Fluid (VOF)
2. Level Set
3. Diffuse Interface (확산 계면)

모드 토글을 VOF, Level Set, Diffuse로 바꾸면 같은 원의 두께와 부호가 달라진다.

VOF는 셀 안의 직선 조각으로 색을 가둔다. Level Set은 영점 등고선이 매끈하고, Diffuse는 띠가 몇 칸으로 번진다.

1. Volume of Fluid (VOF)#

원리#

각 셀에서 특정 유체가 차지하는 체적 분율(volume fraction) $\alpha$를 따라간다.

$$\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = 0$$

$\alpha = 1$이면 그 셀은 유체 A로 가득하다. $\alpha = 0$이면 유체 B다. $0 < \alpha < 1$이면 계면이 셀 안을 지난다.

계면 재구성 (Interface Reconstruction)#

VOF가 하는 일은 $\alpha$에서 셀 안의 면 위치를 다시 그리는 일이다.

• SLIC (Simple Line Interface Calculation): 계면을 좌표축에 평행하게 근사
• PLIC (Piecewise Linear Interface Calculation): 계면을 임의 기울기의 직선으로 근사

PLIC에서 계면의 법선 벡터 $\mathbf{n}$은 $\alpha$ 필드의 기울기로 추정합니다:

$$\mathbf{n} = \frac{\nabla \alpha}{|\nabla \alpha|}$$

장단점#

• 장점: 질량이 정확히 보존됨 (보존형)
• 단점: 계면 재구성이 복잡하고 비용이 큼. 3D 확장이 어려움. 곡률 계산 부정확

2. Level Set Method#

원리#

계면은 부호 거리 함수(signed distance function) $\phi$의 영점 등고선이다.

$$\frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi = 0$$

• $\phi > 0$: 유체 A
• $\phi < 0$: 유체 B
• $\phi = 0$: 계면

Re-initialization#

이류(advection)를 지나면 $\phi$는 부호 거리의 성질($|\nabla\phi| = 1$)을 잃는다. 되돌리려고 re-initialization 방정식을 푼다.

$$\frac{\partial \phi}{\partial \tau} + \text{sign}(\phi_0)(|\nabla\phi| - 1) = 0$$

여기서 $\tau$는 가상 시간(pseudo-time)입니다.

곡률과 표면장력#

Level Set은 곡률 같은 기하량을 그대로 꺼내기 쉽다.

$$\kappa = -\nabla \cdot \left(\frac{\nabla\phi}{|\nabla\phi|}\right)$$

표면장력은 CSF(Continuum Surface Force) 모델로 체적력으로 변환됩니다:

$$\mathbf{f}_s = \sigma \kappa \delta(\phi) \nabla\phi$$

장단점#

• 장점: 곡률/법선 계산 용이. 위상 변화(topology change) 자연스럽게 처리
• 단점: 질량 보존 실패. Re-initialization 과정에서 계면 위치가 이동 가능

CLSVOF: 두 방법의 결합#

VOF의 질량 보존 + Level Set의 기하학적 장점을 합친 Coupled Level Set and VOF (CLSVOF) 방법이 많이 사용됩니다.

3. Diffuse Interface Method (확산 계면)#

철학의 전환#

앞의 두 방법이 계면을 날카로운(sharp) 것으로 보는 반면, 확산 계면법은 계면이 본질적으로 유한한 두께를 갖는다고 봅니다.

이 접근에서는 별도의 계면 추적 방정식이 필요 없고, 체적 분율이 보존 법칙 시스템의 일부로 포함됩니다.

5-equation Model (Allaire et al., 2002; Kapila et al., 2001)#

가장 널리 사용되는 reduced model입니다:

$$\frac{\partial (\alpha_1 \rho_1)}{\partial t} + \nabla \cdot (\alpha_1 \rho_1 \mathbf{u}) = 0$$ $$\frac{\partial (\alpha_2 \rho_2)}{\partial t} + \nabla \cdot (\alpha_2 \rho_2 \mathbf{u}) = 0$$ $$\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + p\mathbf{I}) = 0$$ $$\frac{\partial E}{\partial t} + \nabla \cdot ((E + p)\mathbf{u}) = 0$$ $$\frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 = 0$$

마지막 줄이 갈린다. $\alpha_1$의 이류가 **비보존형(non-conservative)**이다. 계면 압력이 흔들리지 않게 하는 자리가 여기다.

왜 비보존형이 필요한가?#

$\alpha_1$을 보존형으로 풀면:

$$\frac{\partial \alpha_1}{\partial t} + \nabla \cdot (\alpha_1 \mathbf{u}) = \alpha_1 \nabla \cdot \mathbf{u}$$

우변의 $\alpha_1 \nabla \cdot \mathbf{u}$가 이산에서 0이 아니면 압력 진동이 난다.

비보존 이류 $\partial\alpha/\partial t + \mathbf{u}\cdot\nabla\alpha = 0$을 전개하면 $\partial\alpha/\partial t + \nabla\cdot(\alpha\mathbf{u}) - \alpha\nabla\cdot\mathbf{u} = 0$이다. 우변으로 넘긴 것이 위 식이다. 접촉면에서 압력과 속도가 일정하면 $\nabla\cdot\mathbf{u}$는 0이어야 한다. 이산 발산이 그 0을 못 내면 우변이 남는다. $\alpha=0.5$인 계면 셀에서 잔여 발산을 $\delta$라 하면 우변은 $0.5\delta$다. $\alpha$가 그 양만큼 움직이면 혼합 밀도가 바뀌고, 상태방정식은 같은 압력으로 돌아가지 못한다. 진동의 씨앗이 이 항이다.

비보존형은 그 우변을 만들지 않는다. $\alpha$는 $\mathbf{u}$를 따라 움직이고, 셀 면의 플럭스 차분이 $\alpha$를 직접 증감시키지 않는다. 5방정식의 마지막 줄이 비보존인 이유가 여기다. 질량은 위의 두 줄, $\partial(\alpha_k\rho_k)/\partial t + \nabla\cdot(\alpha_k\rho_k\mathbf{u})=0$이 따로 지킨다. $\alpha$를 보존형으로 쓰지 않아도 각 상의 질량은 그 줄이 맡는다.

재초기화 쪽의 숫자도 이미 식에 있다. $|\nabla\phi|=1$이 깨지면 $\mathrm{sign}(\phi_0)(|\nabla\phi|-1)$이 가상시간 $\tau$로 $\phi$를 민다. $|\nabla\phi|$가 1보다 크면 그 항이 $\phi$를 깎고, 1보다 작으면 반대로 민다. 멈추는 곳은 $|\nabla\phi|=1$뿐이다. 그 과정에서 영점 $\phi=0$이 한 셀이라도 미끄러지면 질량이 샌다. VOF의 $\alpha$ 적분에는 그 미끄러짐이 없다. 대신 $\mathbf{n}=\nabla\alpha/|\nabla\alpha|$로 뽑은 법선은 $\alpha$가 0과 1 사이로 번진 칸에서만 기울기가 있다. 곡률이 그 칸의 노이즈를 탄다. $\kappa=-\nabla\cdot(\nabla\phi/|\nabla\phi|)$는 $\phi$가 거리함수일 때 같은 법선을 매끈하게 준다. 그래서 CLSVOF가 둘을 붙인다. 질량은 $\alpha$, 곡률은 $\phi$.

7-equation Model (Baer-Nunziato type)#

각 상에 압력과 속도를 따로 줄 수 있는 완전 비평형 모델이다. 그 길은 Baer–Nunziato 7-방정식이 풀어 둔다.

$$\frac{\partial \alpha_k \rho_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k \mathbf{u}_k) = 0, \quad k = 1, 2$$

이완(relaxation) 과정을 통해 압력/속도 평형으로 수렴시킵니다. 수학적으로 가장 완전하지만 계산 비용이 높습니다.

비교 요약#

보존형 α를 쓰면 접촉면의 압력이 흔들린다#

압축성에서 표를 고르면 아래 두 줄이다. 5방정식은 질량을 $\partial(\alpha_k \rho_k)/\partial t + \nabla\cdot(\alpha_k \rho_k \mathbf{u})=0$으로 지키고, $\alpha$ 이류만 비보존으로 남긴다. 이산 발산이 $\alpha\nabla\cdot\mathbf{u}$를 0으로 만들지 못하면, $\alpha=0.5$인 계면 셀이 혼합 밀도를 밀어 압력을 흔든다.

번진 띠를 매 스텝 되돌릴지는 다른 결정이다. 그 조작은 안티확산으로 계면을 다시 세우는 구현에 이어진다.