Performance Evaluation of OGS-FLAC Simulator for Coupled Thermal-Hydrological-Mechanical Analysis

Dohyun Park; Chan-Hee Park

doi:10.7474/TUS.2022.32.2.144

Preview

Original Article

Tunnel and Underground Space. 30 April 2022. 144-159
https://doi.org/10.7474/TUS.2022.32.2.144

Performance Evaluation of OGS-FLAC Simulator for Coupled Thermal-Hydrological-Mechanical Analysis

열-수리-역학적 연계해석을 위한 OGS-FLAC 시뮬레이터의 성능 평가

Dohyun Park¹

Chan-Hee Park¹^*

박 도현¹

박 찬희¹^*

¹Principal Researcher, Deep Subsurface Storage and Disposal Research Center, Korea Institute of Geoscience and Mineral Resources

¹한국지질자원연구원 심층처분환경연구센터 책임연구원

^{*Corresponding Author}

ABSTRACT

The present study developed a sequential approach-based numerical simulator for modeling coupled thermal-hydrological-mechanical (THM) processes in the ground and investigated the computational performance of the coupling analysis algorithm. The present sequential approach linked the two different solvers: an open-source numerical code, OpenGeoSys for solving the thermal and hydrological processes in porous media and a commercial code, FLAC3D for solving the geomechanical response of the ground. A benchmark test of the developed simulator was carried out using a THM problem where an analytical solution is given. The benchmark problem involves the coupled behavior (variations in temperature, pore pressure, stress, and deformation with time) of a fully saturated porous medium which is subject to a point heat source. The results of the analytical solution and numerical simulation were compared and the validity of the numerical simulator was investigated.

Keywords

Coupled thermal-hydrological-mechanical behavior

Sequential coupling analysis

Open-source simulator

Benchmark test

OGS-FLAC

본 연구에서는 지반의 열-수리-역학적 복합거동을 모델링하기 위한 순차적 접근법 기반의 시뮬레이터를 개발하고 적용된 연계해석 알고리즘의 계산성능을 분석하였다. 본 연구의 순차적 연계해석에서는 다공성 매질의 열 및 유체거동 분석을 위한 오픈소스 기반의 OpenGeoSys 수치코드와 역학해석을 위한 상용 소프트웨어 FLAC3D가 연동되었다. 해석해가 주어진 열-수리-역학적 복합거동 문제를 토대로 개발된 시뮬레이터에 대한 벤치마크 테스트가 수행되었다. 적용된 벤치마크 문제는 완전포화된 지반 내 점열원 작용 시 지반거동(시간에 따른 온도, 간극수압, 응력, 변형 변화)과 관계된다. 해석해와 수치해석 시뮬레이션 결과를 비교 분석하고 연계해석 시뮬레이터의 적정성을 조사하였다.

키워드

열-수리-역학적 복합거동

순차적 연계해석

오픈소스 시뮬레이터

벤치마크 테스트

OGS-FLAC

국문키워드2

국문키워드3

국문키워드4

MAIN

1. 서 론
2. 벤치마크 문제
3. OGS-FLAC 기반 순차적 연계해석의 적정성 분석
3.1 해석모델 및 조건
3.2 연계해석 방식
3.3 해석 결과
4. 결 론

1. 서 론

지표 하에서 수행되는 에너지 개발 및 활용, 폐기물 처분과 관련된 활동들은 지반 자체의 지압, 지열, 지하수 그리고 에너지 및 폐기물 처분과 관련된 열원, 유체압력, 화학작용 등으로 인해 열-수리-역학-화학적 반응이 서로 연계된 다중물리(multiphysics) 복합거동이 발생할 수 있다. 다중물리 복합거동이 수반되는 에너지 개발 및 활용 관련 활동들에는 가스하이드레이트 함유 퇴적층으로부터 자원 생산, 셰일가스와 같은 비전통 탄화수소 자원의 개발, 심부 지반의 지열 회수, 대수층 및 시추공 등을 활용한 지하 에너지 저장, 에너지 말뚝(energy pile) 및 토목섬유로 보강된 옹벽과 같은 지반구조물을 활용한 열에너지 저장 등이 있다. 지반을 대상으로 한 폐기물 처분에는 이산화탄소 및 고준위방사성폐기물의 지층 처분이 포함되며, 이들 분야는 환경 및 안전 이슈와 관련되어 있어 신뢰도 높은 장기적인 복합거동 분석기술이 필요하다.

지반 내 여러 가지 물리현상의 복합적인 상호작용이 매질의 유체흐름, 열전달, 역학적 변형과 같은 지반거동에 미치는 영향을 분석하기 위해 실험적 접근법과 해석적 접근법이 활용된다. 해석적 접근법 중 수치해석의 경우 현장의 다양한 기하학적 조건 및 복잡한 지층 조건을 모델링할 수 있는 장점이 있다. 열-수리-역학적 복합거동 분석을 위한 수치해석적 접근법은 CODE_BRIGHT (UPC, 2021), FRACON (Nguyen and Selvadurai, 1995)과 같은 일체식 연계해석법(monolithic coupling method)과 순차적(sequential or staggered) 연계해석법(Rutqvist et al., 2002, Lim et al., 2007, Taron et al., 2009, Graupner et al., 2011, Kim et al., 2011, Kim et al., 2012, Kim, 2015, Elyasi et al., 2016, Yoon et al., 2020, Kim et al., 2021)으로 크게 분류할 수 있다. 일체식 연계해석법은 다중물리 거동을 하나의 계(system)로 간주하여 계산하는 형태로 서로 다른 물리영역 간의 상호작용을 한번에 고려하며 이를 위해 생성된 미분방정식들이 단일 행렬로 구성된다. 순차식 연계해석법은 해석대상의 물리영역을 하나씩 순차적으로 풀고 사전에 정해진 수치오차의 허용치를 만족할 때까지 반복하여 계산하는 방식이다. 일반적으로 계의 변화를 유발하는 물리현상을 먼저 풀고 그 결과를 다른 물리영역의 초기값으로 넘겨주며 반복 계산이 진행된다. 기 개발된 수치코드들을 활용하는 경우, 코드 커플링 인터페이스 작성이 주된 작업이므로 일체식 연계해석법에 비해 시뮬레이터 개발시간을 줄일 수 있다.

지반의 다중물리 복합거동 해석기술 연구는 국내외적으로 활발히 진행 중이며, 국내의 경우 원천적인 기술은 해외에 의존하는 경향이 크므로 독자적인 연계해석 시뮬레이터 개발이 필요한 상황이다. 해외 의존도를 줄이기 위한 시뮬레이터 개발의 일환으로 Park et al.(2019)은 오픈소스 수치코드인 OpenGeoSys (OGS; OpenGeoSys Community, 2021)와 상용 프로그램인 FLAC3D (Itasca, 2021)를 연동시킨 OGS-FLAC 시뮬레이터 초기버전을 개발하였다. 이 연구의 시뮬레이터는 지반의 온도영향을 고려하지 않은 수리-역학적 연계해석 기반의 알고리즘으로 설계되어 시뮬레이터의 활용도가 제한적일 수 있다. 본 연구에서는 지반의 열적 거동을 고려한 열-수리-역학적 복합거동 시뮬레이션이 가능하도록 OGS-FLAC 해석기능을 업그레이드하고, 벤치마크 테스트를 통해 시뮬레이터의 계산성능을 조사하였다. 유체거동 지배방정식을 구성하는 여러 항(term)들 중 역학적 부분(지반의 체적변화로 인한 유체의 질량변화, 즉 지반의 역학적 변형에 따른 유체의 유입 및 유출)을 반영하여 유체영역의 지배방정식이 완전한 형태로 구현되도록 연계해석 알고리즘이 작성되었으며, 해석해가 주어진 열-수리-역학적 벤치마크 문제를 토대로 작성된 알고리즘의 적정성을 분석하였다.

2. 벤치마크 문제

점열원(point heat source)이 작용하는 지반의 열-수리-역학적 복합거동 문제에 대한 해석해를 제안한 Booker and Savvidou(1985) 연구를 토대로 본 연구의 연계해석 알고리즘의 적정성을 분석하고자 하였다. Booker and Savvidou(1985)는 점열원이 작용하는 주변지반의 시간경과에 따른 지반의 온도, 간극수압, 변위, 응력 변화에 대한 해석적 해를 제안하였다. 대상지반은 물로 완전히 포화되고 역학적으로 탄성 거동을 보이는 것으로 가정하였다. 이 문제에서는 중력으로 인한 영향은 고려하지 않았다. 지반 내 열원은 매질의 온도를 상승시키고 상승된 온도로 인해 간극수와 지반 매질을 팽창시킨다. 간극수의 부피증가량이 매질 내 간극보다 큰 경우에는 간극수압의 증가로 이어지며 이로 인해 매질의 유효응력이 감소하게 되는 열-수리-역학적 복합거동을 보이게 된다. 본 벤치마크 문제에서는 열원에 의한 재료의 열적 거동변화는 유체의 수리적 거동 및 지반의 역학적 거동에 영향을 미치나, 반대로 수리적 및 역학적 거동 변화는 열적 거동에 영향을 미치지 않는다고 가정하였다. 유체의 수리 거동과 지반의 역학 거동은 상호 영향을 받는 것으로 고려하였다. 해석해를 구하기 위한 열, 수리, 역학 분야의 관계식들을 정리하면 다음과 같다. 본 논문에서 언급된 수식들의 매개변수 단위는 사용자가 설정한 단위체계에 따라 달라질 수 있으며, 본 연구의 단위체계는 SI 단위를 기반으로 하고 다음과 같다. 길이 m, 시간 sec, 힘 N, 온도 K, 밀도 kg/m³, 응력 Pa, 압력 Pa, 투수계수 m², 열전도도 W/m-K, 체적 열팽창계수 1/K, 점성도 Paㆍsec, 열용량 J/kg-K, 열원 W이다.

❑ 매질의 열에너지 균형방정식(thermal energy balance equation)

(1)

m \frac{\partial θ}{\partial t} + \nabla^{T} h = q

(2)

m = n ρ_{w} c_{w} + (1 - n) ρ_{s} c_{s}

(3)

h = - K \nabla θ

여기서 m은 체적 열용량(volumetric heat capacity), θ는 온도, t는 시간, h는 열유속(heat flux) 벡터, q는 열원 세기(intensity), n은 공극률, ρ_w와 ρ_s는 각각 물과 지반의 밀도, c_w와 c_s는 각각 물과 지반의 열용량, K는 지반의 열전도도이다.

❑ 유체 질량 균형방정식(mass balance equation)

(4)

\int_{0}^{t} \nabla^{T} v d t = a_{u} θ - ε_{v}

(5)

v = - \frac{k}{γ_{w}} \nabla p

(6)

a_{u} = n a_{w} + (1 - n) a_{s}

여기서 v는 단위시간당 유출량(flux) 벡터, ε_v는 체적 변형률, k는 지반의 투수계수, γ_w는 물의 단위중량, p는 유체압력, a_u는 매질 전체의 체적 열팽창계수(overall volumetric thermal expansion coefficient), a_w와 a_s는 각각 물과 지반의 체적 열팽창계수, n은 공극률이다.

❑ 지반의 운동량 균형방정식(momentum balance equation)

(7)

\partial^{T} σ = 0

(8)

σ_{i j | i = j} = p + b' θ + λ ε_{v} + 2 G ε_{i j}

(9)

σ_{i j | i \neq j} = G γ_{i j}

(10)

b' = (λ + \frac{2 G}{3}) a'

여기서 σ는 유효응력, ε은 변형률, ε_v는 체적 변형률, λ는 Lamé 계수, G는 전단계수, γ는 전단변형률, a′은 지반의 체적 열팽창계수이다.

상기의 열, 수리, 역학 분야의 방정식들은 지반 내 열전달은 열전도(heat conduction)에 의해서만 발생하고, 작용 응력으로 인한 지반입자의 체적변화와 간극수압으로 인한 간극수의 체적변화는 없는 것으로 가정하였다. 지반입자와 물의 비압축성 가정에 의해 유체 유출량은 응력으로 인한 지반의 체적변화와 열에 의한 체적변화에 의해서만 산출된다. 유체 질량 균형방정식과 운동량 균형방정식을 단순화하여 정리하면 다음과 같다.

(11)

\frac{\partial ε_{v}}{\partial t} - a_{u} \frac{\partial θ}{\partial t} + \nabla^{T} v = 0

(12)

κ \nabla^{2} θ = \frac{\partial θ}{\partial t} - \frac{q}{m}

(13)

κ = \frac{K}{m}

(14)

G \nabla^{2} u_{x} - (λ + G) \frac{\partial ε_{v}}{\partial x} = \frac{\partial p}{\partial x} + b' \frac{\partial θ}{\partial x}

(15)

G \nabla^{2} u_{y} - (λ + G) \frac{\partial ε_{v}}{\partial y} = \frac{\partial p}{\partial y} + b' \frac{\partial θ}{\partial y}

(16)

G \nabla^{2} u_{z} - (λ + G) \frac{\partial ε_{v}}{\partial z} = \frac{\partial p}{\partial z} + b' \frac{\partial θ}{\partial z}

여기서 u_x, u_y, u_z는 각각 x, y, z방향 변위이다.

위에서 언급한 지배방정식들에 대한 해를 구하기 위해 푸리에(Fourier) 변환과 라플라스(Laplace) 변환과 같은 적분변환(integral transform)이 적용되었으며 최종 해석해는 다음과 같이 제시되었다.

(17)

θ = \frac{Q}{4 π K R} f (\frac{κ t}{R^{2}})

(18)

p = \frac{X}{1 - c / κ} \frac{Q}{4 π K R} [f (\frac{κ t}{R^{2}}) - f (\frac{c t}{R^{2}})]

(19)

u_{x_{i}} = a_{u} x_{i} \frac{Q}{4 π K R} g^{*}

(20)

σ_{x_{i} x_{i}} = 2 G a_{u} \frac{Q}{4 π K R} [f^{*} - g^{*} + \frac{x_{i}^{2}}{R^{2}} (3 f^{*} - g^{*})]

(21)

σ_{x_{i} x_{j} | i \neq j} = 2 G a_{u} \frac{Q}{4 π K R} \frac{x_{i} x_{j}}{R^{2}} (3 f^{*} - g^{*})

여기서 Q는 열원, x_i는 열원으로부터 이격된 위치, i는 x, y, z축을 나타낸다. 상기 해석해에서 사용된 변수들의 정의는 다음과 같다.

(22)

κ = K / m

(23)

c = k (λ + 2 G) / γ_{w}

(24)

R = \sqrt{x^{2} + y^{2} + z^{2}}

(25)

X = a_{u} (λ + 2 G) - b'

(26)

Y = \frac{1}{λ + 2 G} [\frac{X}{(1 - c / κ) a_{u}} + \frac{b'}{a_{u}}]

(27)

Z = \frac{1}{λ + 2 G} [\frac{X}{(1 - c / κ) a_{u}}]

(28)

f (\frac{κ t}{R^{2}}) = e r f c (\frac{R}{2 \sqrt{κ t}})

(29)

g (\frac{κ t}{R^{2}}) = \frac{κ t}{R^{2}} + (\frac{1}{2} - \frac{κ t}{R^{2}}) e r f c (\frac{R}{2 \sqrt{κ t}}) - \sqrt{\frac{κ t}{π R^{2}}} e^{- R^{2} / 4 κ t}

(30)

f^{*} = Y f (\frac{κ t}{R^{2}}) - Z f (\frac{c t}{R^{2}})

(31)

g^{*} = Y g (\frac{κ t}{R^{2}}) - Z g (\frac{c t}{R^{2}})

발생응력에 대한 해석해는 Chaudhry et al.(2019)에 의해 보완되었으며 다음과 같이 제안되었다.

(32)

σ_{x_{i} x_{i}} = \frac{Q a_{u}}{4 π K R} \{2 G [g^{*} (1 - \frac{x_{i}^{2}}{R^{2}}) + x_{i} g_{i}^{*}] + λ [x_{i} g_{i}^{*} + 2 g^{*}]\} - b' Δ θ

(33)

σ_{x_{i} x_{i} i \neq j} = \frac{Q a_{u}}{4 π K R} \{G (x_{i} g_{i}^{*} + x_{j} g_{i}^{*}) - 2 g^{*} \frac{x_{i} x_{j}}{R^{2}}\}

상기와 같이 본 벤치마크 문제에 대한 Booker and Savvidou(1985)와 Chaudhry et al.(2019)의 해석해는 응력부분에서 차이가 있다. 열-수리-역학적 연계해석 수치코드의 검증을 위해 본 벤치마크 문제를 활용한 Chaudhry et al.(2019) 이전의 과거 연구사례로는 Nguyen and Selvadurai(1995), Blanco-Martín et al.(2017)이 있다. Nguyen and Selvadurai(1995)는 일체식 연계해석 코드인 FRACON을 검증했고, Blanco-Martín et al.(2017)은 순차식 연계해석에 기반한 TOUGH-FLAC 시뮬레이터(Rutqvist et al., 2002)의 성능을 분석하였다. 이 연구사례들에서는 지반온도 및 간극수압에 대해서는 해석해와 시뮬레이션 결과를 비교한 반면, 역학적 물리값(변위 및 응력)에 대해서는 검증결과를 제시하지 않은 것으로 파악되었다. Fig. 1은 두 연구사례에서 사용된 수치해석 모델들을 나타내고, Fig. 2는 지반온도 및 간극수압에 대해 해석해와 시뮬레이션 결과를 비교한 그림으로 역학적 결과는 포함되지 않았다. Chaudhry et al.(2019) 이전에 수행된 상기 선행연구들에서 역학적 항목에 대한 검증결과를 제시하지 않은 것은 Booker and Savvidou(1985)의 응력 해석해와 시뮬레이션 결과의 차이 발생과 관련되었을 것으로 판단된다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-02/N0120320205/images/ksrm_32_02_05_F1.jpg

Fig. 1.

Schematic representation of Booker and Savvidou’s problem in (a) Nguyen and Selvadurai (1995) and (b) Blanco-Martín et al. (2017)

https://static.apub.kr/journalsite/sites/ksrm/2022-032-02/N0120320205/images/ksrm_32_02_05_F2.jpg

Fig. 2.

Comparison between analytical and modeling results for temperature and pore pressure in (a) Nguyen and Selvadurai (1995) and (b) Blanco-Martín et al. (2017)

3. OGS-FLAC 기반 순차적 연계해석의 적정성 분석

3.1 해석모델 및 조건

Booker and Savvidou(1985)가 제안한 해석해는 완전포화 지반에 점열원 작용 시 열원으로부터 일정거리 이격된 위치에서 시간경과에 따른 지반의 온도, 간극수압, 변위, 응력변화 정보를 제공한다. 본 연구에서 설계한 OGS-FLAC 연계해석 알고리즘의 적정성을 평가하기 위해 Booker and Savvidou(1985)의 해석해와 OGS-FLAC 시뮬레이션 결과를 비교 분석하였다. Fig. 3은 수치해석 모델 및 해석 그리드를 나타낸다. 좌표 (0, 0, 0)에서 열원(Q = 300 W)이 작용하고 모델영역은 30 × 30 × 30 m로 설정하였다. 수치해석 시간을 단축하기 위해 해석요소(셀) 개수를 줄일 수 있는 1/8 대칭모델을 적용하였다. 지반거동 계측위치는 임의로 2개 지점(P1, P2)이 설정되었고 열원으로부터 이격된 거리는 각각 0.88, 1.75 m이었다. 계측지점들은 해석요소의 중심점이 되도록 설정하였으며 좌표는 P1 = (0.3346, 0.8054, 0.109), P2 = (0.2226, 1.5096, 0.853) 이었다. 해석해 조건과 동일하게 중력가속도는 고려되지 않았다. 해석해의 적용조건을 반영하여 간극 내 유체는 물, 매질은 물로 완전히 포화된 다공성 탄성지반으로 모델링되었다. 해석해 조건과 동일하게 초기 지반온도는 273.15 K, 초기 간극수압과 지반응력은 0(zero)으로 설정되었다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-02/N0120320205/images/ksrm_32_02_05_F3.jpg

Fig. 3.

Numerical model domain and the model grid

앞서 언급한대로 Booker and Savvidou(1985)에서는 응력으로 인한 지반입자의 체적변화와 간극수압으로 인한 간극수의 체적변화는 없는 것으로 가정하였다. 따라서 아래 식 (34)의 일반화된 유체 질량 균형방정식에서 시간경과에 따른 유체압력 변화를 나타내는 좌측 첫 번째 항이 생략되어 앞서 기술한 식 (11)과 같이 적용되었다. 그러나 서로 다른 수치코드들을 연계하는 OGS-FLAC에서는 식 (34)의 첫 번째 항이 고려되지 않으면, 시간경과에 따른 유체압력 변화를 구할 수 없는 조건이 되어 이 항을 생략하기 어렵다. 즉, OGS-FLAC에서는 매질 전체의 압축률(overall compressibility)을 고려해서 수치해석을 수행해야 한다. 이는 유체영역에서는 해석해의 지반입자 및 물의 비압축성 가정과는 다른 조건으로 OGS-FLAC 연계해석이 이루어짐을 의미한다.

(34)

\frac{1}{M} \frac{\partial p}{\partial t} + α \frac{\partial ε_{v}}{\partial t} - a_{u} \frac{\partial θ}{\partial t} + \nabla^{T} v = 0

(35)

\frac{1}{M} = C = \frac{n}{K_{f}} + \frac{(α - n)}{K_{s}} = n C_{f} + (α - n) C_{s}

(36)

n = n_{0} + α (ε_{v} - ε_{v 0})

여기서 M은 Biot’s modulus, C는 매질 전체의 압축률(또는 저류계수(storativity)), p는 유체압력, α는 Biot coefficient(=1), ε_v는 체적 변형률, θ는 온도, v는 단위시간당 유출량(flux), n은 공극률, K_f와 K_S는 각각 유체와 고체입자의 체적 변형계수, C_f와 C_S는 각각 유체와 고체입자의 압축률이다.

Booker and Savvidou(1985)의 해석해와 OGS-FLAC 연계해석의 재료 압축성에 대한 기본가정의 상이함으로 인한 결과 차이를 분석하기 위하여 두 가지 해석 case (Case1, Case2)를 고려하였다. Case1은 재료의 압축성이 작은 지반조건으로 설정되었고, Case2 해석에서는 Case1보다 상대적으로 압축성이 큰 지반조건을 설정하였다. Table 1은 Case1과 Case2 해석에 적용된 재료물성을 나타낸다. Case2의 경우 지반의 영률은 Case1의 1/10 수준으로 Case2의 압축성이 Case1에 비해 상대적으로 큰 조건이다.

Table 1.

Parameters and values used in the present analysis

Material	Parameter	Unit	Value
Ground	Density	kg/m³	2,700
(Case1)	Young’s Modulus	MPa	50,000
	Poisson’s ratio	-	0.2
	Thermal conductivity	W/m-K	1.83
	Specific heat capacity	J/kg-K	780
	Vol. thermal expansion coefficient	1/K	2.0×10^-5
	Porosity	-	0.015
	Permeability	m²	2.0×10^-19
	Biot coefficient	-	1
Water	Density	kg/m³	1000
	Viscosity	Paㆍsec	0.001
	Specific heat capacity	J/kg-K	4,200
	Thermal conductivity	W/m-K	0.6
	Vol. thermal expansion coefficient	1/K	2.0×10^-4
Ground	Density	kg/m³	2,290
(Case2)	Young’s Modulus	MPa	5,000
	Poisson’s ratio	-	0.3
	Thermal conductivity	W/m-K	1.83
	Specific heat capacity	J/kg-K	920
	Vol. thermal expansion coefficient	1/K	1.5×10^-5
	Porosity	-	0.16
	Permeability	m²	2.0×10^-20
	Biot coefficient	-	1

3.2 연계해석 방식

OGS-FLAC 연계해석을 위해 OGS에서 FLAC3D를 제어하는 방식을 적용하였다. 즉, OGS가 주(master)로 동작하면서 FLAC3D를 종(slave)으로 작동시키는 방식이다. Fig. 4는 적용된 순차적 연계해석 방식을 나타내며, 기존의 순차적 접근법 기반의 시뮬레이터와 유사한 형태로 구성되었다. 유체거동 지배방정식에 포함되는 역학적 변화로 인한 영향(식 (34)의 좌측 두 번째 항)은 초기버전의 OGS-FLAC에서는 고려하지 못했던 부분이다. 본 연구에서는 이 부분이 모사될 수 있도록 알고리즘을 보완하였으며, 이를 통해 유체영역의 지배방정식이 완전한 형태로 구현가능한 것이 초기버전과의 주요한 차이점이다. 식 (34)의 두 번째 항은 단위시간당 매질의 역학적 변형으로 인한 유체의 질량변화를 나타내며, 수리해석 시 FLAC3D의 체적변형률 자료를 전달받아 소스텀(source term)으로 모델링하여 유체거동이 모사된다. 본 연구의 순차적 연계해석에서는 열 및 유체유동이 먼저 계산되고 이후 역학해석이 진행된다. OGS에서 열-수리 해석을 수행하여 매질 내 물리거동 자료(간극수압, 온도 등)를 FLAC3D에 전달하고, 이 자료를 토대로 FLAC3D에서는 역학해석(유효응력 변화 고려)을 실시한다. FLAC3D에서 분석된 역학해석 자료는 OGS로 전달되어 매질의 체적변형률, 공극률 등의 변화를 반영하여 열-수리 해석을 재수행한다. 이러한 각 솔버의 독립적인 계산과 솔버 간 자료교환 과정은 수렴기준을 만족할 때까지 반복적으로 진행된다. 두 개 수치코드 간의 계산자료 교환은 해석요소의 중심을 기반으로 한다. OGS에서는 물리량이 해석절점 기준으로 계산되므로 요소를 구성하는 절점들의 결과값을 내삽하여 해당요소를 대표하는 물리량이 산출된다. 절점기준으로 계산된 결과가 요소중심으로 변환되면서 해당 물리량이 평균효과(averaging)에 의해 과소 평가된 상태로 역학해석이 수행될 수 있으므로 주의가 필요할 것으로 판단된다. 평균효과를 줄이기 위해 해석요소의 크기를 최대한 작게 설정하는 것이 대안이 될 수 있으나, 3차원해석을 수행하는 경우 시뮬레이션 시간과 계산 정확도와의 절충점을 찾을 필요가 있을 것이다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-02/N0120320205/images/ksrm_32_02_05_F4.jpg

Fig. 4.

Sequential coupling scheme used in the present study

3.3 해석 결과

두 가지 해석 case에 대한 OGS-FLAC 시뮬레이션 결과는 Figs. 5~8와 같다. 이 그림들은 계측지점 P1, P2에서 시간경과에 따른 지반의 온도, 간극수압, 변위, 응력 변화를 나타낸다. 열원에 의한 시간경과에 따른 지반의 온도 상승, 간극수압의 발생 및 지반의 투수성으로 인한 간극수압 소산(dissipation), 지반 변위 및 응력의 변화가 모사되었다. 앞서 언급한대로 Case1은 지반의 압축성이 작은 조건이고, Case2는 Case1에 비해 상대적으로 압축성이 큰 지반조건이다.

지반의 압축성이 작은 Case1 해석결과를 보면, OGS-FLAC 계산결과가 해석해와 정확하게 일치하지 않지만 해석해에 상당히 근접하는 것을 알 수 있다(Figs. 5~6). 시간경과에 따른 지반의 온도는 해석해와 차이가 크지 않고 다소 작게 예측되는 경향을 보였으며, 결과적으로 간극수압과 지반변위 발생이 작게 평가된 것으로 판단된다. 응력부분에서는 σ_yy와 σ_xy는 해석해와 잘 일치하는 결과를 보인 반면, σ_xx는 전체적인 변화 경향은 유사하나 해석해보다 작게 예측되는 결과를 보였다. 이러한 해석해와 시뮬레이션 결과의 차이는 앞서 언급한 재료의 압축성 가정조건의 차이, 이종의 수치코드들 간의 자료교환 방식, 수치해석 결과값 추출방법 등의 영향을 받았을 것으로 판단된다. Fig. 9는 계산시간 10⁵ sec에서 온도, 간극수압, 변위, 응력분포를 나타내고, 열원 주변에서 지반거동 변화가 집중됨을 알 수 있다.

지반의 압축성이 상대적으로 큰 Case2 해석결과를 보면, 지반 온도를 제외하고는 Case1에 비해 해석해와의 차이가 비교적 크게 발생하는 것을 알 수 있다(Figs. 7~8). 앞서 해석해 부분에서 기술한대로 Booker and Savvidou(1985) 연구에서는 열원에 의한 재료의 열적 거동변화는 유체의 수리적 거동과 지반의 역학적 거동에 영향을 미치나, 반대로 수리적 및 역학적 거동변화가 열적 거동에는 영향을 미치지 않는다고 가정하였다. 따라서 열해석은 수리-역학적 해석과는 별개로 독립적으로 수행되고, 재료의 압축성 전제조건과는 관계가 없으므로 해석해와 잘 일치하는 결과를 보였다. 간극수압의 경우, 해석해에서는 재료의 압축성이 무시되므로 식 (34)의 좌측 첫 번째 항이 고려되지 않는 반면, Case2에서는 압축률이 해석해보다 큰 상태로 계산이 수행되어 유체압력에 차이가 발생한 것으로 판단된다. 수치해석의 간극수압과 해석해의 편차로 인해 역학적 거동(변위, 응력발생)에 차이가 발생하고, 결과적으로 시간경과에 따라 수리-역학적 연계거동이 반복되면서 해석해에 근접하지 못하는 것으로 추정된다. σ_xx의 경우, peak 발생시점이 해석해와 다소 차이가 있는데 이는 수치코드 간 자료교환 방식, 해석결과 추출방법, 요소형상 등이 복합적으로 관련되어 있을 것으로 판단되며, 추가 연구를 통해 차이발생 원인을 파악할 필요가 있다.

상기의 분석결과로부터 Case1의 경우 OGS-FLAC 연계해석에 의해 해석해가 비교적 잘 예측된 반면, Case2에 대해서는 해석해와 차이가 비교적 크게 발생하는 것을 알 수 있었다. 해석해와 차이가 있는 Case2 결과로 인해 OGS-FLAC의 계산정확도가 떨어지는 것으로 평가하기보다는 해석해의 기본가정과 상이한 조건에서 연계해석이 수행되어 나타난 결과로 판단하는 것이 합리적일 것이다. 즉, Case1에서는 지반의 압축성이 상당히 작아 해석해의 비압축성 가정조건과 유사하여 해석해로부터 편차가 작게 발생하였으나, Case2의 지반조건은 Case1에 비해 지반의 압축성이 상대적으로 커 해석해의 비압축성 조건과 상이하여 계산결과에 차이가 발생한 것으로 판단된다. 두 가지 해석 case에 대한 분석결과로부터 본 연구에서 작성한 열-수리-역학적 연계해석 알고리즘의 적정성이 기본적으로 확보된 것으로 판단할 수 있다. 실제 현장의 유체와 지반입자는 완전한 비압축성 조건이 아니므로 향후 현장의 열-수리-역학적 복합거동 자료에 대한 사례연구를 통해 본 연구에서 설계한 OGS-FLAC 연계해석 알고리즘의 적용성을 추가적으로 분석할 필요가 있다.

https://static.apub.kr/journalsite/sites/ksrm/2022-032-02/N0120320205/images/ksrm_32_02_05_F5.jpg

Fig. 5.

Coupled THM responses with time at point P1 (Case1)

https://static.apub.kr/journalsite/sites/ksrm/2022-032-02/N0120320205/images/ksrm_32_02_05_F6.jpg

Fig. 6.

Coupled THM responses with time at point P2 (Case1)

https://static.apub.kr/journalsite/sites/ksrm/2022-032-02/N0120320205/images/ksrm_32_02_05_F7.jpg

Fig. 7.

Coupled THM responses with time at point P1 (Case2)

https://static.apub.kr/journalsite/sites/ksrm/2022-032-02/N0120320205/images/ksrm_32_02_05_F8.jpg

Fig. 8.

Coupled THM responses with time at point P2 (Case2)

https://static.apub.kr/journalsite/sites/ksrm/2022-032-02/N0120320205/images/ksrm_32_02_05_F9.jpg

Fig. 9.

Distributions of the ground temperature, pore pressure, displacement, and stress at 10⁵ sec (Case1)

4. 결 론

본 연구에서는 오픈소스 기반의 열-수리-역학적 연계해석 시뮬레이터를 개발하고, 해석해가 주어진 벤치마크 문제를 토대로 시뮬레이터의 계산성능을 조사하였다. 분석 결과, 해석해의 적용조건과 유사한 case의 경우 해석해에 근접하는 지반거동 결과를 제공하는 것으로 검토되어 개발된 연계해석 시뮬레이터의 적정성을 확인할 수 있었다. 벤치마크 문제에서 가정한 유체 및 고체의 비압축성 조건은 실제 자연현상과는 차이가 있으므로 실내 및 현장실험 자료에 기반하여 연계해석 알고리즘의 적용성을 추가 분석할 필요가 있을 것으로 판단된다. 한편, 본 연구에서 검토한 벤치마크 문제는 완전포화 상태의 단상 유체거동만 다루고 있으므로 향후 다상 다성분 유체유동 관련 열-수리-역학적 복합거동에 대한 현장자료 기반의 검증을 통해 개발 시뮬레이터의 기술 활용도 및 확장성을 높일 계획이다.

현재까지 연구자에 따라 다양한 접근법의 다중물리 복합거동 해석법들이 제시되어 왔으나, 대상문제별로 연계거동 특성이 다양하고 이를 모델링하기 위한 수식, 해석기법이 정형화되지 않은 상황이다. 오픈소스 기반의 시뮬레이터는 개발자 및 사용자 커뮤니티를 통해 문제점과 해결방법을 공유할 수 있고, 필요 시 해석기능 추가 및 보완 등 업데이트를 신속하게 할 수 있는 장점이 있다. 따라서 복잡한 형태의 물리 방정식 및 모델링 조건을 요구하는 연계해석 특성에 대응하기 위한 툴로서 오픈소스 기반의 시뮬레이터 적용을 고려할 필요가 있을 것으로 판단된다. 다만, 오픈소스 소프트웨어는 사용자의 요구사항에 부합되지 않은 기술개발, 갑작스런 개발 중단이 발생할 수 있으므로 사용자는 해당 소프트웨어에 대해 개발자 수준의 코딩 역량을 갖출 필요가 있다. 지반의 복합거동 연계해석 기술개발은 해외 선진국 위주로 진행되고 있는 상황으로 국내에서도 독자적인 기술력 확보 및 기술 선도를 위한 연구개발이 적극적으로 수행될 필요가 있을 것으로 판단된다.

Acknowledgements

본 연구는 한국지질자원연구원의 기본사업인 ‘심지층 개발과 활용을 위한 지하심부 특성평가 기술개발(과제코드 GP2020-010)’의 일환으로 수행되었습니다.

References

Blanco-Martín, L., Rutqvist, J., and Birkholzer, J.T., 2017, Extension of TOUGH-FLAC to the finite strain framework, Computers and Geosciences, 108: 67-71. 10.1016/j.cageo.2016.10.015

Booker, J.R. and Savvidou, C., 1985, Consolidation around a point heat source, International Journal for Numerical and Analytical Methods in Geomechanics, 9: 173-184. 10.1002/nag.1610090206

Chaudhry, A.A., Buchwald, J., Kolditz, O., and Nagel, T., 2019, Consolidation around a point heat source (correction and verification), International Journal for Numerical and Analytical Methods in Geomechanics, 43: 2743-2751. 10.1002/nag.2998

Elyasi, A., Goshtasbi, K., and Hashemolhosseini, H., 2016, A coupled thermo-hydro-mechanical simulation of reservoir CO₂ enhanced oil recovery, Energy and Environment, 27: 524-541. 10.1177/0958305X16665545

Graupner, B., Li, D., and Bauer, S., 2011, The coupled simulator ECLIPSE - OpenGeoSys for the simulation of CO₂ storage in saline formations, Energy Procedia, 4: 3794-3800. 10.1016/j.egypro.2011.02.314

Itasca, 2021, FLAC3D (Fast Lagrangian Analysis of Continua in 3 Dimensions) software, http://www.itascacg.com/software/flac3d [Accessed January 4, 2022].

Kim, G.I., Lee, C., and Kim, J.S., 2021, A numerical study of the performance assessment of coupled thermo-hydro-mechanical (THM) processes in improved korean reference disposal system (kRS⁺) for high-level radioactive waste, Tunnel and Underground Space, 31: 221-242.

Kim, H.M., 2015, Numerical analysis for geomechanical deformation of sea bed due to gas hydrate dissociation, Journal of the Korean Society of Mineral and Energy Resources, 52: 148-157. 10.12972/ksmer.2015.52.2.148

Kim, H.M., Rutqvist, J., Ryu, D.W., Park, D., and Song, W.K., 2012, Numerical study on the impact of air tight liner of underground lined rock cavern for compressed air energy storage (CAES), Journal of the Korean Society of Mineral and Energy Resources. 49: 147-156.

Kim, H.M., Rutqvist, J., Ryu, D.W., Synn, J.H., and Song, W.K., 2011, Geomechanical stability of underground lined rock caverns (LRC) for compressed air energy storage (CAES) using coupled thermal-hydraulic-mechanical analysis, Tunnel and Underground Space, 21: 394-405.

Lim, J., Jung, S., Ki, S., Kim, Y., and Choe, J., 2007, Solute transport modeling using streamline simulation with multiple contaminant sources coupled with MODFLOW, Journal of the Korean Society of Mineral and Energy Resources, 44: 357-363.

Nguyen, T.S. and Selvadurai, A.P.S., 1995, Coupled thermal-mechanical-hydrological behaviour of sparsely fractured rock: implications for nuclear fuel waste disposal, International Journal of Rock Mechanics and Mining Sciences, 32: 465-479. 10.1016/0148-9062(95)00036-G

OpenGeoSys Community, 2021, OpenGeoSys - Open-source multi-physics, http://www.opengeosys.org [Accessed January 4, 2022].

Park, C.H., Kim, T., Park, E.S., Jung, Y.B., and Bang, E.S., 2019, Development and verification of OGSFLAC simulator for hydromechanical coupled analysis: single-phase fluid flow analysis, Tunnel and Underground Space, 29: 468-479.

Rutqvist, J., Wu, Y.S., Tsang, C.F., and Bodvarsson, G., 2002, A modeling approach for analysis of coupled multiphase fluid flow, Heat transfer, and deformation in fractured porous rock, International Journal of Rock Mechanics and Mining Sciences, 39: 429-442. 10.1016/S1365-1609(02)00022-9

Taron, J., Elsworth, D., and Min, K.B., 2009, Numerical simulation of thermal-hydrologic-mechanical-chemical processes in deformable, fractured porous media, International Journal of Rock Mechanics and Mining Sciences, 46: 842-854. 10.1016/j.ijrmms.2009.01.008

UPC (Polytechnic University of Catalonia), 2021, CODE_BRIGHT, http://deca.upc.edu/en/projects/code_bright [Accessed January 4, 2022].

Yoon, S., Um, E.S., and Kim, J., 2020, Development of an integrated simulator for coupled flow, geomechanics, and geophysics based on sequential approaches, Journal of the Korean Society of Mineral and Energy Resources, 57: 309-317. 10.32390/ksmer.2020.57.3.309

Tunnel and Underground Space ISSN:1225-1275(Print) 2287-1748(Online) 터널과 지하공간

Preview

Performance Evaluation of OGS-FLAC Simulator for Coupled Thermal-Hydrological-Mechanical Analysis

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

Fig. 1.

Schematic representation of Booker and Savvidou’s problem in (a) Nguyen and Selvadurai (1995) and (b) Blanco-Martín et al. (2017)

Fig. 2.

Comparison between analytical and modeling results for temperature and pore pressure in (a) Nguyen and Selvadurai (1995) and (b) Blanco-Martín et al. (2017)

Fig. 3.

Numerical model domain and the model grid

(34)

(35)

(36)

Table 1.

Parameters and values used in the present analysis

Fig. 4.

Sequential coupling scheme used in the present study

Fig. 5.

Coupled THM responses with time at point P1 (Case1)

Fig. 6.

Coupled THM responses with time at point P2 (Case1)

Fig. 7.

Coupled THM responses with time at point P1 (Case2)

Fig. 8.

Coupled THM responses with time at point P2 (Case2)

Fig. 9.

Distributions of the ground temperature, pore pressure, displacement, and stress at 105 sec (Case1)

Acknowledgements

References

Distributions of the ground temperature, pore pressure, displacement, and stress at 10⁵ sec (Case1)