THE ADAPTATION OF IMPROVED NUMERICAL SCHEMES FOR PRESSURE BASED COMPRESSIBLE SOLVER BASED ON OPENFOAM

T.W. Kim

doi:10.6112/kscfe.2026.31.1.041

Preview

Original Article

Journal of Computational Fluids Engineering. 31 March 2026. 41-53
https://doi.org/10.6112/kscfe.2026.31.1.041

THE ADAPTATION OF IMPROVED NUMERICAL SCHEMES FOR PRESSURE BASED COMPRESSIBLE SOLVER BASED ON OPENFOAM

OpenFOAM 기반의 압력 기반 압축성 유동 해석자를 위한 개선된 수치기법의 적용

T.W. Kim¹^†

김 태우¹^†

¹Dept. of Mechanical and Control Eng., Seoul Digital Univ.

¹서울디지털대학교 기계로봇항공학부 기계제어공학전공

^{*Corresponding Author}

License (open-access, http://creativecommons.org/licenses/by-nc/4.0):

This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0) which permits unrestricted noncommercial use, distribution, and reproduction in anymedium, provided the original work is properly cited.

ABSTRACT

This study integrates density-based numerical schemes into a pressure-based solver to enhance compressible flow analysis. The AUSM⁺-up and HLLC schemes were modified for the continuity equation flux and implemented within the OpenFOAM platform. The implementation was validated using 1D shock tube problems, confirming the successful integration of the schemes. Further simulations of 2D airfoils and 3D wings under transonic conditions demonstrate that the enhanced solver is highly effective and reliable. These results indicate that the adapting these schemes significantly enhances the performance of pressure-based solvers for complex compressible flows.

Keywords

Open source CFD

Pressure-based solver

Numerical scheme

OpenFOAM

키워드

오픈소스 전산유체역학

압력 기반 해석자

수치기법

오픈폼

MAIN

1. 서 론
2. 연구 방법
2.1 압력 기반 압축성 유동 해석자
2.2 압축성 유동 해석을 위한 수치 기법 적용
3. 연구 결과
3.1 SOD problem
3.2 RAE2822 익형 해석 결과
3.3 Onera M6 날개 해석 결과
4. 결 론

1. 서 론

유체의 거동에 대한 최신의 연구에서 전산 유체역학(Computational Fluid Dynamics, CFD)은 유체 유동을 해석하고 예측하는 데 있어 실험적 방법과 더불어 핵심적인 도구로 사용되고 있다. 특히 항공우주분야와 같이 마하수가 0.3 미만인 아음속 유동부터 극초음속에 이르는 광범위한 속도 영역을 아우르는 복잡한 유동 현상에 대해 모사해야 할 필요성이 있는 분야에서 전산유체역학을 통한 해석적 연구가 주요한 연구방법으로 발전하고 있다. 이러한 다양한 유동 속도 영역의 유동 해석을 위한 해석자 종류는 크게 밀도 기반(density-based)과 압력 기반(pressure-based)으로 구분된다[1,2,3,4,5].

밀도 기반 해석자는 보존 변수(conservative variables)를 직접 풀어 고속 압축성 유동에 존재하는 충격파와 같은 불연속면을 효과적으로 포착하기 위해 개발되었다[2]. 하지만 마하수가 0에 가까워지는 저속 유동 영역에서는 선속(flux)을 계산하기 위한 행렬계산 과정에서 계산 해가 수렴하지 못하는 계산 체적들이 다수 발생하게 되어 수치적 불안정성 문제가 발생하며, 이로 인해 수렴성이 저하되고 해의 정확도가 낮아지는 한계를 보인다[4]. 이를 해결하기 위해 preconditioning 기법 등이 제안[2]되었지만, 근본적인 문제 해결 방법은 아니라고 알려져 있다.

반면, 압력 기반 해석자는 PISO(Pressure-Implicit with Splitting of Operators), SIMPLE(Semi-Implicit Method for Pressure Linked Equations)과 같은 속도와 압력의 연계 알고리즘을 기반으로 원시 변수(primitive variables)를 주요 변수로 사용하며, 비압축성 및 저속의 아음속 유동에서 효율적이고 강건한 해석 결과를 도출하는 것으로 알려져 있다[3]. 하지만 이런 방식의 해석자를 고속의 압축성 유동 영역에서 발생하는 충격파와 같은 비물리적인 불연속 현상에 적용하는 경우에는 수치 진동이 발생하거나 과도한 수치적 소산(dissipation)으로 인해 불연속면이 명확하게 나타나지 않는 해석 결과를 보여주기도 한다. 특히, 유한체적기법에서 주로 사용되는 collocated grid에서 압력-속도 연계를 위해 많이 사용되는 Rhie-Chow 보간법에 사용되는 격자면에서의 선속 처리가 이러한 수치 소산의 주요 원인 중 하나로 알려져 있다[5].

전 마하수 영역의 유동 해석을 위한 수치 기법의 개발은 밀도 기반 해석자와 압력 기반 해석자가 가진 각각의 장점을 유지하고 단점에 대한 보완하는 방향으로 발전해 오고 있다. 최근에는 압력 기반 해석자의 속도 및 압력 연계 알고리즘은 유지하며 대류 항 계산에 고해상도 수치 플럭스 기법을 접목하는 하이브리드(hybrid) 접근법에 대해 활발한 연구가 진행되고 있다[5,6,7,8]. 이 방식은 압력 기반 해석자의 저속 안정성과 밀도 기반 해석자의 고속 정확도를 동시에 확보하는 것을 목표로 한다. Xisto 등은 PISO 알고리즘에서 Rhie-Chow 보간법을 AUSM 계열 기법[9,10]으로 대체하여 고속 유동에서의 충격파 포착 능력을 향상시키는 연구를 성공적으로 수행했습니다[11,12]. Kraposhin 등은 또한 PISO 알고리즘에 Kurganov-Tadmor central(KT) 기법을 도입하여[13], 복잡한 리만 문제(Riemann problem)를 풀지 않고도 밀도 변화가 큰 유동을 안정적으로 해석하는 방법을 제시하였다. 신재렬 등은 압력 기반 해석자를 중심으로 최근 널리 사용되는 오픈소스 CFD 플랫폼인 OpenFOAM에서 압력 기반의 Coupled 방식 압축성 유동 해석자를 개발하는 연구[14]를 수행하였다. 이처럼 압력 기반 해석자에 AUSM 계열이나 KT와 같은 검증된 플럭스 기법을 결합하는 것은 전 마하수 유동 해석자 개발의 핵심적인 연구 방향으로 제시되고 있다[15].

밀도 기반 해석자에서 사용되는 플럭스 기법 중 HLL(Harten-Lax-van Leer) 계열의 근사 리만 해법은 강건성과 효율적인 해석이 가능한 기법이다[4]. 하지만 초기 HLL 기법은 불연속면을 식별할 수 없어 불연속면에서 과도한 수치적 소산이 발생하는 한계가 있었으며, 이 문제를 해결하기 위해 제안된 기법이 HLLC(Harten-Lax-van Leer-Contact) 모델이다. HLLC 모델은 불연속면에서의 수치적 소산을 획기적으로 감소시킬 수 있으며, 전단파(shear wave)나 다상 경계면과 같은 복잡한 유동 구조를 매우 정확하게 포착할 수 있는 특징이 있다. 그러나 HLLC 기법 역시 고차 정확도 기법과 결합하거나 3차원 고해상도 격자를 사용하여 충격파 등을 해석하는 경우에는 충격파 면에서 수치적 진동과 같은 수치적 불안정성을 나타낼 수 있다[16]. 이를 해소하기 위해 기존의 해결책들은 인위적인 소산을 추가하여 안정성을 확보했지만, 이에 따라 해의 정확도가 낮아지는 특징이 있었다. Fleischmann 등의 연구[17]에서는 이러한 불안정성이 충격파 전파의 수직 방향에서 발생하는 국소적인 저속 유동으로 인해 수치 소산의 스케일링이 부정확해져 발생함을 규명하고, 저속 영역에서의 수치적 소산을 줄여 불안정성을 제거하는 HLLC-LM 기법을 제안하였다.

본 연구에서는 선행연구들에서 확인된 압력 기반의 해석자에 적용되었던 수치 기법들을 개선하여, 대표적인 해석자 검증 문제에 적용하고자 한다. 본 연구에서 사용되는 해석자는 김태우 등의 연구[8]에서 사용된 OpenFOAM 기반의 압력 기반 segregated 해석자를 사용하였다. 본 연구를 통해 기존의 압력 기반 해석자에서 중요하게 사용되는 KT 기법에서 사용되는 격자면에서의 속도 계산을 AUSM Family 계열에서 사용되는 격자면의 수치적 선속을 사용하는 방식으로 개선하고자 한다. 또한 기존 연구에서 제시된 전 마하수 영역에서 안정적인 해를 도출하는 HLLC-LM 기법을 추가적으로 압력 기반 해석자에 적용하고자 한다. 이러한 방식으로 개선된 해석자를 사용하여 불연속적인 유동 현상이 발생하는 검증 문제들을 해석하고 이를 통해 개선된 수치 기법의 효과를 확인하고자 한다.

2. 연구 방법

본 논문에서 설명할 압력 기반의 외부 유동 해석자는 현재 Keysight를 통해 배포되는 OpenFOAM-V2206 및 V2312 버전을 기반으로 개발하였다. OpenFOAM 플랫폼에서 기본적으로 제공하는 압력 기반 비정상 점성 유동 해석자는 PISO 및 PIMPLE(PISO+SIMPLE) 기반의 해석자이며 압축성 및 비압축성 유동을 위한 별도의 해석자가 제공되고 있다[18]. 본 연구에서는 PIMPLE 기반의 알고리즘을 기반으로 압력 기반 압축성 유동 해석자를 개발하였다. PIMPLE 기반 해석자는 기본적으로 비정상 유동의 해석을 위한 알고리즘이나, 각 셀별로 시간을 적용하는 local time 기법을 적용하면 정상 상태 유동 해석도 가능하며 큰 CFL 수를 사용할 수 있는 장점이 존재하여 압력 기반 해석자를 개발할 때 주로 사용되는 해석자이다.

2.1 압력 기반 압축성 유동 해석자

일반적인 유동 해석을 위한 지배방정식은 아래와 같은 식으로 표현된다.

(1)

\frac{\partial ρ}{\partial t} + \nabla ∙ (ρ \vec{u}) = 0

(2)

\frac{\partial ρ \vec{u}}{\partial t} + \nabla ∙ (ρ \vec{u} ∙ \vec{u}) = \nabla ∙ μ [(\nabla \vec{u} + (\nabla \vec{u})^{T}) - \frac{2}{3} (\nabla ∙ \vec{u}) I] - \nabla p

(3)

\frac{\partial ρ h}{\partial t} + \nabla ∙ (ρ \vec{u h}) - \nabla ∙ (α_{e f f} \nabla h) = \frac{\partial p}{\partial t} + \vec{u} ∙ \nabla p

압력 기반으로 지배 방정식을 구성하기 위해서 식 (4)와 같이 밀도와 압력에 관한 관계식을 사용하여 지배방정식을 일부 변경하여야 한다. 압력 기반에서 사용되는 대부분의 속도 및 압력 연계 알고리즘의 경우에는 연속 방정식을 기반으로 반복적인 계산을 통해 수렴해를 구한다. 연속 방정식을 구성하는 밀도를 식 (4)의 관계식을 통하여 식 (5)와 같이 변경하여 연속방정식을 사용한다.

(4)

ρ = \frac{p}{R T} = ψ p

(5)

\frac{\partial ψ p}{\partial t} + \nabla ∙ (ψ p \vec{u}) = 0

식 (5)에서 속도는 SIMPLE이나 PISO 알고리즘에서 사용되는 Rhie-Chow 내삽 기법으로 격자면에서의 면속도를 구하게 된다. 계산된 면속도를 통해 연속 방정식을 만족하는 격자 중심에서 압력을 획득하게 되고 이를 바탕으로 격자 중심점에서의 속도 값을 보정하게 된다[1].

2.2 압축성 유동 해석을 위한 수치 기법 적용

본 장에서는 기존 연구에서 수행한 밀도 변화가 큰 압축성 유동 해석을 위해 압력 기반 해석자에 적용한 KT 기법을 소개하고, 본 연구에서 이를 개선한 내용에 대해 설명하고자 한다. 기존 KT 기법은 리만 솔버를 사용할 필요가 없는 장점이 있으나, 유동 속도의 변화가 큰 영역이나 강한 충격파가 발생하는 경우 수치적 소산이 발생하는 문제가 있는 것으로 알려져 있다[19]. 본 연구에서는 AUSM⁺-up에서 사용되는 수치적 선속 계산 기법을 사용하여 기존 KT 기법의 단점을 개선하고자 하였다. 더하여, Fleischmann등의 연구[17]에서 제시된 저마하수에서의 수치적 소산 문제를 개선한 HLLC-LM 기법을 압력 기반 해석자에 적용하고자 하였다. 압력 기반 해석자의 경우에는 압력을 변수를 사용하는 변형된 연속 방정식을 속도와 압력을 계산하는 알고리즘이 사용되기 때문에, 해당 기법들을 소개했던 선행 연구에서 연속 방정식에 사용되는 질량 선속과 관련된 내용을 사용하였다.

2.2.1 Kurganov-Tadmor 기법의 적용

고속 유동 해석을 위한 선속를 구성하기 위해 사용된 KT(Kurganov-Tadmor) 기법에 대해서는 기존의 연구[8]에서 설명되어 있다. KT 기법은 앞서 서술한바와 같이 복잡한 근사 리만 솔버를 풀지 않고 밀도 변화가 큰 유동 해석을 효율적으로 계산할 수 있는 방법이다. 이 방법은 핵심적인 요소는 Fig. 1과 같이 격자면의 선속을 positive와 negative 방향의 가중합으로 나타내는 것이다. OpenFOAM은 FVM(Finite Volume Method)을 기반으로 개발되었으므로 KT 기법을 FVM의 수치적 선속 계산에 적용하기 위한 방식은 아래와 같다.

(7)

Ψ_{f} ϕ_{f} = Ψ_{f}^{P} (α_{f}^{P} ϕ_{f}^{P} + α_{f}^{P} α_{f}^{\min}) + Ψ_{f}^{N} (α_{f}^{N} ϕ_{f}^{N} - α_{f}^{P} α_{f}^{\min})

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-01/N0500310103/images/jkscfe_2026_311_041_F1.jpg

Fig. 1.

Cell system of FVM to apply the KT scheme

식 (7)에서 사용되는 각 기호의 정의는 아래와 같이 표현된다.

(8)

α_{f}^{P} = \frac{α_{f}^{\max}}{α_{f}^{\max} + α_{f}^{\min}}

(9)

α_{f}^{N} = \frac{α_{f}^{\min}}{α_{f}^{\max} + α_{f}^{\min}}

(10)

α_{f}^{\max} = \max (c_{f}^{P} |\vec{S_{f}}| + ϕ_{f}^{P}, c_{f}^{N} |\vec{S_{f}}| + ϕ_{f}^{N})

(11)

α_{f}^{\min} = - \min (- c_{f}^{P} |\vec{S_{f}}| + ϕ_{f}^{P}, - c_{f}^{N} |\vec{S_{f}}| + ϕ_{f}^{N})

(12)

c_{f}^{P} = \sqrt{γ_{f}^{P} R T_{f}^{P}}, c_{f}^{N} = \sqrt{γ_{f}^{N} R T_{f}^{N}}

$R$ 은 기체 상수이며, 𝛾는 기체의 비열이다. $Ψ_{f}^{P}$ 와 $Ψ_{f}^{N}$ 은 2차 정확도를 가지는 수치적 제한자들(Minmod, vanLeer)를 사용하여 각 격자면에 내삽하여 구한다.

2.2.2 AUSM⁺-up 기법의 적용

AUSM⁺-up은 대류항과 압력항을 나누어 계산하는 플럭스 분할 기법의 한 종류로, 특히 저마하수 영역부터 고마하수 영역까지 광범위한 유동 조건에서 견고하고 정확한 성능을 제공하도록 설계되었다. 압력 기반 해석자는 압력 변수를 사용하는 연속 방정식의 반복적 계산을 통해 속도와 압력을 구하는 방식이 주로 사용된다. 이에 따라 본 연구에선 다음과 같은 AUSM⁺-up 기법에 대한 연구[9]에서 사용되는 선속 계산에서 질량 선속( $\dot{m}$ )을 계산하는 부분만을 적용하였다.

(13)

M_{L / R} = \frac{u_{L / R}}{a_{1 / 2}}

(14)

{\bar{M}}^{2} = \frac{(u_{L}^{2} + u_{R}^{2})}{2 a_{1 / 2}^{2}},

(15)

M_{o}^{2} = \min (1, \max ({\bar{M}}^{2}, M_{\infty}^{2})) \in [0, 1],

(16)

f_{a} (M_{o}) = M_{o} (2 - M_{o}) \in [0, 1],

(17)

M_{1 / 2} = M_{(4)}^{+} (M_{L}) + M_{(4)}^{-} (M_{R}) + M_{d i f f}

(18)

M_{(4)}^{\pm} (M) = \{\begin{cases} M_{(1)}^{\pm} M_{(1)}^{\pm} if | M | \geq 1 \\ M_{(2)}^{\pm} (1 \mp 16 β M_{(2)}^{\mp}) otherwise \end{cases}

(19)

M_{(1)}^{\pm} = \frac{1}{2} (M \pm | M |)

(20)

M_{(2)}^{\pm} = \pm \frac{1}{4} (M \pm 1)^{2}

(21)

M_{d i f f} = - \frac{K_{p}}{f_{a}} \max (1 - σ {\bar{M}}^{2}, 0) \frac{p_{R} - p_{L}}{ρ_{1 / 2} a_{1 / 2}^{2}}

여기서, $K_{p}$ =0.25, 𝛽=1/8을 적용하였다. 압력 기반 해석자에서 사용되는 속도 및 압력 연계 알고리즘에서 사용되는 질량 선속은 아래의 정의를 통해 계산하였다.

(22)

Ψ_{f} ϕ_{f} = Ψ_{f}^{L} a_{f} α_{f}^{L} M_{1 / 2} + Ψ_{f}^{R} a_{f} α_{f}^{R} M_{1 / 2}

(23)

α^{L} = \{\begin{array}{l} 1 if M_{1 / 2} > 0, \\ 0 otherwise \end{array}, α^{R} = 1 - α^{L}

2.2.3 HLLC 기법의 적용

HLLC 기법은 충격파 및 접촉면을 명확하게 포착하는 데 효과적인 특성을 가지나, 저마하수 유동 영역에서는 수치적 확산 증가와 압력 예측의 부정확성 등 한계가 존재한다. 이러한 문제를 개선하기 위해, 본 연구에서는 Fleischmann[17]등이 제안한 shock-stable 및 low-Mach 보정 기법을 도입하였으며, 이 기법은 기존 HLLC 기법의 파동 속도 $S_{L / R}$ 에 감쇠 함수 𝜙를 곱하여 낮은 마하 수 안에서의 대류 소산과 음향 소산의 균형을 전체 소산의 감소를 통해 이루어지도록 식을 수정하였다. 새로운 파동 속도 $S_{L / R}^{H L L C - L M}$ 은 식 (24)과 같다.

(24)

S_{L / R}^{H L L C - L M} = ϕ ∙ S_{L / R}

감쇠 함수 𝜙는 마하수 기반의 함수로, 저마하수 영역에서는 수치 소산을 효과적으로 억제하고, 고마하수 영역에서는 기존 HLLC의 성능을 유지할 수 있도록 설계되었다. 감쇠 함수 𝜙는 국소 마하수 $M_{l o c a l}$ 에 따르고, 식 (25)과 같다.

(25)

ϕ = \sin (\min (1, \frac{M_{local}}{M_{limit}}) ∙ \frac{π}{2})

여기서 국소 마하수 $M_{local}$ 는 셀 경계면 양쪽에 위치한 유동 상태를 고려하여 정의되며, 좌측 상태의 속도 $u_{L}$ 와 음속 $c_{L}$ , 우측 상태의 속도 $u_{R}$ 와 음속 $c_{R}$ 를 기반으로 식 (26)과 같다.

(26)

M_{local} = \max (|\frac{u_{L}}{c_{L}}|, |\frac{u_{R}}{c_{R}}|)

위 방식은 저마하수 조건에서 비정상적인 대류 및 음향 소산 간의 불균형을 방지하고, 전체 수치 소산을 부드럽게 감소시켜 해석의 안정성과 정량적 정확도를 향상시키는 데 효과적이다. 이를 통해 고마하수 영역에서는 기존 HLLC 기법의 성능을 유지하면서도 저마하수 영역에서는 수치적 안정성과 정확성이 향상되도록 개선되었다. 압력 기반 해석자에서 격자면에서 수치적 선속은 아래와 같이 정의할 수 있다.

(27)

Ψ_{f} ϕ_{f} = \{\begin{matrix} Ψ_{f}^{L} ϕ_{f}^{L} & if S_{L} \geq 0, \\ Ψ_{f}^{* L} ϕ_{f}^{* L} = Ψ_{f}^{L} ϕ_{f}^{L} + S_{L} (U_{*_{L}} - Ψ_{f}^{L}) & if S_{L} < 0 \cap S_{*} \geq 0, \\ Ψ_{f}^{* R} ϕ_{f}^{* R} = Ψ_{f}^{R} ϕ_{f}^{R} + S_{R} (U_{*_{R}} - Ψ_{f}^{R}) & if S_{R} > 0 \cap S_{*} \leq 0, \\ Ψ_{f}^{R} ϕ_{f}^{R} & if S_{R} \leq 0, \end{matrix}

중간상태를 나타내는 $U_{*_{L}}$ , $U_{*_{R}}$ 과 $S_{L}$ , $S_{R}$ 은 아래의 식을 통해서 결정된다.

(28)

U_{*_{K}} = \frac{S_{K} - u_{K}}{S_{K} - S_{*}} Ψ_{K}, Ψ = ρ, K = L, R

(29)

S_{L} = \min (u_{L} - c_{L}, \hat{u} - \hat{c}), S_{R} = \max (u_{R} - c_{R} \hat{u} + \hat{c})

여기서, $\hat{u}$ 와 $\hat{c}$ 는 아래와 같은 Roe 평균 정의를 통해 결정된다.

(30)

\hat{u} = \frac{u_{L} \sqrt{ρ_{L}} + u_{R} \sqrt{ρ_{R}}}{\sqrt{ρ_{L}} + \sqrt{ρ_{R}}}

(31)

{\hat{c}}^{2} = \frac{c_{L}^{2} \sqrt{ρ_{L}} + c_{R}^{2} \sqrt{ρ_{R}}}{\sqrt{ρ_{L}} + \sqrt{ρ_{R}}} + \frac{1}{2} \frac{\sqrt{ρ_{L}} \sqrt{ρ_{R}}}{{(\sqrt{ρ_{L}} + \sqrt{ρ_{R}})}^{2}} {(u_{R} - u_{L})}^{2}

$S_{*}$ 는 contact wave speed로 아래의 식을 통해서 구하게 된다.

(32)

S_{*} = \frac{p_{R} - p_{L} + ρ_{L} u_{L} (S_{L} - u_{L}) - ρ_{R} u_{R} (S_{R} - u_{R})}{ρ_{L} (S_{L} - u_{L}) - ρ_{R} (S_{R} - u_{R})}

3. 연구 결과

본 연구를 통해 개선된 압력 기반 압축성 유동 해석자의 해석 결과 검증을 위해 불연속적인 유동현상이 발생하는 3가지 해석 케이스에 대한 해석을 수행하였다. 아래에 설명되는 비교 해석 케이스들의 수치적 정해나 실험적 결과와 개발된 해석자의 해석 결과를 비교하여 개발된 해석자의 유용성을 검토하였다.

3.1 SOD problem

압력 기반 압축성 해석자의 개선된 수치기법의 결과 확인을 위해 해석자 검증 문제에 널리 사용되는 Sod[20] 문제의 해석을 수행하였다. 유차원화 된 변수값을 사용하는 OpenFOAM을 사용하기 위해, 288 K의 기준온도와 101325 Pa의 기준압력을 사용하여 아래의 무차원 초기조건에 따라 초기 조건을 설정하였다.

(33)

(ρ, u, p) = \{\begin{cases} (1, 0, 1) if x \leq 0 \\ (0.125, 0, 0.1) if x > 0 \end{cases}

양끝 단의 길이는 $x$ =±0.5 m이고 끝단의 경계조건으로 외삽조건을 부여하였다. 설정된 유차원 계산시간은 8.23×10^-4 sec이며 해당 유차원 계산 시간에 대한 무차원화 시간은 $t_{f}$ =0.14이다. 단순화된 사각 격자를 사용하였고, 사용된 격자수는 $N$ =100개 이다.

Fig. 2는 무차원 계산시간은 $t_{f}$ =0.14에서 개발된 해석자의 결과를 수치적 정해와 함께 나타낸 그림이다. 개발된 해석자(Present)의 해석 결과는 기존에 사용되는 KT 기법외에 본 연구를 통해 적용된 AUSM⁺-up 및 HLLC 결과를 각각 표시하였다. 해석 결과들과 같이 표시한 정해와 비교하여 세가지 방식의 수치적 기법의 해석 결과가 팽창 영역 및 충격파의 접촉면에서 거의 유사한 수준의 해를 도출하고 있음을 알 수 있다. 이를 통해 본 연구에서 사용된 세가지 수치 기법은 1차원 충격파 해석 결과에서 압력 및 속도, 온도등의 중요한 유동 변수에 대해 동일한 수준의 해를 도출할 수 있음을 보여주고 있다. 다만, 본 해석 케이스의 경우에는 비점성 가정을 기반으로 격자의 형상이 단순한 1차원 유동 해석을 수행하였기 때문에 본 연구에서 사용된 기법에 의한 차이를 확인하기에는 한계가 있는 것으로 판단된다. 그러므로 본 해석 케이스를 통해서는 본 연구에서 사용되는 압력 기반 전산유체역학 해석자에 세 기법이 적절히 적용되었음만을 확인할 수 있었다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-01/N0500310103/images/jkscfe_2026_311_041_F2.jpg

Fig. 2.

Density, Velocity, Pressure, Temperature, Mach No. results of 1D SOD problem

3.2 RAE2822 익형 해석 결과

RAE2822 익형은 NASA가 공개된 웹페이지[21]를 통해 천음속 유동 영역(M_∞ = 0.729)의 해석 및 실험 조건과 결과를 제공되어 해석자의 천음속 유동 해석 능력을 검증하기에 적절한 케이스로 널리 알려져 있다. 본 연구에서는 Fig. 3과 같은 격자를 사용하였으며, 해당 격자는 이전에 수행된 2차원 익형의 해석 연구[8]에서 사용된 격자와 동일하게 35,892개의 복합격자로 구성하였다. 해석에 사용된 격자는 익형의 표면 위쪽의 경계층 부근에서는 사각 형태의 격자를 구성되고 그 이외에는 삼각형의 격자를 구성하는 방식으로 작성되었다. 난류 모델은 k-omega SST 난류 모델을 사용하였고, 익형의 표면은 no-slip 경계조건과 함께 최소 격자 높이(y⁺)의 평균값은 30이상으로 하여 wall function 조건을 부여하였다. 원방 경계 조건은 이전의 연구에서 사용된 외부 유동 해석을 위해 원방 경계 부근에서의 Riemann 조건식을 기반으로 별도로 개발된 characteristic 경계 조건을 사용하였다. 정상 상태 해석을 위한 local time Stepping의 시간 전진 기법을 사용하였다. 해석된 결과는 실험 결과 및 다른 전산유체역학 해석자의 결과와 비교하여 개발된 해석자의 수치기법의 해석 유용성을 검토하였다. 유동의 해석 조건은 M_∞ = 0.729, Reynolds Number = 6.5×10⁶, P_∞ = 101325pa, T_∞ = 255.55 K, AoA(Angle of Attack) = 2.92°이다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-01/N0500310103/images/jkscfe_2026_311_041_F3.jpg

Fig. 3.

Grids system for RAE2822 airfoil

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-01/N0500310103/images/jkscfe_2026_311_041_F4.jpg

Fig. 4.

Comparison of pressure coefficient results on airfoil surface

Fig. 4는 본 연구에서 적용한 수치기법(KT, AUSM⁺-up, HLLC)을 사용하여 획득한 RAE2822 익형 표면의 압력 계수 결과를 실험 결과 및 타 수치해석자의 해석 결과와 비교한 그림이다. 본 연구에 의해 적용된 세가지 수치기법의 표면 압력결과는 대체로 실험 결과 및 타 수치해석과 유사한 결과를 보여주고 있으나, 기법의 특징에 따라 익형 윗면에 해당하는 흡입면에서 발생하는 충격파에 따른 압력 변화의 차이가 확인된다. KT 기법대비 AUSM⁺-up 및 HLLC 기법을 사용한 경우, 충격파에 의한 압력 계수의 기울기가 다소 완만하게 예측되는 것으로 확인된다. 이는 본 연구에서 사용된 central 및 upwind 기법의 수치해석적 차이에 의해 발생하는 것으로 판단된다. 대표적인 central scheme으로 분류되는 KT 기법은 충격파등의 불연속적인 유동 현상에 대해 비교적 강건하게 예측하는 대신 불연속적인 면 주변으로 수치적 불안정성을 보이는 특징을 가지고 있다. Upwind scheme으로 분류되는 두 기법은 AUSM⁺-up 및 HLLC 불연속적인 유동 현상을 강건하게 예측할 수 있으나, 다소간의 수치적 확산을 내포하고 있는 알려져 있다. 이 두 방식의 차이에 따른 충격파 전후에서의 유동 현상이 다소간의 차이를 보여 주고 있는 것으로 판단된다. 다만, 그 해석 차이의 정도가 크지 않기 때문에 본 연구에서 적용된 세가지 수치기법 모두 천음속 영역에서의 불연속적인 유동 현상을 비교적 잘 예측하는 것으로 확인된다.

3.3 Onera M6 날개 해석 결과

본 연구에서 적용된 수치기법의 3차원 해석 결과를 확인하기 위해서 3차원 날개 케이스는 Onera M6 날개 케이스[21]해석을 수행하였다. 앞서 2차원 익형 해석 케이스와 마찬가지로 NASA의 웹사이트[21]의 해석의 유동 조건 및 Fig. 5의 해석 격자를 사용하였다. Structured 형식의 hexahedral 격자를 OpenFOAM에서 사용가능한 unstructured 형식으로 변환하여 해석에 사용하였다. 원방 경계 조건은 앞서 2차원 익형에서 사용한 characteristic 경계 조건을 사용하였고 난류 모델은 k-omega SST 모델을 사용하였고, 날개의 표면은 no-slip 경계조건과 함께 최소 격자 높이(y⁺)의 평균값은 30이상으로 하여 wall function 조건을 부여하였다. 본 해석 케이스의 경우에는 날개 표면에서 예측되는 충격파의 비정상적인 현상이 예측되지 않는 것으로 알려져 있으므로 빠른 계산 해의 획득을 위해서 정상유동 해석을 위한 local time stepping 기법을 적용하여 해석을 수행하였다. 유동의 해석 조건은 M_∞ = 0.8395, Reynolds Number = 11.72×10⁶, P_∞ = 315978pa, T_∞ = 255.55 K, AoA = 3.06°이다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-01/N0500310103/images/jkscfe_2026_311_041_F5.jpg

Fig. 5.

Grids system for Onera M6 wing

Fig. 6은 날개 표면의 압력 계수 예측 결과를 날개 스팬 위치(y/b)에 따라 실험결과 및 타 수치해석 결과와 비교한 그림이다. 날개 스팬 전체 위치에서 본 연구를 통해 세가지 수치기법을 통해 획득한 결과들이 Fig. 6에 제시된 대부분의 스팬 위치에서 실험 결과 및 타 수치해석 결과와 비교적 유사한 결과를 보여주고 있음을 확인할 수 있었다. 다만, y/b가 0.8인 부분에서의 날개 표면의 앞전과 뒷전에서 각각 발달하는 충격파가 합쳐지는 영역으로 타 수치해석 결과와는 다소 차이를 보이고 있으나, 실험 결과와 유사한 경향을 보이고 있어 본 연구를 통해 수행된 수치기법의 적용이 적절하게 이루어졌음을 확인할 수 있다. 또한, 전체적인 충격파에 의한 표면 압력 계수의 예측 경향이 앞서 2차원 익형 해석 결과와 유사하게 KT 기법과 나머지 2개의 기법간의 다소간의 차이를 보이고 있으나, 전체적인 차이는 크지 않아 본 연구를 통해 적용된 수치기법이 3차원 해석에서도 유용하게 활용 될 수 있음을 확인할 수 있었다.

https://cdn.apub.kr/journalsite/sites/kscfe/2026-031-01/N0500310103/images/jkscfe_2026_311_041_F6.jpg

Fig. 6.

Comparison of pressure coefficient results on wing surface according to spanwise ratio

4. 결 론

본 연구에서는 압축성 유동 해석을 위해 압력 기반 해석자에 밀도 기반 해석자에서 널리 사용되는 수치 기법을 적용하고, 수치기법이 적용된 압력 기반 해석자의 압축성 유동 해석자의 해석 유용성을 검토하고자 하였다. 이를 위해서 본 연구에서 개발된 해석자를 사용하여 1차원 충격파 케이스 및 익형 및 날개 케이스의 해석을 수행하였다.

1차원 충격파 해석을 통해 기존에 사용되었던 수치기법과 본 연구를 통해 새롭게 적용된 수치기법들의 해석 결과가 유사하게 나타나고 있음을 확인하여 본 연구에서 수행한 수치기법의 적용이 적절하게 이루어졌음을 확인할 수 있었다. 충격파의 예측 정확도가 중요한 천음속 영역의 익형 및 날개 해석을 통해서 본 연구에서 개발된 압력 기반 해석자가 천음속 영역의 유동에 영향을 받는 물체에 대한 해석이 적절하게 수행되고 있음을 확인하였다.

앞서 제시된 익형 및 날개에서의 유동 해석에서 압력 기반의 압축성 유동 해석자에 적용된 수치기법에 따라 충격파의 강도 및 전후의 유동 현상 예측에는 다소간의 차이가 발생하고 있음을 확인할 수 있었다. 이와 같은 결과의 차이는 본 연구에서 적용된 수치기법 자체가 가지고 있는 차이에 따라 발생한 것으로 판단된다. 이를 통해서도 본 연구에서 수행한 밀도 기반 해석자에서 주로 사용되던 수치기법을 압력 기반에 적용한 연구의 결과가 적절하게 이루어졌다고 판단할 수 있었다.

본 연구를 통해 압축성 유동 해석에 유용하게 활용될 수 있는 압축성 유동 해석자의 개발이 이루어졌음을 확인할 수 있었다. 해당 연구의 결과를 바탕으로 향후에는 고정익 및 회전익, 고속철도와 같은 밀도 기반의 해석자가 주로 사용되는 분야에 대한 추가적인 적용 연구를 수행하며 본 연구를 통해 적용된 수치기법을 속도와 압력을 내재적으로 연계할 수 있는 압력 기반 coupled 해석자등에 적용하는 추가적인 연구를 수행하고자 한다.

Acknowledgements

본 연구는 2025년도 서울디지털대학교 대학연구비의 지원으로 수행되었습니다.

References

2002, Ferziger, J.H. and Peric, M., Computational Methods for Fluid Dynamics, 3^rd edition, Springer, New York.

2007, Hirsch, C., Numerical Computation of Internal & External Flows, 2^nd edition, B-H publications, Burlington.

2007, Versteeg, H.K. and Malalasekera, W., An Introduction to Computational Fluid Dynamics, 2^nd edition, Pearson Prentice Hall, Harlow.

1997, Toro, E.F., Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer Berlin, Heidelberg.

2002, Jemcov, A., Stephens, D.W. and Wahono, S., “A High-Resolution Primitive Variable Solver for Compressible Flow Simulation,” AIAA Scitech 2020 Forum, AIAA 2020-0567, Orlando, Florida.

10.2514/6.2020-0567

2015, Kraposhin, M., Bovtrikova, A. and Strijhak, S., “Adaptation of Kurganov-Tadmor Numerical Scheme For Applying in Combination With the PISO Method in Numerical Simulation of Flows in a Wide Range of Mach Numbers,” Procedia Comput. Sci., Vol.66, p.43.

10.1016/j.procs.2015.11.007

2016, Ong, K.C. and Chan, A., “A pressure-based Mach-uniform method for viscous fluid flows,” Int. J. Comput. Fluid Dyn., Vol.30, No.7-10, p.516.

10.1080/10618562.2016.1245417

2017, Kim, T.W. and Gill, J.H., “Development of Improved Pressure-Based Solver Algorithm for Compressible Flow,” J. Comput. Fluids Eng., Vol.22, No.4, p.28.

10.6112/kscfe.2017.22.4.028

2006, Liou, M.S., “A sequel to AUSM, Part II: AUSM⁺-up for all speeds,” J. Comput. Phys., Vol.214, No.1, p.137.

10.1016/j.jcp.2005.09.020

2022, Krishnan, M., Ray, A. and Peetala, R., “Numerical Analysis of High Speed Flow Applications using Various Flux Schemes,” Trends Sci., Vol.19, No.18, p.5813.

10.48048/tis.2022.5813

2013, Xisto, C.M., Páscoa, J.C. and Oliveira, P.J., “A pressure-based method with AUSM-type fluxes for MHD flows at arbitrary Mach numbers,” Int. J. Numer. Methods Fluids, Vol.72, No.11, p.1165.

10.1002/fld.3781

2014, Xisto, C.M., Páscoa, J.C. and Oliveira, P.J., “A pressure-based high resolution numerical method for resistive MHD,” J. Comput. Phys., Vol.275, p.323.

10.1016/j.jcp.2014.07.009

2000, Kurganov, A. and Tadmor, E., “New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations,” J. Comput. Phys., Vol.160, Issue 1, p.241.

10.1006/jcph.2000.6459

2021, Shin, J.R. and Kim, T.W., “Enhanced Pressure Based Coupled Algorithm to Combine with Pressure-Velocity-Enthalpy for all Mach Number Flow,” Int. J. Aeronaut. Space Sci., Vol.22, p.489.

10.1007/s42405-020-00337-9

2019, Ajaero, C.U., “Validation of the Kurganov-Tadmor/Kurganov-Noelle-Petrova Scheme for Rotating Detonation Engine Simulations using OpenFOAM,” Master’s Thesis, Concordia University, Montreal, Quebec, Canada.

2020, Marzouk, O.A., “The Sod gasdynamics problem as a tool for benchmarking face flux construction in the finite volume method,” Sci. Afr., Vol.10, e00573.

10.1016/j.sciaf.2020.e00573

2020, Fleischmann, N., Adami, S. and Adams, N.A., “A shock-stable modification of the HLLC Riemann solver with reduced numerical dissipation,” J. Comput. Phys., Vol.423, 109762.

10.1016/j.jcp.2020.109762

http://www,openfoam.com, KEYSIGHT, Bracknell, UK.

2002, Leveque, R.J., Finite Volume Methods for Hyperbolic Problems, Cambridge University Press.

1978, Sod, G.A., “A survey of several finite difference methods for systems of non-linear hyperbolic conservation laws,” J. Comput. Phys., Vol.27, Issue 1, p.1.

10.1016/0021-9991(78)90023-2

https://www.grc.nasa.gov/www/wind/valid/archive.html, NPARC Alliance Verification and Validation Archive, NASA.

Journal of Computational Fluids Engineering ISSN:1598-6071(Print) 3022-6252(Online) 한국전산유체공학회지

Preview

THE ADAPTATION OF IMPROVED NUMERICAL SCHEMES FOR PRESSURE BASED COMPRESSIBLE SOLVER BASED ON OPENFOAM

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(7)

Fig. 1.

Cell system of FVM to apply the KT scheme

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

Fig. 2.

Density, Velocity, Pressure, Temperature, Mach No. results of 1D SOD problem

Fig. 3.

Grids system for RAE2822 airfoil

Fig. 4.

Comparison of pressure coefficient results on airfoil surface

Fig. 5.

Grids system for Onera M6 wing

Fig. 6.

Comparison of pressure coefficient results on wing surface according to spanwise ratio

Acknowledgements

References