1. 서 론

항공기가 비행 시 고려해야 할 사항 중 착빙현상이 있다. 착빙현상은 과포화된 액적이 항공기 표면에 부딪혀 얼어붙는 현상을 말한다. 착빙이 지속되면 표면의 얼음이 성장하므로 항공기의 형상이 변하게 된다. 특히 양력을 발생시키는 날개에서 착빙이 발생하면 양력이 감소하고 항력이 증가하는 등 항공기 공력성능이 저하된다. 극단적으로는 실속으로 인한 추락사고까지 발생할 수 있다는 문제점이 있다[1].

이러한 문제를 방지하기 위해 착빙현상에 관련한 연구가 지속적으로 수행되고 있다. 선행연구에서는 착빙풍동시험을 통해 익형의 착빙형상과 공력계수를 모두 측정하였다[2]. 또한 착빙형상을 직접 제작하여 풍동시험을 수행한 바 있다[3,4]. 그리고 익형이나 비행체 형상에 대한 전산유체역학(CFD) 해석을 수행하여 유동 특성을 분석하였다[5,6].

착빙풍동시험 수행 시 얼음으로 인한 날개 형상과 공력계수의 변화를 계측할 수 있다. 그러나 날개 양쪽 끝의 endplate와 날개 표면에서 발생한 경계층의 상호작용으로 인해 3차원 유동구조가 형성되므로 목표했던 2차원 익형의 착빙 형상과 소폭 달라지며[7,8], 공력계수와도 차이가 발생한다[9]. 이러한 Endplate 효과에 대해 분석하는 연구[10]는 이전부터 있었으나, 일반적인 풍동시험 보정 방법론[11,12]을 적용하기가 쉽지 않다. 또한 착빙 형상의 공력을 근사적으로 계산하는 방법론을 제안한 선행연구[13,14]가 있지만 3차원 날개 형상에 패널법을 이용하므로 본 연구와 같이 경계층의 영향이 큰 경우에는 방법론을 적용하기 어렵다. 따라서 착빙풍동시험 환경에서 벽면 경계층 상호작용을 반영하는 간단한 보정 방법론이 필요하다고 판단하였다. 해당 보정 방법론을 통해 endplate 효과가 제거된 2차원 착빙 익형의 공력계수를 얻는 것을 목표로 하였다.

본 연구에서는 착빙풍동시험을 통해 얻은 착빙형상의 공력성능을 보정하는 효율적인 방법론을 제안하였다. 먼저 CFD로 계산한 2차원 익형과 3차원 날개의 공력계수를 풍동시험 결과와 비교하여 CFD 해석 신뢰성을 확인하였다. 그리고 보정 방법론의 검증을 위해 3차원 착빙 날개의 항력계수를 보정하고 풍동시험 결과와 비교하였다. 3차원 착빙 이전 날개의 항력계수에 CFD로 구한 익형의 착빙 전·후 항력계수 변화량을 더하여 성능을 검증하였다. 최종적으로 보정 방법론을 적용하여 착빙으로 인한 날개의 3차원 항력계수 변화량을 2차원 익형 항력계수에 더해 보정하였으며, 보정 방법론으로 도출한 착빙 익형의 공력계수를 CFD 해석 결과와 비교하였다.

2. 수치해석기법

2.1 유동해석자

유동해석자는 인하대학교 공력해석 및 설계 연구실의 In-house 해석자인 UMSAPv[15]을 사용하였다. RANS(Reynolds Averaged Navier-Stokes) 방정식을 해석하였으며, 공간이산화는 유한체적법을 사용하였고 MUSCL 외삽법을 통해 2차 정확도를 구현하였다. 비점성 수치유속은 Roe의 근사리만 해법으로 계산하고 점성 수치유속은 중앙차분법으로 계산한다. 낮은 마하수 유동 해석을 위해 예조건화 기법을 적용하였다. 정상상태 시간적분은 LU-SGS, 비정상 시간 적분은 이중 시간 적분법을 사용하였다. 난류모델은 Menter의 k-ω SST 모델[16]을 사용하였다. 난류모델 방정식은 식 (1)과 (2)와 같다. 𝜌는 밀도, k, 𝜔는 2변수 난류모델의 난류변수이고 각각 난류운동에너지와 소산율이다. v→는 속도, τij는 shear stress, μm은 층류점성계수, μt는 난류점성계수, νt는 난류동점성계수이다. F1은 blending function이고 β*, 𝜎 등은 모델 계수로 참고문헌[16]에서 제시된 값을 사용하였다.

착빙 이전의 경우에는 표면을 매끈하다고 가정하였으나 착빙 이후의 경우에는 착빙에 의한 표면이 거칠어진다. 이를 고려하기 위해 착빙 후 공력 해석시 날개 표면에서는 Knopp이 제안한 표면 거칠기 모델[17]을 적용하였다. 이 경우 거친 표면에서의 난류 변수 kw와 ωw는 식 (3)과 (4)로 계산한다. 여기서 τw는 wall shear stress, uτ는 friction velocity, y는 벽면에서의 최소거리이다. 식 (5)부터 식 (7)은 식 (3)과 (4)를 계산하기 위한 parameter이다. 표면 거칠기 ks에 따라 kw와 ωw의 값이 정해진다.

3. 착빙풍동시험 모사

3.1 착빙풍동시험

착빙풍동시험은 Austria Vienna에 위치한 RTA의 Icing Wind Tunnel[18] 에서 수행하였다. 풍동 시험부의 가로, 세로, 높이는 각각 2.9 m, 3.0 m, 1.7 m이다. 풍동의 구조와 시험부 사진을 Fig. 1에 나타내었다. 풍동에서 시험한 날개 형상의 시위길이는 1 m, 스팬길이는 2 m이다. Fig. 1의 (b)에서 확인할 수 있듯이 날개의 양 끝단에 endplate가 존재한다.

Fig. 1.

Icing wind tunnel: (a) schematic view of icing wind tunnel, and (b) front view of test section

착빙시험 전 Dry air 조건(Clean 조건)에서 시험한 날개의 공력계수는 Fig. 2와 같다. 유동조건은 마하수 0.14, 레이놀즈수 3,500,000이고, 받음각은 –4도부터 10까지 2도 간격이며 공력 계수를 계측하였다. 보안사항으로 인해 구체적인 수치는 공개할 수 없어, 모든 그래프의 공력계수를 Clean 조건 받음각 10도의 값으로 정규화하여 나타내었다.

Fig. 2.

Icing wind tunnel airfoil aerodynamic coefficients: (a) α-CL, and (b) α-CD

Table 1에 비교 및 검증 목적으로 사용할 시험조건을 나타내었다. Rime ice가 1가지, glaze ice가 2가지로 나타난 조건 중에서 공력계수 변화량이 가장 큰 3가지 조건이다. LWC(Liquid Water Contents)는 단위부피당 액적의 양을, Spray time은 착빙시간을, MVD(Mean Volume Diameter)는 액적의 지름을 나타낸다. Fig. 3은 착빙 전·후의 양항곡선(Drag polar)를 비교하였다. 착빙풍동의 가동 범위로 인해 착빙 장치를 작동하지 않은 조건(Clean)과 장치 작동을 고려한 착빙이 발생하기 전인 조건(Before icing) 의 차이가 소폭 있지만 공력 계수는 거의 같다. 현재 관찰되는 미세한 차이는 착빙 장치 작동으로 인한 습도 차이 등 시험조건 차이가 원인인 것으로 판단된다. 착빙된 이후(After icing) 에는 형상 변화에 의해 양력계수는 감소하고 항력계수는 증가하였다. 시험조건별 Before icing 대비 After icing의 공력계수 변화량을 Table 2에 정리하였다.

Fig. 3.

Drag polar for icing wind tunnel iced airfoil

착빙풍동시험 후에 착빙된 날개 표면을 레이저 스캔하여 고해상도의 3차원 형상 데이터를 생성하였다[18]. 보안상의 이유로 전체적인 형상은 나타낼 수 없어 스캔한 형상의 일부분만을 Fig. 4에 나타내었다. Fig. 4의 중앙지점이 스팬의 중심이고 삼각형 격자로 얼음 형상이 세밀히 표현되었다. 스팬방향으로 고르지 않은 표면이 계속되는 등 일반적인 착빙형상의 특징이 관찰된다.

Fig. 4.

Iced airfoil configuration scan image

3.2. CFD 해석

착빙 전·후의 공력 변화를 모사하기 위한 CFD 해석을 수행하였다. 착빙 전 조건(Clean) 에서 2차원 익형 단면(Clean 2D)과 풍동 시험부를 고려한 3차원 날개(Clean 3D)를 해석하였으며, 착빙 시험 이후 생성된 얼음 형상을 고려한 2차원 익형 단면(Ice 2D) 주위에서 유동을 해석하였다.

우선 Clean 조건에서 2차원 익형을 해석하기 위해 익형을 가운데 배치하고, 시위길이의 40배의 정사각형 원방을 계산영역으로 설정하였다. 시위길이 방향으로 위, 아래 각각 400개의 격자점이 있으며 고받음각 조건에서 wake를 포착하기 위해 익형 후류에 시위길이의 0.01배에 해당하는 조밀한 격자를 구성하였다. 경계층 모사를 위해 벽면에서 수직 방향 높이인 y+가 1에 근접하도록 격자를 구성하였으며, 전체 격자수는 약 47,000개이다. 익형 표면에는 No-slip wall 경계 조건, 원방 경계에는 Riemann invariant를 이용한 Far-field 경계 조건을 부여하였다.

풍동 시험부에서 날개의 3차원 공력은 참고문헌[18]의 풍동 시험부 정보와 사진을 바탕으로 Fig. 5와 같이 계산 영역을 구성하였다. 특히 날개 양 끝단에 장착된 측면 벽의 영향을 고려하도록 경계조건을 구분하였다. 비정렬격자를 구성하였으며, 익형 주위에서 시위방향으로 위, 아래 각각 200개의 격자점을 위치시켰으며 스팬방향으로 140개의 격자점이 존재하도록 하여, 총 5,900,000 개의 격자를 구성하였다. 날개 표면과 측면벽의 수직한 방향으로 발달하는 경계층을 포착하기 위해 수직 방향 높이인 y+가 1에 근접하도록 격자를 구성하였다. 풍동 입/출구에는 Far-field 경계 조건을 부여하였다. 날개표면과 양 끝단에 있는 end plate가 포함된 측면 벽면을 회색으로 표시하고 No-slip wall 경계 조건을 부여하였다. 노란색으로 표시한 나머지 영역은 공동으로 Slip wall 경계 조건을 부여하였다.

Fig. 5.

Clean 3D grid: (a) entire domain, (b) airfoil upper surface grid

착빙 후 날개의 공력 변화를 모사하기 위해 착빙된 날개 단면의 형상 모델링을 진행하였다. 참고문헌[8]에 따르면 착빙 형상에 대한 endplate의 영향은 날개 양 끝의 경계층에 제한적이다. 따라서 Fig. 4와 같은 실제 착빙형상에 대해서 스팬방향으로 32등분된 지점의 양끝을 제외한 단면을 종합하여 평균한 단면과 가장 유사한 단면을 형상으로 하였다. 정확한 계산을 위해서는 각 단면의 공력을 계산 후에 평균하는 것이 옳으나 계산 효율성을 중시하여 단면을 평균화하였다. 계산영역과 경계조건은 Clean 2D와 동일하나 앞전 부근에 격자를 집중시켜 얼음으로 인한 형상 변화를 세밀히 모사하였다. 시위방향으로 위, 아래 각각 560개의 격자점이 배치하였으며 전체 격자수는 약 84,000개인 비정렬 격자를 구성하였다. Fig. 6이 얼음 영역 격자를 확대한 그림이다. Clean 2D와 달리 Ice 2D는 No-slip wall 경계조건과 함께 표면 거칠기 모델까지 추가로 적용하였고 표면 거칠기 ks는 Shin의 방법[19]으로 계산한 값인 1 mm이다.

Fig. 6.

Ice 2D grid

착빙 전 조건에서 2차원 및 3차원 해석인 Clean 2D와 Clean 3D는 정상 RANS 유동 해석을 수행하였다. 착빙 후 2차원 익형 유동 해석인 Ice 2D의 경우 앞전에서 발달한 얼음에 의한 비정상 유동 구조를 포착하기 위해 비정상 RANS 유동 해석을 수행하였다. 시위 길이(c)와 자유류 음속(a∞) 으로 무차원화한 시간(t∞=c/a∞)을 기준으로 0.01 무차원시간 간격으로 200 무차원시간 동안 계산을 수행하였으며 마지막 20 무차원시간 동안 평균 공력을 산출하였다.

4. 결 과

4.1 공력 해석 결과 분석

착빙 전 익형에 대해 풍동시험과 CFD 결과를 비교하였다. 2차원 CFD 해석의 경우 측면 벽에서 발달하는 경계층과 풍동시험부와의 간섭을 모사하지 못해서 풍동시험과 공력 특성이 다를 수 있다.

해석결과, 양력이 0인 받음각(αL=0)과 항력이 최소인 받음각이 풍동시험과 CFD 해석 결과가 2도만큼 차이가 발생하였다. 시험을 수행한 착빙풍동 소개자료[18]에서 언급되었듯이 착빙모사를 위해 시험부의 유동조건에 제약이 발생하는 점을 고려할 때, 시험부에 Flow angularity가 존재하며 받음각을 약 2도 만큼 보정하는 작업이 필요하다고 판단하였다. 2도만큼의 받음각 보정을 수행한 풍동시험 결과(Corrected)와 2차원 익형(Clean 2D) 및 3차원 날개(Clean 3D) 에 대한 CFD 결과를 비교하면 Fig. 7과 같다. 풍동 시험부를 모사한 Clean 3D 해석의 공력 계수가 풍동시험과 매우 잘 부합한다. 양력계수는 거의 일치하며, 항력계수는 약간 낮게 예측하였다. Clean 2D, Clean 3D 및 풍동 시험을 비교한 결과 풍동의 벽면을 모사한 Clean 3D가 풍동시험과 유사하였다. Fig. 7을 통해 풍동의 endplate에서 발달한 경계층과 날개 표면 경계층의 상호작용으로 인해 받음각에 따른 양력 기울기가 감소한 것을 확인할 수 있다.

Fig. 7.

Clean airfoil aerodynamic coefficient after correction: (a) α-CL, (b) α-CD, (c) CD-CL

CFD 해석 결과에서도 경계층 상호작용으로 인한 공력계수 변화를 확인할 수 있다. Fig. 8은 받음각 4도에서 익형 윗면의 압력계수분포를 나타낸 그림이다. (a)는 Clean 3D의 half span이며 (b)는 비교대상인 Clean 2D이다. Clean 3D의 압력분포가 스팬의 중심에서는 endplate와 멀리 떨어져 있어 Clean 2D의 압력분포와 유사한 반면 양 끝은 경계층으로 인해 상당한 차이를 보인다.

Fig. 8.

Comparison of airfoil upper surface Cp contour: (a) Clean 3D, (b) Clean 2D

착빙풍동시험 결과와 Clean 2D 및 각 조건별 Ice 2D CFD 해석 결과를 비교하면 Fig. 9와 같다. 풍동시험과 마찬가지로 2차원 CFD 해석에서도 착빙된 얼음에 의해서 양력계수가 감소하고 항력계수가 증가하는 것을 확인할 수 있다. 착빙조건은 다르지만 양항곡선은 거의 같은데, 3가지 조건에서 앞전의 얼음뿔 형태가 유사한 것이 원인으로 여겨진다. 풍동시험 결과와는 차이가 발생했으며 가장 주요한 원인은 2차원 해석이므로 endplate에 의한 공력계수 변화가 반영되지 않았기 때문이다. 선행연구[20]와 같이 고정밀 해석기법을 사용하여 난류 천이 및 박리 유동을 정밀하게 해석하는 방법도 고려할 수 있으나 계산비용이 매우 많이 소모되는 문제가 있다. 본 연구에서는 벽면 효과를 효과적으로 포착하기 위해 비정상 RANS 해석을 수행하였다.

Fig. 9.

Drag polar for experiment and numerical calculation

4.2 보정 방법론

본 연구에서는 착빙 전·후 공력 차이를 효율적으로 추정하기 위한 공학적 근사 모델로서, 착빙 형상의 공력계수를 착빙 전 공력계수에 착빙 전·후 공력계수 차이로 보정하는 방법론을 개발하였다. 일반적으로 물체의 형상과 착빙의 영향을 독립적으로 계산할 수 없지만, 간단한 방법론 적용을 위해 그 상관관계가 작다는 가정을 하였다. 즉, Clean 2D 공력 해석 결과에 착빙으로 인한 변화에 해당하는 Clean 3D와 Ice 3D의 공력 차이를 더해 보정한다. 이 방법은 복잡한 공력 모델링을 요구하지 않고 양항곡선만으로도 계산이 가능하다는 장점이 있다.

개발한 방법을 수식으로 표현하면 식 (8)과 같다. 여기서 CdICE 는 착빙 익형의 항력이고, CdCLEAN는 착빙 전 익형의 항력이다. 항력계수는 양력계수의 함수로 표현할 수 있다. 착빙 효과를 모사하기 위해 3차원 착빙 유무의 차이인 ΔCD 를 보정하며, 식 (9)와 같이 계산한다. 식 (8)과 (9)의 2차원 익형의 양력 계수 (Cl)와 3차원 날개의 양력계수(CL)가 같다. 정리하면, 2차원 양력 계수 Cl을 기준으로 해서 Clean 2D의 항력에, 같은 값을 가지는 3차원 날개의 앙력계수 CL을 기준으로 착빙 유무에 따른 항력 변화량(Ice 3D와 Clean 3D의 차이)를 더하여 보정한다.

이 과정을 Fig. 10과 같이 도식화할 수 있다. Fig. 10의 (a)는 식 (8)과 (9)를 나타낸 것이며 (b)는 방법론 검증을 위해 반대의 과정을 고안한 것이다. 여기서 반대의 과정은 착빙 전·후 익형 항력계수 변화량(Ice 2D와 Clean 2D의 차이)을 착빙 전 날개의 항력계수(Clean 3D)에 더하여 3차원 착빙 날개의 항력계수를 도출하는 과정이며 식 (10) 및 (11)과 같다.

먼저 검증을 위해 보정 방법론을 적용하여 풍동에서 착빙된 3차원 익형의 공력을 도출하였다. 식 (10)의 CDICE(CL) 계산 시 Clean 3D CFD 해석 결과를 CDCLEAN(CL)로 사용한 경우를 Corrected(CFD)로 하였다. 또한 Fig. 7에서 알 수 있듯이 Dry air 조건에서 풍동시험 결과와 Clean 3D가 거의 같다. 따라서 풍동시험 결과값을 CDCLEAN(CL)로 사용하여 착빙 효과를 도출할 수 있으며, 이 경우를 Corrected(Exp)로 하였다. 식 (10), (11) 계산 과정에서 양항 곡선은 스프라인 보간법으로 내삽하여 계산하였다. 보정한 공력을 풍동시험 결과와 비교하면 Fig. 11의 (a)와 같다. CDCLEAN(CL)의 경우 풍동시험과 CFD의 차이가 매우 작았기 때문에 Corrected(Exp)와 Corrected(CFD)의 차이는 거의 없었다. 낮은 양력계수 구간에서의 일부 오차를 제외하면 풍동시험 결과와 전반적으로 잘 일치한다. 그리고 같은 방식으로 식 (8), (9)를 적용하여 보정한 착빙 익형 공력 계수를 CFD 해석 결과(Clean 2D, Ice 2D)와 비교하여 Fig. 11의 (b)에 나타내었다. 검증 케이스와 마찬가지로 낮은 양력계수 구간에서 오차가 발생하였지만 대부분의 구간에서 Ice 2D와 유사한 경향성을 가진다.

Fig. 10.

Schematics of proposed method: (a) 2D iced airfoil, (b) 3D iced wing

Fig. 11.

Corrected drag polar: (a) 3D aerodynamics including endplate effect, (b) 2D aerodynamics of the airfoil

5. 결 론

본 연구에서는 착빙풍동시험에서 얻은 착빙 형상의 공력계수에 보정 방법론을 적용하여 2차원 착빙 익형의 공력계수를 도출하는 방법론을 제시하였다. 그리고 방법론을 적용하여 도출한 결과와 착빙풍동시험 결과를 비교 및 검증하였다.

보정 방법론을 요약하면, 착빙 형상의 공력계수를 착빙 이전의 공력계수와 착빙으로 인한 변화량으로 분리하고 각각을 계산하여 보정된 양항곡선을 도출하는 방법이다. 검증 목적으로 보정 방법론을 적용해 3차원 착빙 날개의 양항곡선을 도출하고 풍동시험과 비교하였으며 이후에 같은 방식으로 2차원 착빙 익형의 양항곡선을 계산하였다. 검증과정에서 착빙 전·후 익형의 공력계수가 필요하므로 CFD 해석을 수행하였다. 또한, 검증에 필요한 착빙 전 날개의 공력계수에 대해서는 풍동시험 데이터와 CFD 해석 결과를 비교하여 신뢰성을 확보하였다. 그 결과, 보정된 3차원 착빙 날개와 풍동시험의 양항곡선이 대부분의 양력계수 구간에서 일치하였다. 그리고 착빙 전·후 날개의 공력계수와 착빙 전 익형의 공력계수를 결합하여 착빙 익형의 양항곡선을 얻었다. 2차원 착빙 익형의 양항곡선도 착빙 익형의 CFD 해석 결과와 매우 유사하였다. 본 연구에서 제안한 보정 방법론은 양항곡선만을 사용하므로 일반적인 유동조건을 가진 착빙풍동시험 데이터에 적용하였을 때에도 효율적으로 공력성능을 파악할 수 있을 것으로 예상한다.

다만, 본 방법론은 받음각에 따른 결빙 조건이 유사하고, 착빙 형상이 스팬 방향으로 큰 비대칭 없이 비교적 균일하게 형성되는 경우에 한하여 적용 가능하다. 만약 복잡한 3차원 착빙 형상이 발생하거나 형상 변화가 급격한 경우에는, 2차원 단면 기반의 보정 기법만으로 전체 유동 특성을 충분히 대표하기 어려우며, 이에 따른 적용 한계가 있다. 향후에는 보다 다양한 착빙 조건에 본 보정 방법론을 적용하고, 그 결과를 풍동시험과 비교함으로써 제안한 기법의 확장 가능성을 평가하고 신뢰성을 확보하고자 한다.