Метод дискретних елементів (DEM): симуляція гранульованих матеріалів

3D Simulations

⚙️ Фізика · Чисельні методи

📅 Червень 2026⏱ 16 хв🔴 Просунутий рівень · Останнє оновлення: 3 липня 2026 р.

Пісок, зерно в силосі, каміння в зсуві, фармацевтичні порошки — ці гранульовані матеріали поводяться ні як тверде тіло, ні як рідина. Вони утворюють лавини, будують стійкі купи з певним кутом природного укосу, передають силу вздовж невидимих ланцюгів контакту й непередбачувано заклинюються. Метод дискретних елементів моделює кожне зерно окремо, описуючи всі ці явища з перших принципів.

Автор: Команда MySimulator · Редакційна перевірка: Редакція MySimulator

1. Чому гранульовані матеріали складні

Купа піску має ~10¹⁰ зерен. Навіть малий лабораторний зразок із 10 000 зерен не можна описати класичною механікою суцільного середовища (яка трактує матеріали як диференційовні поля), бо:

Контакти однобічні: зерна штовхають, але ніколи не тягнуть. Сили лише стискальні, що створює асиметричну, залежну від історії поведінку.
Тертя напрямлене й не залежить від швидкості: кулонівська модель тертя повʼязує дотичну силу з нормальною, а не зі швидкістю.
Немає термодинамічної рівноваги: гранульовані системи атермічні — випадкова теплова енергія за кімнатної температури мізерна порівняно з вагою зерна. Система завжди перебуває в нерівноважному метастабільному стані.
Заклинені стани: струшені зерна течуть; зупинені зерна жорстко заклинюються. Перехід заклинення погано описується стандартною теорією фазових переходів.

Метод дискретних елементів (DEM), запроваджений Кандаллом і Страком у 1979 році, трактує кожне зерно як тверде тіло й обчислює всі контактні сили явно. Колективна поведінка зерен виникає з мільйонів окремих двотільних взаємодій — жодного конститутивного закону для об'ємного матеріалу не припускається.

DEM проти SPH проти FEM для гранульованих течій: DEM точний, але затратний (кількість частинок обмежена ~10⁶–10⁷). Згладжена гідродинаміка частинок (SPH) може опрацьовувати 10⁸ частинок, але використовує рідинні припущення. Континуальний FEM з реологією Друкера–Прагера чи μ(I) — найшвидший, але пропускає ефекти масштабу зерна, як-от силові ланцюги та сегрегацію за розміром. Вибір залежить від масштабу задачі та потрібної деталізації.

2. Виявлення контактів

Перед обчисленням сил ми маємо визначити, які пари частинок перебувають у контакті. Для сфер (переважна форма зерна в DEM) контакт простий:

// Виявлення контакту «сфера-сфера» Частинки i та j мають положення xᵢ, xⱼ і радіуси rᵢ, rⱼ Відстань: d = |xᵢ - xⱼ| Перекриття: δ = rᵢ + rⱼ - d Контакт існує тоді й лише тоді, коли δ > 0 (зерна перекриваються) Одиничний вектор нормалі контакту: n = (xⱼ - xᵢ) / d (вказує від i до j)

Несферичні зерна (еліпсоїди, многогранне каміння, заокруглені куби) потребують складнішого виявлення контактів — алгоритмів GJK або SAT — і у 10–100× затратніші на контактну пару, тому дослідження DEM переважно використовують сфери.

Широка фаза: пошук кандидатних пар

Перевірка кожної пари з N частинок коштує O(N²). За N = 100 000 це 10¹⁰ перевірок на крок часу — неможливо. Просторова хеш-сітка або ієрархія обмежувальних обʼємів (BVH) зводить це до O(N) з малою константою:

// Просторовий хеш на рівномірній сітці (для однакових розмірів зерен) cell_size = 2 \times max_radius Для кожної частинки i: cell = floor(xᵢ / cell_size) // 3D цілочисловий індекс комірки hash(cell) \to bucket // хеш у 1D-масив Пошук сусідів: перевірити частинку проти всіх частинок у 3\times3\times3 = 27 сусідніх комірках. Середня кількість частинок на комірку \approx частка заповнення \times 27 \approx 15 Кандидати на контакт скорочуються з O(N²) \to O(15N) на кадр

Сітку потрібно перебудовувати щокроку (частинки рухаються). Для дуже динамічних симуляцій структуру «звʼязаний список на комірку», що підтримує вставлення/вилучення за O(1), віддають перевагу перед перебудовою з нуля.

3. Моделі контактних сил

Лінійна модель «пружина-демпфер» (LSD)

Найпростіша контактна модель подає кожен контакт як пружину (пружне відштовхування, пропорційне перекриттю) і демпфер (вʼязке згасання, пропорційне відносній швидкості):

// Нормальна контактна сила (лише стиск, без розтягу) Перекриття: δₙ = rᵢ + rⱼ - |xᵢ - xⱼ| (додатне в контакті) Відносна швидкість уздовж нормалі: vₙ = (vᵢ - vⱼ) \cdot n Нормальна сила: Fₙ = kₙ \cdot δₙ - γₙ \cdot vₙ (пружина - демпфер) kₙ = нормальна жорсткість пружини [Н/м] γₙ = нормальний коефіцієнт згасання [Н\cdotс/м] n = одиничний вектор нормалі контакту Застосувати до частинки i: force += Fₙ \cdot n Застосувати до частинки j: force -= Fₙ \cdot n (3-й закон Ньютона)

Жорсткість пружини kₙ — критичний параметр. Вона має бути достатньо великою, щоб перекриття залишалися фізично малими (δ ≪ r), але більше kₙ потребує менших кроків часу (умова стійкості: Δt ≤ π × √(m/kₙ)). На практиці kₙ обирають так, щоб отримати перекриття ~0,5–2% радіуса частинки.

Контактна модель Герца–Міндліна

Фізично точніша модель виводиться з контактної теорії Герца, яка передбачає, що площа контакту зростає з навантаженням, а сила масштабується як δ^(3/2), а не лінійно за δ:

// Модель нормальної сили Герца (нелінійна) Ефективний радіус: R* = (rᵢ \times rⱼ) / (rᵢ + rⱼ) Ефективний модуль: E* = ((1-νᵢ²)/Eᵢ + (1-νⱼ²)/Eⱼ)⁻¹ Жорсткість (залежна від навантаження): kₙ_hertz = (4/3) \times E* \times \sqrtR* \times \sqrtδₙ (змінюється з δₙ) Сила Герца: Fₙ = kₙ_hertz \times δₙ = (4/3) \times E* \times \sqrtR* \times δₙ^(3/2) Нелінійне співвідношення сила-перекриття відображає справжню фізику пружної деформації на поверхні контакту зерна.

Герц–Міндлін потребує параметрів матеріалу (модуль Юнга E, коефіцієнт Пуассона ν), які можна виміряти експериментально, що робить калібрування під реальні матеріали простим. LSD потребує вибору kₙ дещо довільно.

4. Інтегрування за часом

DEM використовує другий закон Ньютона як для поступального, так і для обертального руху кожної частинки:

// Рівняння руху для частинки i Поступальне: mᵢ \cdot ẍᵢ = Σ Fᵢ_contact + mᵢ \cdot g Обертальне: Iᵢ \cdot ω̇ᵢ = Σ Tᵢ_contact mᵢ = маса (сфера: (4/3)π rᵢ³ ρ) Iᵢ = момент інерції (сфера: (2/5) mᵢ rᵢ²) g = прискорення вільного падіння Tᵢ = момент від дотичних контактних сил Інтегрування (явний Ейлер або швидкісний Верле): vᵢ(t+Δt) = vᵢ(t) + (Fᵢ/mᵢ) \cdot Δt xᵢ(t+Δt) = xᵢ(t) + vᵢ(t+Δt) \cdot Δt (симплектичний Ейлер) ωᵢ(t+Δt) = ωᵢ(t) + (Tᵢ/Iᵢ) \cdot Δt

Кроки часу в DEM крихітні: зазвичай Δt ≈ 10⁻⁶–10⁻⁵ секунди (набагато менші за тривалість симуляції від секунд до хвилин). Критичний крок часу задається власним періодом коливань системи «пружина-демпфер». Його порушення спричиняє чисельну нестійкість — перекриття зростають необмежено.

Критичний крок часу: Δt_crit = π × √(m_min / kₙ_max). Для зерен піску діаметром 1 мм (m ≈ 10⁻⁶ кг) з kₙ = 10⁵ Н/м, Δt_crit ≈ 10⁻⁵ с. Симуляція 1 секунди фізичного часу потребує ~100 000 кроків часу. За 10 000 частинок із 5 контактами кожна це 5×10⁹ обчислень сил на секунду часу симуляції.

5. Тертя, кочення та когезія

Кулонівське дотичне тертя

Зерно чинить опір ковзанню в дотичному напрямку. Класична кулонівська модель тертя обмежує дотичну силу часткою нормальної сили:

// Дотична (зсувна) сила з кулонівським обмеженням тертя Відстежувати дотичне зміщення δₜ: інтегрувати відносну дотичну швидкість δₜ(t+Δt) = δₜ(t) + vₜ_rel \cdot Δt Дотична пружна сила: Fₜ_spring = kₜ \cdot δₜ (kₜ \approx (2/7) kₙ для Герца-Міндліна) Кулонівська межа: if |Fₜ_spring| > μ \times Fₙ: Fₜ = μ \times Fₙ \times sign(Fₜ_spring) // тертя ковзання скинути δₜ = Fₜ / kₜ // обрізати пружину, щоб уникнути накопичення else: Fₜ = Fₜ_spring // статичне тертя μ = коефіцієнт тертя ковзання (зазвичай 0,3-0,7)

Опір коченню

Ідеальні сфери котяться без опору коченню — вони утворювали б нереалістично пласкі купи. Момент тертя кочення протидіє обертанню в контакті, параметризований коефіцієнтом тертя кочення μᵣ ≪ μ:

// Момент опору коченню (моделі, обмежені контактом) Величина моменту: |Tᵣ| \leq μᵣ \times rᵢ \times Fₙ Напрям: протидіє ω_relative у контакті (відносна кутова швидкість, спроєктована на площину контакту)

Когезія (Ван дер Ваальса / капілярна)

Дрібні порошки (d < 100 мкм) зазнають значного міжчастинкового притягання від сил Ван дер Ваальса. Вологий пісок набуває когезії від рідинних містків (капілярна адгезія). Це моделюють, додаючи притягальну складову сили, за допомогою розширення контактної механіки Джонсона–Кендалла–Робертса (JKR) чи Дерягіна–Мюллера–Топорова (DMT), або спрощеного підходу «густини енергії когезії», де зерна трохи поза контактом усе ще притягуються.

6. Емерджентні явища: силові ланцюги та кут природного укосу

Силові ланцюги

Вдарте по купі піску — і сила удару не розходиться рівномірно: вона мандрує вздовж силових ланцюгів: вузьких розгалужених шляхів сильно навантажених зерен, оточених майже ненавантаженими «глядацькими» зернами. Неоднорідність вражає: найбільш навантажені зерна несуть у 10× більше за середнє навантаження.

Симуляції DEM передбачили силові ланцюги до того, як їх зобразили експериментально. Експерименти з фотопружними дисками (із застосуванням двозаломних дисків, що оптично показують напруження в поляризованому світлі) підтвердили передбачення DEM у 1990-х роках. Силові ланцюги критично важливі для розуміння передачі напружень у гранульованих середовищах і раптового початку течії (лавини).

Кут природного укосу

Насипте пісок у купу: вона зростає до максимального кута схилу — кута природного укосу θᵣ — і зупиняється. Додайте ще — і малі лавини підтримують цей схил. DEM природно це відтворює:

Коефіцієнт кулонівського тертя μ → tan(θᵣ) ≈ μ (для однакових сфер на пласкій поверхні)
Тертя кочення збільшує θᵣ (кутасті частинки мають вище ефективне тертя кочення)
Розподіли розмірів зерен впливають на θᵣ (ширші розподіли → вище θᵣ через взаємне зчеплення)

Емпіричне співвідношення: θᵣ \approx arctan(μ_eff) μ_eff = ефективне тертя з урахуванням форми зерна, розподілу розмірів і опору коченню Для гладких скляних сфер: θᵣ \approx 22°-25° (μ \approx 0,4-0,45) Для сухого піску: θᵣ \approx 30°-35° (включає зчеплення за формою) Для вологого піску: θᵣ \approx 45°+ (капілярна когезія) Для кутастого гравію: θᵣ \approx 40°-45°

Сегрегація за розміром (ефект бразильського горіха)

Струсіть контейнер зі змішаними розмірами зерен: більші зерна спливають угору — ефект бразильського горіха. Симуляції DEM відтворюють це через механізм конвекції/заповнення порожнин: малі зерна падають у проміжки під великими зернами під час вібрації, поступово підіймаючи великі зерна вгору. Ефект змінюється на протилежний (зворотний ефект бразильського горіха) для щільних великих зерен і чутливий до співвідношення розмірів зерен, співвідношення густин та параметрів вібрації.

7. Продуктивність: просторове хешування та DEM на GPU

Продуктивність DEM на CPU

Добре оптимізований код DEM на CPU опрацьовує ~10⁵–10⁶ частинок на інтерактивних швидкостях (кілька кроків часу на секунду для повільної динаміки). Ключові оптимізації:

Списки сусідів Верле: перебудовувати кандидатів на контакт лише раз на кілька кроків. Щокроку перевіряти, чи перемістилися частинки більш ніж на «товщину оболонки» від останньої перебудови.
Звʼязаний список на комірку: вставлення/вилучення з просторової сітки за O(1) під час руху частинок.
SIMD: векторизувати обчислення сил по пакетах контактних пар за допомогою AVX2/AVX-512.
OpenMP: розпаралелити обчислення сил по частинках (стережіться станів гонитви на спільних масивах сил — використовуйте потокові локальні акумулятори або атомарні операції).

DEM на GPU

DEM на GPU сягає 10⁷–10⁸ частинок, виконуючи ядра контактних сил з одним GPU-потоком на контактну пару. Виклики, властиві DEM на GPU:

Атомарне накопичення сил: кілька потоків одночасно додають сили до тієї самої частинки. Використовуйте CUDA atomicAdd або розбийте на фази, де кожен потік «володіє» частинкою.
Просторове хешування на GPU: сортуйте частинки за ключем комірки за допомогою порозрядного сортування на GPU (бібліотека CUB). Потім скануйте, щоб знайти межі комірок. Перебудовуйте щокроку, але кожен крок — лише ~1 мс на сучасному GPU.
Розкладка памʼяті: AoS (масив структур) неефективний для кешу GPU. SoA (структура масивів) пакує x-координати разом, y-координати разом — значно краще коалесування.

// Структура CUDA-ядра для контактної сили (псевдокод) __global__ void computeForces(float* pos, float* vel, float* force, int* cellStart, int* cellSize, ...) { int i = blockIdx.x * blockDim.x + threadIdx.x; // індекс частинки float3 fi = {0,0,0}; // перебрати сусідні комірки (27 комірок) for each cell in neighbors(cell_of(i)): for j in particles_in(cell): if (i == j) continue; float3 delta = pos[j] - pos[i]; float dist = length(delta); float overlap = radii[i] + radii[j] - dist; if (overlap > 0): fi += contact_force(delta, dist, overlap, vel[i], vel[j]); atomicAdd(&force[i*3], fi.x); // записати результат atomicAdd(&force[i*3+1], fi.y); atomicAdd(&force[i*3+2], fi.z); }

8. Промислові застосування

DEM застосовують усюди, де гранульовані матеріали відіграють ключову роль:

Гірнича справа та переробка мінералів: подрібнення в кульовому млині (руда + сталеві кулі), проєктування конвеєрів, течія в силосах і швидкість вивантаження, режими лійкової течії проти масової. DEM передбачає закупорювання — коли арка зерен блокує бункер, — що критично для запобігання дорогим зупинкам виробництва.
Фармацевтика: покриття таблеток в обертових барабанах, однорідність змішування порошків, наповнення капсул. Моделі руйнування гранул передбачають різницю у швидкості розчинення. FDA дедалі частіше вимагає валідації процесів на основі DEM.
Геотехнічна інженерія: прогнозування дальності зсувів, селеві потоки, спричинене землетрусом розрідження ґрунту, проєктування підпірних стін, механіка підземних обвалів.
Харчова переробка: поводження з крупами, рисом, кавовими зернами — прогнозування сегрегації за розміром і густиною, зносу конвеєрів, оптимізація обʼємного зберігання.
Адитивне виробництво: сплавлення в порошковому шарі металу (SLM, EBM) — DEM передбачає густину й однорідність порошкового шару, що визначає звʼязок лазерної енергії та остаточну пористість деталі.
Планетологія: механіка реголіту на астероїдах (місія OSIRIS-REx використовувала моделі DEM для планування дотику та збору зразків на Бенну), утворення кратерів, динаміка кілець у кільцях Сатурна.

Комерційні коди: EDEM (Altair), Rocky DEM (ESSS), PFC (Itasca), LIGGGHTS (з відкритим кодом на основі LAMMPS). Відкриті MercuryDPM і YADE широко використовуються в дослідженнях. Усі реалізують модель Герца–Міндліна й підтримують прискорення на GPU. Мінімальний DEM з нуля на Python (NumPy + Numba) може опрацьовувати ~10 000 частинок у реальному часі для навчальних цілей.