ML 연구자를 위한 연속체 역학과 시뮬레이션

지배 방정식

왜 필요한가. 연속체 역학의 방정식은 종류가 많아 보이지만, 대부분은 같은 문장 하나로 요약된다. "어떤 양이 한 영역 안에서 늘어나는 속도는, 경계를 통해 들어오는 양과 안에서 만들어지는 양의 합이다." 이 문장, 즉 보존 법칙(conservation law)을 알아보면 처음 보는 방정식도 무엇을 보존하고 무엇이 그것을 옮기는지 읽을 수 있다. 학습된 모델이 깨뜨리면 안 되는 성질도 여기서 나온다.

정의. 밀도 $q(\mathbf{x}, t)$(단위 부피당 양), 유속(flux) $\mathbf{F}$(단위 면적·단위 시간당 경계를 지나는 양), 생성항 $s$에 대해 고정된 영역 $V$에서

가 성립한다. $\mathbf{n}$은 경계의 바깥 방향 단위 법선이다. 해가 매끄러우면 발산 정리로 미분형을 얻는다.

$q$와 $\mathbf{F}$에 무엇을 넣느냐에 따라 익숙한 방정식이 된다.

표에서 $\rho$는 밀도, $\mathbf{u}$는 속도, $p$는 압력, $c_p$는 비열, $T$는 온도, $k$는 열전도도, $\boldsymbol{\sigma}$는 응력 텐서다. 유속은 두 종류로 나뉜다. 흐름이 양을 실어 나르는 이류(advection, $q\mathbf{u}$)와, 기울기를 따라 퍼뜨리는 확산(diffusion, $-k\nabla T$ 같은 항)이다. 방정식이 어떤 성질을 갖는지는 둘 중 무엇이 지배적인지에 달려 있고, 이것이 다음 절의 무차원수다. 방정식 하나를 닫으려면 응력과 변형률, 열유속과 온도 기울기를 잇는 구성 관계(constitutive relation)가 필요하다. Newton 유체, Hooke 탄성, Fourier 법칙이 그런 관계이고, 이것은 보존 법칙과 달리 물질에 대한 모델이다.

직관: 왜 보존형이 중요한가. 미분형 $\partial_t q + \nabla\cdot\mathbf{F} = s$는 해가 매끄러울 때만 적분형과 같다. 충격파나 접촉 불연속처럼 해가 끊기는 곳에서는 적분형만 의미가 있고, 불연속이 움직이는 속도는 적분형에서 나온다(Rankine–Hugoniot 조건). 1차원 Burgers 방정식 $u_t + (u^2/2)_x = 0$은 매끄러운 곳에서 $u_t + u\,u_x = 0$과 같지만, 두 식을 각각 그대로 이산화하면 결과가 다르다. 왼쪽 $u = 1$, 오른쪽 $u = 0$인 계단에서 충격파는 속도 $(1 + 0)/2 = 0.5$로 움직여야 한다.

import numpy as np

N, L_m = 200, 2.0
dx_m = L_m / N
x_m = (np.arange(N) + 0.5) * dx_m                  # (N,) 셀 중심
u0 = np.where(x_m < 0.5, 1.0, 0.0)                 # (N,) 왼쪽 1 m/s, 오른쪽 0 m/s
dt_s, steps = 0.4 * dx_m, 200                      # CFL 0.4, t = 0.8 s

def step_conservative(u):                          # u_t + (u²/2)_x = 0, upwind 유속
    F = 0.5 * u**2                                 # (N,) 셀 j의 오른쪽 면 유속 (u >= 0)
    return u - dt_s / dx_m * (F - np.concatenate([[0.5], F[:-1]]))   # 왼쪽 유입 u = 1

def step_nonconservative(u):                       # u_t + u u_x = 0, 같은 upwind
    du = u - np.concatenate([[1.0], u[:-1]])
    return u - dt_s / dx_m * u * du

for name, step in [("보존형", step_conservative), ("비보존형", step_nonconservative)]:
    u = u0.copy()
    for _ in range(steps):
        u = step(u)
    front = x_m[np.argmax(u < 0.5)]                # u가 처음 0.5 아래로 내려가는 위치
    print(f"{name:5s} 충격파 위치 {front:.2f} m (정확한 값 {0.5 + 0.5 * dt_s * steps:.2f} m)")

보존형   충격파 위치 0.91 m (정확한 값 0.90 m)
비보존형  충격파 위치 0.51 m (정확한 값 0.90 m)

보존형은 셀 경계의 유속을 주고받으므로 한 셀이 잃은 양을 이웃 셀이 정확히 얻는다. 그래서 충격파가 이론 위치(0.90 m) 근처에 있다. 비보존형은 같은 upwind 차분인데 충격파가 출발점에 멈춰 있다. 충격파 바로 앞 셀의 $u$가 0이라 $u\,u_x$ 항이 영원히 0이기 때문이다. 격자를 줄여도 이 오류는 사라지지 않는다. 보존형 방법이 수렴한다면 올바른 약해(weak solution)로 수렴한다는 결과가 Lax–Wendroff 정리이고, 비보존형 방법에는 그런 보장이 없다(LeVeque, 2002).

연구에서의 의미. 보존 법칙은 학습된 모델을 검사하는 가장 값싼 도구다. 정답 해가 없어도 예측된 장(field)의 전체 질량·에너지가 경계 유속과 맞는지는 계산할 수 있다(26.5장). 또 모델 구조를 정할 때 셀 사이 유속을 예측하게 하면 보존이 구조로 보장되고, 셀 값을 직접 예측하면 보장되지 않는다. 같은 정확도라면 앞의 것이 장기 예측에서 덜 표류한다.

무차원수와 영역

왜 필요한가. 같은 Navier–Stokes 방정식이 꿀이 흐르는 문제와 비행기 날개 주위의 공기를 모두 기술한다. 둘의 행동이 완전히 다른 것은 방정식이 아니라 항들의 상대적 크기가 다르기 때문이다. 이 상대적 크기를 재는 것이 무차원수다. 무차원수를 모르면 학습 데이터가 어떤 영역을 덮는지, 새 조건이 내삽인지 외삽인지 판단할 수 없다.

정의. 대표 속도 $U$, 대표 길이 $L$, 동점성계수 $\nu$(m²/s), 확산계수 $D$(m²/s)로

를 정의한다. Reynolds 수 Re가 작으면 점성이 교란을 빨리 감쇠시켜 흐름이 매끄럽고(층류, laminar), 크면 작은 교란이 자라 불규칙한 흐름(난류, turbulent)이 된다. 원관 유동에서는 Reynolds가 1883년 실험에서 이 전이를 하나의 무차원 조합으로 기술했다(Reynolds, 1883). 흔히 Re ≈ 2300 근처를 원관 유동의 전이 기준으로 쓰지만, 입구 조건과 교란의 크기에 따라 실제 전이 지점은 크게 달라진다. 그 밖에도 열전달의 Prandtl 수($\nu/\alpha$, 운동량 확산 대 열 확산), Mach 수(유속 대 음속), 고체 열전달의 Biot 수(내부 전도 저항 대 표면 열전달 저항) 등이 문제의 영역을 정한다.

직관: 무차원수가 같으면 해가 같다. 길이를 $L$, 속도를 $U$, 시간을 $L/U$로 나누어 방정식을 무차원화하면, 점성 Burgers 방정식 $u_t + u u_x = \nu u_{xx}$는

가 된다. 별표는 무차원량이다. 남는 파라미터는 Re 하나다. 따라서 $U, L, \nu$를 바꾸어도 Re가 같으면 무차원 해는 똑같다. 축소 모형으로 하는 풍동 실험이 이 원리(상사 법칙, similarity)에 기댄다.

def burgers(U_ms, L_m, nu_m2s, N=128, t_star=0.3):
    """u_t + u u_x = nu u_xx, 주기 경계, u0 = U sin(2 pi x / L)를 t = t_star L/U까지 푼다."""
    dx = L_m / N
    u = U_ms * np.sin(2 * np.pi * np.arange(N) * dx / L_m)    # (N,)
    T_s = t_star * L_m / U_ms
    n = int(np.ceil(T_s / (0.2 * min(dx / U_ms, dx**2 / nu_m2s))))
    dt = T_s / n                                               # explicit 안정 한계의 0.2배 이하
    for _ in range(n):
        up, um = np.roll(u, -1), np.roll(u, 1)
        u = u - dt * u * (up - um) / (2 * dx) + dt * nu_m2s * (up - 2 * u + um) / dx**2
    return u, dx

for U, L, nu in [(1.0, 1.0, 0.01), (4.0, 0.5, 0.02), (1.0, 1.0, 0.04)]:
    u, dx = burgers(U, L, nu)
    grad = np.max(np.abs(np.roll(u, -1) - np.roll(u, 1)) / (2 * dx)) * L / U   # 무차원 최대 기울기
    print(f"U={U} m/s L={L} m nu={nu} m²/s  Re={U * L / nu:5.1f}  무차원 최대 기울기 {grad:.3f}")

U=1.0 m/s L=1.0 m nu=0.01 m²/s  Re=100.0  무차원 최대 기울기 39.640
U=4.0 m/s L=0.5 m nu=0.02 m²/s  Re=100.0  무차원 최대 기울기 39.640
U=1.0 m/s L=1.0 m nu=0.04 m²/s  Re= 25.0  무차원 최대 기울기 7.557

첫째와 둘째 조건은 속도·길이·점성이 모두 다르지만 Re = 100으로 같아서 무차원 결과가 같다. 셋째 조건은 Re = 25로, 점성이 상대적으로 커서 같은 무차원 시각의 기울기가 다섯 배 이상 작다.

연구에서의 의미. 대리모델의 입력을 $U, L, \nu$ 세 개로 주면 모델은 세 변수가 실제로는 Re 하나를 통해서만 작용한다는 사실을 데이터에서 다시 배워야 한다. 무차원수를 입력으로 쓰면 입력 차원이 줄고, 학습 범위 밖의 단위 조합에도 일반화될 가능성이 커진다. 다만 문제에 여러 무차원수가 관여하면(예: 열유동의 Re와 Prandtl 수, 압축성 유동의 Re와 Mach 수) 그중 하나만 맞추는 것으로는 상사가 성립하지 않는다. 학습 데이터의 범위는 원래 단위가 아니라 무차원수 공간에서 기술하는 것이 좋다.

시뮬레이션 파이프라인

왜 필요한가. "시뮬레이션 결과"는 한 번의 계산이 아니라 여러 단계의 선택이 쌓인 결과다. 각 단계의 선택이 데이터에 흔적을 남기고, 단계마다 비용 구조가 다르다. 대리모델이 무엇을 대신하는지 말하려면 이 단계를 알아야 한다.

단계.

1. 형상(geometry): CAD 모델이나 매개변수화된 형상. 설계 최적화(26.4장)에서는 이것이 설계 변수다.
2. 메시(mesh): 영역을 셀이나 요소로 나눈다. 구조 격자와 비구조 격자가 있고, 벽 근처의 경계층이나 응력 집중부처럼 해가 급하게 변하는 곳에 셀을 몰아 둔다. 메시 생성은 사람이 개입하는 경우가 많아 자동화가 가장 어려운 단계다.
3. 경계조건과 초기조건: 유입 속도 분포, 벽면 조건, 온도나 열유속 조건, 고정·하중 조건. 같은 방정식과 메시라도 경계조건이 바뀌면 다른 문제다.
4. solver: 유한차분·유한체적·유한요소 이산화, 시간 적분, 선형·비선형 solver와 그 tolerance(9장).
5. 후처리(post-processing): 장에서 관심 출력량(quantity of interest, QoI)을 계산한다. 항력·양력 계수, 최대 응력, 평균 Nusselt 수, 스펙트럼 등.

비용의 구조. 3차원에서 한 방향 격자 수를 $N$이라 하면 셀 수는 $N^3$이다. explicit 시간 적분은 CFL 조건 때문에 $\Delta t \propto \Delta x$이므로 step 수가 $N$에 비례하고, 전체 비용은 대략 $N^4$이다. 확산이 지배적이면 $\Delta t \propto \Delta x^2$이라 더 빨리 늘어난다. implicit 방법은 step 수를 줄이는 대신 step마다 큰 선형계를 풀어야 한다. 그래서 해상도를 두 배로 올리는 비용은 두 배가 아니라 수십 배다. 난류처럼 작은 규모까지 해상해야 하는 문제에서는 필요한 해상도 자체가 Re와 함께 늘어나므로(26.3장), 고해상도 시뮬레이션은 소수의 조건에서만 할 수 있다. 대리모델이 필요한 이유이면서, 학습 데이터가 적은 이유이기도 하다.

시뮬레이션 데이터

왜 필요한가. 9.7장은 시뮬레이터를 데이터 생성기로 쓸 때의 원칙(이산화 오차는 체계적이고, 해상도가 다른 데이터는 같은 대상이 아니다)을 다뤘다. 여기서는 그 원칙을 물리 문제에 적용할 때 자주 걸리는 두 가지를 본다. 해상도가 부족하면 오차가 작아지는 대신 비물리적인 해가 나올 수 있다는 것, 그리고 같은 시뮬레이션에서도 출력량마다 수렴 속도가 다르다는 것이다.

격자 Péclet 수. 정상 이류-확산 방정식 $U c_x = D c_{xx}$($c(0) = 0$, $c(1) = 1$)을 2차 중심차분으로 풀면, 격자 간격으로 정의한 격자 Péclet 수 $\mathrm{Pe}_\Delta = U\Delta x/D$가 2를 넘을 때 해가 진동한다(Patankar, 1980). 농도는 0과 1 사이에 있어야 하는데, 음수 농도가 나온다.

Pe = 50.0                                          # U L / D, 정상 이류-확산 U c_x = D c_xx, c(0)=0, c(1)=1
for N in [10, 20, 40, 80]:
    Pc = Pe / N                                    # 격자 Péclet 수 U dx / D
    A = (np.diag(np.full(N - 1, 2.0))
         + np.diag(np.full(N - 2, -(1 + Pc / 2)), -1)
         + np.diag(np.full(N - 2, -(1 - Pc / 2)), 1))   # (N-1, N-1) 중심차분
    b = np.zeros(N - 1)
    b[-1] = 1 - Pc / 2                              # 오른쪽 경계값 c = 1
    c = np.linalg.solve(A, b)                       # (N-1,) 내부 격자점의 농도
    neg = f"있음 (최솟값 {c.min():.3f})" if c.min() < -1e-9 else "없음"
    print(f"N={N:3d}  격자 Péclet {Pc:5.2f}  음수 농도 {neg}")

N= 10  격자 Péclet  5.00  음수 농도 있음 (최솟값 -0.429)
N= 20  격자 Péclet  2.50  음수 농도 있음 (최솟값 -0.111)
N= 40  격자 Péclet  1.25  음수 농도 없음
N= 80  격자 Péclet  0.62  음수 농도 없음

전체 Péclet 수는 50으로 같다. 격자를 바꾼 것뿐인데 $\mathrm{Pe}_\Delta = 5, 2.5$에서는 물리적으로 불가능한 음수 농도가 나오고, 2 아래로 내려가면 사라진다. 이런 해도 solver는 "수렴했다"고 보고한다. 선형계는 정확히 풀렸기 때문이다. 수렴 판정은 이산화된 방정식을 잘 풀었는지를 말할 뿐, 이산화가 물리를 담고 있는지는 말하지 않는다. upwind 차분으로 바꾸면 진동은 사라지지만 대신 수치 확산이 해를 뭉갠다(9.7장). 같은 문제가 26.2–26.3장의 Burgers 데이터에서는 격자 Reynolds 수 $U\Delta x/\nu$로 나타난다.

출력량마다 다른 수렴. 점성 Burgers 방정식(Re = 100)을 격자 다섯 개로 풀고, 두 출력량을 비교한다. 국소량인 최대 기울기 $\max|u_x|$와, 영역 전체의 적분량인 운동에너지 소산율 $\varepsilon = \nu\int u_x^2\,dx$다. 가장 고운 세 격자로 관측 차수와 Richardson 오차 추정을 구하는 방법은 9.7장 그대로다.

rows = []
for N in [32, 64, 128, 256, 512]:
    u, dx = burgers(1.0, 1.0, 0.01, N=N)
    ux = (np.roll(u, -1) - np.roll(u, 1)) / (2 * dx)   # (N,)
    rows.append((N, np.max(np.abs(ux)), 0.01 * np.sum(ux**2) * dx))
    print(f"N={N:3d}  격자 Re {100 / N:5.2f}  최대 기울기 {rows[-1][1]:7.3f}  소산율 {rows[-1][2]:.4f}")
(_, g1, e1), (_, g2, e2), (_, g3, e3) = rows[-3:]
for name, f1, f2, f3 in [("최대 기울기", g1, g2, g3), ("소산율", e1, e2, e3)]:
    p = np.log2(abs(f2 - f1) / abs(f3 - f2))
    print(f"{name}: 관측 차수 {p:.2f}, N=512의 Richardson 오차 추정 {abs(f3 - f2) / (2**p - 1):.1e}")

N= 32  격자 Re  3.12  최대 기울기  37.518  소산율 0.6033
N= 64  격자 Re  1.56  최대 기울기  39.889  소산율 0.5189
N=128  격자 Re  0.78  최대 기울기  39.640  소산율 0.4962
N=256  격자 Re  0.39  최대 기울기  39.475  소산율 0.4903
N=512  격자 Re  0.20  최대 기울기  39.434  소산율 0.4889
최대 기울기: 관측 차수 2.01, N=512의 Richardson 오차 추정 1.3e-02
소산율: 관측 차수 2.07, N=512의 Richardson 오차 추정 4.4e-04

두 출력량 모두 가장 고운 세 격자에서는 2차로 수렴한다. 그러나 거친 격자에서 행동이 다르다. 소산율은 격자를 줄일 때마다 한 방향으로 줄어들어, $N = 64$에서도 오차의 방향과 크기를 짐작할 수 있다. 최대 기울기는 $N = 64$에서 오히려 수렴값보다 커졌다가 다시 내려온다. 격자 Re가 2 근처인 거친 격자에서는 국소량이 단조롭게 수렴하지 않아, 두 격자만 비교하면 오차를 잘못 추정한다. 학습 데이터의 "정답"으로 쓰려는 출력량마다 수렴을 따로 확인해야 하는 이유다.

데이터에 기록할 것. 시뮬레이션 데이터셋의 각 항목에는 결과뿐 아니라 그것을 만든 조건을 함께 저장한다.

• 지배 방정식과 구성 관계, 무차원수와 그 범위
• 형상, 경계조건, 초기조건(무작위로 만들었다면 시드와 분포)
• 메시(격자 수, 최소 셀 크기, 격자 Re·Péclet 수), 시간 간격 또는 tolerance
• solver와 이산화 방식, 코드 이름과 버전
• 수렴 검사 결과: 대표 조건에서의 격자 수렴 연구, 각 출력량의 추정 이산화 오차
• 실패하거나 제외한 실행의 목록과 제외 기준

저장 형식은 배열 데이터를 압축해 담는 HDF5·NetCDF·Zarr 같은 형식이 흔히 쓰인다(집필 시점 예시). 형식보다 중요한 것은 메타데이터가 데이터와 함께 다닌다는 것이다. 파일 이름에 조건을 적어 두는 방식은 데이터가 복사되고 합쳐지는 순간 깨진다.

참고 문헌

• Reynolds, O. (1883). An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels. Philosophical Transactions of the Royal Society of London, 174, 935–982. doi:10.1098/rstl.1883.0029 — 층류·난류 전이를 하나의 무차원 조합으로 기술한 원 논문
• LeVeque, R. J. (2002). Finite Volume Methods for Hyperbolic Problems. Cambridge University Press. doi:10.1017/CBO9780511791253 — 보존 법칙, 보존형 유한체적법, Lax–Wendroff 정리와 비보존형 방법의 실패
• Patankar, S. V. (1980). Numerical Heat Transfer and Fluid Flow. Hemisphere. — 이류-확산 이산화와 격자 Péclet 수
• Oberkampf, W. L., Roy, C. J. (2010). Verification and Validation in Scientific Computing. Cambridge University Press. ISBN 978-0-521-11360-1. — 검증·확인 절차와 이산화 오차 추정