OpenFOAM 기반 밀도 기반 해석자를 활용한 단열 막냉각 효율 비교 분석

2.1 막냉각 홀 형상 및 해석 영역

막냉각 성능은 냉각 홀 형상과 배치에 의해 결정된다. 특히, 홀 형상은 표면 막냉각 효율 분포에 직접적인 영향을 미치는 주요 인자로 작용한다. 본 연구에서는 Kim⁽⁶⁾ 등이 실험한 7-7-7 팬 형상(fan-shaped) 막냉각 홀을 반영하였다. 해당 홀은 유동방향과 횡방향으로 각각 7°의 확장각을 갖는 형상이며, 해당 홀의 형상 및 제원에 대한 자세한 정보는 Table 1과 Fig. 3에 제시되어 있다.

Table 1

7-7-7 Fan-shaped hole geometry

Fig. 3

Configurations of 7-7-7 fan-shaped hole geometry(6)

해석영역의 경계조건을 Fig. 4와 같이 설정하였다. 막냉각 홀의 직경(D)을 기준으로 높이는 20D, 주유동이 통과하는 방향으로 –30D부터 40D까지 총 70D, 폭은 6D로 설정하였다. 이는 막냉각 홀에서 분사되는 냉각 유체가 하류에서 충분히 발달할 수 있도록 하기 위한 조건이다. 막냉각 홀이 위치한 평판 표면은 열 손실을 방지하고 막냉각 효과를 명확히 평가하기 위해 수치해석에서 단열 조건으로 지정하였다. 주유동과 냉각 유체의 난류 강도는 Kim⁽⁶⁾ 등의 실험 조건과 동일하게 입구 조건으로 적용하였으며, 세부 유동 조건은 Table 2에 정리하였다. 이러한 조건은 실제 가스터빈 운용 환경과는 차이가 있으나, 본 연구는 막냉각 유체 특성의 경향을 파악하고 수치해석 모델을 검증하는 데 목적을 두고 있다.

Fig. 4

Boundary conditions

Table 2

Flow inlet conditions

주유동과 냉각 유체의 분사비(Blowing ratio, M)는 식 (1)과 같다.

식 (1)에서 (ρU)_coolant는 냉각 유체이고 (ρU)_Main은 주유동이다.

분사비(M)는 냉각 유체의 분사 강도와 벽면 부착 특성을 결정하는 주요 인자로 분사비가 달라짐에 따라 막냉각 효율이 변화한다. 이를 비교하기 위해 냉각 유체 입구의 경계 조건은 질량 유량을 기준으로 분사비(M) 0.5, 1.0, 2.0의 세 가지 조건을 설정하였고, 분사비에 따른 유량 조건은 Table 3과 같다.

Table 3

Mass flow rate settings for blowing ratio cases

본 연구에서의 수치해석을 위해 적용한 격자는 Fig.5에 제시되어 있으며, Ansys Mesh를 기반으로 생성하였고, 압력 기반의 막냉각 홀 예측에 사용된 바가 있다⁽⁷⁾. 본 연구에서는 적용한 해석격자의 영향을 제거하기 위해 선행연구⁽⁷⁾와 동일한 격자를 활용하였다. 이에 적용한 전체 격자수는 약 1.9 × 10⁷개이다. 해석 정확도 확보를 위해 홀 및 벽면 주변의 격자의 밀집도가 높은 부분(Fine region)에는 홀 직경(D)을 기준으로 0.1D의 격자 크기를 적용하였다. 또한, 벽면 경계층 해석의 정밀도를 확보하기 위해, 벽면 격자의 수는 15개로 설정하였다.

Fig. 5

CFD mesh configuration

평판 막냉각 홀 벽면에서의 유동을 정밀하게 해석하기 위해서는 해석영역에 대해 y⁺ 값이 1 이하가 되도록 격자계를 구성할 필요가 있다. Fig. 6은 분사비(M) 1.0 조건에서의 결과로, 전체 영역에서 y⁺ 값이 1 이하가 유지됨을 보여주며, 본 연구의 격자 해상도가 벽면 열유동 해석에 충분함을 확인하였다.

Fig. 6

y⁺ Value at M=1.0

2.2 해석 기법

본 연구에서는 가스터빈 내에 밀도차가 큰 막냉각 현상을 정밀하게 해석하기 위해 dbturb를 활용한 수치해석을 수행하였다. dbturb는 OpenFOAM 기반의 공개된 밀도 기반 해석자 개발 라이브러리인 HiSA(High-Speed Aerodynamics)⁽⁸⁾을 보다 효과적으로 터빈해석에 적용하기 위하여 수정한 해석자이다. 기존 HiSA⁽⁸⁾는 충격파나 압축성 효과가 두드러지는 고마하수 영역에서는 우수한 성능을 보이나, 밀도차가 크고 유속이 낮은 저마하수 영역에서는 해석 안정성이 저하되는 경향을 보였다. 이러한 특성으로 인해 평판 막냉각 홀과 같이 복합 열유동이 발생하는 영역에서는 유동의 세부 거동을 안정적으로 예측하기 어렵다. 또한 OpenFOAM 내에서 일반적으로 사용되는 압력 기반 해석자는 고속 압축성 유동에서 수치적 불안정성이 발생한다. 한편, 밀도 기반의 segregated solver는 고속 압축성 유동조건에서 높은 해석 정확도를 제공하더라도 CFL(Courant Friedrichs Lewy) 수 제한으로 인해 계산 효율이 떨어지는 문제와 기존 OpenFOAM 라이브러리와의 연계가 어렵다는 문제가 있다. 이러한 한계를 개선하기 위해, 본 연구에서는 HiSA⁽⁸⁾를 기반으로 높은 CFL 수 적용을 위한 Coupled solver를 개발하고 고속 압축성 유동뿐만 아니라 저마하수 영역에서도 수치적 안정성과 해석 정확도를 확보할 수 있도록 Low-Mach 보정 기법을 적용하였다. 이를 통해 기존 압력 기반 해석자의 불안정성을 개선하고 밀도 기반 해석자의 효율 저하 문제를 동시에 보완할 수 있도록 개선하였다.

압력 기반 RANS 해석자는 압력(p), 속도(U), 온도(T)를 주요 변수로 사용하며 OpenFOAM의 기존 라이브러리와 높은 호환성을 가진다는 장점이 있다. 그러나 이러한 접근법은 고속 유동이나 급격한 압력 변화가 존재하는 영역에서 충격파 및 불연속의 정확한 처리에 한계를 가지며, 수치적 불안정성 또는 정확도 저하를 유발할 수 있다. 이를 보완하기 위해 dbturb에서는 밀도(ρ), 운동량(ρU), 에너지(ρE)를 직접 해석 변수로 채택하는 밀도 기반 해석 방식을 적용하였다. 비정상 압축성 점성 유동의 지배방정식인 연속 방정식, 운동량 방정식, 에너지 방정식은 다음의 식과 같다.

여기서 u_i는 직교 좌표계에서의 속도, $$ E=e+{\left|\overrightarrow{u}\right|}^2/2$$는 총 에너지, e는 내부 에너지이다. 연속성, 운동량 및 에너지 방정식을 사용하여 비정상 압축성 점성 흐름을 모델링하며, 이 방정식들은 완벽한 기체 상태 방정식을 사용하여 닫히고, 점성 응력 텐서의 구성 요소를 포함한다.

해석의 정확성을 결정짓는 핵심 요소 중 하나는 flux 분할 기법이며, dbturb에서는 기존에 사용되던 HLLC⁽⁹⁾ 기법을 중심으로 수치 진동을 억제하고, 충격파 및 급격한 유동 변화 영역에서 해석의 안정성을 확보하기 위해 Minmod limiter를 적용하였다. 해당 limiter는 좌우 인접 셀 간의 경사 기울기를 비교하여 수치적으로 가장 안정적인 값을 선택하는 보수적 방법으로, 고차 기법 적용 시 발생할 수 있는 불안정을 효과적으로 제어할 수 있는 장점이 있다. 반면에, HLLC 기법은 충격파 및 접촉면을 명확하게 포착하는 데 효과적인 특성을 가지나, 저마하수 유동 영역에서는 수치적 확산 증가와 압력 예측의 부정확성 등 한계가 존재한다. 이러한 문제를 개선하기 위해, 본 연구에서는 Fleischmann⁽⁹⁾ 등이 제안한 shock-stable 및 Low-Mach 보정 기법을 도입하였으며, 이 기법은 기존 HLLC 기법의 파동 속도 S_l에 감쇠 함수 ϕ를 곱하여 저마하수 안에서의 대류 소산과 음향 소산의 균형을 전체 소산의 감소를 통해 이루어지도록 하여 해석의 안정성과 예측 정확도를 향상시키도록 구성하였다. 새로운 파동 속도 $$ S_l^{HLLC-LM}$$은 식 (5)와 같다.

감쇠 함수 ϕ는 마하수 기반의 함수로, 저마하수 영역에서 수치 소산을 효과적으로 억제하고, 고마하수 영역에서는 기존 HLLC의 성능을 유지할 수 있도록 설계되었다. 감쇠 함수 ϕ는 국소 마하수 Ma_local에 따르고, 식 (6)과 같다.

여기서 국소 마하수 Ma_local는 셀 경계면 양쪽에 위치한 유동 상태를 고려하여 정의되며, 좌측 상태의 속도 u_l와 음속 c_l, 우측 상태의 속도 u_r와 음속 c_r를 기반으로 식 (7)과 같다.

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

수치적 효율성과 안정성을 향상시키기 위해 본 연구에서는 시간 이산화 방식으로 matrix-free implicit solver⁽¹⁰⁾를 채택하였다. 지배 방정식은 시간에 대해 이산화되며, 식 (8)과 같다.

여기서 W는 보존 변수 벡터, R은 지배 방정식의 우변 항, V는 제어 체적을 의미한다. 시간 적분 방식은 계수 β에 따라 결정되며, β=1일 경우는 1차 implicit 방식, β=0은 1차 explicit 방식, β=1/2일 경우는 Crank–Nicolson 방식으로 분류된다. 이 중 explicit 방식은 계산 구조가 단순하고 구현이 용이하다는 장점이 있지만, 수치적 안정성을 위해 CFL (Courant–Friedrichs–Lewy) 조건을 만족해야 하므로 시간간격(∆t)의 제한이 따른다. 이때 시간 간격(∆t)은 식 (9)와 같다.

여기서 λ는 시스템의 고유값이고, ∆x는 격자 길이 척도이며 일반적으로 $$ \mathrm{\Delta }x=\sqrt{\mathrm{\Delta }{x}_i\mathrm{\Delta }{x}_i}$$ 로 계산된다. β=0인 explicit 계산의 경우에는 안정성 확보를 위해 CFL < 1 조건을 만족해야 한다. β=1인 implicit 방식은 안정성 조건에 따른 시간 간격 제한이 다소 완화되기 때문에 상대적으로 큰 ∆t를 사용할 수 있어 전체 계산 효율이 향상될 수 있다. 그러나 implcit 계산의 경우에는 선형 연립 방정식의 반복 풀이가 요구되므로 계산 자원과 메모리 사용량이 증가하는 단점이 존재한다. 이러한 단점을 보완하기 위해 본 연구에서는 matrix-free implicit 방식을 도입하였다. 이 방식은 Jacobian 행렬을 명시적으로 구성하지 않고, 잔차 벡터에 대한 Jacobian 연산만을 반복적으로 수행함으로써, 대규모 시스템에서도 상대적으로 적은 메모리 자원으로 계산을 수행할 수 있다. 이 방법은 특히 비정방형(nonsymmetric) 또는 양의 정부호가 아닌 시스템에서 효과적으로 작동하며, 시간 이산화 항과 선형화된 Jacobian 항을 포함한 좌변 계수 행렬은 식 (10)과 같다.

여기서 좌변의 계수 행렬 A는 시간 이산화 V/∆t로부터의 대각항을 포함하며, 야코비안 행렬 ∂Rⁿ/∂W의 항을 함께 포함한다.

OpenFOAM의 특성상 경계 조건을 구현하기 위한 ghost cell을 사용할 수 없으므로, 밀도 기반 해석자에서 주로 사용되는 ghost cell 기반의 압력 경계 조건 개발이 제한된다. 본 연구에서는 내부 유동장에서 선형화된 관계식을 바탕으로 압력 경계 조건을 구성하였다. 벽면 경계 조건 구현을 Riemann 문제 해결 방식의 특성 관계를 이용하여 내부점 데이터를 외삽하였다. 압력 경계 조건의 경우, zero velocity 조건에서의 적용을 위하여 Liou⁽¹¹⁾가 제안한 선형화된 압력 경계 조건을 사용하였다. 벽면에서의 압력 조건은 식 (11)과 같다.

dp는 압력의 변화, du는 속도의 변화를 의미한다. ρ는 밀도 a는 음속이다. 이를 통해 유도된 경계면 압력 외삽 조건은 식 (12)와 같다.

식 (12)는 경계면에서 압력 p_w를 계산하기 위한 외삽 조건으로, 내부 셀의 특성값을 바탕으로 압력 값을 예측하는 방식이다. 여기서 i는 내부 셀을, w는 벽 또는 경계를 나타낸다. $$ p_i^n,u_i^n,{\rho }_i^n,a_i^n{u}^{n_i}$$는 각각 시간 n에서 내부 셀의 압력, 속도, 밀도, 음속을 의미한다. 이 방식은 경계면에서 압력을 내부 셀의 특성값을 기반으로 외삽하여 계산함으로써, 수치 경계에서 발생할 수 있는 단기 진동(numerical oscillation) 문제를 억제하고, 동시에 계산 안정성과 정확성을 확보할 수 있다.

시간 이산화 기법으로는 RANS 기반의 정상 상태 해석 (steady state)을 위해 dual-time stepping⁽¹²⁾ 기법을 적용하였다. 난류 모델은 k–ω SST⁽¹³⁾ 모델을 사용하였다.

3.1 실험 및 수치해석결과 비교

본 연구에서는 분사비 변화에 따른 단열 막냉각 효율 특성과 유동 구조의 차이를 분석하기 위하여, 동일한 격자계(Fig. 5)를 기반으로 밀도 기반 해석자 dbturb를 이용하여 수치 해석을 수행하였다.

막냉각 효율 예측에 앞서 유동 해석의 적합성을 확인하기 위해 막냉각 홀의 영향이 없는 상류에서 속도 경계층을 실험 결과와 비교하였다. Fig. 7은 Kim⁽⁶⁾ 등이 x/D = -12에서 열선유속계를 이용해 측정한 유동의 경계층 속도 프로파일을 dbturb의 해석 결과와 비교한 결과이다. 이를 통해 본 연구에서 사용된 경계 조건과 해석 조건을 통해 획득한 유동 결과가 실험 결과와 동일한 결과를 보여줌을 확인하였다.

Fig. 7

Velocity profile

Fig. 8은 Kim⁽⁶⁾ 등이 팬 형상 막냉각 홀을 대상으로 압력 감응 페인트(Pressure Sensitive Paint, PSP)를 이용하여 측정한 결과와 dbturb 결과의 단열 막냉각 효율 분포를 분사비 조건별로 제시하였다. 막냉각 홀 벽면의 표면 막냉각 효율의 정의는 식 (13)과 같다.

Fig. 8

Contour comparison of adiabatic film cooling effectiveness from experiment(6) and dbturb at different blowing ratios

T_∞ 는 주유동 온도이며, T_c는 냉각 유체동 온도이다. T_aw는 벽면에서의 단열 상태에서 형성되는 벽면 온도를 의미한다.

Fig. 8에서 dbturb는 초기 영역(x/D ≤ 2.5)에서 실험과 dbturb 간의 평균 막냉각 효율 차이가 비교적 크게 나타났으나, 이후 하류 영역에서는 막냉각 효율이 유사한 감소 추세를 보였다. 이는 dbturb가 냉각 유체의 초기 벽면 부착 및 확산 유동을 과대평가하여 초기에 막 형성 효과를 상대적으로 높게 예측하기 때문으로 해석된다. 또한 막냉각 효율은 홀 중심부에서 낮고 홀 양 끝부분에서 높은 양상을 보인다. 이는 k–ω SST 난류 모델을 적용한 경우, 홀 내부에서 유동이 가속되는 구간과 확장부에서 박리에 따른 감속이 발생하면서 속도 분포가 균일하지 않게 형성되기 때문이다.

Fig. 9은 dbturb 해석자에서 M=2.0일 때 팬 형상 막냉각 홀 내부 z/D=0에서 속도장을 나타낸 것이다. 이러한 유동 특성은 Saumwebere⁽¹⁴⁾ 등의 연구에서도 확인된 바 있다.

Fig. 9

Dbturb velocity vectors and profile at center plane at M=2.0

M=0.5에서는 냉각 유체의 상대 운동량이 작아 주유동과의 혼합이 빠르게 일어나고, 냉각 유체는 벽면에 밀착된 상태에서 단기간 고효율을 유지한다. 그러나 냉각 유체는 빠르게 확산되어 하류로의 냉각 성능은 빠르게 감소하는 양상을 Fig. 8(a)와 8(b)에서 확인할 수 있다.

M=1.0에서 전체 영역에 걸쳐 평균 막냉각 효율이 높게 예측되는 경향이 나타났으며 하류로 갈수록 막냉각 효율이 감소하는 경향은 실험과 수치해석 간에 일치하는 경향을 보였다. 이는 Fig. 8(c)와 8(d)에서 확인할 수 있다. 위 조건에서 주유동에 대한 상대적인 냉각 유체의 관통력이 증가하고, 냉각막이 벽면을 따라 길게 유지되며, 하류 방향으로의 냉각 효율이 개선된다.

M=2.0에서 평균 막냉각 효율을 전체적으로 실험보다 과대 예측되는 양상이 나타났다. 이는 Fig. 8(e)와 8(f)에서 확인할 수 있다. 분사비가 증가함에 따라 냉각 유체의 운동량이 커져 주유동을 강하게 관통하게 되며, 그 결과 중심축을 따라 길게 유지되는 고효율 영역은 유지되지만, 횡방향으로의 분포 폭은 오히려 감소하는 경향을 보인다. 일반적으로 분사비가 증가하면 막냉각 효율은 향상되는 경향을 보이지 않지만, 일정 수준 이상의 분사비에 도달하면 효율이 감소하는 현상이 나타난다. 이는 분사비 증가에 따라 냉각 유체의 운동량이 커지면서, 냉각 유체가 주유동을 관통해 벽면으로부터 멀어지기 때문이다. 따라서 최적의 냉각 효과를 얻기 위해서는 적절한 분사비를 고려해야 한다.

Fig. 10은 분사비에 따른 x/D = 25 단면에서의 냉각 효율 및 속도 벡터장을 통해 유동 구조의 변화를 나타낸다. 분사비가 증가할수록 Fig. 9에서 관찰된 유동의 가속과 불균일한 속도장이 두드러지며 후류 방향에서 CRVP(Counter Rotating Vortex Pairing)현상이 나타난다.

Fig. 10

Velocity vectors and film cooling effectiveness contour plots at x/D=25

M=0.5일 때 유동 내 와류 형성이 미약하여 냉각막이 장거리 유지되기 어려운 구조를 가지며, Fig. 10(a)에서 확인된다.

M=1.0일 때 유동 구조에서 냉각 유체 주입 후 유동 내에 형성되는 CRVP 구조가 냉각 유체와 주유동 간의 급격한 혼합을 억제하고, 냉각 유체가 벽면을 이탈하지 않게 하는 흐름이 유지된다. 와류의 구조는 Fig. 10(b)에서 확인할 수 있다.

M=2.0일 때 CRVP 구조를 Fig. 10(c)에서 가장 명확히 확인할 수 있다. CRVP는 벽면 상부에 형성되어 냉각 유체를 들어 올리는 유동을 유도하기 때문에, 냉각 유체를 벽면에서 멀어지게 하며, 중심축 방향으로 상승하는 흐름을 유도하여 냉각 유체의 벽면 접촉을 약화시킨다. 그 결과 하류 방향으로 갈수록 막냉각 효율의 횡방향 확산을 억제하는 효과를 나타낸다.

Fig. 11은 분사비에 따른 횡방향 평균 단열 막냉각 효율을 비교한 것으로, 실험 결과⁽⁶⁾와 CFX⁽¹⁵⁾, dbturb 결과를 비교하였다.

Fig. 11

Comparison of laterally averaged film cooling effectiveness between CFD results and experiment results(6)

Fig. 11(a)에서는 dbturb 해석 결과가 후류영역으로 갈수록 실험값과의 오차율이 점진적으로 증가하는 양상이 나타난다. 이는 후방 영역에서 냉각 유체의 확산 및 유지가 실제 실험보다 상대적으로 크게 모사하기 때문으로 해석할 수 있다.

Fig. 11(b)에서는 dbturb 해석 결과와 실험값의 결과가 전체 x/D 구간에 걸쳐 일관된 차이를 보인다. 이는 dbturb의 결과가 후방 영역까지 예측 오차가 크게 나타나지 않고 있으며, dbturb가 비교적 안정적으로 구현되고 있다.

Fig. 11(c)에서 dbturb는 전체적으로 실험 결과보다 높은 막냉각 효율 값을 예측하고 있으며, 후방 영역으로 갈수록 실험값과의 차이가 증가하는 양상을 보였다. 그러나 전반적으로 분사비가 증가할수록 CFX⁽¹⁵⁾와 dbturb 두 수치해석 결과는 실험값보다 높은 평균 막냉각 효율을 예측하는 경향을 보였다. 이는 일반적으로 RANS 기반 해석의 난류 모델이 와류 구조를 평균화하여 처리하는 영향으로 판단된다⁽¹⁶⁾. Fig. 12는 분사비에 따른 x/D=3에서 x/D=25 범위의 막냉각 효율 면적 평균값을 나타냈다. dbturb와 CFX⁽¹⁵⁾ 수치해석 값은 실험과 유사한 양상을 보이지만, dbturb의 평균 오차율은 약 25.8%, CFX⁽¹⁵⁾ 평균 오차율은 약 46.1%로 나타나 dbturb의 예측 정확도가 상대적으로 우수하게 나타났다.

Fig. 12

Comparison of area-averaged film cooling effectiveness between CFD and experiment results(6)

막냉각 효율의 실험적 측정은 다양한 요인으로 인해 상당한 어려움이 따른다. 특히, 면적 평균 효율값을 산출할 때 시간 변화나 주변 환경 조건에 따라 오차가 발생하는 것으로 알려져 있다. 예를 들어, Kim⁽¹⁷⁾ 등이 실험한 PSP 기법을 활용한 실험에서는 특정 효율값에서 ±5.2%에서 최대 ±16.1%의 불확실성이 보고된 바 있다. 따라서 막냉각 효율의 절대적인 정확도보다는 분사비에 따른 전체적인 경향성과 분포 패턴을 정밀하게 예측하는 것이 더욱 중요하다. 본 연구에서 사용된 dbturb 해석자는 분사비 변화에 따라 실험과 매우 유사한 경향성을 안정적으로 재현하였으며, 이는 해석자의 신뢰성을 입증하고 향후 다양한 연구 분야에 적용할 수 있는 가능성을 보여주는 결과라 할 수 있다.