지구 온난화 현상으로 북극 해빙의 속도가 가속화되면서 북극 자원 개발이 본격화되고 있다. 최근 러시아는 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)를 사용하여 모델링과 결과 관찰을 수행하므로서 가시화가 뛰어나다. 반면, 그 동안 많은 연구자들이 유체를 모델링하여 유빙의 부력이나 동유체력을 구현하여 왔지만, 매우 방대한 계산 시간으로 인하여 위에서 언급한 여타 방법에 비하여 주목을 받지 못하였다. 본 연구에서는 유체를 모델링하지 않는한 실시간으로 유빙에 작용하는 하중을 생성할 수 없다는 기존의 통념을 벗어나, 유체 모델 없이 유빙에 작용하는 유체력을 실시간으로 생성 및 제어할 수 있는 사용자-서브루틴을 개발하였다. 이 사용자-서브루틴을 양해법 유한 요소 해석에 적용하므로서, 획기적 계산 시간 단축과 정도 높은 빙 저항의 추정이 가능하게 되었다.

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

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

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

Fig. 1 
Subroutine process

요소가 자유 수면보다 높게 위치할 경우 중력만 경험하고, 자유 수면보다 낮게 위치할 경우 중력과 부력을 동시에 경험해야 한다. 요소에 작용하는 중력 가속도는 변동하지 않는다고 가정할 때 식 (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도 이내에 존재하는 경우에만 한다면 양의 항력이 발생하고, 이를 초과할 경우 음의 항력이 발생하므로 각도에 대한 제한 조건이 필요하다.

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

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

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

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