Comparison of Numerical Response Characteristics of DEM Contact Formulations under Varying Timesteps: Spring Model and Impulse-based Method

Je-Kyum Lee; Myung-Kyu Song; Chang-Yong Kim; Wan-Kyu Yoo; Ji-Seok Yun; Sean Seungwon Lee

doi:10.7474/TUS.2026.36.3.257

Preview

Original Article

Tunnel and Underground Space. 30 June 2026. 257-276
https://doi.org/10.7474/TUS.2026.36.3.257

Comparison of Numerical Response Characteristics of DEM Contact Formulations under Varying Timesteps: Spring Model and Impulse-based Method

시간 간격 변화에 따른 DEM 접촉 해석 기법별 수치 응답 특성 비교: Spring model 및 Impulse-based method

Je-Kyum Lee¹

Myung-Kyu Song²

Chang-Yong Kim³

Wan-Kyu Yoo⁴

Ji-Seok Yun¹

Sean Seungwon Lee⁵^*

이 제겸¹

송 명규²

김 창용³

유 완규⁴

윤 지석¹

이 승원⁵^*

¹Postdoctoral Researcher, Dept. of Geotechnical Engineering Research, Korea Institute of Civil Engineering and Building Technology

²Research Professor, Department of Earth Resources and Environmental Engineering, Hanyang University

³Senior Research Fellow, Dept. of Geotechnical Engineering Research, Korea Institute of Civil Engineering and Building Technology

⁴Senior Researcher, Dept. of Geotechnical Engineering Research, Korea Institute of Civil Engineering and Building Technology

⁵Professor, Department of Earth Resources and Environmental Engineering, Hanyang University

¹한국건설기술연구원 지반연구본부 박사후연구원

²한양대학교 자원환경공학과 특임교수

³한국건설기술연구원 지반연구본부 선임연구위원

⁴한국건설기술연구원 지반연구본부 수석연구원

⁵한양대학교 자원환경공학과 교수

^{*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

The discrete element method (DEM) is widely used to simulate the mechanical behavior of discontinuous rock masses; however, the numerical stability, reproducibility of contact history, and computational efficiency of the analysis can vary depending on the contact calculation method adopted. In this study, two representative DEM contact formulations, the spring model and the impulse-based method, were comparatively analyzed in terms of their numerical response to variations in the timestep under identical simulation conditions. To this end, static contact, quasi-static shear, rotational, and dynamic impact conditions were considered, and the contact force, shear resistance, velocity, position, and energy responses were compared. The results show that the spring model effectively reproduces the gradual mobilization of contact force and shear resistance under static contact and quasi-static shear conditions, but exhibits contact force oscillation, distortion of the energy response, and divergence as the timestep increases. In contrast, the impulse-based method maintains a stable kinematic response even at relatively large timesteps under rotational and dynamic impact conditions, yet shows limitations in reproducing the quasi-static shear behavior that relies on contact history. The difference in the contact calculation schemes of the two methods manifests as a difference in their sensitivity to the timestep, which in turn determines the suitability of each method depending on whether the analysis focuses on contact history or momentum change. These findings suggest that, if future studies address the limitations of the impulse-based method, its application can be extended to jointed rock mass analysis while preserving its computational efficiency.

Keywords

개별요소법

불연속 암반

접촉 해석 기법

Spring model

Impulse-based method

불연속 암반의 거동을 모사하기 위해 개별요소법(Discrete Element Method, DEM)이 널리 활용되고 있으나, 접촉 응답을 계산하는 방식에 따라 해석 결과의 수치적 안정성, 접촉 이력 재현성 및 계산 효율성은 달라질 수 있다. 본 연구에서는 DEM의 대표적인 접촉 해석 기법인 Spring model과 Impulse-based method를 대상으로, 동일한 해석 조건에서 시간 간격 변화에 따른 수치 응답 특성을 비교·분석하였다. 이를 위해 정적 접촉, 준정적 전단, 회전 및 동적 충돌 조건을 설정하고, 접촉력, 전단저항, 속도, 위치 및 에너지 응답을 비교하였다. 해석 결과, Spring model은 정적 접촉 및 준정적 전단 조건에서 접촉력 변화와 전단저항의 점진적 발현을 재현하는 데 유리하였으나, 시간 간격 증가에 따라 접촉력 진동, 에너지 응답 왜곡 및 발산이 나타났다. 반면, Impulse-based method는 회전 및 동적 충돌 조건에서 상대적으로 큰 시간 간격에서도 안정적인 운동 응답을 보였으나, 접촉 이력에 기반한 준정적 전단 거동을 재현하는 데에는 제한이 있었다. 두 기법의 접촉 계산 방식의 차이는 시간 간격 변화에 대한 민감성의 차이로 나타났으며, 이는 해석의 초점이 접촉 이력에 있는지, 운동량 변화에 있는지에 따라 각 기법의 적합성을 가르는 요인이 된다. 향후 Impulse-based method의 한계를 보완하는 연구가 수행된다면, 계산 효율성을 유지하면서 절리 암반 해석으로 적용 범위를 확장할 수 있을 것으로 기대된다.

키워드

Discrete element method

Discontinuous rock mass

Contact formulation

Spring model

Impulse-based method

MAIN

1. 서 론
2. DEM 접촉 해석 기법의 이론적 배경
2.1. Spring model의 접촉력 계산
2.2. Spring model의 운동방정식 및 감쇠기법
2.3 Impulse-based method의 접촉력 계산
2.4 Impulse-based method의 운동방정식 및 감쇠기법
2.5 두 기법의 이론적 차이에 따른 역학적 모사 특성 비교
3. 수치 해석 환경 및 비교 기준
3.1 해석 프로그램 및 수치 실험 환경
3.2 대표 해석 조건 및 시간 간격 설정
3.3 절리 모델(Coulomb-Slip model)
4. 조건별 거동 특성 비교
4.1 정적 상태에서의 수렴성 비교
4.2 동적 충돌 시의 거동 및 시간 간격 민감도 분석
4.3 회전 거동 시의 수치 진동 특성 비교
4.4 시간 간격 민감도 및 접촉 응답 특성 고찰
5. 결 론

1. 서 론

암반은 절리, 단층 등 다양한 불연속면의 기하학적 배열과 역학적 특성에 의해 거동이 지배된다. 이러한 불연속 시스템에서는 블록의 미끄러짐, 회전, 분리와 같은 블록 단위의 거동이 연속체의 변형보다 지배적으로 나타날 수 있으며(Hoek and Brown, 1980), 연속체 기반의 수치해석 기법으로 이를 직접적으로 재현하는 데에는 한계가 있다. 이에 대한 대안으로 불연속체의 운동과 접촉 상호작용을 명시적으로 모사할 수 있는 개별요소법(Discrete Element Method, DEM)이 제안되었다(Cundall and Strack, 1979; Cundall, 1988; Hart et al., 1988). DEM은 암반을 독립적인 블록 또는 입자의 집합으로 모델링하고, 시간 간격(timestep)에 따라 블록 간 접촉 관계와 운동을 반복적으로 계산함으로써 불연속체의 거동을 해석한다. 이러한 특성으로 인해 DEM은 터널 굴착, 사면 안정성, 지하 구조물 주변의 블록 거동 해석 등 암반공학과 지반공학의 다양한 문제에 활용되어 왔다(Jing and Hudson, 2002; Lemos, 2012; Lisjak and Grasselli, 2014).

DEM 해석에서 접촉 응답을 계산하는 방식은 해석 결과의 정확도와 수치적 안정성을 결정하는 핵심 요소이며, 대표적으로 Spring model과 Impulse-based method가 있다. Spring model은 접촉을 블록 간의 지속적인 변형 및 상호작용 과정으로 해석하고, 접촉면에 가상의 스프링을 배치하여 상대 변위와 접촉 강성으로부터 접촉력을 계산하는 방법이다. 이 방식은 Cundall and Strack(1979)에 의해 DEM의 기본적인 접촉 처리 방식으로 제안된 이후, Coulomb-slip 및 Barton-Bandis 등 마찰, 점착, dilation, degradation 등 절리면의 역학적 특성을 고려하기 위한 다양한 절리 구성 모델과 결합되어 암반 절리 및 불연속체 거동 해석에 널리 활용되어 왔다(Barton and Choubey, 1977, Barton, 1982, Plesha, 1987). 반면, Impulse-based method는 접촉을 지속적인 변형 과정으로 해석하기보다는, 강체 간 충돌 과정에서 발생하는 운동량 교환을 통해 해석하는 방법이다. 이 방법은 접촉 시점의 상대 속도와 충돌 질량을 이용하여 충격량을 계산하고, 이를 통해 충돌 후 속도를 갱신한다(Hahn, 1988, Baraff, 1989, Mirtich, 1996). Impulse-based method는 본래 animation 및 computer graphics 분야에서 다수의 강체 충돌과 반발을 효율적으로 처리하기 위해 발전하였으며, 이후 토사와 같이 입자의 충돌 분리 및 재배열이 중요한 문제에서도 활용 가능성이 검토되어 왔다(Izadi and Bezuijen, 2018; He et al., 2020).

두 기법은 접촉 응답을 계산하기 위한 물리적 가정과 사용하는 물리량이 다르지만, 모두 블록의 위치와 속도를 시간 간격에 따라 갱신하는 명시적 시간 적분 과정 안에서 구현된다. 따라서 두 기법의 차이는 시간 간격이 변화할 때의 안정성, 정확도 및 계산 효율성 차이로 나타날 수 있다. 시간 간격을 작게 설정하면 접촉 발생과 블록 운동을 세밀하게 추적할 수 있으나 계산 cycle 수가 증가하고, 반대로 시간 간격을 크게 설정하면 계산 효율성은 향상될 수 있으나 접촉 발생 시점이나 접촉 중의 운동 변화를 충분히 반영하지 못할 가능성이 있다. 그러므로 두 접촉 해석 기법의 수치적 특성을 평가하기 위해서는 시간 간격을 주요 변수로 설정하고, 서로 다른 접촉 조건에서 나타나는 수치 응답의 차이를 함께 검토할 필요가 있다.

Spring model과 Impulse-based method는 서로 다른 해석 목적과 적용 분야를 중심으로 개별적으로 활용되어 왔으나, 최근에는 Impulse-based method를 지반공학 분야의 토사 거동 해석에 적용하기 위한 수치 실험적 연구가 수행된 바 있다(Lee and Hashash, 2015; Park et al., 2021). 이러한 연구들은 Impulse-based method가 다수 입자의 충돌, 분리 및 재배열이 빈번한 입상재 문제에서 계산 효율성 측면의 이점을 가질 수 있음을 보여준다. 그러나 절리암반에서는 단순한 블록 충돌뿐만 아니라 접촉 이력과 전단저항의 점진적 발현이 중요하므로, Impulse-based method를 암반공학 문제로 확장하기 위해서는 기존 Spring model과의 기본적인 응답 차이를 먼저 검토할 필요가 있다. 최근 두 기법의 응답 특성을 일부 비교한 연구가 보고된 바 있으나(Lee et al., 2025), 해당 연구는 통합 프레임워크의 제안과 검증에 초점을 두었으며, 두 접촉 해석 기법 자체의 시간 간격 민감도 및 접촉 조건별 응답에 대한 분석은 제한적으로 수행되었다.

이에 본 연구에서는 Spring model과 Impulse-based method의 접촉 응답 차이를 분리하여 검토하기 위해, 동일한 해석 조건에서 시간 간격을 변화시키며 정적 접촉, 준정적 전단, 회전 및 동적 충돌 조건에 대한 수치 실험을 수행하였다. 각 조건에서 나타나는 접촉력, 전단저항, 속도, 위치 및 에너지 응답을 비교하고, 시간 간격 변화에 따른 두 기법의 수치 응답 특성을 분석하였다. 이를 통해 두 접촉 해석 기법의 장점과 한계를 동일한 기준에서 정리하고, 불연속 암반 해석에서 접촉 해석 기법 선택과 시간 간격 설정 시 고려해야 할 주요 수치적 특성을 정리하고자 하였다.

2. DEM 접촉 해석 기법의 이론적 배경

본 장에서는 불연속 암반의 거동을 모사하기 위해 DEM에서 널리 활용되는 두 가지 접촉 해석 기법인 Spring model과 Impulse-based method의 이론적 배경을 고찰한다. 두 기법은 접촉을 해석하는 물리적 가정에 차이가 있으며, 이에 따라 접촉력의 형성, 절리 거동의 표현, 그리고 감쇠 처리 방식에서 서로 다른 수치적 특성을 나타낸다.

2.1. Spring model의 접촉력 계산

Spring model은 Cundall and Strack(1979)에 의해 처음 제안되었으며, 블록 간의 접촉을 Spring-damper 시스템으로 가정하여 해석한다. 두 블록이 접촉하여 겹침(Penetration depth)이 발생하는 경우, 식 (1)과 같이 상대 변위를 이용하여 접촉력을 계산하게 된다(Cundall, 1988, Itasca, 2000, Itasca, 2024). 매 시간 간격마다 상대 변위와 접촉 강성을 고려하여 힘의 증분을 계산하고, 이를 기존 접촉력에 누적함으로써 최종 접촉력을 결정한다(식 (2)-(3)). 이러한 방식은 겹침 정도에 따른 접촉력의 점진적인 변화를 모사할 수 있다는 장점이 있다.

(1)

∆ F = K \cdot ∆ u_{r e l} = K \cdot {\dot{u}}_{r e l} \cdot ∆ t

(2)

F : = F + ∆ F

(3)

M = u_{A r m} \times F

이때, $F$ 는 접촉력 벡터이며 $M$ 은 모멘트를 나타낸다. $∆ u_{r e l}$ 는 단일 시간 간격 동안 발생하는 두 블록 간의 상대 변위이며, ${\dot{u}}_{r e l}$ 는 충돌 시의 상대 속도이다. $u_{A r m}$ 은 접촉점에 대한 각 블록의 모멘트암, $K$ 는 접촉 강성, $∆ t$ 는 시간 간격을 의미한다. 기호 ‘ $: =$ ’는 대입 연산자를 나타낸다.

2.2. Spring model의 운동방정식 및 감쇠기법

접촉력이 산정된 이후에는 운동방정식에 기반한 시간 증분을 통해 블록의 가속도, 속도 및 변위를 갱신한다(식 (4), (5), (6), (7)). 이때 접촉 과정에서 발생하는 에너지 소산을 고려하고 해석의 수치적 안정성을 확보하기 위해 감쇠항이 도입된다. 본 연구에서는 Cundall이 제안한 Local damping 기법을 활용하였다(Cundall, 1987). 이 방식은 접촉력의 크기에 비례하여 감쇠력이 결정되므로 블록의 속도가 0에 가까워지는 준정적 조건에서도 충분한 감쇠 효과를 제공할 수 있어, 잔류 진동을 억제하고 정적 평형 상태로의 수렴을 효과적으로 유도할 수 있다.

(4)

{\dot{u}}_{B}^{(t + ∆ t / 2)} = {\dot{u}}_{B}^{(t - ∆ t / 2)} + \sum_{C o n t s .} (F_{C o n t .} - {(F_{d})}_{C o n t .}) \frac{∆ t}{m_{B}}

(5)

{\dot{θ}}_{B}^{(t + ∆ t / 2)} = {\dot{θ}}_{B}^{(t - ∆ t / 2)} + \sum_{C o n t s .} (M_{C o n t .} - {(M_{d})}_{C o n t .}) \frac{∆ t}{I_{B}}

(6)

{(F_{d})}_{i} = α |\sum_{C o n t s .} F_{C o n t .}| s g n ({\dot{u}}_{B}^{(t - ∆ t / 2)})

(7)

{(M_{d})}_{i} = α |\sum_{C o n t s .} M_{C o n t .}| s g n ({\dot{θ}}_{B}^{(t - ∆ t / 2)})

이때 ${\dot{u}}_{B}$ 와 ${\dot{θ}}_{B}$ 는 각각 블록 속도의 병진 성분과 회전 성분이다. $(t - ∆ t / 2)$ 와 $(t + ∆ t / 2)$ 는 단일 시간 간격 동안 값의 변화를 의미하며, $F_{C o n t .}$ 및 $M_{C o n t .}$ 는 한 블록에 발생된 모든 접촉점에서의 접촉력과 모멘트를 의미한다. $F_{d}$ 와 $M_{d}$ 는 접촉력과 모멘트에 대한 감쇠력이며, $m_{B}$ 와 $I_{B}$ 는 각각 블록의 질량과 관성모멘트이다. $α$ 는 감쇠계수이며, $s g n ()$ 은 방향성분만을 나타내는 부호함수이다.

2.3 Impulse-based method의 접촉력 계산

접촉력이 누적되는 Spring model과 달리 Impulse-based method는 접촉을 단일 시간 간격 내 종료되는 순간적인 충돌 사건으로 간주한다. 이 방법은 충격량-운동량 원리(Impulse-momentum principle)에 기반하며, 충돌 시 발생하는 충격량은 식 (8)-(9)와 같이 충돌 전 두 블록의 상대 속도, 반발 계수(Restitution coefficient), 및 충돌 질량(Collision mass)으로 정의된다(Hahn, 1988; Baraff, 1989; Baraff, 1992; Mirtich, 1996). 계산된 충격량을 이용하여 블록의 충돌 후 속도를 직접적으로 계산할 수 있다. 이와 같이 겹침 정도에 직접적으로 의존하지 않는 계산 방식은 시간 간격의 크기에 대한 민감도가 상대적으로 낮다. 이러한 특성으로 인해 고속 충돌이 지배적인 접촉 조건에서 우수한 수치적 안정성과 높은 연산 효율을 나타낸다.

(8)

J = J {\hat{n}}_{J} = m_{c} ∆ {\dot{u}}_{r e l} \cdot {\hat{n}}_{J} = - (R + 1) m_{c} {\dot{u}}_{r e l} \cdot {\hat{n}}_{J}

(9)

m_{c} = 1 / \{(\frac{1}{m_{B 1}} + \frac{\hat{n_{J}} (u_{B 1 \to C} \times \hat{n_{J}}) \times u_{B 1 \to C}}{I_{B 1}}) + (\frac{1}{m_{B 2}} + \frac{\hat{n_{J}} (u_{B 2 \to C} \times \hat{n_{J}}) \times u_{B 2 \to C}}{I_{B 2}})\}

이때, $J$ 는 충격량 벡터, $J$ 와 ${\hat{n}}_{J}$ 는 각각 충격량의 크기와 방향벡터이다. $R$ 은 반발계수(restitution coefficient)이고, $m_{c}$ 는 충돌질량(Collision mass)을 의미한다. $u_{B 1 \to C}$ 와 $u_{B 2 \to C}$ 는 각각 블록 B1과 B2의 질점으로부터 접촉점까지의 모멘트암(Moment arm)을 의미한다.

2.4 Impulse-based method의 운동방정식 및 감쇠기법

Impulse-based method에서는 접촉 시 발생하는 충격량을 통해 충돌 후 속도 변화를 직접 계산한다(식 (10)-(11)). 따라서 감쇠 또한 접촉력에 기반하기보다는 블록의 속도에 비례하는 형태로 적용하게 된다(식 (12)-(13)). 이 방식은 상대 속도가 큰 고속 충돌에서 과도한 속도 성분을 효과적으로 저감하여 응답을 안정화하는 데 유리하다. 그러나 감쇠의 크기가 속도에 의존하기 때문에, 시스템이 정적 평형 상태에 가까워져 블록의 속도가 작아질수록 감쇠 효과는 감소하게 된다. 따라서 준정적 조건에서는 잔류 진동을 충분히 억제하거나 평형상태로의 수렴을 촉진하는데 한계가 있을 수 있다.

(10)

{\dot{u}}_{B}^{+} = {\dot{u}}_{B}^{-} + \frac{J {\hat{n}}_{J}}{m_{B}}

(11)

{\dot{θ}}_{B}^{+} = {\dot{θ}}_{B}^{-} + \frac{J (u_{B \to C} \times {\hat{n}}_{J})}{I_{B}}

(12)

{\dot{u}}_{B}^{+} : = {\dot{u}}_{B}^{+} (1 - α_{d} ∆ t)

(13)

{\dot{θ}}_{B}^{+} : = {\dot{θ}}_{B}^{+} (1 - α_{d} ∆ t)

이때, ${\dot{u}}_{B}$ 와 ${\dot{θ}}_{B}$ 의 첨자 $-$ 와 $+$ 는 각각 충돌 전과 후의 속도를 의미한다. $α_{d}$ 는 속도기반 감쇠를 위한 감쇠계수(viscous global damping factor)이며 시스템의 고유진동수를 고려하여 계산된다( $α_{d} = 2 ξ_{g} ω_{g}$ , $ξ_{g}$ : global damping ratio, $ω_{g}$ : average natural frequency).

2.5 두 기법의 이론적 차이에 따른 역학적 모사 특성 비교

앞서 살펴본 바와 같이 Spring model과 Impulse-based method는 접촉을 해석하는 기본 가정에서 뚜렷한 차이가 존재한다. Spring model은 접촉을 여러 시간 간격에 걸쳐 지속되는 상호작용으로 간주하고, 상대 변위의 누적에 따라 접촉력을 갱신한다. 따라서 접촉력의 점진적 발현, 정적 평형으로의 수렴, 전단 항복 이전의 탄성 변형과 같이 이력 의존적인 접촉 거동을 모사하는 데 유리하다.

반면, Impulse-based method는 접촉을 단일 시간 간격 내에서 발생하는 순간적인 충돌 사건으로 처리하고 충격량-운동량 관계를 통해 충돌 후 속도를 갱신한다. 이 방식은 접촉 변위나 접촉력 이력을 세밀하게 추적하기보다는, 충돌 후 운동 상태를 안정적이고 효율적으로 계산하는 데 초점을 둔다. 따라서 고속 충돌, 낙석, 붕락과 같이 블록의 분리와 충돌이 지배적인 조건에서는 수치적 안정성과 계산 효율 측면에서 장점을 가지게 된다.

다만, 이러한 장점은 모든 해석 조건에서 동일하게 발현되는 것은 아니며, 접촉 조건과 시간 간격 설정에 따라 각 기법의 수치적 특성이 달라질 수 있다. Spring model은 접촉 강성에 민감하여, 시간 간격이 충분히 작지 않을 경우 과도한 겹침으로 인해 비정상적인 접촉력 증폭 및 수치 진동에 대한 위험성이 있다. 반면, Impulse-based method는 상기와 같은 겹침에 의한 발산에는 상대적으로 덜 민감하지만, 준정적 접촉 이력이나 전단 변위의 점진적 누적을 정밀하게 재현하는 데에는 한계가 있다. 또한 시간 간격이 과도하게 커질 경우 접촉이 누락되는 문제를 근본적으로 배제할 수 없다.

결국 두 접촉 해석 기법은 어느 하나가 보편적으로 우수하다기보다는, 해석하고자 하는 접촉 조건에 따라 서로 다른 적합성을 가진다고 볼 수 있다. 따라서 정적 접촉, 전단 미끄러짐, 회전 및 동적 충돌 조건에서 두 기법의 거동을 동일한 기준으로 비교하고, 특히 시간 간격 변화에 따른 정확도와 안정성의 차이를 정량적으로 검토할 필요가 있다. 본 연구에서는 이러한 관점에서 Spring model과 Impulse-based method의 접촉 조건별 거동 특성과 시간 간격 민감도를 비교·분석하였다.

3. 수치 해석 환경 및 비교 기준

3.1 해석 프로그램 및 수치 실험 환경

두 접촉 해석 기법을 동일한 조건에서 비교하기 위해 자체 개발한 C++ 기반의 2차원 DEM 해석 코드를 활용하였다. 해석 코드는 Fig. 1과 같이 크게 1) 블록 간의 접촉을 탐색하는 단계, 2) 접촉력 혹은 충격량을 계산하는 단계, 3) 도출된 힘 또는 속도 변화를 바탕으로 블록의 위치를 갱신하는 단계로 구성된다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-03/N0120360304/images/ksrm_2026_363_257_F1.jpg

Fig. 1.

Flowchart of the DEM computational procedure

본 연구에서는 접촉 시의 역학적 계산 및 속도 갱신 단계(2단계 및 3단계 일부)에서만 Spring model과 Impulse-based method를 선택적으로 변화시켰으며, 접촉 탐색 알고리즘을 포함한 이외의 계산 과정은 동일하게 유지하였다. 이를 통해 접촉 처리 방식 외의 영향을 배제하고, 물리적 가정과 수치적 구현 방식의 차이에 따른 해석 결과를 직접 비교할 수 있도록 하였다.

블록 자체의 변형이 해석 결과에 영향을 미치지 않도록 하기 위해 해석 대상은 내부 변형을 고려하지 않는 2차원 강체 블록으로 단순화하여 가정하였다. 모든 해석에서는 Table 1과 같이 블록의 밀도, 접촉 강성, 마찰각, 점착력, 반발계수 및 감쇠계수를 동일한 조건으로 설정하여 두 기법 간 비교의 일관성을 확보하였다. 한편 Impulse-based method에서 계산되는 충격량은 결과 비교를 위해 시간 간격으로 나눈 등가의 접촉력으로 환산하였다.

Table 1.

Model properties for numerical tests

Categories	Value
Categories	Static test	Rotation test	Dynamic test
Density ( $Kg / m^{3}$ )	$2, 700$
Joint normal stiffness ( $N / m$ )	$7 \times 10^{9}$
Joint shear stiffness ( $N / m$ )	$5 \times 10^{8}$
Friction angle ( $°$ )	$30$	30	15
Cohesion ( $Pa$ )	$1 \times 10^{4}$
Damping coefficient ( $α$ )	$0.8$	$0.8$	$0.6$
Restitution ( $R$ )	$0.1$	$0.5$	$0.5$

감쇠 관련 계수의 경우, 반발계수( $R$ )는 식 (8)에서와 같이 Impulse-based method 분석에 적용하였으며, 감쇠계수( $α$ )는 식 (6)-(7)에서와 같이 Spring model 분석에 적용되었다. 여기서 두 계수는 에너지 소산을 서로 다른 방식으로 제어하므로, 동적 충돌 조건에서 공정하게 비교하기 위해서는 충돌 전후의 에너지 변화가 서로 대응되도록 설정해야 한다. 특히 Local damping의 권장값인 0.8(Itasca, 2024)을 동적 충돌 조건에 그대로 적용할 경우 과도한 에너지 소산으로 인해 Spring model의 응답이 이론해에서 크게 벗어나고, 결과적으로 Impulse-based method와의 비교에서 불리한 기준이 적용된다. 이에 본 연구에서는 동적 시험에서 Spring model의 감쇠계수를 0.6으로 조정하여, Impulse-based method의 반발계수 0.5와 유사한 수준의 에너지 소산이 발생하도록 하였다. 이를 통해 동일한 기준의 충돌 응답을 기반으로 시간 간격 변화의 영향을 비교할 수 있도록 하였다.

3.2 대표 해석 조건 및 시간 간격 설정

Spring model과 Impulse-based method의 접촉 조건별 거동 특성을 비교하기 위해 불연속체 해석에서 발생할 수 있는 주요 접촉 유형을 단순화하여 정적 평형조건, 준정적 전단 조건, 회전 지배 조건, 동적 충돌 조건의 네 가지 해석 조건을 설정하였다.

구체적으로, 1) 정적 평형 조건에서는 단일 블록이 자중에 의해 고정된 바닥면과 접촉한 후 수직 평형 상태에 도달하는 과정을 모사하였다. 2) 준정적 전단 조건에서는 수직 평형이 이루어진 이후 수평 방향 외력을 부여하여 전단력 증가와 sliding 발생 과정을 비교하였다. 3) 회전 지배 조건에서는 하나의 접촉점을 중심으로 블록이 회전하는 상황을 모사하여 접촉점 주변의 수치 진동과 회전 거동의 재현성을 평가하였다. 4) 동적 충돌 조건에서는 블록의 낙하 및 반발 과정을 통해 동적 접촉 조건에서의 속도 변화와 에너지 소산 특성을 분석하였다.

두 접촉 해석 기법의 시간 간격 의존성을 정량적으로 비교하기 위해 시간 간격 민감도 해석을 수행하였다. 시간 간격은 식 (14)와 같이 Spring model이 안정적으로 수렴할 수 있는 임계 시간 간격(critical time step, $∆ t_{c r}$ )을 기준으로 하여, 0.1배에서 100배까지 혹은 경우에 따라 발산하기까지 증가시켜 적용하였다. 이 범위에서는 동일한 시간 간격 조건에서 두 기법의 정확도, 안정성 및 오차 증가 양상을 직접 비교하였다. 추가적으로 각 기법의 안정성 한계와 failure mode를 확인하기 위해 일부 조건에서는 시간 간격을 더 크게 증가시켰다.

(14)

∆ t = F R A C ∆ t_{c r} = 0.2 \sqrt{\frac{m_{\min}}{2 K_{\max}}}

3.3 절리 모델(Coulomb-Slip model)

암반 불연속면의 전단 거동은 절리 모델에 큰 영향을 받는다. 본 연구의 목적은 절리 구성모델의 차이가 아니라 접촉 계산 방식의 차이를 비교하는 것이므로, dilation, degradation 등 복잡한 이력 의존 거동은 고려하지 않고, 식 (15)와 같은 가장 기본적인 Coulomb-slip 모델을 적용하였다.

(15)

F_{\max} = |F_{N}| \tan φ + c A

이때, $F_{\max}$ 는 최대 전단 마찰력, $F_{N}$ 은 수직 접촉력, $φ$ 는 마찰각, $c$ 는 점착력, $A$ 는 접촉 면적이다. 접촉면 사이에 작용하는 전단력( $F_{S}$ )이 $F_{\max}$ 보다 작은 경우에는, 전단 강성에 비례하여 전단력이 증가하며, $F_{\max}$ 보다 큰 경우에는 전단력을 $F_{\max}$ 로 제한하여 미끄러짐을 모사한다.

Impulse-based method에서도 두 기법의 비교 조건을 유지하기 위해 동일한 Coulomb-slip 모델을 적용하였다. 다만, 접촉 응답이 충격량의 형태로 계산되므로, 식 (16)과 같이 양변에 시간 간격을 곱하여 최대 전단저항 충격량으로 환산하였다.

(16)

J_{\max} = F_{\max} ∆ t = (|F_{N}| \tan φ + c A) ∆ t = J_{N} \tan φ + c A ∆ t

$J_{\max}$ 는 최대 전단저항 충격량이며, $J_{N}$ 은 수직 충격량이다. 전단 충격량 $J_{S}$ 가 $J_{\max}$ 를 초과하는 경우, $J_{S}$ 를 $J_{\max}$ 로 제한하여 전단 미끄러짐을 모사하게 된다.

4. 조건별 거동 특성 비교

3장에서 설정한 수치 해석 조건을 바탕으로, 2차원 강체 블록의 운동을 병진 운동 성분과 회전 운동 성분으로 구분하여 시간 간격 변화에 따른 Spring model과 Impulse-based method의 거동 특성을 비교·분석하였다. 병진 운동에 대해서는 접촉이 지속적으로 유지되어 접촉력 이력과 절리 거동의 모사 정확도가 중요한 정적 접촉 조건과, 접촉 시간이 짧고 충돌 전후의 운동량 변화가 지배적인 동적 충돌 조건을 각각 검토하였다. 또한 회전 운동에 대해서는 제한된 접촉점을 중심으로 접촉력과 모멘트가 동시에 변화하는 조건을 설정하여, 병진 운동만으로는 확인하기 어려운 접촉력 진동 및 각운동 응답의 차이를 평가하였다.

각 조건에서 도출된 해석 결과는 해석해(Analytical solution)와 비교하여 정확도를 평가하였으며, 시간 간격 변화에 따른 정확도 변화는 기본적으로 평균 제곱근 오차(Root Mean Square Error, RMSE)를 활용하여 정량화하였다. 이때, Impulse-based method는 접촉 응답이 충격량의 형태로 계산되므로, Spring model의 접촉력과 직접 비교를 위해 산출된 충격량을 시간 간격으로 나눈 등가 접촉력으로 환산하여 비교에 사용하였다. 또한 각 조건에서 정의한 정량 지표를 통해 시간 간격 선정이 해석의 정확도와 안정성에 미치는 영향을 종합적으로 평가하였다.

4.1 정적 상태에서의 수렴성 비교

단일 블록이 자중에 의해 바닥면과 접촉하여 정적 평형 상태로 수렴하는 과정과, 이후 수평 방향 전단력을 부가하여 전단저항이 발현되는 과정을 해석하였다(Fig. 2). 이를 통해 접촉력의 점진적 발현과 평형 상태로의 수렴과정을 비교하였다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-03/N0120360304/images/ksrm_2026_363_257_F2.jpg

Fig. 2.

Model description for static contact

Fig. 3은 정적 접촉 조건에서 시간 간격 변화에 따른 두 기법의 수직력 및 전단력 발현을 비교한 것이다. Spring model의 수직력 응답(Fig. 3(a))은 $0.1 ~ 1 t_{c r}$ 와 같이 작은 시간 간격 조건에서는 초기 진동 이후 이론값에 수렴하는 경향을 보였다. 시간 간격이 증가할수록 초기 진동의 진폭이 커지고 수렴에 필요한 시간이 증가하였으며, 일정 수준 이상의 시간 간격( $8 t_{c r}$ )에서는 과도한 접촉 반력으로 인해 접촉 상태가 안정적으로 유지되지 못하는 거동이 나타났다. 이는 Spring model에서 접촉력이 단일 시간 간격 동안 발생한 상대 변위와 접촉 강성에 의해 계산되므로, 시간 간격이 커질수록 한 계산 단계에서 발생하는 겹침량과 이에 따른 접촉 반력이 증가하기 때문으로 판단된다. 반면, Impulse-based method의 수직 응답(Fig. 3(b))은 $0.1 ~ 10 t_{c r}$ 범위뿐만 아니라 비교적 큰 시간 간격( $50 t_{c r}$ 및 $100 t_{c r}$ )에서도 발산 없이 약 2 cycle 만에 수직 평형 반력에 도달하는 경향을 보였다. 다만, 시간 간격이 매우 커질 경우 접촉을 인식하는 시점이 지연되면서 수직 반력의 발현 시점에 차이가 발생하였다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-03/N0120360304/images/ksrm_2026_363_257_F3.jpg

Fig. 3.

Contact responses of static contact (For the impulse-based method, the contact forces represent equivalent contact forces converted from impulses, impulse divided by timestep)

전단저항 발현 과정에서는 두 기법의 차이가 보다 뚜렷하게 나타났다. Fig. 3(c)에서 Spring model은 접촉면을 따라 누적되는 상대 전단변위에 따라 전단력이 점진적으로 증가하며, Coulomb-slip 기준에 의해 최대 전단저항까지 도달하게 된다. 작은 시간 간격( $0.1 - 1 t_{c r}$ ) 조건에서는 이론해와 거의 일치하는 전단력-전단변위 관계를 보였으나, 시간 간격이 증가함에 따라 최대 전단저항에 도달하는 전단변위의 오차가 점차 누적되었다. 특히 $8 t_{c r}$ 조건까지는 전단저항의 점진적 발현 과정이 유지되었으나, 이론해와의 차이가 점차 증가하였으며, 그 이상의 시간 간격에서는 접촉력 계산이 불안정해져 발산 거동이 나타났다. 반면, Fig. 3(d)에서 Impulse-based method는 접촉 이력을 누적하지 않고 각 시간 간격에서의 상대 속도에 따라 전단 충격량을 계산하므로, 전단력이 전단변위에 따라 점진적으로 증가하기보다는 최대 전단저항이 즉시 발현되는 형태의 응답을 보였다. 시간 간격이 커질수록 전단저항 발생 시기가 크게 지연되어, 특히 $100 t_{c r}$ 조건에서는 최대 마찰력이 발생하기 이전에 이론값 대비 약 27배에 달하는 1.4mm의 전단 변위가 발생하였다. 이는 정적 전단조건에서 시간 간격 증가가 접촉 탐지 시점 및 전단 응답 발현 위치의 오차로 나타날 수 있음을 보여준다.

Table 2는 정적 접촉 조건에서 수직력 오차와 최대 마찰력 도달 시 전단변위의 오차를 정량적으로 비교한 것이다. 수직력의 경우 Simulation time 기준 0-0.03초 사이의 이론값 대비 RMSE (Root mean square error)를, 전단력의 경우 최대 마찰력이 발생한 시점의 전단 변위를 이론값과 비교한 오차를 평가 지표로 사용하였다. Spring model의 수직력 오차는 시간 간격이 $1 t_{c r}$ 이하로 작은 경우에도 Impulse-based method보다 상대적으로 크게 나타났는데, 이는 Fig. 3(a)에서 확인되는 초기 접촉력 진동의 영향이 오차 산정에 포함되었기 때문이다. 그러나 최대 마찰력 도달 시 전단변위 오차는 매우 작게 발생하였으며, 시간 간격이 증가함에 따라 점진적으로 증가한 뒤 $10.0 t_{c r}$ 조건에서는 발산하였다. 반면, Impulse-based method는 $0.1 - 10.0 t_{c r}$ 범위에서 수직력 오차가 비교적 작고 발산 없이 안정적인 응답을 보였으나, 전단변위 오차는 약 0.05 mm 수준으로 나타났다. 이는 최대 전단저항이 이론적으로 요구되는 전단변위까지 누적된 후 발현된 것이 아니라, 접촉이 인식된 직후 변위가 거의 0에 가까운 시점에 최대 전단저항이 발생하였음을 의미한다. 또한 $50.0 t_{c r}$ 및 $100.0 t_{c r}$ 조건에서는 접촉은 안정적으로 유지되었으나, 전단변위 오차는 크게 증가하는 등, 최대 전단저항이 발현되기까지의 전단 탄성 변형 과정은 재현되지 않았다.

Table 2.

Numerical accuracy comparison for static contact condition

$∆ t / t_{c r}$	Normal force RMSE (kN)		Shear displacement error at max friction (mm)
$∆ t / t_{c r}$	Spring model	Impulse-based method	Spring model	Impulse-based method
0.1	1.3966	0.0181	9.246E-05	0.05058
0.2	1.3706	0.0678	4.046E-05	0.05058
0.4	1.2888	0.0595	1.620E-04	0.05056
0.6	1.1906	0.0520	4.945E-04	0.05053
0.8	1.1723	0.0808	4.449E-04	0.05049
1.0	1.1380	0.0971	6.315E-04	0.05044
2.0	1.0740	0.2191	2.076E-03	0.05000
4.0	1.6122	0.3721	4.192E-03	0.04827
6.0	2.5974	0.5257	1.004E-02	0.04536
8.0	(Diverged)	0.6856	0.01181	0.04128
10.0		0.8644	(Diverged)	0.03608
50.0	-	4.2533	-	0.31080
100.0	-	7.3884	-	1.38870

4.2 동적 충돌 시의 거동 및 시간 간격 민감도 분석

낙석이나 붕괴와 같이 블록 간의 고속 충돌이 지배적인 동적 환경에서는 접촉 기법에 따라 수치적 안정성이 크게 좌우된다. 이를 확인하기 위해 Fig. 4와 같이 바닥면이 고정되어 있으며 단일 블록이 낙하하는 해석 모델을 구성하고, 감쇠 계수 및 반발 계수로 정의된 에너지 소산 조건하에서 시간 간격 변화에 따른 충돌 응답을 비교하였다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-03/N0120360304/images/ksrm_2026_363_257_F4.jpg

Fig. 4.

Model description for dynamic contact

Fig. 5는 시간 간격 변화에 따른 Spring model과 Impulse-based method의 총 역학적 에너지, 연직 방향 속도, 연직 방향 위치를 비교한 것이다. Spring model의 결과(Fig. 5(a), (c), (e))는 $0.1 - 1 t_{c r}$ 범위에서 이론해와 일치하는 안정적이고 일관된 거동을 보였으며, 충돌이 발생할 때마다 일정한 비율로 에너지가 소산되며 점진적으로 평형 상태에 도달하였다. 그러나 $2 t_{c r}$ 이상부터는 에너지 소산량과 충돌 후 반발 속도에서 차이가 발생하기 시작하여, $6 t_{c r}$ 에서는 충돌 후 블록이 비정상적으로 높게 반발하고 에너지가 과도하게 유지되는 현상이 관찰되었으며, $8 t_{c r}$ 에서는 거동이 수렴하지 않고 발산하였다.

반면 Impulse-based method는 $0.1 ~ 100 t_{c r}$ 범위 전반에 걸쳐 에너지 소산 패턴, 속도 이력, 위치 이력이 이론해와 거의 일치하였다(Fig. 5(b), (d), (f)). Spring model이 $1 t_{c r}$ 를 초과하면서 정확도가 저하된 것과 달리, Impulse-based method는 $100 t_{c r}$ 조건에서도 이론적 거동과 큰 차이 없는 안정적인 결과를 보였다. 이는 고속 충돌이나 낙하와 같이 순간적인 운동량 변화가 중요한 조건에서 Impulse-based method가 시간 간격 변화의 영향을 상대적으로 적게 받음을 보여준다. 다만, Fig. 5(f)에서 확인되는 것처럼, $50 t_{c r}$ 및 $100 t_{c r}$ 조건에서는 첫 충돌 시점이 미세하게 빨라지는 경향이 나타났다. 이는 시간 간격이 커질수록 접촉 탐지 시점의 오차가 누적된 결과로, 본 조건에서는 수치적 발산을 유발하지 않았으나, 접촉 발생 시점과 위치에 따른 물리적 정확도는 일부 저하될 수 있음을 보여준다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-03/N0120360304/images/ksrm_2026_363_257_F5.jpg

Fig. 5.

Velocity, position and energy change during collisions

Table 3은 처음 세 번의 충돌 이벤트를 대상으로 충돌 시점 오차와 에너지 소산 오차를 시간 간격 변화에 따라 비교한 것이다. 충돌 시점 오차는 각 충돌 이벤트에서 이론해의 충돌 발생 시각과 수치해석상의 충돌 인식 시각 사이의 오차를 의미하며, 에너지 소산 오차는 충돌 직전 및 직후의 운동에너지 변화량으로 산출한 에너지 소산값에 대한 이론값과 해석값 사이의 상대오차를 의미한다.

Table 3.

Collision timing error and energy dissipation error by timestep ratio

$∆ t / t_{c r}$	Timing error (ms)		Energy dissipation error (%)
$∆ t / t_{c r}$	Spring model	Impulse-based method	Spring model	Impulse-based method
0.1	1.7111	0.0179	0.0417	0.0061
0.2	1.8179	0.0389	0.2282	0.0254
0.4	3.7959	0.0787	1.1641	0.0426
0.6	11.1204	0.1289	4.0260	0.0922
0.8	4.5587	0.1603	1.2885	0.1022
1.0	3.6141	0.1824	2.1003	0.1254
2.0	31.0251	0.4173	8.8021	0.1763
4.0	6.4579	0.7203	29.2654	0.2440
6.0	123.7930	1.1829	94.0084	0.7637
8.0	(Diverged)	1.5340	(Diverged)	0.8628
10.0	-	1.8880	-	1.0753
50.0	-	10.2209	-	6.3973
100.0	-	20.0333	-	8.8442

두 기법 모두 $0.1 t_{c r}$ 조건에서는 에너지 소산 오차가 $0.006 - 0.04 %$ 수준으로 이론값과 잘 일치하였다. 그러나 시간 간격이 증가함에 따라 두 기법 간 뚜렷한 차이를 확인할 수 있었다. Impulse-based method는 $100 t_{c r}$ 조건에서도 에너지 소산 오차 $8.84 %$ 수준으로 오차가 비교적 완만하게 증가한 반면, Spring model은 $4 t_{c r}$ 에서 에너지 소산 오차가 $29.3 %$ 에 달하였으며, $8 t_{c r}$ 에서는 에너지가 오히려 증가하면서 발산하는 결과를 보였다. 충돌 시점 오차의 경우에도 두 기법 간 차이가 명확하게 나타났다. Impulse-based method의 충돌 시점 오차는 시간 간격에 비례하여 완만하게 증가하여, $100 t_{c r}$ 에서도 $20 ms$ 수준을 유지하였다. 반면, Spring model의 충돌 시점 오차는 불규칙하게 증가하였으며, $6 t_{c r}$ 에서 $123.8 ms$ , $8 t_{c r}$ 에서는 발산하면서 1초 이상의 오차가 발생하는 등 급격히 증가하였다.

4.3 회전 거동 시의 수치 진동 특성 비교

불연속 암반 사면의 붕괴나 지하 공간 천단부의 낙반 현상은 단순한 미끄러짐뿐만 아니라, 블록의 무게중심 편심에 의한 회전 거동이 복합적으로 발생한다. 회전이 지배적인 상황에서의 두 기법 간 차이를 비교하고 평가하기 위해, Fig. 6과 같이 하나의 모서리를 중심으로 하여 회전하는 삼각형 블록에 대한 해석을 수행하였다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-03/N0120360304/images/ksrm_2026_363_257_F6.jpg

Fig. 6.

Model description for rotational test

회전 시험에서 두 기법의 수치 응답 정확도를 평가하기 위한 기준으로, 고정된 한 모서리를 회전축으로 하는 강체의 중력 회전 운동에 대한 해석해를 사용하였다. (1) 접촉점이 운동 중 고정되어 있고, (2) 미끄러짐 없이 순수 회전이 발생하며, (3) 접촉이 분리 없이 유지되고, (4) 모든 블록이 강체이며, (5) 에너지 소산이 없다는 이상화된 가정 하에서, 회전각 $θ$ 에서의 각속도는 에너지 보존 법칙으로부터 식 (17)과 같이 유도된다.

(17)

ω (θ) = \sqrt{\frac{2 m g r (\sin θ_{0} - \sin θ)}{I_{O}}}

이때, $m$ 은 블록의 질량, $g$ 는 중력가속도, $r$ 은 회전축으로부터 무게중심까지의 거리, $θ_{0}$ 는 초기 회전각, $I_{O}$ 는 회전축에 대한 관성모멘트이다. 무게중심의 병진속도, 접촉점의 수직력 및 전단력은 모두 회전각의 함수로 유도된다. 이들 식은 표준적인 강체 동역학에 기반하므로, 본 연구에서는 생략하며 구체적인 유도 과정은 참고문헌에 제시되어 있다(Lee, 2026). 다만, 본 해석에서는 일시적인 접촉 분리가 발생하고 수치적 안정성을 위한 감쇠가 적용되므로, 이론식의 이상화 가정이 엄밀하게 유지되지 않는다. 따라서 해석해와의 완전한 일치를 기대하기보다는, 두 기법 간 회전 응답의 상대적 정확도와 안정성을 비교하기 위한 기준으로 활용하였다.

Fig. 7은 회전 조건에서 시간 간격의 변화에 따른 Spring model과 Impulse-based method의 수직력, 전단력, 병진 속도, 각속도에 대하여 비교한 것이다. Spring model의 결과(Fig. 7(a), (c), (e), (g))에서는 정적 및 동적 충돌 조건과 달리 $0.1 ~ 4 t_{c r}$ 의 비교적 작은 시간 간격에서도 수직력 및 전단력 전반에 걸쳐 심한 수치 진동이 관찰되었다. $6 t_{c r}$ 조건에서는 발산이 발생하였다. 이러한 거동은 회전 조건에서 블록이 단일 모서리를 중심으로 회전하면서 접촉 스프링이 매 시간 간격마다 압축과 이완을 반복하기 때문이다. 정적 또는 충돌 후 평형 상태와 같이 접촉이 안정적으로 유지되는 조건에서는 접촉 스프링이 일정한 상태를 유지할 수 있으나, 회전 조건에서는 접촉점 자체가 계속 미세하게 이동하면서 스프링의 변형 상태가 지속적으로 변화하게 된다. 그 결과 접촉력뿐만 아니라 접촉 모멘트도 불규칙하게 증폭되며, 작은 시간 간격에서도 수치 진동을 완전히 억제하기 어렵다.

Impulse-based method의 접촉 응답(Fig. 7(b), (d), (f), (h))은 모든 결과가 이론해와 유사하게 발생하여, 서로 겹쳐서 보이지 않도록 마커의 표시 간격을 조정하였으며, 마커의 간격이 계산값의 불연속이나 이상치를 의미하지 않는다. 응답은 $0.1 ~ 100 t_{c r}$ 범위 전반에 걸쳐 이론해와 비교적 일치하는 매끄러운 응답을 보였으며, 수치 진동 없이 안정적인 접촉력 및 운동 응답이 유지되었다. 다만, 시간 간격이 $50 t_{c r}$ 이상으로 증가함에 따라 수직력 및 전단력 이력에서 이론해와의 미세한 차이가 확인되었으며, 이는 접촉 탐지 시점의 이산화 오차가 접촉력 계산에 영향을 미치기 때문으로 판단된다. 속도 이력에서는 $0.1 ~ 100 t_{c r}$ 모든 범위에서 결과가 거의 일치하였으며, 이론해와도 잘 부합하였다.

https://cdn.apub.kr/journalsite/sites/ksrm/2026-036-03/N0120360304/images/ksrm_2026_363_257_F7.jpg

Fig. 7.

Comparison of rotational motion responses obtained by different DEM formulations (For the impulse-based method, the contact forces represent equivalent contact forces converted from impulses, impulse divided by timestep)

Table 4는 회전 조건에서 시간 간격 변화에 따른 두 기법의 수치적 정확도를 비교한 것이다. 각속도, 병진속도, 수직력, 전단력에 대해 이론해 대비 RMSE를 평가 지표로 사용하였다. Spring model의 각속도 RMSE는 $0.1 - 4 t_{c r}$ 범위에서 $0.936 - 1.032 rad / s$ 로, 시간 간격이 증가하더라도 거의 일정하게 유지되는 경향이 확인되었다. 수직력 및 전단력 RMSE 또한 동일한 범위에서 각각 $3.38 - 3.41 kN$ , $1.826 - 1.883 kN$ 으로 유사한 수준을 유지하였다. 이는 다른 해석 조건에서 시간 간격 증가에 따라 RMSE가 점진적으로 증가하는 것과는 다른 경향이다. 이러한 결과는 회전 조건에서 Spring model의 오차가 시간 간격 변화보다 수치 진동에 의한 영향이 더 크게 나타난 것으로, $0.1 - 4 t_{c r}$ 는 응답에 영향을 줄 정도로 유의미한 오차를 유발하지는 않은 것으로 판단된다.

Table 4.

Numerical accuracy comparison for rotational condition

$∆ t / t_{c r}$	Angular velocity RMSE (rad/s)		Translational velocity RMSE (m/s)		Normal Force RMSE (kN)		Shear Force RMSE (kN)
$∆ t / t_{c r}$	Spring model	Impulse-based method	Spring model	Impulse-based method	Spring model	Impulse-based method	Spring model	Impulse-based method
0.1	0.9359	0.1312	0.7453	0.0899	3.4099	0.4446	1.8832	0.1725
0.2	1.0310	0.1327	0.8163	0.0911	3.3814	0.4371	1.8297	0.1694
0.4	1.0311	0.1371	0.8164	0.0945	3.3803	0.4411	1.8295	0.1623
0.6	1.0311	0.1420	0.8164	0.0983	3.3793	0.4368	1.8292	0.1552
0.8	1.0312	0.1558	0.8164	0.1091	3.3782	0.4352	1.8289	0.1440
1.0	1.0312	0.1659	0.8165	0.1170	3.3772	0.4383	1.8287	0.1444
2.0	1.0315	0.1742	0.8167	0.1234	3.3729	0.4352	1.8278	0.1527
4.0	1.0321	0.1758	0.8172	0.1243	3.3663	0.4394	1.8259	0.1558
6.0	(Diverged)	0.1777	(Diverged)	0.1254	(Diverged)	0.4395	(Diverged)	0.1607
8.0	-	0.1794	-	0.1264	-	0.4357	-	0.1650
10.0	-	0.1813	-	0.1274	-	0.4375	-	0.1666
50.0	-	0.2144	-	0.1477	-	0.4460	-	0.2910
100.0	-	0.2582	-	0.1754	-	0.5784	-	0.4848

반면 Impulse-based method는 $100 t_{c r}$ 조건에서도 발산 없이 안정적으로 수렴하였으며, 각속도 RMSE는 $0.1 t_{c r}$ 의 $0.131 rad / s$ 에서 $100 t_{c r}$ 의 $0.258 rad / s$ 로 완만하게 증가하였다. $1 t_{c r}$ 조건에서 Impulse-based method의 각속도 RMSE는 $0.166 rad / s$ 로 Spring model의 $1.031 rad / s$ 대비 약 6배 낮은 수준이며, 수직력 및 전단력 RMSE 또한 각각 $0.438 kN$ , $0.144 kN$ 으로 Spring model 대비 낮은 수준의 오차를 보였다. 이러한 차이는 $100 t_{c r}$ 까지도 유지되어, Impulse-based method는 모든 조건에서 Spring model의 오차보다 낮은 수준을 보였다.

4.4 시간 간격 민감도 및 접촉 응답 특성 고찰

수치해석의 오차는 접촉 해석 기법과 시간 간격을 포함하여 블록의 형상과 질량, 접촉의 기하적 관계, 감쇠조건, 마찰조건, 초기 운동상태 및 하중 경로 등에 종합적으로 영향을 받는다. 따라서 본 연구에서 제시한 시간 간격 조건에서 나타난 결과는 모든 문제에 보편적으로 적용 가능한 안정성 한계나 정확도 기준으로 해석하기보다는, 본 시험에서 설정한 블록 형상, 초기 조건, 접촉 물성 및 경계조건 하에서 나타난 정성적 응답 특성으로 해석하는 것이 적절하다. 이에 따라 본 연구에서는 특정 시간 간격에서의 발산 여부나 오차 크기 자체보다는, 시간 간격 증가에 따라 각 접촉 해석 기법에서 나타나는 응답 변화의 양상과 그 원인을 해석하고자 하였다.

Spring model에서 시간 간격의 영향은 접촉 강성에 기반한 반력 계산과 직접적으로 관련된다. 시간 간격이 증가하면 단일 시간 간격 내에서 발생하는 겹침량 및 전단 변위 증분이 커지고, 이에 따라 스프링 반력이 과도하게 증가할 수 있다. 이러한 반력 증가는 접촉력 진동, 전단저항 발현 시점의 변화, 에너지 응답 및 감쇠 거동의 왜곡으로 이어질 수 있으며, 결과적으로 물리적 정확도를 저하시킬 수 있다. 따라서 Spring model에서 임계 시간 간격은 수치적 발산을 억제하기 위한 기본적인 안정성 기준으로 작용한다. 다만 다양한 접촉 조건에서 해석의 안정성과 접촉 응답의 재현성을 함께 확보하기 위해서는 보다 보수적인 시간 간격 설정이 필요하다. 특히 접촉 안정화 이후 상대 변위의 누적에 따라 접촉력이 변화하는 과정을 충분히 해상할 수 있도록 시간 간격을 설정하는 것이 중요하다.

반면, Impulse-based method는 접촉을 순간적인 충격량으로 처리하므로, Spring model에서 나타나는 반력의 과도한 증폭이나 그에 따른 발산 현상은 발생하지 않았다. 이에 따라 블록 낙하, 충돌 및 회전과 같이 접촉 시간이 짧고 운동량 변화가 지배적인 조건에서는 상대적으로 큰 시간 간격에서도 안정적인 운동 응답을 나타낼 수 있었다. 그러나 접촉 시점의 상대 속도에 기반하여 충격량을 계산하므로, 시간 간격이 커질수록 실제 접촉 발생 시점과 해석상 접촉 인식 시점 사이의 오차가 증가할 수 있다. 이 경우 블록은 실제보다 더 이동하거나 회전한 상태에서 접촉으로 인식되므로, 실제와 계산상 접촉 위치, 접촉 방향, 접촉면의 상대적 배치가 달라질 수 있으며, 이러한 차이는 충돌 시점과 이후 운동 경로의 오차로 이어질 수 있다.

또한 Impulse-based method는 접촉 변위 이력을 직접 누적하지 않기 때문에, 준정적 접촉 조건에서 주요하게 고려되어야 하는 경로 의존적 거동을 재현하는 데에는 한계가 있다. 예를 들어 절리면에서 수직력과 전단력이 하중 경로에 따라 연속적으로 변화하는 거동은 접촉 이력의 누적을 통해 표현된다. Spring model은 초기 수렴 과정에서 접촉력 진동이 발생할 수 있으나, 접촉이 안정화된 이후에는 상대 변위의 누적에 따라 수직력과 전단력이 변화하는 과정을 시간 이력으로 추적할 수 있다. 따라서 하중 경로에 따른 수직 반력의 변화와 이에 수반되는 전단저항의 점진적 발현을 모사하는 데 유리하다. 반면, Impulse-based method는 수직 반력을 비교적 빠르게 모사할 수 있더라도, 접촉 이력을 기반으로 한 전단저항 발현이나 수직-전단 응답의 연속적인 변화를 정밀하게 재현하는 데에는 제한이 있다.

따라서 두 기법의 적용성은 단순히 정적 문제와 동적 문제로 구분하기보다는, 해석에서 중요하게 평가하고자 하는 접촉 응답의 성격에 따라 판단할 필요가 있다. 접촉이 지속적으로 유지되고, 상대 변위의 누적에 따른 접촉력 변화, 전단저항의 점진적 발현, 응력 및 변형 이력이 주요 관심 대상인 경우에는 Spring model이 더 적합하다. 이 경우 시간 간격은 발산 여부만을 기준으로 설정하기보다는, 접촉력 진동과 전단저항 이력의 왜곡을 충분히 억제할 수 있도록 보수적으로 설정되어야 한다. 반면, 접촉 시간이 짧고 충돌 전후의 속도 변화, 반발 거동, 운동 경로와 같이 운동량 변화가 주요 관심 대상인 경우에는 Impulse-based method가 안정성 및 계산 효율성 측면에서 유리할 수 있다. 다만 Impulse-based method를 큰 시간 간격 조건에서 적용할 경우에는 발산 여부만으로 안정성을 판단하기보다는, 접촉 인식 시점, 접촉 위치 및 충돌 이후 운동 응답의 오차를 함께 검토할 필요가 있다.

5. 결 론

본 연구에서는 불연속 암반 거동 해석에 사용되는 Spring model과 Impulse-based method를 대상으로, 정적 접촉, 준정적 전단, 동적 충돌 및 회전 조건에서 시간 간격 변화에 따른 접촉 응답을 비교·분석하였다. Spring model은 겹침 또는 상대 변위와 접촉 강성을 이용하여 접촉력을 계산하는 반면, Impulse-based method는 접촉 시점의 상대 속도와 충돌 질량을 이용하여 운동량 변화를 계산한다. 수치 실험 결과, 이러한 접촉 응답 산정 방식의 차이는 시간 간격 증가에 따른 수치 안정성, 접촉 이력 재현성 및 접촉 인식 정확도의 차이로 나타났다. 본 연구의 주요 결론은 다음과 같다.

1) Spring model과 Impulse-based method는 접촉 응답 계산에 사용되는 물리량이 다르며, 이에 따라 시간 간격 변화에 대한 수치 응답도 서로 다른 양상으로 나타났다. Spring model은 접촉면의 상대 변위와 접촉 강성을 기반으로 접촉력을 누적하므로, 접촉이 지속되는 조건에서 접촉력의 변화와 전단저항의 발현 과정을 연속적으로 재현할 수 있었다. 그러나 시간 간격이 증가할수록 단일 계산 단계에서 발생하는 겹침량 및 전단 변위 증분이 커졌으며, 접촉력 진동, 수렴 지연, 에너지 응답 왜곡 및 발산이 발생하였다. 반면 Impulse-based method는 접촉 시점의 상대 속도에 기반하여 운동량 변화를 직접 계산하므로, 동적 충돌 및 회전 조건에서 상대적으로 큰 시간 간격에서도 안정적인 운동 응답을 보였다. 다만 접촉 이력을 누적하지 않는 계산 특성으로 인해, 전단저항이 점진적으로 발현되는 과정이나 이력 의존적인 접촉 거동을 정밀하게 재현하는 데에는 한계가 있었다.

2) 두 기법에서 시간 간격은 서로 다른 수치해석적 의미를 갖는다. Spring model에서 시간 간격은 접촉 강성에 의해 계산되는 반력의 안정성과 직접적으로 관련되며, 임계 시간 간격은 수치적 발산을 억제하기 위한 기본적인 안정성 기준으로 작용한다. 그러나 블록 간 충돌 속도가 커질수록 발산 위험성도 커지기 때문에, 접촉 응답의 재현성을 확보하기 위해서는 보다 보수적인 시간 간격 설정이 필요하다. 반면 Impulse-based method에서 시간 간격은 접촉 강성 기반의 안정성 조건이라기보다, 실제 접촉 발생 시점과 해석상 접촉 인식 시점 사이의 차이를 결정하는 시간적·기하학적 이산화 기준에 가깝다. 따라서 큰 시간 간격에서도 반력 증폭으로 인한 발산은 발생하지 않지만, 접촉 위치, 접촉 방향 및 충돌 이후 운동 경로의 오차는 발생할 위험이 있다.

3) 해석 기법의 적용성은 단순히 정적 문제와 동적 문제로 구분하기보다는, 해석에서 중요하게 평가하고자 하는 접촉 응답의 성격에 따라 판단할 필요가 있다. 접촉이 지속적으로 유지되고, 상대 변위의 누적에 따른 접촉력 변화, 전단저항의 점진적 발현, 응력 및 변형 이력이 주요 관심 대상인 경우에는 Spring model이 적합하다. 반면, 접촉 시간이 짧고 충돌 전후의 속도 변화, 반발 거동 및 운동 경로와 같이 운동량 변화가 주요 관심 대상인 경우에는 Impulse-based method가 안정성 및 계산 효율성 측면에서 유리할 수 있다. 따라서 두 기법은 어느 하나가 모든 조건에서 우수하다기보다는, 해석 대상에서 지배적인 접촉 응답과 요구되는 정확도에 따라 선택적으로 적용될 필요가 있다.

4) Impulse-based method의 절리암반 해석 적용성에 대해서는 추가적인 검토가 필요하다. 본 연구는 단순화된 블록 형상과 접촉 조건을 대상으로 두 접촉 해석 기법의 기본적인 수치 응답 특성을 비교한 것이므로, 그 결과를 실제 절리암반 조건으로 직접 일반화하기에는 한계가 있다. 본 연구에서 확인된 특성으로 볼 때, Impulse-based method는 접촉 응답을 누적하지 않기 때문에 절리암반에서 중요한 수직-전단 응답, 전단저항의 점진적 발현, dilation, 강도 저하 및 반복 전단에 따른 응력 재분배와 같은 경로 의존적 거동을 모사하는 데에는 불리할 것으로 예상된다. 그럼에도 불구하고 Impulse-based method는 큰 시간 간격 조건에서도 안정적인 운동 응답을 제공할 수 있으며, 계산 효율성 측면에서 상당한 이점을 기대할 수 있다. 따라서 향후에는 보다 실제적인 절리 모델을 대상으로, Impulse-based method를 절리암반 해석에 적용할 때 발생할 수 있는 한계를 체계적으로 평가하고, 이러한 한계를 완화하면서도 안정성 및 계산 효율성의 장점을 유지할 수 있는 방법을 연구할 필요가 있다.

Acknowledgements

본 연구는 국토교통부/국토교통과학기술진흥원의 지원으로 수행되었습니다(과제명 : 초장대 K-지하고속도로인프라 안전 및 효율 향상 기술 개발, 과제번호 RS-2024-00416524). 이에 감사드립니다.

References

Baraff, D., 1989, Analytical methods for dynamic simulation of non-penetrating rigid bodies, Computer Graphics, 23(3), 223-232.

10.1145/74334.74356

Baraff, D., 1992, Dynamic simulation of non-penetrating rigid bodies, Ph.D. Dissertation, Cornell University.

10.1145/133994.134084

Barton, N., and Choubey, V., 1977, The shear strength of rock joints in theory and practice, Rock Mechanics and Rock Engineering, 10, 1-54.

10.1007/BF01261801

Barton, N., 1982, Modelling rock joint behavior from in situ block tests: implications for nuclear waste repository design, ONWI-308, Terra Tek, Inc., Battelle Memorial Institute, Columbus, OH.

Cundall, P.A., and Strack, O.D.L., 1979, A discrete numerical model for granular assemblies, Géotechnique, 29(1), 47-65.

10.1680/geot.1979.29.1.47

Cundall, P.A., 1987, Distinct element models of rock and soil structure, in Analytical and Computational Methods in Engineering Rock Mechanics, Springer, Vienna, 129-163.

Cundall, P.A., 1988, Formulation of a three-dimensional distinct element model – Part I: A scheme to detect and represent contacts in a system composed of many polyhedral blocks, International Journal of Rock Mechanics and Mining Sciences, 25(3), 107-116.

10.1016/0148-9062(88)92293-0

Hahn, J.K., 1988, Realistic animation of rigid bodies, Computer Graphics, 22(4), 299-308.

10.1145/378456.378530

Hart, R., Cundall, P.A., and Lemos, J., 1988, Formulation of a three-dimensional distinct element model – Part II: Mechanical calculations for motion and interaction of a system composed of many polyhedral blocks, International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 25(3), 117-125.

10.1016/0148-9062(88)92294-2

He, H., Zheng, J., and Li, Z., 2020, Comparing realistic particle simulation using discrete element method and physics engine, Geo-Congress 2020, GSP 317, 464-472.

10.1061/9780784482803.050

Hoek, E., and Brown, E.T., 1980, Underground excavations in rock, Institution of Mining and Metallurgy, London.

Itasca Consulting Group (Itasca), 2000, UDEC – Universal Distinct Element Code: Theory and Background. Itasca Consulting Group, Inc. Minneapolis, MN.

Itasca Consulting Group (Itasca), 2024. 3DEC – Three-Dimensional Distinct Element Code, Version 7.0: Online Documentation. Itasca Consulting Group, Inc., Minneapolis, MN.

Izadi, E., and Bezuijen, A., 2018, Simulating direct shear tests with the Bullet physics library: A validation study, PLOS ONE, 13(4), e0195073.

10.1371/journal.pone.019507329672557PMC5908395

Jing, L., and Hudson, J.A., 2002, Numerical methods in rock mechanics, International Journal of Rock Mechanics and Mining Sciences, 39(4), 409-427.

10.1016/S1365-1609(02)00065-5

Lee, J.K., 2026, Development of a hybrid discrete element method with adaptive contact treatment approach for modeling discontinuous rock masses, Ph.D. Dissertation, Hanyang University.

10.1016/j.compgeo.2025.107527

Lee, J.K., Loy-Benitez, J., Song, M.K., and Lee, S.S., 2025, A hybrid discrete element method with adaptive contact treatment approach for modeling discontinuous rock masses, Computers and Geotechnics, 188, 107527.

10.1016/j.compgeo.2025.107527

Lee, S.J., and Hashash, Y.M.A., 2015, iDEM: an impulse-based discrete element method for fast granular dynamics, International Journal for Numerical Methods in Engineering, 104(2), 79-103.

10.1002/nme.4923

Lemos, J.V., 2012, Recent developments and future trends in distinct element methods: UDEC/3DEC and PFC codes, in Advances in Discontinuous Numerical Methods and Applications in Geomechanics and Geoengineering, Taylor and Francis Group, London, 17-29.

10.1201/b11600-4

Lisjak, A., and Grasselli, G., 2014, A review of discrete modeling techniques for fracturing process in discontinuous rock masses, Journal of Rock Mechanics and Geotechnical Engineering, 6(4), 301-314.

10.1016/j.jrmge.2013.12.007

Mirtich, B.V., 1996, Impulse-based dynamic simulation of rigid body systems, Ph.D. Dissertation, University of California, Berkeley.

Park, E.H., Kindratenko, V., and Hashash, Y.M.A., 2021, Shared memory parallelization for high-fidelity large-scale 3D polyhedral particle simulations, Computers and Geotechnics, 137, 104008.

10.1016/j.compgeo.2021.104008

Plesha, M.E., 1987, Constitutive models for rock discontinuities with dilatancy and surface degradation, International Journal for Numerical and Analytical Methods in Geomechanics, 11(4), 345-362.

10.1002/nag.1610110404

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

Preview

Comparison of Numerical Response Characteristics of DEM Contact Formulations under Varying Timesteps: Spring Model and Impulse-based Method

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

Fig. 1.

Flowchart of the DEM computational procedure

Table 1.

Model properties for numerical tests

(14)

(15)

(16)

Fig. 2.

Model description for static contact

Fig. 3.

Contact responses of static contact (For the impulse-based method, the contact forces represent equivalent contact forces converted from impulses, impulse divided by timestep)

Table 2.

Numerical accuracy comparison for static contact condition

Fig. 4.

Model description for dynamic contact

Fig. 5.

Velocity, position and energy change during collisions

Table 3.

Collision timing error and energy dissipation error by timestep ratio

Fig. 6.

Model description for rotational test

(17)

Fig. 7.

Comparison of rotational motion responses obtained by different DEM formulations (For the impulse-based method, the contact forces represent equivalent contact forces converted from impulses, impulse divided by timestep)

Table 4.

Numerical accuracy comparison for rotational condition

Acknowledgements

References