이산요소법

시뮬레이션 계산물리 물리 마지막 수정: 2026-07-16 04:14:02

1. 개요[편집]

이산요소법 Discrete Element Method
약칭	DEM
제안	Cundall & Strack (1979)
기술 방식	라그랑주 (입자 추적)
핵심 모델	스프링-대시팟, 헤르츠 접촉
주 사용처	분체, 골재, 제약 공정, 광산

모래알 하나하나를 다 풀겠다는 발상. 문제는 모래알이 아보가드로 수만큼 있다는 것.

이산요소법(Discrete Element Method, DEM)은 입자 하나하나를 독립된 강체로 보고, 입자 사이의 접촉력을 계산해 뉴턴의 운동 방정식을 시간 적분하는 라그랑주 기법이다. 연속체 가정이 통하지 않는 분체(粉體, granular material) — 모래, 곡물, 알약, 광석, 토너 가루 — 를 다루기 위해 태어났다.

컨티뉴엄으로 뭉개면 안 되는 이유는 명확하다. 모래는 쌓으면 고체처럼 안식각을 유지하고, 기울이면 액체처럼 흐르고, 흔들면 기체처럼 튄다. 이 삼중인격을 하나의 구성방정식으로 기술하려는 시도는 수십 년째 진행 중이지만 아직 완승하지 못했다.¹ 그래서 그냥 입자를 다 세는 무식한 길을 택한 것이 DEM이다.

시초는 1979년 Cundall과 Strack의 논문 A discrete numerical model for granular assemblies. 원래는 암반역학에서 절리(joint)로 쪼개진 블록의 거동을 보려는 목적이었다.

2. 연접촉 모델[편집]

DEM의 핵심 아이디어는 연접촉(soft-sphere) 가정이다. 실제 입자는 접촉 시 국소적으로 변형되지만, DEM은 입자를 강체로 두는 대신 겹치는 것을 허용하고 그 겹침량 $\delta$ 를 변형의 대리 변수로 쓴다. 겹침에 비례하는 반발력을 스프링으로, 에너지 소산을 대시팟으로 모사한다.

가장 널리 쓰이는 선형 스프링-대시팟(LSD) 모델의 법선 방향 접촉력은 다음과 같다.

F_n = -k_n \delta_n - \gamma_n \dot{\delta}_n

여기서 $k_n$ 은 법선 강성, $\gamma_n$ 은 감쇠 계수다. 접선 방향은 스프링-대시팟에 쿨롱 마찰 한계를 씌운다.

F_t = \min\left( \left| -k_t \delta_t - \gamma_t \dot{\delta}_t \right|,\; \mu |F_n| \right)

즉 접선력은 미끄러지기 전까지는 탄성적으로 버티다가, $\mu |F_n|$ 을 넘으면 미끄러진다. 이 한 줄이 안식각과 아칭(arching)을 만들어내는 원흉이자 영웅이다.

2.1. 헤르츠 접촉[편집]

선형 스프링은 계산이 싸지만 물리적으로는 거짓말이다. 실제 두 구가 눌릴 때 접촉 면적이 겹침량에 따라 커지므로, 힘은 비선형으로 증가한다. 헤르츠(Hertz) 이론은 이를 $\delta^{3/2}$ 의존성으로 준다.

F_n = \frac{4}{3} E^{*} \sqrt{R^{*}}\, \delta_n^{3/2}

$E^{*}$ 는 등가 영률, $R^{*}$ 는 등가 반경( $1/R^{*} = 1/R_1 + 1/R_2$ )이다. 접선 방향까지 엄밀히 다루면 Hertz–Mindlin 모델이 되고, 대부분의 상용 DEM 코드는 이것을 기본값으로 깔아둔다. 정확도는 좋아지지만 강성이 겹침에 따라 변하므로 시간 스텝 산정이 더 까다로워진다.

응집성 분체(미분말, 습한 모래)에는 JKR이나 액체 브릿지 모델 같은 부착력 항을 추가로 얹는다. 알약 코팅이나 토너처럼 반데르발스 힘이 중력을 이기는 스케일에서는 이걸 빼면 결과가 통째로 틀린다.

3. 시간 적분과 레일리 시간[편집]

DEM은 명시적 시간 적분(대개 베를레 적분이나 그 변형)을 쓴다. 즉 시간 스텝이 안정성 한계에 묶인다. 기준이 되는 것이 레일리 시간(Rayleigh time) — 입자 표면을 따라 전파하는 레일리 파가 입자를 한 번 가로지르는 시간이다.

t_R = \frac{\pi R \sqrt{\rho / G}}{0.1631\nu + 0.8766}

실무에서는 $\Delta t$ 를 $t_R$ 의 10~20% 수준으로 잡는 것이 국룰이다. 문제는 이게 살인적으로 작다는 것. 반경 1 mm 강철 구슬이면 $t_R$ 이 마이크로초 단위라, 1초짜리 공정을 보려면 수십만 스텝을 돌려야 한다.

그래서 현장에서는 온갖 편법이 동원된다.

강성 낮추기: 영률을 실제 값의 1/100~1/1000로 깎아 $\Delta t$ 를 키운다. 벌크 거동이 강성에 둔감한 준정적 문제에서는 통하지만, 충돌 반발이 중요한 문제에서는 배신당한다.
입자 조대화(coarse-graining): 입자를 실제보다 크게 잡고 개수를 줄인다. 스케일링 법칙이 성립하는 범위 안에서만 유효하다.
GPU: GPU 컴퓨팅으로 수천만 입자를 돌리는 것이 이제는 특별한 일이 아니다.

4. 이웃 탐색[편집]

입자 $N$ 개의 모든 쌍을 검사하면 $O(N^2)$ . 천만 입자면 그냥 우주가 식을 때까지 기다려야 한다. 그래서 공간을 격자 셀로 나누고 인접 셀만 뒤지는 링크드 셀(linked-cell) 방식이나, 컷오프보다 약간 큰 반경으로 후보 목록을 만들어 몇 스텝간 재사용하는 버렛 리스트(Verlet list)를 쓴다. 둘 다 $O(N)$ 수준으로 떨어뜨린다. 이 부분은 분자동역학에서 그대로 빌려온 자산이고, 더 일반적인 논의는 공간 분할 자료구조와 충돌 감지를 보면 된다.

비구형 입자(쌀알, 자갈, 정제)를 다룰 때는 구 여러 개를 붙여 만드는 클럼프(clumped sphere), 다면체 직접 접촉, 초타원체(superquadric) 방식이 쓰인다. 다면체 접촉 판정에는 GJK 알고리즘이 동원되기도 한다.

5. CFD-DEM 연성[편집]

입자가 공기나 물 속을 날아다니면 유체를 무시할 수 없다. CFD-DEM은 유체를 유한체적법으로 오일러 격자에서 풀고, 입자는 DEM으로 라그랑주 추적하면서 둘을 항력·부력·압력 구배 힘으로 주고받는 연성 기법이다.

비해상(unresolved) CFD-DEM: 격자 셀이 입자보다 크다. 입자의 존재는 셀의 다공도와 항력 소스항으로만 반영된다. 유동층 반응기, 공압 수송(pneumatic conveying), 사이클론 분리기의 표준. 다공도 개념 자체는 다공성 매질 유동과 뿌리가 같다.
해상(resolved) 방식: 격자가 입자보다 작아 입자 표면 유동을 직접 푼다. 침지경계법이나 격자 볼츠만 방법과 조합하는 경우가 많다. 정확하지만 입자 수백~수천 개가 한계.

항력 상관식은 Ergun, Wen–Yu, Di Felice, Gidaspow 등이 쓰이는데, 어느 걸 고르느냐로 유동층 팽창 높이가 수십 % 달라진다. “상관식이 다 했다”는 말이 나오는 지점.

6. 비슷해 보이는 것들과의 차이[편집]

입자를 추적한다는 점에서 SPH나 분자동역학과 헷갈리기 쉬운데, 물리적 위상이 전혀 다르다.

기법	입자의 정체	상호작용	목적
DEM	실제 알갱이	접촉·마찰	분체 거동
SPH	연속체의 보간점	커널 가중 평균	연속체 유동
분자동역학	원자·분자	퍼텐셜 기울기	열역학적 물성

핵심 차이는 이것이다. SPH 입자는 물리적 실체가 없다 — 그냥 나비에-스토크스를 이산화한 계산점일 뿐이라 개수를 늘리면 정확도가 올라간다. DEM 입자는 실제 알갱이라 개수가 물리적으로 정해져 있어, 마음대로 늘리면 다른 물질이 된다. 격자 세분화라는 개념이 DEM에는 없는 이유다.² MD와의 차이는 스케일과 소산 — MD 힘은 보존적 퍼텐셜에서 나오지만, DEM 접촉력은 마찰과 비탄성 충돌로 에너지를 끊임없이 갉아먹는 비보존계다.

7. 활용 분야[편집]

제약: 정제 코팅, 블렌딩 균일도, 캡슐 충전. FDA가 QbD(Quality by Design)를 밀면서 DEM 수요가 폭발했다.
광산·골재: 분쇄기(SAG mill), 컨베이어 이송, 슈트 마모 예측. 라이너 마모를 DEM으로 예측해 교체 주기를 잡는 것이 이 바닥 표준.
농업: 콤바인 내부 곡물 거동, 종자 파종기.
토목·지반: 사면 붕괴, 토석류, 말뚝 관입.
에너지: 유동층 보일러, 화학적 순환 연소 — 여기서는 연소 시뮬레이션과도 손을 잡는다.

대표 코드로는 상용 Rocky DEM, Altair EDEM, 오픈소스 LIGGGHTS(LAMMPS에서 갈라져 나왔다), Yade, 그리고 ANSYS Fluent의 내장 DEM 등이 있다.

8. 현업에서의 현실[편집]

파라미터가 너무 많다. 반발계수, 정지 마찰, 구름 마찰, 영률… 이걸 다 실측할 수는 없으니 결국 안식각 실험 하나 놓고 역으로 튜닝한다. 이걸 점잖게 “캘리브레이션”이라 부른다.³
구름 마찰(rolling friction)은 물리적으로 가장 의심스러우면서 결과에는 가장 크게 영향을 미치는 파라미터다. 완벽한 구는 굴러가면 멈추지 않기 때문에, 비구형 입자의 형상 저항을 이 한 숫자에 몰아넣는다.
입자 수 100만은 이제 노트북에서도 돈다. 문제는 실제 사일로 안에 곡물이 100억 개 들어 있다는 것.
“실험이랑 안 맞아요”의 원인 1순위는 늘 입자 형상이다. 완벽한 구는 자연에 존재하지 않는다.⁴

9. 관련 문서[편집]

분체를 연속체로 기술하려는 시도 중 가장 성공적인 것이 $\mu(I)$ 유변학 모델인데, 이것도 준정적 영역과 충돌 지배 영역 사이의 중간 지대에서 잘 작동한다. 물리학자들이 “모래는 왜 이러는가”로 논문을 쓰는 동안 공학자들은 그냥 DEM을 돌린다. ↩
이걸 헷갈려서 “DEM 격자 수렴성 검토를 하라”는 리뷰 코멘트를 받아본 사람이 의외로 많다. 정답은 “입자 크기는 물성이지 이산화 파라미터가 아닙니다”인데, 코스-그레이닝을 쓰는 순간 이 대답도 반쯤 거짓말이 된다. ↩
검증 및 확인 관점에서 보면 이건 확인(validation)이 아니라 캘리브레이션이다. 캘리브레이션에 쓴 실험으로 다시 검증하면 안 되는데, 실험이 하나뿐이라 그냥 한다. 모두가 알지만 아무도 말하지 않는다. ↩
그래서 구 여러 개를 붙여 쌀알 모양을 만드는데, 접촉 판정 비용이 구 개수에 비례해서 늘어난다. 형상을 정확히 할수록 계산이 느려지는 정직한 세계. ↩