1. 서 론

전산유체역학(computational fluid dynamics, CFD)을 활용하여 선박 성능 평가를 수행하는 경우, 모형 시험보다 적은 시간과 비용으로 평가가 가능하다는 이점이 있어 이전부터 많은 연구가 진행되어 왔다[1, 2]. 이러한 연구들을 바탕으로 현재 CFD를 이용한 저항 추정은 충분히 신뢰할 만한 수준에 도달했다. 하지만 선행 연구들의 대상선은 주로 낮은 프루드 수의 일반 상선이기 때문에 프루드 수가 높은 고속선에 이러한 기법을 그대로 적용하기에는 어려움이 있다. 일반 상선의 경우, 마찰 저항이 전체 저항의 많은 부분을 차지하는 반면, 고속 활주선은 침수 표면적이 작아 마찰 저항의 비율이 낮고, 대신 압력 저항과 조파 저항의 비율이 상대적으로 크다. 이는 고속선의 선형 특성상 속도가 증가함에 따라 선체가 수면 위로 부양되면서 선저 면에 작용하는 양력이 변화하기 때문이다. 이로 인해 선속에 따라 항주 자세가 크게 변화하며, 이러한 특성을 고려한 성능 추정이 필요하다.

고속선의 성능 평가를 위한 기존의 연구들을 살펴보면, Kim et al.[3]은 모형 시험을 통해 얻은 항주 자세를 적용하여 비선형 포텔션 유동과 점성 유동 해석을 수행하였다. Oh and Yoo[4]는 포텐셜 유동 계산 프로그램을 이용하여 선체에 발생하는 부양력을 계산하고 항주 자세를 추정하여 이를 점성 유동장 해석에 적용하였다. Park et al.[5]은 다양한 선속에서 모형시험과 수치해석을 수행하였다. 저속에서는 모형 시험 결과와 매우 높은 일치도를 보여주는 반면, 고속에서는 저항값의 차이가 커지는 것을 확인하였다. 이는 CFD 해석에서의 트림(trim)과 침하량(sinkage)이 실험 결과와 차이를 보임으로써 오차가 발생한 것으로 분석하였다.

고속선 성능을 보다 엄밀하게 추정하기 위해서는 모형 시험에서의 항주 자세와 유사한 결과를 얻는 것이 매우 중요하다. 이를 위해 항주 자세를 결정짓는 선저면의 압력 분포를 정밀하게 모사해야 한다. 하지만 고속 유동에서는 자유수면 경계면에서 수치 확산(numerical diffusion)이 발생하기 쉽다. 이로 인해 상경계면(phase interface)이 부정확해지면서 침수 표면적을 정확하게 구하기 어렵다. 이러한 문제는 선체에 작용하는 압력 분포에 영향을 미쳐 항주 자세 변화에 민감하게 작용할 수 있으므로 자유수면 경계면에서 발생하는 수치 확산을 최소화 하는 것이 중요하다. 이와 관련하여, 이전부터 날카로운 상경계면을 얻기 위한 수치 기법들이 많이 연구되어 왔다[6. 7. 8. 9. 10. 11]. Ubbink[6]와 Muzaferija et al.[7]은 NVD (Normalized Variable Diagram)에 기반한 수치 기법인 CICSAM(compressive interface capturing scheme for arbitrary meshes)과 HRIC(high resolution interface capturing) 기법을 제안하였다. NVD는 도식의 유계성을 명확히 나타내기 위해 변수들을 정규화하는 방법으로 CBC(Convective Boundedness Criterion)로 나타낸 영역 내에서 수치 기법이 정의된다. CBC 영역은 정규화된 값이 상류 기법과 하류 기법 사이의 값을 갖는 영역을 의미하며, 이 영역 내에서 정의된 수치 기법은 유계성을 만족함으로써 안정성을 높일 수 있다. 또한 상류 및 하류 기법의 장점을 결합하여 새로운 수치 기법을 개발할 수 있다는 장점이 있다. 하지만 이 두 기법의 정확도는 Courant number(Co)에 따라 플럭스를 제한하여 국소적으로 결정되며, 격자 면에서의 물리량을 계산할 때 1차 정확도인 상류 기법 또는 하류 기법이 적용될 수 있다. 이로 인해 CICSAM 기법과 HRIC 기법을 사용할 경우, 격자면의 물리량이 1차 정확도로 계산되어 상경계면에서 수치 확산이 심하게 발생할 수 있다. Rusche[8]는 인위적인 속도를 강제로 부여하여 상경계면 근처에서의 수치 확산을 억제하는 경계면 압축 기법(interface compression method)을 제안하였다. 경계면 압축 기법은 날카로운 형상을 포함한 문제에서 수치 확산에 대한 문제를 비교적 잘 해결한다고 검증된 바 있다[12]. 하지만 부여하는 속도 및 격자 형상에 따라 비물리적인 경계면이 나타날 수 있다고 알려져 있다.

수치 확산이 발생하는 주요 원인은 수치 기법의 정확도, 시간 간격, 격자 간격이다. 속도가 빠르면 경계면에서의 계산이 불안정해지고, 이로 인해 낮은 정확도로 계산된 결과는 수치 확산을 유발하게 된다. 고속 유동은 시간에 따른 변화가 빠르게 일어나기 때문에, 충분히 작은 시간 간격을 적용하지 않으면 수치적 불안정성이 발생하여 유체 간 경계가 부정확하게 계산될 수 있다. 또한 고속 유동에서는 경계면에서 빠른 변화가 일어나기 때문에, 경계면이 여러 격자 셀에 걸쳐 확산되어 정확한 경계면을 유지하기가 어렵다. 따라서 충분히 조밀한 격자를 사용하지 않으면 수치 확산이 발생할 가능성이 크다.

본 연구에서는 상용 CFD 소프트웨어인 STAR-CCM+를 사용하여 정수 상태의 선체 저항 해석을 수행하였으며, 자유수면을 고려한 다상 유동 해석자를 이용하였다. 대상선은 고속 활주선으로 속도 조건은 20, 25, 30, 35, 40 knot이고 프루드 수는 각각 0.64, 0.80, 0.95, 1.11, 1.27이다. 여기서 프루드 수는 배수량을 기준으로 무차원화 하였다. 본 논문의 목적은 수치 기법의 정확도, 시간 간격, 격자 간격이 고속 활주선의 침수 표면적에 미치는 영향을 분석하는 것이다. 따라서 격자계와 시간 간격에 대한 검증을 수행 하였으며, 수치 기법의 정확도가 침수표면적에 미치는 영향을 더 자세히 분석하고자 정확도에 따라 압력 분포, 벽 전단 응력, 물의 체적 분율이 어떤 변화를 보이는지 확인하였다.

2. 대상선 및 수치해석 기법

2.1 대상선

본 연구에서는 단동 활주선에 대한 전산유체해석을 수행하였으며, Fig. 1에 측면부와 전면부의 형상을 나타냈다. 모형선의 길이(L_PP), 폭(B), 흘수(D) 및 모형시험에서 계측한 선속별 자세 변화는 아래 Table 1과 같으며, 실선의 선속은 V_S, 모형선의 선속은 V_M으로 나타냈다. 모형시험에서 계측한 침하량과 트림을 적용하여 자세를 고정한 상태로 수치 해석을 수행하였다.

Fig. 1.

Geometry of planing hull

2.2 수치해석 기법

수치 해석에 사용된 해석 영역의 크기를 Fig. 2에 나타내었다. 좌현 반폭에 대해 유동 해석을 수행하였으므로 중앙면의 경계 조건은 대칭 경계 조건을 적용하였다. 계산 영역에서 아랫면과 측면은 벽면 경계 조건을 적용하였다. 시간과 공간에 대해 2차 정확도의 차분 기법을 사용하였으며, k-ω SST 난류 모델을 사용하였다. 다상 유동 해석을 위한 VOF(volume of fraction)에 대해 HRIC 기법을 사용하였다. 자유 수면을 고려한 수치 해석의 경우 격자면(face)에서 VOF 변수값을 결정하는 기법에 의해 상경계면의 정확도가 결정되는데, HRIC 기법은 상류 기법과 하류 기법의 장점을 혼합한 대표적인 기법이다. 특히 식 (1)로 표현되는 Co와 Fig. 3에 나타낸 자유 수면의 법선 벡터(n→i)와 격자 경계면의 법선 벡터(a→f)가 이루는 각도(θ)를 매개 인자로 갖는 가중치 함수를 통해 격자면에서 VOF 변수값을 계산하는 것이 HRIC 기법의 특징이다. Fig. 3에서 P_U, P_D는 현재의 중앙점(P_C)에 대해 상류점, 하류점을 의미하며, v→는 유동의 속도 벡터를 의미한다. 식 (1)에서 v→f는 격자면에서의 유동 속도, V_PC는 중앙점의 셀 체적을 의미한다.

Fig. 2.

Computaional domain

Fig. 3.

Main parameters in HRIC scheme, where v→, n→i, a→f and θ represent velocity, normal vector to the interface, surface area vector and angle between n→i and a→f, respectively

격자면에서 VOF 변수값을 계산할 때, HRIC 기법이 작동하는 방식을 살펴보면, 먼저 식 (2)와 같이 F(θ)값을 통해 각각 HRIC 기법과 상류 기법에 해당하는 ξfCompressive와 ξfBlend에 가중치를 결정한다. F(θ)값의 계산식은 식 (3)과 같으며, Cθ 는 각도 계수이다. 본 논문에서는 상류 기법에 대한 항을 배제하기 위하여 Cθ는 0을 사용하였다. 그 다음 식 (4)와 같이 정규화된 중앙점의 물리량(normalized cell value, ξC)값에 따라 정규화된 격자면의 물리량(normalized face value, ξf)에 해당하는 값이 결정된다. 식 (4)에 대한 NVD를 도식화 하면 Fig. 4와 같다. 구해진 ξfCompressive는 Co에 따른 수정 계산이 이루어지는데 Co가 작동하는 원리를 자세히 살펴보면, 식 (5)에서 보는 바와 같이 ξC과 ξfCompressive사이에서 내삽의 가중치를 결정하기 위해 Co가 사용되는 것을 확인할 수 있다. 여기서 Co가 하한 한계치(Co_l)보다 작은 경우 격자면의 물리량을 그대로 사용하지만, 상한 한계치(Co_u)보다 커질 경우 ξC를 사용하는 상류 기법이 적용된다. 하한 한계치와 상한 한계치 사이에서는 Co에 따른 선형 내삽이 적용된다.

Fig. 4.

Normalized Variable Diagram(NVD) for the HRIC scheme at different courant numbers

식 (5)를 통해 알 수 있듯이 Co에 대한 하한 한계치와 상한 한계치의 설정을 통해 격자면에서의 VOF 변수값에 대한 계산 정확도를 결정할 수 있다. 본 연구에서는 상한 및 하한 한계치를 동일하게 설정함으로써 선형 내삽 구간의 영향을 배제하였다. 따라서 각 셀의 유동 속도 및 격자 크기, 계산 시간 간격 등의 인자에 의해 HRIC 기법 혹은 1차 상류 기법이 선택적으로 적용된 결과를 얻을 수 있다.

3. 수치해석 결과

3.1 격자계 수렴도 평가

계산 영역에 대해 Fig. 5에서 보는 바와 같이 트리머 격자를 구성하고, 선체 주변에 프리즘 층을 구성하였다. 선체 표면에 사용한 격자수는 62,177개이며, 프리즘 층은 첫 번째 격자 크기를 벽단위 y+ 기준 30으로 설정하고 연신율 1.3에 대해 6개의 층을 생성하였다. 상기 설정을 통해 생성한 총 격자수는 1,588,939개이다. 상기 격자 구성의 타당성을 검증하기 위해 격자 수렴도 인자(grid convergence index, GCI)를 평가하였다. 본 연구에서 사용한 격자를 중간 격자(medium mesh)로 두고, 총 격자수를 2배와 1/2배로 변경한 격자를 조밀 격자(fine mesh), 성긴 격자(coarse mesh)로 칭한다. 이러한 격자를 사용하여 수치해석을 수행하였으며, GCI 평가를 위해 선체에 작용하는 전저항(R_TM)에 대하여 모형시험 계측값(R_TM(EFD))과의 백분위값을 계산하였다.

Fig. 5.

Volumetric grids for planing hull

Table 2에서 보듯이 본 연구에서 사용한 격자 구성은 단조 수렴 구간에 존재하며, 격자 밀집도 증가에 따른 전저항은 모형시험 계측값 기준 약 1.6%의 오차에 수렴한다. 그리고 격자 불확실도는 비보정 기준 약 1.2%, 보정 기준 약 0.4%이다. 본 격자계를 통해 얻어진 물리량의 정성적 타당성을 확인하기 위해 Fig. 6에서 선체 바닥면의 압력 분포와 VOF 변수의 분포를 확인하였다. 그림에서 보듯이 VOF 변수의 수치적 확산 등과 같은 비물리적 현상은 발생하지 않았다.

Fig. 6.

Pressure coefficient and volume fraction of water on the hull surface at the medium mesh

3.2 시간 간격 수렴도 평가

2.2절에서 언급한 바와 같이 상경계면의 차분은 Co수에 의한 영향을 받게 되며, Co수는 식 (1)에 의해 계산 시간 간격(time step, Δt)의 영향을 받는다. 따라서 본 연구에서는 격자의 수렴도 평가와 더불어 계산 시간 간격에 대한 영향을 평가하였다. 계산 시간 간격을 4배씩 증가시키면서 해석을 수행하였으며, 격자 수렴도 인자와 마찬가지로 모형시험의 전저항 계측값(R_TM(EFD))에 대한 전저항(R_TM)의 백분위값을 평가 물리량으로 선정하였다. 평가 방식은 격자 수렴도 인자와 유사하게 차분 기법의 정확도와 리차드슨 외삽(Richardson extrapolation)을 이용한 방식을 적용하였다.

Table 3에서 보듯이 본 연구에서 얻은 전저항은 계산 시간 간격에 대해 단조 수렴 구간에 존재하며, 계산 시간 간격 감소에 따른 전저항은 모형시험 계측값 기준 약 1.8%의 오차에 수렴한다. 그리고 Δt에 따른 불확실도는 비보정 기준 약 2.3%, 보정 기준 약 0.7%이다. 여기서 주목해야 할 점은 2차 정확도의 시간 간격을 적용하였음에도 불구하고 해석 결과의 수렴성을 통해 얻은 정확도(p)가 약 1.09 정도로 도출되었다는 점이다. 수치 해석 과정에서 시간 차분에 따른 오차가 정상적으로 작동했다면 이론상 정확도(pest)에 해당하는 2에 근사할 것으로 기대할 수 있으나, 실제 해석 결과의 정확도가 1 정도에 머물렀다는 것은 시간 차분에 따른 오차 외에 다른 오차의 원인이 작용했음을 의미한다.

3.3 시간 간격 영향

시간 간격 수렴도 평가에서 확인한 바와 같이 계산 시간 간격이 증가함에 따라 선체의 전저항 오차가 급격히 증가하는 특성을 보인다. 이러한 해석 결과의 원인을 분석하기 위해 Table 4에서 보듯 Δt=0.000625 s부터 0.1 s까지 계산 시간 간격을 증가시키면서 선체의 압력, 마찰력 그리고 VOF 변수의 분포를 확인하였다. VOF 변수 분포에서 확인할 수 있듯이 Δt=0.000625 s에서는 선체 바닥면에서 물과 공기의 상경계면이 명확하게 보이며, 수치 확산에 의해 선저로 공기가 유입되는 현상이 발생하지 않는다. Δt=0.00125 s로 증가함에 따라 선체 바닥면의 중앙 부분에서 수치 확산에 의한 공기 유입이 다소 발생하였으나, 선체 바닥면에서의 압력과 마찰력은 Δt=0.000625 s인 경우와 매우 유사한 분포를 보여준다. Δt=0.0025 s부터 선체 바닥면의 중앙 부분에서 수치 확산에 의한 공기 유입이 많아지고, 이로 인해 마찰력 및 압력이 낮게 나타나는 현상이 발생하며, Δt=0.005 s로 커지면 마찰력과 압력 감소가 더 뚜렷하게 관찰된다. 특히 Δt=0.01 s부터는 선저에서 심각한 수치 확산이 발생함으로써 공기와 물의 계면이 불분명해짐에 따라 밀도 및 점성의 감소로 인해 압력과 마찰력이 매우 낮아짐을 확인할 수 있다. 수치 확산에 의한 저항 감소는 Δt가 증가함에 따라 점차 심해지다가 Δt=0.08 s 이상이 되면서 일정하게 유지되는 것을 볼 수 있다. 이상의 결과에서 확인할 수 있듯이 계산 시간 간격이 증가함에 따라 시간 차분에 따른 오차뿐만 아니라 상경계면에서 발생하는 수치 확산이 수치 결과의 해에 큰 영향을 미칠 수 있다.

Table 4.

Pressure coefficient, skin friction coefficient, and volume fraction of water over time steps

Δt	Cp	Cf	Volume fraction of water
0.000625
0.00125
0.0025
0.005
0.01
0.02
0.04
0.08
0.1
Legend

계산 시간 간격에 따른 저항 변화를 Fig. 7(a)에 나타냈다. 그림에서 세로 축은 모형시험의 전저항에 대한 백분위이며, 가로축은 계산 시간 간격이다. 그림에서 보듯이 Δt=0.01 s를 기준으로 선체에 작용하는 전저항이 모형시험 대비 80% 이하로 급격히 감소하고, 이후에는 시간 간격에 따른 변화가 크지 않다. 이러한 급격한 변화의 원인은 표 3에서 언급한 바와 같이 수치 확산에 의한 VOF 변수값의 오차에 기인하므로, 이를 정량화하기 위해 각 해석 결과에 대해 선체의 침수 면적을 계산하여 Fig. 7(b)에 나타냈다. 여기서 A 방법(method A)은 VOF 변수값을 연속 가중치로 갖는 침수 면적이며, B 방법(method B)은 VOF 변수값에 대해 임계치를 적용하여 침수 면적의 반영 여부를 결정하는 방법이다. 정확한 식은 아래 식 (6)과 같다.

여기서 α_i는 선체 격자면에서의 물에 대한 VOF 변수값을 의미하며, S_fi는 선체 격자면의 면적 벡터를 의미한다. 그리고 α₀는 침수 여부를 결정하는 임계치에 해당하는데 보통 0.5를 기준으로 한다. A 방법과 B 방법을 비교할 때 약간의 차이는 있으나 계산 시간 간격이 증가함에 따라 침수 면적이 감소하는 특성은 동일하게 나타난다. 선체에 작용하는 전저항과 침수 면적을 볼 때, 계산 시간 간격으로 사용 가능한 최대값은 Δt=0.005 s임을 확인할 수 있다. 계산 시간 간격의 증가는 식 (1)에 의해 Co의 증가를 야기하므로 식 (5)에서 국소적으로 상한 한계치를 초과하는 셀이 발생할 가능성을 높인다. 즉 계산 시간 간격이 증가함에 따라 VOF 변수에 대해 1차 상류 기법이 적용되는 셀이 증가함으로써 수치 확산에 의한 공기 유입이 발생할 수 있다.

Fig. 7.

Wetted surface area and R_TM/R_TM(EFD) over time steps

3.4 Co의 하한 및 상한 경계치 영향

VOF 변수에 대한 수치 기법이 상경계면의 수치 확산에 미치는 영향을 평가하기 위해 Co의 하한 및 상한 경계치 제어를 통해 VOF 변수의 수치 기법을 결정하였다. 식 (5)에서 Co_l과 Co_u를 동일하게 0으로 설정하면 모든 격자에서 VOF 변수에 대해 1차 상류 기법이 적용되는데, 이러한 설정으로 해석한 결과를 Table 5에 나타냈다. 표에서 보듯이 선속이 증가하거나 계산 시간 간격이 증가함에 따라 수치 확산에 의한 공기 유입이 증가한다. 선속이 20노트인 경우 계산 시간 간격이 증가하더라도 수치 확산이 심각하지는 않지만, 선속이 40노트인 경우 계산 시간 간격이 가장 작은 Δt=0.000625 s에서도 수치 확산이 심각하게 발생하는 것을 확인할 수 있다. 이는 VOF 기반의 이송 방정식에서 대류항을 상류 기법으로 처리할 때 나타나는 수치 확산에 의한 해의 번짐 현상으로, 고속선과 같이 자유 수면 및 선체 주변의 유동 속도가 높은 경우, 대류항의 기여도가 높은 상태에서 상류 기법이 심각한 해의 왜곡을 초래할 수 있음을 보여주는 결과이다.

Table 5.

Volume fraction of water according to the time step changes under various speed conditions when Co_u and Co_l are set to 0

VS Δt	20 knot	30 knot	40 knot
0.000625
0.00125
0.0025
0.005
0.01
0.02
0.04
0.08
0.1

한편 식 (5)에서 Co_l과 Co_u를 동일하게 매우 높은 값(본 연구에서는 1000)으로 설정하면 모든 격자에서 VOF 변수에 대해 고차 수치 기법이 적용되는데, 이러한 설정으로 해석한 결과를 Table 6에 나타냈다. 결과에서 보듯이 선속 및 계산 시간 간격에 상관없이 VOF 변수에 의한 수치 확산이 발생하지 않았다. 결국 식 (5)는 수치 안정성 및 해의 수렴성 등을 위해 1차 상류 기법과 고차 수치 기법의 가중치를 반영하는 의미를 지니고 있으므로 식 (5)에 의해 나타날 수 있는 수치 확산의 가능성을 평가한 결과를 Fig. 8에 나타냈다. 그림에서 세로축은 고속선의 침수 면적을 A 방법와 B 방법으로 추정한 결과를 나타냈으며, 가로축은 계산 시간 간격의 변화를 의미한다. 실선은 고차 차분 정확도의 결과이며, 점선은 1차 상류 기법이 적용되었을 때의 결과이다. 그리고 선속이 20, 25, 30, 35, 40 knot인 경우에 대해 각각 검은색, 붉은색, 푸른색, 녹색, 노란색으로 구분하였다. 각 색으로 채워진 부분은 식 (5)에 의해 1차 상류 기법과 고차 정확도의 수치 기법이 내삽될 때 수치 해가 나타날 수 있는 예상 영역에 해당한다. 영역의 크기가 좁다는 것은 VOF를 고려한 대류항의 수치 기법이 수치 확산에 미치는 영향이 적다는 것을 의미하고, 반대로 영역의 크기가 넓은 것은 수치 확산이 VOF의 수치 기법에 민감하다는 것을 의미한다. 결과에서 보듯이 단동활주선의 유동 해석에 있어 20노트 이하의 속도에서는 계산 시간 간격에 상관없이 민감도가 낮은 반면, 속도가 증가함에 따라 계산 시간 간격보다 VOF의 수치 기법이 더 민감한 것을 알 수 있다.

Table 6.

Volume fraction of water according to the time step changes under various speed conditions when Co_u and Co_l are set to 1000

VS Δt	20 knot	30 knot	40 knot
0.000625
0.00125
0.0025
0.005
0.01
0.02
0.04
0.08
0.1

Fig. 8.

Wetted surface area variations according the the time step changes under various speed conditions(black;20 knot, red;25 kts, blue;30 knot, green;35 knot, yellow; 40 knot)

4. 결 론

본 연구에서는 단동활주선에 대한 유동 해석을 통해 계산 시간 간격 및 VOF 변수에 대한 수치 기법의 정확도가 상경계면의 확산에 미치는 영향을 분석하였다. 격자 수렴도 인자를 통해 단동활주선 주변의 격자 구성의 불확실도를 정량적으로 평가하였으며, 계산 시간 간격에 대한 수렴도 인자 평가를 통해 시간 차분 기법의 정확도뿐만 아니라 계산 시간 간격 증가에 따른 수치 확산이 수치 해의 오차에 영향을 미칠 수 있음을 보였다. 수치 확산에 의한 영향을 보다 체계적으로 분석하기 위해 다양한 선속 및 시간 간격에 대한 해석을 수행하였으며, 상경계면의 수치 확산이 점성과 밀도의 변화를 초래하고 이를 통해 선체에 작용하는 전저항에 영향을 미침을 확인하였다. VOF 변수의 수치 기법 정확도가 1차 상류 기법과 고차 수치 기법의 내삽으로 결정되는 식에서 Co의 하한 및 상한 한계치 설정을 통해 1차 상류 기법 혹은 고차 수치 기법이 적용된 결과를 분석하였다.