지질 및 자원 공학 분야에서 수압 균열(hydraulic fracturing)에 의한 균열 전파 메커니즘은 석유 및 천연가스 채굴, 지열 에너지 개발 등 다양한 응용에서 핵심적인 역할을 한다. 특히 드릴 홀 주변에서 발생하는 균열의 거동을 정확히 예측하는 것은 공학적 설계와 안정성 확보에 필수적이다. 이 문제를 효과적으로 다루기 위해 최근 주목받고 있는 방법 중 하나가 상장법(phase-field method)이며, 이를 다물리장 시뮬레이션 플랫폼인 COMSOL Multiphysics와 결합하면 보다 현실적인 암석 손상 해석이 가능하다.

상장법의 기본 개념

상장법은 금속이나 암석과 같은 고체 재료 내에서 균열의 발생 및 전파를 연속적으로 모델링하는 수치 기법이다. 기존의 균열 추적 방식과 달리, 균열의 경로를 사전에 정의하거나 동적으로 재격자화할 필요 없이, 균열 면 자체를 하나의 필드 변수로 처리한다. 이 변수를 흔히 \( d \)로 표기하며, 다음의 값을 가진다:

• \( d = 0 \): 무손상 매질 (균열 없음)
• \( d = 1 \): 완전 균열 영역

중간 값(\( 0 < d < 1 \))은 균열 전단부의 미세 손상 영역을 나타내며, 이는 정수형 균열 모델보다 물리적으로 더 자연스러운 전이를 제공한다. 상장 변수의 시간 발전은 다음과 같은 형태의 편미분 방정식으로 기술된다:

\[ \frac{\partial d}{\partial t} = L \left( -\frac{\partial \psi_0}{\partial d} + G_c \left( \frac{1}{\ell} d - \ell \nabla^2 d \right) \right) \]

여기서 \( L \)은 이동도(mobility), \( \psi_0 \)은 탄성 자유 에너지 밀도, \( G_c \)는 파괴 인성, \( \ell \)은 상장 전이 영역의 특성 길이를 각각 의미한다. 이 방정식은 균열이 에너지 최소화 원칙에 따라 자발적으로 전파되는 과정을 수치적으로 재현한다.

COMSOL을 통한 모델 구축

COMSOL은 사용자 정의 PDE(Partial Differential Equation) 기능과 내장된 구조역학 모듈을 통해 상장법 기반 균열 해석을 직관적으로 수행할 수 있다. 아래는 2차원 평면 변형 조건 하에서 드릴 홀 주변 균열 전파를 시뮬레이션하기 위한 일반적인 절차이다.

우선, 암석 매질 내에 위치한 원형 드릴 홀을 포함한 기하 모델을 생성한다. 예를 들어, 1m × 1m 크기의 사각 영역 중심에 반지름 0.05m의 원을 배치하고, 불리안 차감 연산을 통해 구멍을 형성한다. 이는 COMSOL의 GUI 기반 작업 흐름으로 쉽게 구현되며, 스크립트 기반 자동화도 지원된다.

// MATLAB 스타일의 COMSOL LiveLink 스크립트 예시
model.geom.create('geom1', 2);
model.geom('geom1').create('blk', 'Block');
model.geom('geom1').feature('blk').set('size', {'1[m]', '1[m]'});
model.geom('geom1').create('cir', 'Circle');
model.geom('geom1').feature('cir').set('radius', '0.05[m]');
model.geom('geom1').create('dif', 'Difference');
model.geom('geom1').feature('dif').selection.source({'blk'});
model.geom('geom1').feature('dif').selection.domains({'cir'});
model.geom('geom1').run();

기하 생성 후, "Solid Mechanics" 인터페이스와 함께 "Phase Field Fracture" 모듈을 추가하여 물리장을 연결한다. 상장 모델에서는 유효 응력을 손상 변수 \( d \)로 감쇠시키며, 일반적으로 다음과 같은 수정 응력-변형률 관계를 사용한다:

\[ \sigma = (1 - d)^2 \cdot \mathbb{C} : \varepsilon \]

여기서 \( \mathbb{C} \)는 탄성 계수 텐서이며, 제곱항은 균열 영역에서 강도가 급격히 감소하는 현상을 반영한다.

경계조건 및 수압 부하 적용

수압 균열 시뮬레이션의 핵심은 드릴 홀 내벽에 내압을 가하는 것이다. 이는 COMSOL에서 "Prescribed Pressure" 또는 "Boundary Load" 조건으로 구현할 수 있다. 예를 들어, 초기 5 MPa의 압력을 원형 경계에 단계적으로 적용하면, 국부적인 응력 집중이 발생하고, 특정 임계점을 초과할 경우 균열이 활성화된다.

model.physics('solid').feature.create('pr1', 'PrescribedDisplacement');
model.physics('solid').feature('pr1').selection.set('edge1'); // 드릴 홀 경계
model.physics('solid').feature('pr1').set('p', '5[MPa]');

또한, 대칭 조건을 활용해 계산 효율을 높이고, 시간 적분을 위해 적응형 시간 스텝을 설정함으로써 급격한 균열 전파 시에도 수렴성을 유지할 수 있다.

결과 해석 및 응용

시뮬레이션 결과를 통해 상장 변수 \( d \)의 공간 분포를 시각화하면, 균열의 시작 위치와 전파 경로를 명확히 확인할 수 있다. 일반적으로 드릴 홀 상하부에서 균열이 발생하여 수평 방향으로 전파되며, 이는 암석 내 최대 인장 응력 방향과 일치한다. 또한, 응력 클라우드 플롯을 통해 균열 선단부에서의 응력 집중 현상을 관찰할 수 있으며, 이는 전통적인 파괴 역학 이론과 잘 부합한다.

이러한 수치 해석은 실제 현장 조건을 반영하기 위해 비등방성 암석 특성, 다공성 유체 확산, 온도 변화 등의 요소와 결합될 수 있다. COMSOL의 다물리장 연성 기능을 활용하면, 유체-구조-상장 동시 해석도 가능해져 더욱 정교한 모델링이 실현된다.

결론적으로, 상장법과 COMSOL의 통합은 복잡한 균열 거동을 직관적이고 물리적으로 타당하게 재현할 수 있는 강력한 도구이다. 이 접근법은 암석 파손 메커니즘 이해뿐 아니라, 자원 회수 효율 극대화 및 지반 안정성 평가를 위한 핵심 시뮬레이션 기법으로 자리매김하고 있다.