v 항과 v² 항을 더했더니 양 끝은 맞고 가운데서 42% 빗나갔다

선형 저항과 제곱 저항을 더하면 Re≪60과 Re≫60에서는 맞지만, 두 항이 비슷해지는 Re 10~300에서는 후류가 만드는 추가 저항이 빠져 종단속도가 30~40% 커진다.

같은 물, 같은 유리구슬, 다른 저항 법칙#

물속에 유리구슬 두 개를 떨어뜨린다. 하나는 지름 0.05 mm, 다른 하나는 10 mm다. 작은 구슬이 받는 저항은 속도에 비례한다. 큰 구슬이 받는 저항은 속도의 제곱에 비례한다. 물도 같고 재질도 같다. 크기만 다르다.

18세기 말 프랑스 학자들은 이 두 법칙을 따로따로 발견했다. 그리고 둘을 더한 식 $F = a v^2 + b v$ 로 정리했다. 이 글은 그 합이 어디서 맞고 어디서 틀리는지 숫자로 확인한다. 독자는 두 항의 주도권이 바뀌는 Reynolds 수(관성력과 점성력의 비)와, 단순한 합이 가장 크게 빗나가는 구간을 얻어 간다.

다섯 개의 질문으로 나눠 따라간다.

Q1. 1775년부터 1800년까지, 저항은 v였나 v²였나#

1752년 달랑베르는 이상유체에서 물체가 받는 저항이 0이라는 역설에 부딪혔다. 이론이 막히자 실험이 길을 냈다. 1775년 프랑스 아카데미는 유체 저항 실험팀을 꾸렸다. 실무는 Bossut이 맡았다.

Bossut은 모형 배를 수조에서 끌었다. 정밀하게 재보니 저항은 속도의 제곱에 비례했다. 한편 Du Buat는 파이프 안의 흐름을 쟀다. 그가 본 마찰 저항은 속도에 그냥 비례했다. 두 실험은 서로 다른 답을 냈다.

Prony는 더 정밀한 실험으로 둘을 하나로 묶었다. 저항은 $a v^2 + b v$ 로 쓸 수 있다. Coulomb은 이를 두 형태의 저항으로 정리했다. 속도 제곱에 비례하는 것은 형상저항(form drag, 물체 모양이 만드는 압력 차이), 속도에 비례하는 것은 마찰저항(friction drag)이다. 이 구분은 뒤에 Navier와 Stokes의 점성 유동 연구로 이어진다.

누가 맞았을까. 둘 다 맞았다. 모형 배는 크고 빨랐다. 파이프 속 흐름은 느리고 좁았다. 두 실험은 같은 물리의 양 끝을 본 것이다. 그 "끝"을 가르는 숫자가 다음 질문이다.

Q2. 두 항은 어디서 나오는가#

구 하나에 작용하는 저항을 두 항으로 쓰면 이렇다.

F = \underbrace{3\pi\mu d\, v}_{b\,v} + \underbrace{C_\infty \frac{\pi}{8}\rho d^2 v^2}_{a\,v^2}

$\mu$ 는 점성계수, $\rho$ 는 유체 밀도, $d$ 는 구 지름, $C_\infty$ 는 고속에서의 항력계수(구는 약 0.4)다.

첫 항은 Stokes 저항이다. 관성이 전혀 없는 느린 흐름(크리핑 흐름)의 정확해다. 둘째 항은 Newton 저항이다. 앞면에 부딪힌 유체의 운동량이 그대로 압력으로 바뀐다고 보면 나온다.

한 가지 짚을 점이 있다. Stokes 저항 $3\pi\mu d v$ 는 순수한 표면 마찰이 아니다. 그중 2/3가 마찰이고 1/3은 압력 저항이다. 그래서 Coulomb의 "마찰 대 형상" 구분보다 "점성 지배 대 관성 지배"가 더 정확한 이름이다.

아래 시뮬레이션에서 직접 조작해보자. 구슬을 떨어뜨리면 속도가 오르며 두 저항이 자란다. 둘의 합이 무게와 같아지는 점이 종단속도(terminal velocity, 가속이 멈추는 속도)다.

d = 0.30 mm fluid:

v = 0.00 mm/s · v_t = 1000 mm/s · Re = 0.00 · linear share 100.0% · quadratic 0.0%

지름 슬라이더를 0.02 mm에서 20 mm까지 움직여 보자. 작은 구슬에서는 파란 직선(b·v)이 무게선과 만나고, 큰 구슬에서는 주황 포물선(a·v²)이 만난다. 유체를 glycerol로 바꾸면 10 mm 구슬도 Du Buat 쪽으로 넘어간다.

Q3. 주도권은 언제 넘어가는가#

두 항이 같아지는 속도를 구하면 된다.

3\pi\mu d\, v = C_\infty \frac{\pi}{8}\rho d^2 v^2 \quad\Longrightarrow\quad Re \equiv \frac{\rho v d}{\mu} = \frac{24}{C_\infty}

$Re$ 는 구 지름 기준 Reynolds 수다. $C_\infty = 0.4$ 이면 $Re = 60$ 이다.

이 결과가 Bossut과 Du Buat를 화해시킨다. 저항 법칙을 정하는 것은 유체 종류도, 속도 자체도 아니다. $Re$ 하나다. 파이프 속 느린 흐름은 $Re$ 가 작았다. 모형 배는 $Re$ 가 수십만이었다. 앞의 시뮬레이션에서 물과 glycerol을 바꿔 본 것도 결국 $Re$ 를 옮긴 것이다.

같은 모양의 "두 극한 사이를 잇는" 문제는 다른 곳에도 있다. 입자 열전달의 Knudsen 전이영역에서는 연속체와 자유분자 두 극한을 조화평균으로 이었다. Wood 음속에서는 압축률을 부피분율로 더했다. 극한을 어떻게 섞느냐가 결과를 바꾼다. Prony의 선택은 단순 덧셈이었다.

Q4. 그럼 더하기만 하면 되는가#

확인해보자. 물속 유리구슬( $\rho_p = 2500$ kg/m³) 여섯 개를 정지 상태에서 놓는다. 두 항 저항으로 운동방정식을 적분한다.

m \frac{dv}{dt} = (\rho_p - \rho) V g - b v - a v^2

$m$ 은 구의 질량, $V$ 는 부피, $g$ 는 중력가속도다. 부가질량(added mass, 구가 끌고 가는 주변 유체의 관성)은 무시했다.

종단속도를 기준값과 비교한다. 기준은 Schiller–Naumann 상관식 $C_D = \frac{24}{Re}(1 + 0.15 Re^{0.687})$ 이다. $Re < 1000$ 에서 실험과 수 % 안에서 맞는다고 알려진 식이다.

import math
 
# 물속에 가라앉는 유리구슬: Du Buat(선형) 대 Bossut(제곱) 저항
RHO_P, RHO, MU, G = 2500.0, 1000.0, 1.0e-3, 9.81
C_INF = 0.4          # 두 항 법칙의 제곱(형상) 저항 계수
 
def weight(d):
    return (RHO_P - RHO) * math.pi * d**3 / 6 * G
 
def prony_drag(d, v):
    lin = 3 * math.pi * MU * d * v                 # Du Buat / Stokes:  b v
    quad = C_INF * math.pi / 8 * RHO * d**2 * v**2  # Bossut / Newton:  a v^2
    return lin, quad
 
def vt_prony(d):
    a = C_INF * math.pi / 8 * RHO * d**2
    b = 3 * math.pi * MU * d
    return (-b + math.sqrt(b * b + 4 * a * weight(d))) / (2 * a)
 
def cd_schiller(re):
    return 24 / re * (1 + 0.15 * re**0.687) if re < 1000 else 0.44
 
def vt_reference(d):
    lo, hi = 1e-9, 10.0
    for _ in range(200):
        v = 0.5 * (lo + hi)
        re = RHO * v * d / MU
        f = cd_schiller(re) * 0.5 * RHO * v * v * math.pi * d * d / 4
        lo, hi = (v, hi) if f < weight(d) else (lo, v)
    return 0.5 * (lo + hi)
 
def settle(d, steps_per_tau=400):
    m = RHO_P * math.pi * d**3 / 6
    vt = vt_prony(d)
    tau = vt / ((1 - RHO / RHO_P) * G)       # 가속의 시간 척도
    dt = tau / steps_per_tau
    v, t = 0.0, 0.0
    while v < 0.99 * vt:
        lin, quad = prony_drag(d, v)
        v += dt * (weight(d) - lin - quad) / m
        t += dt
    return t, tau, vt
 
print(" d[mm]   Re_t    vt[mm/s]  ref[mm/s]  err%   lin%   t99/tau  Stokes/vt  Newton/vt")
for dmm in [0.05, 0.1, 0.3, 1.0, 3.0, 10.0]:
    d = dmm * 1e-3
    t99, tau, vt = settle(d)
    lin, quad = prony_drag(d, vt)
    ref = vt_reference(d)
    v_st = weight(d) / (3 * math.pi * MU * d)
    v_nw = math.sqrt(weight(d) / (C_INF * math.pi / 8 * RHO * d * d))
    re = RHO * vt * d / MU
    print(f"{dmm:6.2f} {re:7.2f} {vt*1e3:10.2f} {ref*1e3:10.2f} {100*(vt-ref)/ref:6.1f} "
          f"{100*lin/(lin+quad):6.1f} {t99/tau:8.2f} {v_st/vt:10.2f} {v_nw/vt:10.2f}")
re_x = 24 / C_INF  # 3 pi mu d v = C_inf pi/8 rho d^2 v^2
print(f"crossover Re where lin = quad: 24/C_inf = {re_x:.1f}")

 d[mm]   Re_t    vt[mm/s]  ref[mm/s]  err%   lin%   t99/tau  Stokes/vt  Newton/vt
  0.05    0.10       2.04       1.98    2.9   99.8     4.59       1.00      24.27
  0.10    0.81       8.07       7.29   10.6   98.7     4.55       1.01       8.68
  0.30   17.16      57.21      40.21   42.3   77.8     3.92       1.29       2.12
  1.00  193.49     193.49     145.95   32.6   23.7     2.93       4.22       1.14
  3.00 1121.20     373.73     365.75    2.2    5.1     2.70      19.69       1.03
 10.00 6973.63     697.36     667.76    4.4    0.9     2.65     117.23       1.00
crossover Re where lin = quad: 24/C_inf = 60.0

위 코드를 브라우저에서 직접 돌려보자. 같은 구슬 두 개가 나란히 가라앉는다. 하나는 Prony 두 항 저항, 다른 하나는 Schiller–Naumann 기준 저항을 받는다.

d = 0.30 mm C∞ = 0.40 ρp = 2500 kg/m³ steps/frame 8

step 0 · Prony v = 0.00 · ref v = 0.00 mm/s · t99/τ = …
Re_t 17.16 · vt 57.21 · ref 40.21 · err 42.3% · lin 77.8% · Stokes/vt 1.29 · Newton/vt 2.12

기본값 d = 0.30 mm에서 err 42.3%, lin 77.8%, t99/τ 3.92가 출력 블록과 같게 나온다. 오른쪽 곡선은 지름 전체에 대한 오차다. C∞를 0.6으로 올리면 큰 구슬 쪽 오차가 음수로 돌아서고, ρp를 낮추면 같은 지름에서도 Re_t가 내려가 봉우리가 오른쪽으로 밀린다.

출력에서 세 가지가 읽힌다. 첫째, 양 끝은 잘 맞는다. 0.05 mm에서 2.9%, 3 mm에서 2.2%다. 둘째, 한쪽 항만 쓰면 반대편 끝에서 무너진다. Stokes만 쓰면 10 mm 구슬이 117배 빨라진다. Newton만 쓰면 0.05 mm 구슬이 24배 빨라진다. 셋째, 가운데가 문제다. 0.3 mm( $Re_t = 17$ )에서 42%, 1 mm( $Re_t = 193$ )에서 33% 빗나간다.

가속 시간도 법칙을 따라 바뀐다. 선형 저항만 있으면 속도는 지수함수로 다가가므로 $t_{99}/\tau = \ln 100 = 4.61$ 이다. 제곱 저항만 있으면 $\tanh$ 로 다가가므로 $t_{99}/\tau = \operatorname{artanh} 0.99 = 2.65$ 다. 출력의 4.59와 2.65가 두 극한이고, 가운데 지름은 그 사이를 지난다.

Q5. 왜 가운데서 빗나가는가#

두 항 법칙을 항력계수로 고쳐 쓰면 원인이 보인다.

C_D^{\text{Prony}} = \frac{24}{Re} + C_\infty, \qquad C_D^{\text{S–N}} = \frac{24}{Re}\left(1 + 0.15\, Re^{0.687}\right)

두 식 모두 $Re \to 0$ 에서 $24/Re$ 로 간다. 차이는 괄호 안의 둘째 항이다.

$Re = 17$ 을 넣어 보자. Prony는 $1.41 + 0.40 = 1.81$ 을 준다. Schiller–Naumann은 $1.41 \times 2.05 = 2.90$ 을 준다. 기준보다 저항이 38% 작다. 무게는 같으므로 Prony 구슬은 더 빨리 떨어져야 균형이 맞는다. 그 결과가 42% 큰 종단속도다.

모자란 저항은 어디서 오는가. $Re$ 가 1을 넘으면 앞뒤 대칭이던 흐름이 깨진다. 뒤쪽에 느린 후류가 생기고, $Re \approx 20$ 부근에서 표면 흐름이 떨어져 나가 재순환 영역이 생긴다. 이 후류 압력 저항은 점성 항과 관성 항이 섞여서 생긴다. $v$ 에도 $v^2$ 에도 깔끔하게 속하지 않는다. Prony의 덧셈에는 이 섞인 항이 들어갈 자리가 없다.

그래서 두 항이 비슷한 크기일 때, 바로 $Re \approx 60$ 근처에서 덧셈이 가장 크게 틀린다. 역사적 두 실험이 각자 맞았던 이유도 같다. 둘 다 이 구간 밖에 있었다.

두 항 법칙이 지금도 살아 있는 자리#

틀리는 구간이 있는데도 $a v^2 + b v$ 꼴은 사라지지 않았다. 다공성 매질의 Ergun 식과 Darcy–Forchheimer 식이 정확히 이 모양이다. 파이프 손실의 Prony 식도 같은 뿌리다. 다공성 매질에서는 공극 크기 기준 $Re$ 가 넓은 구간에 걸쳐 있어도 두 항 합이 꽤 잘 맞는다. 계수를 실험에 맞춰 넣기 때문이다.

Euler–Lagrange 입자 코드에서는 다르다. 입자 하나의 $Re$ 가 바로 이 문제 구간(1~~1000)에 들어오는 경우가 많다. 분무 액적, 유동층 입자, 침전 입자가 그렇다. 여기에 $24/Re + 0.4$ 를 쓰면 종단속도가 30~~40% 틀린다. 침전 시간이나 체류 시간이 그만큼 어긋난다.

코드의 항력 모델을 고를 때 이 순서로 확인하면 된다.

입자 $Re$ 분포의 범위를 먼저 뽑는다. 1 미만이면 Stokes, 1000 초과면 상수 $C_D$ 로 충분하다.
1~1000 구간에 질량이 많이 걸려 있으면 Schiller–Naumann처럼 중간 항이 있는 상관식을 쓴다.
단일 구슬 침강 테스트를 하나 돌려 위 표의 0.3 mm와 1 mm 행과 비교한다. 오차 30%가 나오면 두 항 덧셈이 숨어 있다는 신호다.

Bossut의 배와 Du Buat의 파이프는 같은 물을 다른 $Re$ 에서 본 것이었다. 둘 사이의 빈칸은 덧셈만으로 메워지지 않는다.

v 항과 v² 항을 더했더니 양 끝은 맞고 가운데서 42% 빗나갔다 — Prony 두 항 저항법칙과 Re=60