1. 서 론
항공기 결빙은 항공기가 과냉각 액적 구름을 통과할 때 표면에 얼음이 형성되는 현상이다. 실제로 항공기 표면에 결빙이 발생하면 이는 공기역학적 성능 감소를 일으킨다. 따라서 일반적인 항공기는 결빙보호장치를 사용하여 표면을 보호하거나, 결빙보호장치가 없는 항공기의 경우, 결빙 조건을 피해 운항해야 한다[1]. 그러나 전기 추진 항공기의 경우 결빙보호장치를 위한 에너지를 추가적으로 고려해야 하고, 이는 곧 배터리의 중량 증가를 야기하기 때문에 결빙에 특히 취약하다는 단점이 있다[2].
특히 항공기에 발생하는 결빙은 결빙 조건에 따라 Rime ice 및 glaze ice로 구분된다[1]. Rime ice는 낮은 온도 조건에서 항공기 표면에 액적이 부딪히는 순간 바로 얼게 되는 얼음으로 불투명한 얼음이 특징이다. 또한 이 얼음은 들어오는 액적의 분포에 따라 얼음 형상이 나타나며, 이 얼음을 예측하기 위해서는 들어오는 액적의 양을 정확히 예측하는 것이 중요하다. Glaze ice는 상대적으로 높은 온도 조건에서 항공기 표면에 액적이 부딪히면 바로 얼지 않고, 액체 상태로 표면을 따라 흐르면서 어는 얼음으로 투명한 얼음이 특징이다. 이 얼음은 액적의 분포보다 대류열전달로 인해 발생하는 응결 현상이 형상에 더 크게 영향을 미친다. 따라서 에어포일 특성에 따라 대류열전달이 급격히 증가하는 구간에서 더 많이 얼게 되어 뿔 형태로 자라나는 것이 특징이다.
전기 추진 항공기는 추력 장치로 프로펠러 및 로터를 주로 사용하는데, 특히 프로펠러는 회전으로 인한 빠른 속도 및 짧은 코드 길이 특성으로 결빙에 더 취약하다는 단점을 갖는다. 프로펠러 및 로터와 같이 회전하는 물체에 발생하는 결빙의 경우, 온도에 관계없이 root 영역은 속도가 낮아 들어오는 물의 양이 적기 때문에 대부분 rime ice 특성을 갖지만, tip 영역에서는 들어오는 물의 양이 많고 공력 가열이 발생하는 경우에 특히 glaze ice가 나타날 수 있다. 따라서 영역에 따라 발생하는 결빙 형상이 다르며 이를 정확히 예측하는 것은 어렵다. 또한 프로펠러에 발생하는 결빙은 공기역학적 성능 감소로 인한 효율 저하를 야기한다. 앞전에 발생하는 결빙 형상에 의해 회전에 필요한 토크가 증가할 수 있다. 게다가 표면에 발생한 얼음은 뿔 형태를 띄고 있어 같은 rpm 상에서 추력이 감소하는 경향을 보이며 이는 stall로 이어질 수 있다.
이번 연구에서는 프로펠러에 발생하는 결빙을 예측하고, 발생한 결빙이 주어진 결빙 조건에서 얼마나 효율 감소를 야기하는지 예측한다. 기존에는 회전하는 물체에 대한 결빙 연구는 주로 로터 블레이드에 발생하는 결빙 연구들이 진행되었다[3,4,5]. 회전하는 로터에 대한 2D 결빙 예측으로는 BEMT(Blade Element Momentum Theory)을 활용한 방법이 있다. 이 방법은 3차원 해석에 비해 정확하지 않으나 계산속도가 빠르다는 장점이 있다. 따라서 정확한 결과를 위해서는 3D 결빙 예측 방법으로 MRF(Multiple Reference Frame) 방법을 활용한다.
따라서 이번 연구에서는 결빙 조건에서의 프로펠러의 효율을 분석하기에 앞서 프로펠러에 대해 2D 결빙 예측 방법과 3D 결빙 예측 방법에 대한 결과를 비교한다. 이를 위해 기존 결빙 예측 코드인 ICEPAC[6]을 활용하였고, 회전하는 물체에 대한 해석 결과의 신뢰성을 검증하고자 실험 데이터와 결빙 형상을 비교하였다. 또한 예측된 결빙에 의한 프로펠러 형상 변화에 따른 효율 감소를 확인하는 연구를 진행한다. 주어진 결빙 조건에서의 얻은 결빙 형상 데이터를 바탕으로 이를 다시 수치해석을 통해 성능 변화를 확인한다. 또한 전체 결빙 조건에서의 프로펠러 성능 변화를 알아보기 위해 Kriging 기법을 활용하고, 이를 통해 결빙 파라미터들의 영향성을 판단한다. 결과적으로 전체 결빙 조건 중에서 프로펠러 성능에 크게 영향을 미치는 조건들을 분석한다.
2. 본 론
2.1 수치해석 기법
본 연구에서는 프로펠러에 발생하는 결빙 형상을 예측하고 성능 저하에 미치는 영향을 분석하기 위해 오픈소스 기반 CFD(Computational Fluid Dynamics) 해석 툴인 OpenFOAM[7]을 사용한 ICEPAC[6]를 이용하였다. 이 OpenFOAM은 N-S 방정식(Navier-Stokes’ equation)을 이용한 CFD 해석 뿐만 아니라 FVM(Finite Volume Method) 기반의 다른 물리적 해석에 대한 추가 계산이 용이하므로, ICEPAC과 같은 유동장 해석을 기반으로 하는 솔버 개발이 가능하다[7].
본 연구에 사용한 회전하는 프로펠러 결빙 해석 기법은 크게 2D BEMT 방법과 3D MRF 방법 두 가지로 나뉜다. 각각의 방법에서 ICEPAC을 기본적으로 활용하나, 회전하는 물체에 접근하는 방식이 다르다. 따라서 ICEPAC 및 프로펠러 결빙에 대한 각각의 해석 방법은 다음과 같다.
2.1.1 ICEPAC 결빙 예측 코드
먼저 ICEPAC은 네 가지 모듈로 구성되어 있다: 1) 공기역학, 2) 액적장, 3) 열역학, 4) 얼음 성장 모듈. 각 모듈은 준안정 가정을 적용하며, 각 모듈의 결과는 다음 모듈의 입력으로 순차적으로 사용된다.
1) 공기역학 모듈
공기역학 모듈에서는 압축성 레이놀즈 평균 나비에-스토크스 방정식을 OpenFOAM의 rhoPimpleFoam 솔버를 사용하여 해결한다. 이번 연구에서는 표면 위 발생하는 결빙 및 수막에 의한 완전히 발달된 난류를 가정하였고, 이를 해석하기 위해 SA 모델을 이용하였다. 또한 결빙이 발생하기 때문에 표면 거칠기 모델을 적용하였다.
2) 액적장 모듈
액적장 모듈은 오일러 접근법을 사용하며 액적에 작용하는 항력, 중력, 부력을 고려하여 액적의 거동을 모사한다. 액적에 대한 연속 방정식과 운동량 보존식은 다음 식 (1), (2)와 같다[6]. 이를 통해 얻은 액적의 속도 및 밀도장을 바탕으로 다음 식 (3)과 같이 물체 표면으로 충돌하는 액적의 양을 계산한다[6]. , 는 액적의 밀도와 속도를 나타내며, 는 항력계수, 는 액적의 레이놀즈수를 의미한다.
3)열역학 모듈
열역학 모듈에서는 SWIM[6]을 사용하여 수막의 질량과 에너지 보존을 해결한다. 액적장 모듈을 통해 얻은 표면에 충돌하는 액적에 의해() 생성되는 수막의 열전달 및 상태변화로 증발된 수막의 양()표면에 남아있는 수막의 양()과 축적된 얼음의 양() 및 표면 온도()가 결정된다.
4) 얼음 성장 모듈
마지막으로 얼음 성장 모듈에서는 축적된 얼음의 질량을 계산하고 그 결과 얼음 두께를 결정한다. 이 때 곡면을 따라 자라는 결빙 형상의 경우, 사다리꼴 형태를 갖기 때문에 높이에 대한 보정이 요구된다. 따라서 곡면에서 각 셀에서의 얼음 두께에 보정된 얼음 두께 값을 계산하여 반영한다.
2.1.2 2D BEMT 방법
프로펠러 결빙을 효율적으로 분석하기 위해서는 3차원 시뮬레이션이 아닌 2차원 접근 방식을 활용할 수 있다. 이 방법에서는 Fig. 1과 같이 프로펠러 블레이드를 스팬 방향을 따라 여러 개의 2차원 단면으로 나눈다. 이후 BEMT과 결빙 예측 코드를 결합하여 각 단면에서의 결빙 형상을 예측하고 이를 interpolation하며 최종 결빙 형상을 얻는다.
결빙 예측 코드의 입력 조건은 BEMT를 사용하여 계산된다. Fig. 2는 블레이드 요소의 속도를, Fig. 3은 유입 속도 성분 계산을 위한 그림을 나타낸다.
먼저 식 (7)을 통해 Prandtl의 팁 손실 계수를 사용하여 유입 비율을 계산한다.
다음으로, 식 (8), (9), (10)를 이용하여 유입각과 유입 비율을 상관시켜 유효 유입각과 유도 유입 속도를 결정한다. 결정된 입력 조건은 각 구간별로 결빙 형상을 예측하는 데에 활용된다.
2.1.3 3D MRF 방법
Moving Reference Frame(MRF)는 블레이드의 회전 효과를 좌표계를 변환함으로써 블레이드와 주변의 영역을 정지된 격자계로 해석하여 빠른 시간 내에 정상상태의 해석 결과를 획득할 수 있는 유동장 해석 기법이다. 이는 좌표 변환 시 발생하는 Coriolis와 Centrifugal force가 모멘텀 보존식에 고려됨으로써 3차원 및 회전 효과를 고려할 수 있다.
회전하는 블레이드의 결빙 해석을 위하여 MRF를 활용할 경우, 유동장과 액적장 모두 Coriolis와 Centrifugal forces를 적용해야한다. 또한 열역학 모듈에서는 회전 효과가 고려된 유동장의 wall shear stress 를 적용하여 3차원 프로펠러와 회전 효과를 반영할 수 있기에 BEMT보다 정밀한 열역학 모듈로 구성된다. 따라서 MRF를 이용한 결빙 해석에서는 기존 ICEPAC에 MRF 방법을 추가하여 결빙 예측을 했던 선행 연구 방법을 적용하였다[5]. MRF가 적용된 유동장 해석을 기반으로, 바뀐 액적장 모듈에서의 지배방정식은 원심력 및 코리올리 힘이 이 적용된 다음 식 (11), (12), (13) 과 같다.
2.2 수치해석 검증
회전하는 프로펠러 결빙 해석 기법 검증은 실험 결과와의 비교를 통해 수행되었다. 이는 기존 회전하는 물체에 대한 결빙 실험을 기반으로 수행되어야 하므로 NASA가 수행한 헬리콥터의 로터 블레이드에 대한 결빙 비행 시험(HIFT) 결과[9]를 이용하였다. 일반적으로 프로펠러는 로터보다 낮은 마하수 영역에서 작동한다. 따라서 낮은 마하수 영역에서 로터 블레이드의 root 부근의 결빙 형상에 대해 비교 검증을 수행하였다. 실험에 사용된 로터 블레이드는 코드 길이 0.5334m, 반경 7.3152m의 NACA 0012 에어포일 단면을 가진다. 이에 대한 외기 및 결빙 조건은 Table 1과 같다.
Table 1.
Icing condition of hovering rotor blade[9]
| Property | Value |
| Temperature[℃] | -30 |
| MVD[㎛] | 15 |
| LWC[g/㎥] | 0.2 |
| Radius[m] | 7.3152 |
| RPM | 324 |
| Time[s] | 180 |
먼저 3D MRF 방법에 대한 유효성은 Son에 의해 검증되었다[5]. 따라서 이번 연구에서는 2D BEMT 방법에 대한 유효성 검증을 나타내었다. 2D BEMT 방법에 대한 해석 결과는 다음과 Fig. 4와 같다.
해석된 결빙 형상은 FENSAP-ICE[10]로 계산된 얼음 모양과 실험 데이터[9]를 같이 비교하였다. Fig. 4에서 볼 수 있듯이 해석 결과가 실험 결과와 잘 일치하는 것을 확인했다. 다만 r/R=0.292에서 결빙 형상이 실험값에 비해 두께가 다소 작게 도출되었다. 이는 반경이 작은 영역에서는 rime ice가 나타나는데, 이 rime ice의 경우 밀도가 작다는 특징이 있다. 따라서 이를 반영하면 얼음 두께가 늘어나 실험 결과에 더 가까워질 것으로 보인다. 이번 연구에서는 이러한 얼음 밀도에 대한 데이터가 부족하여 얼음 형상에 따른 밀도 변화는 고려하지 않았다.
3. 해석 결과
결빙 시뮬레이션 환경을 설정하기 위해 한국항공우주연구원(KARI)에서 개발한 EAV-3[11,12]로 알려진 HALE UAV의 프로펠러를 선택했다. 이 기체의 제원은 Table 2에 요약되어 있다.
Table 2.
| Property | Value |
| MTOW[kg] | 66 |
| Rate of climb[m/s] | 0.75 |
| Cruise speed(Vcr)[m/s] | 7.8 |
| RPM | 925 |
| 𝑊battery[kJ] | 37,440 |
| Time[s] | 180 |
결빙 조건은 미국 연방항공규정(FAR) Appendix C[13]의 실제 항공기가 만날 수 있는 결빙 조건을 참고하였다. 이 규정은 Liquid Water Content(LWC), Median Volume Diameter(MVD) 및 온도와 같이 결빙에 영향을 미치는 파라미터 간의 상관관계를 간략하게 설명하고 있다. 항공기가 과냉각 구름을 만난다고 가정하였으며, Fig. 5와 Table 3에 표시된 대로 전체 결빙 조건 중 9개 지점의 LWC, MVD 및 온도 값을 이번 시뮬레이션 해석에 사용했다.
Table 3.
9-cases in icing conditions
| Case | LWC(g/m3) | T(℃) | MVD(μm) |
| 1 | 0.735 | -5 | 15 |
| 2 | 0.456 | -14.25 | 15 |
| 3 | 0.2 | -30 | 15 |
| 4 | 0.352 | -5 | 27.5 |
| 5 | 0.21 | -14.25 | 27.5 |
| 6 | 0.083 | -30 | 27.5 |
| 7 | 0.122 | -5 | 40 |
| 8 | 0.08 | -14.25 | 40 |
| 9 | 0.04 | -30 | 40 |
다음으로 고려해야 할 사항은 착빙 노출 시간이다. 이 역시 FAR Appendix C[13]와 임무 형상을 준수하여 결정된다. 결빙은 대기 온도가 0°C 이하일 때 발생하므로 결빙 조건에 대한 최소 고도는 표준 대기 온도가 0°C에 도달하는 고도인 2.3 km로 설정되었다. 또한 결빙 조건의 최대 고도는 HALE UAV가 마주칠 것으로 예상되는 성층 구름의 최고 높이를 고려하여 6.7 km로 설정했다. 항공기가 0.75 m/s의 일정한 상승 속도를 유지하면서 2.3 km에서 6.7 km 사이의 결빙 조건에 노출된다는 점을 감안하여, Fig. 6와 같이 총 결빙 노출 시간은 약 1.63시간으로 결정되었다.
위와 같이 주어진 결빙 조건에서 프로펠러에 대한 결빙 예측 해석을 수행하였다. 결빙 형상은 2D BEMT 방법과 3D MRF 방법을 비교하여 나타낸다. 먼저 3D MRF의 경우 Fig. 7과 같이 3차원 결빙 형상이 나타나므로 분석시 각 섹션별 결빙 형상을 비교한다. Fig. 7의 Glaze 및 Rime 결빙 조건은 각각 Case 1과 Case 3에 대한 형상이다. 두 조건은 같은 MVD를 가지고 있기 때문에 같은 collection efficiency를 갖는다. 그러나 Case 1의 Glaze 결빙 조건에서는 온도가 높고 LWC가 크기 때문에 뿔의 형상 및 결빙 형상의 크기가 Case 3의 Rime 조건에 비해 크게 나타난다.
Case 1,3의 각각의 결빙 조건에 대한 프로펠러 섹션의 결빙 형상에 대해, 각 결빙 예측 방법으로 얻은 결빙 형상을 비교한 그림은 Fig. 8, 9와 같다. 특히 프로펠러는 동체나 날개에 비해 코드 길이가 짧고 속도가 상당히 빠른 것이 특징이며, 이로 인해 결빙 형상이 비교적 크게 형성되는 것을 알 수 있다. 특히 LWC가 큰 Case 1의 경우, 앞전에 결빙이 집중되고 얼음이 많이 쌓이는 것을 볼 수 있다.
또한 각각의 조건에 대한 결빙 형상에 대해 2D BEMT 방법과 3D MRF 방법의 해석 결과가 비교적 일치함을 보였다. 다만 Fig. 9의 Case 3와 다르게 Fig. 8의 Case 1에서는 반경이 큰 위치(r/R = 0.9)에서 뿔의 위치에 차이를 볼 수 있는데, 이는 얼지 않은 수막이 원심력에 의해 회전 반지름 방향으로 넘어가 얼어붙는다는 것을 알 수 있다. 그러나 코드 길이가 상대적으로 작아지는 위치인 만큼 전체 결빙 형상에 큰 영향을 주지 않는다. 결국 이번 해석에 진행된 프로펠러의 경우, 상대적으로 작은 스팬 길이를 갖고 있기 때문에 회전하는 물체에 대한 3D 효과가 크지 않음을 알 수 있다.
다음으로 위와 같이 얻은 결빙 형상을 이용해 프로펠러 성능을 비교 분석하였다. 먼저 프로펠러 성능은 다음 Table 4와 같다.
Table 4.
Performance analysis of propeller
| Property | Value |
| Advance ratio() | 0.42 |
| Thrust coefficient() | 0.065 |
| Power coefficient() | 0.036 |
| Efficiency(𝜂) | 0.747 |
결빙이 발생한 프로펠러의 성능은 앞에서 얻은 결빙 형상을 통해 격자를 재생성하여 공력 해석을 수행하여 얻었다. 결빙 예측을 진행한 OpenFOAM을 통해 유동장 해석을 수행하였고, 총 9개의 조건에 대한 공력 해석 수행 결과는 Fig. 10과 같다.
전체적으로 프로펠러에 결빙이 발생하게 되면 추력이 감소하고 토크가 증가함을 알 수 있다. 추력은 Case 1을 제외하고 크게 감소하지 않는 반면, 토크는 많은 변화를 보였다. 특히 Case 1,2,4의 경우에 가장 많은 효율 변화를 보이며 이는 LWC가 모두 높은 조건에 해당한다.
Kriging 기법을 이용하여 전체 결빙 조건에 대한 프로펠러 효율 변화는 다음 Fig. 11과 같다. 이를 통해 각 결빙 파라미터의 영향을 알 수 있는데, 먼저 LWC가 높은 영역에서 낮은 프로펠러 효율을 확인할 수 있고, MVD의 영향이 크지 않음을 알 수 있다. 결국 LWC 0.2 g/㎥ 이하인 조건 영역에서만 기존 프로펠러 효율의 25% 이상의 효율을 가짐을 확인하였고, 대부분의 영역에서 큰 효율 감소를 확인할 수 있었다.
4. 결 론
본 연구에서는 프로펠러 결빙 영향성을 판단하기 위해 수치해석 방법을 적용하였다. 2D 및 3D 해석 방법을 통해 프로펠러 결빙 예측 해석 결과를 비교하였다. 2D 및 3D 결빙 해석 결과에는 tip 쪽에서 뿔의 위치에 대한 미세한 차이가 있었으나, 크기가 작기 때문에 결빙 두께 및 양에는 큰 차이를 보이지 않았다. 예측된 결빙 형상으로 격자를 재생성한 후 수치해석을 통해 프로펠러 성능 저하를 분석하였다.
결빙 조건에서 비행하는 전기 추진 항공기는 프로펠러 효율에 대한 결빙 영향성 분석이 중요하다. 이번 연구를 통해 결빙으로 인한 프로펠러 성능 저하 원인을 확인하였다. 결빙 형상을 결정하는 결빙 파라미터 중에서도 특히 LWC가 높은 조건일 경우, 프로펠러의 효율이 크게 낮아짐을 확인하였다. 빠른 회전 속도를 갖는 프로펠러의 특성으로 인해 비교적 온도가 높았음에도 결빙이 잘 형성되고, 작은 스팬 길이를 갖기 때문에 MVD 영향이 미비하기 때문이다. 이번 연구를 통해 프로펠러가 결빙 조건에 취약한 특성을 갖는다는 것을 수치해석을 통해 확인하였다.
