지구온난화로 인해 극지 자원개발에 대한 기대가 커지면서, 극지에서의 석유자원 시추를 위한 드릴쉽(Drillship)이나 반잠수식 시추선(Semi-submersible rig)에 대한 관심은 나날이 증가하고 있다. 특히 반잠수식 시추선은 선박 형태의 드릴쉽에 비해 수선면적이 적어 파랑 및 빙하중을 적게 받기 때문에 극지의 시추작업을 위한 해양구조물로 각광을 받고 있다.

극지용 반잠수식 시추선은 그 환경의 특수성으로 인해 일반 해역에 설치되는 것과 달리 추가로 고려해야할 사항들이 있는데 그 중 가장 중요한 요소가 바로 빙하중이다. 특히 유빙(Pack ice)에 의한 빙하중은 그 발생 빈도가 높고 선체와 계류라인에 영향을 줄 수 있어 설계단계에서 정확한 예측이 필요하다.

빙하중을 예측하기 위한 수치해석법의 개발은 컴퓨터의 계산 성능이 발달하면서 활발히 진행되고 있다. Metrikin et al.(2015)은 최근의 빙하중 예측용 수치해석법을 유한요소법(Finite element method), Particle-in-cell(PIC) 그리고 개별요소법(Discrete element method)으로 크게 세 가지로 나누었다. Kim et al.(2013), Wang and Derradji-Aouat(2011) 그리고 Millan and Wang(2011)은 유한요소법의 ALE(Arbitrary lagrangian eulerian) 기법을 적용하여 얼음과 물과 구조물의 상호작용을 구현하였다. Particle-in-cell 기법은 Sayed(1997)에 의해 개발된 방법으로 Semi-lagrangian이라고 할 수 있다. 즉, 이 기법에서는 얼음을 개별입자로 묘사하고, 이들 입자들이 Eulerian 도메인으로 구현된 바다 위에서 움직이도록 하였다. 개별요소법은 Cundall and Strack(1978)에 의해 개발되어 수많은 입자들의 움직임과 충돌을 구현하기 위한 수치해석법으로 널리 알려져 있다. Sun and Shen(2012)은 개별요소법을 이용하여 Pancake ice에 대한 빙하중을 계산하였다. 이 밖의 방법으로, Alawneh(2014)은 GPU(Graphics processing unit)를 이용하여 빙하중을 계산하는 GPU-based event mechanics을 개발하였다. 이 방법은 2D로 구현되는 단점이 있지만 실시간 시뮬레이션이 가능할 정도로 매우 빠르게 계산할 수 있는 장점이 있다.

본 연구에서는 유한요소법을 이용하여 유빙(Pack ice) 조건에서 계류된 반잠수식 시추선에 대한 빙하중과 거동을 계산하였다. 해석시간의 절감을 위해, 얼음과 해수, 얼음과 구조물 그리고 얼음들 간의 상호작용을 간략화된 방법(Simplified method)을 통해 구현하였고, 정확성의 향상을 위해 간략화 방법에 사용되는 주요 계수(Coefficient)들은 별도의 상세 해석을 통해 계산된 후 적용되었다. 최종적으로 계산된 빙하중은 러시아의 Krylov 연구소에서 수행된 해석의 결과와 비교하여, 제안된 방법의 적용성을 평가하였다.

계산에 사용된 프로그램은 상용유한요소해석 프로그램인 ABAQUS/Explicit이며, 자동 모델링 기법을 개발하여 임의의 형상을 가진 유빙을 모델링하였다. 유빙의 거동은 User subroutine을 이용하여 모사하였으며, 유빙의 충돌로 인해 발생하는 얼음의 국부적 손상(Local failure)은 압력-침투량(Pressure-penetration) 관계를 통해 구현하였다.

얼음은 대부분 해수에 부분적 또는 전체적으로 잠겨있기 때문에 이로 인해 발생하는 Drag force, 부가질량(Added mass), 관성력(Inertia) 그리고 부력(Buoyancy force)과 같은 유체역학적 하중을 받게 된다. 본 연구에서는 이러한 하중으로 인한 얼음의 거동을 구현하기 위해 Fig. 1과 Fig. 2 그리고 식 (1)-(3)에서 나타낸 바와 같이 6자유도를 고려한 식을 적용하였다. 먼저, 식 (1)과 식 (2)의 Drag force(Fd→)와 Drag moment(Md→)는 Surge, Sway, Yaw와 같은 평면운동(Planar motion)에 적용되었다 (Hopkins and Shen, 2001). 반면, Heave, Roll, Pitch는 식 (3)에서 나타낸 바와 같이 얼음의 표면에 부력(P)의 형태로 적용하였다. 또한, Roll과 Pitch와 같이 얼음이 기울어진 상황을 고려하기위해 해당 면(Surface)에 국부좌표계 x^'와 y^'를 따라 압력을 계산하여 적용하였고, 대변형을 고려하기 위해 얼음의 모든 면에 압력이 적용되도록 하였다. 부가질량의 구현을 위해 기존 질량에 1+c_m을 곱한 값을 질량으로 대체하였으며, 여기서 c_m은 Newman(1977)에 따라 0.41을 적용하였다.

여기서 CdF는 Drag force 계수를 나타내며, ρ_ω는 해수의 밀도, A_s는 투영면적을 의미한다. V→는 얼음의 속도, 그리고 Vω→는 조류의 속도를 나타낸다. r은 얼음의 중심점과 꼭지점들 사이의 거리를 평균한 값이며, ω→는 각속도를 나타낸다. g와 h는 각각 중력가속도와 해수면으로 부터의 거리를 의미한다. Drag force 계수(CdF)는 Coupled Eulerian-Lagrangian 방법을 이용한 상세 해석을 통해 얻었다. Fig. 2은 Drag force 계수를 얻기 위한 상세해석(a)과 본 연구에서 제안한 방법을 이용한 해석(b)의 모습을 나타내고 있다. 상세해석에서는 얼음을 Lagrangian 도메인, 그리고 해수를 Eulerian 도메인으로 구현하여 얼음을 해수에 뜨도록 한 후 강체로 구현된 선박모델을 전진시켜 충돌하게 하였다. 선박모델로는 대한민국 최초의 연구용 쇄빙선인 아라온호를 적용하였다. 제안된 방법을 이용한 해석모델에서는 얼음과 선박모델은 동일하게 생성한 후 해수를 모델링하는 대신 식 (1)-(3)에서 설명한 Drag force와 부력, 부가질량을 적용하여 얼음과 해수와의 상호작용을 구현하였다. Drag force 계수는 두 해석모델에 대해 일련의 파라메트릭 연구를 수행하여 유사한 얼음들의 거동을 만드는 값으로 정하였다. 최종적으로 얻어진 Drag force 계수는 10.5이며, 선박과 가장 가까운 곳에 위치한 4개의 얼음의 거동을 비교한 결과는 Fig. 3에 나타내었다. 결과적으로 회전운동을 주로 하는 얼음들(1번과 3번)에 비해 직선운동을 주로 하는 얼음들(2번과 4번)에 대한 정확도가 더 높은 경향을 보였다.

일반적으로 얼음과 구조물 또는 얼음들 간의 충돌 시에 발생하는 국부적 손상(Local failure)은 일반적으로 얼음의 다양한 물성치를 고려한 상세한 해석을 통해 구현된다. 하지만, 이러한 상세 해석은 막대한 해석 시간을 요구하기 때문에 수많은 얼음의 충돌을 구현하는 유빙 상태에서의 해석 시에는 적합하지 못하다. 대신, 본 연구에서는 얼음의 충돌 시 발생하는 국부적 손상을 충돌하는 요소들 간의 침투(Penetration)를 통해 구현하도록 하였다. Fig. 4에서는 본 연구에서 적용한 압력-침투 관계(Pressurepenetration relation)를 설명하고 있다. 즉, 얼음의 충돌 시 발생하는 압력에 따라 침투량을 다르게하여 국부적 손상에 의한 하중과 얼음 거동의 변화를 구현하였다. 하중이 아닌 압력과 침투 간의 관계를 사용하는 것은 다양한 충돌의 상황을 고려하기 위해서이다. 즉, 압력은 하중과 접촉면으로 계산되기 때문에 충돌하는 면의 형상에 상관없이 일관된 상관관계를 얻을 수 있다. 얼음의 압력-침투 관계는 여러 가지 충돌 시나리오를 고려한 탄소성해석을 통해 구하였다. 본 연구에서는 Table 1에서 나타낸 10가지 상황에 대해 탄소성해석을 수행하였다. 1년생 빙(First-year ice)의 경우 최대 두께가 2m이지만 1m이하가 대부분을 차지하므로 두께는 0.3m, 0.6m, 0.9m를 선택하였다. 얼음의 크기는 각 얼음두께에 대해 Lindseth(2013)의 기준에 따라 더이상 깨지지 않는 얼음의 최대 크기와 그 최대 크기의 반으로 선택하였다. 마지막으로, 각 경우에 대해 정면충돌과 45° 충돌 경우를 추가하였다.

최종적으로, 해석을 통해 얻은 각 Time step별 압력과 침투량과의 관계를 Fig. 5와 같이 나타낸 후 Kärnä(1993)에 따라 선형회귀(Linear regression)하여 압력-침투 관계를 얻었다. 도출된 압력-침투 곡선의 기울기는 35,367Pa/m이다.

본 연구에 적용된 대상 구조물은 Fig. 6과 같이 4개의 칼럼(Column)을 가진 반잠수식 시추선이다. 주요제원은 Table 2에 나타내었으며, 좌표계는 Fig. 7과 같이 얼음의 진행방향을 x축으로 놓았다.

Fig. 8에서는 대상 시추선에 적용된 계류라인을 나타내고 있다. 각 칼럼 당 2개씩 총 8개의 계류라인이 45° 간격으로 적용되었으며, 라인의 총 길이는 1,920m이고, 수중단위중량(Submerged weight per unit length)은 1,511kg/m이다. 계류라인의 강성(Stiffness)은 모든 방향에 동일하게 Fig. 9와 같이 적용하였다.

본 논문에서 제안한 방법을 이용하여 유빙 조건에서의 수치해석을 수행하였다. 본 연구에 적용된 해석 조건은 Table 3과 같다. 다른 조건은 모두 동일하며 얼음의 속도만 다른 두 개의 경우에 대한 해석을 수행하였다.

먼저, 반잠수식 시추선과 계류라인을 Fig. 10과 같이 단순화하여 모델링하였다. 4개의 칼럼을 강체(Rigid body)로 모델링한 후, 이 칼럼들이 시추선의 중심(Center of gravity, COG)을 따라 움직이도록 연결(Coupling) 하였다. 시추선의 질량과 관성모멘트는 COG에 적용하였으며, 계류라인은 각 칼럼에 연결된 X, Y방향 스프링을 통해 구현하였다. 얼음은 2장에서 소개한 방법을 이용하여 모델링하였다. 임의의 형상을 가진 사각형으로 구현하였으며, 대상 시추선 주변에 80%의 밀집도(Concentration)로 배치하였다. 모델의 가장자리 얼음들의 경계조건을 구현하기 위해 강체의 벽(Rigid wall)을 생성하였으며, 경계조건 효과(Boundary effect)를 최소로 하기 위해 유빙과 구조물의 의 크기를 고려하여 채널의 폭은 500m로 정하였다. Fig. 11에서는 전체모델의 모습을 나타내고 있다.

1,400초 동안 해석을 수행하였으며, 해당시간 동안 대상 시추선이 받는 빙하중과 움직임의 변화를 측정하였다. Fig. 12는 해석의 절반이 진행되었을 때 얼음의 분포와 시추선의 움직임을 보여주고 있다.

해석을 통해 도출된 결과는 러시아 크릴로프 연구소(Krylov State Research Centre, KSRC)에서 최근 수행한 결과와 비교하였다. KSRC에서는 동일한 조건에 대해 개별요소법(Discrete element method, DEM)을 이용하여 계산을 수행하였으며, 각 얼음은 2차원의 원형으로 모델링하였다. 또한, 현실적인 조류(Current)의 적용을 위해 정지상태부터 서서히 속도를 올려 0.5m/s 조건의 경우는 900초, 0.8m/s 조건의 경우는 500초에 목표속도에 도달하도록 하였다(Karulin and Karulina, 2013). Fig. 13에서는 Case 1의 조건에 대해 각각의 방법으로 계산된 대상 시추선의 거동과 하중을 비교하였다. 본 논문에서 소개된 방법을 통해 계산된 결과는 ‘FEM’, 그리고 KSRC에서 계산한 결과는 ‘KSRC'로 각각 표시하였다. 먼저, X방향(빙하중을 받는 방향)을 살펴보면, 두 경우 모두 목표속도에 도달하게 되는 900초가 지난 이후부터 일정한 주기를 가지고 진동하는 것을 볼 수 있다. 대상 시추선의 질량과 계류라인의 강성을 고려해 계산된 고유주기는 약 170초이며, 두 결과 모두 고유주기와 유사한 주기로 진동하는 것을 확인할 수 있다. 해석의 초기에 900초까지의 결과가 다소 차이를 보이는데, 이것은 두 해석의 얼음의 초기 분포와 형상이 다르기 때문인 것으로 판단된다. 하지만 일단 목표 속도에 도달하면 두 경우가 서로 유사한 진폭으로 거동하는 것을 볼 수 있다. 반면, Y방향(빙하중을 받는 방향과 직각인 방향)으로는 움직임이 매우 적은 것을 알 수 있다. 다음으로 X방향의 하중을 살펴보면, 해석시간 동안 발생한 최대 빙하중이 유사한 것을 확인할 수 있다. 다만, ‘KSRC’의 경우가 하중의 변동성이 더 큰 모습을 보였다. 마찬가지로 Y방향의 하중은 X방향에 비해 무시할 정도로 적은 값을 보였다. 이는 대상 구조물의 칼럼(Coulmn)의 경사각이 0°이기 때문으로 판단된다. 경사각이 큰 구조물의 경우는 Y방향의 움직임도 증가할 수 있을 것으로 예상된다. Fig. 14에서는 Case 2의 조건에 대해 각각의 방법으로 계산된 대상 시추선의 거동과 하중을 비교하였다. 먼저, X방향을 살펴보면, 두 경우가 모두 목표속도에 도달하게 되는 500초가 되는 시점에서 최대 변위를 보이고, 그 값 역시 유사한 수준을 보였다. 진동하는 주기는 유사한 경향을 보이나 그 값은 Case 1에 비해 다소 큰 차이를 보였다. Y방향 거동의 경우 ‘KSRC’는 500초가량에서 최대 변위를 보인 이후에는 적은 거동을 보인데 반해 ‘FEM’의 경우는 계속해서 일정한 주기를 가지고 비교적 큰 거동으로 움직이는 것을 알 수 있다. X방향의 하중은 Case 1과 유사하게 최대 빙하중이 유사한 반면 변동성에는 비교적 큰 차이를 보였다. 역시 Y방향의 하중은 X방향에 비해 무시할 정도로 적은 값을 보였다.

본 연구에서는 유한요소법을 이용하여 유빙 상태에서의 빙하중을 계산하는 새로운 방법을 제안하고, 이를 이용하여 계류된 반잠수식 시추선의 빙하중 계산 및 모션 해석에 적용하였다. 상기의 연구를 토대로 다음과 같은 결론을 도출하였다.

(1) 자동 모델링 기법을 적용하여 임의의 형상 및 크기를 가지는 얼음을 모델링하였다. 이로 인해, 유빙상태의 모델링에 필요한 시간을 대폭 절감하였다.

(2) 얼음과 해수, 얼음과 구조물, 얼음들 간의 상호작용은 간략화 기법(Simplified method)을 적용하여 모델링하였다. 이 때 필요한 여러 중요한 계수들은 상세한 해석을 통해 도출하였다.

(3) 얼음의 충돌로 인해 발생하는 국부손상(Local failure)은 압력-침투 관계(Pressure-penetration relation)를 통해 구현하였다. 얼음과 구조물의 압력-침투 관계는 발생 가능한 상황에 대해 상세해석을 수행한 후 계산된 결과들을 선형회귀하여 도출하였다. 상세해석에는 얼음의 탄소성이 고려된 물성치를 사용하였으며, Coupled Lagrangian-Eulerian 기법을 적용하였다.

(4) 얼음과 해수의 상호작용 모델링에 중요한 역할을 하는 Drag force coefficient 또한 Coupled Lagrangian-Eulerian 기법을 이용한 상세해석을 통해 얻었다. 상세해석에서는 얼음을 Lagrangian 도메인, 그리고 해수를 Eulerian 도메인으로 구현하여 얼음을 해수에 뜨도록 한 후 강체로 구현된 선박모델을 전진시켜 충돌하게 하였으며, 본 연구에서 제안된 방법을 이용한 결과와의 비교를 통해 Drag force coefficient를 최종 도출하였다.

(5) 제안된 방법을 유빙 조건에 계류된 반잠수식 시추선에 적용하였다. 해석을 통해 도출된 결과는 러시아 크릴로프 연구소에서 최근 수행한 결과와 비교하였다.

(6) 동일한 조건에 대해 해석을 수행하여 대상 시추선의 거동과 빙하중을 비교한 결과, 두 경우 모두 고유주기와 유사한 주기로 진동하는 것을 확인할 수 있었으며 움직임의 최대값은 유사하나 변동폭은 다소 차이가 있었다. 빙하중도 마찬가지로 해석 기간 동안 계측된 최대값은 두 경우가 유사한 결과를 보였으나 변동폭에서는 차이를 보였다. 전체적으로 0.8m/s보다는 0.5m/s 속도 조건에서 더 유사한 결과를 나타내었다.