Current Issue

Journal of the Society of Naval Architects of Korea - Vol. 61 , No. 1

[ Original Article ]
Journal of the Society of Naval Architects of Korea - Vol. 57, No. 2, pp. 80-87
Abbreviation: J. Soc. Nav. Archit. Korea
ISSN: 1225-1143 (Print) 2287-7355 (Online)
Print publication date 20 Apr 2020
Received 30 Oct 2019 Revised 20 Jan 2020 Accepted 4 Feb 2020
DOI: https://doi.org/10.3744/SNAK.2020.57.2.080

물성치의 공간분포를 고려한 빙 시험편의 확률론적 강도평가
김호준1 ; 김유일1,
1인하대학교 공과대학 조선해양공학과

Probabilistic Strength Assessment of Ice Specimen considering Spatial Variation of Material Properties
Hojoon Kim1 ; Yooil Kim1,
1Department of Naval Architecture and Ocean Engineering, College of Engineering, INHA University, Korea
Correspondence to : yooilkim@inha.ac.kr

Funding Information ▼

Abstract

As the Arctic sea ice decreases due to various reasons such as global warming, the demand for ships and offshore structures operating in the Arctic region is steadily increasing. In the case of sea ice, the anisotropy is caused by the uncertainty inside the material. For most of the research, nevertheless, estimating the ice load has been treated deterministically. With regard to this, in this paper, a four-point bending strength analysis of an ice specimen was attempted using a stochastic finite element method. First, spatial distribution of the material properties used in the yield criterion was assumed to be a multivariate Gaussian random field. After that, a direct method, which is a sort of stochastic finite element method, and a sensitivity method using the sensitivity of response for random variables were proposed for calculating the probabilistic distribution of ice specimen strength. A parametric study was conducted with different mean vectors and correlation lengths for each material property used in the above procedure. The calculation time was about ten seconds for the direct method and about three minutes for the sensitivity methods. As the cohesion and correlation length increased, the mean value of the critical load and the standard deviation increased. On the contrary, they decreased as the friction angle increased. Also, in all cases, the direct and sensitivity methods yielded very similar results.


Keywords: Ice strength, Stochastic finite element method, Drucker-Prager, Random field, Monte Carlo simulation
키워드: 빙 강도, 확률론적 유한요소법, 드러커-프라거, 무작위장, 몬테카를로 해석

1. 서 론

지구온난화 등 여러가지 요인으로 인해 북극 해빙이 감소함에 따라 북극 지역에서 운용되는 선박 및 해양구조물의 수요가 꾸준히 증가하고 있으며 이에 따라 설계를 위한 정교한 빙하중 추정의 중요성이 대두되고 있다.

빙하중 추정을 위한 해빙 재료 물성치에 대한 연구는 과거 1960년대부터 이어져왔다. 해빙의 경우 강철과는 다르게 재료 내부의 불확실성에 의한 이방성을 지니게 되는데 이러한 불확실성은 온도, 얼음의 생성 기간, 공극 그리고 염분 등 다양한 요인에 기인한다. 그러나 지금까지의 대부분의 빙하중 추정 연구에서는 위에서 언급한 불확실성을 정량화하기 어려운 이유로 경험에 입각한 안전계수를 사용하여 결정론적으로 취급해 왔다.

Kujala (1997)는 북극해에서 채취한 얼음시편에 대해 4점 굽힘실험을 수행하여 그것의 강도를 구하였다. Ehlers (2013)는 자신의 수치모델에 파티클 스웜 최적화(Particle swarm optimization) 알고리즘을 사용하여 Kujala (1997)의 실험결과를 만족하는 재료 물성치를 결정론적으로 추정하였으나 이는 해빙의 불확실성을 표현하지 못한 한계를 지녔다. 이와 같이 해빙의 재료 물성치에 대한 연구는 대부분 결정론적 접근법으로 행해졌으며 해빙이 내포한 불확실성을 고려한 연구는 전무한 상황이다.

토목과 같은 타 분야의 경우 확률론적 유한요소법을 적용한 해석사례가 조선분야에 비해 비교적 많다. Cadvar et al. (2010)은 진도대교의 탄성계수와 밀도를 무작위장으로 가정한 확률론적 유한요소법을 수행하여 구조물의 설계를 진행하였다. Noh (2004)는 푸와송 비를 확률변수로 설정하여 그에 따른 철판의 응답을 몬테 카를로 시뮬레이션(Monte-Carlo simulation)을 통해 확률적으로계산한 바 있다.

본 연구에서는 확률론적 유한요소기법을 사용하여 얼음시편의 4점 굽힘 강도해석을 시도하였다. 연구의 주된 목적은 해빙의 탄소성 거동을 고려하여 ANSYS 또는 ABAQUS 등의 상용 소프트웨어를 사용한 확률론적 유한요소 해석을 최대한 빠른 시간 안에 수행하는데 있다. 이는 추후 해빙의 재료 물성치를 확률론적으로 추정함에 있어 확률론적 유한요소 해석 시간이 매우 중요하게 작용하기 때문이다. 간략한 연구 내용은 다음과 같다. 먼저, 빙시편의 항복조건으로 드러커-프라거(Drucker-Prager)를 선정한 후 이 항복조건에 사용되는 물성치들에 대한 공간분포를 상관길이(correlation length)를 고려한 다변수 가우시안 무작위장으로 가정하였다. 그 후 확률론적 유한요소 기법인 몬테 카를로 시뮬레이션을 단순화한 direct method와 확률변수에 대한 응답의 민감도를 이용한 sensitivity method를 제안하여 Drucker-Prager 항복조건에 최초로 도달하는데 필요한 힘을 확률분포로 도출해냈다. 위 과정에서 사용된 각각의 재료 물성치의 가우시안 무작위장의 평균 벡터와 상관길이를 다르게 하여 parametric study를 수행하였다.


2. 이론적 배경
2.1 가우시안 무작위장

본 논문에서는 3차원 공간상에 무작위적으로 분포하는 재료 물성치를 수학적으로 표현하기 위해 무작위장이 사용되었다. 무작위장은 시간영역에서 정의되는 무작위과정을 공간상에 정의한 것일 뿐 확률변수의 집합이라는 기본 개념은 동일하다. 무작위장을 연속적으로 표현한 방법으로는 칼후넨-러브(Karhunen-Loeve) 전개법 (Wang 2008)과 스펙트럼으로 표현한 스펙트럴 표현법 (Schinozuka et al. 1991)이 있으며 3차원 공간상의 특정 위치로 이산화 시킨 방법으로는 중점 방법(Midpoint method), 적분점 방법(Integration point method), 평균 이산화 방법(Average discretization method) 그리고 보간 방법(Interpolation method) 등이 있다. (Papadopoulos & Giovanis, 2018) 본 연구에서는 유한요소 상용 소프트웨어에 적용하기 가장 수월한 중점 방법을 사용하였다. 이는 연속적인 무작위장을 유한요소 모델을 구성하는 각 요소의 도심으로 이산화 시키는 방법이다. 3차원 영역 D3에서 정의된 무작위장 {Z(x):xD}은 유한요소의 도심을 나타내는 위치벡터 xD에 할당된 확률변수의 집합이다.

x1, x2, x3,  , xNDZ=[Z(x1), Z(x2), Z(x3),  , Z(xN)]T(1) 

이때, 각 확률변수가 정규분포를 따를 경우 이 무작위장을 가우시안 무작위장이라고 한다. 즉, 가우시안 무작위장은 공간상에 정의된 다변수 정규분포라고 할 수 있다. 각 확률변수의 평균값은 모두 동일한 값을 지니게 되며 평균 벡터 μ로 표현되고 확률변수간의 상관관계는 공분산 행렬 C로 표현된다.

μi=μ(xi)C=C11C1NCN1CNNcij=C(xi,xj),    i,j=1,  , N(2) 

공간상에 정의된 무작위장의 경우 공분산 행렬의 성분을 공간상 두 확률변수 사이의 거리로 정의할 수 있다. 즉 공간상에 존재하는 두 확률변수의 거리가 가까울수록 두 확률변수의 상관관계가 커지며 반대로 멀어질수록 상관관계가 줄어듦을 식 (3)과 같이 표현할 수 있다. 식 (3)σ2은 분산을 의미하며 ι은 두 확률변수의 거리를 표현한 노름(norm)을 정규화하는 상관길이이다. 상관길이가 0일 경우 두 확률변수는 독립이며 상관길이가 커질수록 상관관계가 커지게 된다.

C(xi,xj)=σ2exp-xi-xjι,  i,j=1.  , N(3) 

가우시안 무작위장의 표본추출에는 여러가지 방법이 있으나 본 논문에서는 공분산 행렬 분해 방법 (Fang et al. 2003)을 소개한다. 우선 가우시안 무작위장의 공분산 행렬을 촐스키 분해(Cholesky decomposition) 방법을 통해 어떤 행렬 L과 그것의 전치행렬 LT의 곱으로 표현한다.

C=LLT(4) 

그 후 평균이 0이고 단위 분산을 갖는 정규분포를 따르는 서로 독립적인 N개의 확률변수의 벡터 G를 도입한다.

G=[g1, g2, g3,  , gi]T,   i=1,  , Ngi=N(0,1),   i=1,  , N(5) 

L, μ 그리고 G를 이용하여 가우시안 무작위장 Z식 (6)과 같이 표현할 수 있다.

Z=μ+LG(6) 

식 (6)은 가우시안 무작위장으로부터 표본을 추출하는 문제가 표준정규분포를 따르는 독립적인 N개의 확률변수들로부터 표본을 추출하는 문제로 단순화되었음을 의미한다. 본 연구에서는 정규분포를 따르는 확률변수로부터 표본을 추출하는데 널리 쓰이는 라틴 하이퍼 큐브 샘플링(latin hypercube sampling)방법을 사용하여 벡터 G를 추출하였다.

2.2 드러커-프라거 항복조건

드러커-프라거 항복조건은 최대 비틀림 에너지를 항복예측의 기준으로 삼는 폰-미세스(Von-Mises) 항복기준을 재료 내부적 마찰 또는 점도가 있는 토양, 암반, 얼음과 같은 물질로 확장시킨 탄소성 모델이다 (Choung et al. 2011). 실제 해빙이 인장으로 작용하는 정수압에도 항복을 일으킨다는 점에서 드러커-프라거 캡(Drucker Prager cap) 모델이 항복조건으로 더 적합하나 연구의 초기단계인 점을 고려하여 본 연구에서는 드러커-프라거를 사용하였다. 이럴 경우 인장으로 작용하는 정수압에 의한 항복이 고려되지 않으므로 드러커-프러거 캡 모델을 사용하는 것에 비해 좀 더 보수적인 결과를 얻게 되는 한계가 있다. 드러커-프라거 선형 항복조건은 식 (7)과 같이 표현할 수 있다.

F=t-p tanβ-d=0t=1q1+1K-1-1Krq3p=-13trace(σ)q=32(S:S)S=σ+pI(7) 

식(7)의 K는 3축 압축에서의 항복 응력에 대한 3축 인장에서의 항복 응력의 비를 의미하며 중간 주응력 값에 대한 항복면의 의존성을 제어한다. 본 연구에서는 K값을 0으로 가정하며 그에 따라 t는 폰-미세스 응력 q와 같아진다. p는 정수압 응력을 의미하며 βp-t응력 평면에서 선형 항복 표면의 기울기이며 일반적으로 재료의 마찰 각도(angle of friction)라고 불린다. d는 재료의 내부 점도(cohesion)를 의미한다. 정리하면 βd는 재료 물성치이며 tp는 유한요소 해석으로부터 얻을 수 있는 결과이다.

2.3 확률론적 유한요소법

전통적인 결정론적 유한요소법과 다르게 시스템의 불확실성을 고려한 유한요소법을 확률론적 유한요소법이라고 한다. 여기서의 불확실성이란 외력, 기하학적 형상 그리고 재료 물성치 등의 불확실성을 의미한다. 본 연구에선 이 중 재료 물성치의 불확실성만을 다룬 확률론적 유한요소법을 사용하였다.

P=KUK=K0+K(8) 

전통적인 유한요소법의 경우 강성행렬 K가 결정론적으로 정의된다. 반면 확률론적 유한요소법의 경우 K는 결정론적으로 정의되는 K0와 그 주변의 불확실성을 나타내는 K으로 나뉜다.

K=Nikeke=k0e+keke=ΩeBTD0eBedΩe+ΩeBTDe(x,θ)BedΩe(9) 

국부 강성행렬로 표기하더라도 위와 같은 원리는 똑같이 적용이 되며 이는 국부 강성행렬을 BD로 표현할 때 쉽게 확인 가능하다. B는 변형률과 절점의 변위 사이의 관계를 나타내는 행렬이며 D는 재료의 구성 속성을 표현하는 행렬이다. D0는 재료 물성치의 평균값으로 이루어진 행렬이며 X(x, θ)D0주변의 유동성을 나타내는 평균이 0인 정상 확률과정을 의미한다 (Papadopoulos & Giovanis, 2018).

D(x, θ)=D0[1+X(x, θ)](10) 

확률론적 유한요소법을 푸는 방법으로는 몬테 카를로 시뮬레이션, 퍼터베이션(perturbation) 등이 있다. 퍼터베이션은 주어진 유한요소 문제를 정식화하여 강성행렬과 힘 벡터를 테일러 전개한 후 그것의 응답인 변위 벡터 역시 테일러 전개로 근사하는 방법이다. 그러나 퍼터베이션의 경우 주어진 문제를 정식화 하지 못하는 경우 적용하기 힘든 문제점이 있다.

몬테 카를로 시뮬레이션은 시스템의 무작위성을 표현하는 확률 분포 또는 무작위장으로부터 표본을 추출하여 각 표본마다의 응답을 구해 응답의 통계적 특성을 얻는 방법을 의미한다. 대수의 법칙에 근간하기 때문에 반복 계산을 필연적으로 수반하며 유한요소 모델이 복잡할 경우 계산량의 부담을 안겨준다. 따라서 표본을 추출함에 있어 시간단축이 필수적으로 고려되어야 하며 하나의 표본으로부터 응답을 구하는데 있어서도 역시 시간단축이 중요한 이슈이다.


3. 해석 절차

몬테 카를로 시뮬레이션의 계산량을 줄이기 위해 해빙의 탄소성 거동을 고려할 필요가 있다. Fig. 1에서 확인할 수 있듯이 빙하중은 탄성구간에서 선형적으로 증가하며 항복을 겪은 이후 얼음의 강한 취성으로 인해 급격하게 감소하게 된다 (Kujala et al. 1990).


Fig. 1 
Result of 4-point bending test of ice specimen (Kujala et al. 1990)

설계적인 측면을 고려했을 때, 빙하중의 시간 이력 전체를 아는 것 보다 더 중요한 것은 항복을 겪을 때의 빙하중이다. 이는 곧 전체 시간 이력을 얻기 위한 탄소성 해석 대신 정적 탄성 해석만을 수행하면 됨을 의미한다. 본 연구에서는 위에서 밝힌 해빙의 탄소성 거동을 고려하여 확률론적 유한요소법을 단순화시킨 방법론을 제시하였으며 방법론의 개략적인 흐름은 Fig. 2와 같다.


Fig. 2 
Overall procedure of analysis

3.1 해석 대상 선정

본 논문에서는 Kujala et al. (1990)이 수행한 얼음 시편의 4점 굽힘 실험을 ABAQUS로 모사한 유한요소 모델을 사용하였다. 빙시편을 총 4개의 강체 지지대가 아래 위로 접하고 있는 모델이다. Fig. 3은 ABAQUS를 통해 모델링 되었다. 얼음 시편의 치수는 길이 4,350mm, 높이 400mm 그리고 폭 400mm이다. 얼음을 모사하기 위한 유한요소로는 C3D8R(Reduced integration 3D solid element)가 사용되었으며 얼음 시편을 아래 위로 지지하고 있는 지지대는 강체로 모델링 되었다. 얼음 시편을 이루는 1,000개의 요소와 2,222개의 노드로 구성되었다.


Fig. 3 
ABAQUS FEM model

3.2 물성치의 무작위장 생성 및 표본추출

본 연구에서는 드러커-프라거 항복 조건에 사용되는 물성치인 마찰각도 그리고 내부 점도에 대한 무작위장을 정의하였다. 각 무작위장은 3차원 가우시안 무작위장이며 서로 독립적이다. ABAQUS의 유한요소 모델로부터 각 요소의 도심의 정보를 가져와 각 물성치의 무작위장을 얼음 시편 요소의 도심으로 이산화시켰다. 3.1에서 선정된 유한요소 모델은 1,000개의 요소로 이산화되었으며 따라서 각 요소의 도심에 정의된 가우시안 무작위장의 확률변수 역시 1,000개이다.

각각의 무작위장에 대한 평균 벡터와 표준편차 값은 Table 1을 기준으로 약간의 변화를 주어 정의되었으며 추후 parametric study에 사용되었다. Table 1은 본 연구에서 수행된 parametric study 중 Zhang et al. (2017)에서 사용한 값을 그대로 사용한 경우를 의미한다. 상관길이, 마찰각도 그리고 내부 점도의 무작위 장의 모수에 대한 기존 연구 및 추정 방법이 정립되지 않은 점을 고려하여 본 연구에서는 임계하중의 확률분포가 양의 영역에서 구해지는 값을 모수로 택하였다. 상관길이의 경우 모든 확률변수가 완전히 독립적인 경우를 의미하는 0m, element의 크기 0.25m보다 작은 0.1m, 그리고 element 크기보다 큰 1m가 사용되었다. 마찰각도의 경우 그 값이 커질 경우 식 (10)을 사용하여 임계하중을 구함에 있어서 분모가 음의 값을 나타나게 되어 임계하중 값 자체가 음수가 나올 수 있으므로 Zhang et al. (2017)에서 사용한 36˚를 기준으로 6˚씩 줄여가며 사용하였다. 내부 점도의 경우 식 (10)에서 확인할 수 있듯이 임계하중의 부호에 영향을 미치지 못하므로 Zhang et al. (2017)에서 사용한 0.58 MPa을 기준으로 50%씩 그 값을 증가시키며 사용하였다. 각 경우에 대하여 표준편차는 평균벡터의 5%로 설정을 하였다.

Table 1. 
Mean and standard deviation of angle of friction and cohesion
Material Property μ(xi) σ
Angle of friction 36˚ 1.8˚
Cohesion 0.58 MPa 0.029 MPa

2.1에서 밝혔듯이 상관길이는 무작위장 위의 두 확률변수 사이의 상관관계를 나타내는 값이다. Fig. 4는 마찰각의 상관길이를 변화하여 그로부터 추출된 표본의 그림이다. 상관길이가 작아질수록 두 확률변수는 독립에 가까워지므로 등고선의 기울기가 가파름을 알 수 있으며 반대로 상관길이가 커질수록 두 확률변수의 상관관계가 커져 등고선의 기울기가 완만해짐을 확인할 수 있다.


Fig. 4 
Sample realization of random field of angle of friction

3.3 파손 조건 정의

위에서 언급한 바와 같이 해빙의 취성을 고려하여 마찰각과 점도의 가우시안 무작위장으로부터 추출된 표본마다 유한요소 중 어느 것이라도 드러커-프라거의 항복면에 최초로 도달하는 순간을 해당 표본의 파손 조건으로 정의하였으며 그 때의 하중을 임계하중으로 정의하였다.

3.4 임계하중 계산

본 연구에서는 각 표본에 대한 임계하중을 구하는 방법으로 두가지 방법론을 제시하였다. 첫번째는 direct method로 일종의 몬테 카를로 시뮬레이션이다. 정적 탄성 해석을 수행하기 위해 모델링 과정에서 탄성계수와 푸와송비를 설정해 주어야 한다. 이때 정적 탄성 해석의 경우 탄성계수와 응력 응답에 영향을 끼치지 못하므로 Table 2Zhang et al. (2017)의 값을 그대로 차용하였다.

Table 2. 
Material property of Zhang et al. (2017)
Young's modulus 4.5 GPa
Poisson's ratio 0.3

그 후 얼음시편 밑에 위치한 두개의 강체 지지대에 얼음시편 길이방향의 수직방향으로 단위 힘을 주어 얼음시편의 굽힘을 유도한 다음 그 결과인 폰-미세스 응력 q와 정수압 p를 모든 요소에 대하여 추출하였다. 단위 힘을 주었을 때의 해석 결과는 Fig. 5와 같다. Fig. 5의 (a)와 (b)는 각각 폰-미세스 응력과 정수압의 결과를 나타낸다. 폰-미세스 응력의 경우 최소 0.163Pa에서 최대 124.998Pa의 분포를 보이며 정수압의 경우 최소 -0.440Pa에서 최대 58.402Pa의 분포를 보인다.


Fig. 5 
Result of static elastic analysis for 4-point bending

이때 qp는 단위 외력에 대한 응답을 의미하므로 정적 탄성 해석에서의 응력 결과가 외력에 비례함을 이용하면 q·xp·xx의 힘으로 외력으로 작용할 때의 응력 결과임을 알 수 있다. 항복함수 식 (7)에 폰-미세스 응력과 정수압에 x를 곱해준 값을, 그리고 마찰각과 점도의 가우시안 무작위장으로부터 추출된 표본을 요소별로 대입하면 식 (10)과 같다.

F=q·x-p·x·tanβ-d=0x=3q-p·tanβ(10) 

마찰각과 점도의 가우시안 무작위장의 특정 표본에 대하여 각 요소마다 식 (10)을 계산함으로써 각 요소의 임계하중을 계산할 수 있으며 이 중 최소값은 곧 특정 표본의 임계하중을 의미한다. 위 과정을 모든 표본에 대해 수행하면 모든 표본에 대한 임계하중을 계산할 수 있다.

두번째 방법은 sensitivity method이다. 본 연구에서 사용된 항복조건인 드러커-프라거의 경우 임계하중을 구하기 위한 식 (10)이 양함수의 형태를 갖지만 음함수의 형태를 갖는 다른 항복조건을 사용하게 될 경우 direct method를 적용하는데 제한이 있다. 음함수를 풀기 위해 수치해를 구해야 하지만 Python 또는 MATLAB에 내장된 수치계산 솔버(numerical solver)의 경우 양함수인 식 (10)에 값들을 대입하여 임계하중을 구하는 시간에 비해 100,000배가량 많은 시간이 소요된다. 이러한 계산시간의 문제를 해결하기 위한 대안이 sensitivity method이다.

Sensitivity method는 테일러 급수를 이용한 방법이다. 식 (10)의 폰-미세스 응력 q와 정수압 p는 마찰각과 점도의 모든 가우시안 무작위장 표본에 대하여 같은 값을 갖으므로 상수로 취급되며 따라서 임계하중 x는 마찰각 β과 점도 d의 함수다. 그러므로 임계하중 x는 마찰각 β과 점도 d의 테일러 전개로 표현 가능하다. 이때 βd에 대한 x의 도함수는 Python 또는 MATLAB을 통해 쉽게 구할 수 있다. 본 연구에서는 3차 도함수까지 근사하였는데 이는 1차까지 근사할 경우 임계하중 값이 음의 값이 나올 정도로 오차가 크게 나오기 때문이며 2차까지 근사할 경우 여전히 약 5% 정도의 오차를 보이기 때문이다. Sensitivity method를 통해 두번째 표본의 첫번째 요소의 임계하중은 다음과 같이 구할 수 있다.

x(d1(2), β1(2))=x(d1(1), β1(1))+xd(d1(1), β1(1)) (d1(2)-d1(1))+xβ(d1(1), β1(1)) (β1(2)-β1(1))+12!xdd(d1(1), β1(1)) (d1(2)-d1(1))2+2xdβ(d1(1), β1(1)) (d1(2)-d1(1)) (β1(2)-β1(1))+xββ(d1(1), β1(1)) (β1(2)-β1(1))2(11) 

식 (11)는 편의 상 2차까지 근사한 식이다. xβxd는 마찰각과 점도에 대한 임계하중의 1계 편미분을 의미하며 동시에 마찰각과 점도에 대한 임계하중의 민감도를 의미한다. 식 (11)의 아래첨자는 요소의 인덱스를 의미하며 윗첨자는 표본의 인덱스를 의미한다. 마찰각과 점도의 가우시안 무작위장의 모든 표본은 3.2에서 가우시안 무작위장 표본 추출로 이미 획득되었다. 따라서 두번째 표본의 첫번째 요소의 임계하중은 첫번째 표본의 첫번째 요소의 임계하중에 표본의 변화에 따른 마찰각과 점도의 변화에 대한 임계하중의 변화를 더하여 구할 수 있다. 결국 첫번째 표본에 대한 계산만 수치계산 솔버로 계산될 수 있다면 그 다음 표본부터는 sensitivity method를 통해 그 해를 구할 수 있다.

본 연구에서 사용한 드러커-프라거처럼 양함수가 아닌 음함수의 형태를 갖는 항복 조건을 사용할 경우 임계하중의 도함수를 구해보면 임계하중을 우변에 포함하는 형태로 도함수를 얻을 수 있다. 이럴 경우 테일러 전개에 사용되는 도함수 안의 임계하중은 이전 표본의 임계하중 이므로 쉽게 처리가 가능하다. 따라서 양함수와 마찬가지로 첫번째 표본의 임계하중만 구할 수 있다면 그 다음 표본의 임계하중부터는 sensitivity method를 통해 손쉽게 얻을 수 있다. 본 연구의 경우 음함수의 경우 MATLAB의 수치계산 솔버로 계산하는 것이 sensitivity method로 계산하는 것보다 약 50,000배의 시간이 소요되는 것을 확인할 수 있었다.

3.5 통계적 후처리

3.4에서 제시한 direct method와 sensitivity method를 통해 얻게 되는 것은 각 표본에 대한 임계하중이다. 이를 사용하여 히스토그램 및 경험적 확률분포를 얻을 수 있으며 간단한 통계적 추론을 사용하여 평균 및 표준편차와 같은 임계하중에 대한 확률분포의 모수를 추정할 수 있다.


4. Parametric study

본 연구에선 앞서 밝힌 두 방법론이 어떻게 작동하는지 보기위해 상관길이, 점도 그리고 마찰각에 대하여 parametric study를 수행하였다. 각 경우에 대해 사용한 표본 수는 100,000이며 그 결과는 Table 3, Table 4, Table 5와 같다.

Table 3. 
Parametric study about correlation length
Reaction Force [N]
Mean Std_dev
Correlation length = 0 Direct 6,904.096 150.850
Sensitivity 6,904.261 150.737
Correlation length = 0.1 Direct 6,970.551 181.126
Sensitivity 6,970.634 181.057
Correlation length = 1 Direct 7,150.651 224.973
Sensitivity 7,150.918 224.737

Table 4. 
Parametric study about cohesion
Reaction Force [N]
Mean Std_dev
Cohesion = 0.58 MPa Direct 6,970.551 181.126
Sensitivity 6,970.634 181.057
Cohesion = 0.87 MPa Direct 10,451.796 272.134
Sensitivity 10,451.828 272.099
Cohesion = 1.16 MPa Direct 3,940.818 363.049
Sensitivity 3,941.094 362.823

Table 5. 
Parametric study about angle of friction
Reaction Force [N]
Mean Std_dev
Friction Angle = 24˚ Direct 7,496.481 218.386
Sensitivity 7,496.505 218.368
Friction Angle = 30˚ Direct 7,236.568 200.615
Sensitivity 7,236.588 200.593
Friction Angle = 36˚ Direct 6,970.551 181.126
Sensitivity 6,970.634 181.057

Fig. 6는 마찰각 36도, 점도 0.58MPa, 상관길이 0.1m일 때의 임계하중의 경험적 확률분포이다. Fig. 5를 포함한 parametric study의 모든 경우에서 100,000번의 표본추출로 충분히 정규분포의 통계적 특성을 찾아낼 수 있었다.


Fig. 6 
PDF of critical force


5. 결 론

본 연구에서는 공간상 랜덤하게 분포하는 재료 물성치를 고려한 얼음시편의 확률론적 강도 해석 방법을 제시하였다. 연구 결과에 대한 분석을 토대로 아래와 같은 결론을 도출할 수 있었다.

• 먼저, 3차원 공간상에 무작위적으로 분포하는 재료 물성치를 가우시안 무작위장을 사용하여 표현할 수 있었다. 이때 재료 물성치는 드러커-프라거에 사용되는 점도와 마찰각이다. 가우시안 무작위장으로부터 표본을 추출하기 위해 무작위장의 표본추출 문제를 독립 동일 표준정규분포를 따르는 확률변수의 표본추출 문제로 바꾸었다. 이를 위해 촐스키 분해와 라틴 하이퍼 큐브 샘플링 기법이 적용되었다.

• 항복 이후 급격한 소성을 겪는 해빙의 탄소성 거동을 고려하여 몬테 카를로 시뮬레이션의 반복계산의 계산량을 단축시킨 direct method을 제안하였다.

• Direct method의 경우 음함수 형태를 갖는 임계하중을 계산하는데 비교적 많은 시간이 걸리는 문제점을 갖는다. 따라서 계산 시간을 단축하기 위해 테일러 전개를 기반으로 하는 sensitivity method를 제안하였다.

• Direct method와 sensitivity method를 사용하여 점도, 마찰각 그리고 상관길이에 대하여 parametric study를 수행하였다. 표본수는 두 방법 모두 표본 수는 100,000이며 계산 시간은 direct method의 경우 10초, sensitivity method의 경우 3분 내외였다. 점도와 상관길이가 증가함에 따라 임계하중의 평균값과 표준편차가 증가함을 확인할 수 있으며 반대로 마찰각의 경우 임계하중의 평균값과 표준편차가 감소함을 확인할 수 있다. 또한 모든 경우에 대하여 direct method와 sensitivity method가 매우 비슷한 결과를 계산해주는 것을 확인할 수 있었다.

• 본 연구는 재료 물성치의 공간분포를 고려하여 단 한 번의 탄성해석을 통해 얼음 시편이 파손을 일으키는 임계하중의 확률분포를 구한다는 점에서 기존 시계열 해석과 차이를 보인다. 평탄빙에 의한 빙하중과 같은 문제에 적용할 경우 전체 시계열을 계산할 필요가 없다는 점에서 계산 시간의 큰 단축을 기대할 수 있다.

• 그러나 선형성이 강한 bending과 같은 해석에만 사용할 수 있다는 한계가 있으며 본 연구에서 필요한 재료 물성치들의 무작위장의 모수를 추정하기 위한 방법론이 정립되지 않은 문제점이 있다.

• 추후 계획으로는 본 연구에서 제시한 direct method와 sensitivity method를 이용한 드러커-프러거 재료 물성치의 가우시안 무작위장 모수 추정이 있다. 이를 위해 얼음시편 굽힘 실험을 계획 중에 있으며 실험결과로부터 얻을 수 있는 임계하중의 확률분포와 direct method 또는 sensitivity method로부터 얻을 수 있는 임계하중의 확률분포의 차이를 줄여 나가는 방식의 베이지안 추론(Bayesian inference)을 통한 모수의 확률분포 추정을 연구 중에 있다.

• 본 연구에서 제시한 임계하중의 확률분포를 얻는 방법론과 추후에 연구할 재료 물성치의 가우시안 무작위장 모수 추정을 통해 평탄빙을 운항하는 빙해선박의 빙하중을 확률적으로 추정하는데 최종적인 목표가 있다.


Acknowledgments

이 논문은 2019년도 정부(산업통상자원부)의 재원으로 한국산업기술진흥원의 지원을 받아 수행된 연구임(No, P0001968, 2019년 산업전문인력역량강화사업). 이 논문은 산업통상자원부 ‘산업전문인력역량강화사업’의 재원으로 한국산업기술진흥원(KIAT)의 지원을 받아 수행된 연구임(2019년 한-영 해양플랜트 글로벌 전문인력 양성사업, 과제번호: N0001287).


References
1. Çavdar, Ö., Bayraktar, A., & Adanur, S., 2010. Stochastic finite element analysis of a cable-stayed bridge system with varying material properties. Probabilistic Engineering Mechanics, 25(2), pp.279-289.
2. Choung, J.M., Nam, J.M., & Kim, K.S., 2011. Comparative study on material constitutive models of ice. Journal of the Society of Naval Architects of Korea, 48(1), pp.42-48.
3. Ehlers, S. & Kujala, P, 2014. Optimization-based material parameter identification for the numerical simulation of sea ice in four-point bending. Proceedings of the Institution of Mechanical Engineers. Part M: Journal of Engineering for the Maritime Environment, 228(1), pp.70-80.
4. Papadopoulos, V., & Giovanis, D. G., 2017. Stochastic finite element methods: An introduction. Springer.
5. Timco, G. W. & Weeks, W. F, 2010. A review of the engineering properties of sea ice. Cold regions science and technology, 60(2), pp.107-129.
6. Wang, L., 2008. Karhunen-Loeve expansions and their applications. Ph.D. London School of Economics and Political Science, United Kingdom.
7. Zhang, N., Zheng, X. & Ma, Q, 2017. Updated smoothed particle hydrodynamics for simulating bending and compression failure progress of ice. Water, 9(11), 882.