# Numerical Analysis for Supercavitation Characteristics around Underwater Vehicle according to Ventilated Gas Temperature

## Article information

J. KIMS Technol. 2022;25(5):487-500
Publication date (electronic) : 2022 October 05
doi : https://doi.org/10.9766/KIMST.2022.25.5.487
1)School of Mechanical Engineering, Pusan National University, Korea
1)부산대학교 기계공학부
* Corresponding author, E-mail: ndtrong@pusan.ac.kr
Received 2022 January 04; Revised 2022 April 04; Accepted 2022 August 10.

## Trans Abstract

Supercavitation is a phenomenon in which the cavity covers the entire underwater vehicle. The purpose of this paper is to compare and analyze the thermal effect on the cavity characteristics by changing the ventilated gas temperature through computational analysis. For this study, a homogeneous mixture model based on the 3D Navier-Stokes equation was used. As a phase change model, it is its own code considering both pressure change and temperature change. A dimensionless number Tm was presented to analyze the numerical results, and as the Tm increased, the cavity length increased by about 3.6 times and the cavity width by about 3.3 times at 393.15 K compared to room temperature. Analyzing these thermal effects, it was confirmed that rapid heat exchange and heat transfer between the gas phase and the liquid phase occurred at the location where the ventilated gas was sprayed, affecting the cavity characteristics. In addition, it can be confirmed that the initial cavity surface becomes unstable as the ventilated gas temperature increases, and it can be confirmed based on the numerical analysis results that the critical temperature at which the cavity surface becomes unstable is 373.15 K.

## 1. 서 론

공동현상(cavitation)은 액체가 빠른 속도로 운동할 때 국부적인 압력 변화(증기압 이하로 하강)로 액체내에 증기 기포가 발생하는 기화 현상을 의미한다. 발생한 증기 기포가 기계 부품 등의 벽에 닿으면 부식이나 소음, 진동을 발생시켜 기계의 성능 저하나 사고의 위험을 수반한다. 한 가지 예로 응축수 펌프를 들 수 있다. 응축수 펌프는 액체를 압송할 목적으로 설계되어 있으며, 액체 압송용 펌프에 기체가 섞이게 되면 공회전 상태가 되어 압송할 수가 없다. 여기서 응축수의 압력이 포화증기 압력보다 낮아질 때, 재증발 증기가 층을 형성하게 되며, 임펠러의 회전이 불안정하게 되어 발생한 펌프 내의 맥동으로 효율이 저하된다. 또한 증기의 공동은 응축하며 임펠러와 케이싱에 손상을 주며 돌발적으로 파손될 수 있다[1]. 이로 인해 공동현상을 줄이거나 피하는 것이 주된 관심이자 연구 대상이었다. 하지만 상당한 항력 저감효과라는 공동현상의 이점을 활용한 초공동(supercavitation) 연구가 국내외에서 활발히 진행되고 있다.

초공동은 공동이 수중운동체(underwater vehicle) 전체를 덮는 현상을 말한다. 어떤 수중운동체가 매질 안에서 운동을 하면 마찰과 저항이 발생하게 되며, 매질의 밀도가 높을수록 마찰과 저항이 커지게 된다. 매질이 해수라면 공기의 밀도보다 800배 이상 높으며, 잠수함이나 어뢰가 수중에서 고속으로 항주하는 데 제한이 되는 주요한 원인이 된다[2]. 이러한 초공동현상이 최초로 적용된 사례는 러시아에서 개발하여 현재까지도 운용하고 있는 쉬크발(shkval) 어뢰(torpedo)이며, 370 km/h의 속도로 진행한다. 기존 어뢰의 평균 진행속도가 90 km/h 임을 생각한다면 4배 이상 빠르다는 것을 알 수 있다[3]. 수중에서 고속으로 진행하기 위해서는 기존 어뢰의 프로펠러 형식에서 로켓추진기관으로 제시되고 발전되어왔다. 현재 초공동 어뢰에 적용하고 있는 로켓추진기관은 해수 흡입형 로켓추진기관(water breathing ramjet)을 사용하고 있다. 이것은 일반적으로 로켓의 추진제에 포함되는 산화제를 외부 해수를 흡입하여 대체함으로써 추진기관의 추진제 량 대비 발생 가능 추력, 즉 비추력(specific impulse)을 크게 높일 수 있는 로켓 추진기관으로, 공기를 흡입하여 산화제로 사용하는 램젯(ramjet) 추진기관의 수중 형태(hydro ramjet)로 볼 수 있다. 해수 흡입형 로켓 추진기관을 사용할 경우 해수로 인한 로켓엔진 효율 향상, 증가한 추력 성능으로 인해서 일반 로켓 추진기관보다 수중 주행 거리를 크게 늘릴 수 있다[4]. 현재 초공동 어뢰 기술 수준은 독일의 바라쿠다(barracuda)가 800 km/h, 한국은 500 km/h에 위치해 있다[5].

공동이 발달하여 초공동 현상을 이루는 공동을 두 가지로 분류할 수 있는데, 그 첫 번째가 자연 공동(natural cavitation)이다. 자연 공동은 말 그대로 자연적으로 발생하는 공동을 말한다. 자연 공동으로 초공동 상태에 도달하기 위해서는 매우 빠른 속도가 요구되기 때문에 현실적으로 불가능하다. 이때 수중운동체의 앞부분에 인위적인 공동(artificial cavitation) 발생기인 공동 발생기(cavitator)를 장착함으로써 초공동 상태로 도달하는 것을 도와준다. 공동 발생기의 형상은 여러 가지가 있으나 콘(cone), 원판(disk) 형태를 가장 많이 사용한다. 원판 형태의 공동 발생기는 자연 공동이 넓게 생기나 에너지가 많이 소모된다. 반면, 콘 형태는 자연 공동이 좁게 형성되나 에너지는 적게 소모된다[6]. 하지만 아직까지도 높은 속도 조건이 요구되기 때문에 낮은 속도에서도 초공동 상태로 도달하기 위해서 연소 가스나 공기와 같은 비응축가스(non condensable gas)를 분사시켜 환기 공동(ventilated cavitation)을 생성시키고, 이것이 초공동 현상을 이루는 두 번째 공동이다[7]. 환기 공동을 이용하면 자연 공동으로 초공동 상태에 도달하지 못하는 낮은 속도 조건에서 자연 초공동 상태에 도달할 때까지 안정적으로 가속할 수 있다. 그렇기에 초공동 현상에 대한 연구나 어뢰에서 분사 및 환기 공동을 빼놓을 수 없다. 앞서 설명한 러시아의 쉬크발 어뢰의 경우 로켓의 뜨거운 배기를 앞쪽 공동 발생기를 통해 분사하고, 독일의 바라쿠다의 경우 콘 형상의 공동 발생기 뒤쪽에서 분사된다[5].

초공동 현상은 액상(liquid)과 기상(vapor)이 혼재하는 다상(multi-phase) 유동에 해당하므로 이를 고려한 수치해석 기법이 필요하다. 대표적인 공동현상과 같은 다상 유동 해석 방법으로 균일 혼상류 모델(homogeneous mixture model)을 많이 사용한다. 액상과 기상 사이에 열역학적, 기계적, 동역학적 평형상태로 가정한 후 혼상류에 대한 연속방정식, 운동량 방정식 그리고 에너지 방정식을 풀어내기 때문에 액상과 기상에 대해 각각 풀어내는 방법보다 방정식의 개수를 줄일 수 있어 계산 시간을 단축할 수 있다. Kunz 등[8]은 균일 혼상류 모델과 예조건화(pre-conditioning), 이중 시간 전진 기법(dual time stepping method)을 사용하여 수중운동체 부분-초공동 해석기법을 개발하였으나 압축성 효과는 반영하지 않았다. Lindau 등[9]은 압축성과 온도 변화까지 모두 고려한 모델을 개발하였다. 밀도 기반(density-based)의 압축성 유동해석에 있어 가장 큰 어려움은 저속 유동에서 특성 속도와 유동 속도의 차이에 기인한 행렬의 강성(stiffness) 문제이다. 이러한 강성 문제를 해결하기 위해 Choi 등[10]은 속도 기반의 스케일링 변수를 사용한 예조건화 행렬을 지배방정식의 시간항에 곱해주었고, Ha 등[11]은 압력 기반의 새로운 스케일 변수를 제안한 압축성 다상 유동 모델을 개발했다.

초공동 현상에 관한 연구는 주로 실험에 의해 이루어졌으나, 전산해석 기술이 발달함에 따라 수치해석 기법을 통한 연구도 늘어나고 있다. Ahn 등[12]은 공동 터널에서 분사 공기 양에 따른 초공동 실험을 수행했다. Fard 등[13]은 원판 공동 발생기를 활용한 수중운동체 환기 초공동 연구를 실험 및 전산 해석을 동시에 수행했다. Kim 등[14]은 3차원 나비에-스톡스(Navier-Stokes) 방정식에 기반한 자체 전산 해석 프로그램을 활용하여 수중운동체 주위 부분 공동 유동 해석을 수행했고, Ha 등[15]은 완전 압축성 효과를 고려한 수중운동체 주위 자연-분사 공동에 대한 다상 유동 전산 해석을 수행했다. Hoang 등[16]은 자유 유동 온도에 따른 수중운동체 주위 공동 유동 해석을 수행했다. Salari 등[17], Moghimi 등[18], Javadpour 등[19]은 수중 터널에서 30도, 45도, 60도 콘 공동 발생기 형상에 따른 자연/환기 초공동 실험을 수행했다.

제시된 선행연구들은 속도, 공동 발생기, 분사 공기의 양 등에 따른 수중운동체 주위 초공동 실험 및 전산해석을 진행하였으나, 속도와 분사 공기의 양을 무한히 증가시키는 것은 불가능하고 또한 공동 발생기만으로 초공동현상을 달성하는 것에는 한계가 있다. 따라서 상온 분사가스에 대한 초공동 실험 결과를 In-house 코드로 검증을 수행하고, 검증이 완료된 코드를 바탕으로 분사가스 온도를 변화시킴에 따라 공동 길이와 두께에 미치는 열적 효과를 비교 및 분석하는 것이 목적이다. 또한 본 연구에서 제시된 결과는 초공동 어뢰가 다른 변수와 달리 효과적으로 초공동현상을 달성할 수 있는 근거 자료로 제시될 수 있다.

## 2. 지배방정식 및 수치해석 방법

### 2.1 균일 혼상류 모델

본 연구에서는 다상 유동 초공동현상의 수치해석을 위해 나비에-스톡스(Navier-Stokes) 방정식 기반의 균일 혼상류 모델을 사용했다. 여기서 혼상류는 3상으로 액상, 기상 그리고 비응축가스로 구성된다. 균일 혼상류 모델의 연속방정식은 각 상이 혼합된 상태로 가정하여 각 상에 대해서 따로 계산된다. 그에 반해 운동량 보존 방정식과 에너지 보존 방정식은 각각의 상 경계에서 열역학적, 동역학적, 기계적 평형 상태로 가정하여 혼상류에 대해 하나의 방정식만으로 계산된다. 본 연구에서 사용될 연속방정식(continuity equation), 운동량 보존 방정식(momentum conservation equation), 에너지 보존 방정식(energy conservation equation)을 텐서(tensor)로 표현하면 아래와 같다.

Continuity Equation :

(1) t(Ylρm)+xj(.Ylρmuj)=[(m˙p++m˙p)+(m˙T++m˙T)]
(2) t(Yvρm)+xj(Yvρmuj)=[(m˙p++m˙p)+(m˙T++m˙T)]
(3) t(Yngρm)+xj(Yngρmuj)=0

Momentum Conservation Equation :

(4) t(ρmui)+xj(ρmuiuj)=βp¯pxj+τijxi+ρmgi

Energy Conservation Equation :

(5) t(ρmhtβp¯p)+xj(ρmhtui)=xi(uiτijqi)h(m˙T++m˙T)

여기서 하첨자 l, v, ng, m은 각각 액상, 기상, 비응축가스, 혼상(mixture)을 의미하며, Y, α는 질량 분율(mass fraction), 체적 분율(volume fraction)을 의미하고 아래와 같이 정의된다

(6) Yl+Yv+Yng=1
(7) αl+αv+αng=1

또한 체적분율은 질량분율과 아래와 같은 상호관계를 가진다.

(8) αi=ViVt
(9) Yi=mimt=ρiρmαi

여기서 Vi, mi는 각각 i상에 대한 체적과 질량을나타내며, Vt는 혼상류의 총 체적(total volume)을 의미한다. 그리고 P, ρm, u, v, g, τij, ht 는 각각 압력, 혼상 밀도, 카테시안 좌표계(cartesian coordinate)에서의 각 방향 속도, 중력, 전단 응력, 전엔탈피(total enthalpy) 를 의미한다. m˙p+, m˙p+는 공동 모델로서 압력변화로 인한 상변화, m˙T+, m˙T는 비등모델로서 온도변화에 의한 상변화 질량전달을 의미하고, h(m˙T++m˙T) 는 잠열 에 의한 상변화를 의미한다. 또한 각 방향의 열유속, 그리고 전엔탈피는 아래와 같이 표현된다[20].

(10) qi=(μmPrm+μtPrt)CPmTxj
(11) ht=[Ylhl+Yvhv+12(ui2)]

여기서 T는 혼상류의 온도이며, CP, Pr, Re, μm, μt는 정압비열, 프란틀 수, 레이놀즈 수, 혼상류 점성계수, 난류 점성계수이다. 또한 난류모델은 비정상상태 혼상류 유동을 다른 방법에 비해 비교적 잘 예측하고, 직접적으로 계산하지 않기 때문에 격자 선정 등에서 자유로울 수 있는 필터 기반 k — ɛ 난류모델을 적용했다[21]. 균일혼상류 모델에서 사용한 여러 변수들의 보조방정식(auxiliary equation)은 아래와 같다.

(12) ρm=αlρl+αvρv+αngρng
(13) μm=αlμl+αvμv+αngμng
(14) Cpm=YlCpl+YvCpv+YngCpng
(15) Prm=αlPrl+αvPrv+αngPrng
(16) km=Ylkl+Yvkv+Yngkng
(17) hm=Ylhl+Yvhv+Ynghng

이러한 나비에-스톡스 기반 균일 혼상류 모델의 경우 방정식의 수를 줄이는 등 많은 장점이 있지만, 액상의 음속, 기상의 음속, 물과 공기의 체적분율이 0.5인 혼상류의 음속이 매우 큰 차이를 가지고 있기 때문에 수렴성이 크게 저하되는 현상이 발생하게 된다. 이를 해결하기 위해 예조건화(pre-conditioning)기법과 의사 시간(pseudo)을 고려한 이중 시간 전진 기법(dual time stepping method)을 적용하여 수렴성과 정확도를 높였다. 식 (1) ~ (5)의 물리적 영역(physical domain)에서의 표현을 하나의 행렬 벡터로 표현하고 이를 다시 해석 공간에서 일반화된 곡선 좌표계(generalized curvilinear coordinate)로 표현하면 아래와 같은 지배방정식으로 표현할 수 있다.

(18) ΓeQ^t+ΓQ^τ+(E^EV^)ξ+(F^FV^)η+(G^GV^)ζ=S^

여기서 Q^ 는 원시 변수 벡터(primitive variable vector), Ê 와 F^, Ĝ 는 대류 선속 벡터(convective flux vector), ÊVF^V, ĜV 는 점성 선속 벡터(viscous flux vector), Ŝ는 소스 벡터(source vector)를 의미한다. Γe, Γ는 자코비안 행렬(flux Jacobian matrix)과 예조건화 행렬이며 아래와 같이 정의된다.

여기서 자코비안 행렬(flux Jacobian matrix)과 예조건화 행렬에서 tτ는 실제 물리적 시간(physical time)과 의사 시간(pseudo time)을 의미한다. ξ, η, ζ는 해석적 영역(computational domain)에서의 좌표를 나타 낸다. 상첨자 ′는 예조건화 계수를 의미하고, 예조건화 항의 경우 실제 존재하지 않은 시간에 대한 매개변수 이므로 이중 시간 기법을 사용하게 되는데, 이때 τ는 예조건화 행렬에 대한 시간을 담당한다. 아래 식 (21) 은 예조건화에 의한 관계식이며, c와 ć은 각각 실제 물리 시간에 대한 음속, 의사 시간에 대한 음속이다.

(21) ρmp=ρmp+(1c21c2)

또한 실제 물리 시간에 대한 음속과 의사 시간에 대한 음속의 관계식을 보면 아래와 같다. Vp는 혼상 밀도에 대한 압력의 비 Va에 기반한 예조건화 기법의 주요 변수이며, 식 (23)에서 저속 구간의 수렴성 저하를 보정하기 위해서 βp값은 1로 부여했다.

(22) c2=ρmhTρmhTρmp+(1ρmhp)ρmT
(23) Vp2=c2=min[c2,max(Va2,βp2)],Va2=p/ρm

### 2.2 공동 모델

공동 모델은 상변화에 따른 질량 전달량을 수학적 평형을 기반으로 한 모델을 사용했다. 균일 혼상류 모델의 경우는 상변화에 의한 질량 전달을 모사하기 위한 모델링이 필요하다. 이때, 압력 변화에 의한 모델링 결과를 공동 모델이라고 한다. 본 논문에서 앞서 언급한 것처럼 상간 질량 이동을 수학적 평형하다고 가정하여 자체 개발한 공동 모델을 사용했다. 연구에 활용된 공동 모델의 경우 수학적으로 증발(evaporation)과 응축(condensation)을 동일하게 모델링 하여 물리적으로 적합한 모델링이고 아래와 같다[22].

(24) m˙+=kvρναLtmin[1,max(pVpkppV,0)]
(25) m˙=klρVαVtmin[1,max(ppVkppV,0)]

여기서 α, ρ, P는 각각 공극분율(void fraction), 밀도, 국소 압력을, PV는 국소 포화 수증기압, t∞ 는 특 성 시간을 의미한다. 또한 kv, kl, kp는 경험적 스케일 상수 및 질량 전달 강도 계수, 증발 모델의 완전 활성 화를 결정하는 스케일링 변수이고, 본 연구에서는 kv = 1,000, kl = 150, kp = 0.1로 적용했다.

### 2.3 상태 방정식

물과 공기가 혼재하는 압축성 다상 유동에서 수치해석을 수행하기 위해서는 다양한 변수들의 상태 방정식이 요구된다. 본 연구에서는 액상과 기상에 대해서는 IAPWS-IF97[23]을 사용했다. 비응축가스에 대해서는 이상기체 상태방정식(ideal gas state equation)을 사용하였고 아래의 식들과 같이 정의되며, Rng, γng는 각각 기체상수와 비열비를 의미한다.

(26) ρng=pRngT
(27) ρngp=1RngT
(28) ρngT=pRngT2
(29) hng=γngγng1RngT
(30) hngp=0
(31) hngT=γngγng1Rng

### 2.4 수치해석 방법

예조건화기법과 이중시간전진기법이 적용된 식 (18)을 수치해석 하기 위해서 암시적 방법(implicit method)을 사용했다. 실제 물리시간 항은 2차 후방 차분(second order accurate backward difference), 의사 시간 항은 1차 후방 차분(first order accurate backward difference)을 사용하여 풀었으며, 이산화 된 식을 정리하면 아래와 같다.

(32) ΔQ^n+1,k[I+ΔτΓ1(Aξ+Bη+1.5ΓeΔtS)]=ΔτΓ1R^n+1,k

여기서 I는 고유 행렬이며, ADI(alternating direction implict) 기법을 잔차(residual)에 사용하기 위해 아래와 같이 정리했다.

(33) R^n+1,k=Γe3Q^n+1,k4Q^n+Q^n12Δt+[(E^ξ+F^η)(Ev^ξ+Fv^η)+S^]n+1,k

대류 선속 벡터 항은 Taylor 급수 전개를 사용하여 선형화(linearization)를 한 후 Upwind Non MUSCL TVD 기법을 사용하여 이산화를 했다. 또한 점성 선속 벡터 항은 2차 중앙 차분(second order central difference)을 사용하여 이산화를 했다. 위의 MUSCL 기법에서 TVD (total variation diminishing) 조건을 만족시키기 위해 각 격자에서의 기울기를 아래와 같은 limiter를 사용하여 제한했다.

(34) ϕ(r,θ)=max[0,min(θr,(1+r)2,θ)],1θ2

## 3. 수치해석 결과

### 3.1 해석 대상 및 경계조건

Fig. 1은 실험 논문[19]에서 수행한 분사 초공동 실험을 본 연구의 수치해석을 위한 해석 모델 개략도이다. 이 실험은 최대 40 m/s의 속도를 생성할 수 있는 반개방형 수중 터널이고, 공동 발생기는 30도, 45도, 60도 대칭 원뿔을 사용했다. 이 공동 발생기 최대 직경 d = 10 mm이고 원뿔 길이는 원뿔 각도에 따라 결정되며, 본 연구에서는 30도 대칭 원뿔 공동 발생기로 선정했다. 경계조건 및 초기조건의 경우 Table 1과 같다. 여기서 Pressure outlet에 해당하는 변수는 정압의 압력조건이고, 이는 outlet 노드들에 고정값으로 적용했다. 또한 Mass flow rate는 분사구에서 나오는 분사가스의 질량 유량, Room temperature는 수중운동체 주위 물 온도인 초기값, Gas temperature는 수중운동체 분사구 위치에 해당하는 노드에 경계조건으로 각 온도에 대한 값을 고정값으로 적용했다. 질량분율은 수중운동체 주위의 Domain을 AL = 1, 분사되는 가스는 AG = 1로 하여 초기값으로 적용했다.

Analysis target schematic diagram and boundary conditions

Boundary and initial conditions

분사가스 분사에 의한 환기 공동의 경우 공동 수 σ 와 공기공급계수(ventilated gas flux coefficient) Cq 와 같은 무차원 수로 유동 특성을 구분할 수 있으며 다음과 같이 정의된다.

(35) σ=PPc12ρu2
(36) Cq=Qud2

여기서 Pu는 자유 흐름(free-stream)의 압력과 속도이며, Pc 는 환기 공동 내부 압력, Q 는 분사가스 의 체적 유량으로 제시된 질량 유량의 단위 변환을 통해 알 수 있고, d 는 공동 발생기의 직경이다. 식 (35)의 공동 수에서 속도 u가 커지거나, 압력차 P - Pc 가 작아지면 공동 수는 감소한다. 여기서 공동 수가 작아질수록 공동의 길이나 두께가 증가하며, 이는 공동이 수중운동체 전체를 감싸는 현상인 초공동현상으로 된다. 보통 공동 수가 0.1 이하인 경우 발 생한다. 또한 Table 1의 조건들로 환기 공동 수와 공기공급계수를 구한 값을 Table 2에 정리했다. 여기서 수치해석 검증의 중요한 결과비교 기준으로 사용할 공동 길이는 공동이 수중운동체를 덮는 X 방향의 길이이고, 공동 두께는 공동이 수중운동체를 덮는 Y 방향의 길이이다.

Ventilated cavitation number and ventilated gas flux coefficient

### 3.2 격자 민감도 검사

격자의 민감도 검사를 위해 정렬 격자 기반의 Coarse Mesh(GRID1)부터 Fine Mesh(GRID3)까지 구성하여 총 격자수를 Table 3에 기술했다. i 방향과 j 방향의 노드 수를 GRID1부터 GRID3까지 25 % 비율로 증가시켰으며, k 방향의 노드 수는 3으로 고정했다. 벽 근처의 최소 격자 크기 y+는 Clustered Size를 부여해 구성하였는데, 이는 Phan 등[16]의 유입되는 유체의 온도에 따른 수중운동체 자연 공동현상 수치해석에서 경험적으로 얻어진 0.005D를 참고했다. 격자 구성은 Fig. 2Fig. 3과 같이 전체 360도에서 5도 만큼에 해당하는 영역에서 Block을 4개로 나누어서 Multi-block으로 구성하여 3D 축대칭(3D axisymmetric)으로 수치해석을 진행했다. 또한 Table 1의 34 m/s Case에 대해 경계조건과 초기값을 부여했다. 이를 바탕으로 격자 민감도 검사 결과를 Fig. 4에 나타내었고, 시간에 따른 공동 길이를 나타낸다. 공동 길이는 액상의 체적분율이 0.5보다 작을 때 측정되는 것을 기준으로 하였으며, 공동이 끊어지거나 비연속적인 공동의 경우는 측정에서 제외했다. 그래프에 실험에서 측정한 값, 즉 목표 값을 검정색 점선으로 강조했다. GRID1은 다른 GRID에 비해 공동 성장이 강하게 일어났다. 하지만 GRID2, 3에서는 목표 공동 길이의 5 % 이내 오차 값에 도달하였고, 시간에 따른 공동 성장에서 큰 차이가 없었다. 따라서 본 연구에서는 GRID1은 적절하지 않은 것으로 판정되고, GRID2, 3에서 비슷하게 발생하였으므로 시간과 계산 비용을 단축시킬 수 있는 GRID2를 최종적으로 선정해 본 수치해석을 진행했다.

Number of grids for grid sensitivity test

Grid construction

Multi-block configuration

Cavity length sensitivity change for each grating with time

### 3.3 수치해석 검증

격자의 민감도 검사를 통해 얻은 GRID2 격자에서 Table 1의 경계조건 및 초기조건, 총 3개의 Case에 대해서 공동 길이 및 공동 두께를 실험 결과와 정성적, 정량적 비교를 함으로써 본 연구에 사용된 In-house 해석 코드를 검증했다. 공동 길이와 공동 두께 측정 기준은 액상의 체적분율이 0.5보다 작을 때 측정되는 것을 기준으로 했다. 수치해석의 Time Step은 CFL 수에 의해 조절되고 결정되며 본 연구에서는 CFL = 0.5로 적용했다. Fig. 5는 공동이 안정화 된 순간의 공동 길이와 두께 그림이며, 그림에서와 같이 유입 속도가 증가함에 따라 점진적으로 성장하며, 실험과 수치해석 결과가 유사함을 확인했다. 또한 정략적 비교 Fig. 6Fig. 7을 보면 수치해석 검증 Case들이 대체적으로 실험 결과와 일치하는 경향성을 나타냈다. 공동 길이는 실험값과 최대 4.1 % 이내, 공동 두께는 최대 3.7 % 이내의 오차를 보였다.

Qualitative comparison of cavity length as a function of speed

Quantitative comparison of cavity length according to inflow velocity

Quantitative comparison of cavity width according to inflow rate

또한 Fig. 8~10은 재진입 유동이 발생하는 지점에서 압력, 속도 분포를 나타낸 그림이다. 내부의 -X 방향의 속도 성분과 외부의 +X 방향의 속도 성분에 의해 순환이 발생하고, 빠른 속도의 재진입 유동으로 인해 재진입 지점에 압력이 높아지고, -X 방향으로의 속도가 빨라지는 것을 직접 확인할 수 있다. 유선의 방향 역시 재진입 유동 발생 지점에서 공동 내부로 들어가는 -X 방향의 속도가 생성되는 것을 확인할 수 있다. 그리고 유입 속도가 증가함에 따라 재진입 유동이 발생하는 지점이 공동 발생기에서 멀어진다는 것과 발생 지점의 압력이 높아진다는 것 역시 확인할 수 있다. 이는 해석 영역으로 유입되는 속도가 증가될수록 주위 유동장의 속도가 증가함을 볼 수 있고, 각 속도에서 공동을 타고 가속되어 재진입 유동지점 유입된다. 또한 유속이 높으면 더 많은 에너지를 가지고 있는 상태에서 수중운동체 재진입 유동 발생 지점을 더 강하게 타격하기 때문에 압력이 높아진다.

Re-entrant jet at 30 m/s

Re-entrant jet at 34 m/s

Re-entrant jet at 37 m/s

### 3.4 분사가스 온도에 따른 초공동현상

분사가스의 온도에 따른 공동 길이 및 두께 등의 특성을 연구하기 위해서 격자 민감도 검사와 수치해석 검증을 마친 U = 34 m/s의 Case에 대해서 분사가스 온도를 바꿔가면서 공동 특성을 분석했다. 기존 실험에서는 분사가스 온도와 유입 유체(물)의 온도는 상온(20 ℃, 293.15 K)이다. 여기서 분사가스의 온도를 50 ℃ (323.15 K), 80 ℃(353.15 K), 100 ℃(373.15 K), 120 ℃(393.15 K)로 증가시켜가면서 공동 특성을 비교했다.

Fig. 11은 서로 다른 분사가스 온도 조건에 따른 공동 길이 및 두께 등의 특성 변화에 대한 수치해석 결과이다. 분사가스의 온도가 증가하면서 공동의 길이나 두께가 증가하는 것을 확연히 확인할 수 있다. 추가적으로 353.15 K에서 373.15 K으로 넘어가는 과정에서 공동의 급격한 변화가 발생하는 이유는 373.15 K이 물의 비등점과 같기 때문일 것으로 예상된다. 물의 비등점과 같은 온도로 분사하면 초기 분사구 주위에서 액상에서 기상으로 급격한 변화가 일어날 것으로 이는 뒤쪽 공동 생성에 영향을 미칠 것이다. 또한 Fig. 12는 온도에 따라 공동 내부에서 어떤 온도 분포를 가지는지에 대한 온도 분포이며, 레전드의 범위는 상온의 온도(293.15 K)부터 물의 비등점(373.15 K)까지로 지정했다. 여기서 분사가스 온도가 분사구 주변부에만 상당한 온도 변화가 존재하고, 그 이후 공동 내부에서는 상온의 온도로 유지되는 것을 확인할 수 있다. 그 이유는 상온의 물과 각기 다른 분사 온도가 분사된 직후 급격한 열교환 때문에 공기 분사구 주변부에만 온도 변화가 있는 것으로 사료된다. 이를 그래프화한 것이 Fig. 13Fig. 14이며, Y 축은 공동 길이와 두께를 공동 발생기 지름으로 나눠주고 X 축도 마찬가지로 식 (37)처럼 온도에 따른 분사가스를 상온의 분사가스로 나누어줘서 무차원화 했다. 그 결과 상온 온도(293.15 K) 대비 393.15 K에서 공동 길이는 약 3.6배 이상, 공동 두께는 약 3.3배 이상 증가한 것을 확인할 수 있다. Fig. 13의 점선은 해석 영역의 한계 때문에 293.15 K, 323.15 K, 353.15 K 3개의 데이터를 로그 스케일로 보간하여 373.15 K, 393.15 K에서의 공동 길이를 예측한 것이다. 최종적으로 총 4 가지 서로 다른 분사가스에 따른 초공동 길이 및 두께를 정성적으로 비교한 것이 Fig. 15이다.

Comparison of cavity characteristics according to ventilated gas temperature

Comparison of temperature distribution according to ventilated gas temperature

Cavity length change with dimensionless number(Thr)

Cavity width change with dimensionless number(Thr)

Comparison of cavity length and width according to ventilated gas temperature

(37) Thr=Hot TemperatureRoom Temperature=Hot Temperature293.15K

이러한 분사가스 온도에 따른 초공동현상 수치해석 결과를 식 (38)처럼 무차원 수(Tm)를 구하여 추가적으로 공동 길이 및 두께의 변화를 비교했다. 여기서 m˙은 상온(293.15 K) 질량 유량, ρ는 밀도, A 는 분사구 면적, V는 분사 속도를 의미하며 ρ만 온도에 따 른 변수이고 나머지는 고정값이다. Table 4는 온도에 따른 밀도 및 Tm을 정량적으로 구한 값이다. Fig. 16Fig. 17은 무차원 수 Tm에 따른 무차원화 된 공동 길이 및 두께의 변화를 나타낸 그래프이다.

Density and dimensionless number(Tm)

Cavity length change with dimensionless number(Tm)

Cavity width change with dimensionless number(Tm)

(38) Tm=ρAVm˙

### 3.5 분사가스 온도가 공동 특성에 미치는 열적 효과

비응축 가스의 물리적 속성은 분사가스 온도에 따라 변하게 된다. 예를 들어, 온도가 올라감으로써 분자 사이의 거리 변화로 인한 밀도 변화는 환기 공동 유동의 구조 및 수중운동체의 유체역학 특성에 영향을 미치게 된다.

Fig. 18은 수중운동체 벽 부근의 온도 변화 곡선을 나타내는 그림이다. 4 가지 다른 분사가스 온도에서 분사구 주변에 상당한 온도 변화가 존재하게 되는 것을 확인할 수 있다. 분사구 위치 이후에는 벽 근처 온도가 일정하게 유지된다. 이는 기상과 액상의 열교환과 열전달이 굉장히 빠르게 일어나고, 분사가스 온도의 영향을 받는 부위는 분사구 주변의 극히 일부에 불과함을 알 수 있다. 여기서 Zhang 등[24]은 고온의 가스를 하나의 포트를 통해 분사하면서 얻은 실험 결과를 바탕으로 전산해석을 수행하였으며, 초기 2,000 K의 연소 가스가 수중으로 나오게 되면서 급격한 온도 강하가 일어남을 확인할 수 있다. 즉, 고온의 분사가스가 수중에 분사되면 기상과 액상 사이의 빠른 열교환과 열전달을 바탕으로 수중운동체 초기 전면 부분의 공동 생성에 큰 영향을 미치고 이는 뒤쪽 공동 특성에도 영향을 미친다는 것을 알 수 있다.

Change of temperature distribution according to wall distance

분사가스 온도가 공동 특성에 미치는 또 다른 열적 효과의 원인을 분석하기 위해서 4 가지 Case 외에 추가적으로 363.15 K에 대해 추가 해석을 수행했다. 이는 물의 비등점 아래인 353.15 K과 363.15 K에서 어떠한 경향이 일어나는지 알아보기 위함이다. Fig. 19에서처럼 온도가 증가함에 따라 초기 공동 표면이 불안정하게 됨을 확인할 수 있다. Nguyen 등[25]이 수행한 쉬크발 형상의 고온 분사가스와 후미부 연소 가스를 고려한 수치해석에서도 분사가스 온도가 올라감으로써 공동 표면이 불안정해지는 현상을 볼 수 있다. 이는 수중운동체 뒤쪽 공동 특성 성장에 영향을 미치는 것으로 예측할 수 있으며, 본 연구에서는 373.15 K, 393.15 K 경우에서 초기 공동 표면이 불안정함을 확인할 수 있다. 이렇게 공동 표면이 불안정해지는 임계온도는 Fig. 19에서 확인할 수 있듯이 373.15 K이며, 그 아래 온도인 363.15 K에서는 확인할 수 없었다. 이러한 분사가스 온도가 공동 특성에 미치는 열적 효과는 실험이나 추가적인 연구로 더 정확한 물리적 거동의 규명이 필요할 것으로 사료된다.

Initial cavity creation and surface shape comparison by ventilated gas temperature

## 4. 결 론

본 연구에서는 Navier-Stokes 방정식 기반의 균일 혼상류 모델을 사용한 In-house 코드를 활용하여 환기 공동을 포함하는 초공동 및 다상 유동 해석을 수행했고, 이를 바탕으로 초공동현상에 도달하기 위한 다양한 변수 중에서 분사가스 온도를 증가시키면서 공동 길이 및 두께의 변화를 분석하고, 수치해석을 수행했다. 따라서 본 논문의 결과 및 결론을 정리하고 요약하면 아래와 같다.

• (1) 공동 주기 특성에 영향을 주는 재진입 유동에 대해 분석한 결과 재진입 유동으로 인해 공동 내부는 불안정해지고, 내부의 -X 방향의 속도 성분과 외부의 +X 방향의 속도 성분에 의해 순환이 발생하여 와류 진동이 발생하는 것을 확인할 수 있다. 빠른 속도의 재진입 유동으로 인해 재진입 지점에 압력이 높아지고, 이로 인한 와류 및 순환 유동으로 공동 이탈이 발생하게 되어 완전히 떨어져 나가는 공동 붕괴로 이어진다. 이후 -X 방향의 속도 성분이 점차 줄어들면서 공동이 안정화되는 것을 확인할 수 있다.

• (2) X 방향의 해석 영역 한계로 상온 대비 393.15 K (무차원 수 Tm이 증가할수록)에서 약 3.6배 이상 증가할 것으로 예상되며, 이를 데이터 간의 오차가 가장 작은 Log 스케일로 보간한 결과 약 4.8배 증가했다. Y 방향의 해석 영역은 한계를 넘어가지 않았고, 공동 두께가 약 3.3배 증가했다.

• (3) 4 가지 다른 분사가스 온도에서 분사구 주변에 상당한 온도 변화가 존재하게 되고, 분사구 위치 이후에는 벽 근처 온도가 일정하게 유지된다. 여기서 분사가스가 분사되는 위치에서 기상과 액상 사이의 빠른 열교환과 열전달이 일어나고, 이는 수중운동체 초기 전면 부분의 공동 생성에 큰 영향을 미치고 이는 다시 뒤쪽 공동 특성에 영향을 미친다는 것을 알 수 있다.

• (4) 분사가스 온도가 증가함에 따라 초기 공동 표면이 불안정하게 됨을 확인할 수 있다. 이는 수중운동체 뒤쪽 공동 특성 성장에 영향을 미치는 것으로 예측할 수 있으며, 본 연구에서는 373.15 K, 393.15 K 경우에서 초기 공동 표면이 불안정함을 확인할 수 있었다. 이렇게 공동 표면이 불안정해지는 임계온도는 373.15 K이다.

## 후 기

이 과제는 부산대학교 기본연구지원사업(2년)에 의하여 연구되었음.

## References

[1]. . Quangnha T., Lee C.. Cavitation in Fuel Pump with 2D Cascade Modeling. The Korean Society for Aeronautical & Space Sciences 37(5):483–489. 2009;
[2]. . Park W. G., Gang T. J., Jeong C. M.. A Study on the Analysis of the Abnormal Flow Field of the Transcendental Cavity Torpedo: Including Full Compressibility and the Effect of Temperature Change. Bulletin of the Society of Naval Architects of Korea 49(3):13–18. 2012;
[3]. . Kim Y. G., Nah Y. I.. Propulsion Technologies of Supercavitating Rocket Torpedo, Shkval The Korean Society of Propulsion Engineers, pp. 383-387. 2011.
[4]. . Um J. Y., Lim H. A., Jin W. S., Choi J. Y.. Conceptual Design of an Underwater Vehicle Powered by Water-Breathing Ramjet. The Korean Society of Propulsion Engineers 18(4):50–60. 2014;
[5]. . Park J. J.. Understanding of Supercavitation Technology and Analysis of Development Cases Defense & Technology. 152–159. 2020.
[6]. . Kim D. H., Park W. G.. Numerical Analysis of Cavity Shape According to Ventilated Condition and Cavitator AoA The Korean Society of Mechanical Engineers. 381–385. 2020.
[7]. . Kunz R. F., Chyczewski T. S., Boger D. A., Stinebring D. R., Gibeling H. J.. Multi-Phase CFD Analysis of Natural and Ventilated Cavitation about Submerged Bodies 3rd ASME/JSME Joint Fluids Engineering Conference. 1–9. 1999.
[8]. . Kunz R. F.. A preconditioned Navier-Stokes Method for Two-Phase Flows with Application to Cavitation Prediction. Computers & Fluids 29(8):849–875. 2000;
[9]. . Lindau J., Kunz R. F., Venkateswaran S., Merkle C.. Development of a Fully-Compressible Multi-Phase Reynolds-Averaged Navier-Stokes Model 15th AIAA Computational Fluid Dynamics Conference. 2648. 2001.
[10]. . Choi Y. H., Merkle C. L.. The Application of Preconditioning in Viscous Flows. Computational Physis 105(2):207–223. 1993;
[11]. . Ha C. T., Park W. G.. Evaluation of a New Scaling Term in Preconditioning Schemes for Computations of Compressible Cavitating and Ventilated Flows. Ocean Engineering 126:432–466. 2016;
[12]. . Ahn B. K., Kim K. S., Jeong S. W., Yoon H. G.. An Experimental Study on Multi-Injected Artificial Supercavitation. The Society of Naval Architects of Korea 58(1):24–31. 2021;
[13]. . Fard M. P., Rashidi I., Nouri N. M.. Numerical and Experimental Study of a Ventilated Supercavitating Vehicle. Journal of Fluids Engineering 136:101301. 2014;
[14]. . Kim D. H., Paramanantham S. S., Park W. G.. Numerical Analysis of Multi-Phase Flow Around Supercavitating Body at Various Cavitator Angle of Attack and Ventilation Mass Flux. Applied Sciences 10:4228. 2020;
[15]. . Jin M. S., Ha C. T., Park W. G.. Numerical Study on Heat Transfer Effects of Cavitating and Flashing Flows based on Homogeneous Mixture Model. International Journal of Heat and Mass Transfer 109:1068–1083. 2017;
[16]. . Phan T. H., Shin J. G., Nguyen V. T., Duy T. N., Park W. G.. Numerical Analysis of an Unsteady Natural Cavitating Flow Around an Axisymmetric Projectile Under Various Free-Stream Temperature Conditions. International Journal of Heat and Mass Transfer 164:120484. 2020;
[17]. . Salari M., Javadpour S. M., Farahat S.. Experimental Study of Fluid Flow Characteristics Around Conical Cavitators with Natural and Ventilated Cavitations. Journal of Marine Science and Technology 25(5):489–498. 2017;
[18]. . Moghimi M., Nouri N. M., Molavi E.. Experimental Investigation on Supercavitating Flow Over Parabolic Cavitators. Journal of Applied Fluid Mechanics 10(1):95–102. 2017;
[19]. . Javadpour S. M., Farahat S., Ajam H., Salari M., Nezhad A. H.. Experimental and Numerical Study of Ventilated Supercavitation Around a Cone Cavitator. Heat Mass Transfer 53(5):1491–1502. 2016;
[20]. . Baron Fourier J. B. J.. The Analytical Theory of Heat The University Press; 1878.
[21]. . Johansen S. T., Wu J., Shyy W.. Filter-based Unsteady RANS Computational. International Journal of Heat and Fluid Flow 25:10–21. 2004;
[22]. . Ha C. T., Park W. G., Merkle C. L.. Multiphase Flow Analysis of Cylinder Using a New Cavitation Model CAV 2009. 662. 2009.
[23]. . Wagner W., Kretzschmar H. J.. IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam International Steam Tables: Properties of Water and Steam based on the Industrial Formulation IAPWS-IF97. 7–150. 2008.
[24]. . Zhang X., Yu Y., Zhou L.. Numerical Study on the Multiphase Flow Characteristics of Gas Curtain Launch for Underwater Gun. International Journal of Heat and Mass Transfer 134:250–261. 2019;
[25]. . Hwang H. S., Nguyen V. T., Kim D. H., Choi D. G., Park S. H., Park W. G.. Analysis of Supercavity Underwater Vehicle Characterized by High Temperature Exhaust Gas at the Rear KIMST Annual Conference Proceedings. 748–749. 2021.

## Article information Continued

### Fig. 1.

Analysis target schematic diagram and boundary conditions

### Table 1.

Boundary and initial conditions

Velocity Inlet Pressure Outlet Gas and room temperature Mass flow rate
30 m/s 1 bar 293.15 K 1.025455×10−4 kg/s
34 m/s
37 m/s

### Table 2.

Ventilated cavitation number and ventilated gas flux coefficient

Velocity Inlet c Ventilated cavitation number Ventilated gas flux coefficient
30 m/s 0.2 0.02837
34 m/s 0.17 0.02503
37 m/s 0.154 0.023

### Fig. 2.

Grid construction

### Fig. 3.

Multi-block configuration

### Fig. 4.

Cavity length sensitivity change for each grating with time

### Table 3.

Number of grids for grid sensitivity test

Division Total number of grids
GRID1 127,680
GRID2 199,500
GRID3 284,232

### Fig. 5.

Qualitative comparison of cavity length as a function of speed

### Fig. 6.

Quantitative comparison of cavity length according to inflow velocity

### Fig. 7.

Quantitative comparison of cavity width according to inflow rate

### Fig. 8.

Re-entrant jet at 30 m/s

### Fig. 9.

Re-entrant jet at 34 m/s

### Fig. 10.

Re-entrant jet at 37 m/s

### Fig. 11.

Comparison of cavity characteristics according to ventilated gas temperature

### Fig. 12.

Comparison of temperature distribution according to ventilated gas temperature

### Fig. 13.

Cavity length change with dimensionless number(Thr)

### Fig. 14.

Cavity width change with dimensionless number(Thr)

### Fig. 15.

Comparison of cavity length and width according to ventilated gas temperature

### Fig. 16.

Cavity length change with dimensionless number(Tm)

### Fig. 17.

Cavity width change with dimensionless number(Tm)

### Table 4.

Density and dimensionless number(Tm)

Temperature Density Tm
293.15 K 1.205 kg/m3 1.0000
323.15 K 1.097 kg/m3 0.9104
353.15 K 1.000 kg/m3 0.8299
373.15 K 0.946 kg/m3 0.7851
393.15 K 0.898 kg/m3 0.7452

### Fig. 18.

Change of temperature distribution according to wall distance

### Fig. 19.

Initial cavity creation and surface shape comparison by ventilated gas temperature