Journal of the Society of Naval Architects of Korea
[ Original Article ]
Journal of the Society of Naval Architects of Korea - Vol. 57, No. 2, pp.88-95
ISSN: 1225-1143 (Print) 2287-7355 (Online)
Print publication date 20 Apr 2020
Received 24 Dec 2019 Revised 4 Feb 2020 Accepted 11 Feb 2020
DOI: https://doi.org/10.3744/SNAK.2020.57.2.088

사용자-서브루틴과 양해법 유한 요소 해석을 이용한 선박의 유빙 저항 추정

한동화1 ; 백광준1 ; 정성엽2 ; 정준모1,
1인하대학교 조선해양공학과
2선박해양플랜트연구소(KRISO)
Ice Floe-induced Ship Resistances using Explicit Finite Element Analyses with a User-subroutine
Donghwa Han1 ; Kwang-Jun Paik1 ; Seong-Yeop Jeong2 ; Joonmo Choung1,
1Dept. of Naval Architecture & Ocean Engineering, Inha University, Incheon, Korea
2Korea Research Institute of Ships & Ocean Engineering (KRISO), Daejeon, Korea

Correspondence to: jmchoung@inha.ac.kr

Abstract

There have been many attempts to predict resistance of vessels in ice floe environment, but they mostly have both strong and weak points at the same time; for instance, simplified formulas are very fast but less flexible to types of ship and ice conditions and other numerical techniques need high computing cost for increased accuracy. A new numerical simulation technique of combining explicit finite element analysis code with a user-subroutine to control real-time forces acting on ice floes was proposed, thereby it was possible to predict ship-to-ice floe resistance with higher convenience and accuracy than other proposed approaches. The basic theory on how real-time hydrostatic and hydrodynamic forces acting on ice floes could be generated using user-subroutine was explained. The heave motion of a single ice floe was simulated using the user-subroutine and the motion amplitudes and periods were almost consistent with analytic values. Towing tests of an icebreaker model ship were simulated using explicit finite element analyses with the user-subroutine. The ice-induced resistance obtained from the towing experiments and simulations showed significant differences. Intentional increase of the drag coefficient to increase the contact duration between the ice floes and rigid model ship leaded the total resistance to be substantially consistent between the model tests and numerical simulations.

Keywords:

Ice floe, Ice resistance, User-subroutine, Drag-force, Buoyancy

키워드:

유빙, 빙 저항, 사용자 서브루틴, 항력, 부력

1. 서 론

지구 온난화 현상으로 북극 해빙의 속도가 가속화되면서 북극 자원 개발이 본격화되고 있다. 최근 러시아는 Yamal 반도에서 생산되는 LNG(Liquefied Natural Gas)를 운반하기 위한 쇄빙 LNG 운반선을 운용하고 있다. 이러한 쇄빙 상선은 북극 항로(Nothern Sea Route, NSR)를 운항하므로 기존 수에즈 운하를 통과하는 항로에 비하여 시간과 비용을 대폭 감소시킬 수 있어, 북극 항로에 대한 관심이 점점 증가하고 있다.

국내 연구 동향을 살펴보면, An et al. (2018)은 아라온호(ARAON)의 북극해 연구 항해를 수행하는 동안 계측된 파랑 하중과 빙 하중 데이터를 바탕으로 선수 프레임 하부의 피로 손상을 평가하는 연구를 수행하였다. Kim (2018b)은 빙 하중 스펙트럼(ice load spectrum)으로부터 시간 영역에서 빙 하중을 고려하는 새로운 접근법을 제안하였다. Cho & Choi (2019)는 쇄빙연구선 아라온호(ARAON)에서 획득한 북극해 현장계측 데이터를 토대로 빙압력-접촉면적의 수정 관계식을 제시하였다.

북극 항로를 운항 가능한 쇄빙선(ice breakers)과 내빙선(ice-strengthened ships)은 해수와의 저항뿐만 아니라 유빙(ice floe)과의 충돌로 인한 부가적인 저항을 경험한다. 따라서 선박의 설계 초기 단계에서 빙 유기 저항을 최소화할 수 있는 선형 설계가 매우 중요하다. 빙 수조(ice basin) 실험을 통하여 가장 신뢰성 있는 빙 유기 저항을 추정할 수 있지만, 비용 및 시간 문제로 인하여 경험식을 이용한 방법, 수치 해석을 이용한 방법 등 다양한 접근법이 개발되고 있다.

Shimansky (1938), Spencer (1992), Poznyak & Ionov (1981)은 경험식을 사용하여 빙 저항을 추정하는 방법을 연구하였다. 위 경험식을 바탕으로 Kim & Lee (2018)은 다양한 빙 입사 조건에 대한 수정된 빙 하중 계산식을 제시하였고 사항 각도(oblique navigation angle)가 빙 하중에 미치는 영향에 대해 연구하였다. 경험식 또는 간이식을 이용한 방법은 빠른 시간에 빙 저항을 예측할 수 있는 장점이 있다. 그러나 이러한 경험식은 선형의 변화, 얼음 조건(두께, 경도, 성분, 염도 등)을 모두 고려하여 개발되지 않았으므로, 범용성에 단점을 가지고 있다.

수치 해석 방법으로는 다물체 동역학(multi-body dynamics, MBD)과 이산화 요소법(discrete element method, DEM) 등이 최근 많이 이용된다. MBD는 유빙과 선박을 강체로 정의하고, 선박의 운항으로 인한 유빙과의 충돌 접촉 마찰 및 충돌 후 유빙의 운동 능력에 기인하는 빙 저항을 추정한다. 여기서 개별 유빙의 운동 능력은 해수 중에서 유빙의 6자유도 거동 특성을 의미하며, 부가수 질량을 무시한다면, 운동 방정식의 감쇠 및 강성 특성에 의하여 결정된다. Kim (2018a)은 단위 얼음의 운동 방정식 상수를 얻기 위하여 CFD(Computational Fluid Dynamics)를 적용한바 있다. Seok et al. (2018)은 CFD를 통한 선박의 6자유도 운동의 해석에서 빙 하중의 영향을 반영하기 위해 개방형 코드(open source code)인 OpenFOAM에 기반한 SNUFOAM을 개발한바 있다. SNUFOAM은 시계열의 빙 하중 데이터를 해석 대상 선형의 무게 중심에 작용하는 외력항으로 적용하여 빙 유기 저항을 추정하였다. Liu et al. (2016)은 평탄빙 및 유빙을 모사할 수 있는 DEM 코드를 개발하였다. 이들은 계산 시간을 감소시키기 위하여 GPU (Graphic Processor Unit) 계산이 가능한 코드를 개발하였다. DEM은 균일한 크기와 형상의 얼음을 가정하는 반면, 불균일한 얼음을 모델링할 수 있는 NDEM(Non-smooth Discrete Element Modeling) 기법도 개발된바 있다 (Read et al., 2018). 이들은 NDEM, 파괴 역학, 내항 해석을 접목하여 선체의 운동을 고려한 평탄빙의 쇄빙 모사를 코드 SAMS(Simulator for Arctic Marine Structures)을 개발한바 있다. 이와 유사한 국내 연구로는 CFD-DEM 연계기법을 이용하여 모노 파일 주위의 입자 거동과 유체 흐름을 해석을 수행한바 있다 (Song et al., 2019). DEM을 위한 LIGGGHTS(LAMMPS Improved for General Granular and Granular Heat Transfer Simulations)과 CFD를 위한 OpenFOAM을 사용하여 수치 해석을 수행하였다.

위 방법들은 비교적 정확한 결과를 얻을 수 있는 장점이 있다. 반면, MBD는 선체와의 접촉을 위한 모사하기 위하여 인하우스 접촉 코드에 의존하였기 때문에 상당히 오랜 계산 시간을 요구되었으며, 새로운 얼음 환경에 대한 빙 저항 추정을 위하여 전처리와 후처리에 상당한 노력이 요구되었다. SAMS는 이산화 모델링에도 상당한 주의가 필요하면서, 여러 기술을 접목한 인하우스 코드(inhouse code)여서 사용자 능력에 매우 민감하고 안정적인 결과를 도출하기 어려운 단점이 있다.

본 연구에서는 상용 유한 요소 코드인 ABAQUS (Simulia, 2018)와 사용자-서브루틴(user-subroutine)를 이용하여 빙 저항을 추정하는 기법을 제시하고자 한다. ABAQUS는 전 세계적으로 매우 넓은 사용자 층을 보유하고 있기 때문에, 사용자가 생성하는 입력에 덜 민감하고, 사용자 입력 오류를 쉽게 검출할 수 있는 강건한 구조 해석 프로그램이다. 유빙의 강체 운동을 허용하기 위하여 양해법(explicit approach) 기반의 유한 요소 해석법을 적용하였다. 이러한 상용 유한 요소 프로그램은 매우 효율적인 접촉 알고리즘을 보유하고 있기 때문에 선체와 유빙과의 접촉으로 인한 계산 시간이 인하우스 코드에 비하여 매우 효율적이다. 또한 전용 전후처리기(pre- and post-processor)를 사용하여 모델링과 결과 관찰을 수행하므로서 가시화가 뛰어나다. 반면, 그 동안 많은 연구자들이 유체를 모델링하여 유빙의 부력이나 동유체력을 구현하여 왔지만, 매우 방대한 계산 시간으로 인하여 위에서 언급한 여타 방법에 비하여 주목을 받지 못하였다. 본 연구에서는 유체를 모델링하지 않는한 실시간으로 유빙에 작용하는 하중을 생성할 수 없다는 기존의 통념을 벗어나, 유체 모델 없이 유빙에 작용하는 유체력을 실시간으로 생성 및 제어할 수 있는 사용자-서브루틴을 개발하였다. 이 사용자-서브루틴을 양해법 유한 요소 해석에 적용하므로서, 획기적 계산 시간 단축과 정도 높은 빙 저항의 추정이 가능하게 되었다.


2. 사용자-서브루틴

ABAQUS 사용 환경에서 사용자-서브루틴은 상용 유한 요소 코드가 지원하지 못하는 특별한 조건이나 환경을 조성할 때 적용된다. 예를 들어 어떤 재료의 점성을 ABAQUS가 지원하지 않는 특별한 변수의 함수로 표현하고 싶다면, 사용자-서브루틴 UVISCO를 이용하여 구현이 가능하다. 일반적인 유한 요소 해석에서 하중은 항상 초기 조건으로 정의하여 작용하므로, 해석 도중에 물체의 변위나 변형에 따른 하중의 제어나 변경은 불가능하다. 이러한 문제도 UDLOAD라는 사용자-서브루틴을 이용하여 실시간으로 사용자가 원하는 하중을 적용할 수 있다.

본 연구에서는 ABAQUS를 이용한 양해법 유한 요소 해석을 통하여 유빙에 작용하는 하중을 실시간으로 제어(생성 및 소멸)하기 위하여 VDLOAD라는 사용자-서브루틴을 사용하였다. 사용자-서브루틴은 포트란(FORTRAN) 컴파일러(compiler)로 코딩되었다.

새로운 해석이 시작되면, ABAQUS는 매 시간 증분 (t)에서의 요소 및 절점 물리량을 계산한다. 사용자-서브루틴은 사용자가 이미 ABAQUS 입력 파일에 정의한 요소에 대한 정보를 ABAQUS 커널(kernal)에 요청한다. VDLOAD 사용자-서브루틴에서 적용한 요소의 요청한 요소 물리량은 요소 중심에서 공간 좌표와 요소의 속도 벡터이다. 접수받은 물리량을 이용하여 사용자-서브루틴은 각 요소별 정수력 (Fb) 및 동수력 (Fd)을 새로이 생성하고 이를 다시 ABAQUS 커널에 되돌려 주므로서, 실시간 하중의 제어가 가능하다 (Fig. 1 참조).

Fig. 1

Subroutine process

요소가 자유 수면보다 높게 위치할 경우 중력만 경험하고, 자유 수면보다 낮게 위치할 경우 중력과 부력을 동시에 경험해야 한다. 요소에 작용하는 중력 가속도는 변동하지 않는다고 가정할 때 식 (1)을 이용하여 실시간 정수력 (부력) 생성이 가능해진다.

Fb(z)=ρsea g z A(1) 

여기서 z은 요소 중심 위치, ρsea는 해수 밀도, g는 중력 가속도, A는 요소의 면적이다.

얼음의 크기가 선박에 비하여 상대적으로 크지 않아서 부가수 질량(added mass) 및 파랑 감쇠(wave damping)로 인한 방사력(radiation force)은 무시될 수 있다고 가정하였다. 또한 유빙 조건에서는 입사파로 인한 파 강제력(wave excitation force) 효과는 없다고 가정하였다. 즉, 형상에 기인한 항력(form-drag force) 만을 동수력으로서 산정하였다. VDLOAD 사용자-서브루틴은 ABAQUS 커널에 요소의 법선 벡터 (n)와 속도 벡터 (v)를 요청한다. 사용자-서브루틴은 이 두 벡터의 내적(dot product)을 통하여 요소에 법선 속도 성분 (vn)을 계산한다 (식 (2)). 이를 식 (3)에 대입하여 동수력 (Fd)을 계산한다. 유빙을 구성하는 요소의 법선 벡터가 얼음의 내부에서 외부를 향한다고 가정할 때, 요소의 법선 벡터와 속도 벡터가 이루는 각이 ±90도 이내에 존재하는 경우에만 한다면 양의 항력이 발생하고, 이를 초과할 경우 음의 항력이 발생하므로 각도에 대한 제한 조건이 필요하다.

vn=nv(2) 
Fd=12Cdρseavn2A(3) 

이러한 알고리즘에 기반한 항력은 유빙이 이동(translation)을 하거나 회전(rotation)을 하는 상황에서도 모두 적용 가능하다. 유빙이 회전하더라도 유빙을 구성하는 특정 요소의 접선 속도만 동유체력을 생성하는데 필요하기 때문이다. 이러한 이유로 요소 항력에 요구되는 항력 계수 (Cd)는 평판에 대한 항력 계수이어야 하며, 얼음의 형상과는 무관하다고 간주될 수 있다.


3. 사용자-서브루틴 교정

3.1 단일 얼음의 운동 주기 및 진폭

사용자-서브루틴의 검증을 위하여 직육면체 형상의 단일 유빙에 대한 수치 해석을 진행하였다. Fig. 2와 같이 요소 분할 정도에 따른 효과를 확인하기 위하여 두 가지 모델(Coarse 및 Fine)이 사용되었다. 단일 얼음의 형상 정보는 Table 1에 제시되었다. 사용된 요소는 감차 적분 4절점 셸 요소(S4R)였으며, 얼음 질량 중심을 참조 절점, 나머지 요소의 절점을 종속 절점으로 하는 변위 커플링 요소를 사용하였다. 참조 및 종속 절점 사이의 변위는 완전 종속으로 가정하였다. 따라서 셸 요소이지만, 강체 요소와 동일한 효과를 유발하였다. 이 참조 절점을 상하 동요를 제외한 나머지 자유도를 모두 구속하여 순수한 상하 동요만을 관찰하였다.

Fig. 2

Single ice models

Information of a single ice

흘수 250mm에 상응하는 흘수를 유지하기 위한 요소 두께 및 질량 밀도를 설정하였다. 수치 해석 모델은 자유 수면으로부터 125mm 잠긴 상태로 모델링 되었다. 이때 적용한 하중은 중력 가속도와 사용자-서브루틴을 이용한 부력 또는/그리고 항력이었다. 부력만 작용한 경우(Case 1) 및 부력과 항력(DNV, 2010)이 동시에 작용한 경우(Case 2)로 나누어 시뮬레이션을 진행하였다. 시뮬레이션 시간은 5.0초 동안 진행되었다.

Case 1은 부력만 작용한 경우이며, 자유 수면(z=0)을 중심으로 감쇠없이 상하 동요가 반복되는 것을 Fig 3.(a)로부터 확인할 수 있다. 즉, 관성력이 부력과 평형을 이루면서, 초기 잠긴 깊이 125mm의 진폭으로 상하 동요가 반복되는 것을 확인할 수 있다. Case 2는 부력과 항력이 동시에 작용한 경우이다(Fig 3. (b) 참조). 항력으로 인하여 자유 수면(z=0)을 중심으로 감쇠하면서 상하 동요가 반복되는 것을 확인할 수 있다. 이러한 감쇠의 효과를 확인하기 위하여 Coarse 모델에 항력 계수를 3.0으로 증가시킨 케이스(Coarse (Cd=3.0))를 Fig 3.(b)에 추가적으로 제시하였다. 항력 계수의 증가는 진폭을 감소시키지만, 상하 동요 주기에 관여하지는 않는 것으로 보인다.

Fig. 3

Heave motion history of single ice models

수치 해석 결과의 타당성을 검증하기 위하여 식 (4)를 이용하여 상하 동요 강성(kh)을 계산하고, 이를 식 (5)에 대입하여 상하 동요 고유 주기(Tn)를 산정하였다.

kh=ρseagLB(4) 
Tn=2πkh/(5) 

여기서, Δ는 얼음의 질량 또는 배수량을 의미한다. 사용자-서브루틴 해석 결과 1.02초로서 이론적 주기 1.0초와 미소한 차이를 보였다. 따라서 사용자-서브루틴은 물리적으로 타당한 고유 주기를 제시하였다.

3.2 단일 얼음의 항력 계수

얼음의 항력 계수를 교정하기 위하여 전후 동요(surge motion)에 대하여 유동 해석을 우선적으로 실시하였다. 유동 해석은 Star-CCM+ (Simens, 2019)을 이용하여 수행되었다. 유체의 지배 방정식으로 식 (6)의 연속 방정식과 식 (7)의 3차원 비압축성 RANS(Reynolds averaged Navier-Stokes) 방정식을 적용하였다. 지배 방정식의 해를 구하기 위하여 각 유한 체적 별로 적분을 취하여 인접하고 있는 유한 체적 내 격자점에서의 물리량들과 상관관계를 구성하여 근사해를 구하는 유한 체적법을 사용하였다. 또한 적용한 난류 모델은 SST k-ω이었다.

uixi=0(6) 
uit+ujuixj=-1ρPxi+ν2uixixj(7) 

유동 해석을 위한 단일 얼음의 밀도와 해수의 밀도를 870 kg/m3 및 1025 kg/m3로 가정하여 얼음의 흘수(0.048m)를 결정하였다 (Table 2). Fig. 4(a)에 보인바와 같은 단일 얼음을 모델링하였으며, 유동 해석에 사용된 격자는 약 200만개로서 이를 Fig. 4(b)에 제시하였다. 단일 얼음 모델에 초기 속도를 부여하므로서, 전진 운동을 모사하였다.

CFD simulation conditions

Fig. 4

Model for computational fluid dynamics simulations

선행 유동 해석을 통하여 얼음 주변의 유동이 너무 심각하게 교란 받지 않는 적절한 수준에서 실제로 가능한 얼음의 초기 속도를 결정하였다. 이렇게 결정된 초기 속도는 최소 0.2m/s부터 최대 1.0m/s까지 0.2m/s 간격으로 총 5수준 속도를 고려하였다.

Fig. 5는 유동 해석으로부터 얻은 속도에 따른 전체 저항 (R)을 나타낸다. 전체 저항은 속도의 증가에 따라 지수 형태로 증가하고 있음을 확인할 수 있다. 이로부터 식 (8)를 사용하여 항력 계수를 산정하여 Fig. 5에 제시하였다. 항력 계수는 운동 속도에 따라 크게 변화하지 않는 것을 확인할 수 있다.

Fig. 5

Total resistances and drag coefficients of an square ice

본 연구에서는 5개 속도 수준에 상응하는 항력 계수의 평균값 (0.9285)을 사용자-서브루틴에 적용하여 전체 저항을 다시 산정하였다. 그 결과, 사용자-서브루틴을 이용한 전체 저항은 유동 해석을 통한 전체 저항과 큰 차이를 보이지 않음을 Fig. 5로부터 확인할 수 있었다.

Cd=R12ρseaAv2(8) 

4. 쇄빙 연구선 모델의 빙 저항 추정

4.1 쇄빙 연구선 모델의 예인 실험 조건

쇄빙 연구선의 예인 실험은 선박해양플랜트연구소 빙해 수조에서 수행되었다. 예인 실험은 실선의 18.667배 축소 모델에 대하여 수행되었다. 예인 실험 조건은 Table 3과 같으며 얼음의 집적도(ice concentration rate) 3수준과 예인 속도(towing speed) 3수준에 따라 총 9회의 빙해 수조 예인 실험이 수행되었다.

Towing test condition

4.2 예인 시뮬레이션 모델링

쇄빙 연구선을 축소한 모델선에 대한 예인 실험을 시뮬레이션하기 위하여 Fig. 6에 보인바와 같이 모델선과 얼음을 유한 요소로 생성하였다. 모델선의 요소 종류는 2차원 강체 요소(R3D4 and R3D3 elements)였다. 실제 예인과 동일한 상황을 재현하기 위하여 질량 중심을 강체 요소의 참조 절점으로 결정하여 이 참조 절점에 예인 방향을 제외한 나머지 방향의 자유도를 모두 구속하였다.

Fig. 6

Finite element model for towing simulation of a research icebreaker

실제 예인 실험에서와 달리 예인 시뮬레이션에서는 정사각형의 얼음 형상을 사용하였다. 따라서 얼음의 간격을 조절하므로서 목표로 하는 얼음 집적도를 달성할 수 있었다. 예인 시뮬레이션을 위하여 적용한 얼음의 크기, 밀도 등은 Table 2를 준용하였다. 얼음에 사용된 요소는 2차원 강체 요소(R3D4 element)였으며, 단일 얼음 질량 중심에 참조 절점을 위치시키고 참조 절점에 질량 요소(MASS element) 및 질량 2차 모멘트 요소(ROTARYI element)를 사용하여 단일 얼음의 질량 정보를 완성하였다.

Simulation conditions

유빙 영역에서의 채널 길이와 폭은 각각 모델선의 2배와 2.5배였다. 얼음의 경계면을 따라 강체 벽을 배치하여 얼음이 채널 내부에서만 이동하고, 채널 외부로 이동할 수는 없었다. 수치 해석은 유빙의 집적도와 선속에 따라 9회 진행되었으며, 항력 계수는 유동 해석의 결과로부터 얻은 0.9285를 사용하였다. 모델선의 참조 절점에 x-방향 변위를 부여하므로서 모델선의 예인을 모사하였다.

4.3 예인 시뮬레이션 결과

모델선의 선미가 유빙 지역에 진입하였을 때 (Fig. 7 (a))부터 선수가 채널 2 m 전방에 도달했을 때 (Fig. 7(b))까지 모델선의 반력 이력을 추출하여 평균값을 빙 저항으로 결정하였다. 이를 유빙 집적도 및 선속에 따라 정리하여 Table 5에 제시하였다.

Fig. 7

Two effective positions for ship resistance calculation

Comparison of towing tests and simulation results

예인 실험과 예인 시뮬레이션으로부터 얻은 빙 저항의 최대 오차율은 집적도가 높을수록 증가하는 경향을 보였으며, 저속일수록 증가하는 경향을 보였다. 따라서 최대 오차율은 얼음 집적도 90% 및 선속 1 노트에서 발생하였다 (Table 5Fig. 8 참조).

Fig. 8

Towing simulation results

예인 실험에서 사용한 얼음은 시뮬레이션에서 가정한 얼음과 달리 비교적 부드럽기 때문에 모델선과의 접촉 시간이 길었을 것으로 추정되며, 이러한 차이가 예인 시뮬레이션의 오차를 크게 발생시켰을 것을 추정하였다. 접촉 시간을 증가시키기 위해서는 얼음을 변형체로 모델링하는 방법과 접촉 하중을 증가시키는 방법이 있을 수 있다. 얼음을 변형체로 모델링할 경우 비교적 실제 상황을 모사할 수 있는 장점이 있는 반면, 계산 시간이 매우 급격히 증가하는 단점이 있다. 양해법 기반 유한 요소 해석의 시간 증분은 변형체 요소 크기에 비례하여 감소하기 때문이다. 반면 마찰 계수를 증가시키거나 항력 계수를 증가시키므로서 접촉 하중의 증가가 가능하다. 마찰 계수의 증가는 수치 해석의 불안정성을 초래할 수 있어서 본 연구에서는 항력 계수를 증가시키는 방법을 선택하였다. 본 연구에서는 시행착오를 거쳐 항력 계수를 2.0으로 증가시킬 경우 예인 실험 결과에 가장 근사하는 시뮬레이션 결과를 얻을 수 있었다 (Fig. 8 참조).

이러한 양해법 기반의 유한 요소 시뮬레이션과 사용자-서브루틴의 조합은 계산 시간의 획기적인 단축을 가져왔다. 선속이 느릴수록 긴 시뮬레이션 시간을 요구하므로 계산 시간이 증가하는 경향을 보였지만, 본 연구에서는 7세대 인텔 i7 10 코어 CPU (central computing unit) 및 32 GB 램 용량 기반의 컴퓨터 환경에서 계산 시간 1시간 이내에 1회의 예인 시뮬레이션을 종료할 수 있었다.


5. 결 론

본 연구에서는 양해법 기반의 상용 유한 요소 해석과 사용자-서브루틴을 이용한 실시간 하중 제어를 통하여 유빙과 모델선의 빙 저항을 추정할 수 있는 가능성을 제시하였다. 사용자-서브루틴의 유효성을 검증하기 위하여 단일 얼음의 상하 동요 진폭과 주기에 대한 시뮬레이션을 실시하였다. 그 결과 정수력(부력)과 동수력(항력)을 모두 포함할 경우 운동 진폭과 운동 주기가 이론해와 거의 일치하는 것을 확인하였다. 전진 운동하는 단일 얼음에 대한 유동 해석을 통하여 항력 계수를 교정한 후에 사용자-서브루틴을 이용한 양해법 유한 요소 해석을 실시한 결과 유동 해석에서 얻은 전체 저항에 근사하는 전체 저항을 얻을 수 있었다. 이로서 새로운 해석 기법의 유효성을 입증할 수 있었다.

쇄빙 연구선의 모형선 예인 실험을 재현하는 시뮬레이션을 실시하였다. 집적도가 높을수록 그리고 선속이 낮을수록 실험과 시뮬레이션의 차이가 큰 경향을 확인할 수 있었다. 항력 계수의 증가를 통하여 접촉 하중의 증가를 유발시킨 후 전체 저항을 비교한 결과, 예인 실험과 매우 근사하는 결과를 얻을 수 있었다. 1회의 시뮬레이션은 일반적인 컴퓨터 환경에서 1시간 이내에 종료될 만큼 혁신적으로 계산 시간을 단축시켰다.

향후 상하 동요 뿐만 아니라 횡 동요(roll) 및 종 동요(pitch)에 대하여도 운동 주기 및 진폭에 대한 검증이 요구된다. 전진 동요에 대하여 도출한 항력 계수(2.0)가 다양한 얼음 형상에 대해서도 적용이 가능한지에 대한 검증이 요구된다. 또한 6자유도 운동에 대해서도 도출된 항력 계수가 유효한지 검증이 필요하다. 변형체 요소를 적용한 얼음 시뮬레이션을 통하여 항력 계수의 증가가 접촉 하중의 증가를 유발한다는 가정에 대한 검증도 요구된다. 균일한 유빙 형상 뿐만 아니라 불규칙한 얼음 형상에 대한 시뮬레이션을 실시하여 얼음 형상의 불규칙성이 전체 저항에 미치는 영향을 검증할 필요가 있다. 유빙뿐만 아니라 평탄빙에도 적용 가능한 사용자-서브루틴의 개발이 요구된다.

Acknowledgments

본 연구(논문)는 산업통산자원부 산업소재핵심기술개발사업 “극한환경용 ICE 내충돌, 고인성 해양플랜트 강재 및 적용 기술 개발” 과제와 산업통상자원부와 한국산업기술진흥원의 “한-영 해양플랜트 글로벌 전문인력양성사업”의 지원을 받아 수행되었습니다.

References

  • An, W.S., Lee, T.K., & Hwang, M.R., 2018. Calculation of fatigue life of bow frame of ARAON considering navigating in ice and open waters. Journal of Ocean Engineering and Technology, 32(6), pp.458-465. [https://doi.org/10.26748/KSOE.2018.32.6.458]
  • Cho, S., & Choi, K., 2019. Modification of local ice load prediction formula based on IBRV ARAON's arctic field data. Journal of Ocean Engineering and Technology, 33(2), pp.161-167. [https://doi.org/10.26748/KSOE.2019.013]
  • Det Norske Veritas (DNV), 2010. DNV-RP-C205 Environmental conditions and environmental loads. DNV.
  • Kim, H., & Lee, J.B., 2018 Estimation method for ice load of managed ice in an oblique condition. Journal of Ocean Engineering and Technology, 32(3), pp.184-191. [https://doi.org/10.26748/KSOE.2018.6.32.3.184]
  • Kim, J., Park, J.C., Kim, S.P., Kim, H.S., & Cho, Y.G., 2018a Multibody dynamics simulation on ice resistance of ship in pack ice environment. Proceedings of the ASME 2018 37th International Conference on Ocean, Offshore and Arctic Engineering, Madrid Spain, 12(S1), pp.88-99.
  • Kim, Y.S., Kim, J.H., Kang, K.J., Han, S., & Kim, J., 2018b Ice load generation in time domain based on ice load spectrum for arctic offshore structures. Journal of Ocean Engineering and Technology, 32(6), pp.411-418. [https://doi.org/10.26748/KSOE.2018.32.6.411]
  • Liu, J., Liu X., & Chen Y., ABS, 2016. Numerical simulations of ice loads on fixed and floating offshore structures using the discrete element method. Proceedings of the Offshore Technology Conference 2016, Newfoundland Canada. [https://doi.org/10.4043/27414-MS]
  • Poznyak, I., & Ionov, B., 1981. The division of icebreaking resistance into components. Proceedings of the 6th STAR Symposium, New York US, pp.249-252.
  • Read, L., Sveinung, L., Wenjun, L., Andrei, T., & Marnix, N., 2018. Simulator for arctic marine structures (SAMS). Proceedings of the ASME 2018 37th International Conference on Ocean, Offshore and Arctic Engineering, Madrid Spain.
  • Seok, W.c., Seo, J.H., & Rhee, S.H., 2018. Application of ice-fluid interaction analysis considering external force term for 6dof motion of an icebreaker. Proceedings of Korean Society for Computational Fluids Engineering, pp.98-99.
  • Shimansky, Y., 1938. Conditional standards of ice qualities of a ship, Northern Sea Route Administration Publishing House, vol. 130.
  • Simens, 2019. Star CCM+ user manual. Simens.
  • Simulia. 2018. ABAQUS user manual. Providence (RI): Dassault Systemes Simulia Corp.
  • Song, S., Jeon, W., & Park, S., 2019. Flow and scour analysis around monopole of fixed offshore platform using method that couples computational fluid dynamics and discrete element method. Journal of Ocean Engineering and Technology, 33(3), pp.245-251. [https://doi.org/10.26748/KSOE.2019.034]
  • Spencer, D., 1992. A standard method for the conduct and analysis of ice resistance model tests. Proceeding of the 23rd American Towing Tank Conference, New Orleans US, pp.301-307.

Fig. 1

Fig. 1
Subroutine process

Fig. 2

Fig. 2
Single ice models

Fig. 3

Fig. 3
Heave motion history of single ice models

Fig. 4

Fig. 4
Model for computational fluid dynamics simulations

Fig. 5

Fig. 5
Total resistances and drag coefficients of an square ice

Fig. 6

Fig. 6
Finite element model for towing simulation of a research icebreaker

Fig. 7

Fig. 7
Two effective positions for ship resistance calculation

Fig. 8

Fig. 8
Towing simulation results

Table 1.

Information of a single ice

Items Unit Value
Length, L mm 1000.0
Breadth, B mm 500.0
Depth, D mm 500.0
Draft, T mm 250.0
Thickness mm 20.0
Shell density ton/mm3 2.5625×10-9
Sea water density ton/mm3 1.025×10-9
Displacement, Δ kg 128.125
Drag coefficient - 2.0

Table 2.

CFD simulation conditions

Item Value
Ice dimensions (m) 0.28×0.28×0.057
Density of ice (kg/m3) 870.0
Density of sea (kg/m3) 1025.0
Draft (m) 0.048
Equation of state RANS
Multiphase model Incompressible
Turbulence model SST k-ω
Prescribed condition Uniform surge velocity

Table 3.

Towing test condition

Item Value
Ice channel width (m) 4.010
Ice thickness (m) 0.057
Density of ice (kg/m3) 870.0
Friction coefficient bet. ice and ship 0.050
Friction coefficient bet. ices 0.081
Flexural strength of ice (kPa) 62.0
Ice concentration (%) 60, 80, 90
Towing speed (knots) 1, 3, 5

Table 4.

Simulation conditions

Item Value
Density of ice (kg/m3) 870.0
Density of sea (kg/m3) 1025.0
Ice dimensions (m) 0.28×0.28×0.057
Length of ice channel (2.5L) (m) 14.125
Sampling rate (s) 0.01
Ice concentration (%) 60, 80, 90
Ship speed (knots) 1, 3, 5

Table 5.

Comparison of towing tests and simulation results

Concent. rate (%) Speed (knots) Ice resistance (N) CPU time (min)
Model test Sim. Rate (%)
60 1 16.33 7.15 56.22 144
3 22.49 13.81 38.60 102
5 31.67 33.75 -6.58 25
80 1 36.92 9.68 73.79 311
3 43.26 22.57 47.82 220
5 49.99 47.01 5.95 47
90 1 42.36 4.67 88.99 436
3 51.04 27.55 46.02 320
5 61.06 64.05 -4.90 48