1. 서 론

고준위 방사성폐기물의 지층처분에서는 장기적인 처분 환경의 변화 속에서 다중방벽시스템(multi-barrier system)이 안정적으로 기능을 유지하도록 하는 것이 주요 과제이다. 특히 공학적 방벽(engineered barrier system, EBS)의 핵심 재료인 벤토나이트는 수화에 따른 팽윤과 투수성 저하로 핵종 이동을 억제하지만, 실제 처분장에서는 불포화 상태, 장기 수화, 열적 영향, 암반–밀봉체 경계 효과 등으로 열–수리–역학(thermo-hydro-mechanical, THM) 상호작용이 누적·교차한다. 이러한 환경 하에서 응력–변형과 항복 메커니즘을 정량적으로 예측하기 위해서는, 불포화 토질 역학에 기초한 구성모델과 검증 가능한 복합거동 해석 기법 개발이 매우 중요하다.

벤토나이트와 같은 팽창성 점토는 불포화 조건에서 흡입력(suction) 변화에 따라 수축과 팽창을 반복하며, 습윤–건조 사이클 동안 비가역적인 체적변형을 나타낸다. 이러한 거동 특성을 규명하고 예측하기 위해 다양한 이론적·수치적 연구가 진행되어 왔다. 초기에는 Bishop(1959)의 유효응력 확장 개념이나 Fredlund and Morgenstern(1977)의 전단강도식처럼 강도 해석에 중점을 둔 단순화된 접근이 있었고, 이후 Prevost(1978)의 cap model(항복면을 모자(cap) 형태로 확장하여 체적 변형을 설명하는 탄소성 모델)이 다축 응력 조건에서의 압밀 및 변형 거동을 설명하기 위해 도입되었다. Roscoe et al.(1963)의 Cam-Clay 모델이 임계상태 이론에 기초하여 포화 점토의 응력–변형 거동을 설명하는 대표적 구성모델로 제안되었으며, Schofield and Wroth(1968)에 의해 수정된 Modified Cam-Clay (MCC)는 항복 조건을 단순화하여 실제 점토의 압밀·전단 응답을 보다 잘 설명할 수 있는 표준 모델로 발전하였다. Wheeler and Sivakumar(1995)는 이러한 MCC 모델을 불포화 조건으로 확장하여 흡입력 변화에 따른 소성 붕괴(collapse) 및 경로 의존적 거동을 포착하려 하였으나, 이들 모델은 불포화 토양의 주요 특성을 부분적으로만 재현하였으며, 벤토나이트의 장주기 수화–건조 과정에서 나타나는 비가역적 체적 변형 거동을 충분히 설명하는 데에는 제약이 있었다.

이러한 한계를 보완하기 위해 Alonso et al.(1990)은 흡입력을 독립 상태변수로 도입하고, 경화법칙과 비연관 유동 법칙을 채택한 Barcelona Basic Model (BBM)을 제안하였다. BBM은 Loading–Collapse (LC)와 Suction Increase (SI) 항복곡면을 통해 건조–습윤 과정에서의 비가역 거동을 설명할 수 있으며, 포화 조건에서는 MCC 모델과 일치하는 이론적 연속성을 갖는다. 이후 Gens and Alonso(1992)는 BBM의 틀을 확장하여 불포화 팽창성 점토의 거동을 설명하기 위한 체계적 이론을 제시하였고, 이를 기반으로 Alonso et al.(1999)은 미세–거시 공극 구조(double structure)를 구분하여 팽윤압과 장기 체적변형을 정교하게 반영할 수 있는 Barcelona Expansive Model (BExM)을 제안하기도 하였다. BBM은 제안된 이후 불포화 점토 거동을 설명하는 대표적 구성모델로 자리매김하였으며, 국제적으로 다양한 수치해석 코드에 구현되어 벤토나이트 완충재의 열–수리–역학적 거동을 모사하는 핵심 구성모델로 활용되어 왔다. 특히 카탈루냐 공대(UPC)에서 개발된 CODE_BRIGHT는 BBM과 BExM의 이론적 기반을 활용하여, 불포화 팽창성 점토의 거동을 정밀하게 분석할 수 있는 주요 모델로 발전하였다(Gens et al., 1998, Sánchez et al., 2012). 또한, Mont Terri Underground Research Laboratory (URL)와 같은 실험과 DECOVALEX 국제 공동연구를 통해 벤토나이트의 장기 거동 분석에 널리 활용되고 있다. 한편, Rutqvist et al.(2011)은 BBM을 FLAC3D (Itasca Consulting Group, Inc., 2025) 내에 사용자 정의 모델(user-defined model, UDM)로 구현하고 이를 다상유동 해석코드 TOUGH2 (Pruess et al., 1999)와 연동하여 불포화 조건에서의 열–수리–역학(thermo– hydro–mechanical, THM) 해석 기법을 제안하였다. 국내에서도 Lee et al.(2020a, 2020b)이 BBM에 기반한 TOUGH–MP/FLAC3D 연동해석을 통해 스위스 Grimsel Test Site의 FEBEX 현장실험을 모사하였으며, 최근에는 Kim et al.(2025)이 FLAC3D에 구현된 BBM을 오픈소스 코드인 OpenGeoSys (OGS, Kolditz et al., 2012)와 연계하기 위한 기초 연구를 수행한 바 있다.

본 연구는 벤토나이트의 수리–역학적 거동과 주변 암반과의 상호작용을 모사하기 위한 수치해석 연구로서, 국제 공동연구인 DECOVALEX (DEvelopment of COupled models and their VALidation against EXperiments) 프로젝트의 일환으로 수행되었다. 이를 위해 BBM 구성식을 FLAC3D의 UDM을 통해 구현하고, OGS와 FLAC을 연동하여 불포화토의 수리–역학적 거동을 해석할 수 있는 기법을 개발하였다. 구현된 해석 모듈의 이론적 정합성과 수치적 신뢰성을 검증하기 위하여 포화 조건에서의 등방 재하–제하 응답을 MCC와 비교하고, Alonso et al.(1990)의 대표적인 응력 경로(Cases 1–3)를 재현하였다. 또한, 개발된 해석모듈을 DECOVALEX–2027 Sandwich Task의 Step 0a 벤치마크 연구에 적용하여 실험실에서 얻어진 불포화 벤토나이트의 거동 특성을 모사하였다.

2. 벤토나이트 탄소성 거동해석을 위한 BBM 모델 개발 및 검증

2.1 Barcelona Basic Model 이론

BBM은 Alonso et al.(1990)에 의해 제안된 구성 모델로, 불포화 토양(unsaturated or partially saturated soils)의 응력 변형률 거동을 체계적으로 설명하기 위해 개발되었다. 기존의 MCC 모델은 포화 점토의 거동을 효과적으로 모사하지만 불포화 조건에서 나타나는 흡입력(suction)에 의한 강도 증가, 붕괴, 팽창(swell), 그리고 경로 의존적 거동(path dependence)을 설명하는 데에는 한계가 있다. BBM은 MCC 모델을 확장하여 흡입력을 독립 변수로 도입함으로써 포화 및 불포화 상태 전반에 걸쳐 적용 가능한 통합 모델을 제시한다. 이 모델은 경화 소성 이론(hardening plasticity)을 기반으로 하며, 두 개의 독립적인 응력 변수인 순평균응력(net mean stress, p)과 흡입력(suction, s)을 핵심으로 한다. 이 변수들은 물리적으로 명확한 의미를 가지며, 실험적으로 독립적 제어가 가능하다는 점에서 모델의 기초 변수로 채택되었다. 또한, BBM에서는 전단응력을 고려한 항복 조건을 기술하기 위해, 편차응력(deviatoric stress, q)을 포함한 삼차원 응력 공간 (p,s,q) 상에서 항복면을 정의한다. 각 변수를 수식으로 표현하면 다음과 같다.

여기서 σm은 평균응력, 𝜎₁, 𝜎₂, 𝜎₃은 각각 최대주응력, ua는 공기압, uw는 수압이다.

BBM은 두 개의 항복 곡면(yield surfaces), 즉 Loading–Collapse (LC) 곡면과 Suction Increase (SI) 곡면을 통해 불포화 토양의 탄소성 거동을 정의한다. Fig. 1은 삼차원 응력 공간 (p,q,s)상에 LC 곡면과 SI 곡면을 도시한 것으로, 두 곡면의 경계 내부에 위치한 응력상태에서는 탄성거동을, 외부에서는 소성 거동을 보이게 된다. 흡입력의 증가는 토양의 강성(stiffness)을 증가시키고 붕괴 거동(collapse behavior)을 유발할 수 있으며, 주로 팽창 또는 수축과 같은 체적 변형에 영향을 미친다. 반면, 흡입력은 전단 변형에는 직접적인 영향을 주지 않는 것으로 가정된다. 이러한 가정 하에서, 탄성 상태에서의 체적변형률 dϵve와 전단 변형률은 dϵvs는 각각 다음과 같은 식으로 표현된다.

여기서 dϵvpe와 dϵvse는 각각 순평균응력에 의한 탄성변형률과 흡입력에 의한 탄성변형률이며, 𝜅는 순평균응력에 대한 탄성 강성 파라미터(elastic stiffness parameter for changes in net mean stress; slope of swelling line), κs는 흡입력에 대한 탄성 강성 파라미터(elastic stiffness parameter for changes in suction), v는 비체적(specific volume), G는 전단탄성계수, pat는 대기압이다. 이때 강성 파라미터는 일반적으로 사용하는 강성(stiffenss) 개념보다는 압축률(compressibility)의 개념에 가까우며, 그 값이 클수록 순평균응력 또는 흡입력의 변화에 따른 변형이 크게 발생한다. 한편, v는 토양의 단위 부피당 전체 부피를 의미하며, 공극비 e와 v=1+e의 관계에 있다. pat는 수치적 안정성을 보장하기 위해 포함된 보정 항이다.

Fig. 1.

Three-dimensional view of the yield surfaces in (p,q,s) stress space (after Alonso et al., 1990)

첫 번째 항복면인 LC 곡면은 외부 응력 증가 또는 흡입력 감소에 따른 붕괴 및 체적 수축과 전단 응력에 의한 항복 조건을 포함하며, 이차원 응력 공간 (p,s) 상에서 다음과 같이 표현된다.

여기서 p0는 현재 흡입력 s에서의 가상 선행압밀응력(preconsolidation stress), p0*는 포화상태(s=0)에서의 선행압밀응력, pc는 기준응력(reference pressure)이다. 이때 λ(s)는 소성 상태(normally consolidated state or virgin state)에서 순평균응력에 대한 강성 파라미터(stiffness parameter for changes in net mean stress for virgin states at suction s)로서 흡입력 증가에 따른 재료의 강성 증가를 반영하며, 식 (4)로 표현된다.

여기서 𝜆(0)는 포화 상태에서의 소성 강성 파라미터(stiffness parameter for changes in net mean stress for virgin states at zero suction; slope of normal consolidation line), r은 최대 강성에 수렴할 때의 비율 상수(0<r<1), 𝛽는 강성 증가 속도를 조절하는 실험적 계수이다.

LC 항복 곡면( f1)은 전단 거동을 고려하기 위해 삼차원 응력 공간 (p,q,s)에서 다음과 같은 타원형 방정식으로 정의된다.

여기서 M은 임계 상태 라인의 기울기, ps=ks는 흡입력 s에 비례한 점착력 증가항, k는 점착력 증가계수이다. 흡입력이 증가하면 항복 곡면은 우측으로 이동하게 되며, 동일한 순평균응력에서 흡입력에 따라 전단강도가 증가하는 효과를 정량적으로 반영할 수 있다.

BBM의 두 번째 항복 곡면은 SI 항복곡면( f2)으로, 흡입력 s가 과거 경험한 최대치 s0를 초과할 때 비가역적인 수축 변형이 발생하게 되며, 다음과 같이 단순한 형태로 표현할 수 있다.

BBM에서 정의된 LC 곡면과 SI 곡면은 원칙적으로 상호 독립적인 항복 조건을 규정하지만, 실험적으로는 소성 체적 변형률의 누적에 따라 두 곡면이 서로 연동(coupling)하여 이동하는 경향이 나타난다. 이러한 상호 연계성을 반영하기 위해, BBM은 두 곡면에 대해 총 소성 체적 변형률(dϵvp)을 기반으로 한 경화 법칙(hardening law)을 적용한다. 이 경화 법칙은 주어진 소성 변형 조건에 따라 포화 상태 선행압밀응력(p0*)과 최대 흡입력(s0)의 진화를 유도하며, 이에 따라 항복 곡면의 위치가 결정된다. 각 경화 파라미터의 진화는 소성 체적 변형률(dϵvp)에 대한 함수로 정의되며, 이를 통해 LC와 SI 곡면은 각각의 기하학적 조건을 유지하면서도 경화 상태에 따라 연속적으로 변화하게 된다. 이때 경화 법칙은 다음과 같이 표현할 수 있다.

여기서 λs는 흡입력에 의한 소성 강성 파라미터(stiffness parameter for changes in suction for virgin states)이다.

총 소성 체적 변형률 dϵvp는 LC 곡면과 SI 곡면에서 각각 독립적으로 발생하는 소성 체적 변형률의 합(dϵvp=dϵv,LCp+dϵv,SIp)으로 표현된다. 이때 dϵv,LCp는 응력 변화에 따른 붕락 또는 전단 거동에 의해 발생하며, dϵv,SIp는 흡입력 증가로 유도된 수축 거동에 의해 발생한다. 항복 상태에서의 소성 변형률은 유동 법칙(flow rule)에 의해 결정되며, BBM은 MCC 모델과는 달리 비연관 유동 법칙(non-associated flow rule)을 채택한다. MCC 모델은 연관 유동 법칙(associated flow rule)을 사용함으로써 전단 팽창(shear dilation)을 과대평가하는 경향이 있으며, 정상압밀 상태에서의 측압계수 K0가 Jaky 상태(Jaky, 1948)의 예측값과 일치하지 않는 결과를 초래하는 것으로 알려져 있다. BBM은 이러한 문제를 개선하기 위하여 소성 포텐셜(plastic potential)을 항복 함수(yield function)와 독립적으로 정의하고, 비연관 파라미터 𝛼를 도입하여 측방 변형률이 0인 조건을 만족하도록 유동 방향을 조정하였다.

LC 항복면에서의 소성 변형률 벡터는 체적 변형률(dϵv,LCp)과 전단변형률 (dϵv,LCp)로 구성되며, 비연관 유동 법칙에 따라 소성 변형률 방향이 다음과 같이 정의된다.

여기서 𝛼는 측방 변형률이 0인 조건을 만족하도록 설정되며 다음과 같이 표현된다.

소성 변형률의 방향은 소성포텐셜 g1의 법선 방향을 따르며, 이에 따라 소성 체적변형률과 소성 전단변형률은 다음의 유동 방정식을 만족한다.

여기서 𝛬₁는 소성승수(plastic multiplier)로, 항복 조건이 계속 유지되도록 보장하는 일관성 조건(consistency condition, df1=0)에 따라 결정되며, 주어진 유동 방향에 대한 소성 변형률의 크기를 정량적으로 조절하는 계수이다.

Alonso et al.(1990)의 논문에서는 소성 포텐셜 g1의 명시적 형태를 제공하지는 않았지만, 식 (9)에 따라 암시적으로 다음과 같은 형태가 가정된 것임을 확인할 수 있다.

한편, SI 항복면의 소성 변형은 흡입력 증가로 인한 체적 변형률만을 유발하며, 전단 변형에는 영향을 미치지 않는다. SI 항복에 의한 체적 변형률은 일관성 조건(df2=0)과 s0 경화법칙(식 (8))을 결합함으로써 다음과 같이 계산할 수 있다.

여기서 𝛬₂는 SI 항복에 대한 소성승수이다.

2.2 BBM 해석모듈 개발 및 검증

2.2.1 FLAC UDM을 이용한 해석모듈 개발

본 연구에서는 FLAC3D 소프트웨어의 UDM 기능을 활용하여 BBM의 코드와 알고리즘을 구현하고 그 정확성을 검증하였다. FLAC3D는 암반 및 지반공학 분야에서 복잡한 지반 거동과 안정성 해석을 수행하기 위해 널리 사용되는 유한차분법(Finite Difference Method, FDM) 기반의 수치해석 소프트웨어로, UDM을 통해 특정 재료의 구성 모델을 맞춤형으로 구현할 수 있는 유연성을 제공한다. FLAC3D의 UDM은 동적 링크 라이브러리(dynamic link library, DLL) 형태로 구현된다. BBM의 구현은 C++ 프로그래밍 언어를 이용하여 수행되었으며, 개발된 코드는 BBM의 핵심 구성 방정식, 즉 LC 항복 곡면(식 (5)), SI 항복 조건(식 (6)), 경화 법칙(식 (7), (8)), 비연관 유동 법칙(식 (10))을 포함하도록 설계되었다. 완성된 코드는 플러그인(plugin) 형태로 FLAC3D에 통합되어, 내장된 기존 구성 모델과 동일한 방식으로 작동한다.

본 연구에서는 FLAC3D에서 제공되는 MCC 모델의 소스코드를 기반으로, 불포화 조건에서의 BBM 거동을 반영하기 위하여 입력 파라미터, 항복면 구성, 탄성 및 소성 변형률 계산 알고리즘을 확장하였다. 특히, 흡입력 증가에 따른 강성 및 전단 강도 변화를 모사하기 위하여 k, 𝛽, r을 입력 파라미터로 추가하였고, 𝜅, λ(s), κs, λs, M, p0*, s0 등 주요 파라미터를 모델 상태 변수로 등록하였다. 또한, 유동 법칙을 구현하기 위해 Jaky 상태 조건을 기반으로 한 비연관 파라미터인 𝛼를 계산하여 소성 변형률 계산식에 반영하였다.

FLAC3D의 응력-변형률 계산 알고리즘은 유한차분법에 기반한 명시적 시간 적분(explicit time-stepping scheme)을 사용하며, 각 해석 요소(zone)에 대해 주어진 변형률 증분에 대응하는 응력 증분을 계산하는 구조를 가진다. 그러나 BBM의 경우, 탄성 상태에서 흡입력에 의해 발생하는 체적 변형률(dϵv,se)은 식 (14)와 같이 응력 변화와는 독립적으로 유발되므로, 흡입력 변화에 따른 FLAC3D의 기본 계산 절차에 직접 반영하는 것이 구조적으로 불가능하다.

이러한 문제를 해결하기 위해 본 연구에서는 흡입력 변화에 의해 유도되는 탄성 체적 변형률에 대응하는 평균 응력 증분을 이론적으로 산정하고, 이를 응력 계산 과정에서 주응력 텐서에 추가하는 방식으로 흡입력의 효과를 구현하였다. 이러한 접근법은 BBM의 이론을 만족하면서도 FLAC3D의 계산 구조와의 호환성을 확보할 수 있도록 한다. 또한, FISH 스크립트를 이용하여 흡입력 효과를 구현하였던 선행 연구들(Rutqvist et al., 2011, Lee et al., 2020a)과는 달리, 본 연구에서는 이를 DLL 기반의 UDM 코드에 직접 반영함으로써 계산 효율성과 수치적 안정성을 동시에 향상시켰다.

식 (14)에 체적 탄성계수 K와 비체적의 관계식(K=vp/κ) 및 응력-변형률 관계식(dϵv=dp/K)를 대입하면 흡입력 변화에 의해 발생하는 평균 응력의 증분은 다음과 같이 표현할 수 있다.

한편, MCC 모델은 포화 상태를 가정하고 유효응력 개념만으로 응력–변형률 거동을 정의하므로 별도의 수리–역학 연동해석을 필요로 하지 않는다. 그러나 BBM은 흡입력을 독립 응력 변수로 취급하므로 흡입력의 계산이 필수적이다. 본 연구에서는 흡입력을 OGS와 같은 외부 코드를 통해 계산하여 FLAC3D로 전달받는 시나리오를 전제하여 코드를 구현하였으며, 이에 대한 구체적인 수리–역학 커플링 절차는 제3.2절에서 상세히 기술하였다.

개발된 BBM 해석모듈의 정확성을 검증하기 위하여 다양한 벤치마크 해석을 수행하였다. BBM 해석모듈이 MCC 모델의 확장 모델로서 포화 조건에서 MCC과 동일한 응력-변형률 거동을 재현할 수 있는지 검토하였으며, Alonso et al.(1990)의 연구에서 제시한 대표적인 응력 경로 시나리오(Cases 1, 2–1, 2‒2,3)를 기반으로 수치해석을 실시하여 해석 결과를 이론적 예측치와 비교하였다. 이 시나리오들은 LC 항복면과 SI 항복면이 실제 불포화 점토의 체적 거동 및 경로 의존성을 적절히 설명하는지를 검증하기 위하여 고안된 예제로, 각 시나리오의 목적과 특성은 다음과 같다.

∙ Case 1: 순평균응력 조건을 달리하면서 흡입력을 점진적으로 감소시키는 습윤(wetting) 시나리오로, 응력 경로에 따라 흡입력 감소가 체적 변형에 미치는 영향을 검증한다.

∙ Case 2: 외력(응력)과 흡입력 변화를 교차 적용하여 응력 경로에 따른 경로 의존성(path dependence)을 평가하는 시나리오로, 응력의 증가와 함께 (Case 2‒1) 흡입력 감소, (Case 2‒2) 흡입력 증가 조건을 교차 적용하고 각기 다른 체적 변형 메커니즘을 분석한다.

∙ Case 3: LC 곡면과 SI 곡면의 연계성(coupling)을 검증하기 위한 시나리오로, 건조–습윤(drying–wetting) 사이클 이후의 응력 증가 경로를 통해 소성 체적 변형의 축적이 항복면 이동과 강성 변화에 미치는 영향을 분석한다.

2.2.2 검증 1: 포화조건 등방 압축 하 재하-제하 반복실험 모델링

본 연구에서는 개발된 BBM 해석모듈을 이용하여 포화, 등방 압축(isotropic compression) 조건에서 반복적인 재하–제하(loading–unloading) 테스트를 수치적으로 수행하고, 해석 결과를 MCC 모델(FLAC3D 내장모델)의 거동과 비교하였다. BBM은 흡입력을 독립 변수로 도입함으로써 불포화 토양의 복잡한 거동을 설명할 수 있도록 확장된 구성모델이지만, 흡입력이 0인 포화 조건에서는 MCC 모델과 이론적으로 동일한 거동을 보여야 한다. 본 연구에서는 MCC 모델의 해석 결과를 기준(reference)으로 설정하고, 동일한 조건에서 흡입력을 0으로 설정한 BBM 해석을 수행함으로써, 구현된 BBM의 수치적 타당성을 검증하고자 하였다.

수치모델링은 단일 해석 요소(zone)를 대상으로 수행되었으며, 과압밀비(Over Consolidation Ratio, OCR)가 1.0인 정규압밀 상태(normally consolidated state)의 점토를 가정하였다. 초기응력은 σxx=σyy=σzz=-5.0MPa(압축)의 등방응력 상태로 설정하였고, 총 5단계의 반복적인 재하–제하 경로에 따른 평균 유효응력, 비체적, 체적 탄성계수의 변화 이력 등을 비교 분석하였다. 해석에 사용된 입력 물성값은 Table 1에 제시되어 있으며, 재하 및 제하 과정은 요소 경계면에 일정 속도를 부여하는 방식으로 수행되었다. MCC 모델과 BBM에 대하여 동일한 해석 조건과 하중 경로를 적용하였으며, BBM 해석에서는 흡입력을 0으로 고정하여 포화 상태에서 점토의 거동을 재현하도록 하였다. 두 모델은 모두 비선형 탄소성 모델로서 체적 탄성계수 K가 현 시점의 응력상태와 비체적에 따라 변하게 되며, K=vp/κ의 관계를 갖는다. Table 1에 제시된 Kmax는 체적탄성계수의 상한값으로 계산 중 K>Kmax인 조건이 발생하는 경우 오류 메시지를 출력하도록 하여 수치적 안정성을 확보하기 위해 설정된 값이다.

Fig. 2는 두 모델의 해석 결과를 비교한 그림으로 Fig. 2(a)는 순평균응력과 수직 변형률 간의 관계를 보여주며 Fig. 2(b)는 로그스케일의 순평균응력과 비체적의 관계를 보여준다. 그림에서 확인할 수 있듯이 반복적인 하중 경로 동안 포화 조건 BBM 해석 모듈의 결과는 MCC 모델과 정확히 일치하는 것으로 나타났다. 응력과 변형률, 비체적의 관계가 정확히 일치하며, 제하 구간에서 소성–탄성 영역 전이로 인해 체적 탄성계수가 증가하는 비선형 탄성 거동이 뚜렷하게 나타난다. 또한, 각 하중 사이클이 진행됨에 따라 정규 압밀선(𝜆)과 탄성 팽윤선(𝜅)을 따르는 전형적인 점토의 거동이 적절히 재현됨을 확인할 수 있었다.

Fig. 2.

Comparison of numerical results between the developed BBM (at zero suction) and the built-in MCC model in FLAC.

2.2.3 검증 2: Alonso의 Case 1

Case 1은 Alonso et al.(1990) 논문에 제시된 ‘Volumetric deformation induced by wetting at increasing confining stress’ 시나리오에 해당하며, 불포화 상태의 점토에 대해 다양한 순 평균 응력 조건에서 흡입력을 감소시키는 조건에서 체적의 변화를 분석하는 것이 목적이다. 이 시나리오는 낮은 평균 응력에서는 팽창, 높은 평균 응력에서는 붕괴(소성 압축)가 발생하는 불포화 점토의 전형적인 거동 특성을 재현하기 위하여 제시되었다.

Fig. 3은 Case 1에서 제시된 세 가지 응력 경로에 대한 이론 모델과 수치해석 결과를 비교한 것이다. Fig. 3(a)는 각 응력 경로에 따라 순평균응력과 비체적의 관계를 도식화한 것으로, 각 경로는 흡입력을 감소시키는 시점의 순평균응력 크기에 따라 Path 1(A‒B‒F), Path 2(A‒C‒D‒F), Path 3(A‒E‒F)으로 구분된다. 이들 경로는 모두 동일한 초기 상태점 A(p = 0.15 MPa, s = 0.20 MPa)에서 시작하여, 최종 상태점 F(p = 0.60 MPa, s=0)까지 진행되며, 각 경로는 흡입력 s를 0으로 감소시키는 시점에 따라 거동 특성이 달라진다. 이러한 경로 설정은 BBM에서 정의된 LC 항복면의 위치 변화를 통해, 흡입력 변화에 따른 체적 거동과 항복 조건의 경로 의존성을 재현하기 위해 구성된 것이다. 해석에서 사용한 입력 물성은 논문에 제시된 수치를 그대로 적용하였으며, Table 2에 정리되어 있다. 수치 모델링은 단일 해석 요소를 대상으로 수행되었으며, 초기 상태의 p(0.15 MPa)와 p0*(0.20 MPa)를 통해 확인할 수 있듯 대상 재료는 과압밀 상태(over consolidated state)의 초기 조건을 가정한다. 또한, 흡입력의 변화는 s0 이하에서 발생하므로 Case 1의 경우 SI 항복면의 영향은 배제된다.

Fig. 3.

Comparison between theoretical predictions and numerical results from the BBM for Case 1 scenario

Fig. 3(b)-(d)는 Alonso et al.(1990)의 이론에 따른 응력 경로와 본 연구의 수치해석 결과를 비교한 것이다. Path 1은 낮은 순평균응력 조건(p = 0.15 MPa)에서 흡입력을 0으로 감소시킨 후 평균 응력을 증가시키는 습윤–재하 경로이다. 초기 AB 구간에서는 흡입력 감소에 따라 탄성 영역 내에서 부피 팽창(elastic swell)이 발생하며, 이로 인해 비체적이 증가한다. 이후 응력이 증가하면서 재료는 포화 조건의 LC 항복면에 도달하게 되고 소성 수축(plastic shrinkage)이 유발된다. 수치해석 결과는 이러한 팽창–수축 거동을 정밀하게 재현하였으며, 특히 탄성–소성 전이 구간의 곡선 형태에서 이론 해석과 매우 우수한 일치를 보였다. Path 2는 p를 0.35 MPa까지 먼저 증가시킨 다음 흡입력을 감소시키고, 다시 응력을 증가시키는 재하‒습윤–재하 경로이다. 이 경우에는 p가 0.35 MPa에 이르기 전에 이미 LC 항복면에 도달하여 탄성‒소성 전이가 발생한다. 또한 C 지점에서 흡입력을 감소시키면 동시에 급격한 체적 수축이 유발되어 D 지점에 이르게 되며, 이후 F 지점에 이를 때까지 응력 증가에 따른 소성 수축이 지속적으로 발생한다. Path 1에 비해 더 높은 응력 수준에서 소성 거동이 시작되며, 강성 증가에 따라 소성 구간의 기울기는 상대적으로 완만하게 나타난다. 수치해석 결과는 붕괴 시의 급격한 체적 감소뿐만 아니라, 이후의 소성 거동에 이르기까지 이론 해석 곡선과 매우 높은 정확도로 일치함을 보였으며, 이는 BBM 해석 모듈이 LC 항복 조건 및 경화 거동(hardening behavior)을 정밀하게 구현하고 있음을 보여준다. Path 3은 가장 높은 순평균응력 상태(p = 0.30 MPa)에서 흡입력의 감소를 유도하는 재하‒습윤 경로이다. 이 경우 AC 구간까지는 Path 2와 동일한 거동을 보이지만, E 지점까지 응력 증가로 인한 체적 감소가 지속적으로 발생하게 되며, EF 구간에서 흡입력의 감소로 인한 추가적인 소성 수축이 발생한다. 이 경로에서도 수치 결과는 이론 해석과 매우 유사한 경향을 나타내며, BBM 해석 모듈이 고응력 상태에서도 탄소성 거동을 정확하게 재현하고 있음을 확인할 수 있다.

수치해석 시 각 경로 별로 A, C, E 지점에서 흡입력의 변화를 모사할 때, 비선형 응력–변형률 거동을 정밀하게 구현하기 위하여 흡입력 변화 구간을 총 100회의 점진적 증분(incremental step)으로 분할하고, 각 단계마다 평형 조건을 만족하도록 계산하였다. 이 과정에서 C–D 구간과 E–F 구간에서는 해석 요소 내부에서 소성 수축 변형률이 발생하므로, 내부 응력의 재분배(stress redistribution)가 일어나고 이에 따라 일정한 경계 응력 조건을 만족시키기 위한 미세한 응력 변동(stress variation)이 관찰된다. 본 해석에서는 해석 요소의 경계면을 따라 절점(node)에 일정한 외력(traction)을 부여하는 방식으로 일정 응력 경계 조건을 적용하였기 때문에, 이러한 응력 변동은 수치적 오차가 아닌 흡입력 증가에 따른 새로운 평형상태로의 점진적 이행을 반영하는 타당한 현상이다.

2.2.4 검증 3: Alonso의 Case 2

Case 2는 “Effect of Alternate Application of Load and Suction Changes” 시나리오로, 흡입력 변화와 하중 증가를 교차 적용했을 때 나타나는 응력 경로별 비체적 변화 특성을 규명하는 것을 목적으로 한다. 이 시나리오는 흡입력 변화 방향에 따라 흡입력을 감소시키는 경우(Case 2‒1, wetting)와 증가시키는 경우(Case 2‒2, drying)로 구분되며, 각 조건에서의 경로 의존적 탄성–소성 전이 거동과 체적 변형 메커니즘을 정량적으로 분석하는 데 중점을 둔다.

Fig. 4는 Case 2‒1의 세 가지 응력 경로에 대한 이론해와 수치해석 결과를 비교한 것이다. Fig. 4(a)는 세 응력 경로에 대한 이론적 순평균응력–비체적 관계를 보여주며, Fig. 4(b)-(d)에는 각 경로에 대한 수치해석 결과를 이론해와 함께 도시한 것이다. 세 응력 경로는 모두 초기 상태점 A(p = 0.15 MPa, s = 0.20 MPa)에서 시작하여 최종 상태점 F(p = 0.60 MPa, s = 0)에 도달한다. Path 1(A–B–F)은 흡입력 0.20 MPa를 유지한 채 순평균응력이 0.60 MPa가 될 때까지 하중을 가하고 BF 구간에서 흡입력을 0으로 감소시키는 재하‒습윤 경로이며, Path 2(A–C–D–F)는 습윤(s = 0.10 MPa)–재하–습윤(s = 0) 경로이다. Path 3(A–E–F)은 흡입력을 0으로 감소시킨 후 하중을 가하는 습윤‒재하 경로이다. 각 경로는 초기 흡입력 조건과 변동 폭이 달라 LC 항복면의 위치 변화, 탄성–소성 전이 시점, 소성 구간의 체적 변화 기울기에 차이를 보인다. Path 1은 높은 흡입력 조건에 의해 가장 큰 항복강도를 보이며 이에 따라 가장 높은 순평균응력 상태에서 탄성–소성 전이가 발생한다. 한편, Path 2와 Path 3은 초기 습윤 과정에서 흡입력 감소에 의한 탄성 팽창이 발생한다. 경로에 따라 체적의 감소 패턴은 다르지만 세 경로 모두 최종 지점 F에 도달했을 때 동일한 비체적에 수렴하게 된다. 흡입력이 낮은 상태에서 하중 재하가 이루어질수록 강도와 강성이 감소하게 되므로, 탄성-소성 전이가 더 낮은 응력 조건에서 발생하게 되고 소성 체적 변형이 더 뚜렷이 나타난다.

Fig. 5는 Case 2‒2의 세 가지 응력 경로에 대한 이론해와 수치해석 결과를 비교한 것이다. Fig. 5(a)는 세 응력 경로에 대한 이론적 순평균응력–비체적 관계를 보여주며, Fig. 5(b)‒(d)에는 각 경로에 대한 수치해석 결과를 이론해와 함께 도시한 것이다. Case 2‒2는 흡입력의 증가와 하중 재하가 교차 적용되는 시나리오로서, 모두 동일한 초기 상태점 A(p = 0.15 MPa, s = 0)에서 시작하여 최종 상태점 F(p = 0.60 MPa, s = 0.2 MPa)에 도달한다. Path 1(A–B–F)은 흡입력이 0인 포화 상태에서 순평균응력을 0.60 MPa까지 증가시킨 후, B‒F 구간에서 흡입력을 0.2 MPa로 증가시키는 건조‒재하 경로이다. Path 2(A–C–D–F)는 건조(s = 0.10 MPa)–재하–건조(s = 0.20 MPa) 경로를 따르며, Path 3(A–E–F)은 흡입력을 0.20 MPa까지 증가시킨 후 하중을 가하는 건조‒재하 경로이다. Path 1은 흡입력이 0인 포화 상태에서 하중이 가해지므로 가장 낮은 응력 수준에서 탄성–소성 상태의 전이가 일어난다. 항복 이후 B 지점까지는 응력 증가에 의한 소성 수축이 발생하지만, B‒F 구간에서는 흡입력 증가에 의한 탄성 수축이 발생하게 된다. Path 2에서는 첫 번째 건조 단계(A–C)에서 탄성 수축이 발생하고, LC 항복면의 이동으로 인해 재하 구간(C–D)에서 Path 1보다 더 높은 항복강도를 보인다. 또한, 두 번째 건조 단계(D–F)에서 소성이 아닌 탄성적인 수축이 발생하게 된다. Path 3에서는 높은 흡입력 증가량에 의해 초기 탄성 수축량이 상대적으로 크게 발생하게 되며, 강도와 강성이 증가한 상태에서 재하가 이루어지므로 탄성‒소성 전이 응력이 높고 비체적의 감소가 비교적 완만한 기울기로 발생하게 된다. Case 2‒1의 경우 세 응력 경로 모두 최종 상태점 F에서 동일한 비체적 값으로 수렴하지만, Case 2‒2에서는 서로 다른 비체적으로 수렴하게 되는데, 이는 습윤 과정이 소성 변형을 가속화하는 반면 건조 과정은 흡입력 증가를 유도하여 응력 경로를 탄성 영역 쪽으로 이동시키는 작용을 하기 때문이다. 따라서 건조 과정이 포함되는 경우, 시편이 동일한 응력 변화량과 흡입력 변화량을 경험한다 하더라도 이력에 따라 체적 변형 메커니즘이 확연히 달라지는 뚜렷한 경로 의존성을 나타내게 된다.

Fig. 4.

Comparison between theoretical predictions and numerical results for BBM under the Case 2-1 (wetting and loading) scenario

Fig. 5.

Comparison between theoretical predictions and numerical results for BBM under the Case 2-2 (drying and loading) scenario

Fig. 4와 Fig. 5의 결과를 종합적으로 살펴볼 때, 본 연구의 BBM 해석모듈은 Case 2‒1과 Case 2‒2의 모든 응력 경로에 대해 탄성–소성 전이 시점, 소성 구간 기울기 변화, 습윤에 의한 소성 붕괴, 건조에 의한 탄성 수축, 체적 변형의 경로 의존성을 매우 정밀하게 재현하였다.

2.2.5 검증 4: Alonso의 Case 3

Case 3은 “coupling of the SI and LC yield curves” 시나리오로, 흡입력의 변화가 항복 흡입력 s0이하에서 발생하는 Case 1, Case 2-1 및 Case 2-2와 달리, 초기의 s0보다 큰 흡입력이 작용하여 SI 항복면과 LC 항복면이 동시에 작용하는 조건을 다룬다. 이 경우에는 흡입력 변화에 따른 SI 항복면의 이동이 모델 예측에서 핵심적인 역할을 하며, 그 결과 LC 항복면의 위치와 형상도 함께 변화하게 된다.

Fig. 6은 Case 3에 대한 두 가지 응력 경로에 대한 이론 해석 결과와 수치해석 결과를 비교한 것이다. Fig. 6(a)는 이차원 응력 공간(p, s) 상에 표현된 두 경로와 각 경로의 순평균응력–비체적 관계를 함께 보여주며, Fig. 6(b)와 Fig. 6(c)는 각 경로별 수치해석 결과와 이론 해석의 상세 비교를 나타낸다. Path 1(A–B)은 흡입력이 0인 포화 조건에서 순평균응력을 0.15 MPa에서 0.6 MPa 증가시키는 단순 재하 경로이다. Path 2(A–C–D–E)는 먼저 흡입력을 0.3 MPa로 증가시키는 건조 구간(A–C)과 이를 다시 0 MPa로 감소시키는 습윤 구간(C–D)을 거친 뒤, 순평균응력을 증가시키는 재하 구간(D–E)으로 구성된다. 이때 A와 D 지점, B와 E 지점은 각각 동일한 순평균응력과 흡입력 값에 대응한다. Case 3에서 초기 s0가 0.025 MPa이므로, Path 2의 건조 구간(A–C)에서 재료는 SI 항복면에 도달하게 된다. 이때 탄성 수축과 소성 수축이 동시에 발생하여 비가역적인 체적 감소가 나타난다. 이어지는 습윤 구간(C–D) 구간에서는 흡입력 감소에 따른 탄성 팽창이 발생하지만, 이전 단계에서 발생한 소성 압축량은 회복되지 않는다. 따라서 D 지점 (p = 0.15 MPa, s = 0.0 MPa)에서 비체적은 동일한 응력 상태인 A 지점에 비해 작게 나타난다. Fig. 6(a)를 살펴보면, 이 과정에서 SI 항복면은 SIi에서 SIc로 상향 이동하며, 경화법칙에 따라 소성변형이 p0*의 증가를 유도하므로 LC 항복면 또한 초기 위치 LCi에서 LCD로 우측으로 이동하게 된다. 따라서 Path 2의 재하 구간(D–E)에서 탄성–소성 전이가 Path 1보다 높은 응력 수준에서 발생하게 된다. Path 2에서 흡입력을 변화시키는 A–C–D 구간에서 나타나는 응력 변동의 원인은 제 2.2.3절에서 설명한 바와 같다.

Fig. 6.

Comparison between theoretical predictions and numerical results for BBM under the Case 3 (SI and LC yield curve coupling) scenario

Fig. 6(b)와 6(c)에서 확인할 수 있듯이, 본 연구의 BBM 해석모듈은 건조-습윤 이력에 따른 비가역적 변형, 소성변형에 기인한 SI–LC 항복면의 커플링 및 이동, 그리고 최대 흡입력과 선행압밀응력의 변화를 합리적으로 모사하였다. 수치해석 결과는 이론 해석 곡선과 거의 일치하는 것으로 나타났으며, 이를 통해 본 해석모듈이 불포화 점토의 경로 의존적 거동과 체적 변형 메커니즘을 정확하게 예측할 수 있음을 확인하였다.

3. OGS-FLAC 연동해석을 통한 Sandwich Task Step0a 모델링

3.1 DECOVALEX-2027 Sandwich Task

DECOVALEX (DEvelopment of COupled models and their VALidation against EXperiments) 프로젝트는 1992년부터 시작된 국제 공동연구로, 고준위 방사성폐기물의 지층처분을 고려하는 각국의 방사성폐기물 관리·규제기관과 연구기관이 참여하여, 처분장 내 열-수리-역학-화학적(thermo-hydro-mechanical-chemical, THMC) 연계거동의 이해와 해석기법 개발을 목표로 한다(Birkholzer et al., 2019, DECOVALEX, 2025). DECOVALEX 프로젝트는 분담금을 납부하는 정식 회원기관(Funding Organization, FO)과 연구기관(Research Organization, RO)으로 구성되며, 9번째 단계인 DECOVLAEX-2027에는 프랑스 ANDRA, 독일 BASE/BGE/BGR, 중국과학원(CAS), 캐나다 CNSC/NWMO, 네덜란드 COVRA, 미국 DOE, 독일 DynaFrax, 스페인 ENRESA, 스위스 ENSI, 영국 NWS, 스웨덴 SSM, 체코 SURAO, 대만 Taipower, 그리고 한국지질자원연구원(KIGAM), 한국원자력연구원(KAERI) 등 총 18개 정식 회원기관이 참여하고 있다. 이들은 재정 분담과 함께 각국의 연구기관·대학과 협력하여, 현장시험 모사, 수치해석, 이론 모델 개발 및 성능 검증을 수행한다. 프로젝트는 4년 주기로 진행되며, 각 주기마다 새로운 연구 과제(Task)가 설정되고, 참여기관들은 각자의 수치해석 기법으로 동일한 실험 또는 현장시험 조건을 모사한 뒤 결과를 상호 비교·평가한다. 현재는 DECOVALEX‒2027 단계가 진행 중이며(2024‒2027년), 총 7개의 Task가 수행되고 있다. 한국지질자원연구원은 정식 회원기관으로서 Sandwich Task와 Safenet2 Task에 참여 중이며, 본 연구는 Sandwich Task에 참여하여 수행하고 있는 연구 내용의 일부를 소개한 것이다.

Sandwich Task는 스위스 Mont Terri Rock Laboratory에서 수행된 대규모 Sandwich 실험(Fig. 7)을 기반으로, 수직 유압 밀봉 시스템(Sandwich sealing system)의 수리–역학적(HM) 연계 거동을 정밀하게 모사하고 해석 코드를 검증하는 것을 주요 목표로 한다(Friedenberg and Hinze, 2024, Friedenberg et al., 2025). 이 시스템은 독일 Karlsruhe Institute of Technology (KIT)에서 제안된 개념으로(Nüesch et al., 2002, Emmerich et al., 2019), 벤토나이트 기반의 밀봉 세그먼트(DS: binary bentonite mixture segment)와 높은 투수성을 갖는 사질토 기반의 등전위 세그먼트(ES: equipotential sand segment)를 교차 배치함으로써, 수분의 균일한 침투와 벤토나이트의 고른 팽윤을 유도하는 것을 목표로 한다. 이러한 설계는 굴착손상대(excavation damaged zone, EDZ)에서 발생하는 비균질 팽윤과 우선 유동(preferential flow) 문제를 완화하여 밀봉 성능의 신뢰성을 높일 수 있는 장점이 있다. Sandwich 밀봉 시스템의 성능을 평가하고 수리-역학적 연계거동을 규명하기 위하여 다양한 규모의 실험이 수행되었다. MiniSandwich 및 Semi-technical scale 실험(HTV-7 실험)에서는 벤토나이트의 팽윤 거동, 투수성 저감, 지화학적 변화 등에 관한 장기 데이터가 축적되었으며(Emmerich et al., 2019, Emmerich et al., 2021), 현재 스위스 Mont Terri Rock Laboratory에서 대규모 현장 실험이 진행 중이며, 이를 통해 얻어진 데이터는 수치해석 모델 검증 및 물성 매개변수 보정에 활용되어 복잡한 수리-역학적 연계 거동 해석의 중요한 기초자료를 제공하고 있다.

Sandwich Task는 독일 Gesellschaft für Anlagen- und Reaktorsicherheit (GRS) gGmbH가 주도하고 있으며, Larissa Friedenberg와 Matthias Hinze가 Task lead로서 실험 데이터 관리와 모델링 지침을 총괄하고 있다. GRS는 팽윤압, 투수성, 수분함량 등 핵심 데이터를 조정·배포하고 있으며, 각 참여 연구팀은 이를 바탕으로 블라인드 예측 –수치해석–보정 과정을 반복 수행한다. 이를 통해 각자의 해석 코드를 비교·검증하는 공동연구가 이루어지고 있으며, Table 3은 Sandwich Task의 참여 연구팀과 적용 기법을 요약한 것이다.

Fig. 7.

Scheme of the Sandwich sealing system and the in-situ Sandwich experiment at Mont Terri (after Wieczorek et al., 2024)

Sandwich Task는 4년간의 연구 기간 동안 단계적 과정을 거쳐 벤토나이트–사질토 복합 밀봉 시스템의 수리–역학적(HM) 거동을 규명하고, 이를 모사할 수 있는 수치해석 기법을 개발·검증하는 것을 목표로 한다. Fig. 8은 Sandwich Task의 단계별 모델링 과정을 도식화한 것이다. 첫 번째 단계인 Step 0은 실험실 규모의 시뮬레이션으로, 벤토나이트의 물성 보정과 단순화된 조건에서의 HM 연계 거동을 분석하는 데 목적이 있다. 이 단계는 다시 세부 단계로 구분되며, Step 0a에서는 Calcigel 벤토나이트의 팽윤 압력 실험 결과를 활용하여 기본 물성을 보정한다. Step 0b에서는 MiniSandwich 실험을 모사하여 수화 과정에서의 팽윤, 수분 분포, 투수성 변화를 재현하며, Step 0c에서는 준대형(semi-technical) 규모의 HTV-7 실험을 대상으로 3차원 모델을 적용해 모래 렌즈와 굴착손상대(EDZ)를 포함한 복잡한 유체 흐름 조건을 분석한다. Step 1은 스위스 Mont Terri Rock Laboratory에서 수행된 Shaft 1 실험을 대상으로 하며, 벤토나이트–사질토 Sandwich 시스템과 모암인 Opalinus Clay 사이의 상호작용을 정량적으로 평가한다. Step 1a에서는 순수 수화 과정에 따른 압력 및 응력 변화를 모사하고, Step 1b에서는 굴착 및 환기 이력까지 고려하여 보다 실제적인 현장 조건을 재현한다. 이때 산출되는 액체 압력, 반경 응력, 축응력 등이 주요 비교 지표로 활용된다. 마지막으로 Step 2는 Mont Terri 현장의 전체 Sandwich 시스템(Shaft 1과 Shaft 2)을 대상으로 하는 대규모 현장 시뮬레이션 단계이다. 단순화된 3차원 모델을 적용하여 장기적인 수화 과정을 재현하고, 벤토나이트와 모암의 HM 상호작용을 종합적으로 평가한다. Step 2a에서는 기본적인 수화 이력을 고려하고, Step 2b에서는 발굴, 환기, 암반의 비선형 거동까지 반영한 고도화된 해석을 수행한다. 이를 통해 Sandwich 밀봉 시스템이 장기적으로 지하수 유입과 방사성 핵종 이동을 효과적으로 차단할 수 있는지를 검증하게 된다.

Fig. 8.

Modeling steps of the Sandwich Task (Friedenberg and Hinze, 2024)

3.2 OGS-FLAC 연동을 통한 수리-역학 커플링

앞서 제2.2.1절에서 언급한 바와 같이, FLAC3D의 수리유동 해석은 포화 상태의 유동을 전제로 하므로, BBM의 주요 응력 변수인 흡입력을 계산하기 위해서는 외부의 수리유동 해석 모듈이 필요하다. 본 연구에서는 OpenGeoSys (OGS)를 활용하여 BBM의 흡입력을 계산하고, 이를 FLAC3D와 연동하는 방식을 채택하였다. 이 접근법은 Park et al.(2019)의 연구에서 제안된 OGS–FLAC 해석기법으로, 열–수리–역학(thermal–hydrological–mechanical, THM) 복합 거동을 모사하기 위해 개발되었다. Fig. 9는 OGS-FLAC의 계산과정 및 연동방식을 보여주는 개념도이다(Kim et al., 2024). OGS-FLAC은 오픈소스 해석 코드인 OGS와 상용 해석 코드 FLAC3D를 연동하여, 검증된 두 코드의 상호보완적 기능을 통해 단일 코드로 다루기 어려운 THM 연계 거동을 정량적으로 평가할 수 있다는 장점이 있으며, CO₂ 저장 대수층의 수리–역학적 안정성 분석(Kim et al., 2021)과 Mont Terri 지하연구시설에서 수행된 FE 시험의 벤토나이트 완충재 및 암반의 THM 거동 해석(Kim et al., 2024)에 적용된 바 있다.

OGS–FLAC은 순차적 연동(sequential coupling) 해석 방식을 기반으로 하며, 매 시간 증분(time step) 또는 iteration 단계에서 각 코드가 순차적으로 독립 실행되고 결과를 교환하는 구조를 따른다. 구체적으로, OGS에서 계산된 간극수압과 온도 분포는 FLAC3D로 전달되어 변형 및 응력 해석에 반영되고, 이어서 FLAC3D에서 산출된 체적 변형률 및 변형 정보는 다시 OGS로 전달되어 다음 단계의 열–수리–유동 해석에 활용된다. 이러한 약결합(loose coupling) 접근은 모든 방정식을 단일 시스템으로 통합하여 동시에 해석하는 일체형 연계(monolithic coupling) 방식에 비해, 구성방정식 구현과 모델 개발에 필요한 비용과 시간을 절감할 수 있으며, 동시에 검증된 상용 코드와 오픈소스 코드를 효과적으로 활용할 수 있다는 장점을 지닌다.

Fig. 9.

OGS-FLAC simulation scheme based on a file-based sequential coupling with grid generation (Kim et al., 2024)

본 연구의 OGS 해석에서는 불포화대에서 지하수 유동 해석에 널리 사용되는 Richards 방정식을 지배방정식으로 적용하였다. Richards 방정식은 Darcy 법칙과 연속 방정식을 결합한 형태로, 포화대와 불포화대 모두에서 적용가능하며 다음과 같이 표현된다.

여기서 p는 유체압력, C(p)=∂θ/∂pθ는 압력에 대한 비체적수분용량(specific moisture capacity), k는 고유투수계수, 𝜇는 유체 점성계수, krl(S(p))는 상대투수도, 𝜌는 유체 밀도, g는 중력가속도, z 위치수두, θ=θ(p)는 체적수분함량을 의미한다.

불포화 조건에서는 모세관압(capillary pressure, Pc​)과 포화도(Sl) 간의 상관관계가 유동 특성을 결정한다. 모세관압은 기체압(pg)과 액체압(pl)의 차이로 정의된다. 불포화 조건에서는 기체압이 충분히 빠르게 소산되므로 기체 압력은 대기압으로 유지된다고 가정한다.

본 연구에서는 van Genuchten(van Genuchten, 1980) 경험식을 이용하여 모세관압-포화도 관계를 기술하였다. 주어진 모세관압에 대한 유효 포화도(Se)는 다음 식으로 표현된다.

여기서 Sl, Srl, Ssl은 각각 액체 포화도, 잔류 포화도, 최대포화도이며, P0는 기체 진입 압력(air entry pressure), m는 경험적 계수이다.

상대투수도(relative permeability, krl)는 Mualem–van Genuchten 관계식(Mualem, 1976, van Genuchten, 1980)에 의해 계산하였으며, 불포화 조건에서 포화도 감소에 따른 투수도의 저하는 다음과 같이 나타낼 수 있다.

OGS–FLAC 해석의 시간 증분(time step)은 OGS의 수리해석을 통해 결정되며, 각 단계에서 OGS의 해석 결과(포화도, 흡입력, 간극수압 등)가 FLAC3D에 전달된다. FLAC3D는 이를 바탕으로 모델 변수와 상태 변수를 갱신한 뒤 준정적 평형 해석(quasi-static equilibrium analysis)을 수행하여 응력과 변형을 계산한다. 다만, 본 연구는 흡입력 변화가 역학적 거동에 미치는 영향을 규명하는 데 주안점을 두었으므로, FLAC3D에서 계산된 역학적 결과는 다시 OGS 해석에 반영하지 않는 일방향 커플링(one-way coupling) 방식을 적용하였다.

3.3 Step 0a 모델링 및 결과 분석

앞서 언급한 바와 같이, Sandwich Task의 Step 0a는 Sandwich 실험에 사용된 벤토나이트 재료의 수리–역학적 거동을 재현하기 위한 구성모델 구현과 실험 기반 물성 보정을 목적으로 수행되었다. 이 단계의 주요 목표는 Calcigel 벤토나이트의 팽윤 및 수리적 응답을 수치모델로 신뢰성 있게 모사하기 위해 기본 물성을 조정하는 데 있으며, 이를 위해 다양한 건조밀도 조건에서 수행된 팽윤압 시험(swelling pressure test) 자료가 활용되었다. 총 19개의 실험 자료 중 건조밀도(ρdry)가 1,401 kg/m³, 1,500 kg/m³, 1,700 kg/m³일 때의 팽윤압 시험 결과가 기준 데이터로 선정되었으며, 시간 경과에 따른 팽윤압과 포화도(liquid saturation) 변화 데이터를 통해 모델 파라미터를 보정하고 연구팀 간의 결과를 비교 검증하였다. 각 시험편의 건조 밀도와 입자밀도(2,680 kg/m³, 105 〬C 건조 기준)에 기반하여 계산된 초기 공극률은 각각 0.477, 0.440, 0.366이며, 이를 BBM의 초기 비체적을 산정하는 주요 입력값으로 활용하였다.

팽윤압 시험은 직경 50 mm, 높이 30 mm의 원통형 시편을 대상으로 Oedometer 셀을 이용해 수행되었으며, Fig. 10은 팽윤압 시험의 개요와 수리–역학적 경계 조건을 나타낸 모식도이다. 시편은 건조 상태에서 하부로부터 0.2 bar의 과압 조건으로 수화되며, 시험 중 부피는 일정하게 유지된다. Step 0a에서 각 연구팀 간 비교·검증을 위해 요구되는 출력 결과는 (1) 포화도 증가에 따른 시편의 팽윤압, (2) 수분 포화도 변화, (3) 수직 변형률, (4) 주입 유체량이다.

Fig. 10.

Constant-volume swelling test setup and axisymmetric numerical model for Step 0a (Friedenberg and Hinze, 2024)

Fig. 11은 본 연구의 수치해석에 사용된 모델 형상과 경계조건을 나타낸 것으로, 축대칭 조건을 고려하여 실제 시험편의 1/6에 해당하는 30° 부채꼴 형상만을 모델로 구현하였다. 수리적 경계조건으로는 바닥을 제외한 모든 경계면에 대해 무유동(no-flow) 조건을 부여하였으며, 바닥에는 수압을 가하여 수화가 진행되도록 하였다. OGS의 Richards flow 모듈은 압력을 기반으로 편미분 방정식을 해석하므로, 초기 포화도 조건은 van Genuchten 관계식을 통해 상응하는 압력 조건으로 변환하여 적용하였고, 역학적 경계조건으로는 모든 경계면에서 수직방향 속도를 0으로 고정하여 일정 체적(constant volume) 조건의 팽윤압 시험을 재현하였다. 수치해석 결과의 해석 시, 팽윤압은 시험편의 수직방향 전응력 증가량의 평균을 통해 산정하였으며, 이 때 벽면 마찰이 없는 이상적인 조건을 가정하여 전 구간에서 수직방향 응력이 균일하다고 간주하였다. 수직 변형률은 시편 상하부 경계면에서 각각 5 mm 떨어진 두 지점에서 평가하였고, 유입수량은 초기 포화도 및 시간에 따른 포화도의 변화를 반영하여 계산하였다.

Fig. 11.

OGS-FLAC3D model geometry and applied boundary conditions for Step 0a swelling test simulation

각 건조밀도 조건별 해석 케이스에 대한 파라미터 보정은 수치모델이 해당 조건의 실내 실험 결과를 충실히 재현할 수 있을 때까지 반복적인 시행착오법(trial-and-error method)을 통해 수행하였다. 수리적 파라미터 보정은 고유투수계수(intrinsic permeability)를 조정하여 실험에서 측정된 액상 포화도의 시간 변화가 실험 결과와 일치하도록 하였으며, 역학적 파라미터 보정은 시간에 따른 팽윤압 측정 결과를 적절히 재현하는 BBM 파라미터 조합을 찾는 것을 목표로 진행하였다. 최종적으로 산출된 보정 물성치는 Table 4에 정리하였다.

Fig. 12는 P307-1/9 시편에 대한 수치해석 결과로, 시간 경과에 따른 흡입력, 순평균응력, 체적변형률의 공간 분포 변화를 보여준다. 시편 하부에서 수화가 시작되므로 흡입력의 감소는 공간적으로 순차적이고 불균질하게 나타난다. 이에 따라 초기 단계(20 hours)에서는 하부에서 구속조건과 부피 팽창에 따른 팽윤압 발생과 응력 증가가 먼저 나타나는 반면, 상부는 여전히 높은 흡입력이 유지되어 응력 상태와 변형이 공간적으로 불균질한 양상을 보인다. 시간이 경과함에 따라(222 hours) 수분 확산이 상부로 진행되면서 흡입력이 점차 소산되고, 그에 따라 순평균응력과 체적변형률도 하부에서 상부로 전달되어 비교적균질한 분포로 전환된다. 약 500시간 경과 시(504 hours)에는 시편 전체가 완전 포화 상태에 도달하여 흡입력이 0으로 수렴하였으며, 이 시점부터 더 이상의 새로운 팽윤압이 발생하지 않는다. 동시에 순평균응력과 체적변형률 분포 역시 전 영역에서 안정화되는 양상을 보인다. 또한, 순평균응력과 체적변형률의 분포가 유사한 경향을 나타내는 점은 주어진 팽창압 조건에서 시편이 탄성 영역 내에서 거동하고 있음을 의미하며, 이에 따라 비가역적인 소성 변형은 발생하지 않은 것으로 판단할 수 있다.

Fig. 12.

Simulation results of P307-1/9: temporal evolution of suction, net mean stress, and volumetric strain distributions

Fig. 13은 수치해석에서 도출된 포화도 결과를 실험 측정치와 비교한 것이다. 세 가지 건조밀도 조건에서 수치해석 결과는 측정치와 전반적으로 잘 일치하였으며, 특히 초기 수화 단계의 급격한 포화도 증가와 장기적으로 안정화되는 경향이 재현되었다. 여기에서는 복잡한 양방향 연계 해석을 적용하지 않았음에도, 투수계수의 보정만으로 초기 포화도 증가율을 적절히 설명할 수 있었으며, 이는 수리 물성 보정의 타당성을 보여준다. 한편, Fig. 14는 시뮬레이션을 통해 예측된 주입수량을 나타낸 것으로 약 500시간 이내에 대부분의 수분 유입이 이루어지고 이후에는 일정한 값으로 수렴하여 시편이 충분히 포화되었음을 확인할 수 있다. 또한 건조밀도가 높을수록 유효 공극률이 작아 최종 유입수량이 감소하는 경향이 관찰되었다.

Fig. 13.

Comparison between experimental measurements and numerical simulation results of degree of saturation for three specimens (P307-1/5, P307-1/9, and P307-1/18)

Fig. 14.

Numerical simulation results of cumulative injected water volume for three specimens (P307-1/5, P307-1/9, and P307-1/18)

Fig. 15는 수치해석에서 도출된 팽윤압(축응력 평균) 결과를 실험 측정치와 비교한 것이다. 세 가지 건조밀도 조건 모두에서 해석 결과는 측정치와 전반적으로 잘 일치하였으며, 초기 약 500시간 이내에 급격한 팽윤압 상승이 발생하고 이후 일정한 값에 수렴하는 경향을 재현하였다. 건조밀도가 높을수록 최종 팽윤압이 크게 나타났으며, 이는 실험 결과와 동일한 경향을 보였다. 다만 초기 단계에서는 수치해석과 실험 간 차이가 나타났는데, 이는 실험의 경우, 초기부터 상대적으로 높은 응력상태가 계측된 반면, 시뮬레이션에서는 일괄적으로 초기 조건을 0.1 MPa의 순평균응력으로 설정하였고, 수화가 시편 하부에서 시작되므로 초기에는 팽윤압 증가가 완만하게 진행되다가, 전반적인 포화도가 높아진 이후에는 급격히 증가하는 양상이 나타났기 때문이다.

Fig. 16은 동일한 조건에서 계산된 시편 상·하부의 축방향 변형률 분포를 나타낸 것이다. 초기에는 수화로 인한 팽윤압 발생에 따라 변형률이 증가하였으나, 포화도가 증가하면서 흡입력과 팽윤압이 감소함에 따라 변형이 점차 복원되는 경향을 보였다. P307- 1/9와 P307-1/18 시편에서는 변형률이 거의 완전히 회복되어 탄성 거동을 보인 것으로 판단된다. 반면 P307-1/5 시편에서는 포화도 수렴 이후 상부에서 양(+)의 변형률, 하부에서 음(–)의 변형률이 나타나는 수축 경향이 관찰되었으며, 이는 국부적인 파괴와 이로 인한 시료의 소성 수축이 발생했음을 시사한다. 실험에서는 변형률 분포가 직접적으로 계측되지 않아 정량적 비교는 어려우나, 수치해석 결과에 따르면 건조밀도가 높은 시편일수록 전체 변형의 크기가 상대적으로 작게 나타났다.

Fig. 15.

Comparison of experimental measurements and numerical simulation results of swelling pressure (axial stress) for three specimens (P307-1/5, P307-1/9, and P307-1/18)

Fig. 16.

Numerical simulation results of axial strains at the top and bottom of the specimens (P307-1/5, P307-1/9, and P307-1/18)

4. 요약 및 결론

본 연구에서는 방사성폐기물 심층처분장의 공학적 방벽 재료로 널리 활용되는 벤토나이트의 장기적 거동을 정량적으로 규명하기 위해, BBM의 구성식을 FLAC3D에서 구현하고, OGS와 연계한 수리–역학적 연동 해석기법을 개발하였다. 구현된 BBM 모듈은 먼저 포화 조건에서 MCC 모델과의 등가성을 검증함으로써 기본적인 이론적 정합성을 확보하였다. 또한 Alonso et al.(1990)이 제시한 대표적 벤치마크 시나리오(Case 1: 흡입력 감소에 따른 팽창/붕괴, Case 2–1/2–2: 하중–흡입력 교차 경로, Case 3: SI–LC 항복면 커플링)를 재현하여, 불포화 토양의 전형적 거동 특성과 경로 의존성을 충실히 모사할 수 있음을 입증하였다. 특히 흡입력 변화에 따른 체적 팽창과 수축, 건조–습윤 이력에 따른 소성 변형의 축적, 그리고 항복면 이동을 통한 강성 변화가 합리적인 수준에서 구현되었음을 확인하였다. 한편, DECOVALEX-2027 Sandwich Task의 Step 0a 팽윤압 시험 조건에 본 해석기법을 적용하여 Calcigel 벤토나이트의 세 가지 건조밀도(1,401, 1,500, 1,700 kg/m³)에 대해 팽윤압과 포화도 변화를 재현하였으며, 이를 위해 수리적·역학적 물성치를 보정하였다. 수치해석 결과는 실험 자료와 높은 수준의 일치를 보였으며, 이를 통해 본 해석모듈이 실제 실험을 정량적으로 설명할 수 있음을 확인하였다.

본 연구에서 개발한 BBM 기반 OGS–FLAC 연계 해석기법은 벤토나이트의 경로 의존적 변형 메커니즘, 불포화 상태에서의 수리–역학적 상호작용, 장기적 팽윤 및 투수성 저감 과정을 단계적으로 모사할 수 있다는 점에서 학문적·실용적 의의가 크다. 본 연구에서 확보한 물성치와 해석기법을 바탕으로, MiniSandwich, HTV–7, Mont Terri URL 등 보다 대형·복잡한 현장실험 조건에 적용하여 해석기법의 확장성과 예측력을 추가적으로 검증할 필요가 있다. 나아가 열적 영향(thermal effect)과 암반–밀봉재 경계에서 발생하는 비선형 상호작용까지 고려한 THM 통합 해석으로 발전시킬 경우, 실제 처분장 환경에서 다중방벽 성능평가와 안전성 확보에 기여할 수 있을 것으로 기대된다.

Fig. 17은 국제 공동연구에 참여한 여러 연구팀 간의 수치해석 결과를 실험 자료와 비교한 것으로 상단은 전응력의 변화를, 하단은 포화도의 변화를 나타낸다(Friedenberg et al., 2025). 이는 DECOVALEX-2027 제3차 워크숍에서 발표된 결과로 본 연구에서 제안한 모델의 결과는 ‘KIGAM’으로 표시되어 있다. 비록 본 연구에서 제시한 최종 보정 물성과 일부 다른 조건으로 계산된 결과이나, 본 연구의 모델이 벤토나이트 거동을 재현함에 있어 국내외 유수 참여팀들의 해석 코드와 유사한 수준의 성능을 확보하고 있음을 보여준다. 본 연구는 DECOVALEX-2027 Sandwich Task에 참여하여 국내외 기관들과의 협력 및 교류를 통해 지속적으로 개선, 확장할 예정이다.

Fig. 17.

Comparison of simulation results from different participants with experimental data for three tests: (a) P307-1/5, (b) P307-1/9, and (c) P307-1/18. The upper row shows the evolution of total stress, and the lower row shows the degree of saturation (Friedenberg et al., 2025)