물성치의 공간분포를 고려한 빙 시험편의 확률론적 강도평가

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차원 영역 $$ D\subset {ℝ}^3$$에서 정의된 무작위장 $$ \left\{Z\right(x):x\in D\}$$은 유한요소의 도심을 나타내는 위치벡터 $$ x\in D$$에 할당된 확률변수의 집합이다.

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

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

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

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

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

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

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

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

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

2.3 확률론적 유한요소법

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

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

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

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

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

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

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

Fig. 4

Sample realization of random field of angle of friction

3.3 파손 조건 정의

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

3.4 임계하중 계산

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

Table 2.

Material property of Zhang et al. (2017)

그 후 얼음시편 밑에 위치한 두개의 강체 지지대에 얼음시편 길이방향의 수직방향으로 단위 힘을 주어 얼음시편의 굽힘을 유도한 다음 그 결과인 폰-미세스 응력 $$ 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

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

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

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

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

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

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

3.5 통계적 후처리

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