Total Lagrangian SPH Numerical Framework for Large Deformation and Collision Simulations

Min Jun Kim; Gil Ho Yoon

doi:10.7734/COSEIK.2026.39.5.319

Preview

Research Paper

Journal of the Computational Structural Engineering Institute of Korea. 31 October 2026. 319-326
https://doi.org/10.7734/COSEIK.2026.39.5.319

Total Lagrangian SPH Numerical Framework for Large Deformation and Collision Simulations

대변형 및 접촉 시뮬레이션을 위한 Total Lagrangian SPH 수치 해석 프레임워크 구축

Min Jun Kim¹

Gil Ho Yoon²^*

김 민준¹

윤 길호²^*

¹Graduate Student, Department of Mechanical Engineering, Hanyang University, Seoul, 04763, Korea

²Professor, Department of Mechanical Engineering, Hanyang University, Seoul, 04763, Korea

¹한양대학교 융합기계공학과 석사과정

²한양대학교 기계공학부 교수

^{*Corresponding Author}

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

In this study, a Total Lagrangian Smoothed Particle Hydrodynamics (TL-SPH) framework was developed for collision and large-deformation problems of elastic solids. The total Lagrangian formulation was applied to stress evaluation and time integration, whereas contact forces evaluated in the current configuration were coupled to the reference-frame equations of motion. To handle collisions, a penalty contact model with a surface-aware penetration depth and Coulomb friction model were coupled to the TL-SPH equations of motion. The surface-aware measure computes the contact forces only when the surfaces penetrate, removing the non-physical repulsive force that conventional SPH produces from kernel overlap before contact. The framework was verified through three benchmarks—a 2D oscillating cantilever plate, an elastic cylinder on an inclined plane, and the collision of two rubber rings—verifying convergence, agreement with analytical solutions, and energy conservation.

Keywords

total lagrangian SPH

contact force model

large deformation

energy conservation

본 논문에서는 탄성체의 접촉과 대변형 거동을 해석하기 위한 Total Lagrangian SPH(TL-SPH) 기반 수치 해석 프레임워크를 제안한다. 응력 평가와 시간 적분에는 초기 기준 형상 기반 정식화를 적용하고, 현재 형상에서 평가한 접촉력을 이에 결합하였다. 접촉 처리는 표면 인식(surface-aware) 침투 깊이 기반의 페널티 접촉 모델과 쿨롱 마찰 모델로 수행한다. 이 표면 인식 방식은 표면이 침투할 때만 접촉력을 계산하여, 기존 SPH에서 접촉 이전의 커널 중첩으로 발생하던 비물리적 반발력을 제거한다. 2D 진동판, 경사면 접촉･마찰, 고무 링 충돌의 세 가지 벤치마크 문제를 통해 프레임워크의 정확성, 수렴성, 에너지 보존 특성을 검증하였다.

키워드

total lagrangian SPH

접촉력 모델

대변형

에너지 보존

MAIN

1. 서 론
2. 본 론
2.1 SPH 기법의 수학적 정식화
2.2 Total Lagrangian 운동량 보존 방정식
2.3 SPH 접촉 동역학
2.4 시간 적분
3. 수치 예제
3.1 2D 진동판 (2D Oscillating plate)
3.2 경사면 탄성체의 접촉･마찰 해석
3.3 고무 링의 대변형 충돌 해석
4. 결 론

1. 서 론

본 연구에서는 탄성체의 대변형 및 접촉 시뮬레이션을 위한 Total Lagrangian smoothed particle hydrodynamics(TL-SPH) 기반 수치 해석 프레임워크를 제안한다. 제안하는 프레임워크는 응력 평가와 시간 적분을 기준 형상에서 수행하며, 현재 형상에서 평가한 접촉력을 기준 형상 운동방정식에 결합한다. 이를 통해 대변형 해석의 수치적 안정성과 접촉 거동의 물리적 정확성을 확보하고자 한다.

대변형 및 고체 간 접촉 해석을 위해 무격자 방법인 SPH가 활용된다(Monaghan, 2005). 유체 해석에서 출발한 고체 SPH는 주로 Updated Lagrangian SPH(UL-SPH) 방식을 따랐으나, 이는 인장 응력 상태에서 입자 분포가 비물리적으로 군집되는 인장 불안정성(tensile instability)을 발생시켜 수치적 안정성을 저하시켰다(Swegle et al., 1995).

이를 극복하기 위한 접근으로 TL-SPH는 모든 연산을 초기 기준 형상(reference configuration)에서 수행한다. 이 방식은 인장 불안정성을 회피하여, 고체 대변형 해석의 표준 정식화로 자리잡았다(Belytschko et al., 2000). 최근에는 무격자 기법을 이용한 구조 동역학 해석 등 관련 연구가 이루어지고 있다(Lee et al., 2023). 반면 TL-SPH는 단일 물체의 대변형 해석에는 강건하나, 고체 간 접촉에서는 접촉 이전 단계의 커널 간섭이 비물리적 반발력을 일으킨다(Adepu and Ramachandran, 2024).

본 연구는 이 문제를 표면 인식(surface-aware) 접촉 모델로 해결한다. 침투 깊이 기반의 페널티 접촉 모델과 쿨롱 마찰 모델을 TL-SPH 운동 방정식에 결합하였다. 이 모델은 표면이 실제로 침투할 때만 접촉력을 계산하므로, 커널 중첩에 의한 비물리적 반발력이 나타나지 않는다. 이로써 대변형과 접촉을 단일 프레임워크에서 처리하였고, 세 가지 벤치마크로 정확성과 에너지 보존성을 검증하였다.

논문의 구성은 다음과 같다. 2장에서는 SPH 정식화와 TL-SPH 지배 방정식, 표면 인식 페널티 접촉력 모델, 시간 적분 기법을 서술한다. 3장에서는 2D 진동판, 경사면 접촉･마찰, 고무 링 충돌의 세 가지 벤치마크를 통해 프레임워크의 정확성, 강건성, 에너지 보존 특성을 검증한다. 4장에서는 결론과 향후 연구 방향을 제시한다.

2. 본 론

2.1 SPH 기법의 수학적 정식화

SPH는 연속체 변수 𝜙(X)를 입자들의 가중합으로 근사하는 무격자 방법이다. 위치 X에서의 함숫값은 커널 함수 $W$ 를 이용하여 다음 식 (1)과 같이 보간된다.

(1)

ϕ (X_{i}) \approx \sum_{j} ϕ (X_{j}) W_{ij} V_{j}^{0}

여기서, $V_{j}^{0}$ 는 이웃 입자 $j$ 의 초기 부피, $W_{i j}$ 는 입자 $i$ 와 $j$ 간의 초기 거리에 따른 커널 함숫값이다. 본 연구에서는 2차원 Wendland 5차 커널을 사용하였다(Wendland, 1995). 여기서 커널의 영향 반경을 결정하는 평활화 길이(smoothing length) $h$ 는 초기 입자 간격 $∆ p$ 에 대해 $h = 1.15 ∆ p$ 로 설정하였다.

한편, 자유 표면이나 경계 입자 부근에서는 표준 SPH 미분 연산자가 선형 변화를 정확하게 계산하지 못한다. 이를 보완하기 위해 기준 형상에서 정의되는 그래디언트 보정 행렬(gradient correction matrix, $B_{i}^{0}$ )(식 (2))을 도입한다(Bonet and Kulasegaram, 2002).

(2)

B_{i}^{0} = [\sum_{j} V_{j}^{0} (X_{j} - X_{i}) \otimes \nabla_{i}^{0} W_{ij}]^{- 1}

여기서, $\nabla_{i}^{0} W_{i j}$ 는 기준 좌표계에 대한 커널 함수의 미분이다. 이 보정 행렬을 커널 그래디언트에 곱하면 경계 부근에서도 변형 구배 텐서를 정확히 계산할 수 있다.

2.2 Total Lagrangian 운동량 보존 방정식

TL-SPH는 초기 기준 형상(𝛺₀)에서 정의된 입자 위치와 커널 그래디언트를 시뮬레이션 전 기간 동안 고정하고, 변형 정보를 변형 구배 텐서(Deformation gradient tensor, F)로 표현하는 정식화이다. 변형 구배 텐서 F는 현재 위치 x를 기준 형상 X에 대한 그래디언트( $\nabla_{0} x = \nabla_{0} u + I$ )로 정의된다. 여기서 u=x-X는 변위 벡터, $\nabla_{0}$ 는 기준 형상에 대한 그래디언트 연산자, 그리고 I는 단위 텐서이다. 변형 구배 텐서의 변화율은 입자 이산화를 통해 다음과 같이 계산된다.

(3)

{\dot{F}}_{i} = \sum_{j} V_{j}^{0} ({\dot{u}}_{j} - {\dot{u}}_{i}) \otimes \nabla_{i}^{0} W_{i j} \cdot B_{i}^{0}

여기서, $\dot{F}$ 는 변형 구배 텐서의 시간 변화율, $\dot{u}$ 는 입자의 속도 벡터이다. 변형 구배 텐서로부터의 응력 계산에는 대변형 탄성 거동에 적합한 Neo-Hookean 재료 모델을 사용한다. 부피 변형부와 편차 변형부를 분리한 변형 에너지 함수로부터 유도되는 Kirchhoff 응력 텐서는 다음 식 (4)와 같다.

(4)

τ = \frac{K}{2} (J^{2} - 1) I + Gdev (\bar{b})

여기서, $K$ 는 체적계수, $G$ 는 전단 계수, $J = \det (F)$ 이며, $\bar{b} = J^{- 2 / d} {FF}^{T}$ 는 modified left Cauchy-Green 변형 텐서, $d$ 는 공간 차원이다. 운동량 보존 방정식은 이로부터 얻어지는 제1 Piola-Kirchhoff 응력 텐서( $P = τ F^{- T}$ )를 통하여 표현되며, 이를 이산화한 기준 지배 방정식은 다음과 같다.

(5)

ρ_{0} {\ddot{u}}_{r, i} = \sum_{j} (P_{i} B_{i}^{0^{T}} + P_{j} B_{j}^{0^{T}}) \nabla_{i}^{0} W_{ij} V_{j}^{0}

여기서, 𝜌₀는 초기 밀도이다. 식 (5)의 표준 정식화에서는 대변형 시 입자가 지그재그 패턴으로 엇갈리는 아워글래스 모드(hourglass mode)가 발생하여 수치적 불안정성을 유발한다. 이를 억제하기 위해 표준 가속도 ${\ddot{u}}_{r, i}$ 에 전단 보정 가속도 ${\ddot{u}}_{s, i}$ 를 추가하며 다음 식 (6)과 같이 계산된다(Wu et al., 2025).

(6)

ρ_{0} u ¨_{s, i} = \sum_{j} ({\tilde{G}}_{i} F_{i}^{- T} B_{i}^{0^{T}} + {\tilde{G}}_{j} F_{j}^{- T} B_{j}^{0^{T}}) \frac{\partial W_{i j}}{\partial r_{i j}^{0}} V_{j}^{0} e^_{i j}^{0}

여기서, $\partial W_{i j} / \partial r_{i j}^{0}$ 는 기준 형상에서의 입자 간 거리 $r_{i j}^{0}$ 에 대한 커널 함수의 도함수이다. $\tilde{G}$ 는 전단 거동과 관련된 재료 강성 계수이고, 보정 벡터 ${\hat{e}}_{i j}^{0}$ 는 변형 구배로 기술되지 않는 입자 간 비물리적 상대 변위, 즉 아워글래스 오차를 나타내며, 두 항의 구체적 정의는 Wu 등(2025)을 따른다. 최종 가속도는 표준 항과 전단 보정 항의 합으로 구한다.

2.3 SPH 접촉 동역학

TL-SPH는 단일 물체 해석에는 효과적이나 접촉을 다루려면 별도의 알고리즘이 필요하다. 이를 해결하기 위해 본 연구에서는 표면 인식 페널티 접촉력 모델(Adepu and Ramachandran, 2024)을 도입하여 대변형 해석을 위한 TL-SPH 운동 방정식에 결합한다.

2.3.1 수직 접촉력 모델

두 물체가 접촉할 때, 한 쪽을 주(primary) 물체로 다른 쪽을 부(secondary) 물체로 지정한다. 주 물체의 입자 $a$ 가 부 물체 $b$ 로부터 받는 수직 접촉력은 침투 깊이에 비례하는 페널티 형태로 계산된다(식 (7)).

(7)

F_{a}^{N} = k_{r} δ_{n^{c}}^{a} n_{a}^{c}

여기서, $k_{r}$ 은 수직 강성 계수(normal spring stiffness coefficient), $δ_{n^{c}}^{a}$ 는 입자 $a$ 와 부 물체 표면 사이의 침투 깊이, 그리고 $n_{a}^{c}$ 는 접촉면 단위 법선 벡터이다. 이 모델의 핵심은 침투 깊이를 입자 단위가 아닌 표면 단위로 정의하는 데 있다. Fig. 1은 이러한 표면 단위 접촉의 기하를 나타낸 것으로, 부 물체가 다수의 입자로 이산화되어 있음에도 주 물체 입자 $a$ 가 개별 입자가 아닌 가중평균화로부터 매끄러운 표면 전체와 상호작용함을 보인다. 표면 단위의 침투 깊이는 부 물체 입자들의 가중 평균을 이용하여 다음 식 (8)과 같이 계산된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-05/N0040390505/images/Figure_jcoseik_39_05_05_F1.jpg

Fig. 1.

Schematic of the surface-aware contact model

(8)

d_{a} = \frac{\sum_{b} (n_{a}^{c} \cdot r_{ab}) \frac{m_{b}}{ρ_{b}} W_{ab}}{\sum_{b} \frac{m_{b}}{ρ_{b}} W_{a b}}

여기서, r_a는 입자 $a$ 의 현재 위치 벡터이고 r_ab=r_a-r_b이며, $m_{b}$ 와 $ρ_{b}$ 는 각각 부 물체 입자 $b$ 의 질량과 밀도, $W_{a b}$ 는 입자 $a$ 와 $b$ 사이의 커널 함숫값이며, 현재 형상의 입자 간 거리로 평가한다. Fig. 1에서 $d_{a}$ 는 입자 $a$ 로부터 접촉 법선 $n_{a}^{c}$ 방향을 따라 부 물체 표면까지 측정한 거리이며, 침투 깊이 $δ_{n^{c}}^{a}$ 는 초기 입자 간격 $∆ p$ 에서 $d_{a}$ 를 뺀 값으로 정의된다. 접촉 법선 벡터는 인접 부 입자들의 위치 정보로부터 SPH 가중평균과 정규화 과정을 거쳐 계산한다(식 (9)).

(9)

{\hat{n}}_{a}^{c} = \frac{\sum_{b} \frac{r_{ab}}{‖ r_{ab} ‖} \frac{m_{b}}{ρ_{b}} W_{ab}}{\sum_{b} \frac{m_{b}}{ρ_{b}} W_{ab}}

이를 정규화하여 단위 법선 벡터 $n_{a}^{c}$ 를 얻는다. Fig. 1과 같이, $n_{a}^{c}$ 는 인접 부 입자들이 주 입자 $a$ 를 향하는 방향 벡터 $r_{ab}$ 들의 가중평균으로 결정되어 부 물체 표면에 수직한 방향을 가리킨다. 이 표면 인식 방식은 실제 표면이 침투할 때만 정확한 접촉력을 계산하므로, 커널 반경의 중첩만으로 생기던 비물리적 반발력을 막을 수 있다.

2.3.2 마찰을 고려한 접선 접촉력 모델

접촉면의 마찰은 입자와 상대 물체 표면 사이에 가상의 접선 스프링을 두어 표현한다. 접선력은 접촉 이후 누적된 미끄러짐 거리에 비례하므로, 매 단계 접선 방향의 상대 속도를 시간 적분하여 스프링의 총 변형량으로부터 시험(trial) 접선력을 계산한다(식 (10)).

(10)

{\tilde{F}}_{a}^{t, n + 1} = - k_{f} [(∆ l_{a}^{n} + v_{ab}^{n + 1} ∆ t) \cdot t_{a}^{n + 1}] t_{a}^{n + 1}

여기서, $k_{f}$ 는 접선 강성 계수(tangential spring stiffness), $v_{a b}$ 는 입자 $a$ 가 부 물체에 대해 갖는 상대 속도이다. $∆ l_{a}$ 은 이전 단계까지 누적된 접선 스프링 변형이다. 또한, $∆ t$ 는 시간 증분을 나타낸다. 접선 단위 벡터는 다음 식 (11)과 같다.

(11)

t_{a} = \frac{(v_{ab} - (v_{ab} \cdot n_{a}^{c}) n_{a}^{c})}{‖ v_{ab} - (v_{ab} \cdot n_{a}^{c}) n_{a}^{c} ‖}

쿨롱 마찰 법칙(Coulomb friction law)에 따라 마찰력의 최댓값은 수직력에 비례한다(식 (12)).

(12)

F_{a}^{t, n + 1} = \min (μ ‖ F_{a}^{N} ‖, ‖ {\tilde{F}}_{a}^{t, n + 1} ‖) \cdot ({\tilde{F}}_{a}^{t, n + 1} ∕ ‖ {\tilde{F}}_{a}^{t, n + 1} ‖)

여기서, 𝜇는 마찰 계수이다. 이 방법은 정지(stick)와 미끄러짐(slip) 상태를 자동으로 구분하며, 미끄러짐이 시작되면 접선 스프링의 누적된 변형을 마찰력 최댓값까지만 허용하여 마찰 거동을 재현한다. 입자 a에 작용하는 총 접촉력은 수직력과 접선력의 합으로 구해지며, 식 (5)의 우변에 외력으로 더해진다.

이 접촉력은 작용-반작용의 법칙(Newton’s third law)에 따라 부(secondary) 물체에서 가장 가까운 입자에 같은 크기, 반대 방향으로 가해진다.

2.4 시간 적분

본 연구에서는 에너지 손실을 줄이기 위해 위치 기반 Verlet(position-based Verlet) 시간 적분 기법을 사용한다. 먼저 $n$ 번째 단계에서 가속도로 변형 구배 텐서, 밀도, 입자 위치를 중간 시각으로 갱신한다. 여기서 r는 입자의 현재 위치 벡터이다(식 (13)).

(13)

\{\begin{cases} F^{n + 1 / 2} = F^{n} + (\frac{∆ t}{2}) {\dot{F}}^{n} \\ ρ^{n + 1 / 2} = ρ_{0} \frac{1}{J} \\ r^{n + 1 / 2} = r^{n} + (\frac{∆ t}{2}) {\dot{r}}^{n} \end{cases}

갱신된 중간 형상에서 2.2절의 방법으로 Kirchhoff 응력을 구하고, 이로부터 가속도 ${\ddot{u}}^{n + 1 / 2}$ 를 계산한다. 이 가속도로 입자 속도를 갱신한다(식 (14)).

(14)

{\dot{u}}^{n + 1 / 2} = {\dot{u}}^{n} + (∆ t) {\ddot{u}}^{n + 1 / 2}

갱신된 속도로 식 (3)의 변형 구배 변화율 ${\dot{F}}^{n + 1}$ 을 재평가한 뒤, 다음 단계( $n + 1$ )의 변형 구배 텐서, 위치, 밀도를 확정한다(식 (15)).

(15)

\{\begin{cases} F^{n + 1} = F^{n + 1 / 2} + (\frac{∆ t}{2}) {\dot{F}}^{n + 1} \\ ρ^{n + 1} = ρ_{0} \frac{1}{J} \\ r^{n + 1} = r^{n + 1 / 2} + (\frac{∆ t}{2}) {\dot{u}}^{n + 1} \end{cases}

안정적인 계산을 위해 시간 간격은 CFL(Courant-Friedrichs-Lewy) 조건으로 정한다(식 (16)).

(16)

∆ t = CFL \cdot \min (\frac{h}{c_{v} + | \dot{u} |_{\max}}, \sqrt{\frac{h}{| \ddot{u} |_{\max}}})

여기서, $| \dot{u} |_{\max}$ , $| \ddot{u} |_{\max}$ 는 전체 입자 중 최대 속도와 최대 가속도의 크기, $c_{v} = \sqrt{K / ρ}$ 는 체적 음속이다. 모든 예제에서 CFL=0.2를 사용하였다.

3. 수치 예제

본 장에서는 세 가지 벤치마크 문제를 통해 제안된 TL-SPH 기반 접촉 및 대변형 해석 프레임워크의 정확성, 강건성, 그리고 에너지 보존 특성 등을 정량적으로 검증한다. 모든 수치 예제에서 커널과 평활화 길이(smoothing length)는 2.1절과 동일하게 Wendland 커널과 $h = 1.15 ∆ p$ 를 사용하였다.

3.1 2D 진동판 (2D Oscillating plate)

본 절에서는 동적 거동 해석의 정확성과 공간 해상도(Spatial resolution)에 따른 수렴성을 검증하기 위해, Fig. 2와 같이 좌측 끝단에 고정(fixed) 조건을 부여한 형태의 2D 진동판 문제를 해석한다. 해석에 적용된 대상의 물성치는 Table 1과 같으며, Gray 등(2001) 이래 널리 쓰이는 표준 벤치마크 설정을 따른다. 자유 진동을 위해 전체 영역에 식 (17)의 초기 속도 분포를 적용하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-05/N0040390505/images/Figure_jcoseik_39_05_05_F2.jpg

Fig. 2.

Initial configuration of the 2D oscillating plate

Table 1.

Material properties of the 2D oscillating plate

Property	Symbol	Value
Length	$L$	0.2m
Thickness	$H$	0.02m
Density	𝜌₀	1000.0kg/m³
Young’s modulus	$E$	2.0MPa
Poisson’s ratio	𝜈	0.4
Initial velocity factor	$v_{f}$	0.05

(17)

v_{y} (x) = v_{f} c_{v} [\frac{f (x)}{f (L)}]

여기서, $v_{f} = 0.05$ =0.05는 초기 속도 계수, $c_{v}$ 는 체적 음속이며, $f_{(x)}$ 는 다음 식 (18)과 같이 주어진다.

(18)

f (x) = (\sin (k L) + \sinh (k L)) (\cos (k x) - \cosh (k x)) - (\cos (k L) + \cosh (k L)) (\sin (k x) - \sinh (k x))

여기서, $k$ 는 1차 굽힘 모드의 파수(wavenumber)로 $k L = 1.875$ =1.875를 만족하며, 이론 진동 주파수는 다음과 같이 주어진다.

(19)

ω^{2} = \frac{E H^{2} k^{4}}{12 ρ (1 - ν^{2})}

해석 결과 검증에 앞서 TL-SPH 정식화의 적용 효과를 분석한다. 공간 해상도 $H / ∆ p = 10, 20, 40$ =10, 20, 40에 대해 관측점(observation point)의 수직 방향 변위 시간 이력을 Fig. 3에 도시하였다. 해상도가 증가함에 따라 진동의 위상과 진폭이 수렴하며, $H / ∆ p$ =40에서 계산된 진동 주기는 식 (19)의 이론해와 약 5% 오차로 근접하였다. 또한, 장시간 모사에서도 진폭의 저하없이 거동이 안정적으로 유지되는 점을 통해 TL-SPH 기법의 정밀도 및 강건성을 확인하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-05/N0040390505/images/Figure_jcoseik_39_05_05_F3.jpg

Fig. 3.

Time histories of the vertical displacement at the observation point for three spatial resolutions ( $H / ∆ p = 10, 20, 40$ =10, 20, 40), compared with the analytical solution

에너지 보존 특성을 평가하기 위해 운동 에너지(Kinetic energy), 변형 에너지(Strain energy) 및 총 에너지(Total energy)의 시간 이력을 Fig. 4에 비교하였다. $H / ∆ p = 40$ =40에서 총 에너지 변동이 초기 값 대비 1% 이내로 유지되어 에너지 보존 성능을 보였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-05/N0040390505/images/Figure_jcoseik_39_05_05_F4.jpg

Fig. 4.

Time histories of total, kinetic, and strain energies ( $H / ∆ p = 40$ =40)

3.2 경사면 탄성체의 접촉･마찰 해석

마찰 모델을 검증하기 위해 Fig. 5와 같이 경사각 60°의 경사면 위에 놓인 반경 0.5m의 탄성 원통(elastic cylinder)이 자중에 의해 경사면을 따라 내려가는 문제를 해석한다. 재료 물성은 대변형 탄성 벤치마크에서 널리 쓰이는 고무 물성(Wu et al., 2025)을 사용하였고, 수치 파라미터와 함께 Table 2에 정리하였다. 페널티 강성을 포함한 수치 파라미터는 Adepu와 Ramachandran(2024)을 참고하여 선정하였다. 선정된 강성은 접촉 과정에서 최대 침투 깊이가 초기 입자 간격 $∆ p$ 보다 충분히 작게 유지되도록 하는 값이다. 두 마찰 계수(𝜇=0.3, 𝜇=0.6)에 대해 미끄러짐(slip)과 구름(rolling)으로 분기되는 양상을 평가하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-05/N0040390505/images/Figure_jcoseik_39_05_05_F5.jpg

Fig. 5.

Sliding elastic cylinder on an inclined plane (𝜃=60°)

Table 2.

Parameters used in the sliding-on-slope problem

Property	Symbol	Value
Young’s modulus	$E$	10MPa
Poisson’s ratio	𝜈	0.4
Density	𝜌₀	1200kg/m³
Gravity	$g$	9.81m/s²
Slope angle	𝜃	60°
Friction coefficient	𝜇	0.3, 0.6
Normal stiffness	$k_{r}$	1 × 10⁷N/m
Tangential stiffness	$k_{f}$	1 × 10⁵N/m

경사면 위 원통은 마찰력이 부족하면 미끄러지며 내려가고, 충분하면 미끄러짐 없이 굴러 내려간다. 두 거동의 전이는 경사각과 마찰계수의 관계로 정해지며, 알짜 외력만 가속에 기여한다는 가정에서 질량중심의 경사 방향 변위에 대한 이론해는 다음과 같이 분기 형태로 주어진다.

(20)

x_{cm} (t) = \{\begin{cases} x_{0} + \frac{1}{2} g t^{2} (\sin θ - μ c o s θ, & \tan θ > 3 μ \\ x_{0} + \frac{1}{3} g t^{2} \sin θ & \tan θ \leq 3 μ \end{cases}

여기서, $x_{cm}$ 은 실린더 질량중심의 경사 방향 위치, $x_{0}$ 는 초기 위치이며 본 해석에서는 $x_{0} = 0$ =0으로 둔다. 식 (20)에서 미끄러짐(slip) 분기는 마찰이 부족한 경우로, 마찰 항 𝜇cos𝜃가 감속에 기여한다. 구름(rolling) 분기는 마찰이 충분해 굴러 내려가는 거동으로, 회전 관성에 의한 감속이 계수 1/3에 반영된다. 𝜃=60°에서는 𝜇=0.3이 미끄러짐, 𝜇=0.3이 구름 조건에 해당하므로, 본 예제는 경사를 고정한 채 마찰계수만 바꿔 두 운동 모드를 함께 검증한다.

Fig. 6과 Fig. 7은 질량중심의 경사 방향 변위 시간 이력을 이론해와 비교한 결과이다. 미끄러짐 조건(𝜇=0.3, Fig. 6)에서는 t = 0.6s에서 약 1.28m의 변위가 나타나 이론해와 상대 오차 1.53 %로 근접하였다. 구름 조건(𝜇=0.6, Fig. 7)에서는 회전 관성에 의한 감속으로 약 1.08m의 변위를 보여 이론해보다 5.9 % 큰 변위가 나타났다. 두 경우 모두 변형률이 회전 성분의 약 2 %이하에 머물러, 변형이 없는 강체에 가까운 회전을 유지하였다. 질량중심의 이동 속도와 표면의 회전 속도의 비는 미끄러짐 조건에서 3.02로 접촉점에서 미끄러지며 내려간 반면, 구름 조건에서는 1.14로 약 12 %의 미끄러짐이 남았다. 다만 미끄러짐은 식 (20)의 슬립 분기에 따라 변위를 오히려 줄이므로, +5.87 %의 오차는 평균 접선력이 구름 유지에 필요한 마찰력에 미달한 결과이다. 이는 본 조건(𝜃=60°, 𝜇=0.6)의 마찰력이 전이 임계(tan𝜃≈3𝜇) 부근이어서 접촉력 변동에 마찰 포화가 민감하기 때문으로 실제로 접선 강성 $k_{f}$ 를 10배 높여도 미끄러짐은 줄지 않은 반면, 임계에서 먼 𝜇=0.8에서는 미끄러짐이 3.6% 감소하였다. 마찰계수에 따라 미끄러짐과 구름으로 분기되는 거동이 재현되어 본 마찰 모델의 타당성을 확인하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-05/N0040390505/images/Figure_jcoseik_39_05_05_F6.jpg

Fig. 6.

Time history of the center-of-mass position along the incline for the cylinder with 𝜇=0.3 (slip case)

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-05/N0040390505/images/Figure_jcoseik_39_05_05_F7.jpg

Fig. 7.

Time history of the center-of-mass position along the incline for the cylinder with 𝜇=0.6 (rolling case)

3.3 고무 링의 대변형 충돌 해석

본 프레임워크의 강건성과 에너지 보존 특성을 평가하기 위해, Fig. 8과 같이 내반경 0.03m, 외반경 0.04m의 동일한 두 고무 링이 각각 $v_{0} = 0.06 c_{v}$ 의 균일 속도로 정면 충돌하는 대변형 문제를 해석하였다. 재료 물성과 수치 파라미터는 Wu 등(2025)을 참고하여 Table 3에 정리하였으며, 페널티 강성은 3.2절과 동일한 기준으로 선정하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-05/N0040390505/images/Figure_jcoseik_39_05_05_F8.jpg

Fig. 8.

Initial configuration of the colliding rubber rings problem

Table 3.

Material parameters of the colliding rubber rings

Property	Symbol	Value
Young’s modulus	$E$	10MPa
Density	𝜌₀	1200kg/m³
Poisson’s ratio	𝜈	0.4
Initial velocity of each ring	$v_{0}$	0.06 $c_{v}$
Friction coefficient	𝜇	0.0
Normal stiffness	$k_{r}$	1 × 10⁷N/m
Tangential stiffness	$k_{f}$	1 × 10⁵N/m

Fig. 9는 $t = 0, 2.4, 4.8, 7.2, 9.6, 12.0$ =0, 2.4, 4.8, 7.2, 9.6, 12.0ms의 여섯 시점에서 두 링의 변형 형상과 von Mises 응력 분포를, Fig. 10은 좌･우 링 각각의 총･운동･변형 에너지의 시간 이력을 도시하였다. 충돌 초기에 접촉 영역에서 큰 변형이 발생하지만, 입자 분포의 불규칙화 없이 안정적인 응력장을 유지하였다. Fig. 10의 곡선에 나타나듯이 최대 변형 시점에서 운동 에너지의 대부분이 변형 에너지로 전환된 뒤 두 링이 분리되며 다시 운동 에너지로 회복되었다. 두 링의 좌･우 에너지 곡선이 구분되지 않을 정도로 일치하였으며, 이는 표면 인식 접촉 모델이 작용-반작용에 따라 충돌 과정 전반에서 대칭성을 보존함을 뜻한다. 또한 총 에너지의 상대 변동이 초기 운동 에너지 대비 2 % 이내에 유지되어, 수치적 소산 없이 탄성 충돌의 에너지 보존이 재현되었다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-05/N0040390505/images/Figure_jcoseik_39_05_05_F9.jpg

Fig. 9.

Snapshots of the colliding rubber rings colored by von Mises stress at $t = 0, 2.4, 4.8, 7.2, 9.6, 12.0$ =0, 2.4, 4.8, 7.2, 9.6, 12.0ms

https://cdn.apub.kr/journalsite/sites/jcoseik/2026-039-05/N0040390505/images/Figure_jcoseik_39_05_05_F10.jpg

Fig. 10.

Time histories of total, kinetic, and strain energies of the left and right rings

4. 결 론

본 연구에서는 Total Lagrangian SPH 정식화를 적용한 접촉･대변형 수치 해석 프레임워크를 구축하였다. 기준 형상에 기반한 TL-SPH 정식화에 전단 보정을 결합하여 대변형 시 발생하는 아워글래스 모드를 억제하였으며, 접촉 처리에는 표면 인식 페널티 모델과 쿨롱 마찰 모델을 TL-SPH 운동 방정식에 결합하여 충돌과 마찰을 단일 프레임워크에서 다루었다.

세 가지 벤치마크를 통해 프레임워크를 검증하였다. 2D 진동판 문제에서는 관측점 변위의 시간 이력이 이론해와 근접하였고, 공간 해상도에 대한 수렴성과 장시간 모사에서의 진폭 보존을 확인하였다. 경사면 문제에서는 두 마찰 계수(𝜇=0.3, 𝜇=0.6)에 대한 질량중심 변위가 이론해와 부합하였고, 마찰 계수 크기에 따라 미끄러짐과 구름으로 분기되는 거동이 재현되었다. 고무 링 충돌 문제에서는 von Mises 응력장이 입자 분포의 불규칙화없이 유지되었으며, 좌･우 링의 에너지 이력이 대칭적으로 일치하며 충돌 전반에서 총 에너지가 보존되었다.

본 프레임워크는 동적 거동 해석, 접촉･마찰 처리, 대변형 시 에너지 보존을 단일 정식화에서 다룰 수 있음을 보였다. 다만 표면 인식 접촉 모델은 법선의 가중평균인 식 (9)에 기반하므로, 날카로운 모서리와 같이 법선이 불연속인 형상에서는 접촉력 추정이 불안정해 질 수 있어 별도의 보정이 필요하다. 또한 본 연구는 2차원 검증에 한정되므로, 향후 3차원으로 확장하여 곡면 접촉 및 대규모 접촉 문제에 적용하고자 한다.

Acknowledgements

본 연구는 기후에너지환경부(MCEE)와 한국에너지기술평가원(KETEP)의 지원을 받아 수행한 연구 과제입니다(No. RS-2026-25527655).

References

Adepu, D., Ramachandran, P. (2024) Improved Collision Handling of Elastic Solids in SPH Using a Contact Force Model, Comput. Part. Mech., 11(2), pp.545~563.

10.1007/s40571-023-00638-2

Belytschko, T., Guo, Y., Liu, W.K., Xiao, S.P. (2000) A Unified Stability Analysis of Meshless Particle Methods, Int. J. Numer. Methods Eng., 48(9), pp.1359~1400.

10.1002/1097-0207(20000730)48:9<1359::AID-NME829>3.0.CO;2-U

Bonet, J., Kulasegaram, S. (2002) A Simplified Approach to Enhance the Performance of Smooth Particle Hydrodynamics Methods, Appl. Math. Comput., 126(2-3), pp.133~155.

10.1016/S0096-3003(00)00143-0

Gray, J.P., Monaghan, J.J., Swift, R.P. (2001) SPH Elastic Dynamics, Comput. Methods Appl. Mech. Eng., 190(49-50), pp.6641~6662.

10.1016/S0045-7825(01)00254-7

Lee, S.W., Heo, J.Y., Kim, K.H. (2023) Development of Meshless Method Using Least-Squares Method with Geometric Conservation Law for Structural Dynamic Analysis, J. Comput. Struct. Eng. Inst. Korea, 36(1), pp.67~74.

10.7734/COSEIK.2023.36.1.67

Monaghan, J.J. (2005) Smoothed Particle Hydrodynamics, Rep. Prog. Phys., 68(8), pp.1703~1759.

10.1088/0034-4885/68/8/R01

Swegle, J.W., Hicks, D.L., Attaway, S.W. (1995) Smoothed Particle Hydrodynamics Stability Analysis, J. Comput. Phys., 116(1), pp.123~134.

10.1006/jcph.1995.1010

Wendland, H. (1995) Piecewise Polynomial, Positive Definite and Compactly Supported Radial Functions of Minimal Degree, Adv. Comput. Math., 4(1), pp.389~396.

10.1007/BF02123482

Wu, D., Tang, X., Zhang, S., Hu, X. (2025) Unified Non-Hourglass Formulation for Total Lagrangian SPH Solid Dynamics, Comput. Mech., 75(3), pp.1081~1113.

10.1007/s00466-024-02552-7

Journal of the Computational Structural Engineering Institute of Korea ISSN:1229-3059(Print) 2287-2302(Online) 한국전산구조공학회 논문집

Preview

Total Lagrangian SPH Numerical Framework for Large Deformation and Collision Simulations

ABSTRACT

MAIN

(1)

(2)

(3)

(4)

(5)

(6)

(7)

Fig. 1.

Schematic of the surface-aware contact model

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

Fig. 2.

Initial configuration of the 2D oscillating plate

Table 1.

Material properties of the 2D oscillating plate

(17)

(18)

(19)

Fig. 3.

Time histories of the vertical displacement at the observation point for three spatial resolutions (H/∆p=10,20,40=10, 20, 40), compared with the analytical solution

Fig. 4.

Time histories of total, kinetic, and strain energies (H/∆p=40=40)

Fig. 5.

Sliding elastic cylinder on an inclined plane (𝜃=60°)

Table 2.

Parameters used in the sliding-on-slope problem

(20)

Fig. 6.

Time history of the center-of-mass position along the incline for the cylinder with 𝜇=0.3 (slip case)

Fig. 7.

Time history of the center-of-mass position along the incline for the cylinder with 𝜇=0.6 (rolling case)

Fig. 8.

Initial configuration of the colliding rubber rings problem

Table 3.

Material parameters of the colliding rubber rings

Fig. 9.

Snapshots of the colliding rubber rings colored by von Mises stress at t=0,2.4,4.8,7.2,9.6,12.0=0, 2.4, 4.8, 7.2, 9.6, 12.0ms

Fig. 10.

Time histories of total, kinetic, and strain energies of the left and right rings

Acknowledgements

References

Time histories of the vertical displacement at the observation point for three spatial resolutions ( $H / ∆ p = 10, 20, 40$ =10, 20, 40), compared with the analytical solution

Time histories of total, kinetic, and strain energies ( $H / ∆ p = 40$ =40)

Snapshots of the colliding rubber rings colored by von Mises stress at $t = 0, 2.4, 4.8, 7.2, 9.6, 12.0$ =0, 2.4, 4.8, 7.2, 9.6, 12.0ms