해석 대상 선형은 저항성능과 관련한 비교검증용 데이터로 널리 활용되고 있는 3,600 TEU급 컨테이너선인 KRISO Container Ship (KCS)이다 (Kim et al., 2001; Hino et al., 2020). Fig. 1과 Table 1은 각각 해석 대상 선형과 주요 치수이다.

Fig. 1 
Lines of KCS

Table 2에 CFD 해석 조건의 프루드 수(Fr=V/gLPP)와 레이놀즈 수(Reynolds number, Re=VLPP/ν)를 나타내었다. Hino et al. (2020)에서는 실선 기준 10 knots에서 26 knots까지의 선속 조건에 대한 저항시험 결과를 CFD 검증 자료로 제공하고 있다. 대상 선형은 Kim et al. (2001)의 연구를 통해 설계속도인 24 knots에 대응되는 Fr = 0.260에서 수면파형을 계측한 결과 또한 공개되어있으므로, 이를 CFD의 파형 결과와 비교하거나 실험 결과를 이용해 직접 파형저항을 구할 수도 있을 것으로 판단하였다.

본 연구에서는 KCS 선형의 저항시험이 수행된 속도 영역 중 조파저항이 유의미하게 식별될 것으로 기대되는 실선 기준 14 knots에서 26 knots까지의 속도영역에 대한 저항 해석을 수행하였다. 선체 자세는 저항시험 중 계측된 값으로 고정하였다. CFD 해석은 수면이 없는 이중모형(double body)에 대해서도 같은 Re 영역에 대해 수행하여 선체 저항에 대한 수면의 영향을 비교하였다.

본 연구에서는 오픈소스 CFD 라이브러리인 OpenFOAM을 사용하여 저항시험 조건의 KCS 모형 주위의 유동장 해석을 수행하였다. 유동의 운동 해석을 위해 식 (1), (2)와 같이 비압축성 유동의 지배방정식인 연속 방정식과 RANS(Reynolds-averaged Navier-Stokes) 방정식을 사용하였다. 그리고 공기-물의 경계면인 자유수면을 고려하기 위해 VOF(volume of fluid) 기법을 식(3)과 같이 적용하였다.

여기서 ρ는 밀도, u¯는 속도벡터, p는 정압을 나타낸다. s는 유체에 가해지는 체적력의 소스항을 나타내는데, 여기에는 중력이 포함된다. α는 각 격자에서 물의 체적 비율을 의미하는데, α는 격자가 물로 완전히 채워진 경우 1, 공기만 채워진 경우 0의 값을 갖는다. 본 연구에서는 α=0.5를 자유수면의 위치로 정하였다. τ는 점성응력텐서를 나타내는데, 비압축성 유동에서 식 (4)와 같이 표현된다.

여기서 S는 평균 변형률 텐서, μ_eff는 유효 점성계수로 유체 점성계수와 난류 점성계수(μ_t)의 합에 해당한다. 난류 점성 계수를 계산하기 위해 RANS 방정식에 기반한 k-ω SST 난류 모델을 사용하였다. 이는 식 (5)와 같다.

여기서 a₁는 상수이고 F₂는 혼성함수이며, 난류 운동에너지(k)와 난류 진동수(ω)는 각각의 수송방정식으로부터 계산하였다 (Yeo et al., 2019).

수치해석을 위한 계산 영역 및 경계 조건은 Fig. 2와 같다. 유동장 해석은 반폭의 공간에 대해 이뤄졌다. 좌표축의 원점은 수선면에서 선체의 길이, 폭방향 중심이며, 선미에서 선수 방향을 +x, 선체 중심에서 좌현 방향을 +y, 중력 가속도의 반대 방향을 +z로 설정하였다.

Fig. 2 
Schematic diagram of the computational domain

계산 영역의 크기는 통상적인 저항 해석을 위한 조건을 참고하여 정했다. 선수에서 입구 경계(inlet)까지의 거리는 수선간길이(L_PP)의 1.0배, 선미에서 출구 경계(outlet)까지의 거리는 L_PP의 2.5배, 중심선에서 폭 방향 바깥 경계까지는 L_PP의 1.5배로 설정하였다.

경계면 중 입구 경계에서의 속도 및 압력 경계조건은 각각 Dirichlet과 Neumann 경계조건을 적용하였으며, 출구 경계에서의 속도 및 압력 경계조건은 각각 Neumann과 Dirichlet 경계조건을 사용하였다. 그리고 선체 표면의 속도 경계조건은 속도값은 0으로 고정된 Dirichlet 경계조건, 압력은 Neumann 경계조건을 적용하였다. 모형에서 멀리 떨어진 상부 경계(top), 하부 경계(bottom), 측면 경계(side)는 Dirichlet 경계조건을 적용하였다. 선체 중심선을 따라 놓인 대칭 중심면(symmetry)에서의 속도 및 압력 경계조건은 대칭 경계조건을 사용하였다. 그리고 자유수면이 없는 이중모형의 해석에서는 수면 위치에 경계면을 설정하고, 대칭 경계조건을 적용하였다,

곡률이 심한 벽면에서 격자 간격을 일정하게 유지하기 위해 선체 주변에서는 5개의 프리즘 층을 생성하였으며, 벽함수를 사용하기 위해 첫 번째 격자점은 30 ≤ y+ ≤ 40을 기준으로 위치하였다. 자유수면의 영향을 정확히 고려하기 위해 자유수면 근처에 격자를 밀집시켰으며, 선박으로부터 멀어질수록 격자의 밀집도를 감소시켰다. 본 연구에서 사용된 격자 수는 1,183,658개이다.

시간과 공간 차분은 2차 정확도의 Crank-Nicolson 기법과 중앙 차분법을 각각 적용하였다. 수치 안정성을 위해 속도, 압력, 난류성분들에 대해 OpenFOAM에서 제공하는 linear-upwind 제한자(linearUpwind limiter)를 반영하였다. 또한, VOF 지배방정식의 이송항은 안정된 자유수면 유동 결과를 얻을 수 있는 total variation diminishing 특성을 지닌 VanLeer 제한자를 사용하였다. 계산 시간 간격에 따른 수치 안정성을 위해 Courant 수가 1보다 작도록 계산 시간 간격을 설정하였으며, 압력과 속도의 연성을 계산을 위해 PIMPLE(Merged PISO-SIMPLE) 알고리즘을 사용하였다.

실험 결과로부터 설계 속도인 Fr = 0.26의 파형저항을 구하고, 이를 CFD의 결과와 비교하였다. 실험, CFD의 조파저항과 파형저항의 비교를 통해 파형저항 해석법의 특성을 파악하고 CFD 해석으로 파형저항의 수식을 검증할 수 있는지를 고찰하였다.

Table 3은 CFD를 통한 전저항, 조파저항 해석과 Hino et al. (2020)의 저항 시험 결과의 비교이다. 전체적으로 CFD를 통해 얻은 전저항계수(C_T)는 실험 대비 2% 내외의 차이를 보였기 때문에, CFD 해석의 세팅은 적절히 이뤄진 것으로 판단하였다.

Fig. 4는 실험과 CFD 해석의 전체 파형 비교이다. 파고는 k₀z로 무차원화 되었다. 실험에서는 파형이 먼 하류까지 유지되는데 반해, CFD 해석에서는 선체에서 멀어질수록 급격히 파형의 크기가 줄어드는 것을 알 수 있다. 이러한 문제는 점성을 고려하지 않은 포텐셜 유동 해석의 조파저항 분석에서도 보고된 바 있으나 (Kim et al., 2000), 그 정도는 본 연구의 유한체적법 해석에 비하면 작은 편이므로 수치감쇠와 난류모델의 영향이 함께 작용하는 것으로 보인다. 뒤에서는 이런 CFD 해석 특성의 영향을 정량적으로 분석하였다.

Fig. 4 
Wave elevation field of experiment and CFD analysis

ITTC에서는 파형의 추출 위치가 L_PP의 21%에서 167% 사이에서 정하도록 한다. 이는 선속에 따른 켈빈파형의 파장 변화를 고려하지 않은 결과이므로, 저속에서는 선체 주위에서만 짧은 파장으로 발달되는 켈빈파형을 추출하기가 어렵다는 문제가 있다. 본 연구에서는 Table 4에 보인대로 설계속도에서 y/L_PP 위치를 선정하고, 이에 대응하는 k₀y값 기준으로 종방향 추출 위치를 선정하였다. k₀y의 형식으로 무차원화된 파형 추출 위치는 선속 조건이 반영되었으므로 낮은 선속에서는 선체에 더 가까운 위치의 파형을 추출하게 된다.

Table 4에 파형 추출 위치를 정리하여 나타내었다. k₀y는 선체 중앙선으로부터의 거리이기 때문에 모형선의 반폭인 y/L_PP = 0.07보다 안쪽에 대해서는 수면파형을 얻을 수 없다. Table 4에서 굵은 글씨로 표시한 구역이 파형 추출이 가능한 영역이다.

ITTC 파형저항 해석법에서는 실험을 통해 계측한 수면파 영역을 3가지로 구분한다. 이를 Fig. 4의 실험 결과 부분에 표시하였다. A 영역은 선체의 영향을 받지 않은 정수면으로, 본 연구에서는 x/L_PP = -0.5인 FP 위치로 정하였다. B영역은 실제 파형 계측이 이뤄지는 영역으로, 파형저항 해석에서 가장 중요한 영역이다.

Fig. 5의 실험, CFD 해석의 파형 비교에서 x/L_PP = 1.0까지는 서로의 파형이 일치하는 것으로 보았다. 그 이후로는 CFD 해석의 파형이 감쇠하여 비교할 수가 없게 되어, 이를 B’영역으로 따로 정의하였다. C 영역은 켈빈 파형이 수조 벽면을 맞고 반사되는 영역으로, 파고 계측 결과가 간섭을 받으므로 ITTC에서는 C영역의 데이터는 잘라내고(truncation) 식(28)의 계수(c₁, c₂, c₃)를 잘 조정하여 B 영역의 파형을 잘 추종하는 단일 파장의 파형을 추정하여 덧붙이도록 하고 있다 (ITTC, 2021b).

Fig. 5 
Comparison of longitudinal wave cuts of CFD analysis (solid lines) and experiments (dotted lines) at Fr = 0.26

CFD에서는 수조 벽면을 구현하지 않을뿐더러 Fig. 5에서 보인 바와 같이 생성된 켈빈파가 이미 B’영역에서 감쇠하고 있어 C 영역의 파형을 식별할 수는 없었다. 따라서 trunctaion 개념을 적용해 B’영역과 C영역의 파고를 외삽하여 얻을 수는 없었다.

실험에서 조파저항계수(C_W)는 전저항계수(C_T)에서 점성저항 성분((1+k)C_F)을 뺀 값으로 정의된다. 형상계수(k)는 ITTC 가이드라인을 따라 0.1 < Fr < 0.2 영역에 대해서 Fr⁴과 전저항 계수의 회귀분석을 통해 구하였다 (ITTC, 2021a). 이렇게 구한 조파저항계수를 실험의 파형저항계수(C_WP)과 비교한 결과를 Fig. 6에 나타내었다. 파형저항계수는 B(실선), B’(점선), C(파선) 영역에 대한 성분을 구분하여 나타내었고, k₀y = 0의 파형저항을 외삽하여 표시하였다.

Fig. 6 
Comparison of C_W and C_WP of the experimental and CFD results at Fr = 0.26

Fig. 6에서 B+B’ 영역의 파형저항은 k₀y가 증가함에도 별 차이를 보이지 않았으므로, 쇄파저항의 영향은 무시할만한 것으로 보인다. B 영역의 파형저항은 선체에 가까워질수록 증가하는 경향이 나왔는데, 이는 선체 주위에서 켈빈파형의 파고가 높게 형성되는 것을 감안하면 타당한 것으로 판단된다. 특이하게도 B와 B+B’의 파형저항 추이를 각각 외삽하여 구한 k₀y = 0에서의 파형저항은 비슷한 값으로 추정되었다. 이는 B영역의 파형저항 해석만으로도 k₀y = 0의 파형저항을 유효하게 추정할 수 있음을 의미한다.

Fig. 5에서 보인 바와 같이, CFD와 실험 모두 B 영역에서는 파형이 유사하였다. 따라서 B 영역의 파형저항은 실험과 CFD에서 비슷한 경향을 보였다. 다만 B’영역은 CFD에서 해석하지 못하였으므로 본 연구에서 사용된 CFD 해석 기법만을 이용해 B+B’영역의 파형저항을 추정할 수는 없었고, CFD 결과의 B 영역 파형저항을 k₀y = 0에서 외삽하고 이를 통해 B+B’영역의 파형저항을 역으로 구하는 접근법은 CFD 해석에서도 유효해 보인다.

C 영역의 파형저항은 특이하게 k₀y가 증가할수록 그 크기가 증가하는 경향을 보였다. 이러한 특성은 선미파에서 발달하는 발산파가 B’과 C 영역의 경계에 가까워짐에 따른 해석의 오차로 보인다. Fig. 4의 파형 결과에서 알 수 있듯, 종방향 파형의 횡방향 추출 위치가 선체에서 멀어질수록 선미 발산파가 B‘과 C 영역의 경계에 가깝게 위치하게 된다. 따라서 Fig. 7의 k₀y = 1.5와 4.5의 파형 비교에서 보이듯 C 영역의 외삽된 파고의 진폭(c₁, c₂)은 k₀y가 증가할수록 크게 나타난다. 만약 실험이 수행된 예인수조의 폭이 더 컸다면 C 영역의 경계가 더 하류로 이동하여 이러한 문제는 발생하지 않았을 것으로 생각된다.

Fig. 7 
Comparison of the longitudinal wave pattern at k₀y = 1.5 and 4.5

실험에서 B+B’영역의 파형저항 추이로부터 k₀y = 0에서 외삽해 추정한 파형저항의 값은 0.249로, 1978 ITTC 방법으로 구한 모형시험의 조파저항에 비하면 55% 정도의 크기를 가졌다. 이는 파형저항 해석식에서는 다루지 않는 비선형성이나 다른 물리현상이 실제 조파현상에 영향을 미침으로 이해할 수 있다. 이런 경향은 일반적인 선형을 이용한 파형저항 연구에서도 지적된 사항이다 (Moran and Landweber, 1972; Kang and Lee, 1981; Kang and Kim, 1989).

4.2장에서는 CFD 해석에서의 파형저항 해석 결과를 Fr 변화와 함께 설명하였다. 4.1장에서 확인한대로 B’, C 영역의 파형저항을 CFD를 통해서는 구할 수 없었기 때문에 B영역의 파형저항 결과로부터 k₀y = 0에서 추정된 파형저항을 Fr 변화와 연관하여 설명하였다.

CFD 해석에서 조파저항으로 볼 수 있는 저항 성분은 실험과 같이 1978 ITTC 방법으로 구한 조파저항계수 외에도 수면을 포함한 모형의 전저항과 수면이 없는 이중모형의 저항 차이로도 구할 수 있다. 선체 표면의 법선방향 압력의 적분값은 점성압력저항 성분을 포함하므로 엄밀한 의미의 조파저항은 아니나, CFD 해석에서 바로 구할 수 있는 저항성분이므로 비교에 포함하였다.

본 연구에서는 조파저항이 수면과 관련한 모든 물리현상으로 인한 것으로 여겨 이중모형과 이상유동 조건의 전저항 차이를 CFD의 조파저항으로 정의하여 Table 3의 오른쪽 열에 나타냈고, CFD의 파형저항도 이 값과 비교하였다. 이렇게 구한 CFD의 조파저항계수는 실험의 1978 ITTC 방법을 통해 구한 조파저항계수보다 큰 값을 보였다.

CFD 해석 결과의 파형을 Fig. 8에 나타내었다. 앞서 설계속도의 파형 비교와 같이, 켈빈파형은 선체 근처에만 한정하여 나타나고 후류로 전파되는 양상은 거의 나타나지 않아서 파형은 B 영역에 한정하여 식별할 수 있었다. 모든 선속의 CFD 해석 결과에서 Table 4의 추출 위치별 파형저항을 구해 Table 3의 가장 오른쪽 열의 CFD를 통해 구한 C_W와 비교해 FIg. 9에 나타냈다.

Fig. 8 
Wave elevation (k_oz) around the model with Fr variation of CFD analyses

Fig. 9 
Comparison of C_WP/C_W of CFD results

낮은 Fr일수록 파형은 선형파에 가까워져 파형저항의 해석 정확도가 향상될 것으로 예상하였으나, Fr = 0.152, 0.195 조건에서는 파형저항에 대해 유의미한 결과를 확인할 수 없었다. 이는 격자 해상도의 한계로 낮은 Fr의 작은 수면파형을 정밀하게 얻어내지 못한 탓으로 판단된다. 그보다 높은 선속에서는 C_W와 C_WP의 비율은 실험의 55%와 비슷한 수준으로 나타났기 때문에 선형 수면파 가정에 근거한 파형저항 해석은 실험에서와 마찬가지로 CFD에서도 조파저항과 차이가 있음을 확인하였다.