J. Korean Soc. Hazard Mitig Search

CLOSE


J. Korean Soc. Hazard Mitig. > Volume 17(2); 2017 > Article
비선형 점성유체의 다상유동 모형을 이용한 용암류 전산해석

Abstract

Due to high temperature and viscosity, lava flow can cause large-scale damage on surrounding terrain and humanity. It is formed by molten rocks when magma in volcano is erupted out of its chamber beneath earth surface. The physical property of lava is the nonlinear behavior of viscosity of molten rocks with the influence of temperature-varying properties. Because of this, experimental estimation in the laboratory is nearly impossible. Thus many efforts have been carried out to numerical approach but mostly numerical method is little availability. This study presents a numerical approach to simulate the unsteady flow of molten rocks effusing from a vent using open source code of computational fluid dynamics utility called OpenFOAM with two phase modeling for lava and surrounding atmosphere. The viscosity of lava flow is simulated the Bingham plastic or Herschel-Bulkley model. In order to validate the numerical method, we carried out the quantitative evaluations including the extent of damage, travel distance and thickness of the lava flow at Mt. Etna, Italy.

요지

용암은 용융된 암석이 분출되어 지표면을 흐르는 현상으로써 높은 온도와 점도 등에 의해 주변 지형과 인류에게 막대한 피해를 유발할 가능성이 있는 자연재해다. 이때 용암의 물리적 성질은 온도에 따라서 점성이 변하는 비선형적 거동 특성을 보인다. 이러한 극한 온도 상태에서의 유동흐름은 실험으로 추정하는데 한계가 있으므로 수치적 접근 방법에 의해 다수 연구가 수행되었다. 이에 본 연구에서는 소스가 공개된 전산유체해석 기법인 OpenFOAM을 이용하여 용융된 암석이 시간에 따라서 분출되는 비정상 흐름현상에 대해서 수치해석을 수행하였다. 여기서 사용된 수치해석 기법은 이상유동으로 유체의 자유표면을 추적하는 VOF기법을 적용하였으며, 용암의 점성은 비선형 점성 모델인 Herschel-Bulkley 유체를 적용하였다. 이와 같은 수치해석에 의한 용암의 거동을 정량적으로 평가하기 위해서 2001년 분화한 이탈리아의 에트나 화산의 용암류의 관측치를 대상으로 용암의 이동거리, 두께, 확산범위 등을 비교하였다.

1. 서론

용암류는 높은 온도에서 용융된 암석이 지표로 분출되어 흐르는 현상이다. 암석의 조성 물질에 의해 용암의 점도, 유속, 두께 등이 다양하게 나타난다. 일반적으로 이산화규소(SiO2)의 양이 적은 현무암질 용암이 유동성이 크고 유속이 비교적 빠르므로 광역적인 피해를 일으키며, 용암의 성분에 따른 주요 특징은 Table 1에 나타내었다. 그리고 용암의 온도로 인해서 암석이 경화 혹은 용융현상이 발생하며, 특히 외부와 직접 맞닿아 있는 용암의 표면과 지표면에서의 전도 및 대류를 통한 열손실로 인해 온도 분포가 결정되며, 이때 용암이 냉각되는 과정에서 용암의 점도, 속도에 영향을 미쳐 유동형태에 변화를 준다(Cordonnier et al., 2015). 또한, 밀도가 높은 고온의 암석이므로 지형이 변하거나 산불이 발생하는 등 2차 피해의 우려가 있다.
Table 1
Feature of Lava Flow According to Composition(Kilburn, 2000)
Composition SiO2<55% (Basalt, Basaltic andesites) SiO2>55% (Andesites, dacites, trachytes)
Length(km) <10 (basalt to 50) <5 (some to 15)
Mean thickness(m) 3~20 20~300
Mean discharge rate(m3/s) 10~100 (basalt to 1,000) 1~10 (andesites to 100)
Flow volume(km3) 0.01~0.1 0.01~1.0
Temperature(°C) 1,050~1,200 (basalt) 950~1,170 (andesite)
Fig. 1은 2007년 이탈리아의 에트나(Etna) 화산에서 분출된 용암으로 이로 인해 집이 매몰되는 등 지표면에 상당량이 퇴적된 바 있다. 따라서 용암의 이동 경로와 확산 범위를 사전에 예측함으로써 피해를 줄이는 방안이 필요하다. 그러나 비선형 점성 거동 특성으로 인해 용암의 유동 특성을 파악하는 데 어려움이 많다.
Fig. 1
House Buried in Lava from Etna Volcano, Sicily
KOSHAM_17_02_117_fig_1.gif
특히, 실험으로 재현하기에 여러 제약이 있으므로 주로 댐 붕괴 실험과 같이 단순한 모형이나 수치적 접근에 의해 연구되고 있다. 댐 붕괴 실험의 경우 실리콘, 시럽 등 용암과 유사한 높은 밀도의 유체를 이용하여 수평 이동거리를 분석하는 등 용암을 분석하기 위한 기초 연구로 다양한 조건으로 수행되고 있다(Castruccio et al., 2010; Minussi et al., 2012; Saramito et al., 2013; Castruccio et al., 2014). 또한, 용암의 특성을 반영하기 위해 고온에서 암석을 용융시킨 후 평판 혹은 경사면에서의 흐름 실험을 수행하였지만 실제 흐름에 적용하기 위해서 지형조건을 고려한 검증이 요구된다(Lev et al., 2012). 이와 같이 이상의 모형실험은 지형조건 등의 복잡한 거동 현상을 재현하는데 한계가 있다. 이러한 이유로 용암류를 분석하기 위해서 수치적 접근 방법이 많이 발전하였다. 주로 CA(Cellular Automata) 기법을 사용하며, 이는 지형 경사를 반영하여 고도가 낮은 방향으로 확산시키는 방법으로써 SCIRA, MAGFLOW 등이 있다. Crisci et al. (2004)는 SCIRA를 이용해서 1991~2001년 에트나 화산의 용암류의 이동범위를 정성적으로 비교한 바 있다. Vicari et al.(2007)은 2차원 CA 기법의 MAGFLOW를 이용해서 2001년 에트나 화산의 용암류를 재현하였고 이를 관측 값을 기준으로 이동거리, 면적에 대해서 정성적·정량적인 비교를 통해 유효성을 검증하였지만 2차원 기법이기 때문에 수직분포에 대한 고려가 배제되었다. Cordonnier et al.(2015)은 MAGFLOW, VOLCFLOW 등 여러 용암류 수치해석 프로그램을 사용해서 평판과 경사면에서의 유동, 온도 확산 등 다양한 경우에서 검증을 수행한 후 2001년 에트나 화산의 이동 범위를 비교·분석한 바 있다. 그러나 이상의 수치해석 코드는 2차원 영역에서 수행되어 빠른 계산속도와 정성적인 검증을 통해 널리 사용되었지만 용암의 두께와 열에너지 등 3차원 영역에서의 수직분포를 반영하지 못하는 단점이 있다. 따라서 본 연구에서는 3차원 전산유동해석의 다상유동 기법을 이용해서 용암류의 수치해석에 접근하였다.
용암류는 화산이 분화하는 경우 발생하는 근접재해로써 분화구와 인접한 지역에서 막대한 피해가 발생한다. 특히, 백두산은 2002년 이후 분화의 전조 현상이 나타난 바 있으며, 946년경으로 추정되는 화산폭발지수 7에 해당하는 대규모 분화 기록 등으로 인하여 동북아시아 지역의 대규모 화산재해 가능성이 연구되고 있다(Liu et al., 2011; Yun et al., 2012). 또한 최근 연구에 의하여 백두산 지하에 존재하는 마그마 방이 확인되고 있고 20억 톤의 수량을 갖는 천지로 인해서 화산 분화시 대규모 용암류나 화산이류가 발생할 가능성이 높을 것으로 판단된다(Yun et al., 2012). 본 논문에서는 비선형성을 지닌 용암류에 대해서 공개 소스 코드를 적용하여 백두산 분화시에 발생할 수 있는 근접재해에 대한 분석뿐만 아니라 국내에서 빈번하게 발생하는 토석류와 이류 등에 적용할 수 있는 수치해석 방법을 제기하고자 하였다.
본 연구는 화산의 근접재해 중에서 용암류를 대상으로 하였고, 유한체적법 기반으로 Navier-Stokes 지배방정식을 해결하는 전산유동해석을 활용하였다. 특히, 기존 수치해석 방법을 보완할 수 있도록 3차원 공간을 구현함으로써 수직분포의 영향을 고려할 뿐만 아니라 비선형 점성유체의 물성을 더욱 상세하게 반영한 수치해석을 수행함으로써 용암의 특성을 정량적으로 검증할 수 있는 방안을 제시하였다. 이때 소스가 공개된 전산유체해석 라이브러리인 OpenFOAM을 사용하였고 이는 사용자 임의로 코드를 수정하거나 추가할 수 있는 장점이 있다. 본 연구에서 사용한 OpenFOAM의 표준 솔버인 interFoam은 해석 영역 내에서 용암과 대기 영역으로 가정한 이상(Two-phase) 유동을 고려하여 두 유체의 상호작용을 동시에 처리하는 VoF(Volume of Fluid) 기법이 적용되었다. 그러나 등온 상태의 유체를 가정하기 때문에 열에너지방정식을 결합함으로써 온도의 영향을 고려한 용암류를 구현하였다.

2. 지배방정식과 수치해석 기법

2.1 지배방정식

본 연구에서는 3차원 지형을 기반으로 지표면에서 흐름이 발생하는 비압축성, 비정상유동을 가정한 용암류를 모의하였다. 이를 구현하기 위한 수치모형은 Ubuntu 14.04의 운영체제에서 공개 코드인 OpenFOAM-v2.4.0을 사용하였다.
전산영역 내에서 용암류를 모사하기 위한 방법으로 두 유체의 상 경계를 추적하여 유동 현상을 동시에 처리하는 이상(Two-phase)유동해석 중에서 VoF기법을 적용하였다. VoF기법은 검사체적 내에서 유체의 체적분율(Volume fraction, α)이 0 혹은 1인 경우로 두 유체의 상을 구분하며, 자유표면은0< α <1에 위치한다. 여기서 유체의 점성과 밀도는 물성이 다른 각 유체를 고려하여 Eqs. (1)-(3)과 같이 산정된다.
(1)
ρ=α1ρ1+α2ρ2
(2)
μ=α1μ1+α2μ2
(3)
α2=1α1
여기서μ는 점성계수(μ= ρ·v)이고 아래첨자 1과 2는 상이 다른 두 유체를 의미하며 본 연구에서는 공기와 용암이다.
이와 같이 이상 유동해석 기법을 이용하여 유동 문제를 해결하기 위한 지배방정식은 Eqs. (4)-(6)과 같이 연속방정식, 운동량 방정식, 그리고 수송방정식이다.
(4)
·U=0
(5)
 ρUt +·(ρUU)·(νU)=ρgP+F
(6)
 αt +·(Uα)=0
여기서U는 속도, ρ는 유체의 밀도, v는 동점성계수, g는 중력가속도, P는 압력, 그리고F는 표면장력을 의미한다.
한편, Eq. (6)에서 유체의 체적분율에 의해 구분된 두 유체사이의 경계에서 발생하는 상대속도를 고려하여 수정된 수송방정식은 Eq. (7)과 같다.
(7)
 α1t +·(Uα1)+·(α2Uc)=0
여기서Uc는 두 유체의 상대속도의 벡터성분을 의미하며, 통상 Compression velocity라고 한다.
본 연구는 공개소스 라이브러리 OpenFOAM에서 Eqs. (1)-(7)과 같은 지배방정식을 이용한 표준 수치해법인 interFoam에 Eq. (8)의 열에너지 방정식을 결합하여 수정한 솔버를 사용하였다.
(8)
ρCP( Tt +U·T)=·(κT)+ϕ
여기서CP는 열용량이고k는 열확산율을 의미한다. 이 때 열확산율은k = ρCP/Pr로 정의되고Pr은 유체의 동점성과 온도 전도의 관계를 의미하는 프란틀 수(Prandtl number)다.

2.2 비선형 점성모델

현무암질 암석이 1000°C 이상의 고온에서 용융되어 유동성을 갖는 용암류의 전단응력은 전단변형률과 비선형 관계를 갖는 비뉴턴유체이다. 이러한 성질을 갖는 용암류를 모사하는 방법에는 Bingham plastic 모델과 Herschel-Bulkley 모델 등이 있다(Griffiths, 2000). 이 두 모델은 항복응력과 점성이 중요한 매개변수로 작용하며, 항복응력에 도달하기 전에는 고체와 같이 거동하지만 항복응력에 도달하는 경우 유체와 같이 유동성을 지니는 특징이 있다. 한편, Bingham plastic 모델은 항복응력 이후의 거동은 Newtonian의 성질을 띠어 전단응력과 전단변형률이 선형관계를 갖는 반면 Herschel-Bulkley 모델의 경우 의가소성유체(Shear-thinning) 혹은 팽창성유체(Shear-thickening)의 성질을 갖는다.
본 연구에서는 이상의 비뉴턴 유체의 모델 중에서 비교적 상세한 물성을 반영할 수 있는 Herschel-Bulkley 모델을 이용하였으며, 이는 Eq. (9)와 같이 정의된다.
(9)
τ=τ0+kγ˙n
여기서τ는 전단응력, τ0는 항복응력, ϒ는 전단변형률, k는 비례상수를 의미하는 컨시스턴시 지수(Consistency index) 그리고n은 유체의 상태를 나타내는 유동지수(Flow index)다.
이와 같은 비선형 점성모델을 수치기법에 적용하기 위해서 Eq. (10)과 같이 비뉴턴유체의 겉보기점성(η)을 도입하여 뉴턴유체와 같이 선형의 형태로 식을 변환한 후η는 Eqs. (11)-(12)와 같이 적용한다.
(10)
τ=η(γ˙)γ˙
(11)
η={ν0,|γ˙|γ˙0τ0|γ˙|+k|γ˙|n1|γ˙|γ˙0
(12)
η=min(ν0,τ0/γ˙+kγ˙n1)
이 때 겉보기점성은 Eq. (11)과 같이 전단변형률에 의해 점성이 결정된다. 이와 같이 정의된 비선형 점성모델의 매개변수(v0, τ0, k, n)는 대상 유체인 용암의 물성에 맞게 적용된다.

2.3 수치해석 기법

용암류의 수치해석을 위해서 사용된 식(4-8)의 지배방정식을 해석영역 내에서 전산 격자를 통해 물리량을 산정하기 위한 이산화방법에는 유한체적법이 사용되었다. 그리고 각 항에 적용된 기법은 다음과 같다. 먼저 시간항은 1차 정확도의 Euler 기법, 대류항은 1차 정확도의 풍상측 차분법(Upwind differential), 유체의 상은 2차 정확도의 vanLeer 그리고 압력구배는 선형보간법을 사용하였다. 또한, 지배방정식은 압력기반 해법을 사용함으로서 운동량과 압력이 주요 변수로 작용되며, 이를 계산하기 위해서 하나의 지배방정식으로 결합하는 PISO(Pressure implicit with splitting of operators) 알고리즘을 사용하였다.

3. 댐 붕괴 현상을 이용한 비선형 점성유동 해석 검증

3.1 개요

복잡한 지형과 물성 등을 고려한 용암류를 분석하기 전에 유동의 비선형성과 수치모형의 검증을 위해서 비교적 단순한 형태의 댐 붕괴 모형을 이용하여 유체의 거동을 분석하였다. 본 연구에서는 비선형 점성유체를 대상으로 댐 붕괴 해석을 수행하여 실험과 수치해석 결과를 비교함으로써 본 연구에서 사용한 수치기법의 유효성을 검증하였다.

3.2 비선형 점성유동의 해석영역과 경계조건

댐 붕괴 모형을 수치적으로 접근하기 위해서 Fig. 2와 같이 1,900mm × 145mm × 5mm의 해석영역을 설정하였고 11,020개의 정육면체 격자를 사용하였다. Fig. 2에서 초기 설정된 유체영역이 중력(-y 방향)의 영향으로 붕괴되는 것을 모사하기 위해서 상단은 대기압 조건이며, 이외의 각 단면은 점착조건을 설정하였다. 이때 유체영역의 높이(H)는 각각 70, 100, 130mm를 초기조건으로 사용하여 각 경우에 대해서 유체의 이동거리를 비교하였다. 여기서 유체영역은 Minussi et al.(2012)의 실험에서 사용된 비선형 유체인 Carbopol 940을 모사하기 위해서 본 연구에서는 Herschel-Bulkley 모델을 사용하였고 매개변수는 Table 2와 같다.
Fig. 2
Computational Domain of Dam-break Model
KOSHAM_17_02_117_fig_2.gif
Table 2
Herschel-Bulkley Parameters Of Carbopol 940(Minussi et al., 2012)
 Properties  Samples  Mean value 
1 2 3
τ0  17.79   18.83   18.10  18.24
k 1.74 1.85 2.12 1.90
n 0.54 0.54 0.52 0.53
ρ = 1000kg/m3

3.3 비선형 점성유체의 이동거리 분석

Minussi et al.(2012)는 비선형 점성유체의 이동거리를 분석하기 위해서 댐 붕괴 실험과 더불어 상용 소프트웨어 Ansys CFX의 Herschel-Bulkley 모델을 사용하였다. 본 연구에서는 동일한 조건에서 OpenFOAM을 이용하여 선행연구와의 비교를 통해서 수치기법을 검증하였다. Fig. 3에 선행연구와 본 연구의 결과를 시간에 따른 유체의 이동거리로 나타내었다. 그리고 유체의 초기 높이에 따라서 (a)~(c)에 각각 나타내었으며, 유체의 이동 양상이 매우 유사한 것을 확인하였다. 특히, (b)와 (c)의 경우 실험의 이동 거리와 잘 일치하였다. 이상 실험과 수치해석 간의 비교를 통해서 본 연구에서 사용한 수치 기법이 유체의 거동을 잘 모사하고 있음을 알 수 있다.
Fig. 3
Front Horizontal Position Versus Time-varying
KOSHAM_17_02_117_fig_3.gif

4. 열 유체 분출 현상을 이용한 열대류⋅전도 해석 검증

4.1 개요

비선형 점성유체의 다상유동 수치기법을 댐 붕괴 모형을 이용해서 유효성을 확인한 후 고온의 유동체인 용암류를 구현하기 위해서 이상의 수치기법과 열에너지 방정식을 결합하였다. 이와 같이 수정된 수치기법을 적용하여 온도에 의한 유체의 거동을 검증하기 위해서 열 유체의 분출 현상을 모의하여 열대류와 전도에 관한 수치해석을 수행하였다.

4.2 열 유동의 해석영역과 경계조건

Fig. 4와 같이 뜨거운 유체가 유입되어 평판 위에서 흐름이 발생하는 실험이며, 유체의 이동과 온도의 확산에 관해서 정량적 수치를 비교하였다(Garel et al. 2012). 이때 대기 온도는 20°C, 유입되는 유체의 온도는 42°C이고 8,100초 동안 2.2 × 10-8 m3/s로 유입된다. 여기서 Eq. (13)의 푸리에법칙(Fourier’s law)에 따라서 대기와의 온도 차에 의해 열 이동이 발생한다.
Fig. 4
Computational Domain of Thermal Diffusion Model
KOSHAM_17_02_117_fig_4.gif
(13)
q=kT
여기서 q는 단위시간에 단위면적을 통하는 열량(Heat flux), ∇T는 온도구배를 의미한다. 그리고v는 3.56 × 10-3 m2/s, ρ는 954 kg/m3이며, CP는 1,500 J/kg·K, k는 0.15W/m·K이다.

4.3 열 유동의 확산 분석

열에너지 방정식이 결합된 다상유동 기법을 검증하기 위해서 Garel et al.(2012)의 실험과 비교하였다. 유체가 평판을 이동함에 따라서 발생하는 열의 확산을 평가하였으며, Fig. 5는 본 연구의 수치해석에 의한 확산 결과이다. 이를 정량적으로 비교하기 위해서 Fig. 6과 같이 1,200초일 때 유체가 이동한 거리에 따른 온도의 분포를 실험과 비교하여 나타내었다. 본 연구 결과에서 온도의 감소가 다소 느리지만, 유체가 이동하면서 발생하는 열의 이동이 전반적으로 유사한 경향을 보이는 것을 확인할 수 있다. 따라서 본 연구에서 사용한 열에너지 방정식을 결합한 수치해석 기법의 타당성을 확인하였고 이를 용암류의 수치해석에 적용할 수 있다고 판단하였다.
Fig. 5
Thermal Diffusion Versus Time-varying
KOSHAM_17_02_117_fig_5.gif
Fig. 6
Surface Temperature Versus Distance to the Source
KOSHAM_17_02_117_fig_6.gif

5. 용암류 수치해석

5.1 개요

본 연구에서 용암류를 수치적으로 해결하기 위해서 Fig. 7과 같이 이탈리아의 시칠리아(Sicily, Italy) 섬에 위치하고 있는 에트나 화산(북위 37° 45´ 18´´, 동경 14° 59´ 42´´, 고도 3,350m)을 대상으로 정량적 평가를 수행하였다. 에트나 화산은 세계에서 가장 활발한 활화산 중 하나이며, 빈번하게 용암이 분출되고 있다. 특히, 현무암질 용암류가 분출되어 멀리 이동하며, 정상화구 이외에 지각의 틈에서도 분출되어 다량의 용암이 멀리 확산된 바 있다(Coltelli et al., 2007). 본 연구에서는 Fig. 8에 나타낸 바와 같이 2001년 분화가 발생한 에트나 화산에서 ‘LFS1’(Lower Fissure System 1)에서 분화지점을 대상으로 수치해석을 수행하였다. Coltelli et al.(2007)은 LFS1에서 분출한 용암을 매일 관측하였고 이로 부터 분출량, 이동경로, 확산 범위를 비롯하여 분화 전·후의 지형 고도를 측정하여 두께를 산정하는 등 상세한 분석이 이루어졌다. 이러한 관측 자료를 바탕으로 여러 수치모델을 활용하여 수치해의 유효성을 검토하는 연구가 수행된 바 있다(Cordonnier et al., 2015). 여기서 사용된 수치기법은 CA기법, 확률적 모델(Stochastic models) 혹은 2차원 영역으로 두께 산정과 수직분포의 영향을 고려하는 데 한계가 있다. 본 연구는 이를 보완하여 정량적인 수치를 산정하기 위해서 OpenFOAM에 의해 용암류를 모사하였다. 그리고 3차원 영역에서 수치기법의 유효성을 검증하기 위해서 현무암질 용암 특성이 잘 나타난 에트나 화산을 대상으로 선정하였고 화산 주변 지형에 대하여 90m 간격의 수치표고모형(Digital Elevation Model, DEM)을 이용해서 전산 해석에 적합한 형태로 적용하였다.
Fig. 7
Etna Volcano in Sicily, Italy(DEM) and Computational Domain(White line)
KOSHAM_17_02_117_fig_7.gif
Fig. 8
Computational Grid
KOSHAM_17_02_117_fig_8.gif

5.2 해석영역과 경계조건

본 연구의 전산해석 영역은 Fig. 7에 표기한바와 같이 대상 용암류와 인접한 영역으로 축소하였다. LFS1에서 분출된 용암의 이동경로를 중심으로 약 80여만 개의 육면체 격자를 생성하였다(Fig. 8). 수치해석에 사용된 경계조건은 대기 중에 노출된 해석영역을 고려하기 위해서 상단은 중력의 영향을 반영한 대기압 조건과 지표면은 점착조건을 설정하였다. 용암의 분출구는 시간에 따라서 체적유량(m3/s)이 유입되도록 Table 3과 같이 설정하였다.
Table 3
Flow Rate Inlet Boundary Condition(Coltelli et al., 2007)
Date  Time (s)   Cumulative volume (×106m3)   Daily Effusion Rate(m3/s) 
 07.18, 03:00  0 0.00 0.00
07.18, 13:00 36,000 0.37 10.28
07.19, 16:00 133,200 1.70 13.68
07.20, 13:00 208,800 3.50 23.81
07.22, 11:00 374,400 8.58 30.68
07.26, 12:00 723,600 14.98 18.33
07.28, 16:00 910,800 16.99 10.74
용암은 현무암질이 용융된 형태로서 그 물성은 유체의 점성, 밀도, 열용량 그리고 열전달계수와 비선형 점성 모델 Herschel-Bulkley의 매개변수를 Table 4와 같이 설정하였다. 여기서 설정한 매개변수는 초기조건에 해당하며, 이는 에트나 화산의 용암 시료와 경험식을 근거로 추정한 물성이다(Castruccio et al., 2010, 2014). 한편, 점성계수의 경우 비선형성과 온도의 영향으로 매 시간 갱신되어 적용된다.
Table 4
Properties of the Basaltic Lava Flow(Castruccio et al., 2010, 2014)
 Properties   Parameter  Unit Value
Lava Air
Fluid ν m2/s 1.51 1.51e-05
ρ kg/m3 2,700 1.205
CP  J/kg-K  1,200 1,005
κ W/m.K 1-2 0.0257
Viscous ν m2/s  9.63e+04   Newtonian μ = const
τ0 m2/s 0.07
κ m2/s2 1,163
n - 0.45

5.3 용암류 거동 분석

본 연구의 비교 대상으로 LFS1에서 분화한 용암류는 약 6,400m 이동하였으며, 최대 폭 545m, 평균 두께 11m로 관측되었다. 그러나 수치해석에 의한 이동거리는 약 3,400m로서 관측 값에 미치지 못하는 것으로 나타났다.
앞서 수행한 상대적으로 단순한 형태의 비선형 점성유동과 열 확산 등에서 수치기법의 타당함을 확인하였지만 해석영역이 확대되고 물리적 현상이 복잡한 용암을 반영하는 데 한계가 있는 것으로 판단된다. 이를 파악하기 위해서 Fig. 9에 나타낸 바와 같이 각 단면에서의 폭과 두께를 분석하였다(Figs. 10, 11). 수치해석에 의한 용암의 폭과 두께는 관측 값보다 과도하게 산정 되었다. 특히, 분출 초기에 이동성이 저하되고 측면으로 확산되거나 퇴적되는 경향을 보이고 있다.
Fig. 9
LFS1 Lava Flow Field(Coltelli et al., 2007)
KOSHAM_17_02_117_fig_9.gif
Fig. 10
Lava Width of Each Section
KOSHAM_17_02_117_fig_10.gif
Fig. 11
Lava Thickness of Each Section
KOSHAM_17_02_117_fig_11.gif
이러한 문제는 용암의 물리적 성질에 기인한 것으로 판단되며, 이는 용암이 약 1,000~1,300°C에서 용융된 암석으로 이루어진 특성과 지표면을 흐르면서 점차 온도가 낮아지고 점성이 높아짐에 따라서 이동성이 저하되어 오차가 발생한 것으로 판단된다. Fig. 12는 이동거리가 3,400m에 도달했을 때 용암의 종단면에 대한 온도-점도분포이며, 일반적으로 용암의 고상선(Solidus line)의 온도는 985~1,010°C지만
Fig. 12
Viscosity Versus Temperature at Lava and Ground Surface
KOSHAM_17_02_117_fig_12.gif
(Helz et al., 1987; Kauahikaua et al., 1998), 본 연구에서는 다소 높은 1,150°C정도에서 점도가 증가하여 이동성이 저하된 것으로 판단된다. 또한, 지표면에서의 점도보다 용암의 표면에서 점도가 높은 것으로 보아 공기와 맞닿은 면에서 먼저 경화가 발생할 것으로 사료된다. 이때 시간상 먼저 분출된 용암이 공기와 맞닿은 표면부터 경화되면서 단열작용을 하여 내부에서 흐르는 용융상태의 용암의 온도를 유지함으로써 유동성을 증가시키는 역할을 한다(Hon et al., 2008; Peterson et al., 1994). 그러나 본 연구에서는 이와 같은 경화-용융현상을 배제하였기 때문에 유동성에 대해 오차가 발생한 것으로 보아 향후 상전이 현상을 수치적으로 구현해야할 것으로 판단된다.
이처럼 현재 열에너지 방정식이 상전이(Phase change) 현상을 배제하고 있기 때문에 경화 및 용융현상이 이루어지지 않으며, 따라서 Eqs. (14)-(15)와 같이 외부 유동 조건과 온도의 영향을 반영하여 엔탈피(Enthalpy)가 고려된 에너지 방정식이 도입이 요구된다.
(14)
ρCP(Ht+UΔH)=ϕ+(κT)
(15)
H={CPTforT<TsSolidCPT+TTsTlTsLforTsTTlMushyCPT+LforT>TlLiquid
또한, 직교 격자를 사용함에도 불구하고 수치 해가 불안정한 원인은 수직방향의 격자가 다소 큰 것으로 판단된다. 현재 해석 영역에서 수직방향의 최소 격자는 약 6-7m이며, 이는 지표면의 최소 격자 크기가 7m인 것과 전 영역의 격자수를 고려하여 설정한 수치다. 그러나 LFS1에서 분출된 용암의 평균 두께는 11m로서 전산 격자가 약 1-2개 포함되어 해의 정확성이 현저히 감소한 것으로 판단된다. 이처럼 실제 현상과 비교해서 오차가 발생한 원인은 용융상태의 암석으로 이루어진 용암이 온도에 따라서 발생하는 상전이 현상에 대해서 고려되지 않은 것과 전산 격자에 기인한 것으로 판단된다.

6. 결론

본 연구에서는 다상유동 기법을 활용하여 비선형 점성유체를 대상으로 공개 소스 라이브러리인 OpenFOAM을 이용하여 3차원 전산유동해석을 수행하였다. 화산의 근접재해 중 용암류를 다루기 위해서 비교적 단순한 형태의 댐 붕괴 해석으로 비선형 점성유동의 이동거리 특성을 검증하였다. 댐 붕괴 실험과 비교하여 본 연구의 수치기법이 유체의 거동을 잘 모사하고 있음을 확인하였다. 그리고 온도의 영향을 고려하기 위해서 기존의 다상유동 기법과 열에너지 방정식을 결합하였으며, 유체의 이동에 따라서 열이 이동하는 양상을 실험과 비교하여 검증하였다. 이때 수치해석에서 온도의 감소는 다소 느리지만 전반적으로 유체가 이동하며 열이 이동하는 경향은 유사한 것으로 확인하였다. 따라서 이상의 수치기법을 이용해서 에트나 화산에서 발생한 용암류에 적용하였고 Herschel-Bulkley 유체의 매개변수와 열 물성을 통해서 수치해석을 수행하였다. 그러나 관측에 의한 이동거리와 수치해석과 비교했을 때 각각 6,400m와 3,400m로 큰 차이가 발생하였으며, 이러한 원인을 파악하기 위해서 용암의 각 단면에서 폭과 두께를 분석하였다. 여기서 분화 초기에 폭과 두께가 관측 값보다 과도하게 산정된 것으로 보아 용암의 온도가 감소하면서 점성이 높아짐에 따라서 이동성이 저하된 것으로 판단되었고 실제 용암의 온도가 변할 때 발생하는 경화현상과 용융현상 등 상전이 효과가 반영되지 않아서 발생한 오차로 판단되었다. 이를 해결하기 위해서 향후 상전이 영향이 고려된 에너지 방정식이 도입될 필요가 있을 것으로 판단된다.

감사의 글

본 연구는 정부(국민안전처)의 재원으로 자연재해저감기술개발사업단의 지원을 받아 수행된 연구임[MPSS-자연-2015-81].

References

Castruccio, A, Rust, A.C, and Sparks, R.S.J (2010) Rheology and flow of crystal-bearing lavas: Insights from analogue gravit currents. Earth and Planetary Science Letters, Vol. 297, pp. 471-480. 10.1016/j.epsl.2010.06.051.
crossref
Castruccio, A, Rust, A.C, and Sparks, R.S.J (2014) Assessing lava flow evolution from post-eruption field data using Herschel-Bulkley rheology. Journal of Volcanology and Geothermal Research, Vol. 275, pp. 71-84. 10.1016/j.jvolgeores.2014.02.004.
crossref
Coltelli, M, Proietti, C, Branaca, S, Marsella, M, Andronico, D, and Lodato, L (2007) Analysis of the 2001 lava flow eruption of Mt. Etna from three-dimensional mapping. Journal of Geophysical Research, Vol. 112, pp. 1-18. 10.1029/2006JF000598.
crossref
Cordonnier, B, Lev, E, and Garel, F (2015). Benchmarking lava-flow models. Geological Society. London: Special Publications, Vol. 426: p 1-21.
crossref
Crisci, G.M, Rongo, R, Gregorio, S.D, and Spataro, W (2004) The simulation model SCIARA: the 1991 and 2001 lava flows at Mount Etna. Journal of Volcanology and Geothermal Research, Vol. 132, pp. 253-267. 10.1016/S0377-0273(03)00349-4.
crossref
Garel, F, Kaminski, E, Tait, S, and Limare, A (2012) An experimental study of the surface thermal signatrue of hot subaerial isoviscous gravity currents: Implications for thermal monitoring of lava flows and domes. Journal of Geophysical Research, Vol. 117, pp. 1-18. 10.1029/2011JB008698.
crossref
Griffiths, R.W (2000) The Dynamics of Lava Flows. Annual Reviews of FLuid Mechanics, Vol. 32, pp. 477-518. 10.1146/annurev.fluid.32.1.477.
crossref
Kilburn, C.R.J (2000). Lava flows and flow fields. In: Sigurdsson H, ed. Encyclopedia of Volcanoes. Academic Press, San Diego: p 291-305.
crossref
Lev, E, Spiegelman, M, Wysocki, R.J, and Karson, J.A (2012) Investigating lava flow rheology using video analysis and numerical flow models. Journal of Volcanology and Geothermal Research, Vol. 247-248, pp. 62-73. 10.1016/j.jvolgeores.2012.08.002.
crossref
Liu, G, Yang, J, Wang, L, and Sun, J (2011). Analysis of Tianchi Volcano Activity in Changbai Mountain. NE China: Global Geology, Vol. 14: No. 1, p 44-53.
crossref
Minussi, R.B, and Maciel, G.F (2012) Numerical experimental comparison of dam break flows with non-Newtonian fluids. Journal of the Brazilian Society of Mechanical Sciences and Engineering, Vol. 34, No. 2, pp. 167-178. 10.1590/S1678-58782012000200008.
crossref
Saramito, P, Smutek, C, and Cordonnier, B (2013) Numerical modeling of shallow non-Newtonian flows: Part I. The 1D horizontal dam break problem revisited. International Journal of Numerical Analysis & Modeling, Series B, Vol. 4, No. 3, pp. 283-298.
crossref
Vicari, A, Alexis, H, Negro, C.D, Coltelli, M, Marsella, M, and Proietti, C (2007) Modeling of the 2001 lava flow at Etna volcano by a Cellular Automata approach. Environmental Modelling & Software, Vol. 22, pp. 1465-1471. 10.1016/j.envsoft.2006.10.005.
crossref
Yun, S.H, and Lee, J.H (2012) Analysis of Unrest Signs of Activity at the Baegdusan Volcano. The Petrological Society of Korea, Vol. 21, No. 1, pp. 1-12. 10.7854/JPSK.2012.21.1.001.
crossref


ABOUT
ARTICLE CATEGORY

Browse all articles >

BROWSE ARTICLES
AUTHOR INFORMATION
Editorial Office
1010 New Bldg., The Korea Science Technology Center, 22 Teheran-ro 7-gil(635-4 Yeoksam-dong), Gangnam-gu, Seoul 06130, Korea
Tel: +82-2-567-6311    Fax: +82-2-567-6313    E-mail: master@kosham.or.kr                

Copyright © 2024 by The Korean Society of Hazard Mitigation.

Developed in M2PI

Close layer
prev next