Verification of Thermo-Hydro-Mechanical Coupling in the Sequential OGS-FLAC3D Simulator with a Fixed-Stress Split Stabilization

Taehyun Kim; Chan-Hee Park; Changsoo Lee; Jin-Seop Kim

doi:10.7474/TUS.2026.36.5.545

Preview

Original Article

Tunnel and Underground Space. 31 October 2026. 545-562
https://doi.org/10.7474/TUS.2026.36.5.545

Verification of Thermo-Hydro-Mechanical Coupling in the Sequential OGS-FLAC3D Simulator with a Fixed-Stress Split Stabilization

고정응력 분할 안정화를 적용한 OGS-FLAC3D 순차 연계해석의 열-수리-역학 복합거동 검증

Taehyun Kim¹

Chan-Hee Park²^*

Changsoo Lee³

Jin-Seop Kim³

김 태현¹

박 찬희²^*

이 창수³

김 진섭³

¹Senior Researcher, Korea Atomic Energy Research Institute

²Principal Researcher, Korea Institute of Geoscience and Mineral Resources

³Principal Researcher, Korea Atomic Energy Research Institute

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

²한국지질자원연구원 책임연구원

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

^{*Corresponding Author}

License (open-access, http://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 (http://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

Reliable prediction of coupled thermo-hydro-mechanical (THM) behavior induced by decay heat is essential for the safety assessment of deep geological repositories for high-level radioactive waste. OGS-FLAC3D, a sequential coupling simulator that combines OpenGeoSys for thermo-hydraulic analysis with FLAC3D equipped with various geomechanical constitutive models, offers the practical advantages of a modular framework; however, the convergence and stability of sequential iterations may deteriorate under strongly coupled conditions. In this study, a fixed-stress split stabilization term was introduced into the fluid mass balance equation, and the convergence and accuracy of the coupled analysis were verified against the analytical solution of the Booker-Savvidou point heat source problem. The numerical results for temperature, pore pressure, displacement, and stress showed good agreement with the analytical solution, with peak relative errors in pore pressure of −2.73% and −1.60% at two monitoring points located at different distances from the heat source. In the unstabilized fixed-strain split analysis, in which a small specific storage coefficient was introduced, the prescribed convergence criterion was not satisfied within the given number of coupling iterations, and the pore pressure was overestimated by approximately 67% relative to the analytical solution. This result demonstrates that fixed-stress stabilization plays an important role in achieving convergence within a practical number of coupling iterations for the strongly coupled benchmark considered in this study. For the stabilized cases, the final converged solutions showed only negligible differences among the tested stabilization coefficients and the T-H-M and H-T-M coupling sequences. These results confirm that the fixed-stress- stabilized OGS-FLAC3D framework can reliably reproduce the coupled THM response of the present benchmark problem.

Keywords

OGS-FLAC3D

Thermo-hydro-mechanical coupling

Fixed-stress split

Booker- Savvidou problem

Sequential coupling

고준위 방사성폐기물 심층처분장의 안전성 평가에는 붕괴열에 의한 열-수리-역학(THM) 복합거동의 신뢰성 있는 예측이 필수적이다. 열-수리 해석에 강점을 지닌 OpenGeoSys와 다양한 지반역학 구성모델을 갖춘 FLAC3D를 결합한 순차 연계해석 시뮬레이터(OGS-FLAC3D)는 모듈성의 장점을 지니나, 물리 과정 간 결합이 강한 조건에서는 순차 반복의 수렴성과 안정성이 저하될 수 있다. 본 연구에서는 이러한 한계를 해소하기 위해 고정응력 분할(fixed-stress split) 기반의 안정화항을 유체 방정식에 도입하고, Booker-Savvidou 점열원 문제의 해석해를 이용하여 연계해석의 수렴성과 정확성을 검증하였다. 수치해석 결과 온도· 간극수압·변위·응력이 해석해와 잘 일치하였으며, 간극수압의 첨두 상대오차는 점열원으로부터 거리가 다른 두 계측점에서 각각 −2.73%, −1.60%로 나타났다. 안정화항을 제거하고 작은 비저류계수를 부여한 고정변형 분할 해석에서는 주어진 반복 횟수 내에서 수렴 기준에 도달하지 못하고 간극수압이 해석해보다 약 67% 크게 나타났다. 이는 본 강결합 벤치마크에서 고정응력 분할 안정화가 실용적인 반복 횟수 내에서 연계 수렴을 확보하는 데 중요한 역할을 함을 보여준다. 또한 안정화항을 적용한 경우, 본 벤치마크 조건에서 시험한 안정화항 계수와 T-H-M 및 H-T-M 연계 순서에 따른 최종 수렴해의 차이는 매우 작게 나타났다. 이상의 결과로부터 고정응력 분할 안정화를 적용한 OGS-FLAC3D가 본 벤치마크의 THM 복합거동을 안정적으로 재현할 수 있음을 확인하였다.

키워드

OGS-FLAC3D

열-수리-역학 연계해석

고정응력 분할

Booker-Savvidou 문제

순차 연계해석

MAIN

1. 서 론
2. 지배방정식 및 커플링 기법
2.1 지배방정식
2.2 순차 커플링 구조
2.3 고정응력 분할(fixed-stress split) 안정화
2.4 요약
3. 검증: Booker-Savvidou 점열원 문제
3.1 문제 정의
3.2 해석 결과
3.3 안정화항 적용 여부에 따른 비교
3.4 안정화항 계수에 따른 영향
3.5 커플링 순서에 따른 비교
4. 고 찰
5. 결 론

1. 서 론

고준위 방사성폐기물 심층처분장에서는 폐기물에서 발생하는 붕괴열에 의해 주변 암반의 온도가 상승하고, 이로 인한 간극수압과 응력의 변화가 장기간에 걸쳐 발생한다. 이러한 열-수리-역학(thermo-hydro-mechanical, THM) 복합거동은 처분장의 안전성 평가에 핵심적인 요소이므로, 이를 신뢰성 있게 예측하기 위한 수치해석 기법이 요구된다. THM 복합거동은 온도·간극수압·변형이 서로 강하게 연동되므로, 각 물리 과정을 정확히 결합하는 커플링 기법이 해석의 정확성과 안정성을 좌우한다.

THM 커플링 해석은 크게 모든 지배방정식을 하나의 연립방정식으로 동시에 푸는 일체형(monolithic) 방식과, 각 물리 과정을 개별 해석코드로 번갈아 푸는 순차(sequential) 방식으로 구분된다. 일체형 방식을 적용하는 프로그램으로는 CODE_BRIGHT, OpenGeoSys 등이 있으며(CODE_BRIGHT, 2004; Grunwald et al., 2022), 순차 방식을 적용하는 프로그램으로는 TOUGH- FLAC, OGS-FLAC3D 등이 있다(Rutqvist et al., 2011; Park et al., 2019). 이 중 순차 방식은 검증된 기존 해석 코드를 그대로 재사용할 수 있고 각 물리 과정을 독립적으로 다룰 수 있는 모듈성의 장점을 가진다. 다만 이러한 모듈성은 주로 검증된 기존 솔버(solver)를 재사용하고 서로 다른 물리 모듈을 조합할 수 있는 구현·유지관리 측면의 실용적 장점이며, 강결합 문제를 수치적으로 다루는 측면에서는 일체형과 다른 특성을 가진다. 순차 방식에서는 커플링 지연(coupling lag)이나 반복 수 증가가 나타날 수 있으며, 본 연구에서 도입한 고정응력 분할 안정화는 이러한 순차 방식의 수렴 안정성을 보완하기 위한 것이다.

이러한 장점에 따라, 유한요소법(finite element method) 기반으로 열-수리 해석에 강점을 지닌 OpenGeoSys(Kolditz et al., 2012)와 유한차분법(finite difference method) 기반으로 다양한 지반역학 구성모델을 갖춘 상용 코드인 FLAC3D(Itasca, 2012)를 결합한 순차 연계해석 시뮬레이터(OGS-FLAC3D)가 개발되어(Park et al., 2019, 2020; Kim et al., 2021a), 단상 수리-역학 해석에서 다상 유동을 포함한 열-수리-역학(THM) 복합거동 및 실규모 현장 실험 해석으로 확장·적용되어 왔다(Kim et al., 2021b, 2022, 2024; Park and Park, 2022).

그러나 순차 방식은 각 물리 과정을 분할하여 풀기 때문에, 물리 과정 간 결합이 강한 강결합(strong coupling) 조건에서 반복의 수렴이 저하되거나 해의 정확도가 제한되는 한계를 가진다. 특히 유체와 고체 알갱이가 비압축에 가까워 비저류계수가 작은 조건에서는, 역학적 변형이 간극수압에 미치는 영향이 커져 순차 반복이 불안정해지기 쉽다. 비저류계수가 작을수록 압력 변화에 대한 직접적인 저류 효과가 감소하므로, 현재 압력 변화에 대응하는 골격 변형이 직전 값으로 늦게 반영되는 순차 연계의 지연 효과가 상대적으로 커지게 된다. 이러한 한계를 해소하기 위해 유체-역학 연계해석 분야에서는 고정응력 분할 기법(fixed-stress split)과 같은 안정화 기법이 제안되어 왔다(Kim et al., 2011; Mikelić and Wheeler, 2013). 고정응력 분할 기법은 총 평균응력이 고정된다는 가정으로부터 유도되는 안정화항을 유체 방정식에 도입하여, 비저류계수와 무관하게 순차 반복의 안정성을 확보한다.

본 연구에서는 고정응력 분할 기법 기반의 안정화항을 OGS-FLAC3D THM 커플링에 도입하고, 그 정확성과 안정성을 점열원(point heat source) 문제의 해석해로 검증하고자 하였다. 검증 문제로는 포화 다공질 매질 내 점열원에 의한 THM 거동의 해석해가 존재하는 Booker-Savvidou 문제(Booker and Savvidou, 1985; Chaudhry et al., 2019)를 채택하였다. 안정화항의 유도와 순차 커플링 구조를 제시하고, 온도·간극수압· 변위·응력에 대해 해석해와 비교하였으며, 안정화항의 적용 여부와 안정화항 계수, 커플링 순서가 해석 결과에 미치는 영향을 분석하였다. 이를 통해 순차 방식의 강결합 조건에서의 적용성을 검토하였다.

2. 지배방정식 및 커플링 기법

2.1 지배방정식

본 연구에서는 포화 다공질 매질에서 점열원에 의한 열-수리-역학(THM) 복합거동을 선형 열탄성 압밀 이론에 기반하여 해석하였다. 미소변형과 상수 물성을 가정하며, 지배방정식은 온도, 간극수압, 골격 변위에 대한 세 방정식으로 구성된다.

먼저 열 전달은 열전도가 지배하는 것으로 가정하며, 에너지 보존 방정식에 기반하여 다음 식과 같이 표현된다.

(1)

{(ρ_{c})}_{e f f} \frac{\partial T}{\partial t} - \nabla \cdot (λ \nabla T) = Q

위 식에서 ${(ρ_{c})}_{e f f}$ 는 매질의 유효 체적열용량(volumetric heat capacity, J/m³·K), $λ$ 는 유효 열전도도(effective thermal conductivity, W/m·K), $Q$ 는 점열원(point heat source, W)을 의미한다. 각 항을 살펴보면 ${(ρ_{c})}_{e f f} \frac{\partial T}{\partial t}$ 는 매질의 열저장(storage), $\nabla \cdot (λ \nabla T)$ 는 열전도(conduction) 항에 해당한다. OGS-FLAC3D은 열-수리 간 양방향 결합 해석이 가능하나, 본 연구의 Booker-Savvidou 문제 (Booker and Savvidou, 1985)에서는 유체 이류에 의한 열 이송을 무시하고 열전도도 및 열용량을 상수로 가정한다. 이 경우 열 방정식은 압력 및 유동에 의존하지 않는 닫힌 형태가 되므로, 본 문제에 한하여 온도장은 각 시간 단계에서 독립적으로 결정된다.

다음으로 유체의 거동은 질량보존 방정식에 기반하며, 골격의 변형과 온도 변화가 간극수압에 미치는 영향을 함께 고려하여 다음 식과 같이 표현된다.

(2)

S \cdot \frac{\partial p}{\partial t} + α \frac{\partial ε_{v}}{\partial t} - α_{u} \cdot \frac{\partial T}{\partial t} + \nabla \cdot v = 0, v = - \frac{k}{μ} \nabla p

위 식에서 $S$ 는 비저류계수(specific storage coefficient, Pa^-1)로 $S = ϕ / K_{f} + (α - ϕ) / K_{s} = 1 / M$ 로 정의되며, $ϕ$ 는 공극률, $K_{f}$ 와 $K_{s}$ 는 각각 유체와 고체 알갱이의 체적탄성계수(Pa), $M$ 은 Biot 탄성계수를 의미한다. $α$ 는 Biot 상수, $ε_{v}$ 는 골격의 체적변형률, $p$ 는 간극수압(Pa), $T$ 는 온도(K)이다. v는 Darcy 유속(specific discharge, m/s)으로 투수계수 k(m²)와 유체의 점성계수 $μ$ (Pa·s)에 의해 결정된다. $α_{u}$ 는 유체와 고체의 열팽창 차이에서 비롯되는 열-수리 연동계수(thermal expansion parameter)로 다음과 같이 주어진다.

(3)

α_{u} = α_{f} \cdot ϕ + (α + ϕ) \cdot 3 α_{s}

여기서 $α_{f}$ 는 유체의 체적열팽창계수, $α_{s}$ 는 고체 알갱이의 선형열팽창계수(1/K)를 의미한다. 식 (2)의 각 항을 살펴보면 $S \cdot \frac{\partial p}{\partial t}$ 는 유체 자체의 저류(storage), $α \cdot \frac{\partial ε_{v}}{\partial t}$ 는 골격의 체적변형이 간극수압에 미치는 영향, $α_{u} \cdot \frac{\partial T}{\partial t}$ 는 가열에 따른 열팽창(thermal pressurization), ∇·v는 유체 흐름(flux)에 해당한다. 이 중 $α \frac{\partial ε_{v}}{\partial t}$ 항은 역학 해석 결과를 필요로 하는 역학적 연동항이다. Park and Park(2022)은 이 항을 매질의 변형에 의한 소스/싱크(source/sink) 항으로 OGS-FLAC3D 해석 시뮬레이터에 처음으로 반영하였으며, Kim et al.(2021b, 2022, 2024) 또한 동일한 방식을 적용한 바 있다.

마지막으로 역학 해석을 위해서는 운동량 보존 방정식(linear momentum balance equation)을 적용하며, 간극수압과 열팽창에 의한 응력 변화를 함께 고려한다.

(4)

\nabla \cdot σ + ρ g = 0

(5)

σ = σ' - α pI, σ = D : (ε - α_{s} ∆ T m)

위 식에서 $σ$ 는 전응력(total stress), $σ'$ 는 유효응력(effective stress), $ρ$ 는 매질의 밀도, $g$ 는 중력가속도를 의미한다. D는 탄성계수 텐서, $ε$ 는 변형률, $∆ T = T - T_{o}$ 는 응력이 없는 기준 상태에 대한 온도 변화이다. m = [1, 1, 1, 0, 0, 0]ᵀ은 Voigt 표기의 Kronecker delta로, 열변형이 수직 성분에만 등방적으로 작용함을 나타낸다. 유효응력은 총응력에서 간극수압 항( $α pI$ )을 제외한 값으로 정의되며(Biot의 유효응력), 총변형( $ε$ )에서 열변형( $α_{s} ∆ T m$ )을 차감한 뒤 탄성계수를 곱하여 계산된다. 본 연구에서는 중력을 무시하였다.

Biot 상수는 배수 체적탄성계수 $K_{d r}$ 와 고체 알갱이 체적탄성계수 $K_{s}$ 에 의해 다음과 같이 정의된다.

(6)

σ = 1 - \frac{K_{d r}}{K_{s}}

본 연구의 대상 문제는 고체 알갱이가 비압축인 조건( $K_{s}$ → ∞)에 해당하여 α = 1이 된다. 나아가 Booker-Savvidou 문제는 유체와 고체 알갱이를 모두 비압축으로 가정하므로( $K_{f}$ , $K_{s}$ → ∞), 비저류계수는 0이 된다. 이 경우 압력 변화에 대한 직접적인 저류 효과가 사라져 유체 질량보존이 골격의 체적변형과 더욱 강하게 연동된다. 따라서 순차 연계에서는 수리-역학 사이의 반복 지연에 민감한 강결합 조건(strong coupling)이 된다.

2.2 순차 커플링 구조

OGS-FLAC3D은 열-수리 해석을 담당하는 OpenGeoSys와 역학 해석을 담당하는 FLAC3D를 파일 기반 순차적 연계해석 방식(file-based sequential coupling method)으로 결합한 해석 시뮬레이터이다. OpenGeoSys와 FLAC3D는 동일한 격자 형상을 공유하며 절점과 요소가 서로 대응된다. 데이터 교환은 파일 기반으로 이루어지며, OpenGeoSys는 절점에서 계산한 온도와 간극수압을 FLAC3D로 전달하고, FLAC3D는 요소에서 계산한 체적변형률을 OpenGeoSys로 전달한다. OpenGeoSys는 전달받은 체적변형률 증분을 요소 단위의 소스/싱크 항으로 반영하며, 안정화항 계산에 필요한 절점 간극수압은 형상함수로 요소 내에서 보간하여 사용한다.

Park et al.(2019) 및 Kim et al.(2021a)이 단상 유체 조건의 수리-역학(H-M) 연계해석을 대상으로 처음 개발하였으며, 이후 Kim et al.(2021b)이 다상 유동을 포함한 열-수리-역학(THM) 복합거동으로 확장하였다. Park and Park(2022)은 역학적 연동항을 변형률 기반 소스/싱크 항으로 도입하여 단상 유동 조건의 Booker- Savvidou 벤치마크를 검증하였고, Kim et al.(2021b, 2022, 2024)은 이를 다상 유동 해석에 적용하였다. 다만 이들 선행 연구는 안정화(보상)항 없이 변형률 소스텀만을 적용하였으며, 본 연구는 여기에 고정응력 분할 기반의 안정화항을 도입하고 커플링 구조를 개편하였다.

본 연구에서 각 시간 단계의 해석 순서는 T → H → M이며, 이는 Booker-Savvidou 문제의 경우 점열원에 의한 온도 변화가 간극수압을 유발하고, 이로 인한 유효응력 변화가 역학적 변형을 유발한다는 점을 고려하였다. 이러한 물리적 인과 순서를 반영하여, 선행 OGS-FLAC3D 연구들이 H-T-M 순서를 적용한 것과 달리 본 연구에서는 T를 먼저 해석하는 T-H-M 순서를 채택하였다. 또한 선행 연구(Park et al., 2019; Kim et al., 2021a, 2021b; Kim et al., 2022)에서 역학(M) 해석을 시간 단계 내 최종 단계에서 수행한 것과 달리, 본 연구에서는 T-H-M 전체를 매 반복(iteration)마다 수행하도록 구성하였다. 이에 따라 FLAC3D는 매 반복마다 호출되며, 반복은 간극수압이 충분히 수렴하도록 수행된 후 다음 시간 단계로 진행된다. 해석 순서와 각 프로그램 간 데이터 교환 시점은 Fig. 1과 같으며, OpenGeoSys에서 FLAC3D로는 간극수압과 온도가, FLAC3D에서 OpenGeoSys로는 체적변형률이 전달된다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-05/N0120360506/images/ksrm_2026_365_545_F1.jpg

Fig. 1.

Conceptual diagram of the sequential T-H-M coupling procedure in OGS-FLAC3D

2.3 고정응력 분할(fixed-stress split) 안정화

2.3.1 순차 분할의 안정성

순차 분할에서는 유체식 (2)의 역학적 연동항 $α \frac{\partial ε_{v}}{\partial t}$ 가 역학해를 요구하므로, 수리 해석을 푸는 시점에 미지수이다. 이 항을 처리하는 방식에 따라 순차 분할은 고정변형 분할(fixed-strain split)과 고정응력 분할로 구분되며, 두 방식 모두 유동 문제를 먼저 푼다(Kim et al., 2011). 본 절에서는 안정화항이 없는 고정변형 분할을 통해 안정성 문제를 살펴보고, 다음 절에서 고정응력 분할을 유도한다. 또한 본 절에서 다루는 순차 반복은 수리·열(TH) 해석과 역학(M) 해석 사이의 커플링 반복을 의미하며, 각 해석 프로그램 내부의 비선형 반복과는 구별된다.

이후의 반복식은 직전 시간 단계(time step) n에서 다음 시간 단계 n+1로 전진하는 한 시간 단계를 대상으로 하며, 두 종류의 인덱스를 사용한다. n은 시간 단계 인덱스로 $p^{n}$ 은 이미 수렴하여 확정된 직전 단계의 값을, $p^{n + 1}$ 은 이번에 구하려는 값을 의미한다. k는 한 시간 단계 내부의 반복 (iteration) 인덱스로 $p^{k}$ , $p^{k + 1}$ 은 반복 과정의 중간값이며, 반복이 수렴하면 $p^{n + 1}$ 이 된다. 한편 온도 T는 매 시간 단계에서 가장 먼저 해석되므로(T-H-M), 반복 인덱스 없이 n+1 단계에서 이미 갱신된 값 $T^{n + 1}$ 으로 반영된다.

고정변형 분할은 유동 문제를 풀 때 체적변형률의 변화율 $\frac{\partial ε_{v}}{\partial t}$ 를 미지수로 갱신하지 않고 직전 반복값 $ε_{v}^{k}$ 로 고정하여 명시적(explicit)으로 처리한다. 즉 역학적 연동항 $α \frac{\partial ε_{v}}{\partial t}$ 를 직전 반복에서 계산된 값으로 명시적(explicit)으로 평가하여 우변으로 이항한다. 이때 n에서 n+1로 전진하는 시간 단계의 유체 방정식은 다음과 같이 표현된다.

(7)

S \cdot \frac{p^{k + 1} - p^{n}}{∆ t} + \nabla \cdot v^{k + 1} = α_{u} \cdot \frac{T^{n + 1} - T^{n}}{∆ t} - α \cdot \frac{ε_{v}^{k} - ε_{v}^{n}}{∆ t}

식 (7)에서 역학적 연동항은 직전 반복에서 계산된 값이 되어 우변으로 넘어가며, 좌변에 남는 시간미분 항은 $S \frac{\partial p}{\partial t}$ 하나뿐이다. 즉 압력의 시간변화율을 비저류계수 $S$ 가 단독으로 감당하게 된다.

(8)

\frac{\partial p^{k + 1}}{\partial t} = \frac{1}{S} \cdot [- \nabla \cdot v^{k + 1} + α_{u} \frac{\partial T}{\partial t} - α \frac{ε_{v}^{k}}{\partial t}]

식 (8)로부터 S가 작을수록 $1 / S$ 가 커져 반복 중 압력 변동이 증폭되며, S가 0에 가까운 비압축· 강결합 조건에서는 압력 변화에 대한 직접적인 저류 효과가 감소하고, 현재 압력 변화에 대응하는 역학적 변형이 직전 반복값으로 지연되어 반영되기 때문에 연계 반복의 수렴성이 저하된다.

이러한 안정성은 반복의 오차 증폭인자 G로 정량화된다. 반복 오차를 $e^{k} = p^{k} - p^{*}$ ( $p^{*}$ 는 수렴해)로 정의하면, 순차 반복은 $e^{k + 1} = G \cdot e^{k}$ 형태로 표현되며, |G| < 1일 때 오차가 매 반복 감소하여 수렴한다. 고정변형 분할의 증폭인자 $G_{f s}$ 는 다음과 같다(Kim et al., 2011).

(9)

G_{f s} = \frac{α^{2}}{α^{2} + K_{d r} \cdot S}

식 (9)에서 분자 $α^{2}$ 는 역학적 연동의 세기를, 분모 $α^{2} + K_{d r} \cdot S$ 는 비저류계수(S)와 골격 강성 ( $K_{d r}$ )에 의한 흡수 효과를 나타낸다. 이를 결합 강도(coupling strength) $τ = \frac{α^{2}}{K_{d r} \cdot S}$ 로 표현하면 $G_{f s} = τ / (1 + τ)$ 에 해당하며, 결합이 강할수록( $τ$ 가 클수록) 1에 근접한다. 따라서 S > 0이면 G_fs < 1이지만, S가 작아질수록 G_fs가 1에 접근하여 반복 오차의 감소가 매우 느려진다. 특히 S → 0에서는 $G_{f s}$ → 1이 되어 수렴이 보장되지 않는다. 본 연구의 대상 문제는 S = 0에 해당하므로, 식 (8)의 증폭인자 또한 $G_{f s}$ → 1이 되어 반복 오차의 감소율이 매우 작아지고, 실용적인 반복 횟수 내에서 수렴을 확보하기 어려워진다. 이것이 다음 절의 고정응력 분할이 필요한 이유이다.

2.3.2 고정응력 분할 기법 안정화항

고정변형 분할과 고정응력 분할은 모두 Kim et al.(2011)이 제시하고 안정성을 분석한 순차 분할 기법이다. 고정변형 분할이 체적변형률의 변화율을 직전 반복값 $ε_{v}^{k}$ 으로 고정하는 것과 달리, 고정응력 분할은 총 평균응력의 변화율이 고정된다는 가정( $δ {\dot{σ}}_{m} = \partial σ_{m} / \partial t$ = 0)으로부터 체적변형률을 압력의 함수로 표현한다. 본 연구의 기여는 이 고정응력 분할을 OGS-FLAC3D THM 커플링에 적용하고, 안정화항 계수에 자유도를 부여하여 검증하는 데 있다. 총 평균응력의 구성식은 다음과 같다.

(10)

σ_{m} = K_{d r} \cdot ε_{v} - α \cdot p

여기에 $δ σ_{m} = 0$ 조건을 부과하면 체적변형률 증분이 압력 증분에 종속된다.

(11)

d ε_{v} = \frac{α}{K_{d r}} d p

식 (11)을 역학적 연동항에 대입하면 다음을 얻는다.

(12)

α \frac{\partial ε_{v}}{\partial t} = \frac{α^{2}}{K_{d r}} \frac{\partial p}{\partial t}

즉 역학적 연동항이 압력에 대한 항으로 변환된다. 이 변환을 반영하기 위해, 안정화항 계수( $β$ )를 $β = α^{2} / K_{d r}$ 로 두고 식 (2)에 $β \frac{\partial p}{\partial t}$ 를 더한 뒤 동일한 항을 다시 빼면, 원 방정식과 수학적으로 동일한 형태를 유지하면서 다음과 같이 재구성된다.

(13)

(S + β) \frac{\partial p}{\partial t} + [α \frac{\partial ε_{v}}{\partial t} - β \frac{\partial p}{\partial t}] - α_{u} \frac{\partial T}{\partial t} + \nabla \cdot v = 0

여기서 더해진 $β \frac{\partial p}{\partial t}$ 는 좌변의 저장항과 결합하여 $(S + β) \frac{\partial p}{\partial t}$ 을 이루고, 다시 빼진 $- β \frac{\partial p}{\partial t}$ 는 역학적 연동항과 묶여 대괄호 안의 보정항 $[α \frac{\partial ε_{v}}{\partial t} - β \frac{\partial p}{\partial t}]$ 을 이룬다. 좌변의 $S + β$ 는 비저류계수와 안정화항을 함께 포함하며, 대괄호 항은 반복 과정에서 실제 연동항과 그 근사치의 차이를 보정하는 역할을 한다(수렴 시 상쇄, 2.3.3 참조). 안정화항은 새로운 물리적 저류항이 아니라, 순차 수리 해석 시 아직 완전히 갱신되지 않은 역학적 되먹임 효과를 수치적인 저류항 형태로 근사하여 연계 반복을 안정화하는 계수라는 점을 유의할 필요가 있다. 고정변형 분할과 동일한 분석을 적용하면 고정응력 분할의 증폭인자는 다음과 같이 근사된다.

(14)

G_{s s} \approx \frac{α^{2}}{α^{2} + K_{d r} \cdot (S + β)}

식 (9)의 S가 식 (14)에서 $S + β$ 로 대체되므로, S가 0이더라도 $β$ > 0에 의해 분모가 유한한 값을 유지하여 $G_{s s}$ < 1이 보장된다. 따라서 고정응력 분할은 비저류계수와 무관하게 무조건 안정성 (unconditional stability)을 확보한다(Kim et al., 2011). 이는 강결합 조건에서 조건부 안정성만을 갖는 고정변형 분할과 대비된다. 한편 역학적 연동항은 실제 해석에서 FLAC3D가 반환하는 체적변형률 증분 $α \cdot (ε_{v}^{k} - ε_{v}^{n}) / ∆ t$ 로 소스/싱크 항을 구성하여 반영된다.

2.3.3 상쇄에 의한 β 독립성

본 절에서는 식 (13)에 도입한 안정화항이 수렴 시 상쇄되어 원 방정식의 해를 보존하는 과정을 보다 자세히 기술한다. 안정화항의 반복 보정은 앞서 정의한 시간 단계 인덱스 n과 반복 인덱스 k로 표현된다. 이 때 시간 미분의 기준점은 좌변과 우변 모두 직전 시간 단계 $p^{n}$ 이며, 식 (13)의 반복식은 다음과 같이 표현된다.

(15)

(S + β) \cdot \frac{p^{k + 1} - p^{n}}{∆ t} + α \cdot \frac{ε_{v}^{k} - ε_{v}^{n}}{∆ t} - α_{u} \cdot \frac{T^{n + 1} - T^{n}}{∆ t} + \nabla \cdot v^{k + 1} = β \cdot \frac{p^{k} - p^{n}}{∆ t}

좌변의 $(S + β) \cdot \frac{p^{k + 1} - p^{n}}{∆ t}$ 은 안정화항을 포함한 압력 시간항이며, 우변의 $β \cdot \frac{p^{k} - p^{n}}{∆ t}$ 은 직전 반복값을 기준으로 하는 안정화 보정항이다. 반복이 수렴하여 $p^{k + 1} = p^{k} = p^{*}$ 가 되면, 좌변의 안정화 성분과 우변의 보정항이 모두 $β \cdot \frac{p^{*} - p^{n}}{∆ t}$ 이 되어 서로 상쇄된다.

(16)

S \cdot \frac{p^{*} - p^{n}}{∆ t} + α \cdot \frac{ε_{v}^{*} - ε_{v}^{n}}{∆ t} - α_{u} \cdot \frac{T^{n + 1} - T^{n}}{∆ t} + \nabla \cdot v^{*} = 0

식 (16)은 안정화항 $β$ 가 소거된 원 방정식 (2)의 이산 형태와 동일하다. 따라서 수렴해 $p^{*}$ 는 안정화항 계수 $β$ 의 값에 무관하며, 이를 $β$ 독립성이라 한다. 안정 하한으로는 Kim et al.(2011)이 $β \geq \frac{α^{2}}{K_{d r}}$ 을 제시하였으며, Mikelić and Wheeler(2013)는 $β = \frac{α^{2}}{2 K_{d r}}$ 로도 수축(contraction)이 보장됨을 보였다. 본 연구에서는 Mikelić and Wheeler(2013)의 하한을 채택하였다.

다만 $β$ 독립성은 수렴해가 $β$ 에 무관함을 의미할 뿐, 임의로 선택 가능한 값임을 뜻하지 않는다. $β$ 는 수렴해의 위치는 바꾸지 않으나, 식 (14)의 증폭인자를 통해 수렴에 이르는 반복 경로와 속도에는 영향을 미친다. 즉, 고정응력 분할 안정화의 목적은 원 지배방정식의 수렴해를 변경하여 정확도를 직접 보정하는 것이 아니라, 순차 분할 과정에서 안정적인 수렴 경로를 제공하는 데 있다. $β$ 가 안정 하한 미만일 경우 수렴이 느려지거나 보장되지 않으며, 하한 이상에서 무조건 안정성이 확보된다. 따라서 β는 안정 하한 이상에서 수렴 성능을 고려하여 설정되어야 하며, β의 수렴해에 대한 무영향과 반복수에 대한 영향은 3장에서 실증한다.

2.4 요약

본 장에서는 순차 THM 커플링(T-H-M, FLAC3D 매 반복 호출)에 고정응력 분할을 적용한 안정화 기법을 기술하였다. 순차 분할에서 미지인 역학적 연동항( $α \frac{\partial ε_{v}}{\partial t}$ )을, 총 평균응력 고정 가정 ( $δ {\dot{σ}}_{m} = \partial σ_{m} / \partial t = 0$ )으로부터 유도되는 관계(식(11))를 통해 압력에 대한 안정화항 $β \frac{\partial p}{\partial t}$ 로 변환하여 유체 방정식 좌변에 도입하고, 실제 연동항과의 차이를 우변에서 $β \cdot \frac{p^{k} - p^{n}}{∆ t}$ 로 반복 보정하였다. 이 안정화항은 비저류계수와 무관하게 무조건 안정성을 제공하며, 이는 고정변형 분할이 강결합에서 갖는 조건부 안정성과 대비된다. 또한 수렴 시 좌·우변에서 상쇄되어 원 방정식의 해를 보존한다( $β$ 독립성). 안정 하한으로는 Mikelić and Wheeler(2013)의 $β = \frac{α^{2}}{2 K_{d r}}$ 를 채택하였으며, $β$ 가 수렴해에는 영향을 주지 않으나 수렴 성능에는 영향을 미침을 3장에서 실증한다.

3. 검증: Booker-Savvidou 점열원 문제

3.1 문제 정의

Booker-Savvidou 점열원 문제는 포화 다공질 매질 내 점열원에 의한 THM 복합거동에 대한 해석해가 존재하여, 커플링 해석의 검증 벤치마크로 널리 사용된다(Booker and Savvidou, 1985). 본 연구에서는 Booker-Savvidou 원 해석해의 응력 부호 오류를 수정한 Chaudhry et al.(2019)의 해석해를 기준으로 수치해를 검증하였다.

문제의 대칭성을 이용하여 전체 영역의 1/8만을 모델링하였다. 해석 영역은 한 변의 길이가 5 m인 정육면체이며, 점열원은 대칭 모서리인 원점에 위치하고 열원 강도는 100 W이다. 초기 온도와 간극수압은 모두 0으로 설정하였다. 격자는 약 24,700개의 사면체 요소(절점 약 5,500개)로 구성하였으며, 점열원 근방의 급격한 변화를 반영하기 위해 열원 근방에서 약 0.01 m, 외곽 경계에서 약 0.4 m의 요소 크기를 갖도록 국소 세밀화하였다(Fig. 2).

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-05/N0120360506/images/ksrm_2026_365_545_F2.jpg

Fig. 2.

Finite element mesh of the one-eighth symmetric domain for the point heat source problem

경계조건으로 세 대칭면(x = 0, y = 0, z = 0)에는 법선방향 변위를 구속하였으며, 외곽 경계면은 변위를 자유로 두고 온도와 간극수압을 0으로 설정하였다. 점열원으로부터 해석 영역이 충분히 크므로 외곽 경계의 영향은 무시할 수 있다.

해석에 사용한 물성은 Table 1과 같다. 지반의 탄성계수와 포아송비(E = 5,000 MPa, ν = 0.3)로부터 배수 체적탄성계수는 $K_{d r}$ = 4,166.7 MPa이며, Biot 상수 $α$ = 1이다. 본 문제에서는 유체와 고체 알갱이를 모두 비압축으로 가정하므로 비저류계수 S = 0이 된다. 열팽창계수는 유체의 경우 $α_{f}$ = 2.0×10^-4 /K, 고체의 경우 1.5×10^-5 /K(선형열팽창계수 $α_{s}$ = 5.0×10^-6 /K)이다. 안정화항 계수는 Mikelić and Wheeler(2013)의 하한을 채택하여 $β = \frac{α^{2}}{2 K_{d r}}$ = 1.2×10^-10로 설정하였다.

Table 1.

Material properties used in the analysis

Material	Variable	Unit	Value
Water	Density	kg/m³	1,000
	Viscosity	Pa·sec	0.001
	Specific heat capacity	J/kg·K	4,200
	Thermal conductivity	W/m·K	0.6
	Volumetric thermal expansion coefficient	1/K	2.0×10^-4
Ground	Density	kg/m³	2,290
	Young’s modulus	MPa	5,000
	Poisson ratio	-	0.3
	Thermal conductivity	W/m·K	1.83
	Specific heat capacity	J/kg·K	920
	Volumetric thermal expansion coefficient	1/K	1.5×10^-5
	Porosity	-	0.16
	Permeability	m²	2.0×10^-20
	Biot coefficient	-	1

수치해와 해석해의 비교를 위한 계측점은 Table 2와 같이 점열원으로부터 서로 다른 거리에 위치한 두 지점(P1, P2)을 선정하였다. 점열원 특이점의 영향을 함께 살펴보기 위해 점열원으로부터의 거리(R)가 다른 두 점을 사용하였다.

Table 2.

Coordinates of the monitoring points

Monitoring points	x (m)	y (m)	z (m)	R (m)
P1	0.2427	0.1747	0.1461	0.333
P2	0.3779	0.2828	0.1747	0.503

해석의 총 시간은 2×10⁶ 초이며, 가변 시간간격(adaptive time stepping, 최대 Δt = 36,000 s)을 적용하였다. 각 시간 단계 내 순차 반복(Picard)의 수렴 정도는 반복 간 최대 변화량으로 확인하였으며, 간극수압 100 Pa, 온도 0.1 K 미만을 기준으로 충분히 수렴하도록 반복을 수행하였다.

3.2 해석 결과

본 절에서는 2장의 고정응력 분할을 적용한 THM 순차 커플링 해석 결과를 Chaudhry et al.(2019)의 해석해와 비교하였다. 해석해는 Chaudhry et al.(2019)에 제시된 식을 사용하였으며, 그 상세한 유도는 해당 문헌을 참조한다. 안정화항 계수는 3.1에 제시한 바와 같이 $β$ = 1.2×10^-10로 설정하였으며, 각 시간 단계에서 온도·간극수압·변위가 충분히 수렴하도록 순차 반복을 수행하였다.

수치해의 정확도는 첨두 상대오차와 정규화 제곱근평균오차(Normalized Root Mean Square Error, NRMSE)로 평가하였다. 첨두 상대오차는 계측점 최댓값의 진폭 정확도를 나타내나, 곡선의 위상이나 형상 오차는 반영하지 못하는 한계가 있다. 이를 보완하기 위해 전 시간이력에 대한 NRMSE를 함께 사용하였으며, 이는 다음과 같이 정의된다.

(17)

N R M S E = \frac{\sqrt{\sum {(u_{n u m} - u_{a n a})}^{2}} / N)}{u_{a n a_\max} - u_{a n a_\min}} \times 100 (%)

여기서 $u_{n u m}$ 과 $u_{a n a}$ 는 각각 수치해와 해석해, $N$ 은 비교에 사용한 시간 이력 데이터의 개수이다. 분모는 해석해의 범위(range)로 정규화한 것으로, 단위가 서로 다른 물리량 간의 비교를 가능하게 한다. 또한 물리 단위의 절대 잔차를 확인하기 위해 절대 제곱근평균오차(RMSE)를 결과 비교 시 병기하였다.

간극수압의 시간이력을 Fig. 3(a)에 나타내었다. 두 계측점 모두 가열 초기에는 간극수압이 0에 가깝게 유지되다가, 열전달이 계측점에 도달하는 시점부터 유체의 열팽창에 의해 급격히 상승하여 첨두값에 이른 후 배수에 의해 점차 소산되는 거동을 보인다. 첨두 간극수압은 점열원에 가까운 P1에서 약 1.6 MPa, 상대적으로 먼 P2에서 약 1.05 MPa로 나타났다. 수치해는 이러한 상승·첨두·소산의 전 과정을 해석해와 잘 일치하도록 재현하였으며, 첨두 상대오차는 P1에서 −2.73%, P2에서 −1.60%, 전곡선 NRMSE는 각각 1.16%, 1.42%로 나타났다(Table 3). 절대 잔차 또한 RMSE 기준 약 0.02 MPa로, 첨두값(약 1.6 MPa)의 1% 남짓에 불과하여 해석해와 사실상 일치한다. 이는 강결합 조건에서도 고정응력 기반 분할을 적용한 OGS-FLAC3D 해석 결과가 간극수압을 정확히 재현함을 보여준다.

Fig. 3(b)는 온도의 시간 이력을 보여주고 있다. 온도는 점열원으로부터의 열전도에 의해 상승하며, 수치해가 해석해와 잘 일치함을 확인할 수 있다. 첨두 상대오차는 P1에서 −1.40%, P2에서 −1.62%, NRMSE는 각각 1.44%, 1.65%이다.

다음으로 변위와 유효응력의 시간 이력을 Fig. 3(c)~(f)에 나타내었다. 응력은 점열원 문제의 구면 대칭성을 고려하여 반경방향 유효응력( $σ_{r r}$ )과 접선방향 유효응력( $σ_{θ θ}$ ), 그리고 직교좌표계 전단응력( $σ_{x y}$ )에 대해 비교하였다. 구면 대칭 조건에서 주응력 방향은 반경방향 및 접선방향과 일치하므로 반경·접선 유효응력은 각각 FLAC3D의 주응력으로부터 취하였으며, 이는 Chaudhry et al.(2019)이 제시한 구면좌표계 응력 공식에 대응한다. 한편 구면좌표계에서 전단응력은 대칭성에 의해 0이 되므로(Chaudhry et al., 2019), 전단 거동의 검증을 위해서는 직교좌표계 성분 $σ_{x y}$ 를 사용하였다. 변위와 세 응력 성분 모두 수치해가 해석해의 시간이력을 잘 재현하였으며, 첨두 상대오차와 NRMSE는 Table 3과 같다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-05/N0120360506/images/ksrm_2026_365_545_F3.jpg

Fig. 3.

Comparison between the numerical and analytical solutions at P1 and P2

Table 3.

Peak relative error, NRMSE, and absolute RMSE of the numerical solution against the analytical solution at the monitoring points ( $β$ = 1.2×10^-10)

Variable	P1			P2
Variable	Peak(%)	NRMSE(%)	RMSE	Peak(%)	NRMSE(%)	RMSE
Temperature	-1.40	1.44	0.176℃	-1.62	1.65	0.121℃
Pore pressure	-2.73	1.16	0.0186 MPa	-1.60	1.42	0.0150 MPa
Displacement	-2.95	1.37	0.00048 mm	-1.178	1.04	0.00038 mm
$σ_{r r}$	-2.13	1.38	0.0123 MPa	-5.55	1.71	0.0118 MPa
$σ_{θ θ}$	-5.17	1.76	0.0179 MPa	-2.68	1.08	0.0086 MPa
$σ_{x y}$	-4.75	1.30	0.0040 MPa	-10.84	3.13	0.0081 MPa

전반적으로 수치해는 온도·간극수압·변위·응력의 모든 물리량에 대해 해석해와 잘 일치하며, 전곡선 NRMSE는 대부분 2% 이내로 나타났다. 이는 고정응력 분할 기반 안정화항을 적용한 순차 THM 커플링이 강결합 조건에서도 해석해를 정확히 재현함을 의미한다. 일부 응력 성분에서 첨두 상대오차가 5~10%로 크게 나타나나, 동일 구간의 NRMSE가 2% 내외인 점을 고려하면 이는 오차가 첨두 근방에 국소화되어 있고 전 시간이력은 잘 일치함을 의미한다. 특히 P2의 전단응력에서 첨두 상대오차가 −10.84%로 나타났으나, 이는 해당 위치에서 $σ_{x y}$ 의 절대 크기가 작아 상대오차가 확대된 것으로, 절대 오차에 해당하는 RMSE는 0.008 MPa로 작다. 이러한 잔여오차 주요 원인 중 하나로 점열원 특이점(점열원~계측점 간 거리, R → 0에서 발산) 근방의 이산화 한계가 추정되며, 이에 대해서는 4장에서 상세히 논의한다.

3.3 안정화항 적용 여부에 따른 비교

안정화항 적용 여부가 연계 수렴성과 해석 결과에 미치는 영향을 확인하기 위해, 안정화항을 제거한 해석을 수행하여 3.2의 결과와 비교하였다. 구체적으로 식 (15)의 반복식에서 안정화항에 해당하는 항, 즉 좌변의 $β \cdot \frac{p^{k + 1} - p^{n}}{∆ t}$ 와 우변의 안정화항 $β \cdot \frac{p^{k} - p^{n}}{∆ t}$ 를 모두 제거하여 저장항과 역학적 연동항만으로 유체 방정식을 구성하였다. 이는 2.3.1에서 기술한 고정변형 분할에 해당한다. 이 경우 기준 해석과 동일한 S=0 조건에서는 고정변형 분할의 연계 수렴을 확보하기 어려워(해석이 발산), 해석 수행을 위해 S = 1.2 × 10⁻¹⁰ Pa⁻¹ 를 부여하였다. 따라서 두 결과의 차이에는 안정화항 유무와 더불어 비저류계수 차이가 포함되므로 해석 결과의 정량적 차이를 안정화항의 유무만으로 직접 귀속할 수는 없다는 점에 유의할 필요가 있다. 반복 조건은 3.2와 동일하게 설정하였다.

간극수압의 시간이력을 Fig. 4에 나타내었다. 안정화항을 제거하고 작은 비저류계수를 적용한 경우, 간극수압의 첨두 상대오차가 P1에서 약 +67%, P2에서 약 +72%로, 해석해를 크게 과대평가하였다. 이는 2.3.1의 증폭인자 분석과 부합하는 결과이다. 비저류계수가 매우 작은 조건에서는 고정변형 분할의 증폭인자 $G_{f s}$ 가 1에 근접하므로 순차 반복이 동일한 반복 조건에서 수렴에 이르지 못한다. 물리적으로는 안정화항이 없으면 역학적 연동에 의한 압력 완화가 반복 내에서 충분히 반영되지 못하여 간극수압이 과도하게 상승한 상태로 남게 된다. 반면 3.2에서 살펴본 바와 같이 안정화항 적용 시에는 첨두 상대오차가 −2.73%(P1), −1.60%(P2)로 나타나, 안정화항이 강결합 조건에서 순차 반복의 수렴을 실용적인 반복 범위 내로 가속시켜 안정적인 연계 수렴을 확보하는 데 중요한 역할을 하였음을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-05/N0120360506/images/ksrm_2026_365_545_F4.jpg

Fig. 4.

Pore pressure without the stabilization term (fixed-strain split), showing significant overestimation compared to the analytical solution

3.4 안정화항 계수에 따른 영향

2.3.3에서 기술한 바와 같이, 안정화항 계수 $β$ 는 수렴해 자체를 변화시키지 않으나 수렴 과정에는 영향을 미친다. 이를 확인하기 위해 $β$ 를 안정 하한(1.2×10^-10)과 그 이상의 값(2.4, 5.0 × 10^-10)으로 변화시켜 해석을 수행하고, 간극수압의 첨두 상대오차와 NRMSE, RMSE를 비교하였다(Table 4).

Table 4.

Peak relative error, NRMSE, and RMSE of the excess pore pressure for different values of the stabilization coefficient $β$

$β (\times 10^{- 10})$	Iteration (number)	P1			P2
$β (\times 10^{- 10})$	Iteration (number)	Peak(%)	NRMSE(%)	RMSE	Peak(%)	NRMSE(%)	RMSE
1.2	Same (10)	-2.73	1.16	0.0186 MPa	-1.60	1.42	0.0150 MPa
2.4	Same (10)	-2.72	1.14	0.0183 MPa	-1.59	1.40	0.0149 MPa
5.0	Converge (40)	-4.37	1.82	0.0293 MPa	-3.24	2.25	0.0238 MPa
5.0	Converge (40)	-2.77	1.15	0.0185 MPa	-1.64	1.41	0.0149 MPa

* “Same (10)” indicates the case with the maximum iterations limited to 10, identical to other analyses; “Converge (40)” corresponds to the case iterated until satisfying the convergence criteria in Section 3.1 (pore pressure: 100 Pa, temperature: 0.1 K), which required approximately 40 iterations.

Fig. 5 및 Table 4에 나타난 바와 같이 $β$ = 1.2, 2.4 × 10^-10에서 압력 첨두오차는 각각 −2.73%, −2.72%로 동일하게 나타나, 안정화항 계수가 수렴해를 변화시키지 않는 $β$ 독립성을 확인할 수 있다. $β$ = 5.0 × 10^-10의 경우 동일한 반복 조건에서는 첨두오차가 −4.37%로 증가하나, 반복을 수렴에 이를 때까지 수행하면 −2.77%로 $β$ = 1.2의 −2.73%와 거의 일치한다. 이는 $β$ 가 수렴해 자체는 바꾸지 않되 수렴에 이르는 속도에는 영향을 미침을 의미하며, 2.3.3의 이론적 성질과 부합한다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-05/N0120360506/images/ksrm_2026_365_545_F5.jpg

Fig. 5.

Pore pressure for different stabilization coefficients, showing $β$ -independence of the converged solution

이러한 $β$ 의 영향은 양방향으로 나타난다. $β$ 가 안정 하한보다 작으면 역학적 완화가 과소 반영되어 수렴이 저하된다. 그 극단인 안정화항을 완전히 제거한 경우, 비저류계수가 0이면 직접적인 압력 저류 효과가 사라지고 현재 압력 변화에 대응하는 역학적 연동이 지연되어 반영되므로 연계 반복의 수렴성이 크게 저하된다(2.3.1). 또한 이를 보완하기 위해 비저류계수에 작은 값을 부여하더라도 압력이 크게 과대평가된다(3.3). 반대로 $β$ 가 과도하게 크면 유체 방정식 좌변의 유효 저류계수 (S + β)가 증가하여 반복당 압력 보정량이 지나치게 감쇠되고, 동일한 수렴 기준에 도달하기 위한 연계 반복 횟수가 증가할 수 있다. 따라서 $β$ 는 계산 효율 측면에서는 안정성을 확보할 수 있는 범위에서 안정 하한 근방의 값을 사용하는 것이 합리적이며, 무조건 크거나 작은 값이 유리한 것은 아니다. 본 연구에서 채택한 Mikelić and Wheeler(2013)의 하한인 $β = α^{2} / 2 K_{d r}$ 은 이러한 관점에서 실용적으로 타당한 선택이라 판단된다.

3.5 커플링 순서에 따른 비교

2.2에서 물리적 인과 순서에 따라 T-H-M을 채택하였으며, 본 절에서는 커플링 순서가 해석 결과에 미치는 영향을 확인하기 위해 H-T-M 순서로 수행한 결과와 비교하였다. 간극수압의 첨두 상대오차와 NRMSE를 Table 5에 나타내었다.

Table 5.

Comparison of the pore pressure between the T-H-M and H-T-M coupling orders

Coupling order	P1			P2
Coupling order	Peak(%)	NRMSE(%)	RMSE	Peak(%)	NRMSE(%)	RMSE
T-H-M	-2.73	1.16	0.0186 MPa	-1.60	1.42	0.0150 MPa
H-T-M	-2.72	1.15	0.0184 MPa	-1.59	1.41	0.0150 MPa

두 순서 모두 동일한 반복 조건에서 계측점의 간극수압이 사실상 동일한 값으로 나타나(Fig. 6), 본 벤치마크 및 적용 조건에서는 시험한 두 연계 순서에 따른 수렴해의 차이가 매우 작게 나타났다. 다만 본 문제에서는 온도장이 수리 및 역학 과정으로부터 받는 영향이 제한적이므로, 이를 일반적인 비선형 THM 문제로 확대 해석 시에는 확인할 필요가 있다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-05/N0120360506/images/ksrm_2026_365_545_F6.jpg

Fig. 6.

Pore pressure for the T-H-M and H-T-M coupling orders, showing a negligible effect of the tested coupling order under the present benchmark conditions

4. 고 찰

3장의 검증 결과에서 수치해는 온도·간극수압·변위·응력의 전 시간이력에 걸쳐 해석해와 잘 일치하였다. 다만 일부 응력 성분의 첨두 근방에서 상대오차가 다소 크게 나타났으며, 이러한 잔여오차의 주요 원인 중 하나로 점열원 특이점 근방의 공간 이산화 한계가 추정된다. 동일한 Booker-Savvidou 문제를 대상으로 한 OpenGeoSys v5 기반의 일체형(monolithic) 연계해석(Park and Park, 2026)에서도 점열원 근방에서 수치해와 해석해 사이의 잔여오차가 보고된 바 있다. 일체형 방식은 순차 분할에 따른 연계 오차가 원리적으로 존재하지 않음에도 점열원에 가까운 영역에서 오차가 증가하는 경향을 보였으며, 이는 본 연구에서 나타난 근거리 잔여오차가 순차 연계방식 자체뿐만 아니라 점열원 특이점과 공간 이산화의 영향과도 관련됨을 시사한다.

안정화항의 필요성은 선행 연구와의 관계에서도 확인된다. 선행 OGS-FLAC3D 연구(Park and Park, 2022)는 변형률 기반 소스텀을 도입하였으나 안정화항을 적용하지 않았으며, 탄성계수가 작은 강결합 조건에서 간극수압이 해석해 대비 과대평가되는 경향이 있음을 보고한 바 있다. 본 연구에서 안정화항을 제거한 해석 결과(3.3)가 간극수압을 약 67% 과대평가한 결과는 이러한 경향과 정합한다. 다만 강결합 조건에서 안정화항 미적용 해석 시 나타난 큰 과대평가는 안정화항 부재뿐 아니라 서로 다른 비저류계수와 미수렴 상태의 영향을 함께 포함한다는 점을 유의할 필요가 있다.

고정응력 분할 기반 안정화항을 적용한 해석에서는 간극수압 첨두 상대오차가 P1에서 −2.73%, P2에서 −1.60%로 나타났으며, 연동해석의 안정적 수렴과 함께 해석해와 잘 일치하는 결과를 얻었다. 즉, 고정응력 분할 안정화는 강결합 조건에서 연계 반복이 안정적으로 수렴하도록 하였으며, 그 결과 얻어진 수렴해가 해석해와 잘 일치함을 확인하였다.

한편 동일한 Booker-Savvidou 문제를 대상으로 한 OpenGeoSys v5 기반의 일체형(monolithic) 연계해석(Park and Park, 2026)에서는 커플링 분할에 따른 오차가 원리적으로 존재하지 않으며, 간극수압을 해석해 대비 약 1% 수준으로 재현된 것으로 보고되었다. 본 연구의 순차 연계해석에서도 간극수압 첨두 상대오차가 −2.73%와 −1.60%로 나타나, 두 연계방식 모두 해석해를 잘 재현하는 결과를 보였다. 다만 두 연구는 해석 영역, 격자, 열원 조건 및 계측점 위치가 서로 다르므로 오차값의 직접적인 정량 비교보다는 해석해 대비 정확도 수준의 관점에서 비교하는 것이 타당하다.

종합하면, 본 연구의 기여는 Kim et al.(2011)이 제안한 고정응력 분할 기반 안정화항을 OGS-FLAC3D THM 커플링에 적용하고 Booker-Savvidou 해석해로 검증함으로써, 순차 방식으로도 해석해와 잘 일치하는 커플링 충실도를 확보하면서 모듈성이라는 순차 방식의 장점을 유지한 데 있다. 즉, 강결합 조건에서 순차 분할이 갖는 수렴성 한계를 완화하고 모듈형 순차 커플링을 심층처분 THM 해석에 신뢰성 있게 적용할 수 있는 근거를 제시하였다. 본 연구는 단상 유동과 선형 열탄성 조건의 벤치마크를 대상으로 안정화항의 효과를 검증하였으며, 점열원 특이점에 의한 근거리 잔여오차의 저감은 향후 과제로 남는다.

한편 다상 유동 및 비선형 구성모델을 포함한 OGS-FLAC3D 해석은 선행 연구(Kim et al., 2021b, 2022, 2024)에서 이미 수행된 바 있으며, 본 연구에서 검증한 안정화항을 이들 해석에 통합하여 강결합 조건에서의 정확도와 안정성을 확보하는 것이 향후 연구 방향이다.

5. 결 론

본 연구에서는 고정응력 분할(fixed-stress split) 기반의 안정화항을 OGS-FLAC3D 열-수리-역학(THM) 커플링에 도입하고, 점열원 문제의 해석해(Booker-Savvidou 문제)를 이용하여 그 정확성과 안정성을 검증하였다. 주요 결과는 다음과 같다.

1. 총 평균응력 고정 가정으로부터 유도되는 안정화항 $β = α^{2} / 2 K_{d r}$ 을 유체 방정식에 도입하여, 비저류계수가 0인 강결합 조건에서도 순차 반복의 안정성과 수렴성을 향상시켰다. 이 안정화항은 수렴 시 좌·우변에서 상쇄되어 원 방정식의 해를 보존하므로, 안정화항 계수 $β$ 는 수렴해에 영향을 주지 않는다.

2. 고정응력 분할을 적용한 수치해는 온도·간극수압·변위·응력의 전 시간이력에 걸쳐 해석해와 잘 일치하였다. 간극수압의 첨두 상대오차는 계측점에서 −2.73%, −1.60%, 전곡선 NRMSE는 2% 내외로, 강결합 조건에서도 해석해에 근접하는 정확도를 확인하였다.

3. 안정화항을 제거한 고정변형 분할 해석에서는 작은 비저류계수를 사용하였음에도 주어진 반복 횟수 내에서 수렴 기준에 도달하지 못하고 간극수압이 해석해 대비 크게 과대평가되었다. 이를 통해 본 벤치마크에서 고정응력 분할 안정화가 실용적인 연계 수렴을 확보하는 데 핵심적인 역할을 함을 확인하였다. 또한 본 벤치마크 해석 조건에서는 시험한 안정화항 계수와 T-H-M 및 H-T-M 연계 순서에 따른 최종 수렴해의 차이는 매우 작게 나타났다.

이상의 결과로부터, 고정응력 분할 기반의 안정화항을 도입한 OGS-FLAC3D는 강결합 조건의 THM 복합거동을 신뢰성 있게 해석할 수 있음을 확인하였다. 본 방법은 유한요소법(finite element method) 기반의 열-수리 해석 코드(OpenGeoSys)와 유한차분법(finite difference method) 기반의 역학 해석 코드(FLAC3D)의 강점을 결합한 순차 연계해석의 모듈성을 유지하면서, 강결합 조건에서의 수렴성과 안정성을 향상시켰다. 이는 심층처분과 같이 열-수리-역학 결합이 강한 문제에 순차 연계해석을 보다 신뢰성 있게 적용하는 데 기여할 것으로 판단되며, 향후 다상 유동, 비선형 구성모델 및 균열 지배 THM 문제로의 확대 적용에는 추가적인 검증이 필요하다.

Acknowledgements

This research was supported by the Institute for Korea Spent Nuclear Fuel (iKSNF) and National Research Foundation of Korea(NRF) grant funded by the Korea government(Ministry of Science and ICT, MSIT)(RS-2021-NR056198).

References

Booker, J.R., and Savvidou, C., 1985, Consolidation around a point heat source, International Journal for Numerical and Analytical Methods in Geomechanics, 9(2), 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

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

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

Kim, J., Tchelepi, H.A., and Juanes, R., 2011, Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits, Computer Methods in Applied Mechanics and Engineering, 200, 1591-1606.

10.1016/j.cma.2010.12.022

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 CO₂ injection analysis in the Otway Basin, Australia, Environmental Earth Sciences, 80(121), 1-15.

10.1007/s12665-021-09411-1

Kim, T., Park, C.-H., Lee, C., and Kim, J.-S., 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 and 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

Kim, T., Park, C.-H., Lee, C., and Kim, J.-S., 2022, A Thermo-Hydro-Mechanical Coupled Numerical Simulation on the FE Experiment: Step 1 Simulation in Task C of DECOVALEX-2023, Tunnel and Underground Space, 32(6), 518-529.

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

Mikelić, A., and Wheeler, M.F., 2013, Convergence of iterative coupling for coupled flow and geomechanics, Computational Geosciences, 17, 455-461.

10.1007/s10596-012-9318-y

Grunwald, N., Lehmann, C., Maßmann, J., Naumov, D., Kolditz, O., and Nagel, T., 2022, Non-isothermal two-phase flow in deformable porous media: systematic open-source implementation and verification procedure, Geomechanics and Geophysics for Geo-engineering and Geo-resources, 8(107), 1-29.

10.1007/s40948-022-00394-2

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 and 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 and its case study, Tunnel and Underground Space, 30(6), 509-518.

10.7474/TUS.2020.30.6.509

Park, C.-H., and Park, J.-W., 2026, A monolithic finite element formulation for coupled THM analysis with two-phase flow, Computers and Geotechnics, 198, 108282.

10.1016/j.compgeo.2026.108282

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

10.7474/TUS.2022.32.2.144

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

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

Preview

Verification of Thermo-Hydro-Mechanical Coupling in the Sequential OGS-FLAC3D Simulator with a Fixed-Stress Split Stabilization

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(6)

Fig. 1.

Conceptual diagram of the sequential T-H-M coupling procedure in OGS-FLAC3D

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

Fig. 2.

Finite element mesh of the one-eighth symmetric domain for the point heat source problem

Table 1.

Material properties used in the analysis

Table 2.

Coordinates of the monitoring points

(17)

Fig. 3.

Comparison between the numerical and analytical solutions at P1 and P2

Table 3.

Peak relative error, NRMSE, and absolute RMSE of the numerical solution against the analytical solution at the monitoring points (β = 1.2×10-10)

Fig. 4.

Pore pressure without the stabilization term (fixed-strain split), showing significant overestimation compared to the analytical solution

Table 4.

Peak relative error, NRMSE, and RMSE of the excess pore pressure for different values of the stabilization coefficient β

Fig. 5.

Pore pressure for different stabilization coefficients, showing β-independence of the converged solution

Table 5.

Comparison of the pore pressure between the T-H-M and H-T-M coupling orders

Fig. 6.

Pore pressure for the T-H-M and H-T-M coupling orders, showing a negligible effect of the tested coupling order under the present benchmark conditions

Acknowledgements

References

Peak relative error, NRMSE, and absolute RMSE of the numerical solution against the analytical solution at the monitoring points ( $β$ = 1.2×10^-10)

Peak relative error, NRMSE, and RMSE of the excess pore pressure for different values of the stabilization coefficient $β$

Pore pressure for different stabilization coefficients, showing $β$ -independence of the converged solution