STUDY ON BOUNDARY CONDITIONS FOR CFD/FREE WAKE COUPLED HYBRID METHOD

D. Park; M.H. Cho; S.H. Park

doi:10.6112/kscfe.2024.29.4.102

Preview

Original Article

Journal of Computational Fluids Engineering. 31 December 2024. 102-116
https://doi.org/10.6112/kscfe.2024.29.4.102

STUDY ON BOUNDARY CONDITIONS FOR CFD/FREE WAKE COUPLED HYBRID METHOD

CFD/자유 후류 모델을 결합한 하이브리드 기법의 원방 경계 조건 연구

D. Park¹

M.H. Cho¹

S.H. Park²^†

박 다운¹

조 민형¹

박 수형²^†

¹Department of Aerospace Information Engineering, Konkuk University

²Department of Mechanical and Aerospace Engineering, Konkuk University

¹건국대학교 항공우주정보시스템공학과

²건국대학교 기계항공공학부

^{*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

In the aerodynamic analysis of rotors, low-fidelity BEM-wake modeling methods have limited accuracy, while high-fidelity Computational Fluid Dynamics(CFD) poses significant computational challenges. To address these limitations, a hybrid method combining CFD with a free-wake model has been developed. This employs CFD in the vicinity of the rotor blades and the free-wake model for the far-field wake. The specification of CFD boundary conditions at the coupling interface is critical. Therefore, this study investigates suitable boundary conditions for hybrid methods in rotor problems. Each of boundary condition is based on the Bernoulli equation and the specified constant total enthalpy. It is demonstrated that using boundary condition the based on constant total enthalpy provides more accurate aerodynamic analysis results over a wider range of flight speed compared to it based on the Bernoulli equation.

Keywords

Computational fluid dynamics

Rotor aerodynamics

Free wake

Hybrid methods

Boundary condition

키워드

전산유체역학

로터공기역학

자유 후류

혼합연계기법

경계 조건

MAIN

1. 서 론
2. 수치해석기법
2.1 자유 후류모델
2.2 CFD/자유 후류모델 연계
2.3 해석 모델 및 격자
2.4 해석 조건
3. 해석 결과
3.1 저속 회전 로터
3.2 고속 회전 로터
3.3 하이브리드 기법의 효율성 분석
4. 결 론

1. 서 론

일반적인 회전익기 로터에서도 복잡한 유동 현상이 발생하고 후류의 거동으로 인한 블레이드, 동체 간섭 현상이 성능에 중요하게 작용해 로터의 공력 해석이 중요하다. 특히 최근에 연구가 활발히 진행 중인 형상이 복잡한 차세대 복합형 회전익기 및 차세대 수직이착륙기는 일반적인 회전익기보다 더욱 복잡한 유동 현상이 발생해 공력 해석이 필수적이다. 이것을 기존에 사용되던 저 정밀도 BEM(Blade Element Theory)-후류 모델링 기법을 사용해 분석하면 정확도가 감소하고, 고정밀도 전산유체해석으로 분석하면 계산의 난이도와 계산량이 상당히 증가한다. 따라서 두 기법을 보완한 기법이 필요하다. 각 기법의 단점을 보완하기 위해 전산유체해석기법과 자유 후류 모델을 연계 해석하는 CFD/자유 후류모델을 연계한 하이브리드 기법(CFD/free-wake coupled hybrid method)이 개발되었다. 이것은 로터 블레이드 근방 영역은 Euler 방정식 또는 Navier-Stokes 방정식을 풀어 전산유체해석을 수행한다. 원방의 후류는 자유 후류 모델을 이용하여 모델링 한다. 자유 후류 모델로 후류를 모델링하기 위해서 전산유체해석 결과를 이용하고, 자유 후류모델로 모델링된 후류 구조를 비오-사바르(Biot-Savart) 법칙을 통해 전산유체해석에 필요한 유동 조건을 구한다.

자유 후류 모델을 전산유체기법에 보정 하는 방법은 일반적으로 2가지로, 장속도기법(FVA, Field Velocity Approach)과 경계 조건 보정 방법(boundary correction approach)과 이다. 장속도 기법은 모델링 된 후류가 전산유체해석이 이뤄지는 모든 영역에 유도 속도를 보정 하는 방법이다. 이 기법에 관한 대표적인 연구는 Lee[1], Inada[2], Sitaraman[3] 등이 있다. 이 기법은 각 와류 요소가 모든 전산유체해석영역에 유도되는 속도를 계산해야 하므로 계산 시간이 느리고, 계산량이 증가한다. 경계 조건에 유도 속도를 보정 하는 방법은 모델링 된 후류가 전산유체해석 영역 경계면에 유도하는 속도를 비오-사바르 법칙으로 계산한다. Wie[4], Rajmohan[5], Min[6] 등이 이 기법을 사용해 다양한 연구를 진행하였다. 경계 조건 보정 방법은 격자의 경계점에만 유도 속도를 보정 하므로 계산 시간이 빠른 장점이 있어 이번 연구에서는 해당 방식을 채택하였다.

로터 공력 해석을 위한 자유 후류 모델은 각 블레이드 제어점에서 발생하는 모든 후류 요소를 고려하는 방법 [1, 4, 5] 과 하나의 끝단 후류만을 고려하는 방법 [6, 7, 8, 9]이 있다. 끝단 후류만 고려하는 것은 블레이드 스팬 방향에 따른 와류를 자세히 모사하는 것은 어렵다. 하지만 모든 와류 요소를 고려하면 모델링의 복잡성이 증가하고, 와류 필라멘트 요소가 블레이드 패널 수만큼 늘어나기 때문에 계산 간격이 줄어들거나 계산이 진행됨에 따라 와류 요소가 증가할수록 계산 비용이 상당히 증가하는 단점이 있다. Min[6], Leishman[8], Yang[9] 등 많은 선행 연구자들이 끝단 후류만을 고려한 자유 후류 모델링으로 타당한 결과를 얻었기 때문에 끝단 후류만을 고려하는 것은 로터 후류를 예측할 때 문제를 단순화할 수 있는 합리적인 접근법이다. 따라서 본 연구에서는 블레이드 스팬 방향에 따른 제어점에서 발생하는 모든 와류 요소를 개별적으로 고려하지 않고, 이를 끝단 후류로 치환하여 하이브리드 기법을 개발하였다.

하이브리드 기법의 계산 영역은 CFD와 자유 후류 모델의 결합이 발생하는 경계면의 벽과의 거리가 일반적인 아음속 전산유체해석 계산 영역에 비해 매우 가깝다. 특히 중첩 격자 기법을 사용해 배경 격자가 있는 연계 해석 기법보다도 본 연구에서 개발한 해석 기법은 블레이드의 더 좁은 근방 영역만을 해석한다. 따라서 리만 불변량(Riemann invariant) 경계 조건과 같이 일반적인 아음속 원방 경계 조건을 사용하면 계산의 정확도가 낮아진다. 그러므로 특수한 경계 조건을 지정하는 것이 중요하다. 그렇지만 일반적으로 간단한 1차원 파동방정식의 고유값과 고유 벡터를 이용한 리만 불변량 경계 조건을 사용하거나[1, 10, 11] 경계 조건에 대한 자세한 서술은 없으며 경계 조건의 영향에 대한 자세한 공력 특성 분석에 관한 연구는 부족하다. 따라서 본 논문에서는 로터 문제에 적합한 CFD와 자유 후류 모델을 결합한 하이브리드 기법의 원방 경계 조건 연구를 진행하였다. 베르누이 방정식(Bernoull’s equation)을 활용한 경계 조건과 전엔탈피(total enthalpy)가 일정하다는 조건을 사용해 경계 조건을 개발하고 각 경계 조건의 특징과 계산의 정확도에 대해 분석하였다.

2. 수치해석기법

2.1 자유 후류모델

시간과 공간에 대한 후류의 위치의 지배방정식은 식 (1)[12, 13] 로 표현된다.

(1)

\frac{\partial \vec{r} (ψ, ζ)}{\partial ψ} + \frac{\partial \vec{r} (ψ, ζ)}{\partial ζ} = \frac{1}{Ω} [\vec{V_{\infty}} + \sum_{i = 1}^{N_{v}} \vec{V_{i n d}} (\vec{r} (ψ, ζ), \vec{r} (ψ_{i}, ζ))]

위 식에서 $\partial \vec{r} (ψ, ζ)$ 은 후류 계산점(collocation point)의 위치 벡터이다. 이때 𝜓, 𝜁은 각각 방위각과 후류 요소의 wake age이다. 위 식의 우변에서 𝛺는 로터의 회전 속도, $\vec{V_{\infty}}$ 는 자유류의 속도이며, $\sum_{i = 1}^{N_{v}} \vec{V_{i n d}} (\vec{r} (ψ, ζ), \vec{r} (ψ_{i}, ζ))$ 은 각 제어점에서 모든 후류 필라멘트 요소 $N_{v}$ 에 의해 유도되는 속도의 합이다. 위 식을 풀기 위해 Bagai[12, 13]이 제안한 5-계산점 중앙 차분법을 적용한 Pseudo-Implicit 시간 전진 기법을 사용하였다.

후류의 위치를 결정하기 위한 유도 속도는 비오-사바르 법칙을 이용하였다. 이것을 통해 와류강도 $Γ_{v}$ 인 와류 필라멘트의 부분 요소 길이 $d \vec{l}$ 에서 거리 $\vec{r}$ 만큼 떨어진 위치에서 유도 속도의 증분값 $\vec{v}$ 을 식 (2)로 구할 수 있다. 특정 위치에서의 전체 속도는 와류 필라멘트의 길이를 따라 적분해 구할 수 있다. 본 논문에서는 곡선 와류를 직선 와류로 연결된 것으로 가정하였다. 또한 점성에 의한 와류 확산 효과를 고려하고 특이점(singularity) 문제를 해결하기 위해 Squire이 제안한 와류핵성장(vortex core growth)[14, 15] 모델을 사용하였다. 와류 핵 반지름( $r_{c}$ )는 와류 초기 반지름( $r_{i n i t i a l}$ ), 1.25643인 Lamb-Oseen 상수(𝛼), 평균유효점성계수(𝛿), 동점성계수(𝜈), 로터의 회전속도(𝛺), 후류 나이(wake age) 𝜁에 대한 함수로 식 (3)으로 표현된다.

(2)

\vec{d v} = \frac{Γ_{v}}{4 π} \frac{\vec{d} \times \vec{r}}{| \vec{r} |^{3}}

(3)

r_{c} (ζ) = \sqrt{(r_{i ∋ tial}^{2} + 4 α δ v ζ / Ω)}

2.2 CFD/자유 후류모델 연계

CFD와 자유 후류 기법 연계를 위해 모든 격자 셀에 와류 요소로 인한 유도 속도를 직접 반영하는 장속도 기법이 아닌 경계면에만 유도 속도를 보정 하는 방법을 적용하였다. 경계면을 통한 CFD/자유 후류 기법 연계 방법은 Fig. 1을 통해 알 수 있다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F1.jpg

Fig. 1.

Principle of hybrid method through the CFD boundary interface

자세한 연계 절차는 Fig. 2로 알 수 있다. 각 블레이드에서 발생하는 후류의 강도는 전산유체해석을 통해 블레이드에서 발생하는 양력을 통해 구한다. 본 논문에서는 계산의 효율성을 위해 끝단 후류만 고려한다. 각 블레이드의 끝단 후류의 강도는 Kutta-Joukowski 이론[15]과 양력선(lifting line) 이론[15]을 이용해 얻을 수 있다. 반지름 방향으로 블레이드 요소의 양력을 계산해 각 요소에서 발생하는 와류의 강도를 구할 수 있다. Kutta-Joukowski 이론을 통해 단위 반지름당 양력은 식 (4)로 표현된다. 식 (4)에서 양력( $L$ )과 양력계수( $C_{l}$ )는 전산유체해석을 통해 구한다. 𝜌는 밀도, $r$ 은 반지름, c는 각 위치에서의 코드 길이를 의미하며, $V (r)$ 은 각 반지름 위치에서 회전방향 속도와 전산유체해석을 통해 얻은 유입류 속도로 계산된 자유류 속도를 나타낸다. 반지름 방향으로 블레이드 요소에서 발생하는 각 말굽 와류를(Fig. 3) 중첩 시켜 블레이드 전체를 의미하는 하나의 말굽 와류로 나타낼 수 있다. 즉, 양력선 이론을 통해 각 반지름 방향 블레이드 요소에서 발생하는 양력으로 구한 와류 강도를 적분해 하나로 표현된 말굽 와류의 자유끌림 와류 강도를 계산할 수 있다.

(4)

d L (r) = ρ V (r) Γ (r) d r = \frac{1}{2} ρ V (r)^{2} c C_{l} (y) d r

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F2.jpg

Fig. 2.

Procedure of hybrid method with tight coupling method

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F3.jpg

Fig. 3.

Lifting-line model representing near wake trailed from a rotor blade[15]

전산유체해석을 통해 구한 와류강도를 이용해 앞서 2.1에서 서술한 방법대로 다음 후류의 위치를 구한다. 새로 구한 후류 요소의 위치와 강도를 이용해, 모든 와류 요소가 전산유체해석 경계면에 미치는 유도 속도를 계산한다. 이 유도 속도로 인해 경계면 유동의 밀도, 속도, 에너지가 변한다. 새로 계산된 경계 조건을 이용해 전산유체해석을 하고, 이를 반복한다. 방위각마다 CFD와 자유 후류 기법의 연계가 이뤄지도록 강성 결합을 진행한다. 연계 해석이 진행되는 경계면의 경계 조건으로 유도 속도가 전산유체해석에 영향을 미치는데 본 연구에서는 베르누이 경계 조건, 전엔탈피 고정 경계 조건을 사용하였다. 이에 대한 자세한 설명은 아래 2.2.1와 2.2.2에서 서술된다.

2.2.1 베르누이 방정식 경계 조건

로터가 저속으로 회전하는 회전익기의 블레이드는 비압축성 영역에서 비행한다고 가정할 수 있다. 따라서 본 연구에서는 베르누이 방정식을 사용한 경계 조건을 개발하였다. 식 (5)는 베르누이 방정식으로 첫 번째 항은 정압, 두 번째 항은 동압으로 전압이 일정하다는 가정이 성립된다. 경계면에서 등엔트로피 관계를 만족한다고 가정하였고, 식 (6)는 등엔트로피 관계식이다. 이때 $u, v, w$ 는 자유 후류 기법으로 발생한 후류로 인한 유도 속도가 반영된 속도이다. 각 방향의 자유류의 속도, 압력 밀도( $u_{\infty}, v_{\infty}, w_{\infty} p_{\infty}, ρ_{\infty}$ )와 매시간 간격(time step)마다 유도 속도가 반영되는 $\vec{v}$ 를 식 (5), (6)에 적용해 경계면의 유동 조건 $p, ρ$ 를 구한다. 이를 적용해 전산유체해석을 진행해 익형에서 발생하는 와류의 강도를 계산해서 자유 후류 기법에 반영해 매시간 간격마다 발생하는 후류의 위치를 계산한다. 위 2.2절에서 서술한 방식으로 이 과정을 반복한다.

(5)

p + \frac{1}{2} ρ (u^{2} + v^{2} + w^{2}) = p_{\infty} + \frac{1}{2} ρ_{\infty} (u_{\infty}^{2} + v_{\infty}^{2} + w_{\infty}^{2})

(6)

\frac{p}{ρ^{γ}} = \frac{p_{\infty}}{ρ_{\infty}^{γ}}

또한 유동의 방향을 분석하여 유입류, 유출류에 따라 경계 조건을 달리하였다. 격자가 실제로 회전하기 때문에 격자 셀 내부의 값에 격자 속도를 반영해 경계면의 속도를 구하여 유동의 방향을 결정하였다. 만약 유출류인 경우에는 유도 속도가 반영된 경계면의 속도는 고스트 셀과 내부 셀의 평균이라 설정해 고스트 셀의 속도를 계산하였고, 이를 식 (5), (6)에 적용해 밀도, 에너지를 구하였다. 유입류인 경우에는 리만 불변량(Riemann invariants)[16]을 이용해 경계면의 속도를 구하고 식 (5), (6)을 적용해 밀도, 에너지를 계산한다.

2.2.2 전엔탈피 고정 경계 조건

고속으로 로터가 회전하거나 고속 전진 비행에서는 블레이드는 압축성 유동 영역에서 비행한다. 고속으로 회전하는 로터와 고속 비행을 요구하는 복합형 회전익기의 경우에는 베르누이 방정식이 성립하지 않는 영역에서 비행한다. 따라서 본 연구에서는 베르누이 방정식을 사용하지 않고 전엔탈피가 변하지 않는 경계 조건을 개발하였다. 전엔탈피 ( $H_{T}$ )는 전산유체해석에서 계산된 셀 내부의 압력, 밀도, 속도로 $(p, u, v, w)$ 구성된 식 (7)으로 정의된다.

(7)

H_{t} = \frac{p}{ρ} (\frac{γ}{γ - 1}) + \frac{1}{2} (u^{2} + v^{2} + w^{2})

(8)

\frac{T_{b}}{T_{t}} = 1 + \frac{γ - 1}{2} M_{b}^{2}

(9)

\frac{p_{b}}{p_{t}} = {(1 + \frac{γ - 1}{2} M_{b}^{2})}^{\frac{γ}{γ - 1}}

(10)

\frac{ρ_{b}}{ρ_{t}} = {(1 + \frac{γ - 1}{2} M_{b}^{2})}^{\frac{1}{γ - 1}}

전엔탈피 고정 경계 조건도 유동의 방향을 분석하여 유입류, 유출류에 따라 경계 조건을 달리하였다. 유동의 방향을 결정하기 위해 격자 셀 내부의 값에 격자 속도를 반영해 경계면의 속도를 구하였다. 만약 유출류인 경우에는 유도 속도가 반영된 경계면의 속도는 고스트 셀과 내부 셀의 평균이라 설정해 고스트 셀의 속도를 구했고, 식(7), (8), (9), (10)을 이용해 밀도, 에너지를 계산하였다. 유입류인 경우에는 리만 불변량을 이용해 경계면의 속도를 구한뒤 식(7), (8), (9), (10)를 적용해 경계면의 유동 조건 $(p_{b}, ρ_{b})$ 을 구하였다.

2.3 해석 모델 및 격자

본 연구에서 개발한 하이브리드 기법을 단일 로터에 대해서 검증하고, 경계 조건이 미치는 영향을 알아보기 위해 대표적인 단일 로터 Caradonna&Tung 로터[17]를 사용하였다. 코드 길이 대비 반지름이 6인 비틀림각과 테이퍼비가 없는 직사각형 블레이드이다. 익형은 콜렉티브 피치각 8°의 NACA 0012로 구성되어있다.

뒷전 방향으로 해석영역 크기는 일정하며 아래 흐름 방향으로 해석영역 크기를 달리한 3개의 격자를 사용하였다. 제자리 비행에서는 전진 비행에 비해 후류가 아래 흐름 방향으로 이동하기 때문에 Fig. 4와 같이 격자를 구성하였다. 근방 후류(near-wake)를 포함하기 위해 뒷전 방향으로는 후류의 30° 영역으로(코드 길이의 3.5배) 격자를 생성하였다. 또한 아래 흐름 방향으로는 기존에 수행된 실험 결과와 전산유체해석 결과를 기준으로 가장 작은 격자도 후류 요소가 1개 이상 포함되도록 하였다. 따라서 가장 작은 격자인 grid 1의 아래 흐름 방향으로 블레이드부터 원방 격자 경계면까지 거리는 코드 길이의 1배이다. 해석영역의 크기의 영향을 알아보기 위해 아래 흐름 방향으로 코드길이 1배, 1.5배, 1.75배인 격자로 영향도를 분석하였다. 블레이드 윗면 아랫면 격자점 개수는 코드 방향 121개, 스팬 방향 35개다. Gird 1, 2, 3은 아래 흐름 방향으로 코드 길이의 1배, 1.5배, 1.75배이고 격자점 개수는 49개, 55개, 59개이다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F4.jpg

Fig. 4.

Single rotor grids for hybrid method

2.4 해석 조건

공간 이산화를 위해 격자 중심의 유한 체적법(FVM, Finite Volume Method)를 사용하였고, 2차 정확도의 시간 전진을 위해 이중 시간 전진(dual-time stepping) DADI(Diagonalized Alternating Direction Implicit)[18] 기법을 사용하였다.

본 연구에서 전산유체해석영역은 블레이드 해석을 위한 격자로만 구성된다. 즉, 유동장 해석을 위한 추가적인 포텐셜 영역을 구성하지 않았고, 중첩 격자기법을 사용하지 않아 배경 격자가 없다. 블레이드별 격자 영역 해석만으로도 충분히 타당한 결과를 얻을 수 있도록 CFD/자유 후류를 강성 결합한 하이브리드 기법을 구현하였다. 전산유체해석은 압축성 유동 해석이 가능한 건국대학교 in-house 해석자 KFLOW[19]를 사용하였으며 일반적인 로터 공력 문제에서 검증되고 활용된 해석자이다. 하이브리드 기법 검증을 위해 단독으로 진행한 단일 로터 전산유체해석은 비점성 조건에서 Roe 공간 차분법과 2차 정확도 Minmod 제한자를 사용하였다. 본 논문에서 개발한 하이브리드 기법으로 계산 시 같은 공간 차분법을 적용하였다. 초기 와류 반지름은 코드 길이의 1%로 설정하였다. 자유 후류 기법과 전산유체해석 연계가 이뤄지는 블레이드 근방 경계 조건은 2.2에서 설명한 조건을 사용하며, 블레이드는 비점성 벽으로 설정하였다. 단일 로터는 블레이드 끝단 마하수는 0.439, 0.877로 비행할 때 5° 간격으로 7바퀴 회전하도록 설정하였다.

3. 해석 결과

본 기법의 정확한 검증 및 자세한 분석에 앞서 하이브리드 기법 해석으로 나타나는 후류 거동 모습을 확인하였다. Fig. 5는 전엔탈피 고정 경계 조건을 사용한 하이브리드 기법으로 얻은 후류의 결과다. Grid 1을 5°간격으로 7바퀴 회전한 결과로, 선으로 표현된 후류는 하이브리드 기법의 자유 후류 모델로 구한 후류의 모습이다. Iso-surface는 하이브리드 기법에서 전산유체해석 결과에서 얻은 것이며, 자유류 엔탈피 대비 엔탈피가 1.001일 때 결과다. 즉, 하이브리드 기법을 사용하면 자유 후류 모델로만 후류를 예측하는 것이 아니라 전산유체해석 영역 내에서도 후류의 거동을 확인할 수 있다. 전산유체해석 영역 내에서 예측된 후류와 자유 후류 모델로 예측된 후류의 위치가 유사하다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F5.jpg

Fig. 5.

Result of wake configuration with hybrid method

3.1 저속 회전 로터

3.1.1 추력 계수 변화

계산 시간(방위각)에 따라 총 추력 계수의 변화를 확인해 수렴 여부를 확인하였다. 이는 Fig. 6에서 확인 할 수 있으며, 계산 시간에 따른 총 추력 계수 변화를 나타낸 것으로 grid 1(긴 파선), grid 2(실선), grid 3(짧은 파선) 격자마다 파란색 선은 베르누이 경계 조건, 빨간색 선은 전엔탈피 고정 경계 조건을 사용한 결과다. 이를 통해 2.3에서 제시한 해석 조건으로 5도 간격으로 7바퀴 회전하여 방위각 2,520°에 도달하였을 때, 끝단 마하수 0.439인 경우 논문에서 제시한 모든 해석 문제에 대해 충분히 수렴한 결과를 얻을 수 있다고 판단된다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F6.jpg

Fig. 6.

Convergence of total thrust coefficient at tip Mach 0.439

추력 계수 비교를 통해 하이브리드 기법의 정확도를 분석하였다. Table 1은 각 계산에서 구한 총 추력 계수이다. 하이브리드 기법으로 계산한 총 추력 계수의 평균값을 도출해 전산유체해석 결과와 실험값과 비교하면 평균 오차가 5% 미만이므로 우수한 성능을 보이는 해석 기법이라 판단된다. 정확한 분석을 위해 반경 방향으로 추력 계수(sectional thrust coefficient)를 분석하였고 이것은 Fig. 7에 도시하였다. 파란색 선은 베르누이 방정식 경계 조건, 빨간색 선은 일정한 전엔탈피 경계 조건 결과로 긴파선은 grid 1, 실선은 grid 2, 짧은 파선은 gird 3에 대한 것이며 검은색 점 은 실험 결과[17]다. 초록색 긴 파선은 전산유체해석만을 진행한 결과다. 이를 통해 하이브리드 기법으로 해석한 반지름에 따른 추력 분포의 경향성은 실험 및 전산유체해석 결과와 유사함을 알 수 있다. 또한 이것을 통해 계산의 정확도는 원방 경계면과 블레이드와의 격자의 간격에 영향을 받음을 알 수 있다. 일반적인 아음속 유동의 전산유체해석은 물체로부터 원방 격자까지의 거리가 과하게 좁으면 물체 근방에 원방 경계 조건으로 인한 영향이 작용하기 때문에 오차가 증가한다. 본 논문에서도 격자가 가장 작은 grid 1이 grid 2를 사용한 결과에 비해 정확도가 감소한다. 그러나 본 연구에서는 원방 경계면에서 자유 후류 모델과 연계 해석이 수행된다. 특히 모든 격자 셀에 유도 속도가 반영되지 않고, 경계면에 유도 속도를 반영하는 기법을 사용하므로 원방 격자의 거리가 증가하면 자유 후류의 후류 요소로 인한 유도 속도가 정확히 반영되지 않는다. 따라서 격자가 큰 grid 3의 정확도도 grid 2에 비해 낮아진다. 또한 본 연구에서 개발한 하이브리드 기법에서는 로터의 끝단 후류만을 고려한다. 따라서 추후 하이브리드 기법을 사용한 연구 진행 시 해석영역 크기에 대한 분석도 선행돼야 한다. 경계 조건에 따른 영향은 반경 50% 부근에서는 거의 미미하지만 블레이드 끝단 방향으로 경계 조건에 따른 해석 결과 차이가 증가한다. 해당 결과가 끝단 마하수 0.439인 경우이므로 반경 50% 부근은 비압축성 유동장 영역이고, 블레이드 끝단 부근은 압축성 유동장 영역이다. 따라서 로터 운용 조건에 따라 베르누이 경계 조건을 사용하는 것보다 전엔탈피 고정 경계 조건을 사용해야 하며, 고속으로 회전하는 로터 해석 결과(3.2절)에서 추가적인 결과를 확인할 수 있다.

Table 1.

Comparison of total thrust coefficient at tip Mach 0.439

Hybrid method					CFD	Exp
Boundary condition	Grid 1	Grid 2	Gird 3	Mean value	CFD	Exp
Bernoulli’s equation	0.0041	0.0046	0.0049	0.0045	0.0048	0.0046
Total enthalpy	0.0041	0.0046	0.0049	0.0045	0.0048	0.0046

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F7.jpg

Fig. 7.

Comparison of sectional thrust coefficient on radial direction at tip Mach 0.439

Fig. 8(a-c)는 블레이드 반지름의 50%, 80%, 96% 지점에서 표면 압력 계수를 도시한 결과다. 파란색은 베르누이 방정식 경계 조건, 빨간색은 전엔탈피 고정 경계 조건을 사용한 결과를 의미하며, 긴 파선은 grid 1, 실선은 grid 2, 짧은 파선은 gird 3을 사용한 결과이다. 초록색 파선은 전산유체해석결과, 검은색 점은 실험 값[17]이다. 하이브리드기법의 와류 요소와 경계면이 가까울수록 후류로 인한 더 강한 유도속도가(식 (2)) 전산유체해석 경계 조건에 반영된다. 따라서 하이브리드 기법에서 발생한 끝단 후류로 인해 발생한 유도 속도의 영향이 더 많이 고려된 블레이드 끝단 방향의 정확도가 높고, 윗면보다 아랫면의 압력 계수 분포 정확도가 높다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F8.jpg

Fig. 8.

Comparison of pressure coefficient on boundary conditions and grids at tip Mach 0.439

3.1.2 끝단 후류 분석

후류의 자세한 위치는 Fig. 9에서 분석하였다. Fig. 9(a)는 후류가 생성됨에 따라 반지름 방향으로의 끝단 후류 위치 변화이고, Fig. 9(b)는 아래 흐름 방향으로의 끝단 후류의 위치를 나타낸 것이다. 초록색 2점 쇄선은 KFLOW를 사용해 전산유체해석 단독으로 해석한 결과다. 파란색은 베르누이 방정식 경계 조건, 빨간색은 전엔탈피 고정 경계 조건을 사용한 결과를 의미하며, 각각의 긴 파선은 grid 1, 실선은 grid 2, 짧은 파선은 gird 3을 사용한 결과다. 검은색 점은 실험값[17]을 의미한다. 저속으로 회전 로터에서 자유 후류 모델로 생성된 후류에 경계 조건과 격자가 미치는 영향이 미미하다. 로터 허브로 수축하는 경향과 아래 흐름으로 이동하며 생성되는 후류의 경향성은 실험값과 전산유체해석 결과와 유사하나, 후류의 진동을 보인다. 후류의 진동에 미치는 변수를 확인하기 위해 시간 간격과 초기 후류 반지름에 따른 후류의 변화를 Fig. 10을 통해 확인하였다. Grid 2 기준으로, 그래프에 명시된 Rc는 코드길이 대비 초기 와류 반지름, 5deg, 10deg는 시간간격을 의미한다. 와류 초기 반지름이 미치는 영향은 미미하나 시간 간격마다 후류 요소를 선으로 가정하기 때문에 시간 간격을 증가하면 진동을 완화시킬 수 있다. 와류 초기 반지름 크기와 무관하게 총 추력 계수는 시간 간격이 5°이면 0.0046, 시간 간격이 10°이면 0.0045로 시간 간격이 커질수록 오차는 다소 증가한다. 다만, 시간 간격이 증가하면 자유후류로 생성되는 후류 요소가 감소해 계산량이 감소한다. 따라서 시간 간격이 5°에서 약 1시간 30분 소요되던 계산 시간이 시간 간격을 10°로 증가했을 때 계산 시간은 약 30분으로 감소했다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F9.jpg

Fig. 9.

Comparison of wake geometry on wake age

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F10.jpg

Fig. 10.

Comparison of wake geometry on wake age to analyze the wake vibration

3.1.3 유동장 분석

Fig. 11은 가장 높은 정확도로 추력 계수를 예측한(Table 1) grid 2의 해석 결과로, Fig. 11(a)는 베르누이 경계 조건, Fig. 11(b)는 전엔탈피 고정 경계 조건을 사용한 결과이다. 이것을 전산유체해석 결과(Fig. 11(c))와 비교해 물리적으로 타당한 유동장이 발생하는지 확인하였다. 하이브리드 기법을 사용한 결과와 전산유체해석만을 진행했을 때 발생하는 끝단 와류와 유사한 위치에서 끝단 와류가 발생하며 와류 강도도 유사하다. 따라서 비슷한 끝단 와류가 발생하기 때문에 3.1.1절에서 설명한 추력 분포도, 3.1.2절에서 명시한 자유후류모델로 발생한 후류의 위치가 전산유체해석과 유사한 결과를 보임을 설명할 수 있다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F11.jpg

Fig. 11.

Comparison of vorticity distribution between from hybrid Method and CFD at tip Mach 0.439

3.2 고속 회전 로터

3.1절에서는 끝단 마하수 0.439로 상대적으로 저속 회전 로터에 대한 문제로 베르누이 경계 조건을 사용하더라도 로터의 허브 근방 추력 성능에는 영향이 적은 운용 조건이다. 경계 조건에 따른 명확한 차이를 확인하기 위해 블레이드 끝단 속도 0.877로 회전하는 문제에 대해 분석하였다.

계산 결과의 수렴성을 확인하기 위해 계산 시간에 따른 로터의 총 추력 계수 변화를 확인하였다. 이것은 Fig. 12에서 확인 가능하다. Grid 1(긴 파선), grid 2(실선), grid 3(짧은 파선) 격자마다 전엔탈피 고정 경계 조건을 사용한 결과다. 이를 통해 2.3에서 제시한 해석 조건으로 5도 간격으로 7바퀴 회전하여 방위각 2,520°에 도달하였을 때, 끝단 마하수 0.877인 경우 충분히 수렴한 결과를 얻을 수 있다고 판단된다. 다만, 베르누이 경계 조건을 사용하면 발산해 값을 도출할 수 없다. 본 문제에서는 끝단 속도 0.877로 회전해 블레이드에 발생하는 유동이 압축성 유동이다. 따라서 비압축성/비회전성 유동에서만 성립하는 베르누이 방정식을 이용한 경계 조건을 사용하면 비물리적인 유동이 발생해 계산이 발산하므로 정확한 결과를 구할 수 없다.

Table 2는 격자와 경계 조건에 따른 로터 총 추력 계수 결과와 실험 결과[17]를 비교한 것이다. 전엔탈피가 일정한 경계 조건은 비압축성/압축성, 비회전성/회전성 유동에서 모두 성립되므로, 해당 경계 조건을 사용했을 때는 실험값, 전산유체해석 결과와 5% 미만의 오차로 타당한 추력을 구할 수 있다. 그러므로 하이브리드 기법을 사용해 다양한 운용 조건에서 물리적으로 정확한 로터의 공력을 해석하기 위해서는 전산유체해석과 자유 후류 모델이 연계되는 원방 경계면에 전엔탈피가 일정하도록 설정해야 한다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F12.jpg

Fig. 12.

Convergence of total thrust coefficient

Table 2.

Comparison of total thrust coefficient at tip mach 0.877

Hybrid method					CFD	Exp
Boundary condition	Grid 1	Grid 2	Gird 3	Mean value	CFD	Exp
Bernoulli’s equation	Divergence	Divergence	Divergence	-	0.0046	0.0047
Total enthalpy	0.0040	0.0048	0.0054	0.0047	0.0046	0.0047

그리고 반경에 따른 추력 계수 변화를 Fig. 13를 통해 분석하였고 파선은 grid 1, 실선은 grid 2, 짧은 파선은 gird 3, 초록색 긴파선은 전산유체해석 결과, 검은색 점은 실험 결과[17]다. 저속 유동에서 운용하는 것과 유사하게 블레이드와 원방 경계면의 거리에 영향을 받는다. 일반적인 아음속 유동의 전산유체해석은 물체 벽면과 원방 격자의 거리가 지나치게 작으면, 물체 근방에 원방에서 발생하는 유동의 영향이 작용하여 오차가 증가한다. 본 논문에서도 grid 1을 사용한 결과의 정확도가 grid 2보다 낮다. 그러나 본 연구에서는 경계면에 유도 속도를 반영하기 때문에 원방 격자의 거리가 증가하면 끝단 후류로 인해 경계면에 유도되는 속도가 정확히 반영되지 않는다. 이에 따라, grid 3의 정확도도 grid 2에 비해 낮다. 또한, 저속으로 회전하는 경우보다 고속으로 회전하는 경우 하이브리드 기법에 계산 영역의 영향성이 증가한다. 회전 속도가 증가할수록 경계면에서 발생하는 유동이 블레이드 주변 유동에 미치는 영향이 증가하므로 해석 영역의 크기가 중요하다. 또한 원방 경계면에서 전엔탈피가 고정되기 때문에 압축성 유동, 충격파가 발생하는 문제 등에서는 해석 영역의 크기가 중요하다. 문제에 따라 로터 근방에서 발생하는 유동장 및 후류의 거동이 모두 다르기 때문에 본 논문에서 진행한 연구만으로 하이브리드 기법에서 적합한 해석 영역의 크기를 일반화하여 제시하기 어렵다. 따라서 하이브리드 기법으로 로터의 공력 성능을 분석할 때, 반드시 해석 영역의 크기에 따른 결과를 분석해 격자에 따른 수렴성을 확인해야 한다.

추가적으로 Fig. 14을 통해서 반경 50%, 80%, 96% 지점에서의 표면압력을 분석하였다. 실험 결과[17]와 비교했을 때 충격파가 발생하지 않는 반경 50% 지점의 표면 압력 계수뿐만 아니라 충격파가 발생하는 반경 80%, 96% 지점의 표면 압력 분포도 상당히 일치하였다. 위에서 서술한 이유로 격자 원방의 위치에 따라 충격파 발생 위치가 민감하게 변하지만 grid 2를 사용한 결과가 전반적으로 압력 분포의 경향을 잘 따라가고 있음을 알 수 있다. 결론적으로 고속의 압축성 유동에서 자유 후류 모델 연계 해석 기법을 사용할 때 본 연구에서 개발한 전엔탈피 고정 경계 조건이 효과적으로 적용될 수 있음을 확인하였다.

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F13.jpg

Fig. 13.

Convergence of total thrust coefficient

https://cdn.apub.kr/journalsite/sites/kscfe/2024-029-04/N0500290409/images/jkscfe_2024_294_102_F14.jpg

Fig. 14.

Comparison of pressure coefficient on grids at tip Mach 0.877

3.3 하이브리드 기법의 효율성 분석

하이브리드 기법의 효율성을 확인하기 위해 전산유체해석에 요구되는 계산 자원 및 계산 시간을 비교하였다. Table 3에 각 기법 별 격자계 및 격자 개수를 명시하였다. 하이브리드 기법의 격자는 grid 2 기준이고, 전산유체해석 격자는 본 논문에서 제시한 전산유체해석 결과를 도출하기 위한 격자다. I 방향은 블레이드 위/아래 방향이며, j 방향은 블레이드 스팬 방향, k 방향은 익형 코드 방향으로 격자점 수를 나타낸 것이다. 전산유체해석은 로터의 운동을 모사하기 위해 중첩격자기법을 사용한다. 따라서 부격자 계와 배경 격자 계로 구성되고, 총 격자 개수는 25,966,080개다. 반면 하이브리드 기법은 블레이드 격자 계만 구성하면 되므로 4,136,832개의 격자가 필요하다. 즉, 격자수가 약 6배 차이가 발생하며, 이에 따라 계산 자원이 더 많이 요구되며 계산 시간도 증가한다. 정확도를 비교하면, 두 방법 총 추력 계수를 비교했을 때 앞선 3.1, 3.2절에서 확인 할 수 있듯 모두 실험값과 5% 이내다. 하지만 계산 시간은 Table 4에서 확인할 수 있듯, 수렴되기 위해 하이브리드 기법 대비 전산유체해석 시간이 약 9배 더 소요된다. 만약 같은 회전수 조건 7바퀴로 계산한다면 21배 더 소요된다. 따라서 정교한 유동장을 확인하기 위해서는 전산유체해석을 하는 것이 중요하지만, 효율적으로 공력 성능을 얻기 위해서는 하이브리드 기법을 사용하는 것이 유용함을 알 수 있다.

Table 3.

Comparison of grid between CFD and hybrid method

	Grid	The number of cells
CFD (Unsteady)		Subgrid	2 × 2,995,200 (2× (49×201×313)
		Background	19,975,680 (271×273×273)
		Total	25,966,080
Hybrid (Unsteady)		Only blade	2 × 2,068,416 (2× (55×115×337))
Hybrid (Unsteady)		Total	4,136,832

Table 4.

Comparison of wall clock time between CFD and hybrid method

	Revolution	Wall clock time
CFD(Unsteady)	3 rev (7rev)	13.7 h (32.0 h)
Hybrid(Unsteady)	7 rev	1.5 h

4. 결 론

본 논문에서는 전산유체해석기법과 자유 후류 기법을 결합한 하이브리드 기법을 개발하였고, 정보의 연계가 이루어지는 경계면의 경계 조건에 관해 연구하였다. 경계 조건은 압축성 유동을 고려한 베르누이 방정식을 활용한 경계 조건, 전엔탈피가 일정한 경계 조건 2가지를 개발해 단일 로터 제자리 비행 조건에서 비교 분석하였다. 하이브리드 기법으로 구한 끝단 후류의 구조와 총 추력 계수를 실험 결과, 전산유체해석 단독 해석 결과와 비교하였을 때 신뢰할만한 결과를 획득하였다. 하지만 베르누이 방정식은 비압축성 비회전성 유동에서 만족하는 방정식이므로, 이를 압축성 유동 해석의 경계 조건에 적용하면 비물리적인 결과를 얻는다. 반면, 본 연구에서 개발한 전엔탈피가 일정한 조건은 압축성/비압축성, 회전성/비회전성과 관계없이 유동장 영역에서 적용할 수 있고, 이를 압축성 유동 해석방법에 적용하면 물리적으로 타당한 결과를 얻을 수 있다. 기존의 베르누이 관계식에 기반한 경계 조건을 보완하여 제시된 전엔탈피 경계 조건이 더 물리적으로 타당한 결과를 제시할 수 있음을 확인하였다.

Acknowledgements

이 논문은 2021년도 정부(과학기술정보통신부)의 재원으로 한국연구재단의 지원을 받아 수행된 연구입니다(No.2021R1A5A1031868).

References

2016, Lee, J.B., "Numerical Study on Time Efficiency of Rotor Performance Analysis Using the CFD-Freewake Coupling Method," Ph. D. Thesis, Pusan University.

2008, Inada, Y., Yang, C., Iwanaga, N. and Aoyama, T., "Efficient Prediction of Helicopter BVI Noise under Different Conditions of Wake and Blade Deformation," Trans. Jpn. Soc. Aeronaut. Space Sci., Vol.51, No.173, pp.193-202.

10.2322/tjsass.51.193

2006, Sitaraman, J. and Baeder, J.D., "Field Velocity Approach and Geometric Conservation Law for Unsteady Flow Simulations," AIAA J., Vol.44, No.9, pp.2084-2094.

10.2514/1.5836

2010, Wie, S., "Development of Tim-Marching Free-Wake and Navier-Stokes Solver Coupling Method for Rotor Performance," Ph. D. Thesis, KAIST.

2010, Rajmohan, N., "Application of Hybrid Methodology to Rotors in Steady and Maneuvering Flight," Ph. D. Thesis, Georgia Institute of Technology.

2010, Min, B. and Sankar, L.N., "Hybrid Navier-Stokes/Free-Wake Method for Modeling Blade-Vortex Interactions," J. Aircr., Vol.47, No.3, pp.975-982.

10.2514/1.46550

1997, Berkman, M.E., Sankar, L.N., Berezin, C.R. and Torok, M.S., "Navier-Stokes/full potential/free-wake method for rotor flows," J. Aircr., Vol.34, No.5, pp.635-640.

10.2514/2.2240

2002, Leishman, J.G., Bhagwat, M.J. and Bagai, A., "Free-Vortex Filament Methods for The Analysis of Helicopter Rotor Wakes," J. Aircr., Vol.39, No.5, pp.759-775.

10.2514/2.3022

2002, Yang, Z., Sankar, L.N., Smith, M.J. and Bauchau, O., "Recent Improvements to a Hybrid Method for Rotors in Forward Flight," J. Aircr., Vol.39, No.5, pp.804-812.

10.2514/2.3000

2014, Farrokhfal, H. and Pishevar, A.R., "A New Coupled Free Wake-CFD Method for Calculation of Helicopter Rotor Flow-Field in Hover," J. Aerosp. Technol. Manag., Vol.6, No.2, pp.129-147.

10.5028/jatm.v6i2.366

2010, Egolf, T.A., Rajmohan, N., Reed, E. and Sankar, L.N., "A Hybrid CFD Method for Coaxial Rotor Performance Prediction in Forward Flight," AHS Specialist's Conference on Aeromechnics.

1995, Bagai, A. and Leishman, J., "Rotor Free-Wake Modeling Using a Pseudo-implicit Relaxation Algorithm," J. Aircr., Vol.32, No.6, pp.1276-1285.

10.2514/3.46875

1994, Bagai, A. and Leishman, J.G., "Rotor Free-Wake Modeling Using a Pseudo-Implicit Technique - Including Comparisons with Experimental Data," Proceedings of the 50th Annual Forum of the American Helicopter Society, Washington D.C., USA.

1965, Squire, H., "The Growth of a Vortex in Turbulent Flow," Aeronautical Quarterly, Vol.16, pp.302-306.

10.1017/S0001925900003516

2006, Leishman, J.G., Principle of Helicopter Aerdoynamics, Cambridge University.

2011, Carlson, J., "Inflow/Ouflow Boundary Conditions with Application to FUN3D", NASA Langley Research Center, NASA/TM-2011-217181.

1981, Caradonna, F.X. and Tung, C., " Experimental Analytical Studies of a Model Helicopter Rotor in Hover," NASA Technical Memorandum 81232.

2002, Park, S.H., Kim, Y. and Kwon, J.H., "Prediction of Dynamic Damping Coefficients using Unsteady Dual-Time Stepping Method," 40th AIAA Aerospace Sciences Meeting&Exhibition, pp.715-723.

2008, Kim, J.W., Park, S.H., Ko, J.H., Park, I., Jung, S.N., Yu, Y.H., Kim, E. and Kwon, J.H., "Hovering Rotor Computation including Aero-Elastic Effects with an Overlapped Grid Navier-Stokes Solver," 26th AIAA Applied Aerodynamics Conference, pp.6729-6741.

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

Preview

STUDY ON BOUNDARY CONDITIONS FOR CFD/FREE WAKE COUPLED HYBRID METHOD

ABSTRACT

MAIN

(1)

(2)

(3)

Fig. 1.

Principle of hybrid method through the CFD boundary interface

(4)

Fig. 2.

Procedure of hybrid method with tight coupling method

Fig. 3.

Lifting-line model representing near wake trailed from a rotor blade[15]

(5)

(6)

(7)

(8)

(9)

(10)

Fig. 4.

Single rotor grids for hybrid method

Fig. 5.

Result of wake configuration with hybrid method

Fig. 6.

Convergence of total thrust coefficient at tip Mach 0.439

Table 1.

Comparison of total thrust coefficient at tip Mach 0.439

Fig. 7.

Comparison of sectional thrust coefficient on radial direction at tip Mach 0.439

Fig. 8.

Comparison of pressure coefficient on boundary conditions and grids at tip Mach 0.439

Fig. 9.

Comparison of wake geometry on wake age

Fig. 10.

Comparison of wake geometry on wake age to analyze the wake vibration

Fig. 11.

Comparison of vorticity distribution between from hybrid Method and CFD at tip Mach 0.439

Fig. 12.

Convergence of total thrust coefficient

Table 2.

Comparison of total thrust coefficient at tip mach 0.877

Fig. 13.

Convergence of total thrust coefficient

Fig. 14.

Comparison of pressure coefficient on grids at tip Mach 0.877

Table 3.

Comparison of grid between CFD and hybrid method

Table 4.

Comparison of wall clock time between CFD and hybrid method

Acknowledgements

References