본 연구에서는 KVLCC2 선형의 1/100 축척 모형을 대상 선박으로 사용하였다. 방향타(rudder)는 혼(horn)타의 형상이며, 로터세일은 Garenaux et al. (2020)의 실험에서 사용된 풍동 모형을 1/33로 축소하여 적용하였다. 선박과 로터세일의 제원은 각각 Table 1과 Table 2에 정리하였다. 로터세일은 Fig. 1과 같이 선박의 무게중심 x_G로부터 종 방향으로 ±0.3L_PP 위치에 설치하였다. 회전비(spin ratio)는 3.0으로 고정하였으며, 평면도 기준 로터세일 모두 반시계 방향으로 회전한다.

Fig. 1 
Target model

본 연구에서는 Yasukawa and Yoshimura (2015)가 제안한 MMG 모델을 사용하여 로터세일 선박의 3자유도 시뮬레이션을 수행하였다. Fig. 2와 같이 물체고정좌표계의 원점 o는 선체 중앙에 위치하며, 선수 방향을 양의 x축, 우현 방향을 양의 y축으로 정의하였다. 하첨자 E는 지구좌표계 기준값을 나타내며, 초기 선박이 직진하는 상황에서 x축과 X_E축은 일치한다. 선수각 ψ은 X_E축과 x축 사이 각도로 정의하고, 우현 방향 회전을 양의 방향으로 한다. 타각 δ은 방향타 기준 선체 중심선 대비 각도로 정의하며 우현 편향을 양의 방향으로 한다. 실제 풍향 θ_true은 선수를 기준으로 측정하며, 선수풍(head wind)일 때를 0^o로 정의하고 시계 방향으로 증가한다.

Fig. 2 
Coordinate system of a rotor-sail ship

본 연구에서는 전후동요 및 좌우동요, 선수동요 운동을 고려하며 아래의 운동방정식을 계산하였다.

여기서 m_x와 m_y는 전후동요 및 좌우동요에 대한 부가질량(added mass)을 의미하며, J_z는 선수동요에 대한 부가 관성 모멘트(added yaw moment of inertia)를 뜻한다. 하첨자 m은 선체 중앙에서 정의되었음을 나타낸다. 외력항은 아래식과 같이 선체(H) 및 프로펠러(P), 방향타(R), 두 개의 로터세일(Rotor1, Rotor2)의 항들로 구성된다.

로터세일에 의한 외력은 Sauder and Alterskjaer (2022)의 방법을 참고하여 준정상(quasi-steady) 접근법과 스트립 이론 (strip theory)을 이용하여 계산하였다. 아래 식과 같이 먼저 θ_true에 따라 구한 실제 유입풍(true wind) 속도 U_true에서 선박의 속도 U_ship을 뺀 값으로 겉보기 유입풍(apparent wind) 속도 U_A를 계산한다.

로터세일의 각 높이에서 유입풍의 세기가 차이가 나기 때문에 기준 풍속 U_∞에 비균일 유입풍 모델을 곱하여 각 높이에서의 실제 유입풍 속도를 산정하였다. 겉보기 유입풍 속력 U_A의 제곱을 z방향으로 적분하여 겉보기 유입풍 속도 방향에 수직한 성분인 양력과 평행한 성분인 항력을 다음과 같이 계산할 수 있다.

여기서 Hwang (2024)의 대리모델을 이용하여 양력계수 C_L과 항력계수 C_D를 다음과 같이 산정하였다.

이 대리모델은 Marco et al. (2016)의 모델을 참고한 모델로 로터세일의 작동 조건을 고려하여 2≦SR≦4, 4≦AR≦7, 1.5≦DR≦2.5의 범위에서 CFD 단독성능 해석을 수행하고 이를 최소자승법을 이용하여 다항식 형태로 모델링하였다.

양력과 항력을 겉보기 풍향 θ_A을 이용하여 다음과 같이 로터세일에 의한 추력 성분 X_Rotor와 횡력 성분 Y_Rotor를 계산하였다.

유입풍 모델링을 위해 Rotor Fwd 주위에서 실제 풍속 U_true를 CFD 해석 결과로부터 측정하였다. 측정 방향은 Fig. 3과 같이 θ_true와 같은 방향에서 0^o부터 180^o까지 30^o 간격으로 설정하였고, 측정 거리는 Rotor Fwd의 중심으로부터 1.5D부터 4.0D까지 0.5D 간격으로 측정한 유입풍을 비교하였다.

Fig. 3 
Measuring position of true wind speed

Fig. 4에 비교한 바와 같이 1.5D 위치에서는 로터세일의 회전운동으로 인해 유입풍속이 감소하는 것을 확인할 수 있다. 0.5D 간격으로 반경 방향으로 넓혀가며 유입풍속을 확인한 결과, 3.5D에서의 U_true가 4.0D의 유입풍속과 거의 일치하는 것을 확인하였다. 이는 로터세일의 회전유동에 영향을 받지 않는 거리로 판단되어 3.5D에서 측정한 유입풍속 U_true를 유입풍 모델링에 사용하였다.

Fig. 4 
Comparison of wind speed about measurement location (θ_true = 0^o)

각 풍향에 따라 유입풍의 z방향 프로파일을 Fig. 5에 비교하였다. θ_true = 150^o와 θ_true = 180^o의 경우, 상갑판과의 상호간섭이 다른 풍향에 비해 높아 유입풍속 모델링에서 제외하였다.

Fig. 5 
Logistic regression model of inflow wind speed

이후 θ_true = 0^o부터 θ_true = 120^o까지 30^o간격으로 Fwd Rotor의 3.5D 위치에서 측정한 유입풍속을 사용하여 로지스틱 회귀를 통해 아래 식의 매개변수를 도출하였다. 본 연구에서 사용한 유입풍 모델(Wind model)의 매개변수는 Table 3에 정리하였다.

본 연구에서는 상용 소프트웨어인 STAR-CCM+(v23.06)을 사용하여 비압축성 유동에 대한 Reynolds-averaged Navier-Stokes (RANS) 해석을 수행하였다. 자유수면은 volume-of-fluid (VOF) 방법으로 포착하고, 선체의 강체 운동은 중첩 격자 기법으로, 로터세일의 회전 효과는 moving reference frame (MRF) 기법으로 모델링하였다. 일반적인 MRF 기법에서는 비관성좌표계 도입에 따른 가속도를 지배방정식의 외력항으로 처리하므로 별도의 MRF 영역 설정이 요구된다. 반면 로터세일과 같은 원통형 형상은 로터세일 실린더와 톰 디스크 벽면에 회전속도 경계조건(rotating wall boundary condition)을 부여함으로써 비교적 간단하게 회전하는 물체의 모사가 가능하다.

프로펠러의 추력은 계산 효율을 위해 가상 디스크(virtual disk) 방법 중 체적력 프로펠러(body force propeller) 모델을 적용하였다. 체적력 모델은 프로펠러 블레이드의 기하학적 형상을 격자로 구현하지 않고 추력 및 토크를 체적력 소스항으로 운동량 방정식에 직접 부가하는 방식이다(Hough and Ordway, 1964). 반경방향 부하 분포에는 Goldstein optimum 함수를 적용하였으며(Goldstein, 1929), 이를 통해 슬립스트림 내 선회류(swirl) 성분이 시간 평균적으로 모사된다. 가상 디스크 방식은 평균적인 유동 특성만을 고려하기 때문에 방향타 유입류의 비대칭 압력 분포가 간접적으로 반영된다. 난류모델은 realizable k-ε을 적용하고 벽함수로 two layer all y+를 사용하였다. 격자는 trimmer 격자를 사용하였으며 격자 수렴도 해석 시 모든 격자 수준에서 y+는 30-100 범위로 유지하도록 prism layer를 조정하였다. 시간 적분은 암시적 비정상(implicit unsteady) 방법을 적용하였고, 속도와 압력은 SIMPLE 방법을 이용하여 연성하였다.

본 연구에서 적용된 계산영역은 Fig. 6과 같이 선박의 길이(LPP)로 무차원하여 표현하였다. 선수 방향으로 1.5L_PP, 선미 방향으로 파의 반사를 최소화하기 위해 2.5L_PP로 설정하였고, 폭은 1.0L_PP, 깊이는 2.5L_PP 길이로 설정하였다 (ITTC, 2024). 경계조건으로 양측면(side), 입구면(inlet), 출구면(outlet)과 아랫면(bottom)은 속도유입조건(velocity inlet)을 적용하였고, 윗면(top)에는 압력유출조건(pressure outlet)을 적용하였다.

Fig. 6 
Computational domain

3자유도 시뮬레이션에서 로터세일 외력으로 인한 모멘트를 상쇄하기 위해 방향타에 PI 제어를 적용하였다. MMG 모델을 통해 조정된 PI 계수들을 CFD 계산에 동일하게 사용하였고 제어는 아래 식과 같이 선수각에 대한 제어를 적용하였다.

여기서 K_P와 K_I는 각각 비례 및 적분 이득, ψ_target은 목표 선수각, ψ는 현재 선수각이다. 오차를 목표 선수각과 현재 선수각의 차이로 두고, 비례-적분(PI) 제어를 적용하여 목표 조타각 δ_target(t)을 계산하였다. K_P와 K_I는 MMG 모델을 활용하여 선수각 응답의 수렴성과 안정성을 기준으로 시행착오법(trial and error)을 통해 K_P = 5.0, K_I = 0.2로 결정하였다. 실제 조타기 구동계의 물리적 제약을 반영하기 위해 조타각 변화율에 상한을 두었다.

Fig. 7은 본 연구에서 사용한 격자계를 나타내며, Fig. 7(a)와 같이 선박의 Kelvin wave를 모사하기 위해 선체 주위 격자의 해상도를 높였다. 선체 표면으로부터 반사되는 선체 주변 수면파를 예측하기 위해 Kim et al. (2021)을 참고하여 자유수면 부근의 격자를 해상하였고, 추후 6자유도 선박 조종 운동을 해석하기 위해 선체 중첩격자를 생성하였다. 또한 방향타 회전 움직임을 반영하기 위해 방향타 주변에도 중첩격자를 생성하였다. 배경 격자는 선체의 수평면(Surge, Sway, Yaw) 움직임을 따라가도록 설정하였다.

Fig. 7 
Grid configuration of rotor sail ship

격자의 수렴성을 확인하고자 KVLCC2 선박의 자항 시험을 ITTC 권고 절차에 따라 수행하였다. 모든 격자 조건에서 프로펠러 회전수는 17.5 rps로 고정한 후 트림(trim) 및 침하(sinkage)를 고려한 상태의 자항 수렴 선속 결과를 이용하여 격자 수렴도를 확인하였다. 격자수와 수렴도 해석결과는 Table 4에 정리하였으며, 격자 불확실성은 1.57%로 사용한 격자의 수렴성을 나타낸다.

계산의 유효성을 확인하기 위해 로터세일이 없는 상태의 선박을 대상으로 Kim et al. (2019) 모형실험 결과와 자항점을 비교하여 Table 5에 정리하였다. Medium 격자를 사용하여 0.5 rps 간격으로 세 가지 프로펠러 회전수 조건에서 수렴 선속을 측정하였고, 수렴 선속 결과를 최소 자승법(least squares method)을 통해 자항점을 도출하였다. 수렴 선속은 약 0.7979 m/s로 모형시험과 약 0.08% 오차율을 나타냈고, 자항점의 프로펠러 회전수는 2.17% 오차율을 나타냈다.

타력 계산의 정확도를 확인하기 위해 정적 방향타 테스트(static rudder test)를 수행하고 이를 통해 얻어진 유체력 계수를 Yasukawa and Yoshimura (2015)의 모형시험 결과와 비교하였다. 타각 δ를 0°에서 20° 범위에서 변화시키며 선체에 작용하는 종방향력 X' 및 횡방향력 Y', 선수동요 모멘트 N', 타 수직력 F_N'을 비교하여 Fig. 8에 정리하였다. 본 연구의 CFD 해석 결과는 모형시험 결과와 우수한 일치를 보인다.

Fig. 8 
Comparison of CFD and experimental results for static rudder tests

로터세일 단독 성능 테스트를 위해 로터세일 주변 격자의 경우 Hwang (2024)의 설정을 동일하게 구성하였다. 해당 설정에서 양력 및 항력 계수의 격자 불확실성은 각각 1.77%와 0.82%로 나타났으며, 로터세일 단독성능 모형시험 결과(Garenaux et al., 2020)와 양력 및 항력 계수 모두 약 10% 이내의 오차율을 보였다.

1자유도 CFD 해석을 수행하여 로터세일 추진력이 작용하는 조건에서 MMG모델에 로터세일을 결합한 모델을 검증하였다. θ_true = 0^o부터 θ_true = 180^o까지 30^o 간격으로 설정하여 해석하였고, Fig. 9과 같이 로터세일 선박의 선속 수렴 시간인 90초까지 계산을 수행하였다.

Fig. 9 
CFD analysis of rotor-sail ship speed convergence time in a 1DoF simulation

모든 θ_true에 대하여 수렴된 선속을 Fig. 10에 나타내었으며, 그래프의 반경 방향, 각도 방향은 각각 로터세일이 장착된 KVLCC2 선박의 수렴 선속과 유입풍향 θ_true을 나타낸다. 수렴 선속은 로터세일이 장착되지 않은 KVLCC2의 자항선속 U0으로 무차원화 하였으며, θ_true = 0°부터 θ_true = 180°까지의 결과를 대칭시킨 그래프이다. 로터세일 높이에 따라 유입풍속을 균일하게 입력한 MMG 모델과 유입풍 모델이 적용된 MMG 모델을 CFD의 결과와 비교하였다.

Fig. 10 
Converged speed of rotor sail ship in 1DoF simulation (t=90s)

균일풍을 적용한 MMG 모델은 상갑판과 유입풍의 상호작용에 따른 로터세일 성능 저감 효과를 반영하지 못하므로 전반적으로 선속을 과대 예측하였다. 이에 비해 유입풍 모델을 적용한 경우에는 선박의 상갑판 근처 유속 감소 효과가 반영되어 CFD 결과와의 일치도가 전반적으로 유사하였다.

Fig. 11은 길이방향으로 일렬 배치된 두 개의 로터에서 측정한 겉보기 유입풍속 그래프이다. 로터 간 상호작용이 상대적으로 큰 θ_true = 0°및 θ_true = 150^o에서는 Rotor Fwd와 Rotor Aft에서 측정된 겉보기 유입풍속 U_A 차이가 크게 나타났다. 반면, 상호작용이 작은 θ_true = 90^o에서는 두 로터의 U_A가 거의 일치하였다. 이처럼 두 로터세일의 상호작용으로 인해 유입풍속 개형이 달라질 경우 Fig. 9에서 살펴본 바와 같이 선속의 오차가 커질 수 있다.

Fig. 11 
Apparent wind speed U_A (t = 90s)

특히 정면풍인 θ_true = 0^o에서는 Rotor Fwd가 먼저 유입풍을 받고, 그 후류 및 회전 유동의 영향으로 Rotor Aft의 U_A가 감소하였다. 반대로 θ_true = 150^o에서는 Rotor Aft의 회전으로 형성된 유동이 Rotor Fwd 주변 유동에 영향을 미치면서 Rotor Fwd의 U_A가 증가하는 경향이 나타났다. 즉, θ_true = 0^o과 θ_true = 150^o에서는 로터간 상호작용으로 인해 실제 각 로터의 유입 조건이 Rotor Fwd 기준 유입풍 모델과 달라지므로 CFD와 MMG 모델의 결과 차이가 상대적으로 크게 나타난 것으로 해석할 수 있다. 반면 θ_true = 180^o에서는 θ_true = 0° 및 θ_true = 150°에 비해 두 로터 간 U_A 차이가 상대적으로 작게 나타났으며, 이는 선박의 전진 방향과 풍향이 서로 일치하기 때문에 로터 간 상호작용이 비교적 작은 것을 의미한다.

앞선 1자유도 해석을 통해 선속이 수렴된 이후 선수각 오차를 기준으로 방향타 제어를 시작하며, 선수각이 0^o로 수렴하는 시간에 대응하는 150초까지 수평면 3자유도 해석을 수행하였다.

기존의 MMG모델은 방향타 유효 유입각 α_R을 선박 표류각 β에 대한 식을 이용하여 계산하므로 1자유도 조건에서는 α_R이 0으로 예측되어 방향타 모멘트가 발생하지 않는다. 하지만 CFD 해석 결과에서 Fig. 12의 0.5ρU₀²로 무차원화된 동압 분포와 같이 프로펠러의 비대칭 후류로 인해 방향타각 δ = 0^o에서 방향타 모멘트(NR)가 발생한다. 따라서 CFD 결과의 풍향별 평균 방향타 모멘트(NR)를 Yasukawa and Yoshimura (2015)의 모델에 대입하여 식 (10)과 같이 α_R을 미지수로 정리하였다. 이를 통해 α_R = 2.0^o로 산정하여 MMG 모델의 유효 유입각을 결정하였다.

Fig. 12 
Pressure distribution on rudder induced by asymmetric propeller flow (t = 90s, δ = 0°)

여기서 A_R은 방향타의 면적, U_R은 타에 유입되는 축방향 속도 성분, f_α는 타 양력 구배 계수, x_R은 타의 종방향 위치, aH는 타력 증가 계수, x_H는 선체 추가 횡력의 실제 작용 위치를 나타낸다.

Fig. 13는 CFD와 MMG 모델 해석 방법별로 로터세일 선박의 직진 안정성을 확인하기 위해 θ_true = 0^o, 90^o, 180^o 조건에서 유효유입각 보정 전과 보정 후 궤적을 비교한 그래프이다. 그래프에서 원점은 1자유도 해석의 마지막 시점을 기준으로 하였다. 모든 풍향에서 로터세일의 외력으로 발생하는 모멘트로 인해 로터세일 선박의 표류각이 발생하고, 이를 상쇄하기 위해 방향타 제어를 실시한 결과, 선수각은 0^o로 수렴하였으나, 선박은 원래의 직진 선로로 복귀하지 못하고 횡방향으로 이탈한 상태에서 전진하였다. 즉, 본 결과는 선박이 선수각은 유지한 채 일정한 crab angle을 가지고 항주하는 상태임을 의미한다. θ_true = 0°의 경우, 로터세일이 반시계 방향으로 회전할 때 양력이 좌현 방향으로 작용하므로 –y 방향으로 치우치는 것을 확인할 수 있고, 반대로 θ_true = 180°의 경우에는 +y 방향으로 치우치는 것을 확인할 수 있다. 유효유입각이 보정된 MMG 모델의 궤적(실선)이 보정 이전 MMG 모델(굵은 파선)보다 CFD 결과 점선)에 근접한 것을 알 수 있다.

Fig. 13 
Comparison of trajectory in 3DoF simulation

이를 정량적으로 비교하기 위해 3자유도 시뮬레이션의 계산 완료 시점에서 횡방향 편차(cross track error)를 비교하여 Fig. 14에 나타내었다. 유효 유입각 보정 전후 CFD 결과와 횡방향 편차 최대 차이값은 θ_true = 0^o에서 각각 약 0.07 LPP과 0.04 LPP로 나타났으며, 유효유입각 보정 후 CFD 해석 결과에 근접하는 것을 확인 할 수 있다.

Fig. 14 
Cross track error of trajectory in 3DoF simulation

Fig. 15는 θ_true에 대한 타각의 시계열을 비교한 그래프이다. 1자유도 해석에 이어 3자유도 해석을 진행하며 로터세일의 외력으로 인한 표류각을 상쇄하기 위해 ψ_target = 0^o를 제어입력으로 하여 방향타를 제어를 시작하였으며, 모두 동일 PI 제어계수를 적용하였다. 계산 시작 후 10s 이내 오버슛이 발생한 뒤 선수각이 수렴함에 따라 타각도 함께 수렴하는 경향을 보였다.

Fig. 15 
Comparison of rudder angle in time domain

α_R 보정 전에는 3자유도 계산 초기에 프로펠러 비대칭 후류로 인한 유효유입각이 반영되지 않아 중립 타각을 과소 예측하였다. 반면 방향타에 작용하는 모멘트로 α_R을 역산 및 보정한 후에는 CFD와의 일치도가 개선되었고 최대 차이율은 약 40%에서 19%로 감소하였다.

Fig. 16는 3자유도 계산 종료 시점에서 α_R이 보정된 MMG 모델과 CFD 결과의 선체 모멘트(N_H), 방향타 모멘트(N_R), 로터세일 모멘트(N_Rotor)를 무차원하여 비교한 그래프이다. 전반적으로 CFD 해석과 MMG 모델의 모멘트는 서로 잘 일치하는 경향을 보였다. 다만 Fig. 15의 결과를 보면, MMG 모델의 중립 타각이 CFD 해석 결과보다 작게 예측됐지만 타 모멘트는 CFD 결과보다 근소하게 크게 예측되었다. 이러한 차이는 MMG 모델에 적용된 방향타 관련 계수가 CFD가 아닌 모형시험을 통해 도출되었기 때문으로 판단된다. 본 연구에서 확인된 차이는 크지 않지만, 추후 선회나 지그재그 시험과 같이 대각도 타각이 요구되는 조건에서 방향타로 인한 모멘트를 추가적으로 비교 및 검증할 필요가 있다.

Fig. 16 
Comparison of nondimensionalized yaw-moment components (t = 150s)