Implementation of Barcelona Basic Model Module into OGS-FLAC3D for Simulating Elasto-Plastic Behavior of Unsaturated Soils

Taehyun Kim; Gwang-Il Kim; Changsoo Lee; Jin-Seop Kim

doi:10.7474/TUS.2025.35.2.112

Preview

Technical Note

Tunnel and Underground Space. 30 April 2025. 112-124
https://doi.org/10.7474/TUS.2025.35.2.112

Implementation of Barcelona Basic Model Module into OGS-FLAC3D for Simulating Elasto-Plastic Behavior of Unsaturated Soils

불포화토 탄소성 거동 모사를 위한 OGS-FLAC3D 해석 시뮬레이터 내 Barcelona Basic Model 해석 모듈 구축

Taehyun Kim¹

Gwang-Il Kim¹

Changsoo Lee²^*

Jin-Seop Kim²

김 태현¹

김 광일¹

이 창수²^*

김 진섭²

¹Senior Researcher, Korea Atomic Energy Research Institute

²Principal Researcher, Korea Atomic Energy Research Institute

¹한국원자력연구원 선임연구원

²한국원자력연구원 책임연구원

^{*Corresponding Author}

License (open-access, https://creativecommons.org/licenses/by-nc/4.0/):

This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

ABSTRACT

Buffer materials in a high-level radioactive waste disposal system serve as a key component of the engineered barriers by delaying radionuclide migration and protecting the disposal canister. Over extended periods, these buffer materials undergo a transition from an unsaturated to a saturated state under high-temperature conditions induced by radioactive decay heat and groundwater inflow from surrounding rock. The buffer materials exhibit swelling properties upon saturation, making it critical to predict the long-term capacity of swelling performance throughout the disposal period, which spans hundreds of thousands of years. In this study, the Barcelona Basic Model, a representative model for the elastoplastic behavior of unsaturated soils, was modularized within the thermo-hydro-mechanical coupled numerical simulator, OGS-FLAC3D. Following the construction of the numerical code, we conducted validation through comparative analysis with case studies evaluating volumetric changes along the pathways of mean effective stress and suction variations. The numerical simulation results and validation case studies showed a good agreement. Based on the developed module, we expect to use the OGS-FLAC3D simulator in various research studies related to the behavior of buffer materials in disposal systems.

Keywords

OGS-FLAC3D

Barcelona Basic Model

Elasto-Plastic behavior

Buffer material

Unsaturated soil

고준위 방사성폐기물 처분 시스템 내 완충재는 핵종 유출 시 이동을 지연하고 처분 용기를 보호하는 역할을 하는 공학적 방벽의 핵심 요소이다. 완충재에서는 방사성 붕괴열로 인한 고온 및 주변 암반으로부터의 지하수 유입 조건 하에서 장기간에 걸쳐 불포화 상태에서 포화 상태로 변화하는 과정을 거치게 된다. 완충재는 포화 시 팽윤하는 특성을 지니고 있으며 수십만 년에 이르는 처분 기간 동안 팽윤 성능 요건을 만족할 수 있는지 예측하는 것이 중요하다. 따라서 본 연구에서는 대표적인 불포화토의 탄소성 거동 모델인 Barcelona Basic Model을 열-수리-역학 복합거동 해석 시뮬레이터인 OGS-FLAC3D 내에 모듈화 하였다. 해석 코드 구축 후 검증을 위하여 평균 유효응력 및 흡입력 변화 경로에 따른 비체적 변화를 평가하는 예제와 비교 분석을 수행하였으며, 수치해석 결과 및 검증 예제 결과가 잘 일치하는 것을 확인하였다. 개발된 해석 모듈을 기반으로 향후 OGS-FLAC3D 해석 시뮬레이터를 처분 시스템 내 완충재 거동과 관련한 다양한 연구에 활용할 수 있을 것으로 기대된다.

키워드

OGS-FLAC3D

Barcelona Basic Model

탄소성 거동

완충재

불포화토

MAIN

1. 서 론
2. OGS-FLAC3D 해석 시뮬레이터
3. OGS-FLAC3D 내 Barcelona Basic Model 해석 모듈
3.1 Barcelona Basic Model
3.2 Barcelona Basic Model 내 응력 계산
3.3 Barcelona Basic Model 내 변형률 계산
4. Alonso 예제를 활용한 Barcelona Basic Model 해석 모듈 검증
4.1 Alonso 예제 1
4.2 Alonso 예제 2
4.3 Alonso 예제 3
4.4 Alonso 예제 4
5. 결 론

1. 서 론

고준위 방사성폐기물은 주로 원자력 발전 시 핵연료를 사용하면서 필연적으로 발생하는 사용후핵연료를 포함한 고독성 물질들을 의미한다. 고방사성으로 인해 인체에 치명적인 해를 가할 수 있으며, 따라서 방사성 물질의 반감기를 고려하여 인간 생활권으로부터 수십만 년 이상 격리할 필요가 있다. 고준위 방사성폐기물의 안전한 격리를 위하여 전세계적으로 다양한 방안에 대한 논의가 이뤄져 왔으며, 현재까지 지하 수백 미터 깊이의 암반에 처분장을 건설하고 폐기하는 심층처분 방식이 가장 널리 동의를 얻고 있다. 고준위 방사성폐기물의 안전한 격리를 위해 처분 시스템은 크게 두 가지 요소로 구성되는데, 먼저 인공적으로 제조 및 구성된 공학적 방벽과 모암을 의미하는 천연 방벽으로 이루어진다. 이중 공학적 방벽의 경우, 폐기물을 직접 보관하는 처분 용기와 처분공 및 처분 용기 주변을 채우는 완충재, 처분 터널을 채우는 뒷채움재, 처분 터널 주변의 근계암반 등으로 구성된다. 이 중에서 완충재는 처분공의 빈 공간을 메우는 역할을 하며, 혹시 있을 수 있는 핵종 유출 시 핵종 이동을 지연하는 역할을 하게 된다. 또한 처분 용기를 안전하게 보호하는 역할 또한 담당하게 된다. 완충재의 주재료로는 벤토나이트가 가장 널리 활용되고 있으며, 낮은 투수성과 함께 지하수 유입 시 팽윤 특성을 통해 처분공 벽면과의 빈 공간을 채울 수 있어야 한다.

처분 시 완충재 내에서는 처분 용기에서 발생하는 고온의 열과 주변 암반으로부터 유입되는 지하수 유입이 동시에 발생하며, 따라서 포화도 변화가 복잡하게 발생하게 된다. 앞서 언급한 완충재의 팽윤 특성은 포화도 변화에 따라 발생하며, 팽윤 특성은 처분 시스템의 안정적 작동을 위해 필수적인 요소이므로 적절한 평가를 수행할 필요가 있다. 완충재 팽윤 특성 모사를 위해서 oedometer 시험을 비롯해 다양한 실내 및 현장실험이 수행된 바 있다(Komine and Nobuhide, 1994, Huertas et al., 2000, Börgesson et al., 2001, Wang et al., 2012, Lee et al., 2024, Yoon et al., 2024). 하지만 처분 시스템이 수십 만년 이상 인간 생활권으로부터 안전하게 격리되어야 한다는 점을 고려한다면, 최장 수십 년 정도 수행이 가능한 실내 및 현장실험은 장기 예측 측면에서 한계점을 지니고 있다.

따라서 수치해석을 활용한 연구를 병행하는 것이 필수적이며, 불포화 상태에서의 완충재 팽윤 거동 특성 모사를 위한 다양한 해석모델들이 제시되어 왔다(Alonso et al., 1990, Gens, 1995, Mašín, 2005, Shao et al., 2024). 이중 Barcelona Basic Model (BBM)은 Alonso et al.(1990)에 의해 소개된 이후로 CODE_BRIGHT, TOUGH-FLAC, COMSOL 등 다양한 해석 코드들에 적용되며 완충재 팽윤 거동 모사를 위한 연구에 활발히 활용되고 있다(CODE BRIGHT, 2004, Kristensson and Ákesson, 2008, Rutqvist et al., 2011, Lee et al., 2020, Lee et al., 2024). 이때, 각 해석 코드들이 적용하고 있는 수리 모델, 해석 솔버, 해석기법 등의 차이로 인해 세부적으로는 해석 결과 상의 차이가 발생할 수 있으며, 따라서, 보다 안정적인 장기 진화 예측을 수행하기 위해서는 실제 처분 조건 하에서 다상 유동 시 완충재 팽윤 거동 모사를 위한 다양한 해석 코드를 확보하여 해석해 및 해석 코드 간 결과 비교 검증을 수행하는 과정이 필요하다. 본 연구에서는 처분 시스템 내 복합거동 평가를 위한 연구에 적용되고 있는 OGS-FLAC3D 해석 시뮬레이터(Park et al., 2019, Kim et al., 2024) 내에 BBM 모듈을 구축함으로써 완충재 팽윤 거동 모사 연구에 활용하고자 하며, 구축 완료 후 검증을 위해 수행한 검증 예제 해석 결과를 소개하고자 한다.

2. OGS-FLAC3D 해석 시뮬레이터

OGS-FLAC3D 해석 시뮬레이터는 오픈소스 기반의 THM 복합거동 해석 프로그램인 OpenGeoSys (OGS) ver. 5.8 (Kolditz et al., 2012)와 지반 분야 해석에 강점을 갖고 있는 FLAC3D ver. 7.0 (Itasca, 2012)을 연동한 해석 시뮬레이터로서 지반에서의 복잡한 열-수리-역학(T-H-M) 복합거동 해석을 위해 개발되었다(Park et al., 2019, Park et al., 2020). 두 가지 서로 다른 해석 코드의 장점을 활용하고 연동 시 효율성을 고려하여 순차적 연동해석 방식(sequential coupling method)을 적용하고 있으며, OpenGeoSys에서 열-수리 해석을 수행한 후 해석 결과를 FLAC3D에 전달하면 이를 바탕으로 역학 해석을 수행하고 해석된 결과를 다음 스텝의 OpenGeoSys에 전달하는 방식으로 해석이 진행된다(Fig. 1). 수리-역학 해석을 수행한 초기 버전(Park et al., 2019, Kim et al, 2021a) 이후 Park and Park(2022)이 단상 유동 시 T-H-M 복합거동 벤치마크 문제를 검증한 바 있으며, 최근에는 고준위 방사성폐기물 처분장에서의 복잡한 다상 유동을 고려한 T-H-M 복합거동 문제에서도 적용성을 확인한 바 있다(Kim et al., 2021b, Kim et al., 2024).

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-02/N0120350203/images/ksrm_2025_352_112_F1.jpg

Fig. 1.

Sequential coupling scheme of OGS-FLAC for a THM coupled simulation with timestep

OGS-FLAC3D를 활용하여 고준위 방사성폐기물 처분장 내 공학적방벽 해석을 위해서는 앞서 언급한 바와 같이 일반적인 암석과 달리 포화도 변화에 따라 팽윤 특성을 보이는 완충재 거동을 모사할 수 있는 구성 모델을 확보할 필요가 있다. 따라서 본 연구에서는 팽윤 거동 모사가 가능한 Barcelona Basic Model (BBM)을 OGS-FLAC3D 해석 시뮬레이터 내에 구축하였으며, 다음 장에서는 BBM에 대한 간략한 소개 및 검증 결과를 수록하였다.

3. OGS-FLAC3D 내 Barcelona Basic Model 해석 모듈

Lee et al.(2020)은 기존 BBM 관련 연구인 Alonso et al.(1990), Gens(1995), Rutqvist et al.(2011)을 바탕으로 FLAC의 User Defined Model (UDM)과 내장된 Modified Cam Clay Model (MCCM)을 활용하여 BBM 해석 모듈을 개발하고 검증한 바 있다. 따라서 본 연구에서는 Lee et al.(2020)이 소개한 바 있는 TOUGH2-FLAC3D 시뮬레이터에 삽입된 BBM 모듈을 기반으로 하여 OGS-FLAC3D에서 구동 가능한 BBM 모듈을 개발하였으며, 본 장에서는 BBM 및 OGS-FLAC3D에서 활용 가능하도록 개발된 BBM 모듈에 대해 간단히 소개하고자 한다.

3.1 Barcelona Basic Model

본 논문에서는 OGS-FLAC3D 해석 시뮬레이터의 검증 결과 보고가 주 목적이므로 BBM에 대해서는 간단하게 설명하고자 하며, 보다 자세한 설명은 Lee et al.(2019) 및 Lee et al.(2020)에서 확인할 수 있다. BBM은 포화 상태의 토질 거동까지만 모사가 가능한 MCCM을 확장하여 불포화 상태에서 발생하는 거동까지 모사할 수 있도록 Alonso et al.(1990)이 개발하였다. 따라서 포화 상태, 즉 흡입력이 0인 상태에서는 MCCM 모델과 동일한 항복면을 갖지만, 흡입력이 0이 아닌 불포화 상태에서는 Fig. 2(a)에 나타난 바와 같이 항복면이 증가하게 된다. 즉, MCCM은 전응력(total mean stress, p)과 von-Mises 응력(q)으로 이뤄지는 p-q 면에서의 파괴포락선만으로 구성되지만, BBM은 포화도와 관계된 흡입력(suction, s)이 새로운 축으로 추가되며, 흡입력 증가 시 항복면 함수(Loading-Collapse yield surface, LC)를 따라 항복면의 크기가 증가하는 것을 알 수 있다. 따라서 BBM에서는 Fig. 2(b)에 나타난 바와 같이 변형과 관련된 변수인 정규 압밀 압축 지수 Fig. 2(b)에서 전응력-비체적 관계의 기울기인 $λ_{0}$ 가 흡입력에 따라 변화하며( $λ (s)$ ) , 흡입력에 따른 비체적의 변화, 즉 포화도 변화에 따른 불포화토의 팽윤 거동을 모사할 수 있게 된다.

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-02/N0120350203/images/ksrm_2025_352_112_F2.jpg

Fig. 2.

Conceptual digram of Barcelona basic model yield surface and compression curve (modified after Alonso et al., 1990)

3.2 Barcelona Basic Model 내 응력 계산

OGS-FLAC 시뮬레이터 내 BBM 모듈 구동 시 평균 유효응력 $p'$ 과 von Mises 응력 $q$ 는 Lee et al.(2020)과 동일하게 식 (1) 및 (2)를 통해 계산된다.

(1)

p^{'} = p - p_{p} = \frac{1}{3} (σ_{1} + σ_{2} + σ_{3}) - p_{p}

(2)

q = \sqrt{\frac{{(σ_{1}^{'} - σ_{2}^{'})}^{2} + {(σ_{2}^{'} - σ_{3}^{'})}^{2} {(σ_{3}^{'} - σ_{1}^{'})}^{2}}{2}}

위 식에서 $p'$ 은 평균 유효응력, $p_{p}$ 는 공극수압, $σ_{1, 2, 3}$ 은 주응력, $σ_{1, 2, 3}'$ 은 유효주응력을 의미하며, von Mises 응력은 주응력 방향의 벡터량을 단일한 스칼라 값으로 표현한 등가 응력에 해당한다. OGS-FLAC3D 해석 시뮬레이터 내 다상 유동 해석 시 다양한 해석 모델이 있으며, 현재 개발된 모듈은 Richards flow 모델을 기반으로 하고 있다. Richards flow 모델에서는 포화도를 기반으로 압력을 계산하며, 포화도가 1 미만일 때, 즉 불포화 시에는 압력이 음의 값을 갖게 되고 해당 음의 값이 흡입력으로 계산된다. 이때, 공극수압은 대기압과 동일한 것으로 계산한다. 포화도가 1일 때, 즉 완전 포화 시에는 압력이 양의 값을 갖게 되며, 해당 양의 값이 공극수압으로 작용하고 흡입력은 0으로 계산하게 된다(Kim et al., 2024).

3.3 Barcelona Basic Model 내 변형률 계산

OGS-FLAC3D 내 BBM 모듈에서 변형률 계산은 탄성 변형률( $d ϵ^{e}$ ), 소성 변형률( $d ϵ^{p}$ ), 흡입력에 따른 변형률( $d ϵ^{s}$ ) 및 온도 변화에 따른 변형률( $d ϵ^{T}$ )의 합으로 계산된다. 각 변형률의 계산 방식은 Lee et al.(2020)에 자세한 설명이 있으므로 본 논문에서는 간략히 소개하고자 한다.

체적 탄성 변형률의 변화량( $d ϵ_{v}^{e}$ )은 체적변형계수( $K$ ) 및 평균 유효응력 변화량( $d p'$ )을 통해 식 (3)과 같이 계산되며, 체적변형계수는 흡입력에 따른 과압밀압축지수( $κ_{p s} (s)$ ) 및 간극비( $e$ )를 통해 계산된다(식 (4)).

(3)

d ϵ_{v}^{e} = \frac{1}{K} d p'

(4)

K = \frac{(1 + e) p^{'}}{κ_{p s} (s)}

소성 변형은 Fig. 2(a)에 나타난 바와 같이 온도 및 흡입력 변화에 따른 항복면 함수( $L C$ ) 및 흡입력 증가 시 불포화토의 이전 최대 흡입력( $s_{o}$ )을 초과할 때 발생하게 된다(Alonso et al., 1990). 이를 고려하면 체적 소성 변형률 변화량( $d ϵ_{v}^{p}$ )은 평균 유효응력 변화에 따른 항복면 함수 변화로 인해 발생하는 소성 변형( $d ϵ_{p}^{p}$ )과 흡입력 변화로 인해 발생하는 소성 변형( $d ϵ_{s}^{p}$ )의 합으로 계산된다. 이 때, $d ϵ_{p}^{p}$ 와 $d ϵ_{s}^{p}$ 는 다음 식과 같이 정의된다.

(5)

d ϵ_{p}^{p} = d Λ \frac{\partial g}{\partial p^{'}}

(6)

d ϵ_{s}^{p} = \frac{(λ_{s} - κ_{s p} (p^{'}, s))}{v} \frac{d s_{0}}{s + p_{a t m}}

식 (5)에서 $d Λ$ 는 소성 계수(plastic multiplier) 변화량을 의미하며, $g$ 는 소성 포텐셜 함수를 의미한다. 식 (6)에서 $λ_{s}$ 는 초기 상태에서 흡입력 변화에 따른 강성과 관계된 변수이며, $v$ 는 비체적(specific volume), $κ_{s p} (p', s)$ 는 평균 유효응력( $p'$ )과 흡입력( $s$ )에 따라 변화하는 수리적 팽윤지수, $p_{a t m}$ 은 대기압, $d s_{0}$ 는 기존 불포화토의 최대 흡입력의 변화량을 의미한다.

흡입력에 따른 체적 변형률 $d ϵ_{v}^{s}$ 는 흡입력 강성 계수( $K^{s}$ ) 및 흡입력 변화 $d s$ 를 통해 식 (7)과 같이 계산되며, 온도 변화로 인한 체적 변형률 $d ϵ_{v}^{T}$ 는 체적열팽창계수 $α_{T}$ 및 온도 변화를 통해 식 (8)과 같이 계산된다.

(7)

d ϵ_{v}^{s} = \frac{1}{K^{s}} d s

(8)

d ϵ_{v}^{T} = α_{T} d T

OGS-FLAC3D의 BBM 모듈을 통한 수리-역학 연동 해석 시 공극률 변화는 위와 같이 계산된 체적변형률을 통해 계산된다. 공극률과 간극비는 비체적 변화를 통해 계산되며, 비체적(specific volume)은 식 (9)와 같이 정의된다.

(9)

v = \frac{V}{V_{s}} = \frac{V_{s} + V_{ϕ}}{V_{s}}

위 식에서 $V_{s}$ 는 토질 입자의 부피, $V_{ϕ}$ 는 공극의 부피, $V$ 는 전체 부피에 해당한다. 토질 입자의 부피는 압력의 변화에 따라 거의 변화하지 않지만, 전체 부피는 공극의 부피 변화로 바뀌게 된다. 따라서 위 식과 같이 비체적을 표현하면 상대적으로 일정한 고체 입자의 부피를 기준으로 전체 부피의 변화에 대한 공극 변화의 영향을 설명할 수 있으며 압밀이나 팽윤 거동을 해석할 때 유용하게 활용할 수 있다.

연동 해석 시 비체적의 변화는 앞서 획득한 체적 변형률 변화를 통해 다음과 같이 계산된다. 첨자 $i$ 는 현재 해석 스텝, $i$ +1은 다음 해석 스텝을 의미하며, 공극률과 간극비는 해당 변수의 정의에 기반하여 비체적 변화를 적용하여 식 (11) 및 (12)와 같이 계산된다. 계산된 공극률은 OGS-FLAC3D 연동 해석 시 OpenGeoSys로 전달되어 수리 해석을 위한 질량 보존 방정식(mass conservation equation)의 계산에 활용된다.

(10)

v_{i + 1} = v_{i} (1 - d ϵ_{v}) = v_{i} (1 - \frac{d v}{v})

(11)

e_{i + 1} = \frac{V_{ϕ}}{V_{s}} = v_{i + 1} - 1

(12)

ϕ_{i + 1} = \frac{V_{ϕ}}{V} = \frac{v_{i + 1} - 1}{v_{i + 1}}

4. Alonso 예제를 활용한 Barcelona Basic Model 해석 모듈 검증

본 연구에서는 개발된 해석 모듈의 검증을 위해 Alonso et al.(1990)에서 제시한 불포화토의 거동을 대상으로 한 예제를 활용하여 비교 검증을 수행하였다. BBM은 앞서 설명한 바와 같이 흡입력, 평균 유효응력 및 온도 변화에 따른 불포화토의 체적 변화 거동을 설명하기 위한 모델이며, 따라서 Alonso 예제는 평균 유효응력 또는 흡입력을 각각 일정한 조건이나 변화하는 조건, 그리고 상호 작용하는 조건 하에서 비체적 변화가 어떻게 발생하는지를 살펴보기 위한 예제이다.

Alonso et al.(1990)에서는 논문 상에 그래프로만 결과값이 제시되어 있기 때문에 Alonso et al.(1990)에 수록된 그래프를 추출한 후 OGS-FLAC3D 시뮬레이터 해석 결과를 도시한 그래프의 축을 일치시켜 결과를 비교하였다. 해석 시 입력 물성은 Alonso et al.(1990)과 동일한 물성을 적용하였다(Table 1).

Table 1.

Input parameters for validation simulation based on Alonso et al.(1990)

Variable (unit)	Value	Variable (unit)	Value
$α_{α}$ (-)	0.4	$M$ (-)	1.0
$κ_{s}$ (-)	0.6	$p_{a t m}$ (MPa)	0.1
$κ_{p s 0}$ (-)	0.02	$p^{c}$ (MPa)	0.1
$κ_{s p 0}$ (-)	0.008	$p_{r e f}$ (MPa)	0.1
$λ_{p s 0}$ (-)	0.2	$p_{0}^{*}$ (MPa)	0.2
$λ_{s}$ (-)	0.08	$G$ (MPa)	10.0
$r_{λ}$ (-)	0.75	$v_{i n i}$ (-)	1.9
$β_{λ}$ (MPa^-1)	12.5	$v_{c}$ (0)(-)	2.033

4.1 Alonso 예제 1

Fig. 3(a)는 Alonso et al.(1990)의 첫 번째 예제의 응력 경로를 나타내고 있으며, Fig. 3(b)는 Alonso et al.(1990) 결과와 해석 결과의 비교를 나타내고 있다. Fig. 3(b)를 통해 볼 수 있는 바와 같이 Alonso et al.(1990)의 결과(실선)와 해석 결과(원형 기호)가 정확하게 일치하는 것을 확인할 수 있다.

1번 예제에서는 1번(빨간색 실선), 2번(파란색 실선), 3번(검은색 실선) 응력 경로 순으로 해석을 진행하였으며, 1번 응력 경로(ABDF)에서는 일정한 평균 유효응력(0.15 MPa)에서 흡입력이 감소할 때, 즉 포화가 진행된 후(AB), 완전 포화 상태에서 평균 유효응력을 0.6 MPa 까지 증가시키게 된다 (BD, DF). 이 구간에서는 Fig. 3(b)에서 볼 수 있는바와 같이 포화도 증가로 인한 팽윤으로 인해 비체적이 증가한 후(AB), 평균 유효응력의 증가에 따라 비체적이 감소하게 된다(BDF). 이때 Fig. 3(a)에 도시된 초기 항복면 $L C_{i}$ (P = 0.2 MPa 일 때)을 만나게 되면 소성 거동으로 인해 비체적이 더 급격히 감소하는 것을 확인할 수 있다(BD).

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-02/N0120350203/images/ksrm_2025_352_112_F3.jpg

Fig. 3.

Stress paths scenario 1 and comparison results for decreasing suction at increasing net mean stress

2번 응력 경로(ACDF)에서는 Fig. 3(a)에 나타난 바와 같이 흡입력이 일정한 상태에서 평균 유효응력을 0.35 MPa까지 증가시킨 후(AC), 완전 포화시키는 구간(CD)을 지나 다시 평균 유효응력을 0.6 MPa 까지 증가시키게 된다(DF). 2번 응력 경로 적용 시 Fig. 3(b)에 나타난 바와 같이 평균 유효응력 증가에 따른 압밀로 인해 비체적이 감소하며, 초기 항복면 $L C_{i}$ 을 만난 후 비체적이 급격히 감소하게 된다(AC). 소성 영역에 있는 상태이므로 이후 평균 유효응력은 일정하지만 흡입력이 0.2 MPa에서 0 MPa까지 감소하며 완전 포화되는 CD 구간에서는 소성 변형이 발생하면서 비체적이 급격히 감소하는 붕괴 현상이 발생하게 된다. 포화가 진행됨에도 1번 응력 경로의 AB 구간과 다르게 비체적이 감소하는 것은 소성 영역에서 항복면( $L C$ 곡선)이 변화하면서 발생하는 압축 변형이 팽윤으로 인한 비체적 증가보다 크기 때문인 것으로 나타났다(Lee et al., 2019). 다음으로 평균 유효응력이 증가하는 DF 구간에서는 소성 변형 발생으로 인해 비체적이 감소하는 것으로 나타났다.

마지막 Fig. 3(a)의 3번 응력 경로(ACEF)에서는 평균 유효응력을 0.15 MPa에서 0.6 MPa까지 증가시킨 후(ACE), 완전 포화시키는 경로(EF)를 적용하였다. Fig. 3(b)에 나타난 바와 같이 초기 항복면 $L C_{i}$ 을 만난 이후 소성 변형으로 인해 비체적의 감소가 더 급격히 발생하게 되며(ACE), 포화가 진행되는 EF 구간에서는 앞서 2번 응력 경로의 CD 구간과 마찬가지로 포화로 인한 팽윤보다 소성 변형으로 인한 압축 변형이 더 크게 작용하면서 비체적이 감소하는 것으로 나타났다.

4.2 Alonso 예제 2

Fig. 4(a)는 두 번째 Alonso 예제의 응력 경로를 나타내고 있고 Fig. 4(b)는 해석 결과(원형 기호)와 Alonso et al.(1990) 결과(실선)의 비교를 나타내고 있으며, 두 결과가 잘 일치함을 보여주고 있다.

먼저 첫 번째 응력 경로인 ABDF의 경우 예제 1의 ACEF 응력 경로와 일치하며, 따라서 비체적 변화도 동일하게 나타나는 것을 확인하였다. 두 번째 응력 경로인 ACDF의 경우 Fig. 4(a)에 도시된 바와 같이 흡입력이 0.1 MPa까지 감소하며 포화가 진행된 후, 평균 유효응력이 0.6 MPa까지 증가하며(CD), 이후 일정한 평균 유효응력 상태에서 완전 포화가 진행된다(DF). 해당 구간에서는 먼저 Fig. 4(b)에 나타난 바와 같이 탄성 구간에서 포화로 인한 팽윤으로 비체적이 증가한다(AC). 1번 예제의 응력 경로 AB 구간과 비교 시 포화 정도가 낮기 때문에 비체적 증가량도 적은 것을 확인할 수 있다. 이후 평균 유효응력 증가에 따라 압밀이 진행되며, 평균 유효응력이 약 0.24 MPa인 지점에서 초기 항복면 $L C_{i}$ (Fig. 4(a)의 회색 선)을 접한 이후 소성 변형으로 인해 비체적 변화가 더 크게 일어남을 확인할 수 있다. 다음 단계인 DF 구간에서는 앞서 1번 예제와 마찬가지로 포화로 인한 팽윤보다 소성 거동으로 인한 압축 변형으로 인해 비체적 감소가 더 큰 것을 확인할 수 있다. 세 번째 응력 경로인 ACEF의 경우에도 1번 예제의 ABDF와 일치하며, 따라서 비체적 변화도 동일하게 나타남을 확인하였다. Lee et al.(2020)이 언급한 바와 같이 1, 2번 예제를 통해 다른 응력 경로를 갖더라도 최종 지점인 F까지의 평균 유효응력 및 흡입력 변화가 같고, 포화 상태의 선행압밀 응력 최대값이 같을 때는 비체적 변화가 동일하게 나타남을 알 수 있다.

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-02/N0120350203/images/ksrm_2025_352_112_F4.jpg

Fig. 4.

Stress paths scenario 2 and comparison results for decreasing suction at increasing net mean stress

4.3 Alonso 예제 3

Fig. 5(a)는 세 번째 예제의 응력 경로를 나타내고 있으며, 앞서 1, 2번 예제와 다르게 포화 상태에서 불포화 상태로 변화하는 응력 경로를 나타내고 있다. Fig. 5(b)는 비교 결과를 도시하였으며, 해석 결과(원형 기호)와 Alonso et al.(1990) 결과(실선)가 정확하게 일치함을 확인할 수 있다.

먼저 첫 번째 응력 경로인 ABDF 구간에서는 완전 포화 상태에서 평균 유효응력이 증가함에 따라 비체적이 감소하게 되며(Fig. 5(b)의 AB), 초기 항복면(P=0.2 MPa 일 때 $L C_{i}$ )을 만난 이후 비체적 감소가 더 크게 발생한다. ABDF 구간에서의 최종 항복면은 포화 상태의 선행 압밀응력(B 지점에서의 평균 유효응력)이 가장 큰 B 지점에서의 $L C_{B}$ 이므로 BDF 경로에서는 탄성 영역에 위치하게 된다. 따라서 흡입력 증가, 즉 포화도 감소로 인해 수축이 발생함을 확인할 수 있다(Fig. 5(b)의 BF).

두 번째 응력 경로인 Fig. 5(a)의 ACDF 구간에서는 탄성 구간에서 초기 포화도 감소로 인해 수축이 발생하며(Fig. 5(b)의 AC), 이후 평균 유효응력이 증가함에 따라 비체적이 감소하게 된다. 초기 항복면을 지나며 비체적 감소량이 증가하며, D 지점 도달 시 최종 항복면 $L C_{D}$ 를 갖게 된다. 따라서 DF 구간에서는 탄성 영역에 위치한 상태로 포화도 감소로 인한 압축이 발생하며, 흡입력 변화량이 BF 보다 적기 때문에 첫 번째 응력 경로의 BF 구간에 비해 압축량은 더 작게 나타나는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-02/N0120350203/images/ksrm_2025_352_112_F5.jpg

Fig. 5.

Stress paths scenario 3 and comparison results for decreasing suction at increasing net mean stress

마지막 세 번째 응력 경로인 ACEF 구간에서는 흡입력이 증가함에 따라(ACE) 압축이 발생하게 되고(Fig. 5(b)의 AE), 이후 평균 유효응력 증가로 인해 비체적 감소가 발생하게 된다. 앞서와 마찬가지로 초기 항복면 이후 비체적 감소가 더 크게 발생하며, 최종적으로 항복면 $L C_{F}$ 를 갖게 된다.

예제 1 및 2와 다른 점은 최종 지점은 F로 동일하지만 최종적으로 갖는 항복면의 크기가 $L C_{B}$ < $L C_{D}$ < $L C_{F}$ 순으로 결정되기 때문에 최종 비체적 변화, 즉 소성 변형의 크기도 이에 따라 1번 > 2번 > 3번 응력 경로 순으로 나타남을 확인할 수 있다. 즉, 시료가 받게 되는 최종 선행 압밀응력이 클수록 더 많은 압축 변형을 받게 되므로 소성 변형도 더 많이 발생하게 됨을 알 수 있다.

4.4 Alonso 예제 4

네 번째 예제에서는 앞서 초기 불포화토 최대 흡입력이 0.3 MPa이었던 것과 달리 초기 불포화토의 최대 흡입력을 0.025 MPa로 설정하였다. 응력 경로는 Fig. 6(a)에 나타난 바와 같이 완전 포화 상태에서 평균 유효응력을 증가시키는 1번 경로(AB), 이후 일정한 평균 유효응력(0.015 MPa) 하에서 흡입력이 증가 후 다시 포화되는 ACD 구간을 지나 평균 유효응력을 0.6 MPa까지 증가시키는 DE 구간으로 구성되어 있다. 해석 결과와 Alonso et al.(1990)의 결과는 Fig. 6(b)와 같이 잘 일치하는 것을 확인할 수 있다.

먼저 1번 응력 경로에서는 앞선 예제들과 마찬가지로 평균 유효응력 증가로 인해 비체적이 초반 감소하며, 초기 항복면 $L C_{i}$ (P = 0.2 MPa 일 때)를 만난 후 비체적 감소량이 더 크게 증가하게 된다. 2번 응력 경로에서는 먼저 일정 평균 유효응력 하에서 흡입력이 증가하게 되는데, 이로 인해 포화도 감소로 수축이 발생하므로 비체적은 감소하게 된다. 이때 초기 최대 흡입력 0.025 MPa을 넘어서게 되면 소성 변형이 발생하여 3번 예제에서의 ACE 응력 경로보다 더 큰 비체적 감소가 발생하며 새로운 포화 상태 선행압밀 응력을 갖게 되어 항복면을 새롭게 형성하게 된다. 흡입력 0.3 MPa에 도달 후 다시 완전 포화 과정으로 돌아오는 CD 경로에서는 탄성 구간에 위치하기 때문에 포화로 인한 비체적의 증가만 발생하게 된다. 마지막으로 포화 상태에서 평균 유효응력을 증가시키는 DE 경로에서는 압축으로 인한 비체적 감소가 발생하며, 앞서 최대 흡입력 변화로 인해 새롭게 형성된 항복면 $L C_{D}$ 를 지나면서 소성 변형으로 인해 비체적의 급격한 감소가 발생하는 것을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/ksrm/2025-035-02/N0120350203/images/ksrm_2025_352_112_F6.jpg

Fig. 6.

Stress paths and comparison results for decreasing suction at increasing net mean stress (Scenario 4)

5. 결 론

본 논문에서는 고준위 방사성폐기물 처분 시스템의 핵심 구성요소인 완충재의 팽윤 거동을 모사하기 위한 수치해석적으로 모사하기 위한 연구를 진행하였다. 해석 코드로는 순차적 연동해석 방식에 기반한 OGS-FLAC3D 해석 시뮬레이터를 활용하였으며, 팽윤 거동 모사를 위한 모델로는 BBM을 선정하였다. BBM은 탄소성 토질 모델인 MCCM 모델에 흡입력 변화 및 열적 변화에 따른 변형을 추가하여 팽윤 거동을 모사할 수 있도록 개발된 모델이며, 본 연구에서는 기존 연구를 바탕으로 FLAC3D의 UDM 및 FISH 함수를 활용하여 BBM 모듈을 구축하였다.

구축한 모듈의 검증을 위해서 대표적으로 활용되고 있는 Alonso et al.(1990)의 네 가지 검증 예제 결과와 해석 결과를 비교하였다. 각 검증 예제는 평균 유효응력 및 흡입력 변화 경로에 따라 비체적이 변화하는 양상을 모사하고 있으며, 각 예제에서는 평균 유효응력 증가 시 비체적이 감소하고 초기 항복면을 만날 시 소성 변형으로 인해 비체적 변화가 더 급격히 발생하는 것을 확인하였다. 또한 흡입력 증가 시(시료 건조 시) 수축 변형, 흡입력 감소 시(시료 포화 시) 팽창 변형이 발생하는 것을 확인하였으며, 탄성 혹은 소성 구간에 위치하는지에 따라 시료의 변형량이 달라질 수 있음을 확인하였다.

해석 결과는 검증 예제에서 제시된 결과와 정확히 일치하는 것으로 나타났으며, 이를 통해 OGS-FLAC3D 내에 구축한 BBM 모듈이 안정적으로 작동하는 것을 확인하였다. 개발된 연구 결과를 기반으로 향후 OGS-FLAC3D 해석 시뮬레이터를 보다 다양한 고준위 방사성폐기물 처분 관련 연구에 적용할 수 있을 것으로 기대되며, 후속 연구에서 실제 실험 결과와 해석 결과를 비교함으로써 활용성을 확장하고자 한다.

Acknowledgements

본 연구는 과학기술정보통신부의 재원으로 사용후핵연료관리핵심기술개발사업단(2021M2E1A1085193)의 지원을 받아 수행되었습니다.

References

Alonso, E.E., Gens, A., and Josa, A., 1990, A constitutive model for partially saturated soils, Geotechnique, 40(3), 405-430.

10.1680/geot.1990.40.3.405

Börgesson, L., Chijimatsu, M., Fujita, T., Nguyen, T.S., Rutqvist, J., and Jing, L., 2001, Thermo-hydro-mechanical characterisation of a bentonite-based buffer material by laboratory tests and numerical back analyses, International Journal of Rock Mechanics and Mining Sciences, 38(1), 95-104.

10.1016/S1365-1609(00)00067-8

CODE BRIGHT, 2004, A 3D program for thermo-hydro-mechanical analysis in geological media manual, Barcelona, UPC.

Gens, A., 1995, Constitutive laws. In: Gens, A., Jouanna, P., Schrefler, B.A. (Eds.), Modern Issues in Non-saturated Soils, Springer-Verlag, New York, 129-158.

10.1007/978-3-7091-2692-9_2

Huertas, F., Fuentes-Cantillana, J.L., Jullien, F., Rivas, P., Linares, J., Fariña, P., Ghoreychi, M., Jockwer, N., Kickmaier, W., Martínez, M.A., and Samper, J., 2000, Full-scale engineered barriers experiment for a deep geological repository for high-level radioactive waste in crystalline host rock(FEBEX project), Final report, EUR(Luxembourg).

Itasca, 2012, FLAC3D Version 5.0, Fast Lagrangian Analysis of Continua in 3 Dimensions, Constitutive Models, Itasca Consulting Group, Minnesota.

Kim, T., Park, C.-H., Watanabe, N., Park, E.-S., Park, J.-W., Jung, Y.-B., and Kolditz, O., 2021a, Numerical modeling of two‑phase flow in deformable porous media: application to CO2 injection analysis in the Otway Basin, Australia, Environmental Earth Sciences, 80(121), 1-15.

10.1007/s12665-021-09411-1

Kim, T., C.-H. Park, C. Lee and J.-S. Kim, 2021b, A numerical study on the Step 0 benchmark test in Task C of DECOVALEX-2023: Simulation for thermo-hydro-mechanical coupled behavior by using OGS-FLAC, Tunnel & Underground Space, 31(6), 610-622.

Kim, T., Park, C.-H., Lee, C., Kim, J.-S., Park, E.-S., and Bastian, G., 2024, A numerical analysis of Thermo-Hydro-Mechanical behavior in the FE experiment at Mont Terri URL: Investigating capillary effects in bentonite on the disposal system, Geomechanics for Energy and the Environment, 40, 100597.

10.1016/j.gete.2024.100597

Kolditz, O, Bauer, S., Bilke, L., Böttcher, N., Delfs, J.O., Fischer, T., Görke, U.J., Kalbacher, T., Kosakowski, G., McDermott C C. et al., 2012, Opengeosys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (thm/c) processes in porous media, Environmental Earth Sciences, 67(2), 589-599.

10.1007/s12665-012-1546-x

Komine, H., and Nobuhide, O., 1994, Experimental study on swelling characteristics of compacted bentonite, Canadian Geotechnical Journal, 31(4), 478-490.

10.1139/t94-057

Kristensson, O., and M. Ảkesson, 2008, Mechanical modeling of MX-80 - quick tools for BBM parameter analysis, Physics and Chemistry of the Earth, 33, 508-515.

10.1016/j.pce.2008.10.003

Lee, C., Lee, J., Kim, M., and Kim, G.Y., 2020, Implementation of Barcelona Basic Model into TOUGH2-MP/FLAC3D, Tunnel & Underground Space, 30(1), 39-62.

Lee, C., Yoon, S., Lee, J., and Kim, G.Y., 2019, Introduction of Barcelona Basic Model for Analysis of the Thermo-Elasto-Plastic Behavior of Unsaturated Soils, Tunnel & Underground Space, Vol. 29(1), 38-51.

Lee, D.-H., Go, G.-H., Lee, G.-J., and Yoon, S., 2024, Swelling behavior simulation study of KJ-II bentonite buffer blocks under various experimental conditions, Journal of the Korean Geotechnical Society, 40(2), 29-40.

Mašín, D., 2005, A hypoplastic constitutive model for clays, International Journal for numerical and analytical methods in Geomechanics, 29(4), 311-336.

10.1002/nag.416

Park, C.-H., Kim, T., Park, E.-S., Jung, Y.-B., and Bang, E.-S., 2019, Development and Verification of OGSFLAC Simulator for Hydro-mechanical Coupled Analysis: Single-phase Fluid Flow Analysis, Tunnel & Underground Space, 29(6), 468-479.

Park, C.-H., Kim, T., Park, E.-S., Jung, Y.-B., and Bang, E.-S., 2020, Construction of open-source program platform for efficient numerical analysis andits case study, Tunnel & Underground Space, 30(6), 509-518.

Park, D. and Park, C.-H., 2022, Performance evaluation of OGS-FLAC simulator for coupled thermal-hydrological-mechanical analysis, Tunnel & Underground Space, 32(2), 144-159.

Rutqvist, J., Ijiri, Y., and Yamamoto, H., 2011, Implementation of the Barcelona Basic Model into TOUGH-FLAC for simulations of the geomechanical behavior of unsaturated soils, Computers & Geosciences, 37, 751-762.

10.1016/j.cageo.2010.10.011

Shao, H., Rhoelke, C.h., Emmerich, K., Popp, T., Koeniger, F., Schuhmann, R., Wieczorek, K., and Hesser, J., 2024, Coupled hydromechanical modelling of a vertical hydraulic sealing system based on the sandwich principle, Rock Mechanics and Rock Engineering, 57(6), 4355-4365.

10.1007/s00603-023-03404-1

Wang, Q., Tang, A.M., Cui, Y.J., Delage, P., and Gatmiri, B., 2012, Experimental study on the swelling behaviour of bentonite/claystone mixture, Engineering Geology, 124, 59-66.

10.1016/j.enggeo.2011.10.003

Yoon, S. and Lee, G.-J., Chang, S., and Lee, D.-H., 2024, Assessment of thermal-hydraulic-mechanical-chemical (THMC) performances of Ca-type bentonite-graphite mixtures as buffer materials for a high-level radioactive waste repository, Case studies in Construction Materials, 21, e03876.

10.1016/j.cscm.2024.e03876

Tunnel and Underground SpaceISSN:1225-1275(Print) 2287-1748(Online)한국암반공학회

Preview

Implementation of Barcelona Basic Model Module into OGS-FLAC3D for Simulating Elasto-Plastic Behavior of Unsaturated Soils

ABSTRACT

MAIN

Fig. 1.

Sequential coupling scheme of OGS-FLAC for a THM coupled simulation with timestep

Fig. 2.

Conceptual digram of Barcelona basic model yield surface and compression curve (modified after Alonso et al., 1990)

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

Table 1.

Input parameters for validation simulation based on Alonso et al.(1990)

Fig. 3.

Stress paths scenario 1 and comparison results for decreasing suction at increasing net mean stress

Fig. 4.

Stress paths scenario 2 and comparison results for decreasing suction at increasing net mean stress

Fig. 5.

Stress paths scenario 3 and comparison results for decreasing suction at increasing net mean stress

Fig. 6.

Stress paths and comparison results for decreasing suction at increasing net mean stress (Scenario 4)

Acknowledgements

References