[OpenFOAM 튜토리얼] 후방 계단 유동 해석 — simpleFoam으로 재순환 영역 포착하기

하류 1.5 m는 격자 H로 15H다. 주석의 H=0.01 m로 읽으면 같은 구간이 150H가 된다.

후방 계단 바로 뒤에서 유속이 음이 된다. 그 구간이 끝나는 x가 재부착점이다. simpleFoam이 그 점을 몇 H에 찍는지는, 격자의 높이와 nu 주석의 높이가 같을 때만 읽을 수 있다.

Re 슬라이더를 50에서 1500으로 밀면 스텝 뒤 역류 구간이 길어진다. 입자 슬라이더(50–500)는 그 구간을 채우는 궤적 수만 바꾼다.

Re 400입자 220

스텝 직후 재순환 영역(Lr)이 형성된다. Re가 커지면 Lr가 길어지다가 Re ≳ 1000부근부터 비정상이 된다 — Armaly et al. 1983 실험과 일치.

화면 아래 숫자 Lr이 스텝 몇 개 분량으로 늘어나는지만 보면 된다. Re가 낮은 쪽에 있으면 거품이 스텝에 붙어 있고, 슬라이더 끝으로 가면 하류로 늘어진다.

선행 조건#

OpenFOAM v2312 이상. .com이든 .org든 이 케이스 구성은 같다.
터미널과 케이스 세 칸(0, constant, system)은 이미 익숙한 것으로 둔다.
그림은 ParaView 5.x에서 연다.

1단계: 케이스 디렉토리 생성#

빈 디렉터리에 0, constant, system을 만든다.

# OpenFOAM 환경 로드 (설치 경로에 따라 다름)
source /opt/openfoam2312/etc/bashrc
 
# 튜토리얼 디렉토리 복사 또는 직접 생성
mkdir -p backwardFacingStep
cd backwardFacingStep

트리는 이렇게 갈라져야 한다.

tree backwardFacingStep/

backwardFacingStep/
├── 0/
│   ├── U          # 속도 초기/경계 조건
│   ├── p          # 압력 초기/경계 조건
│   ├── k          # 난류 운동에너지
│   └── epsilon    # 난류 소산율
├── constant/
│   ├── transportProperties
│   └── turbulenceProperties
└── system/
    ├── blockMeshDict
    ├── fvSchemes
    └── fvSolution

2단계: 격자 생성 — blockMeshDict#

system/blockMeshDict가 계단을 블록 두 개로 자른다. 꼭짓점의 y 간격은 1이고 scale은 0.1이라, 주석의 H = 0.1 m는 이 격자와 맞다.

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2312                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
scale   0.1;   // 전체 스케일 (m 단위 → 0.1 곱하면 최종 크기)
 
vertices
(
    // 입구 채널 블록 (계단 위쪽)
    (0  1  0)  // 0
    (5  1  0)  // 1
    (5  2  0)  // 2
    (0  2  0)  // 3
    (0  1  1)  // 4
    (5  1  1)  // 5
    (5  2  1)  // 6
    (0  2  1)  // 7
 
    // 계단 하류 블록 (계단 전체 높이)
    (5  0  0)  // 8
    (20 0  0)  // 9
    (20 2  0)  // 10
    (5  2  0)  // 11 (= vertex 2)
    (5  0  1)  // 12
    (20 0  1)  // 13
    (20 2  1)  // 14
    (5  2  1)  // 15 (= vertex 6)
);
 
blocks
(
    // 입구 채널
    hex (0 1 2 3 4 5 6 7)   (50 10 1) simpleGrading (1 1 1)
    // 계단 하류 (더 조밀한 격자)
    hex (8 9 10 11 12 13 14 15) (150 20 1) simpleGrading (1 1 1)
);
 
boundary
(
    inlet
    {
        type patch;
        faces ((0 3 7 4));
    }
    outlet
    {
        type patch;
        faces ((9 10 14 13));
    }
    wall
    {
        type wall;
        faces
        (
            (0 1 5 4)   // 상단 입구 벽
            (1 8 12 5)  // 계단 수직면
            (8 9 13 12) // 하단 바닥
            (3 2 6 7)   // 상단 벽 (입구)
            (11 10 14 15) // 상단 벽 (하류)
        );
    }
    frontAndBack
    {
        type empty;  // 2D 해석
        faces
        (
            (0 1 2 3)
            (4 5 6 7)
            (8 9 10 11)
            (12 13 14 15)
        );
    }
);

# 격자 생성 실행
blockMesh

예상 출력:

Create time
 
Creating block mesh from "system/blockMeshDict"
...
Writing polyMesh with 3600 cells
End

3단계: 경계 조건 설정#

0/U — 속도 경계 조건#

입구는 x 방향 1 m/s, 벽은 noSlip, 출구는 zeroGradient다.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
dimensions      [0 1 -1 0 0 0 0];  // m/s
 
internalField   uniform (0 0 0);   // 초기값: 정지 상태
 
boundaryField
{
    inlet
    {
        type        fixedValue;
        value       uniform (1 0 0);  // x 방향 1 m/s 균일 유입
    }
    outlet
    {
        type        zeroGradient;  // 출구: 자유 유출
    }
    wall
    {
        type        noSlip;        // 벽면: 점착 조건
    }
    frontAndBack
    {
        type        empty;
    }
}

0/p — 압력 경계 조건#

압력은 운동압력이다. 출구만 0으로 고정한다.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
 
dimensions      [0 2 -2 0 0 0 0];  // 운동 압력 (p/rho), m²/s²
 
internalField   uniform 0;
 
boundaryField
{
    inlet
    {
        type        zeroGradient;
    }
    outlet
    {
        type        fixedValue;
        value       uniform 0;  // 출구 기준 압력 0
    }
    wall
    {
        type        zeroGradient;
    }
    frontAndBack
    {
        type        empty;
    }
}

0/k — 난류 운동에너지#

입구 난류강도 5%면 $k = 1.5 \times (0.05)^2 = 0.00375$ 다. 파일에 적힌 값과 같다.

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      k;
}
 
dimensions      [0 2 -2 0 0 0 0];  // m²/s²
 
// 입구 난류강도 5%, 유속 1 m/s → k = 1.5*(I*U)² = 1.5*(0.05*1)² ≈ 0.00375
internalField   uniform 0.00375;
 
boundaryField
{
    inlet
    {
        type        fixedValue;
        value       uniform 0.00375;
    }
    outlet
    {
        type        zeroGradient;
    }
    wall
    {
        type        kqRWallFunction;
        value       uniform 0.00375;
    }
    frontAndBack
    {
        type        empty;
    }
}

4단계: 물성치 및 난류 모델 설정#

constant/transportProperties#

ν는 이 파일에 있다. Re 1000을 만드는 H가 격자 높이와 같은지는 결과 표 뒤에서 맞춘다.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      transportProperties;
}
 
transportModel  Newtonian;
nu              1e-5;   // 운동점성계수 (m²/s), Re = U*H/nu = 1*0.01/1e-5 = 1000

constant/turbulenceProperties#

RAS에 kEpsilon을 켜 둔다.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      turbulenceProperties;
}
 
simulationType  RAS;
 
RAS
{
    RASModel        kEpsilon;   // k-epsilon 두 방정식 모델
    turbulence      on;
    printCoeffs     on;
}

5단계: 수치 기법 설정#

system/fvSchemes#

시간항은 steadyState, 속도 대류는 linearUpwind, k와 epsilon은 upwind다.

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
 
ddtSchemes      { default steadyState; }  // 정상 상태 해석
 
gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;  // 경사도 제한자
}
 
divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);  // 2차 정확도
    div(phi,k)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}
 
laplacianSchemes { default Gauss linear corrected; }
 
interpolationSchemes { default linear; }
 
snGradSchemes { default corrected; }

system/fvSolution#

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
 
solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-7;
        relTol          0.01;
    }
    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
    }
    "(k|epsilon)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0.1;
    }
}
 
SIMPLE
{
    nNonOrthogonalCorrectors    0;
    pRefCell                    0;
    pRefValue                   0;
 
    residualControl
    {
        p           1e-4;
        U           1e-4;
        "(k|epsilon)" 1e-4;
    }
}
 
relaxationFactors
{
    fields      { p 0.3; }           // 압력 완화계수 (발산 방지)
    equations
    {
        U       0.7;
        k       0.7;
        epsilon 0.7;
    }
}

6단계: 해석 실행#

압력 완화 0.3과 U·k·epsilon의 0.7이, 로그를 500회에서 잔차 1e-4 아래로 내리는 결합이다. 그 압력 보간이 격자 진동을 어떻게 감추는지는 SIMPLE의 Rhie–Chow가 따로 적는다.

# 격자 품질 확인 (선택 사항, 권장)
checkMesh
 
# simpleFoam 실행
simpleFoam 2>&1 | tee log.simpleFoam

수렴 로그:

Time = 1
 
smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 0.0234, No Iterations 3
smoothSolver:  Solving for Uy, Initial residual = 0.412, Final residual = 0.00891, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.987, Final residual = 0.00821, No Iterations 15
...
 
Time = 500
 
smoothSolver:  Solving for Ux, Initial residual = 1.23e-05, ...
GAMG:  Solving for p, Initial residual = 8.91e-05, ...
 
SIMPLE solution converged in 500 iterations

멈추는 기준은 모든 변수의 잔차가 1e-4 아래다. 이 로그는 500번째에 그 선을 넘는다.

7단계: 결과 시각화 (ParaView)#

paraFoam으로 연다.

# ParaView용 더미 파일 생성 후 실행
touch backwardFacingStep.foam
paraFoam

ParaView에서 순서는 이렇다.

Coloring에서 U를 고르고 Magnitude로 속도 크기를 칠한다.
Glyph로 벡터를 뿌려 역류 구간을 본다.
Stream Tracer로 유선을 그린다.

VTK로 넘긴 뒤 숫자만 뽑을 때는 이 명령이다.

foamToVTK -latestTime

벽 근처 ux의 부호가 바뀌는 x를 재부착으로 잡는다. 스크립트는 H를 0.01 m로 박아 둔다.

import numpy as np
 
# foamToVTK 출력 파일에서 하단 벽(y=0) 데이터 읽기
# (실제 경로는 환경에 따라 조정)
data = np.loadtxt("VTK/backwardFacingStep_500/internal.vtk",
                  skiprows=10)  # VTK 헤더 제외
 
x = data[:, 0]
ux = data[:, 3]  # x방향 속도
 
# y ≈ 0인 셀 필터링 (하단 벽 근처)
wall_mask = data[:, 1] < 0.005
x_wall = x[wall_mask]
ux_wall = ux[wall_mask]
 
# 부호 변화 위치 = 재부착 지점
sign_changes = np.where(np.diff(np.sign(ux_wall)))[0]
if len(sign_changes) > 0:
    x_reattach = x_wall[sign_changes[0]]
    H = 0.01  # 계단 높이 (m)
    print(f"재부착 길이: {x_reattach/H:.1f}H")

흔한 오류와 해결책#

오류 1: `FOAM FATAL ERROR: Cannot find file "0/epsilon"`#

--> FOAM FATAL ERROR:
    cannot find file
    "backwardFacingStep/0/epsilon"

0/epsilon이 없다. k와 같은 뼈대로 하나 만든다.

// 0/epsilon
dimensions      [0 2 -3 0 0 0 0];  // m²/s³
 
// epsilon = C_mu^(3/4) * k^(3/2) / L, L ≈ 0.1*H
internalField   uniform 0.0009;
 
boundaryField
{
    inlet       { type fixedValue; value uniform 0.0009; }
    outlet      { type zeroGradient; }
    wall        { type epsilonWallFunction; value uniform 0.0009; }
    frontAndBack { type empty; }
}

오류 2: 해석 발산 (Floating point exception)#

Floating point exception (core dumped)

완화가 세거나 격자가 성기면 발산한다. fvSolution의 완화계수를 내린다.

relaxationFactors
{
    fields      { p 0.2; }   // 0.3 → 0.2로 감소
    equations   { U 0.5; }   // 0.7 → 0.5로 감소
}

오류 3: residual이 감소하지 않고 진동#

nNonOrthogonalCorrectors가 0인데 비직교가 크면 잔차가 출렁인다. 보정 횟수를 올린다.

SIMPLE
{
    nNonOrthogonalCorrectors    2;  // 0 → 2
    ...
}

예상 결과 정리#

항목	값
Reynolds 수 ( $Re = UH/\nu$ )	1,000
재부착 길이 ( $x_r / H$ )	6–8
최대 역방향 속도	$\approx -0.1$ m/s
수렴 반복 횟수	300–600회

Armaly et al. 1983의 $Re = 800$ 은 약 6H로, 표의 6–8 옆에 놓인다.

x_r \approx 6\text{–}8 \cdot H \quad (Re \sim 10^3)

격자 꼭짓점에서 계단은 y가 0에서 1까지다. scale 0.1을 곱하면 H는 0.1 m다. nu 주석은 $Re = 1 \times 0.01 / 10^{-5} = 1000$ 이라 H를 0.01 m로 적는다. 표의 1000은 그 주석과 같다. 격자의 H를 넣으면 $1 \times 0.1 / 10^{-5} = 10^{4}$ 다. 표의 6–8H는 0.1 m이면 0.6–0.8 m이고, 0.01 m이면 0.06–0.08 m다.

하류 블록의 x는 5에서 20이다. 계단면은 0.5 m, 출구는 2.0 m, 간격은 1.5 m다. 격자 H로는 15H라 0.6–0.8 m가 안에 들어간다. 주석의 H로 나누면 150H가 되고, 6H는 계단면에서 0.06 m, 좌표 0.56 m다. 추출 스크립트의 H도 0.01 m다. 바닥 칸 높이는 $(2/20) \times 0.1 = 0.01$ m이고, y가 0.005 미만인 칸만 고르면 그 칸의 절반이다.

k는 주석의 $1.5 \times (0.05)^2 = 0.00375$ 와 같다. epsilon 식에 $C_\mu = 0.09$ 를 넣으면 $C_\mu^{3/4} k^{3/2} \approx 3.77 \times 10^{-5}$ 이고, 0.0009로 나눈 길이는 0.042 m다. $L \approx 0.1H$ 는 0.01 m 또는 0.001 m라, 그때 epsilon은 0.00377이나 0.0377이어야 한다. 파일의 0.0009는 둘 다 아니다. 역류 $-0.1$ m/s는 입구의 10%이고, 강도 5%의 0.05 m/s보다 두 배 크다.

6H를 열 배로 늘려 버리는 높이#

격자 H를 지키면 표의 Re 1000은 거짓이다. nu를 $10^{-4}$ 로 올려야 그 1000이 돌아온다. 주석의 H를 지키면 scale을 0.01로 내려야 하고, 6H는 0.06 m로 줄어 하류 1.5 m의 앞에 붙는다.

어느 쪽을 이 케이스의 H로 적을 것인가. 입구 1000을 관 전이로 읽으면, 이미 켜 둔 k-ε이 과한지, 계단 박리가 먼저인지는 원관이 흐트러지는 순간에 질문이 남아 있다.