1. 서 론
유체역학과 암석역학은 오랜 기간 동안 각자의 독자적인 영역에서 전문성을 갖추며 발전해 왔다. 하지만 21세기에 당면한 범국가적인 문제들을 해결하기 위해서는 서로 분리되어 적용 가능했던 한계를 넘어서, 서로 융합하여 유기적으로 해결해야만 하는 상황에 직면하게 되었다. 이렇게 해결해야 할 대형 국책 연구 사업들에는 온실가스의 지중저장, 신재생에너지의 지중연계 효율적 활용, 사용후 핵연료의 심층처분, 각종 석유가스 및 하이드레이트의 개발과 같은 다양한 사업들이 존재한다. 이와 같은 지중 환경을 대상으로 하는 사업들에서는 지중에 지하수, 이산화탄소, 석유가스와 같은 단상 혹은 다상의 유체가 복합적으로 존재할 수 있으며, 유체의 흐름 및 압력과 대수층의 역학적 변형이 상호작용하여 지중의 환경에 영향을 미칠 수 있다.
먼저 온실가스 중 하나인 이산화탄소를 오랜 시간동안 평형을 이루어 온 지중의 대수층에 압력을 가하여 주입할 경우에 지중의 구조적 안정성 문제를 해결해야만 안전한 주입이 가능하다. 또한, 신재생 에너지의 하나인 지열발전에서 유체의 흐름을 원활히 하기 위해서 지중에 수압파쇄를 하게 될 때도 유체와 암석역학의 상호작용의 이해가 필수적으로 요구된다. 사용후 핵연료의 심층처분에서는 고준위방사성폐기물이 붕괴열을 방출하고, 이러한 열에 대한 다상 유체(기체와 액체)의 거동이 핵연료를 보관하는 용기(canister) 및 완충재(buffer material)가 묻혀 있는 지중의 응력 상태와 상호작용하여 구조물의 안정성에 어떠한 영향을 미치는지에 대해 잘 이해하고 분석하는 것이 필수적이라고 할 수 있다.
이러한 문제들을 해결하기 위하여 실내 및 현장 실험, 수치해석 등의 다양한 접근법이 활용되고 있으며, 수치해석의 경우 다양한 현장 조건을 전산해석을 통해 효율적으로 모사할 수 있다는 장점을 지니고 있다. 그러나 앞서 언급한 문제들의 경우 열(T), 수리(H), 역학(M), 화학(C)을 복합적으로 고려해야 하므로 기존과는 다른 접근이 필요하며 전 세계적으로 다양한 기법 및 프로그램들을 활용한 연구가 진행 중에 있다(Gutierrez and Makurat, 1997, Rutqvist et al., 2002, Taron et al., 2009). T-H-M 복합거동 해석의 경우 지중에서의 연계거동이 복잡한 양상을 나타내며 이를 모사하고 예측하기 위한 수식이나 해석기법은 정형화 되어 있지 않다. 국내의 경우 기존에 제시된 T-H-M 연동해석이 가능한 일부 프로그램들을 활용하여 실험실 시험 및 현장시험을 모사하는 연구들이 진행되고 있으나(Lee et al., 2017, Park et al., 2018a, b, Lee et al., 2019a, b) T-H-M 복합거동 해석을 위한 독자적인 시뮬레이터를 확보하려는 노력은 부족한 실정이다.
따라서 본 연구에서는 T-H-M 복합거동 해석을 위한 시뮬레이터 개발의 첫 단계로써 H-M 연동해석이 가능한 시뮬레이터 개발을 수행하였으며, 이를 위해 열-수리 해석 및 다상 유체 거동 해석에 강점을 지니고 있는 OpenGeoSys(Rink, 2015)와 지반 분야의 안정성 검토 및 역학 해석에 널리 쓰이고 있는 FLAC3D(Itasca, 2013)를 연계하고자 하였다. OpenGeoSys는 멀티-플랫폼을 지원하는 해석용 소프트웨어 중 하나로써 다공질 혹은 균열 매질 내의 독립적이거나 복잡한 T-H-M 복합 거동을 해석하기 위하여 개발되었다. OpenGeoSys는 범용성이 높은 C언어를 기반으로 하였으며, 코드가 공개되어 있고 무료로 사용이 가능하다. 따라서 물리, 화학, 수학, 유체역학, 암석역학 등 다양한 전문가들이 각자의 코드를 자신이 관심 있는 영역에서 개발을 지속적으로 수행하고 있다. 특히 최근에는 이산화탄소 지중저장, 암반 내 유체이동, 방사성폐기물에서 발생할 수 있는 핵종 이동과 같은 다상(multiphase) 유체의 유동과 관련된 문제들의 해석에 활용되고 있다(Park et al., 2008, Park et al., 2011). 또한, 연속체 기반 대표 해석 기법인 유한요소해석(Finite Element Method, FEM)에 기반하여 복잡한 3차원 영역의 표현에 강점을 지니고 있다. 반면, 심부 암반의 안정성 평가에 중요한 역학 해석 부분은 기존의 다른 수치해석 프로그램 혹은 시뮬레이터에 비해 역학적인 분야에 대한 검증이 부족하고 제공되는 구성 모델이 한정되어 있다는 약점이 있다.
FLAC3D는 ITASCA에서 개발되었으며 암반 및 터널 등의 안정성 해석에 널리 활용되고 있는 상용 소프트웨어로 유한차분법(Finite Difference Method, FDM)에 기반하여 공학적 계산을 수행한다. 또한, 지반 안정성 해석을 위해 변형 특성에 기초한 다양한 구성 모델이 제공되어 있으며, 필요할 경우 사용자가 직접 내장 언어인 FISH 혹은 사용자 구성 모델인 UDM(User Defined Model)을 통해 수정하거나 구성할 수 있다. 그러나 FLAC3D의 경우 해석 시 Lagrangian 기법을 활용하며, 따라서 지반과 같은 고체를 대상으로 한 해석에는 적합한 특성을 보이나 유체 해석에 있어 단상 유동에 제한되어 있는 단점이 있다.
본 연구에서는 OpenGeoSys와 FLAC3D를 연동하여 각 프로그램의 단점을 보완하고 그리드 생성의 효율성, 정밀한 다상유체 해석 등의 장점을 활용하여 복잡한 H-M 복합 거동 해석을 수행할 수 있는 시뮬레이터를 개발하고자 하였으며, 검증의 첫 단계로써 포화 상태에서의 단상 유체 거동과 관련한 벤치마크 문제를 통해 개발된 시뮬레이터의 검증을 수행하고자 한다.
2. 접근 방법
수리-역학적 복합 거동해석을 위한 해석 방식으로는 하나의 코드 안에서 하나의 지배 방정식을 통해 결과를 도출하는 단일(monolithic) 해석방법과 두 가지 이상의 코드를 연계하여 순차적으로 결과를 주고 받는 순차적(sequential) 해석방법이 있다. 단일 해석방법의 경우 높은 해석 안정성과 정확성을 보장하나 고사양의 컴퓨팅 소스를 필요로 하며 코드 구성이 복잡하다는 단점이 있다. 순차적 해석방법의 경우 단일 해석방법에 비해 해석 여부가 경계조건과 입력변수의 변화에 민감하게 반응한다는 단점이 있으나 기존 코드를 수정함으로써 효율적이고 편리한 코드 구성이 가능하다는 장점을 지니고 있다.
본 연구에서는 코드 구성의 효율성과 각 코드가 지닌 장점을 활용하기 위해 순차적 해석방법을 기반으로 하여 연성해석 알고리즘을 설계하였다. Fig. 1(a)는 OpenGeoSys와 FLAC3D의 연동에 대한 개념도를 나타내고 있다. OpenGeoSys에서 먼저 수리 해석을 수행하여 변화된 압력 정보를 FLAC3D 에 전달하며, FLAC3D 에서는 전달된 정보를 바탕으로 역학 해석을 수행하게 된다. 이 때 압력 변화는 Fig. 1(a)와 같이 유효응력 계산을 통해 반영되며, 역학 해석 시 변화된 응력과 변형을 바탕으로 공극률, 투수계수, 모세관압의 변화를 계산하여 다시 OpenGeoSys 에 전달하며, 이러한 계산이 반복적으로 수행된다.
이 때, 수리-역학적 연동 과정에서 수리 특성의 변화는 공극 혹은 저류층의 응력 변화 또는 변형으로 인하여 발생하며 유체의 흐름은 질량 및 운동량 균형 방정식에 기반하여 다음과 같이 표현할 수 있다.
| $$-\nabla\;\bullet\;q=\lbrack nC_f+(\alpha-n)C_s\rbrack\;\frac{dp}{dt}+\alpha\frac{d\varepsilon}{dt}$$ | (1) |
위 식 (1)에서 q는 저류층을 흐르는 유속을 나타내는 겉보기 속도인 비유량(specific discharge, m/s)에 해당하며, n은 공극률, Cs와 Cf는 공극 및 유체의 압축률(Pa-1), p는 압력(Pa), t는 시간, α는 Biot 상수, ε는 저류층의 변형률을 의미한다. 압축률의 단위는 Pa-1 이며, 공극률, Biot 상수는 무차원으로 표현된다. 식 (1)을 통해 확인할 수 있듯이 유체의 흐름은 시간당 압력 변화()에 따라 발생하는 공극 내 유체 자체의 변화(nCf), 저류층 내 공극의 변화((α-nCs))및 저류층 자체의 변형()에 의해 결정된다. 이 때 [nCf+(α-nCs)]는 비저류율(specific storage, Ss (m-1)), 즉 단위수두의 변화에 따라 공극의 단면적을 통해 유입되거나 유출되는 물의 체적비와 관련된 항이며, 유체의 밀도 ρf (kg/m3), 중력가속도 g(m/sec2)와 함께 다음 식 (2)를 통해 비저류율을 계산할 수 있다. OpenGeoSys에서는 수리 해석 시 비저류율, 투수계수 및 모세관압의 변화를 통해 유체의 흐름이 결정되지만, 본 논문에서는 포화상태의 단상유체 거동에 해당하는 벤치마크 문제의 가정 상 투수계수 및 모세관압의 변화 없이 비저류율의 변화를 통해 유체의 거동이 발생하였으며, 추후 연구를 통해 모세관압, 다상유체에 의한 수리-역학적 복합거동 해석에 대한 검증을 수행하고자 한다.
| $$S_s=\rho_fg\lbrack C_s(1-n)+C_f\;n\rbrack$$ | (2) |
OGSFLAC의 연동과정에서 유체의 밀도와 압축률, 공극의 압축률은 매질의 특성으로써 고정된 입력변수로 활용되며, 따라서 비저류율 계산 시 역학적 특성 변화에 의해 수리 특성에 영향을 미칠 수 있는 변수는 공극률이다. 응력 혹은 변형률에 따른 공극률의 변화와 관련해서 다양한 관계식들이 제안된 바 있지만(David et al., 1994, Davies and Davies, 1999, Dong et al., 2010), 이는 일반적으로 암종, 지역적 특성에 따른 경험식으로 제시된다는 제약이 있기 때문에, 본 연구에서는 Zienkiewicz et al.(1999)이 제안한 다음 식 (3)을 적용하여 변형률 변화에 따른 공극률의 변화를 계산하였다.
| $$n=f(\varepsilon)=n_o+\alpha(\varepsilon_v-\varepsilon_{vo})+\frac1\varphi(p-p_o)$$ | (3) |
식 (3)에서 no는 초기 공극률, εv 및 εvo는 각각 체적변형률 및 초기 체적변형률, p와 po는 압력 및 초기 압력을 의미하며, α와 φ는 상수이다(Biot, 1941). Elyasi et al.(2016)은 저류층 해석시 α와 φ는 일반적으로 각각 1 및 무한대의 값을 갖는 것으로 가정한다고 보고한 바 있으며, 본 연구에서는 이를 적용하였다.
연동 해석 시 해석 순서는 Fig. 1(b)와 같다. 먼저 OpenGeoSys에서 해석이 수행되며, 시간 tk에서 tk+1로 해석이 진행되는 동안 공극률 및 투수계수는 이전 단계(tk-1)에서 전달받은 응력 및 변형률을 바탕으로 계산된 값으로 유지되며(ɸk, kk), 이를 바탕으로 유체의 압력(∆P) 및 모세관압(∆Pc)의 변화를 계산하게 된다(과정 ①). OpenGeoSys 해석이 완료되면 유체 압력 및 모세관압을 통해 계산된 압력 변화량(∆P)을 FLAC3D 로 전달하게 되며(과정 ②), 변화된 압력 정보를 바탕으로 OpenGeoSys 와 동일한 시간에 대해 FLAC3D 에서도 역학 해석을 수행하게 된다(과정 ③). FLAC3D에서 시간 tk+1까지의 역학 해석이 완료되면 변화된 응력(σk+1) 및 변형률(εk+1)을 바탕으로 공극률(∆ɸ), 투수계수(∆k)의 변화를 환산하여 OpenGeoSys 로 전달하게 되며(과정 ④), 이러한 과정이 반복적으로 수행된다.
이러한 경우 열-수리-역학적 복합거동 해석을 위해 활발하게 쓰이고 있는 TOUGH-FLAC과 커플링 구조 및 방식에는 공통점이 있으나, 그리드 생성의 효율성에 있어 차이를 보이게 된다. OpenGeoSys는 유한요소법(Finite Element Method, FEM)을 기반으로 하였기 때문에 복잡한 3차원 도메인의 표현에 유리한 반면, TOUGH는 유한체적법(Finite Volume Method, FVM) 기반으로 개발되었으며, 요소 간 직교연결성을 유지해야 하므로(Croucher and O’Sullivan, 2013) 상대적으로 요소 생성의 유연성에 제한이 있다. 또한, 다상유체 모델의 내부 구조면에서 TOUGH는 습윤성 유체의 압력(wetting fluid pressure, Pw)과 비습윤성 유체의 포화도(non-wetting fluid saturation, Snw) 간의 관계를 풀이하는 Pw-Snw 모듈 중심으로 개발이 되어 있는 반면에, OpenGeoSys는 모세관압(capillary pressure, Pc)과 비습윤성 유체의 압력(non-wetting fluid pressure, Pnw) 간의 관계를 풀이하는 Pc-Pnw 모듈과 Pw-Snw 모듈을 각각 가지고 있다는 장점이 있다(Park et al., 2011). Jha and Juanes(2014)는 단층 내 CO2 주입과 같은 다상 유체 거동 해석 시 모세관 현상(capillary effect)이 강하게 발생하는 경우에는 모세관압을 1차 변수(primary variable)로 설정하여 공극 내 압력을 계산해야만 정밀한 해석이 가능하다고 보고한 바 있으며, 따라서 Pc-Pnw 모듈의 적용이 가능한 OpenGeoSys의 경우 다상 유체 거동 해석에 있어 장점을 지니고 있다고 판단된다.
연동의 매개체는 두 개의 방법을 디자인 단계에서 고려할 수 있는데, 첫 번째는 파일을 구성해서 연동하는 방법이고, 두 번째는 네트워크의 통신을 통해서 구현하는 방법이 있다. Table 1은 각각의 방법에 대한 장단점을 나타내고 있다. 본 연구에서는 코드 수정에 있어 보다 효율적인 파일 구성 기반 연동 방식을 통해 연계모듈을 구성하였다.
Table 1. Comparison between file-based coupling method and TCP/IP-based coupling method
3. OGSFLAC 검증
개발된 모듈 혹은 프로그램의 검증을 위해서는 수학적으로 풀이된 해석적 해(analytic solution)이 존재하는 벤치마크 모델을 통해 검증하는 과정이 필수적이다. 본 연구에서는 수리-역학적 연계해석 검증 모델로 널리 활용되고 있는 Terzaghi 압밀 문제의 해석적 해를 활용하여 개발된 모듈의 검증 과정을 수행하였다.
3.1 벤치마크 모델
Terzaghi의 압밀 문제는 외부 하중에 의해 발생한 과잉간극수압이 시간이 지남에 따라 소산되는 과정을 시간의 함수로 표현한 문제로써 다음의 가정들을 전제로 한다(Terzaghi et al., 1996).
- 외부하중의 재하는 1차원 방향으로 발생
- 유체의 흐름은 Darcy의 법칙에 따름
- 투수계수는 일정
- 2차 압축은 발생하지 않음
- 변형은 무시할 수 있을 만큼 작음
- 재하 대상 매질은 포화되어 있으며 균질함
Terzaghi의 1차원 압밀방정식은 편미분 형태로 주어지며, 식 (4) 및 (5)는 위치 x에서 시간에 따른 압력 및 변위를 나타내는 식으로서 s는 라플라스 변환을 통해 표현된 시간에 대한 항이며, pi는 초기 재하 하중, L은 모델의 초기 길이를 의미한다.
| $$\overline p(x,s)=-p_i\;\frac{\cosh\;(\sqrt{s/\chi}\;x)}{s\;\cosh\;(\sqrt{s/\chi}\;L)}$$ | (4) |
| $${\overline u}_s(x,s)=-\frac{p_i}E\left(1-\frac{2\nu^3}{1-\nu}\right)\;\frac{\sinh\;(\sqrt{s/\chi}x)}{s\sqrt{s/\chi}\;\cosh\;(\sqrt{s/\chi}L)}$$ | (5) |
| $$\chi=\frac{kE}\mu/\left(1-\frac{\nu^2}{1-\nu}\right)$$ | (6) |
식 (4) 및 (5)의 χ는 라플라스 변환 풀이 시 활용되는 변수로써 투수계수 k(m2), 탄성계수 E(Pa), 유체의 점성계수 µ(Pa·sec), 포아송비 ν 간의 관계에 의해 정의되며 식 (6)을 통해 계산된다.
본 연구에서는 총 길이 20 m의 모델을 가정하였으며, Table 2는 해석적 해 및 수치해석을 수행하는 과정에서 적용된 입력변수를 나타낸다. 초기 재하 하중은 0.1 MPa로 설정하였으며 압밀이 이루어지는 시간은 압밀이 충분히 완료되도록 5350초를 설정하였다. 또한, 편미분 방정식의 라플라스 변환 풀이를 위해 오픈 소스 코드인 파이썬을 활용하여 해석적 해를 획득하였다.
Fig. 2는 생성된 수치 모델 과 경계조건을 보여주고 있다. 탄성모델 및 포화 상태를 가정하였으며, 유체는 상온 조건의 물을 가정하여 해당하는 물성을 활용하였다. 초기 조건 모사 시 비배수 상태에서 하중 재하에 따른 과잉공극수압의 발생 과정은 모사가 어려웠기 때문에 모델 전체에 걸쳐 0.1 MPa의 공극 수압이 발생한 상태를 초기 조건으로 가정하였다. Terzaghi의 압밀 문제는 1차원 방향의 하중 재하를 가정하고 있으므로 해석 시 하중 재하 방향(x 방향) 이외의 방향(y, z)에 대해서는 각 노드마다 고정 경계조건을 적용하여 변위 발생을 막고자 하였다. 경계조건의 적용 및 초기 조건 모사가 완료된 후에는 우측면을 배수 상태(drainage condition)로 하여 압밀이 일어나는 과정을 모사하였다.
3.2 OpenGeoSys 및 FLAC3D 해석 결과
OGSFLAC 시뮬레이터에 대한 검증 해석에 앞서 OGSFLAC 시뮬레이터가 OpenGeoSys 및 FLAC3D의 개별 프로그램에 비해 어떠한 장단점을 보이는지 파악할 필요가 있다고 판단하였다. 이를 위해 Table 2에 제시된 바와 같이 두 가지 변수 세트를 설정하여 OpenGeoSys 및 FLAC3D의 개별 프로그램을 통해 벤치마크 모델에 대한 수리-역학적 연동해석을 수행하였다. 또한, 해석 방법에 따른 차이점을 파악하기 위해 OpenGeoSys의 경우 수리-역학적 연동해석 시 하나의 행렬(matrix)을 구성하여 연동해석을 동시에 풀어내는 단일 해석방법과 본 연구의 시뮬레이터에서 적용하고자 하는 방식과 같이 수리해석과 역학해석을 순차적으로 풀어내는 순차적 해석방법을 적용하여 결과를 비교하였다.
해석 수행 결과 OpenGeoSys는 Case A, B의 경우 모두 원활히 해석이 수행되는데 반해서 FLAC3D는 Case B에서만 해석 결과를 획득할 수 있었다. 유한차분법에 기초한 FLAC3D의 경우 유한요소법에 비해 수리 해석 시 요소의 크기, 입력 물성에 따라 정확도 및 해석 속도가 민감하게 반응하기 때문인 것으로 판단되며, Case B의 경우 압밀이 발생할 수 있는 일반적인 연약 지반의 물성과는 차이가 있으나 본 연구에서는 해석적 해의 검증을 목적으로 하므로 물성을 가정하여 활용하였다.
해석을 통해서 x = 8.0 m 인 지점에서의 시간에 따른 공극수압, x방향 변위의 변화와 함께 해석이 완료된 5000 sec일 때 모델의 위치별 변위 결과를 비교 분석하였다. Fig. 3는 각 프로그램 별로 x = 8.0 m 위치에서 시간에 따른 공극수압 및 변위에 대한 해석 결과를 나타내고 있다. 시간 t = 0 sec일 때, 0.1 MPa의 초기 하중이 적용됨에 따라 비배수 상태에서의 초기 과잉공극수압 0.1 MPa이 모든 프로그램에서 잘 모사되었음을 확인할 수 있다. 초기 조건에서는 하중의 대부분이 포화상태인 간극수에 작용하고 있으므로 변위가 발생하지 않으나, 시간경과에 따라 0 MPa의 공극수압조건이 적용된 우측 경계면에서 배수가 진행되며, 이로 인해 모델 전체에 걸쳐 공극수압이 감소하게 되어 탄성체로 가정된 암석에 작용하는 유효수직응력이 증가하게 된다. 각 프로그램별 시간에 따른 공극수압의 변화 결과는 해석적 해와 비교적 잘 일치하는 것을 확인할 수 있다(Fig. 3(a)).
시간이 경과함에 따라 Fig. 3(b)에서와 같이 모델 전체에 걸쳐 유효수직응력의 증가로 x 방향 변위가 증가하는(수축하는) 압밀 현상이 발생하게 되며, 각 프로그램별 해석 결과가 해석적 해와 잘 일치하는 것을 확인할 수 있다. Fig. 3(c)는 5000 sec에서 각 위치별 변위 결과를 나타내고 있다. 변위가 고정된 좌측 경계면(x = 0.0 m)에서는 압밀로 인한 수축현상이 발생하지 않으며, 유효수직응력의 증가가 가장 큰 우측 경계면(x = 20.0 m)에서 가장 큰 변위가 발생한 것을 확인할 수 있다.
해석 속도에 있어서는 OpenGeoSys 단일 해석방법, OpenGeoSys 순차적 해석방법, FLAC3D 순차적 해석방법 순으로 나타났으며, 이는 기존에 보고된 해석방법에 따른 해석 속도 결과와도 동일한 결과를 나타냈다(Kim et al., 2011).
Table 2. Input parameters for benchmark model in Terzaghi’s consolidation problem: Cases A and B were compared to investigate the proper input parameters
3.3 OGSFLAC 해석 결과
개발된 OGSFLAC 시뮬레이터의 포화 상태에서의 단상 유체 거동에 대한 검증을 수행하기 위해 Table 2의 Case B 입력 변수를 활용하여 해석을 수행하였다. 또한, 기 검증되었으며, OGSFLAC 시뮬레이터와 같은 연계방식을 적용하고 있는 TOUGH-FLAC을 활용하여 결과를 비교하였다. OGSFLAC 및 TOUGH-FLAC 해석 시 투수계수는 고정된 것으로 가정하였으며, 연동해석 시 식(2)에 제안된 바와 같이 매 연동해석 단계마다 공극률의 변화에 따라 비저류율의 변화를 반영하여 유체 유동이 발생할 수 있도록 양방향 연동 해석을 수행하였다.
Fig. 4은 시간에 따른 모델 전체에서의 압력 변화를 보여주고 있다. 초기 비배수 상태(0 sec)에서는 과잉공극수압 0.1 MPa이 모델 전체에 걸쳐 형성되어 있으며, 압밀 해석을 위해 우측면에 배수 조건이 설정된 이후부터는 시간 경과시 우측면을 통해 배수가 이루어지기 시작하고 그에 따라 모델 내 공극수압이 점차적으로 낮아짐을 확인할 수 있다. 압밀 시험이 거의 종료된 5000sec 에 이르러서는 모델 내부의 배수가 거의 완료되어 모델 전체에 걸쳐 공극 수압이 0에 가까워진 것을 확인할 수 있다.
Fig. 5는 압밀 과정에서 계측된 압력 및 변위의 변화 양상을 보여주고 있으며, Fig. 5(a) 와 (b)는 각각 계측지점인 x=8.0 m에서 측정된 시간에 따른 압력 및 변화에 대한 그래프이다. Fig 4에서 나타난 바와 같이 시간 경과에 따라 압밀로 인한 배수 발생으로 압력이 감소하는 것을 확인할 수 있으며(Fig. 5(a)), 식 (7)에서 보이는 바와 같이 압력(P) 감소 시 유효 응력(σ′)의 증가로 인해 압축 방향으로의 변위가 증가하는 것으로 나타났다.
| $$\sigma'=\sigma-P$$ | (7) |
Fig. 5(c)는 5000 sec 경과 후 하중 재하 방향(길이 방향)으로 위치별 발생한 변위를 나타내며, 고정 조건이 적용된 x = 0.0 m인 지점에서는 변위가 발생하지 않은 반면에 배수 조건이 적용되어 압력 감소가 가장 크게 발생한 x = 20.0 m인 지점에서는 압축 방향 변위가 가장 크게 발생한 것을 확인할 수 있다.
Fig. 5에서 제시된 바와 같이 개발된 OGSFLAC 시뮬레이터를 통해 해석된 결과는 해석적 해 및 TOUGH-FLAC을 통한 해석 결과와 잘 일치하는 것으로 나타났으며, 이를 통해 포화 상태에서 단상 유체의 수리-역학적 복합거동 해석에 대한 검증이 잘 이루어졌음을 확인하였다. 해석 속도에 있어서는 OGSFLAC 시뮬레이터를 활용한 해석이 단일 해석방법 및 순차적 해석방법을 적용한 OpenGeoSys에 비해 느린 것으로 나타났다. 그러나 OpenGeoSys에서 제공되는 역학적 구성 모델이 제한적이라는 점을 고려할 경우 OGSFLAC 시뮬레이터가 향후 다양한 현장 적용에 있어 효용성을 지닐 수 있다고 판단되며, 후속 연구를 통해 다양한 역학 구성 모델에 대한 적용 및 검증을 수행하고자 한다.
4. 결 론
본 연구에서는 수리-역학적 복합거동 해석의 필요성 증가에 부응하여 독자적인 해석 모듈을 개발하고자 수리 해석을 위한 OpenGeoSys 와 역학 해석을 위한 FLAC3D를 연동하여 시뮬레이터를 개발하였다. 개발된 시뮬레이터는 OGSFLAC으로 명명하였으며, 연동 방식으로 결과 파일 교환 방식에 기반한 순차적 해석 방법을 적용하였다. 연동 과정에서 OpenGeoSys가 주(master)로 작용하여 해석 시간을 결정하며 압력 정보를 전달하게 된다. FLAC3D는 종속(slave)으로 작용하여 OpenGeoSys에서 전달받은 압력 정보를 바탕으로 역학 해석을 수행하여 응력 및 변형률 정보를 전달하며, 이를 바탕으로 공극률 및 투수계수의 변화를 계산하여 반복적인 연동해석이 이루어지게 된다. 개발된 모듈은 그리드 생성, 다상 유체 해석 등에 있어 기존의 열-수리-역학적 복합거동 해석 모듈과 차이점을 보일 수 있을 것으로 판단된다.
개발된 모듈의 검증을 위하여 수학적 해가 존재하는 Terzaghi의 압밀 문제를 검증 모델로 선정하여 해석을 수행하였다. OGSFLAC을 통한 해석에 앞서 OpenGeoSys 및 FLAC3D와 비교하여 OGSFLAC이 갖는 장점을 확인하고자 개별 프로그램을 통한 해석을 수행하였다. FLAC3D는 입력 변수에 따라 해석 결과의 정확성 및 해석 속도가 민감하게 반응하는 것으로 나타났으나 OpenGeoSys는 안정적인 해석 결과를 보였으며, 이는 수리 해석 시 해석 기법에 의한 차이 때문인 것으로 판단된다. 해석 속도에 있어서는 OpenGeoSys의 단일(monolithic) 해석 방법이 가장 빠른 것으로 나타났으나 OpenGeoSys에서 제공하는 역학 구성 모델이 제한적이라는 점을 고려하면 향후 다양한 현장을 대상으로 한 해석에 있어 OGSFLAC 시뮬레이터가 효용성을 지니고 있다고 판단된다.
OGSFLAC 시뮬레이터를 통한 해석 결과, 시간에 따른 압력 및 변위 변화가 수학적 해와 잘 일치하는 것을 확인하였으며, 추가적으로 OGSFLAC과 같은 연동 방식을 사용하는 TOUGH-FLAC의 해석 결과와도 동일한 결과를 보임으로써 포화 상태에서의 단상 유체 거동에 대한 검증을 완료하였다. 그러나 본 연구에서 개발된 모듈을 통해 수행된 해석은 포화 상태에서 단상 유체 거동만을 고려한 압밀 문제이기 때문에 향후 다상 유체 거동, 모세관압에 의한 영향을 고려한 실제 규모에서의 검증 해석을 계획하고 있으며, 이를 통해 OGSFLAC의 단상 및 다상 유체 거동에 대한 수리-역학적 복합거동 해석에 대한 검증을 완료하고자 한다. 또한, 메쉬 생성, 해석, 결과 후처리를 포함한 오픈소스 프로그램 기반의 플랫폼을 구성하여 개발된 시뮬레이터를 효과적으로 활용하고자 계획하고 있으며, 이를 통해 다양한 국내외 현장에 대한 후속 연구를 수행할 수 있을 것으로 기대된다.
