포아즈이유 유동: 해석해 유도와 유한차분 검증

평판 두 장, 간격 $2h$. 압력은 $x$로만 떨어지고, 완전 발달이면 $v=0$, $u$는 $y$만 남는다. 원관도 같은 이름(Poiseuille flow)으로 불리지만, 아래 식은 2차원 채널이다. 해가 이 포물선에서 벗어나면 점성항이나 no-slip이 틀린 것이다.

−dp/dx와 μ 슬라이더를 움직이면 화면 아래 $u_{\max}$만 바뀌고, 곡선 모양은 그대로다.

그려진 곡선은 $u/u_{\max} = 1 - (y/h)^2$다. Re 슬라이더는 구석의 숫자만 옮긴다. 봉우리 값은 $h=1$일 때 그 슬라이더를 $2\mu$로 나눈 것이다.

문제 설정#

형상 및 경계 조건#

두 무한 평행 평판 사이의 완전 발달 층류 유동을 고려한다.

• 채널 반폭: $h$ (평판은 $y = \pm h$에 위치)
• 유동 방향: $x$ (무한히 길다고 가정)
• 경계 조건: 양 벽면에서 no-slip $\Rightarrow$ $u(\pm h) = 0$
• 구동력: 일정한 압력 구배 $\dfrac{dp}{dx} = -G < 0$ ($G > 0$은 상수)

지배 방정식#

완전 발달 유동이므로 $\partial/\partial x = 0$, $v = 0$이다. $x$-방향 나비에–스토크스 방정식:

$\rho \left( u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} \right) = -\frac{dp}{dx} + \mu \frac{\partial^2 u}{\partial y^2}$

완전 발달 조건($\partial u/\partial x = 0$, $v = 0$)을 적용하면:

$\mu \frac{d^2 u}{dy^2} = \frac{dp}{dx} = -G$

해석해 유도#

1단계: ODE로 단순화#

$\frac{d^2 u}{dy^2} = -\frac{G}{\mu} \equiv -C \quad (C > 0)$

2단계: 1차 적분#

$\frac{du}{dy} = -C y + A_1$

3단계: 2차 적분#

$u(y) = -\frac{C}{2} y^2 + A_1 y + A_2$

4단계: 경계 조건 적용#

$u(+h) = 0$:

$u(-h) = 0$:

(1) $-$ (2): $0 = 2 A_1 h \Rightarrow A_1 = 0$ (대칭)

(1) $+$ (2): $0 = -C h^2 + 2A_2 \Rightarrow A_2 = \dfrac{C h^2}{2}$

최종 해석해#

$\boxed{u(y) = \frac{G}{2\mu}\bigl(h^2 - y^2\bigr)}$

채널 중심($y=0$)에서 최대 속도:

$u_{\max} = \frac{G h^2}{2\mu}$

물리적 해석#

속도 프로파일#

해석해는 포물선(parabola) 분포다. 벽면 마찰이 유동 단면 전체에 균일하게 전달되어, 중심에서 최대, 벽에서 0이 된다.

무차원 속도 $u/u_{\max} = 1 - (y/h)^2$는 채널 반폭 $h$나 점도 $\mu$에 무관한 보편적 프로파일이다.

벽면 전단 응력#

$\tau_w = \mu \left.\frac{du}{dy}\right|_{y=h} = \mu \cdot \frac{G}{\mu}(-h) = -G h$

(크기: $\tau_w = G h$) 압력 구배가 클수록, 채널이 넓을수록 벽 마찰이 커진다.

체적 유량#

$Q = \int_{-h}^{h} u \, dy = \frac{G}{2\mu} \int_{-h}^{h} (h^2 - y^2) \, dy = \frac{G}{2\mu} \cdot \frac{4h^3}{3} = \frac{2 G h^3}{3\mu}$

이를 Hagen-Poiseuille 법칙의 2D 버전이라 부른다. 유량이 $h^3$에 비례하므로, 채널이 두 배 넓어지면 유량은 여덟 배 증가한다.

코드 기본값 $h=G=\mu=1$을 같은 식에 넣으면 $u_{\max}=1/2$, $Q=2/3$이다. 높이 $2h$로 나누면 평균 속도는 $1/3$이다. 봉우리와의 비는 $3/2$라서, 평균은 중심의 $2/3$다. $1-(y/h)^2$를 $-h$에서 $h$까지 적분하고 $2h\,u_{\max}$로 나누면 같은 $2/3$가 나온다. $h$와 $\mu$는 비에서 빠진다.

$h$만 두 배로 키우면 거듭제곱이 갈라진다. $u_{\max}$는 $h^2$이라 4배, $Q$는 $h^3$이라 8배, 벽 응력의 크기 $Gh$는 2배다. 유량이 여덟 배가 됐다고 응력도 그 비율일 것이라 두면, 벽 값을 네 배 크게 읽는다. 관 단면으로 넘어가 이 포물선의 운동에너지를 평균 속도로 나누면 계수가 달라진다. 층류 원관에서 그 계수가 2.0000으로 고정되는 계산은 속도수두 앞의 계수에 있다.

Python으로 해석해 시각화#

import numpy as np
import matplotlib.pyplot as plt
 
# 파라미터
h   = 1.0    # 채널 반폭 [m]
G   = 1.0    # 압력 구배 [Pa/m]
mu  = 1.0    # 동점도 [Pa·s]
 
# 격자
y = np.linspace(-h, h, 400)
u_exact = G / (2 * mu) * (h**2 - y**2)
u_max   = G * h**2 / (2 * mu)
 
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
 
# 속도 프로파일 (가로형)
axes[0].plot(u_exact / u_max, y / h, "k-", lw=2)
axes[0].fill_betweenx(y / h, u_exact / u_max, alpha=0.15)
axes[0].axvline(0, color="gray", lw=0.8, ls="--")
axes[0].set_xlabel(r"$u / u_{\max}$")
axes[0].set_ylabel(r"$y / h$")
axes[0].set_title("Poiseuille 속도 프로파일")
axes[0].set_xlim(-0.05, 1.1)
 
# 전단 응력 분포
tau = mu * np.gradient(u_exact, y)
axes[1].plot(tau / G / h, y / h, "r-", lw=2, label=r"$\tau / \tau_w$")
axes[1].axvline(0, color="gray", lw=0.8, ls="--")
axes[1].set_xlabel(r"$\tau / \tau_w$")
axes[1].set_ylabel(r"$y / h$")
axes[1].set_title("전단 응력 분포 (선형)")
axes[1].legend()
 
plt.tight_layout()
plt.savefig("poiseuille_profile.png", dpi=150)
plt.show()

유한차분법으로 수치 검증#

FD 이산화#

$d^2u/dy^2 = -G/\mu$를 2차 중앙차분으로 이산화한다. 균일 격자 $\Delta y = 2h/(N-1)$:

$\frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta y)^2} = -\frac{G}{\mu}$

경계 조건: $u_0 = u_{N-1} = 0$.

이를 삼대각 행렬(tridiagonal) 연립방정식으로 쓰면:

$\begin{pmatrix} -2 & 1 \\ 1 & -2 & 1 \\ & \ddots & \ddots & \ddots \\ & & 1 & -2 \end{pmatrix} \mathbf{u} = -\frac{G (\Delta y)^2}{\mu} \mathbf{1}$

Python 구현#

import numpy as np
 
def solve_poiseuille_fd(N, h=1.0, G=1.0, mu=1.0):
    """
    N : 내부 격자점 수 (경계 제외)
    returns (y, u_num, u_exact, L2_error)
    """
    dy = 2 * h / (N + 1)
    y_inner = np.linspace(-h + dy, h - dy, N)
 
    # 삼대각 계수 행렬
    diag  = -2 * np.ones(N)
    upper =  np.ones(N - 1)
    lower =  np.ones(N - 1)
    A = (np.diag(diag) + np.diag(upper, k=1) + np.diag(lower, k=-1))
 
    rhs = -G * dy**2 / mu * np.ones(N)
 
    u_num = np.linalg.solve(A, rhs)
 
    # 경계 포함 전체 격자
    y_full  = np.concatenate([[-h], y_inner, [h]])
    u_full  = np.concatenate([[0.0], u_num, [0.0]])
    u_exact = G / (2 * mu) * (h**2 - y_full**2)
 
    L2 = np.sqrt(np.mean((u_full - u_exact)**2))
    return y_full, u_full, u_exact, L2
 
# 격자 수렴 테스트
Ns = [4, 8, 16, 32, 64, 128, 256]
errors = []
for N in Ns:
    _, _, _, err = solve_poiseuille_fd(N)
    errors.append(err)
 
# 수렴 차수
for i in range(1, len(Ns)):
    ratio = errors[i-1] / errors[i]
    order = np.log2(ratio)
    print(f"N={Ns[i]:4d}  L2={errors[i]:.3e}  수렴차수={order:.2f}")

수렴 결과#

“

참고: 포아즈이유 유동은 2차 다항식이므로 2차 중앙차분으로 기계적 정밀도(machine precision) 수준의 정확도를 얻는다. 실제 오차는 $\Delta y^2$ 비례가 아니라 부동소수점 반올림 오차 수준에서 포화된다.

낮은 $N$에서도 포물선은 식에 붙는다. 점성항이 죽었는지는 여기서 바로 보인다. 차수 열은 그 판정에 쓸 수 없다.

이산화 문장은 $\Delta y = 2h/(N-1)$라고 적었다. 표의 간격은 $N=4$에서 0.4000, $N=8$에서 0.2222, $N=16$에서 0.1176, $N=32$에서 0.0606, $N=64$에서 0.0308이다. $h=1$이면 $2/(N+1)$과 같다. $2/5=0.4$, $2/9 \approx 0.2222$, $2/17 \approx 0.1176$, $2/33 \approx 0.0606$, $2/65 \approx 0.0308$. 코드는 양 벽을 미지수에서 빼서 내부 $N$칸을 $2h/(N+1)$로 나눈다. $N-1$은 경계점까지 미지수로 센 개수다.

$L_2$ 다섯 칸은 $2.78\times 10^{-16}$, $1.67\times 10^{-16}$, $2.22\times 10^{-16}$, $2.78\times 10^{-16}$, $2.22\times 10^{-16}$이다. 인접 비의 이진로그는 약 0.7, 음수, 음수, 약 0.3이다. 표의 ~2.0과 다르다. $y$의 2차식은 3계 도함수가 0이라 2차 중앙차분의 절단오차가 0이다. 이 열로 차수를 재지 못한다.

해석 결과 비교 플롯#

import matplotlib.pyplot as plt
 
fig, ax = plt.subplots(figsize=(6, 5))
 
colors = plt.cm.viridis(np.linspace(0.2, 0.9, 4))
Ns_plot = [4, 8, 16, 64]
 
for N, c in zip(Ns_plot, colors):
    y, u_num, _, _ = solve_poiseuille_fd(N)[:3]
    _, _, u_exact, _ = solve_poiseuille_fd(N)
    ax.plot(u_num, y, "o--", color=c, ms=4, label=f"FD N={N}")
 
# 해석해
y_fine = np.linspace(-1, 1, 400)
ax.plot(0.5 * (1 - y_fine**2), y_fine, "k-", lw=2, label="해석해")
 
ax.set_xlabel(r"$u$ (무차원)")
ax.set_ylabel(r"$y/h$")
ax.set_title("FD vs 해석해 비교")
ax.legend(fontsize=9)
plt.tight_layout()
plt.savefig("poiseuille_fd_vs_exact.png", dpi=150)
plt.show()

OpenFOAM 설정 힌트#

포아즈이유 유동을 OpenFOAM icoFoam으로 설정할 때 핵심 파라미터:

0/U — 속도 경계 조건#

boundaryField
{
    walls
    {
        type    noSlip;         // 벽면 no-slip
    }
    inlet
    {
        type    fixedValue;
        value   uniform (1 0 0); // 또는 parabolicVelocity
    }
    outlet
    {
        type    zeroGradient;
    }
}

0/p — 압력 경계 조건#

boundaryField
{
    inlet   { type fixedValue; value uniform 1; }
    outlet  { type fixedValue; value uniform 0; }
    walls   { type zeroGradient; }
}

constant/transportProperties#

nu              nu [ 0 2 -1 0 0 0 0 ] 0.01;  // 동점성계수 조정

발달 길이와 중심 속도#

1. 완전 발달 조건: 채널 길이가 입구 발달 길이($L_{dev} \approx 0.05 \, Re \cdot h$)보다 충분히 길어야 함
2. 격자 독립성: $y$-방향 격자 수 16 이상이면 2차 정확도 확보
3. 검증 지표: 중심선 속도 $u_c = G h^2 / (2\mu)$와 벽면 전단 응력 $\tau_w = G h$ 비교

L2가 1e-16이면 차수 칸을 볼 것인가#

$N=4$에서 $L_2$는 이미 $2.78\times 10^{-16}$이다. 격자를 64까지 넣어도 그 칸은 $10^{-16}$대에서 오르내린다. ~2.0은 그 비에서 나오지 않았다. 2차 다항식 위에서 중앙차분은 절단오차가 0이다. 차수 실험으로 합격 도장을 찍는 판단이 여기서 틀린다. 점성항 $\mu \nabla^2 \mathbf{u}$가 살아 있는지는 중심 속도와 벽 응력으로 본다.

볼 숫자는 $Gh^2/(2\mu)$와 $Gh$다. 오더가 갈라지는 데이터는 불연속이 있는 해 쪽에 있다. 같은 코드가 0.79차와 2.00차를 동시에 낸 경우는 차수 칸이 스킴 차수가 아닐 수 있음을 따로 잰다. 잔차가 반올림에 붙어 있으면 격자를 더 넣어도 차수는 나오지 않는다.

“

참고 문헌

• Batchelor, G. K. An Introduction to Fluid Dynamics. Cambridge, 1967.
• Ferziger, J. H., Perić, M. Computational Methods for Fluid Dynamics. Springer, 2002.
• OpenFOAM Documentation — icoFoam tutorial: cavity