해석 대상으로 연구용 선체와 프로펠러를 고려하였으며, 선박해양플랜트연구소에서 수행한 모형 시험과 같은 경계조건을 적용하였다. 전체 계산영역은 대형 캐비테이션 터널 시험부와 동일한 조건으로 구성하였다. Fig. 1에 해석 대상과 영역을 나타내었으며, 관련 제원은 Table 1에 표시하였다. 경계조건은 입구부에 유속 7 m/s, 출구부에 게이지 압력 0 Pa로 설정하여 선체의 운항 속도를 모사하였으며, 추진기의 회전속도는 37.05 rps로 설정하였다. 시험부 벽면, 선체 벽면 그리고 추진기에서의 경계층 유동을 모사하기 위해 no-slip 조건을 설정하였다. 압축성 해석의 경우 출구부에 압축성 파가 반사되지 않도록 비반사 조건을 추가로 설정하였다. 격자 구성은 벽면에서의 경계층 유동을 잘 예측하기 위해 프리즘 격자를 구성하였으며, 추진기로 유입되는 선체의 반류장을 보다 정확히 모사하기 위해 선미 부근에 조밀하게 격자를 구성하였다. Fig. 2에 격자 구성을 나타내었으며, 터널 시험부와 선체의 y⁺는 30, 추진기는 무차원 벽면 거리 변수인 y⁺가 5 이하가 되도록 구성하였다. 총 1400만 개의 격자를 구성하여 수치해석을 수행하였다.

Fig. 1 
Computational domain and geometries of numerical simulation

Fig. 2 
Grids of computational domain, hull and propeller for numerical simulation

이 연구는 상용프로그램인 ANSYS Fluent를 사용하여 균일혼상류모델의 이상유동해석을 수행하였다. 정상상태 RANS해석을 먼저 수행하고, 계산된 정상 유동을 초기조건으로 비정상 DDES 해석을 수행하였다. 식 (1) ~ (5)는 균일혼상류 RANS 지배방정식을 나타내었다.

여기서, 아래첨자 m은 혼상류, l은 액체, υ는 증기를 의미하고, p, ρ, u, μ, α는 각각 압력, 밀도, 속도, 점성계수, 각 상의 체적분율을 나타낸다. ^-는 평균값이며, '는 변동값을 나타낸다. 공동 모델은 기포동역학 Rayleigh-Plesset 방정식을 기반으로 유도한 Schnerr-Sauer 모델을 적용하였으며, 난류 모델은 k-ω shear stress transport(SST)를 적용하였다. k-ω SST모델은 벽면에서는 k-ω 모델을 적용하여 경계층에서의 역압력 구배를 보다 잘 예측하고, 벽면이 아닌 그 외 영역은 k-ω 모델이 적용하여 계산을 수행한다. 선행연구^(15,17)를 통하여 Schnerr-Sauer 모델과 k-ω SST 난류 모델을 사용하여 추진기에서 발생하는 공동과 수중방사소음에 대한 예측결과의 유효성을 확인한 바가 있다. SST기반의 DDES는 벽면 근처의 경계층 내부는 RANS로 계산이 진행되고 경계층 이외의 영역은 large eddy simulation(LES)로 계산이 진행된다. 식 (6) ~ (8)에는 DDES의 이송방정식을 나타내었다.

식 (8)은 소산항의 길이 척도, f_d는 분포 함수, ν_T와 ν는 난류에 의한 동점성 계수, 매질의 동점성 계수, 그리고 k는 von Karman 상수를 의미한다. 그 외 식 (6) ~ (8)에 사용된 기호들은 참고문헌(24)에 정의된 것과 같다. 충분한 샘플링 주파수와 추진기에 의한 비정상 유동을 정확히 모사하기 위해 추진기 0.1° 회전시간인 7.497376 × 10^-6 s로 시간 간격을 설정하였다.

이 연구에서 보다 현실적인 유동 모사와 파수-주파수 분석을 수행하기 위해 압축성 유동해석을 진행하였다. 유체의 압축성 반영은 에너지 방정식과 상태 방정식을 추가로 적용하여 반영하였다. 에너지 방정식은 식 (9)와 같다.

여기서 k_eff는 유효열전도율이고, S_E는 열원을 나타낸다. 물의 압축성은 다음의 상태 방정식인 Tait 방정식을 고려하였다.

여기서 ρ, K는 액체의 밀도와 체적 탄성 계수를 나타낸다. 다음은 증기의 압축성을 고려할 수 있는 이상기체 방정식이다.

식 (12) 기체 방정식을 통해 외부 압력에 의한 증기의 밀도와 체적 변화를 계산할 수 있다. R은 기체상수, M_ω은 분자량 그리고 T는 온도이다. 위의 에너지 방정식과 상태 방정식을 적용하여 압축성 수치해석을 수행하였다.

이 연구에서 소음을 예측하기 위해 Lighthill의 음향사상법을 기반으로 하는 FW-H 방정식을 사용하였으며, 식 (13)에는 FW-H 방정식을 나타내었다.

여기서 x, y는 각각 수음점과 소음원의 위치이며, r은 수음점과 소음원의 거리이다. t, τ는 관찰자 시간과 지연시간이며, τ = t-r/c₀이다. M_r은 소음원에서 수음점 방향으로의 마하수를 나타내며, ρ, u, p는 각각 매질의 밀도와 속도, 압력을 의미한다. 상첨자 ”' ”는 섭동값, 하첨자 “₀”은 평균값을 의미하며, δ_ij는 Kronecker delta 함수, n은 적분면에서의 수직 벡터, c는 음속을 나타낸다. 식 (13)의 오른쪽 첫 번째 항은 시간변화에 따른 유체 체적 기인 소음원인 홀극 소음원을 의미하고 두 번째 항은 유체가 물체의 표면 압력 변화에 의한 소음인 쌍극 소음원, 마지막 항은 난류에 의해 방사되는 사중극자 소음원이다. 일반적으로 유동장 평균 속도가 저 마하수 유동일 때 난류 구조의 기여도는 효과가 미미한 것으로 알려져 있다. 사중극자 항을 무시한 홀극, 쌍극 소음원은 Farassat 1A공식에 의해 다음과 같이 나타낼 수 있다⁽²⁵⁾.

여기서 “^⋅”와 “ret”는 시간 미분과 지연시간을 의미한다. 아래 첨자 n과 r는 수직방향과 수음점 방향을 나타내는 벡터이다.

이 연구에서는 식(16)과 같이 홀극(pT´)과 쌍극자(pL´) 소음원만을 고려하였다.

파수-주파수 분석은 운송체 실내 소음에 대한 외부 유동의 기여도 평가 및 내부 유동에서의 기여도 평가 등 많은 연구에서 다뤄지고 있다^(26~28). 외부 유동으로 인한 벽면에서 형성되는 압력 섭동은 난류 경계층 유동에 의한 비압축성 압력과 외부 소음에 의한 압축성 압력으로 나눌 수 있다. 하지만 압축성 압력인 음파의 경우 대류에 의한 비압축성 압력에 비해 압력변동이 작아 표면 압력 섭동에서 압축성 압력을 바로 분리하기가 어렵다. 이 연구에서는 압축성 압력을 분리하기 위해 선체 선미부의 압력을 이용하여 푸리에 변환인 파수-주파수 변환을 사용하였으며, 시간 정보인 주파수와 공간 정보인 파수의 함수는 다음 식으로 나타낼 수 있다.

여기서 x와 t는 공간과 시간을 의미하며, p, ξ 그리고 τ 는 각각 압력, 공간 지연, 시간 지연을 나타낸다. 파수-주파수 스펙트럼은 주파수와 파수로 표현되며 식 (20)과 같다.

압축성 압력을 3차원 스펙트럼에서 분리하기 위해 선체 선미부 표면을 따라 평행하게 전파하는 음향파의 정보로부터 식 (19) ~ (20)을 정의할 수 있다.

여기서, U₀는 평균 유속이며, M은 평균 유속에 대한 마하수를 나타낸다. 위 식을 통해 디락콘(dirac cone)을 나타낼 수 있다. 선미부 표면에서 음향 파수(acoustic wavenumber) k_a를 통해 디락콘으로 표시되는 영역 내부는 음속보다 빠른 영역인 압축성 압력으로 분리할 수 있으며, k ≤ k_a를 만족한다. 그 이외의 영역은 음속보다 느린 영역으로 k > k_a 이며, 음속이 아닌 영역이기에 비압축성 압력 섭동 영역으로 정의할 수 있다.

먼저, 추진기가 없는 선체(bare hull)에서의 RANS 정상상태 해석을 통해 선체 반류장을 확인하였다. Fig. 3은 선체 반류장을 축방향 유동에 대해 평균 유속으로 무차원화하여 나타내었다. 실험은 예인수조에서 측정된 반류장이며, 수치해석은 터널에서의 반류장이기 때문에 정확한 비교는 어려우나, 실험과 해석 모두 갈고리 형태의 유사한 유동분포를 확인할 수 있다. 반류장 확인 후, 추진기가 설치된 선체의 정상상태 해석을 수행하였고, 정상상태 해석을 초기조건으로 비정상 유동해석을 수행하였다. Fig. 4는 추진기에서 발생하는 공동을 체적 분율 α = 0.1에 대해 가시화하여 실험과 비교한 그림을 나타내었다. 두 결과 모두 터널 시험부 영역에서 수행되었으며, 각각 θ = 10°와 θ = 20°에서 비교하였다. Fig. 5는 선체 선미부의 압력 분포를 나타내었다. 선미부에서는 와류에 의한 압력변동이 명확하게 나타나지 않으며, 방향타 목을 지나면서 와류가 형성되어 비압축성 압력변동이 발생함을 알 수 있다. 이러한 결과를 미루어 볼 때, 선미부에서의 압력은 프로펠러 공동에 의한 압력변동이 지배적임을 알 수 있다.

Fig. 3 
Hull wake-field at propeller plane

Fig. 4 
Comparison of propeller cavitation between experimental result and compressible SST-DDES solver result

Fig. 5 
Temporal pressure perturbation obtained from compressible SST - DDES solver results

유동해석 결과를 입력값으로 FW-H 적분방정식을통해 소음해석을 수행하였다. 투과성 적분면은 Fig. 6과 같이 영역을 선정하여 유동해석과 동일한 시간 간격으로 필요한 정보를 추출하였다. 수음점의 위치는 Fig. 6과 같이 실험과 비교하기 위해 대형 캐비테이션 터널 수중 청음기와 동일한 위치로 설정하여 계산을 수행하였다. 소음 예측결과와 실험과의 유효성을 평가하기 위해 투과성 적분면으로부터 예측된 결과와 터널 시험결과를 1/3옥타브 밴드로 비교하였으며, Fig. 7에 나타내었다. 고주파에서 스펙트럼 경향과 에너지 수준이 잘 일치하는 것을 확인할 수 있으나, 저주파 영역에서는 실험 결과에 비해 소음해석 결과가 낮게 나타나는 결과를 확인할 수 있다.

Fig. 6 
FW-H integral surface and observer location for noise analysis

Fig. 7 
Comparison of predicted one-third octave band spectrum using permeable FW-H in compressible solver and experimental results

Fig. 7의 전체 음압 수준(overall sound pressure level, OASPL)을 Table 2에 표시하였다.

유동해석과 소음해석을 통해 해석의 유효성을 확인하였으며, 유동해석을 바탕으로 압축성 압력을 추출하기 위해 파수-주파수 변환을 수행하였다. 선미부의 압력은 유동해석과 동일한 시간 간격으로 추출되었으며, Fig. 8 (a)에 입력값으로 사용할 선미부의 영역을 나타내었다. 3차원 파수-주파수는 엄밀하게는 사각 평판 2차원 영역에 대해서 적용되기 때문에 이 연구에서는 Fig. 8의 곡면 형상을 1 mm의 등간격의 정렬 격자계 사각 평판으로 근사하였으며, 400 × 500개의 압력 신호를 사용하였다. 파수-주파수를 통해 계산된 파워 스펙트럼 밀도는 Fig. 9 (a)와 같이 파수-주파수 영역에서 나타낼 수 있으며, 파수 분해능은 △k_x = 1/(N_x△x)로 표현이 가능하다. N_x는 압력 신호의 개수이며, △x는 격자간의 거리이다. 디락콘(-c₀+U₀ ≤ ω/k ≤ c₀+U₀) 내부 영역은 압축성 압력이고, 나머지 영역은 비압축성 압력으로 볼 수 있다. 자세히 확인하기 위해 2차원 파워 스펙트럼 밀도 다이어그램(k_x2 = 0)을 Fig. 9 (b)에 나타내었으며, 압축성 압력을 추출하는 기울기의 특성은 |-c₀+U₀|과 |c₀+U₀|이다. 여기서 c₀(= 1500 m/s)는 음속, U₀(= 7 m/s)는 평균 유속이다. 압축성 영역인 |-c₀+U₀|에 큰 레벨의 압력 성분들이 분포함을 확인할 수 있다. 압축성과 비압축성 압력은 명확하게 분리된 영역을 통해 추출할 수 있다. Fig. 10은 디락콘을 통해 분리된 압축성 압력과 비압축성 압력을 각각 역 파수-주파수를 적용하여 시계열 압력으로 계산하였으며, 특정 시간에서의 선미부 압축성 압력과 비압축성 압력을 나타내었다.

Fig. 8 
Configuration of computation domain of wavenumber-frequency transformation

Fig. 9 
Power spectral density diagrams of Wavenumber-frequency transformation

Fig. 10 
Decomposition of the surface pressure field: compressible pressure (top); incompressible pressure (bottom)

다락콘 내부의 압축성 압력과 디락콘 외부의 비압축성 압력의 스펙트럼은 식 (21)을 사용하여 파수영역 전체에 걸쳐 계산하였으며, Fig. 11에 스펙트럼을 나타내었다. 이 연구의 수치해석 조건에서는 압축성 압력이 저주파 영역을 제외한 고주파 영역에서 비압축성 압력에 비해 기여도가 상대적으로 다소 높게 나타남을 볼 수 있다.

Fig. 11 
Comparison of the power spectral density levels by integrating wavenumber