연구에 사용된 추진기 모형은 KP505이다. 이는 선박해양플랜트연구소(Korea Research Institute of Ships and Ocean Engineering, KRISO)에서 연구목적으로 선형과 실험 결과가 공개된 3,600TEU급 컨테이너선 KRISO Container Ship(KCS)를 위해 설계되었다 (Kim et al., 2001). 추진기의 지름(D)은 대상선형을 활용한 연구를 참고하여 선정하였다 (Kim et al., 2021a; Kim et al., 2021b). 연구에 사용한 KP505의 형상과 주요 제원은 Fig. 1과 Table 1에 나타내었다.

Fig. 1 
Design of KP505 propeller blade

본 연구의 CFD 해석에 선체를 포함하지는 않았지만, 선박의 운동 해석과의 연계를 고려하여 좌표계는 가상의 선체에 고정된 것으로 정하였다. 따라서 추진기의 회전을 따라 좌표계가 회전하지 않는다. Fig. 2는 좌표계의 구성이다. 운동 해석의 좌표계는 원점(O)이 선체의 길이방향 중심을 지나는 수직선과 추진기 축의 연장선의 교점에 위치한 오른손 직교 좌표계(O - xyz)이다.

Fig. 2 
Ship-fixed coordinate system

x, y, z축의 양의 방향은 각각 추진기의 축을 따라 상류 방향, 축에 수직인 우현 방향, 아래 방향으로 정의하였다. 추진기는 x축을 중심으로 회전하며 x와 z방향으로 전진하도록 하였다. 추진기 날개의 위상(ϕ)은 날개의 기준선과 z축이 이루는 각도로, x축을 중심으로 하는 회전 방향 상에서 정의된다. 선수파 중 선박 운동을 고려한다면 z축은 선체의 파랑 중 종동요, 상하동요의 영향을 직접 받는 방향, y축은 운동이 없는 방향으로 볼 수 있다. 따라서 본 논문에서는 z방향 운동을 주된 시험 조건으로 정의하였다.

실제 선박의 선수파 중 운동은 x - z평면 상에서 전후동요, 상하동요, 종동요로 이뤄진다. 추진기의 운동 또한 이 세 운동이 복합적으로 나타날 것으로 기대된다. 하지만 추진기 성능에 대해 세 가지 운동의 속도와 가속도의 복합적 영향을 전체적으로 이해하기는 어려운 일이므로, 해석에 적용할 운동 성분을 줄여 문제를 단순화할 필요가 있다.

가속도를 포함한 해석 조건을 정하는데 있어 Park et al. (2020)의 연구를 참고하여 부가질량의 개념으로 추진기의 영향을 가정하였다. 전체 모형선의 원점 기준 부가질량(m_z)과 부가관성모멘트(J_yy)는 각각 식 (1)과 (2)와 같이 나타낼 수 있다. 여기서 x_P는 선체 길이 방향 중심에서 본 추진기의 x방향 위치이다.

만약 선체가 상하동요와 종동요 운동을 한다면 선체의 상하동요 운동에 대한 부가질량(m_zH)과 부가관성모멘트(J_yyH)에 더하여 추진기의 상하동요의 부가질량(m_zP)과 종동요의 부가관성모멘트(J_yyP)가 발생할 것이다. 식 (1)의 우변의 두 성분에 대해서, 선체의 상하동요 중 투영되는 면적이 추진기보다 월등히 크기 때문에 추진기의 부가질량은 무시할만한 수준으로 생각된다. 같은 이유로 식 (2)의 J_yyP 또한 J_yyH에 비할 수준은 아닌 것으로 판단하였다.

반면 식 (2)에서는 추진기가 무게중심에서 가장 먼 위치인 선체의 후방에 위치하면서 xP2mzP항이 전체 J_yy에 영향을 줄 수 있을 것으로 보았다. 따라서 추진기의 상하동요 병진 운동으로 인한 가속도 관련 항이 선체의 종동요 운동에 미치는 영향을 식별하는 방향으로 해석 조건을 한정하였다.

추진기의 작동 조건은 전진 속도(U_A), 회전수(n), 추진기의 축에 수직인 종방향의 운동 속도(W), 같은 방향의 가속도(W˙)의 네 가지로 구성된다. 앞의 두 가지는 추진기 성능과 관련하여 필수적으로 사용되는 전진비(advance ratio, J_A)로 무차원화하여 식 (3)과 같이 나타내었다.

W는 사항각과의 연계를 고려하여 식 (4)와 같이 U_A로 무차원화하여 W'로 나타내었다. 본 연구에서 다룰 시험 조건에서는 W'가 크지 않기 때문에 W'는 추진기 사항각과 거의 비례한다고 볼 수 있다.

W˙의 무차원화를 위해 W˙로 인한 추진기의 유체력이 어떤 양상을 보일지를 먼저 고려해야 한다. 추진기는 양력이 지배적이므로 유체력에 대해 점성의 영향이 작은 편이기 때문에, 해석 조건의 선정에서는 비점성 유동을 가정할 수 있다. 따라서 부가질량의 관점에서 가속도 영향을 받는 성분의 크기는 U_A, E에는 무관하고 W˙에 비례할 것으로 예상하여, 실험 조건 전체에서 고정된 값인 n, D를 이용해 식 (5)와 같이 W˙'으로 무차원화하였다.

Table 2에 전체 해석 조건을 정리하여 나타내었다. 시험은 정속 직진(straightforward), 정속 사항(steady drift), 수직 가속(vertical acceleration)의 세 가지로 구성된다. 정속 직진 조건은 같은 조건의 모형시험과 비교검증을 위해 다양한 J_A 조건에서 해석을 수행하였다. 사항 조건의 J_A 는 본 연구가 참고한 자항시험(Kim et al., 2021b)의 조건에서 사용한 모형선 자항 상태에서 추진기의 작동 조건을 바탕으로 구하였다. 모형선의 전진속도(V_S), 반류비(w_P0), n가 주어진다면 식 (6)을 통해 J_A 를 구할 수 있다.

KP505의 상하 운동과 가속도는 추진기의 장착 선형인 KCS의 선수파 중 부가동력에 대한 CFD 해석 결과를 이용해 그 범위를 구하였다 (Kim et al., 2021a). 파장 별로 얻어진 선체의 운동과 추진기의 수직방향 변위로부터 W와 W˙의 최대값을 구했고, 이를 포함하도록 CFD 해석 조건의 범위를 정하였다.

정속 사항 조건의 해석이 수행된 W'의 최대값은 0.38로, 이를 입사각으로 환산하면 20.8°이다. 정속 사항 조건의 해석은 최대 W'에 대해 25, 50, 75, 100% 조건에서 수행하였다.

Fig. 3은 제안된 수직 방향 가속 운동을 적용한 시계열 속도, 가속도의 예시이다. 가속도 조건의 해석은 다음의 순서로 수행되었다. 우선 J_A = 0.535 조건에서 추진기를 직진한다. 그리고 시험 조건으로 정해진 0.026, 0.052의 W˙'를 적용하여 W'이 0.6이 될 때까지 수직 방향 속도를 증가시킨다. 이후 -W˙'를 적용해 W'이 0까지 감소하도록 한다. 운동 주기 (T)동안 하나의 W' 조건에 대해 +W˙', -W˙'의 두 가지 가속도 조건의 유체력 변화를 확인할 수 있다.

Fig. 3 
Time history of W˙ and W˙' in periodic vertical motion (W˙'=±0.026)

2.3장에서 정한 해석 조건에서 추진기의 중심인 원점 위치에 작용하는 3방향 유체력(F = [X, Y, Z])과 3방향 모멘트(M = [K, M, N])를 얻는다. 일반적인 추진기의 성능 해석과 같은 직진 조건에서는 축방향 추력(T)과 축방향 토크(Q)를 구하여 성능 해석에 활용한다. 본 연구에서는 6분력을 모두 다루는 관계로 변수 표현의 일관성을 위해 T와 Q 대신 X와 K를 이용하여 해석 결과를 표시하였다. 힘과 모멘트 차원은 각각 식 (7), (8)과 같이 X', K'로 무차원화하였다.

X와 K를 제외한 다른 방향의 유체력과 모멘트도 같은 방식으로 무차원화하였다. 이를 무차원화한 유체력, 모멘트 벡터는 각각 F' = [X', Y', Z'], M' = [K', M', N']로 표시하기로 한다. M', N'은 프로펠러 중앙 위치에서의 무차원 모멘트값인 M_P', N_P'에 식 (9), (10)을 이용해 얻은 선체 중심에서의 값인 M_O', N_O'를 더한 결과를 함께 보였다. K'에 대해서는 K^P'와 K^O'의 차이가 없으므로 본문에서 K^P'로만 나타냈다.

수직 가속도를 포함한 운동 조건의 유체력과 모멘트는 식 (11), (12)와 같이 정속 사항(F'(J_A, W'), M'(J_A, W'))과 가속도 관련 항(ΔF'W˙', ΔM'W˙')의 선형 결합으로 가정하였다.

가속 운동 중 추진기의 유체력 해석을 위해 전산유체역학 프로그램인 STAR-CCM+ V2302(Siemens Digital Industries Software, Plano, Texas)를 이용하였다. 해당 프로그램은 유한체적법에 기반한 비정상 레이놀즈 평균 나비에-스토크스 방정식(Unsteady Reynolds-averaged Navier-Stokes equation)을 비정렬격자계에서 이산화하고 난류모델은 k - ω SST 모델을 사용한다. 대류항은 2차 Upwind 방식으로 이산화되며 속도-압력 연성은 SIMPLE 알고리즘을 사용한다. 시간항에 대한 처리의 경우 타임스텝 당 시간간격이 작아 1차 Euler 방식을 사용하였다. 해석에서 시간 간격은 프로펠러 회전수와 해석 격자계의 크기를 고려하여 직진 조건의 경우 0.0005초, 정속사항/수직가속 조건의 경우 0.00025초로 두었다. 이에 대응하는 추진기 위상 변화는 각각 3.6°, 1.8°이다.

본 연구에서는 추진기의 회전, 축방향의 직진운동에 더하여 축에 수직인 방향의 가속 운동을 구현해야 한다. 해석 영역은 크게 배경 유체 영역과 회전운동을 하는 추진기 영역으로 구성된다. Fig. 4는 계산 영역의 구성이다.

Fig. 4 
Computational domain and grid regions

격자계는 영역별로 배경(Background), 격자 재분할 구역(Grid refinement region), 추진기 원반(Propeller disk)으로 나뉘며 각각 구성된 격자수는 Table 3에서 보인 바와 같다. 추진기 원반 영역은 슬라이딩 메쉬(sliding-mesh)로 회전을 구현하였다.

경계 조건은 Table 2에서 보인 직진과 정속 사항/수직 가속 조건에 대해 각각 적용되었다. 직진 조건에서는 Fig. 4의 유입면(Inlet)에 속도 및 압력 경계조건은 각각 Dirichlet과 Neumann 경계조건을 적용하였고, 유출면(Outlet)에는 속도에 대해 Neumann, 압력에 대해 Dirichlet 경계조건을 적용하였다. 배경 영역의 나머지 경계면으로의 유동 유입이나 유출은 고려하지 않았다. 추진기 표면의 속도 경계조건은 속도값은 0으로 고정된 Dirichlet 경계조건, 압력은 Neumann 경계조건을 적용하였다. 정속 사항과 수직 가속 조건에서는 배경 격자를 포함한 전 계산 영역이 주어진 속도나 가속도로 수직운동을 하도록 하여, 유입면과 유출면을 상하부에도 정의하였다.

Fig. 5는 해석을 위한 프로펠러의 표면과 x - z평면 상의 격자 구성이다. 벽면에서의 전단력(τ_w) 보정을 위한 벽함수(Wall function)를 사용하기 위해 추진기 날개면과 허브에는 9개 층의 프리즘 격자를 배치하였다. 벽으로부터 법선 방향의 첫 번째 격자 높이(Δs)는 식 (13)을 통해 y⁺가 30보다 작은 값을 가지도록 설정하였다.(Versteeg and Malalasekera, 2007)

Fig. 5 
Grid distribution on the model surface and x - z plane

여기서 ρ와 ν는 각각 유체의 밀도와 동점성 계수이다.

CFD를 이용한 정속 직진 중의 추진기 성능 해석 결과를 실험과 비교하여 검증하였다. Fig. 6은 같은 조건의 실험을 통해 얻은 추진기 단독 성능 곡선과의 비교이다. 추진기 효율(η_O)은 식 (14)를 통해 정의하였다. 토크에서 약간의 차이는 있으나, 본 연구에서 주요하게 다룰 추력에서는 실험 결과와 잘 일치하는 것을 확인할 수 있다.

Fig. 6 
Comparison of propeller performance for straight forward test

수직 방향 운동에 따른 유체력과 모멘트의 변화를 살펴보기에 앞서 Dubbioso et al. (2013, 2014)를 참고하여 Fig. 7과 같이 KP505 프로펠러의 2차원 날개 단면에서 사항각 중 받음각 변화를 분석하였다.

Fig. 7 
Schematic of the angle of attack of the propeller blade in oblique flow

추진기 날개의 받음각(α)은 식 (15)를 통해 구할 수 있다.

여기서 β, θ는 각각 날개 단면을 포함한 평면에서 추진기 유입류가 이루는 각도와 추진기 날개의 피치각으로, 식 (16), (17)을 통해 반지름 위치(r), 추진기 피치(P)로부터 계산할 수 있다.

Fig. 8은 본 연구의 정속 사항 시험 조건에서 구한 ϕ에 따른 α의 추세로, r = 0.35D 위치에서 구한 값이다. 프로펠러 날개의 중앙선이 +z축에 놓였을 때를 ϕ = 0˚로 정의하여 나타냈다.

Fig. 8 
Angle of attack of the propeller blade with respect to the phase angle

사항 유입류의 영향으로 α는 90˚, 270˚ 위상에서 서로 반대로 나타났다. K' = 0.38에서 그 진폭은 1.3˚로, 이는 K' = 0조건에서 J_A가 0.48에서 0.59까지 ±10% 정도 변화하는 수준의 α 변화에 대응한다. 이런 추진기 날개의 위상에 따른 하중 변화는 유체력의 불균형을 유발하여 진동의 원인이 될 수 있다. 하지만 평균 α는 W에 무관하게 일정하므로 각 날개의 유체력의 총합인 추진기의 추력, 토크는 W 조건에 대해서 큰 변화가 없을 것으로 생각된다.

모형선 자항 조건인 J_A = 0.535에서 W' 의 변화에 따른 6방향 유체력과 모멘트의 변화를 조사하였다. x축 방향 외의 유체력 성분 중에는 수직 방향 유입류의 방향인 z방향의 유체력(Z)과 x - z평면 상의 모멘트인 y방향 모멘트(M)가 주로 발생할 것으로 예상하였다. 본 연구의 해석 범위와 유사한 30° 이내의 입사각에서의 선형적인 유체력과 모멘트의 증가 양상을 보인 Dubbioso et al. (2013)의 연구 결과에서는 사항 유입류 방향의 유체력 성분이 그렇지 않은 방향 성분에 비해 약 5배 더 큰 것으로 확인되었고, 모멘트는 두 방향이 서로 비슷한 크기를 보였다.

Fig. 9는 정속 사항 조건의 유체력, 모멘트 결과이다. 모든 유체력과 모멘트는 J_A = 0.535, W' = 0 조건의 X₀', K_P0'로써 무차원화하여 나타내었다. Dubbioso et al. (2013)의 연구에서 보인 바와 같이 W'에 따른 Z'의 변화는 운동 방향과 반대로 나타났다. 그리고 그 크기는 운동이 없는 방향인 Y'보다 크게 나타났지만, 그 비율은 2.5배 정도로 Dubbioso et al. (2013)의 연구보다는 두 방향 성분 간의 차이가 작게 얻어졌다. 이는 본 연구에서 사용된 추진기의 형상이 스큐(skew)를 가져 더 복잡한 3차원 형상을 가지기 때문에 z방향 유입류도 날개 표면을 따라 y방향으로 흘러나가는 비율이 증가함에 따른 반작용으로 생각된다.

Fig. 9 
Hydrodynamic force and moment of the propeller with the vertical inflow speed variation

유체 모멘트 성분인 M_P'와 N_P'의 크기가 비슷한 것은 Dubbioso et al. (2013)의 연구에서도 지적되었던 내용으로, 유체력 결과와는 다르게 본 연구에서도 이와 같은 경향을 보였다. 좌표계 원점인 선체의 길이방향 중심에 작용하는 모멘트인 M_O'와 N_O'에 대해서는 추진기 자체에서 발생하는 모멘트인 M_P'와 N_P'보다 추진기의 유체력과 추진기 선미 배치로 인한 모멘트 성분(x_P × F)이 더 지배적인 것으로 확인되었다. 그리고 M_O'이 추진기의 수직방향 속도와 반대로 발생하기 때문에 선미 위치에서 작동하는 추진기는 종동요에 대한 일종의 감쇠 작용(Damping)을 수행하는 것으로 볼 수 있다.

Fig. 8에서 설명한대로 W'의 변화에 따른 X', K_P'의 변화는 그렇게 두드러지지 않았다. 이는 선체의 상하동요, 종동요 운동이 발생하더라도 추진기의 추력, 토크는 전진 속도만을 기준으로 산정해도 큰 문제가 없음을 의미한다. 따라서 이전 연구에서 사용해 온 전진비에만 의존한 추진기 모델이 사항각을 포함한 유동 조건에 놓이더라도 x방향 유체력의 해석에 한해서는 타당한 결과를 낸다고 볼 수 있다.

Fig. 9에서 보인 결과를 이용해 J_A = 0.535 조건에서 W'에 따른 F', M_P'을 식 (18)과 같이 근사하여 나타낼 수 있었다. X_P', K_P'은 추세의 개형을 감안하여 W'의 2차 함수로, 나머지 유체력과 모멘트 성분은 W'의 1차 함수로 근사하였다.

Fig. 10은 W˙'=0.026, 0.052 조건의 수직 가속도 운동 중 얻은 W'와 F', M_P'의 시계열로, 식 (18)을 이용해 추정한 W˙'=0 조건의 W' 에 따른 유체력과 모멘트도 함께 나타내었다. 시계열 자료에서는 추진기 날개의 위상 변화로 인한 짧은 주기의 유체력과 모멘트의 요동이 나타났으나, 이는 크기가 작아 가속도 조건의 영향 분석에서 무시하였다. 전체적으로 유체력의 변화는 W'의 변화를 따르는 경향을 보였는데, X'와 K'는 4.1장에서 다룬 대로 W'에 대해 선형적인 반응을 보이지는 않았다.

Fig. 10 
Time series of hydrodynamic force and moment of vertical acceleration test

Fig. 10에서 실제 결과와 식 (18)을 이용해 추정한 결과의 차이가 식 (11), (12)의 ΔF'W˙', ΔMP'W˙'이다. Fig. 11에 이를 추세선과 함께 나타내었다. 앞에서 언급한 유체력과 모멘트의 위상에 따른 요동으로 인해 데이터가 산포되었으나, 전체적인 경향의 확인에는 문제가 없었다.

Fig. 11 
Off-axial components of ΔF' and ΔM' with variation of acceleration (symbol: mean, error bar: 95% confidence interval)

W˙'에 비례하여 나타나는 유체력 성분을 Y', Z', M', N'에 대해서 잘 확인할 수 있었으므로, 부가질량의 개념으로 이를 설명하고자 하는 접근이 타당함을 확인할 수 있었다. Y', Z', M', N'의 네 가지 운동 성분에 대해 가속도의 영향을 받는 성분과 W˙'와의 비례상수를 Table 4에 나타내었다. Y', Z'간의 크기 차이는 약 4배로 얻어져 정속 사항 중 유체력의 추세와 마찬가지로 z방향 유체력의 변화가 지배적임을 알 수 있다.

그 중 수직 방향 운동에 대해 직접적인 영향을 받는 Z'에 대한 부가질량을 해석하였다. 식 (19)와 같이 Table 4의 비례상수에 차원을 적용하면 종동요 가속에 대한 부가질량 m_z를 구할 수 있다.

이를 식 (2)에 대입하면 선체의 가속 운동에서 추진기가 미치는 영향을 구할 수 있다. 식 (2)의 xP2mzP 항은 0.126 kgm²로 계산되었다. 이는 본 연구의 대상 선형인 KCS의 종동요 방향 관성 모멘트(I_yy) (Kim et al., 2021b)인 38.58 kgm²의 0.3% 정도로, 운동에 영향을 줄 만한 크기는 아닌 것으로 판단된다.

가속도로 인한 유체력 변화를 확인하기 위해 Fig. 12에 추진기의 표면 압력 분포를 정리하여 나타내었다. 표면 압력은 식 (20)에서 보인 바와 같이 nD를 이용해 무차원화되었으며, 식 (21)과 같이 직진 상태를 기준으로 사항각의 영향 성분, 가속도 성분을 분리하여 나타냈다.

Fig. 12 
Composition of pressure coefficient (p') on the propeller surface in heave motion

압력 성분 중 사항각에 대한 p'(W')은 위상에 따른 변화가 앞날에서 명확히 드러났으며, 이는 Fig. 8에서 추진기 날개의 사항 중 받음각 변화에 대해 설명한 내용을 뒷받침한다. 위상이 0° - 180°일 때에서는 추진기 날개의 받음각 증가로 인해 압력면의 압력이 더 높아지고 흡입면의 압력은 낮아지는 양상이 발생했다.

p'W˙'은 +W˙', -W˙'에 대해 각각 나타내었다. 가속도 영향에 대해, 앞서 전체 유체력과 모멘트에서는 가속도에 비례한 유체력 변화가 얻어졌으므로 가속도의 방향에 따라 p'W˙'이 반대 부호로 얻어질 것으로 기대하였다. 하지만 국부 압력 분포에서는 이러한 경향이 나타나지 않았으며, 가속도로 인한 압력의 변화는 주로 90° 내외의 위상에서 허브 근처 압력의 감소로 얻어졌다. 이는 Fig. 7과 식 (15)에서 알 수 있듯, 허브 근처는 r이 작아지면서 추진기 날개의 받음각에 W'가 미치는 영향이 증가하기 때문에 가속 운동에서 유선의 급격한 변화는 허브 근처에서 발생하는 것을 의미한다. 추가로, 유동의 받음각이 바뀌면서 정체점의 위치가 변하게 되는 앞날에서도 압력의 국부적인 변화가 나타났다.