I. 서  론

공동(cavitation)의 발생은 선박의 운용에 있어 프로펠러의 추진효율을 저하시키고 선체 기진력을 증가시키는 등의 문제를 야기한다. 특히 고속함정의 프로펠러와 타(rudder) 주위에 발생하는 공동은 수중방사소음의 주요 원인으로 작용하기 때문에 공동의 발생 특성을 정량적으로 평가할 수 있어야 한다. 본 연구에서는 공동의 정량적 평가를 위한 해석법의 개발을 목표로, 먼저 속도 포텐셜을 기저로 하는 판요소법^[1]을 확장하여 부분공동 유동을 해석할 수 있는 수치해석 프로그램을 개발하였다. 현재 날개단면에서 발생하는 부분 공동 문제에 대해서는 정량적으로 비교가 가능한 실험 결과가 부족한 실정이며, 비점성 해석법을 사용한 결과는 공동이 발생한 후에 작용하는 항력 및 비정상 상태에서의 계산에 한계가 있다. 최근 여러 점성해석 상용 프로그램들이 다양한 난류와 공동 모델을 조합하여 공동해석을 수행할 수 있는 단계로 발전하고 있으며,^[2] 본 연구에서는 그 중 FLUENT에서 제공하는 난류와 공동 모델을 조합하여 성능해석을 수행하고 그 결과를 비교, 분석하였다. 이를 통해 부분 공동문제에 대한 난류 및 공동 모델의 조합을 사용함에 따른 점성해석 성능을 평가하고 점성효과를 고려한 비점성 해석법의 성능향상을 위한 방안을 모색하고자 하였다.

II. 수치해석 방법

2.1 비점성 수치해석 방법

입사각을 갖고 균일 유동 중에 일정한 속도로 움직이는 2차원 날개 단면을 가정한다. 날개 단면 표면에 법선 다이폴과 소스를 분포함으로써 경계치 문제를 정식화하였다. 발생한 공동의 표면에 대해 시간에 따른 변화량은 없고, 공동의 표면에서의 압력은 증기압으로 동일하다고 가정한다. 즉, 정상 상태의 무한 유동 가정 하에서 공동 표면에서의 접선 방향 속도는 모두 일정하다. 이러한 공동 표면에서의 동역학적 경계 조건을 먼저 적용하고, 선형 공동 마감 모형의 사용에서 나타나는 오류항을 수정하는 방식으로 공동 표면을 계산하며, 이때의 오류항은 공동 표면에 분포된 소스들의 세기의 합이 0이 되어야 한다는 수치 해석적 조건식으로 표현된다. 발생한 공동의 두께는 얇은 비행익 이론^[3]에 따라 공동 표면에서의 소스의 세기로 나타나며, 그 끝에서의 두께는 0 이라는 마감조건을 바탕으로 반복 계산을 수행하여 그에 따른 날개 두께의 수정 량을 구한다.^[4] 경계면으로 둘러싸인 무한유체 영역에서 라플라스 방정식을 지배방정식으로 하는 유기 속도 포텐셜 는 그린 정리를 만족하며, 이를 바탕으로 법선 다이폴과 소스의 세기를 미지수로 갖는 Eq.(1)을 구한다. 여기서 날개 단면의 경계를 , 발생한 공동의 경계를 라고 하며 G 는 그린 함수를 의미한다.

. (1)

아래 Eq.(2)는 이산화한 적분방정식이며, 이 때 는 번째 판에서의 전체 속도 포텐셜, 는 다이폴의 개수, 는 소스의 개수, 는 전체 미지수의 개수 (), 와 는 각각 , 번째 판에서의 다이폴과 소스의 세기를 나타낸다.

. (2)

또한 와 는 각각 번째 판에 위치한 법선 다이폴과 소스에 의해 번째 판에 유기되는 속도 포텐셜로 아래 Eqs.(3)과 (4)를 통해 형상정보만을 가지고 계산할 수 있다.

, (3)

. (4)

발생한 초월공동의 표면에 위치한 법선 다이폴의 세기는 공동표면에서의 접선방향 속도가 일정하고 공동이 시작되는 부분의 기울기가 날개와 접하는 부분의 기울기와 같다는 점착조건에 의해 계산 가능하다. 따라서 적분방정식을 이산화하여 구하고자하는 미지수의 세기는 전체 날개에 분포한 판의 개수와 같다.

. (5)

 (6)

이때 는 공동발생 지점으로부터 번째 판까지의 공동표면의 길이를 나타낸다. 위 식을  행렬 형태로 정리한 연립방정식으로 풀어 계산한 다이폴과 소스의 세기를 이용하여 앞서 언급하였던 공동 표면에서의 오류 항을 수정해가는 반복 계산을 수행하여 수렴 값을 구한다.

본 해석에 사용된 비점성 해석 방법은 기존의 속도를 기저로 한 판요소법에 의한 선형 및 비선형 계산 결과^[5] 뿐만 아니라 실험 결과^[6] 등과의 비교를 통해 상당한 신뢰성을 확보하였다. 따라서 이를 바탕으로 기본적인 2차원 날개단면 NACA 0012 단면과 NACA 66 단면에서 발생하는 부분 공동의 기하학적 특성에 대한 해석을 수행하였다.

2.2 점성 수치해석 방법

FLUENT에서는 해석을 위해 여러 가지 난류와 공동 모델을 제공하고 있는데, 각 모델의 조합에 의한 결과가 계산 결과에 차이를 주므로^[7] 부분 공동 해석에 보다 적합한 조합을 찾고자 하였다.

다양한 공동 모델이 존재하지만 현재 FLUENT에서는 Schnerr & Sauer, Zwart-Berber-Belamri와 Singhal et al. 모델을 제공하고 있다.^[8] 그 중 앞 선 두 모델을 사용할 경우 계산이 빠르고 안정되게 수렴되며, Singhal et al. 모델의 경우 물리적으로 위의 두 모델과 매우 유사하지만, 수치적으로 다소 불안정하고 사용에 어려움이 있다고 알려져 있다. 따라서 본 연구에서는 Schnerr & Sauer 와 Zwart-Berber-Belamri 모델을 사용하였으며, Schnerr Sauer 모델이 2차원 나비어-스톡스 방정식을 기반으로 하는 반면 Zwart-Berber-Belamri 모델은 비정상 3차원 나비어-스톡스 방정식을 기반으로 만들어졌다. 또한 공동현상과 같은 복잡한 다상유동 현상을 계산하기 위해 사용되는 난류 모델로서 Realizable k-ε 과 Reynolds stress 모델이 있으며, 이러한 각각 두 개의 공동 및 난류 모델을 이용하여 총 4가지 조합으로 계산을 수행하였다.

이에 각 모델들의 해석 방법의 차이가 실제 공동 현상에 어떠한 영향을 미치는지 계산 결과를 통해 확인하고자 하였다. 앞선 비점성 해석에서와 같은 모델과 조건으로 상용 프로그램인 ANSYS FLUENT v.13 을 사용하여 부분 공동이 발생한 정상 상태 유동에 대해서 계산을 수행하였다. 표준 벽함수를 사용하여 계산하였으며, 이때의 레이놀즈수, 에 따른 는 100이고 벽으로부터 떨어진 첫 번째 격자점까지의 거리는 이다.

III. 수치해석 결과 및 분석

3.1 비점성 수치해석 결과

비점성 수치해석 방법을 통해 계산되는 결과 값은 아래 Fig. 1과 같다. 계산하고자 하는 날개 단면의 형상 정보와 초기 부분 공동의 길이를 입력하면 경계조건을 만족하는 공동 형상을 반복 계산을 통해 찾는다. 이를 바탕으로 유입 유동의 입사각에 따른 NACA 0012 날개 단면과 NACA 66 단면에서 발생하는 부분공동의 길이의 변화와 그 때의 공동수를 계산한 결과를 Figs. 2와 3에 정리하였다. 동일한 입사각을 갖는 경우, 공동수가 작아질수록 발생하는 초월공동의 길이가 길어짐을 확인할 수 있다.

3.2 점성 수치해석 결과

점성 해석 도메인의 크기는 유입경계가 날개단면 코드길이의 10배, 유출경계는 20배, 외부경계는 10배이다. Figs. 4와 5는 NACA 0012 날개단면에 대하여 Table 1의 공동수 조건에 대한 각 해석 모델 조합에 따라 계산된 부분 공동의 형상 정보를 나타내었다.

이때 발생하는 부분 공동 유동해석에 적합한 해석 모델을 찾기 위해 각각 두 가지 난류와 공동 모델을 조합한 4가지 방법으로 해석을 수행하였다. 그림에서 나타나는 경계선은 부피율 계수(volume fraction factor)가 0.5 인 경우 생성된 부분공동의 형상을 나타낸 것이다.  이 때, 부피율 계수는 물과 수중기로 이루어진 다상유동장에서 수중기가 차지하는 비율을 나타낸다. 그림에 표시된 3개의 선분은 각 각도 조건별로 3개의 공동수 조건에 대한 결과 값을 나타내며, 오른쪽으로 갈수록 계산한 공동수와 압력 조건이 작을 때의 결과를 나타낸다. Figs. 5 ~ 8은 NACA 66 날개단면에 대하여 Table 2에 나타낸 받음각과 공동수 조건에 대해 각각의 계산 모델별 조합에 따라 계산된 부분 공동의 형상 정보를 보여준다. Fig. 8 (c)는 받음각 5°에 대해 Reynolds stress + Schnerr & Sauer 모델을 사용한 결과로 공동수 1.72에서는 부분 공동이 발생하지 않았다.

Figs. 5 ~  8을 비교한 결과 두 종류의 날개단면에 대해서 난류 모델에 따라서는 Figs. 5 ~ 8에 대해서 각각의 (a), (b) 결과와 (c), (d) 결과를 비교해 보면 Realizable k-ε 모델을 사용하여 계산한 결과 (a)와 (b)가 공동의 길이가 더 길고 두껍게 나타나고 있으며, 압력 회복 구간이 크게 나타나지 않아 공동의 전체적인 면적이 더 크게 나타나고 있다. 또한 같은 난류 모델을 사용하였을 때, 공동 모델에 따른 차이는 Figs. 5 ~ 8에 대해서 각각의 (a), (c) 결과와 (b), (d) 결과를 비교해 보면, 공동수가 큰 경우, 즉 가장 바깥쪽 경계의 크기를 비교했을 때는 앞서 난류 모델에 따라 발생하는 차이에 비하여 크지 않으며, 공동수가 낮아질수록 그 차이는 더 크게 차이나는 것을 확인할 수 있다. 공동 수가 큰 경우, 즉 부분 공동의 길이가 짧은 경우 공동 모델의 차이는 날개 표면에서의 공동 마감의 형태에 미미한 차이를 줄 뿐 발생하는 부분공동의 기하학적 특징에는 크게 영향을 미치지 않는 것으로 보인다. 하지만 Figs. 7과 8 같이 NACA66 날개 형상의 경우 받음각이 커질수록 공동수가 큰 경우에도 부분 공동의 형상이 차이가 나타나는 것을 확인할 수 있다. 이는 받음각이 커지면서 각 모델별로 와류 또는 박리되는 유동을 모사하는 차이에 의한 것으로 판단된다.

3.3 수치해석 방법에 따른 결과 비교

앞서 FLUENT를 사용하여 계산한 점성 해석과 동일한 조건에서 비점성 해석을 수행하여 발생한 부분공동의 길이 정보를 비교해보았다. Figs. 9와 10에서는 각각 NACA 0012 와 NACA 66 날개 단면에서 발생하는 부분 공동의 길이 정보를 상사 계수인 에 따라 나타내었다. 그림에서 나타난 실선은 비점성 해석방법을 이용한 계산 결과를 보이며, 비점성 해석결과의 받음각에 따른 경향성을 보이기 위해 추세선을 표시하였다. 받음각 조건은 명도로 차이를 주었고, 각 난류와 공동 모델의 조합을 사각형과 원형 포인트로 나타내었다. 동일한 받음각(같은 명도) 조건하에 계산된 공동수 조건이 같은 경우(x축 값) 각 해석방법에 따라 발생한 공동의 길이를 비교할 수 있다. 위 계산 결과에 따르면 공동 모델에 상관없이 Realizable k-ε 모델을 사용하여 계산한 공동의 형상이 비점성 해석 방법의 결과와 부분 공동의 길이를 유사하게 나타내고 있다. Reynolds stress 모델의 경우 앞서 언급하였듯이 정상 상태의 계산임에도 불구하고 계산 시간에 따라 발생하는 공동의 마감 형태가 일정하지 않고 발생 길이가 Realizable k-ε 모델과 비교하였을 때 다소 작게 계산된다. 또한 공동수가 작은 경우 공동 모델은 발생하는 공동의 길이에는 크게 영향을 주지 않지만 벽과 가까운 부분에서의 공동 마감 형상에 영향을 미치고 있음을 Figs. 5 ~ 8을 통해 확인할 수 있다.

공동수가 큰 경우에도 그 영향이 난류 모델에 의한 결과와 비교하였을 때 미미하며, 발생한 공동의 길이를 비교하였을 때, 공동 모델 중 Zwart-Gerber- Belamri 모델의 계산 결과가 비점성 해석 결과와 더 유사한 형상을 보여준다. 하지만 실제 유동에서의 경우 공동의 뒷부분에서 압력 회복구간이 존재하며 유입 유동의 반대 방향으로 제트유동이 형성되는데, 비점성 해석방법의 경우 이러한 압력 회복 구간의 존재를 무시하고 모든 공동 표면에서 압력이 일정하다는 선형 마감 모형을 채택하여 사용하고 있다. 비점성 해석 방법은 빠른 시간 안에 비교적 정확한 공동의 길이 정보를 계산해주지만, 두께 및 마감 형태에서는 실제 유동 현상의 재현에 한계가 있다. 앞서 언급한 실제 유동에서 관찰되는 재유입제트유동(Re-entrant Jet)은 공동이 마감됨에 따라 그 끝부분에서 정체점이 형성되고 이에 따라 공동 내부와 외부의 힘의 평형 상태가 깨져 상류 방향으로의 제트유동이 형성되는 현상이다.^[9] 이러한 재유입제트유동 현상이 Reynolds stress 난류 모델을 사용하여 계산한 결과에서는 나타나고 있으며, 이러한 관점에서는 실제 유동의 형상과 가깝게 나타나고 있다고 할 수 있다.

Fig. 11에서는 NACA 0012 날개단면에 대해 입사각이 5도 일 경우 각 공동수 조건별로 계산된 결과의 압력 분포를 보여준다. 비점성 해석방법의 경우 앞서 설명하였듯이 공동이 마감되는 지점에서 강제적으로 압력 회복 구간이 주어졌음을 확인할 수 있다. FLUENT를 이용하여 계산한 각 난류와 공동 모델 조합의 결과를 비교해보면 공동수가 낮을수록, 즉 발생한 초월공동의 길이가 짧을수록 두 해석방법에 의한 결과가 유사하게 나타나며, 이때에 FLUENT의 각 모델 조합에 의한 결과도 크게 차이를 보이지 않는다. 하지만 공동수가 커질수록, 즉 발생한 초월공동의 길이가 길어질수록 해석방법에 의해 발생한 공동의 길이가 다소 차이나며, 앞서 언급 하였듯이 난류 모델에 의한 차이가 공동 모델에 따른 차이보다 상당히 크게 나타난다.

해석 결과 점성해석 방법을 사용한 결과 값이 비점성 해석 방법에 의해 계산한 결과보다 발생한 공동의 길이가 짧게 나타나고 있다. 물론 점성력에 의한 저항에 의해 공동의 크기가 더 작게 계산될 것이라고 예측하였으나, 실제 실험 결과와 함께 그 결과가 비교 검토 될 필요가 있다. 점성해석 방법에 의한 결과는 각 모델 조합에 따라 달라지는 경향성에 대해서 파악하였으며, 특정 모델 조합을 선택하기에는 신뢰할만한 실험 결과와의 검증 단계를 거쳐야 할 것으로 보인다.

IV. 결  론

본 연구에서는 저소음 프로펠러 설계단계에서부터 적용 가능한 부분공동 성능 해석법 개발을 목표로 선박용 프로펠러 날개 단면에서 발생하는 부분공동 다상 유동의 점성 및 비점성 수치해석을 수행하고 그 결과를 비교, 평가하고자 하였다. 비점성 해석 방법에 의한 계산 결과와 상용 프로그램을 사용하여 계산한 결과는 발생한 부분공동의 길이 및 마감형태에서 차이를 보이고 있고, 상용 프로그램에서 제공하는 난류 모델과 공동 모델에 따른 부분 공동 해석 결과는 또한 길이와 마감형태가 서로 다르게 나타난다. 어떠한 결과가 실제 형상과 가깝게 나타내는 지는 추후 실험 및 보다 다양한 날개 단면에 대한 해석을 통해 데이터베이스를 구축, 적합한 조합을 도출해내야 할 것으로 판단된다. 본 연구를 기반으로 부분 공동 해석에 있어 점성 효과를 고려한 비점성 해석법을 개발하고, 프로펠러 설계 단계에서 유용하게 적용할 수 있는 수치해석법으로 확장할 계획이다.