쐐기형 모형선 주위 연속 쇄빙과정에 관한 입자 기반 수치 시뮬레이션

2.1 지배방정식

본 연구에서는 얼음을 연속적이고 균일한 탄성체로 가정하여 선형 운동량 방정식(Linear momentum equation)과 각운동량 방정식(Angular momentum equation)을 지배방정식으로 사용하였다. 먼저, 선형 운동량 방정식은 아래 식과 같다

여기서 $$ \rho $$는 밀도, $$ \overrightarrow{u}$$는 속도, $$ t$$는 시간, $$ I$$는 단위행렬, $$ ϵ$$는 변형률, $$ \overrightarrow{f}$$는 외부 하중, $$ \lambda $$는 라메의 제1계수, $$ \mu $$는 라메의 제2계수, $$ E$$는 영률(Young’s modulus), $$ \nu $$는 푸아송 비(Poisson's ratio)를 각각 의미한다.

한편, 각운동량 방정식은 다음과 같다

여기서 $$ J$$는 관성 모멘트, $$ \overrightarrow{\omega }$$는 각속도, $$ \overrightarrow{r}$$는 위치 벡터, $$ \overrightarrow{F_s}$$는 전단력을 의미한다.

2.2 Moving Particle Semi-implicit (MPS) 법

MPS법(Koshizuka & Oka, 1996)은 원래 비압축성 자유표면 유동 해석을 위해 개발된 무격자 수치 방법 중 하나이며, 지배방정식의 편미분 연산자를 가중치 함수(weight function) 기반으로 이산화된 편미분 연산자로 대체하여 대수 연립 방정식으로 변환한 뒤 시간 전진에 따라 근사적인 해를 구해 나간다. Chikazawa et al.(2001)은 MPS법으로 탄성체 거동을 묘사한 방법을 제안하였으며, 이후 몇 몇 연구자들에 의해 발전해 왔다(Suzuki et al, 2007; Suzuki & Koshizuka, 2008; Khayyer et al., 2018). 본 연구에서는 Hwang et al.(2011), Hwang et al.(2014), Hwang et al.(2016)이 제안한 구조 해석 모델링 기반의 수치 알고리즘을 사용하여 탄성체 거동을 계산하였다.

2.2.1 가중치 함수

본 연구에서 사용한 가중치 함수는 Lee et al.(2011)이 사용한 6승 함수이며, 아래와 같다.

여기서 $$ r$$은 입자 간의 거리, $$ r_e$$는 유효반경을 의미한다.

2.2.2 기울기 모델

초기의 Koshizuka & Oka(1996)가 제안한 MPS법의 기울기 모델은 식(4)와 같다.

여기서, $$ \Phi $$는 물리량, $$ d$$는 해석하려는 문제의 차원(dimension) 수, $$ \overrightarrow{x}$$는 입자의 위치 벡터, $$ W$$는 가중치 함수, $$ n_0$$는 초기 입자수 밀도를 의미한다. 또한, 아래첨자 $$ ij$$는 물리량의 상대적인 차(예를 들어, $$ {\overrightarrow{x}}_{ij}={\overrightarrow{x}}_j-{\overrightarrow{x}}_i$$)를 나타낸다.

식(4)는 초기의 입자배치가 균일할 경우 비교적 정확한 근사해를 얻을 수 있지만, 입자가 이동함에 따라 비균일한 분포가 될 경우에는 많은 수치적 오차를 포함하게 된다. 이러한 문제를 해결하기 위하여 Khayyer & Gotoh(2011)는 수정된 모델을 제안했다. 하지만 이들 식은 비균일한 입자배치 상태에서 상대적으로 정확한 해를 주지만 초기 입자수 밀도의 사용으로 인해 특히 자유표면 부근의 불안정성 유발 및 프로그램 코딩 시 복잡한 연산을 포함하게 된다.

따라서, 본 연구에서는 아래에 서술하는 것과 같이 음적 기반 기울기 모델을 새롭게 제안한다.

일반적으로 위치 $$ x$$에서의 함수 $$ \Phi $$를 입자 $$ i$$의 위치 $$ x_i$$에서 Taylor 급수 전개하여 2차항 이상을 무시하면 다음과 같다.

완벽한 선형 기울기 연산을 위해, 스칼라 양인 $$ \Phi $$의 1차 항을 소거할 목적으로 기울기가 0인 함수 $$ \Phi \text{'}$$를 다음과 같이 정의한다.

식(5a)를 식(5b)에 대입한 뒤 $$ \Phi \text{'}$$에 대한 기울기를 계산하면 다음과 같다.

이 때, 식(5b)와 (5c)로부터 위치 $$ i$$와 $$ j$$에서 각각 $$ \Phi \text{'}$$의 값과 기울기를 구하면 다음과 같다.

위 식에서 $$ \Phi \text{'}$$의 기울기인 $$ \nabla \Phi {\text{'}}_i$$는 식(4)에 대입하여 아래와 같이 근사할 수 있다.

따라서 식(5d), 식(5e), 식(5f)를 식(5g)에 대입하면

위 식에서 1차 미분 성분인 은 식(5f)에 의해 0이 되어 식(4)와 등가가 되며, 다음과 같이 다시 쓸 수 있다.

위 식을 간략화하면 다음과 같다

여기서 $$ \otimes $$는 텐서곱을 의미한다.

최종적으로, 본 연구에서 새롭게 제안한 기울기 모델은 아래 식(7)과 같이 정리되며, 식(4)에 나타나는 초기 입자수 밀도나 차원 수가 필요 없고 프로그래밍 코딩 또한 간편해 진다.

2.2.3 탄성체 모델

본 연구에서 얼음을 연속적이고(continuous), 동질하며(homogeneous), 등방적인(isotropic) 재질로 가정하여 연속체 이론의 지배 방정식에 포함된 공간 편미분 연산자를 MPS법의 상호작용 모델로 이산화한다.

운동량 보존 방정식을 계산하기 위하여 식(1)의 두 항을 분리하여 따로 계산한다. 첫 번째 항에 포함된 변위벡터의 발산은 MPS의 기울기 모델을 이용하여 아래 식과 같이 이산화한다.

그 다음에 변위벡터의 발산에 대한 기울기를 다음과 같이 이산화하고 첫 번째 항을 구한다.

한편, 식(1)에서 우변의 두 번째 항은 다음 식으로 계산한다.

여기서 $$ R$$은 입자의 회전 행렬을 의미하며, 회전 행렬을 사용하여 입자의 회전 변위가 제거될 수 있다.

다음으로, 각운동량 보존 방정식을 계산하기 위해서 식(2)을 아래 형식으로 이산화하였다.

보다 자세한 내용은 Hwang et al.(2013), Hwang et al.(2014), Hwang et al.(2016)에 구체적으로 서술되어 있다.

2.3 탐색 알고리즘

2.3.1 입자-입자 탐색 알고리즘

얼음 구조물은 입자로 구성되며 임의의 입자에서 유효반경 내에 있는 주변 입자의 탐색은 linked-list 법(Liu et al. 2003)을 기반으로 수행하며, 개략적 탐색(rough searching)과 정확한 탐색(exact searching)의 두 단계에 걸쳐 수행한다.

먼저, 개략적 탐색 단계에서는 Fig. 1에 나타낸 것처럼 전체 계산 영역을 정규 배치된 격자(cell)로 나누고 각 격자마다 포함된 입자 번호를 linked-list에 저장한다.

Fig. 1.

Segmented computational domain into cells and created linked-list to store particle numbers contained in cells

다음으로, 정확한 탐색 단계에서는 Fig. 2와 같이 입자마다 자신 속한 격자와 인접한 격자들의 linked-list에 포함된 입자만으로 탐색을 수행하며, 탐색 조건식은 다음과 같다.

단, $$ i$$는 중심 입자 번호, $$ j$$는 linked-list에 포함된 주변의 입자 번호이다.

Fig. 2.

Cell linked-list searching for neighbor particles

2.3.2 입자-다면체 탐색 알고리즘

입자로 구성된 얼음 구조물과 다면체로 구성된 선체 표면과의 접촉 판단을 위해서 입자-다면체 탐색 알고리즘이 필요하다. 기본적으로, 한 입자가 어느 다면체와 접촉하는지 여부를 판단하기 위해 구조물의 형상 파일로부터 계산 영역 내의 해당 다면체 정보를 탐색해 낸 후 그 다면체까지의 법선 거리를 구해야 한다. 하지만, 이 형상 파일에 다면체 정보가 대량 포함될 경우 임의의 입자와 접촉하는 해당 다면체를 탐색하는데 많은 시간이 소요될 수 있다. 또한, 다면체의 면요소(surface element)는 모서리(edge)와 정점(vertex)을 포함하며, 각각의 모서리와 정점에 관한 정보는 인접한 면요소들과도 공유된다. 이 때문에 특히 법선 벡터의 투영점이 모서리나 정점에 위치할 경우, 인접한 요소와 중복 탐색 및 부정확한 부동소수점자리 문제로 인해 인접 요소가 선택되는 등의 탐색오류가 발생하게 된다. 가장 간편하고 확실한 방법은 면요소와 모서리 및 정점에 대한 법선 벡터를 별도로 탐색하는 것이겠지만, 대부분의 형상 파일 형식(*.STL, *.OBJ 등)이 모서리나 정점의 공유 관계에 대한 직접적인 정보를 포함하지 않기 때문에 탐색이 쉽지 않다. 따라서, 본 연구에서는 입자의 위치 정보와 주변 다면체의 면요소 정보를 이용하여 한 입자로부터 인접한 다면체 상으로의 법선 벡터를 구할 수 있도록 신속하고 정확한 탐색 알고리즘을 적용한다.

먼저, 개략적 탐색 단계에서는 Fig. 3과 같이 전술한 입자-입자 탐색 과정과 유사하게 전체 계산영역을 정규 배치된 격자로 나누고 각 격자마다 자신과 인접한 격자 정보 및 해당되는 다면체 번호를 탐색하고 이를 저장한다. 이 때 하나의 다면체 번호가 자신을 포함하는 모든 격자들에 저장되어야 한다.

Fig. 3.

Cell searching for polygon

본 연구에서는 선체를 강체(ridge body)로 가정하기 때문에 개략적 탐색 단계가 시뮬레이션이 시작하기 전에 한번만 수행하면 된다. 다만, 시뮬레이션 과정에서 선체가 이동함에 따라 격자의 위치와 방향을 수정해야 한다.

다음으로, 정확한 탐색 단계에서는 좌표 변환을 통해 입자가 속한 격자 정보를 얻을 수 있고 개략적 탐색 단계의 결과를 참조하여 임의 입자 근처의 다면체 정보를 파악할 수 있다. 이 후 임의 입자의 투영점이 인근 다면체의 면요소 범위 내에 있는지 판단하고, 판단이 성공되면 해당 다면체로의 수직 벡터를 구한다.

2.4 외부 하중

얼음 입자에 가해지는 외부 하중은 충돌력, 유체력과 중력으로 구성된다. 충돌력은 입자-다면체 탐색 알고리즘으로 구한 수직 벡터를 사용하여 아래 식의 구와 평면에 관한 Hertz 모델(Hertz, 1882)로 계산한다.

여기서 $$ E$$는 탄성계수, $$ L_p$$는 입자 크기, $$ {\overrightarrow{V}}_{nor}$$는 입자로부터 다면체까지의 수직 벡터, $$ L_p/2$$는 입자의 반경, $$ L_p-\left|{\overrightarrow{V}}_{nor}\right|$$는 충돌 지점의 수직 변형량, $$ {\overrightarrow{V}}_{nor}/\left|{\overrightarrow{V}}_{nor}\right|$$는 하중의 방향 벡터를 각각 의미한다.

식(13)으로부터 선체와 충돌한 각 얼음 입자의 충돌력이 구해지고 그 총합의 반발력은 선체에 작용하는 순간 빙 저항이 된다.

유체력으로는 부력만 고려하며, 중력($$ F_g$$)과 부력($$ F_b$$)의 합력은 극한값이 있는 선형 함수로 근사하여 다음과 같이 쓸 수 있다.

여기서, $$ z_1=0.5L_p$$, $$ z_2=-0.5L_p$$, $$ F_b=g\left(1-{\rho }_{water}/{\rho }_{ice}\right)$$, $$ F_g=g$$이다.

정리하면 다음과 같다.

2.5 시간 적분

본 연구에서 시간 적분은 4단계의 Runge Kutta 법을 사용하였다(Press et al., 1992). MPS법을 이용한 유동 해석에서는 기름유출 문제에 대해 4단계의 Runge-Kutta 법이 적용된 바 있다(Jeong et al., 2013).

다음은, $$ m$$번째 시간 스텝에서의 물리량을 이용하여 $$ m+1$$번째 시간 적분의 진행 과정에 대해 서술한다.

여기서 $$ \overrightarrow{{\rm X}}$$는 전체 벡터 변수, $$ {\overrightarrow{f}}_L\left(\overrightarrow{\chi }, \left[R\right]\right)$$는 벡터 형태 물리량의 변화율 함수, $$ {\overrightarrow{f}}_A\left(\overrightarrow{\chi }\right)$$는 행렬 형태 물리량의 변화율 함수, $$ ∆t$$는 시간 간격, $$ f_u$$는 속도벡터의 시간 변화율 함수, $$ f_x$$는 위치벡터의 시간 변화율 함수, $$ f_w$$는 각속도벡터의 시간 변화율 함수, $$ f_R$$는 각변위벡터의 시간 변화율 함수를 각각 나타낸다. 또한, $$ Q$$는 4원수(Quaternion) 혹은 Hamilton 함수로, 다음과 같이 정의된다.

단, $$ \theta $$는 각변위의 크기(norm)이며, $$ \overrightarrow{\Theta }$$는 각변위의 회전축을 각각 의미한다.

마지막으로 $$ m+1$$번째 시간 스텝에서의 속도, 위치, 각속도 및 회전은 다음과 같이 정해진다.

2.6 파괴 모델

본 연구에서 파괴 판정은 최대 법선 응력(maximum normal stress)을 파괴 기준으로 사용하였으며, 판별식은 다음과 같다(Juvinall & Marshek, 2006).

여기서 $$ {\sigma }_f$$는 재료의 파괴 응력, $$ {\sigma }_1$$는 최대 주응력, 즉 응력 텐서의 최대 고유값을 의미한다(Gonzalez & Stuart, 2008).

파단이 발생하는 경우 파단면을 기준으로 양쪽에 위치하는 2개의 입자들간의 연결을 차단할 필요가 있으며, 이 작업을 정확히 수행하기 위하여 파단면 상의 방향 벡터를 구할 필요가 있다. 본 연구에서는 내부 마찰각의 영향을 무시하여 파단면에서 법선벡터를 응력 텐서의 제일 고유벡터로 근사하였다. 제일 주응력의 최대 위치에서 파단이 발생하기 때문에 아래 관계가 성립된다.

여기서 $$ \widehat{n_1}$$는 파단면에서의 법선 벡터를 의미한다.

이에 따라서 파단면 위치는 아래와 같이 스칼라 함수(최대 주응력의 기울기와 방향의 내적)로 정의할 수 있다.

파단면의 한 쪽면의 함수 값이 양이라고 한다면 이와 대응하는 반대쪽면은 음이 되는 성질을 이용하여 두 입자가 파단면 양쪽에 위치하는 조건을 다음과 같이 구성할 수 있다.

최종적으로, 식(19)와 식(22)를 동시에 만족할 경우 두 입자간의 가중치 함수를 0으로 설정하여 상호작용을 차단한다.

2.7 MPI 기반 분산 병렬계산

평탄빙과 모형선의 상호 충돌에 의한 쇄빙과정을 포함하는 시뮬레이션은 일반적인 구조해석에 비하여 해석 과정이 복잡하고 대규모의 계산량을 필요로 한다. 따라서 본 연구에서는 MPI(Message Passing Interface)를 이용한 분산 병렬계산을 수행하여 연산시간을 단축하였다.

분산 계산에는 기본적으로 잘 알려진 영역분할법을 적용하였으며, 각 계산 영역에 포함되는 입자의 물리량 정보를 계산 노드별로 따로 저장하여 연산을 수행하였다. 이 때, 각 노드에서 계산을 수행할 때 인접한 노드에 저장된 입자 정보에 대한 경계조건이 필요하게 된다. Fig. 4에 나타낸 바와 같이 각 영역의 경계 근처의 입자들에서 MPI통신으로 정보를 교환한 뒤 경계조건을 구성하였다.

Fig. 4.

Schematic diagram of boundary condition for local domain on different nodes

추가적으로, 빙파괴 직후 파쇄된 입자가 계산영역간 경계면을 통과하는 경우 이러한 입자들의 정보에 대해 계산 노드간 전송 작업도 필요하게 된다. 하지만, 이로 인해 입자 정보의 정렬 순서가 변할 수 있기 때문에, 노드 내부에서 사용하는 로컬 입자 번호와 전체 계산 영역의 모든 입자에 대한 전체 입자 번호의 두 세트를 사용하여 입자 정보를 정확히 인식하게 하였다.

3.1 얼음 보의 3점 굽힘 테스트

제안된 파괴 모델을 검증하기 위하여 얼음 보의 3점 굽힘 시뮬레이션 결과를 KRISO(Korea Research Institute of Ships and Ocean Engineering)의 빙수조에서 수행된 실험 결과와 비교ㆍ검증하였다. Fig. 5은 초기 설정 및 실험 조건을 표시한다. 아랫면에서의 지지점 사이의 거리($$ L$$)는 0.18m, 시편의 폭($$ B$$)은 0.08m~0.11m, 시편의 두께($$ H$$)는 0.04m~0.06m이다. 윗면의 집중 하중점의 변위는 0.025m/sec 속도로 증가한다. 얼음 시편과 아랫면 지지 구조물과의 접촉은 자유 조건을 부여하였다. 얼음의 밀도는 900kg/m³이고 푸아송 비는 0.3를 사용하였다. 한편, 시뮬레이션의 이산화단계에서 사용된 입자 크기는 2x10^-3m이고 시간 간격은 1x10^-6sec이다.

Fig. 5.

Initial setup of 3-point bending test with ice beam

Fig. 6~Fig. 7에 t=0.0810~0.0822sec까지 3x10^-4sec 간격으로 파괴 과정을 나타낸다. Fig. 6에서는 종방향의 수직응력 분포를 나타내며, Fig. 7에서는 표면 입자(파란색)와 내부 입자(빨간색)를 각각 나타낸다. 단, 얼음 보 내부의 상황을 알 수 있도록 보의 1/4면을 제거하여 표시하였다. 시각 t=0.0813sec 부근에서 파괴가 발생하며 파괴선단에 응력집중 현상과 이로 인해 중앙단면에서 표면입자가 내부입자로 침입되는 것을 확인할 수 있다.

Fig. 6.

Time-sequence of longitudinal normal stress

Fig. 7.

Time-sequence of crack path

얼음 보의 위와 아랫면의 중앙점에 작용하는 수평방향의 수직응력에 관한 시계열을 Fig. 8에 나타내며, 시간 축을 확대하여 Fig. 9에 나타낸다. 처짐의 증가에 따른 응력 값은 거의 선형 관계를 만족한다. 파괴는 아랫면에서 먼저 일어나며, 파괴가 발생하는 순간 응력은 0에 가까워지고 심하게 진동하지만 시간이 지남에 따라 진동폭은 감소한다. 한편, t=0.0814sec 부근의 얼음 보 중앙지점에서의 압축응력은 파괴되기 직전에 다소 증가를 하는 경향을 보이는데, 이는 Fig. 10에서 알 수 있듯이 파괴의 진행 과정 중에 중앙 수직 횡단면에서 압축응력이 인장응력에 비해 커져 일순간 토크 밸런스를 유지 못해 나타난 것으로 볼 수 있다.

Fig. 8.

Time history of longitudinal normal stress at two mid-points on top and bottom surface of ice beam

Fig. 9.

Detailed time history of longitudinal normal stress

Fig. 10.

Longitudinal normal stress at t=0.0814sec

실험에서 계측된 하중, 처짐 및 시편의 기학적 정보를 이용하여 식(23)과 같이 얼음의 탄성계수, 굽힘 파괴응력과 굽힘 파괴변형률을 계산하고 시뮬레이션에는 탄성계수와 파괴응력을 직접 적용한다. 최종적으로 시뮬레이션에서 계산된 파괴변형율을 실험과 비교하고 검증하기로 한다.

단, $$ P_{max}$$는 파괴 발생 시점에서의 최대하중, $$ \delta $$는 파괴가 발생한 시점에서의 처짐, $$ {\sigma }_f$$는 파괴응력, $$ {\epsilon }_f$$는 파괴변형율이다.

Table 1은 4가지 실험 조건과 계측된 데이터, 그리고 Table 2는 식(23)을 통해 구한 얼음 시편의 물성치 값을 정리하여 나타낸다. 실험에서 상대적으로 탄성계수가 큰 Case 01은 담수빙(fresh water ice)이며 그 외는 일반적으로 빙수조에서 사용되는 모형빙(model ice)이 사용되었다.

Table 1.

Raw data of three-point bending test cases

Table 2.

Calculated material properties of ice beam

Fig. 11에는 비교 결과를 표시한다. Case 01의 경우 가장 적은 오차율(8.33%)을 보이고, 그 외는 10.20%~11.52%사이의 오차율을 보인다. 모형빙의 경우 상대적으로 큰 오차가 발생하는 이유는 본 연구에서 얼음을 탄성체로 가정하고 파괴 판정 기준을 단순화된 취성파괴 조건으로 적용했기 때문인 것으로 판단된다. 즉, 담수빙의 파괴 과정에 비해 모형빙의 파괴과정에서 소성변형에 기여하는 성분이 더 높은 비율을 차지하고, 본 연구의 가정에 더 가까운 담수빙의 오차가 더 낮게 나타나기 때문이다.

Fig. 11.

Comparison of flexural strain with experiments

3.2 쐐기형 선수 주위 쇄빙 시뮬레이션

Fig. 12과 같이 20°의 선수각을 갖는 쐐기형 모형선 주위의 연속쇄빙 실험이 KRISO의 빙해 수조에서 수행되었다. 모형선의 길이와 높이는 각각 0.722m와 0.262m이고 반폭은 0.495m이다.

Fig. 12.

Schematic view of wedge used for level-ice test

시뮬레이션에서는 Fig. 13에 나타낸 바와 같이 평탄빙의 폭이 5.6m의 조건을 사용하였으며, 쐐기형 모형선이 일정 속도로 전진하는 동안 작용하는 빙저항을 계산하였다. 시뮬레이션에서는 Table 3과 같이 실험에서 계측된 평탄빙의 물성치를 사용하였다.

Fig. 13.

Initial set-up of simulation for level-ice test

Table 3.

Other material properties of ice

평탄빙에 대한 입자 기반 시뮬레이션에서 평탄빙의 두께가 길이와 폭에 비해 상대적으로 매우 얇기 때문에 두께 방향에 배치 가능한 입자 개수가 많지 않게 되며 이로 인해 시뮬레이션 결과에 큰 영향을 미칠 수 있다. 따라서 본 연구에서는 세 가지 입자 크기(두께 h에 각각 2, 3, 4개의 입자를 배치)에 대한 수렴 테스트를 수행하고 타당한 입자 크기를 선정한 다음, 쐐기형 모형선의 세 종류의 선속 변화에 따른 빙저항 추정 시뮬레이션을 수행하였다. 시뮬레이션에 사용된 입자 크기 및 쐐기형 모형선의 속도 조건을 Table 4에 정리하여 나타낸다.

Table 4.

Test matrix of the wedge level-ice simulation

Fig. 14는 입자 크기별 빙저항의 시계열 결과를 나타낸다. 전체적으로 모형선이 얼음과의 충돌로 인해 받는 충격력이 매우 불규칙하게 나타나며, 입자 크기가 가장 큰 Case L1의 빙저항이 다른 경우에 비해 다소 낮게 추정되는 것을 알 수 있다. 정량적 비교를 위해 빙저항의 시간 평균과 오차율을 산정하여 Fig. 15에 나타낸다. 결과적으로 입자의 크기가 감소함에 따라 평균 빙저항 값이 수렴 경향을 보이는 것을 확인할 수 있다. 계산의 정확도과 연산의 효율성을 종합적으로 고려하여 이후의 속도별 시뮬레이션에서는 Case L2(V1)의 입자 크기를 사용하기로 한다.

Fig. 14.

Time-history of Ice resistance on wedge-type model ship for convergence test of particle size

Fig. 15.

Time-averaged ice resistance and error ratio for convergence test of particle size

Fig. 16은 Case L2(V1)의 t=3.8sec에서의 주응력 분포(왼쪽)와 쇄빙 패턴(오른쪽)을 나타낸다. 쐐기형 모형선이 지나간 후류 경로에 얼음의 파쇄 현상이 확인된다. 또한 모형선 부근에서 평탄빙의 외부로 전파되는 환형의(ring-shaped) 응력파(stress wave)가 목격되는데, 이는 얼음을 탄성체로 가정한데서 비롯된 결과라고 할 수 있다. 한편, 1차 주응력 분포와 파괴 패턴을 비교하면 가장 높은 1차 주응력 값이 얼음 조각(ice floe)의 충돌된 부위와 crack의 끝단에서 나타나는 것을 알 수 있다.

Fig. 16.

Snapshots of first principal stress(left) and fracture pattern(right) at t=3.8 sec in case of V1(L2)

다음으로, 쐐기형 모형선의 선속에 따른 세 가지 Case(V1, V2, V3)의 빙저항에 관한 시계열을 Fig. 17에 나타낸다. 시뮬레이션의 초반에 쐐기형 모형선의 빙저항은 가속 구간에서 점차적으로 증가한 뒤, 이후 일정 범위 내에서 불규칙적으로 변동하는 것을 알 수 있다. 또한 모형선의 선속이 클수록 빙저항이 크게 나타나는 것을 알 수 있다.

Fig. 17.

Time-history of Ice resistance on wedge-type model ship with different speeds

본 연구에서는 앞서 서술한 바와 같이 선체 및 얼음과 상호작용하는 유체의 운동은 부력만을 고려하기 때문에 빙저항 값을 실험과 비교하기 위해서는 open water에서의 저항 값이 추가적으로 필요하다. 따라서 상용 CFD(Computational Fluid Dynamics) 소프트웨어인 Star-CCM+을 이용하여 쐐기형 모형선의 open water에서의 저항 값을 추정하였다. Fig. 18는 Case V1(L2)의 시뮬레이션 결과이며, 자유표면의 연직 변위 및 상부의 순간 유선을 나타낸다. 자유표면은 선수 부분에서 상승하고 선미 후방에서 하강하는 것을 확인 할 수 있다. 이 때 모형선에 작용하는 선속별 유체저항의 결과를 Fig. 19에 나타낸다. 유체저항은 선속의 제곱에 정확히 비례하여 증가하고 있으며, 대부분 압력저항이 지배적이고 마찰저항은 전체 유체저항의 5~8%정도를 차지한다.

Fig. 18.

Free-surface elevation and stream line fields around wedge-type model ship for Case L1(V2)

Fig. 19.

Open water resistance with different wedge velocity

마지막으로, open water 저항을 포함한 총 빙저항을 KRISO의 빙해수조에서 수행된 연속 쇄빙의 실험결과와 비교하하여 Fig. 20에 나타낸다. 선속이 증가함에 따라 시뮬레이션 및 실험 모두 빙저항이 증가하는 경향을 보이고 있다. 시뮬레이션의 빙저항은 실험에 비해 약 7~40% 가량 크게 나타나며 선속의 증가에 따라 오차율은 감소하는 것을 알 수 있다. 한편, 총 저항에서 open water 저항이 차지하는 비율은 선속에 비례하여 나타난다.

Fig. 20.

Total ice resistance on wedge-type model ship with different velocity

3.3 분산 병렬계산 성능

앞선 모든 시뮬레이션은 HPC(High Performance Computing) 시스템 상에서 수행되었으며, 사용된 HPC의 사양은 Table 5에 정리하였다. Fig. 21은 MPI 기반 분산 시뮬레이션의 노드별 가속비를 나타낸다. 최대 9노드를 사용할 경우 약 80% 미만의 성능을 보이지만, 총 입자수가 증가함에 따라 전체적으로 계산 효율이 개선되는 경향을 나타내고 있다.

Table 5.

Specification of HPC used in simulation

Fig. 21.

Speedup ratio for MPI-based parallel computing