계면 폭이 2칸이면 표면장력의 39%가 사라진다 — Korteweg 응력의 비등방성

계면을 가로질러 압력계를 옮겨 보았다#

정지한 평평한 기액 계면이 하나 있다. 기체 쪽에서 액체 쪽으로 압력계를 천천히 옮겨 본다. 바늘은 움직이지 않는다. 평평한 계면에는 Laplace 압력 차(곡률이 만드는 압력 점프)가 없기 때문이다. 그런데 표면장력은 분명히 있다. 그럼 표면장력은 어디에 숨어 있을까?

답은 압력계를 90도 돌리면 보인다. 계면에 나란한 방향으로 재면 계면 안에서만 바늘이 내려간다. 이 글은 그 차이가 어디서 나오는지 Korteweg 응력(밀도 기울기가 만드는 응력)으로 유도한다. 그리고 격자 간격 1짜리 LBM 격자에서 이 적분이 계면 폭에 따라 얼마나 사라지는지 직접 잰다. 계면 폭이 2칸이면 39%가 없어진다.

정지한 유체의 응력은 압력 하나가 아니다#

다상 LBM의 지배방정식을 쓸 때 운동량 식의 응력은 세 조각으로 나뉜다.

$p$는 열역학 압력, $\boldsymbol{\tau}_K$는 Korteweg 응력, $\boldsymbol{\tau}_{NS}$는 Navier–Stokes 점성 응력이다. $\mathbf{S}^{d}$는 대각합을 뺀 변형률 텐서(trace-free rate-of-strain, 부피 변화를 뺀 순수 전단 변형률)다.

정지 상태에서는 $\mathbf{S}^{d}=0$이므로 점성 응력이 사라진다. 남는 것은 $p$와 $\boldsymbol{\tau}_K$뿐이다. 균일한 유체에서는 $\boldsymbol{\tau}_K$도 0이다. 밀도 기울기가 있는 계면 안에서만 살아난다. 그리고 이 항은 방향을 가린다. 계면에 수직인 방향과 나란한 방향에 서로 다른 값을 준다.

아래 시뮬레이션에서 직접 조작해보자. 노란 점선이 계면을 가로질러 움직이며 그 위치의 유체 요소에 걸리는 응력을 보여준다.

초록 화살표(수직 성분)는 계면 어디서나 점선 원과 같은 길이다. 주황 화살표(접선 성분)는 계면 한가운데서 짧아진다. κ를 키우면 계면이 넓어지고 주황 면적, 곧 σ가 $\sqrt{\kappa}$에 비례해 커진다.

자유에너지에서 Korteweg 텐서로#

확산 계면 모델은 자유에너지 범함수 하나에서 출발한다.

$\psi(\rho)$는 균일 유체의 자유에너지 밀도, $\kappa$는 기울기 에너지 계수다. 이 범함수의 변분에서 압력 텐서가 나온다.

$p_0=\rho\psi'-\psi$는 상태방정식 압력이다. 마지막 항 $\kappa\nabla\rho\otimes\nabla\rho$가 비등방성의 원천이다. 이 텐서의 발산은 $\nabla\cdot\mathbf{P}=\rho\nabla\mu$를 만족한다. $\mu=\psi'(\rho)-\kappa\nabla^2\rho$는 화학 퍼텐셜이다. LBM에서 이 두 형태를 힘으로 넣는 차이는 비이상 LBM forcing의 두 형태에서 다뤘다.

이 글에서는 반 데르발스 상태방정식을 환산 단위로 쓴다.

$\rho$, $T$, $p$ 모두 임계값으로 나눈 값이다. $T<1$이면 기체 밀도 $\rho_g$와 액체 밀도 $\rho_l$이 공존한다. 두 값은 Maxwell 조건(두 상의 $\mu$와 $p$가 같음)으로 정해진다.

수직 성분은 평평하고, 접선 성분은 꺼진다#

계면 법선을 $z$축으로 두자. 밀도는 $z$에만 의존한다. 텐서의 두 성분을 쓰면 이렇다.

$p_N=P_{zz}$는 수직 성분, $p_T=P_{xx}$는 접선 성분이다. 정지 평형에서는 $dP_{zz}/dz=0$이다. 그래서 $p_N$은 계면 전체에서 벌크 압력 $p_b$와 같다. 압력계가 움직이지 않은 이유다.

두 성분의 차이는 $p_N-p_T=\kappa\rho'^2\ge 0$이다. 접선 방향 압력이 계면 안에서만 낮다. 압력이 낮다는 것은 그 층이 옆으로 당기고 있다는 뜻이다. 이것이 표면장력의 역학적 정의다.

평형 방정식 $\kappa\rho''=\psi'(\rho)-\mu_b$에 $\rho'$를 곱해 적분하면 첫 적분이 나온다. $\frac{\kappa}{2}\rho'^2 = W(\rho) \equiv \psi(\rho)-\mu_b\rho+p_b$다. 이를 대입하면 공간 적분이 밀도 적분으로 바뀐다.

$W(\rho)$는 두 벌크 상태 사이의 자유에너지 언덕이다. 이 식에서 $\sigma\propto\sqrt{\kappa}$가 바로 보인다. 격자와 무관한 기준값이다.

코드: dz=1 격자 위의 반 데르발스 계면#

LBM은 격자 간격을 1로 둔다. 그러면 계면 폭은 κ와 $T$로만 정해진다. 아래 코드는 화학 퍼텐셜 잔차를 0으로 보내는 경사 흐름(gradient flow, 잔차 방향으로 조금씩 내려가는 반복)으로 평형 계면을 만든다. 그 뒤 격자 위에서 $p_N$, $p_T$를 계산해 응력 적분 σ와 첫 적분 공식 σ를 비교한다.

import math
 
# 환산 반 데르발스: rho_c = T_c = p_c = 1
def p0(r, T):  return 8*T*r/(3 - r) - 3*r*r
def mu0(r, T): return 8*T/3*(math.log(r/(3 - r)) + r/(3 - r)) - 6*r
def psi(r, T): return r*mu0(r, T) - p0(r, T)
 
def maxwell(T):
    # 화학 퍼텐셜과 압력이 같아지는 점 (Newton, 수치 야코비안)
    rg, rl = 1 - 2*math.sqrt(1 - T), 1 + 2*math.sqrt(1 - T)
    h = 1e-7
    for _ in range(50):
        f1 = mu0(rl, T) - mu0(rg, T)
        f2 = p0(rl, T) - p0(rg, T)
        a = (mu0(rl+h, T) - mu0(rl, T))/h; b = -(mu0(rg+h, T) - mu0(rg, T))/h
        c = (p0(rl+h, T) - p0(rl, T))/h;   d = -(p0(rg+h, T) - p0(rg, T))/h
        det = a*d - b*c
        rl -= ( d*f1 - b*f2)/det
        rg -= (-c*f1 + a*f2)/det
    return rg, rl
 
def korteweg_profile(T, kappa, n=100, dz=1.0, tol=1e-7):
    rg, rl = maxwell(T)
    mub, pb = mu0(rl, T), p0(rl, T)
    L = n*dz                                         # 격자 단위: LB 격자처럼 dz = 1
    z = [(i + 0.5)*dz - L/2 for i in range(n)]
    r = [0.5*(rl + rg) + 0.5*(rl - rg)*math.tanh(zi*math.sqrt(1 - T)/math.sqrt(kappa)) for zi in z]
    dt = 0.2*min(dz*dz/kappa, 1.0)
    for it in range(1, 30001):
        e = [r[0]] + r + [r[-1]]                     # 양 끝 기울기 0
        res = [mu0(r[i], T) - kappa*(e[i+2] - 2*r[i] + e[i])/dz**2 - mub for i in range(n)]
        r = [r[i] - dt*res[i] for i in range(n)]
        if max(abs(x) for x in res) < tol:
            break
    e = [r[0]] + r + [r[-1]]
    pN, pT = [], []
    for i in range(n):
        d1 = (e[i+2] - e[i])/(2*dz)
        lap = (e[i+2] - 2*r[i] + e[i])/dz**2
        iso = p0(r[i], T) - kappa*r[i]*lap
        pN.append(iso + 0.5*kappa*d1*d1)             # Korteweg 압력 텐서의 수직 성분
        pT.append(iso - 0.5*kappa*d1*d1)             # 접선 성분
    sigma_stress = sum(a - b for a, b in zip(pN, pT))*dz
    # 첫 적분 kappa/2 rho'^2 = W(rho)  ->  sigma = int sqrt(2 kappa W) d rho
    m = 4000
    sigma_exact = 0.0
    for k in range(m):
        rr = rg + (k + 0.5)*(rl - rg)/m
        W = max(psi(rr, T) - mub*rr + pb, 0.0)
        sigma_exact += math.sqrt(2*kappa*W)*(rl - rg)/m
    w = (rl - rg)/max(abs(e[i+2] - e[i])/(2*dz) for i in range(n))
    return dict(w=w, it=it, rg=rg, rl=rl, spread=max(abs(p - pb) for p in pN),
                dip=min(pT) - pb, s_st=sigma_stress, s_ex=sigma_exact)
 
print("  T   kappa  width  iters  max|pN-pb|  min(pT-pb)  sigma_stress  sigma_exact  ratio")
for T in (0.85, 0.95):
    for kappa in (0.1, 0.5, 2.0, 8.0):
        d = korteweg_profile(T, kappa)
        print(f"{T:.2f}  {kappa:4.1f}  {d['w']:5.2f}  {d['it']:5d}   {d['spread']:.2e}   {d['dip']:+.4f}"
              f"     {d['s_st']:.5f}       {d['s_ex']:.5f}    {d['s_st']/d['s_ex']:.4f}")

출력은 다음과 같다.

  T   kappa  width  iters  max|pN-pb|  min(pT-pb)  sigma_stress  sigma_exact  ratio
0.85   0.1   2.08     47   2.09e-02   -0.0309     0.10067       0.16545    0.6085
0.85   0.5   2.63    275   1.08e-02   -0.1550     0.32864       0.36995    0.8883
0.85   2.0   4.35  30000   9.00e-03   -0.2429     0.71074       0.73990    0.9606
0.85   8.0   8.21    414   2.52e-03   -0.2650     1.46491       1.47980    0.9899
0.95   0.1   2.39    254   2.67e-03   -0.0112     0.02759       0.03239    0.8518
0.95   0.5   3.99  30000   9.96e-04   -0.0255     0.06887       0.07243    0.9508
0.95   2.0   7.38    290   3.61e-04   -0.0290     0.14301       0.14485    0.9873
0.95   8.0  14.48   1168   9.36e-05   -0.0298     0.28878       0.28971    0.9968

위 코드를 브라우저에서 직접 돌려보자. 기본값은 출력 둘째 줄($T=0.85$, $\kappa=0.5$)이다.

κ 슬라이더를 0.1로 내리면 오른쪽 아래 점이 비율 0.61 근처로 떨어진다. 8로 올리면 0.99에 붙는다. 실행이 끝날 때마다 점이 하나씩 남으니, $T$를 바꿔 가며 출력 표에 없는 조합을 채워 보자.

2칸짜리 계면에서 σ의 39%가 사라진다#

ratio 열이 핵심이다. 계면 폭(밀도 차를 최대 기울기로 나눈 값)이 2.08칸이면 0.6085다. 표면장력의 39%가 격자에서 사라졌다. 폭 2.63칸에서 0.888, 4.35칸에서 0.961, 8.21칸에서 0.990이다. 온도와 무관하게 비율은 폭 하나로 거의 정해진다. $T=0.95$의 폭 3.99칸은 0.951로 $T=0.85$의 4.35칸과 비슷하다.

min(pT-pb) 열도 같은 이야기를 한다. $T=0.85$에서 κ를 키울수록 접선 압력의 골은 -0.031에서 -0.265로 깊어진다. 연속체 값에 가까워지는 것이다. 얇은 계면에서는 골의 바닥이 격자점 사이에 숨어 버린다.

max|pN-pb| 열은 수직 압력이 얼마나 평평하지 않은지 보여준다. 연속체라면 0이다. 폭 2.08칸에서는 $2.1\times10^{-2}$로 접선 골의 3분의 2나 된다. 이 잔여 기울기는 정지 계면에 가짜 힘으로 작용한다. 스퓨리어스 전류(spurious current, 정지해야 할 계면 주위에 생기는 가짜 흐름)의 한 원인이다.

iters가 30000인 두 줄은 수렴 실패가 아니다. 격자 위 계면은 평행 이동에 대해 대칭이 아니다. 그래서 프로파일이 격자점에 맞물리는 위치로 아주 천천히 미끄러진다. 잔차가 $10^{-7}$ 근처에서 멈춘 동안에도 σ는 다섯 자리까지 변하지 않는다.

√κ와 (1−T)^1.5#

sigma_exact 열을 보면 κ가 네 배씩 커질 때 0.165 → 0.370 → 0.740 → 1.480으로 두 배씩 커진다. 첫 적분 공식의 $\sqrt{\kappa}$ 그대로다. 폭도 $\sqrt{\kappa}$에 비례한다. κ를 키우면 계면이 두꺼워지는 대신 표면장력도 함께 커진다는 뜻이다.

온도 의존성은 평균장 임계 지수를 따른다. κ=0.5에서 $T=0.85$와 $T=0.95$의 σ 비는 0.36995/0.07243 = 5.11이다. $(0.15/0.05)^{3/2}=5.20$과 가깝다. 임계점에 가까울수록 두 상의 밀도 차가 줄고 계면은 넓어진다. 그래서 같은 κ에서도 $T=0.95$ 쪽 비율이 더 좋다.

실무에서는 이 둘이 충돌한다. 원하는 σ와 밀도비를 정하면 κ와 $T$가 정해지고, 그러면 폭도 정해진다. 폭이 4칸 아래로 내려가면 σ가 5% 넘게 줄어든다. 보존형 위상장 모델에서도 같은 폭 제약이 나타난다. 보존형 Allen–Cahn의 곡률 항에서 본 계면 두께 W도 같은 이유로 4~5칸 이상을 권한다.

κ를 고르기 전에 한 번 재 볼 폭#

다상 LBM에서 κ를 정했다면 시뮬레이션을 돌리기 전에 평면 계면 하나를 1D로 이완해 보자. 그리고 두 숫자를 확인한다. 하나는 격자 위 $\kappa\sum\rho'^2$와 $\int\sqrt{2\kappa W}\,d\rho$의 비다. 다른 하나는 $\max|p_N-p_b|$다. 비가 0.95 아래라면 Laplace 시험의 σ도 그만큼 작게 나온다. 그 경우 κ를 키우거나 $T$를 임계점 쪽으로 올리는 편이 낫다. 수직 압력이 평평하지 않다면 정지 액적 시험에서 스퓨리어스 전류를 먼저 의심하자.