Высокий
уровень распространения травмы опорно-двигательного аппарата, достигающий 24% в
структуре временной нетрудоспособности, и высокая встречаемость переломов
бедренной кости указывает на актуальность изучения механизма их образования
[1]. Тем не менее, большинство научных работ описывает механизм образования
переломов при одном виде деформации. Полученные в таких работах критерии
установления механизма образования повреждения используются в
судебно-медицинских экспертизах методом нестрогой аналогии. В реальной жизни
такие простые условия нагружения бедренной кости встречаются редко, поэтому
морфология реального перелома имеет множественные дополнительные признаки, не
характерные для одного вида деформации, что вызывает обоснованное сомнение у
эксперта в его правоте установления механизма травмы и требует разработки новых
методов установления механизма травмы.
Проведенные
ранее нами эксперименты по математическому моделированию методом
конечно-элементного анализа, широко применяемого для решения задач механики
деформируемого твёрдого тела в науке «сопротивление материалов», переломов
бедренной кости позволили валидировать твердотельную параметрическую модель
бедра, созданную по данным компьютерной томографии реальных людей, и применить
её для изучения механизма образования переломов бедренной кости в различных
условиях и обстоятельствах [2].
Целью
настоящей работы является установление
возможности
применения метода конечно-элементного анализа для исследования
сложнонапряженных состояний при переломе бедренной кости с последующей
визуализацией результатов.
В
ходе предварительной части исследования были проанализированы литературные
источники, в которых рассматривалась морфология переломов диафизов длинных
трубчатых костей, вызванных ударным воздействием тупых предметов. Данный анализ
основывался на данных, полученных при проведении экспериментов на биоманекенах
или в ходе практических наблюдений (342 эксперимента на 200 биоманекенах, 56
экспериментов на образцах костей нижних конечностей, 116 экспертных наблюдений)
[3]. Условия нагружения бедренной кости и характер разрушения, использованные в
модели, были взяты из этих работ для проверки достоверности модели. В
экспериментах ударное воздействие на биоманекены осуществлялось при
горизонтальном положении тела на металлической поверхности и ударе тупым
твердым предметом по передней поверхности диафиза длинной трубчатой кости, при
этом фиксировался дистальный отдел бедренной кости и осуществлялось давление на
таз с целью имитации нагрузки вращения верхнего отдела бедра по часовой или
против часовой стрелки. Сила удара, вызывающая появление переломов в
экспериментах, составляла от 1800Н до 1900Н [3].
Для
проверки возможности математического моделирования процесса формирования
повреждения перелома бедренной кости был использован метод конечных элементов
(МКЭ), реализованный с помощью программной среды ANSYS LS-DYNA, являющейся
популярной программой для конечно-элементного анализа, разработанная компанией
Ansys inc., которая используется для решения линейных и нелинейных динамических
задач механики деформируемых твердых тел, включая анализ разрушения (ANSYS inc.,
https://www.ansys.com) [4]. В
данной работе с помощью ANSYS LS-DYNA были воспроизведены эксперименты,
проведенные на биоманекенах и в ходе практических наблюдений при ударных
воздействиях. Для воссоздания условий нагружения модели бедра было выполнено
закрепление в области суставных поверхностей с жестким основанием и образование
упругоподатливой подложки, моделирующей основание Винклера. Сетка конечных
элементов была сгенерирована автоматически с использованием элементов типа
Solid размером 5 мм для моделирования объемного напряженно-деформированного
состояния.
Задано
ограничение движение кости по оси
Z
в проксимальной (верхней) части кости (область тазобедренного сустава) и по в
дистальной (нижней) части кости (область коленного сустава). Преднагрузка
смоделирована приложением силы в 100 Н к поверхности головки бедренной кости и
направленной в первом эксперименте по оси
Y,
а во втором – против оси
Y.
Приложение силы было под углом 90 градусов к оси кости модели бедра в проекции
средней трети передней поверхности бедра и реализовано смоделированным
цилиндрическим стальным индентором с радиусом кривизны 30 мм и длиной 180 мм.
При этом его движения ограничены во всех направлениях, кроме вертикального.
Ударное воздействие моделировалось скоростью воздействия 18 м в сек. Между
частями модели заданы условия воздействия:
-
между
губчатым веществом, компактным веществом и мышечной тканью – условия
неразрывного соединения, что достигается построением конформной сетки конечных
элементов на этих частях;
-
между
мышечной тканью, пуансоном и основанием – контакт с трением.
Учитывая,
что при проведении экспериментов на биоманекенах и в практических наблюдениях
свойства костной и мягких тканей не исследовались, были применены усредненные
механические свойства тканей из разных литературных источников (см. таблицы 1 и
2) [5-8]. Такие характеристики тканей хорошо описывали поведение материалов в
предыдущих экспериментах [2]. Свойства составляющих костной ткани описаны
моделями материала – изотропного упругопластического для губчатого вещества и
анизотропного упругого для компактного вещества. В качестве модели разрушения
принята теория максимального главного напряжения (теория Треска).
Таблица 1.Механические свойства компактного вещества костной ткани
Density
(плотность)
|
2000 kg/m3
(кг/м3)
|
Young’s modulus X
direction (модуль
Юнга
по
оси
X)
|
12 GPa
(ГПа)
|
Young’s modulus Y
direction (модуль
Юнга
по
оси
Y)
|
12 GPa
(ГПа)
|
Young’s
modulus
Z
direction
(модуль Юнга по оси
Z)
|
20 GPa
(ГПа)
|
Poisson’s
ratio
XY
(коэффициент Пуассона по оси
XY)
|
0.38
|
Poisson’s
ratio
YZ
(коэффициент Пуассона по оси
YZ)
|
0.22
|
Poisson’s
ratio
XZ
(коэффициент Пуассона по оси
XZ)
|
0.22
|
Shear
modulus
XY
(модуль сдвига по оси
XY)
|
4.5 GPa
(ГПа)
|
Shear
modulus
YZ
(модуль сдвига по оси
YZ)
|
5.6 GPa
(ГПа)
|
Shear
modulus
XZ
(модуль сдвига по оси
XZ)
|
5.6 GPa
(ГПа)
|
Compressive ultimate
strength (предел прочности при
сжатии)
|
0.205 GPa
(ГПа)
|
Tensile
ultimate
strength
(предел прочности при растяжении)
|
0.133 GPa
(ГПа)
|
Maximum
tensile
stress
(максимальное растягивающее напряжение)
|
52 MPa
(МПа)
|
Maximum
shear
stress
(максимальное сдвиговое напряжение)
|
65 MPa
(МПа)
|
Таблица 2.
Механические свойства губчатого вещества костной ткани
Density
(плотность)
|
127
kg/m3
(кг/м3)
|
Young’s modulus
(модуль
Юнга)
|
0.38
MPa
(МПа)
|
Poisson’s ratio
(коэффициент
Пуассона)
|
0.33
|
Bulk
modulus
(объёмный модуль
упругости)
|
0.37255 MPa
(МПа)
|
Shear modulus
(модуль
сдвига)
|
0.14286 MPa
(МПа)
|
Compressive ultimate
strength
(предел прочности при
сжатии)
|
6.23 MPa
(МПа)
|
Tensile
ultimate
strength
(предел прочности при растяжении)
|
8.4 MPa
(МПа)
|
Maximum
tensile
stress
(максимальное растягивающее напряжение)
|
8.4 MPa
(МПа)
|
Maximum
shear
stress
(максимальное сдвиговое напряжение)
|
7.4 MPa
(МПа)
|
Yield strength
(предел текучести)
|
1.75 MPa
(МПа)
|
Tangent modulus
(касательный модуль)
|
41.8 MPa
(МПа)
|
В
связи с тем, что в нашей модели разрушение мягких тканей не изучалось, а физические
свойства мягких тканей обобщенно использовались только как прокладка между
индентором и основанием, механические свойства мягких тканей взяты из свойств
баллистического геля. В качестве модели разрушения принята модель
Джонсона-Кука, учитывающая изменение критериев прочности от скорости нагружения
и температуры (см. таблицу 3):
Таблица 3.
Механические свойства мягких тканей
Density
(плотность)
|
1030
kg/m3
(кг/м3)
|
Young’s modulus
(модуль
Юнга)
|
59.862
kPa
(кПа)
|
Poisson’s ratio
(коэффициент
Пуассона)
|
0.4956
|
Bulk
modulus
(объёмный модуль
упругости)
|
29 kPa
(кПа)
|
Shear modulus
(модуль
сдвига)
|
20 kPa
(кПа)
|
Johnson
Cook
failure
(критерии разрушения по Джонсону-Куку)
|
|
Damage
constant
D1
(параметр поврежденности
D1)
|
−0.13549
|
Damage
constant
D2
(параметр поврежденности
D2)
|
0.6015
|
Damage
constant
D3
(параметр поврежденности
D3)
|
0.25892
|
Damage
constant
D4
(параметр поврежденности
D4)
|
0.030127
|
Damage
constant
D5
(параметр поврежденности
D5)
|
0
|
Melting temperature (температура
плавления)
|
20 °C
|
Reference
strain
rate
(/sec)
(эталонная скорость
деформации в /сек)
|
1
|
Визуализация
результатов исследования получена в постпроцессоре ANSYS
LS-DYNA.
В
эксперименте 1 при направлении преднагружении в проксимальной части бедренной
кости величиной 100 Н по оси
Y
в
результате моделирования ударного воздействия перпендикулярно к продольной оси
бедренной кости в условной средине диафиза бедренной кости сформировалась косопоперечная
асимметричная линия перелома с выкрашиванием элементов в линии перелома и
началом формирования линии перелома в зоне минимальной толщины компактного
костного вещества в травмируемой зоне проксимальнее места приложения
травмирующей силы бедренной кости на противоположной стороне от нагружения. Максимальные
эквивалентные напряжения достигали 52,958
MPa
(Рис. 1-3). При анализе векторов главных напряжений (Рис. 3) отмечается
динамическое распределение сжимающих (синие стрелки) и растягивающих (красные
стрелки) напряжений, меняющие свое расположение и величину в зависимости от
этапа разрушения кости. В начальном этапе нагружения в проксимальной части
кости возникали асимметричные напряжения, указывающие на винтообразную
деформацию кости, не превышающую пределы разрушения структур, возникающие при
преднагружении (Рис. 4).
|
Рис. 1. Визуализация конечного
этапа моделирования перелома диафиза бедренной кости в сложнонапряженном
состоянии при преднагрузке в проксимальной части бедренной кости по оси
Y.
Анализ по Фон Мизесу.
|
|
Рис. 2. Динамическая визуализация
моделирования перелома диафиза бедренной кости в сложнонапряженном состоянии
при преднагрузке в проксимальной части бедренной кости по оси
Y.
Анализ по Фон Мизесу.
|
|
Рис. 3. Динамическая визуализация
моделирования перелома диафиза бедренной кости в сложнонапряженном состоянии
при преднагрузке в проксимальной части бедренной кости по оси
Y.
Анализ по векторам главных напряжений.
|
|
Рис. 4. Визуализация внутренних
напряжений в проксимальном отделе при деформации кручения верхней части бедренной
кости при преднагрузке в проксимальной части бедренной кости по оси
Y.
|
В
эксперименте 2 при направлении преднагрузки в проксимальной части бедренной
кости величиной 100 Н в противоположном направлении оси
Y
в результате моделирования ударного воздействия перпендикулярно к продольной
оси бедренной кости в условной средине диафиза бедренной кости сформировалась
косопоперечная асимметричная линия перелома с выкрашиванием элементов в линии
перелома и началом формирования линии перелома в зоне минимальной толщины
компактного костного вещества в травмируемой зоне проксимальнее места
приложения травмирующей силы бедренной кости на противоположной стороне от
нагружения. Максимальные эквивалентные напряжения достигали 59,218
MPa
(Рис. 5-7). При анализе векторов главных напряжений (Рис. 7) отмечается
динамическое распределение сжимающих (синие стрелки) и растягивающих (красные
стрелки) напряжений, меняющие свое расположение и величину в зависимости от
этапа разрушения кости. В начальном этапе нагружения в проксимальной части
кости возникали асимметричные напряжения, указывающие на винтообразную
деформацию кости, не превышающую пределы разрушения структур, возникающие при
преднагружении (Рис. 8).
|
Рис. 5. Визуализация конечного
этапа моделирования перелома диафиза бедренной кости в сложнонапряженном
состоянии при преднагрузке в проксимальной части бедренной кости в
противоположном направлении оси
Y.
Анализ по Фон Мизесу.
|
|
Рис. 6. Динамическая визуализация
моделирования перелома диафиза бедренной кости в сложнонапряженном состоянии
при преднагрузке в проксимальной части бедренной кости в противоположном
направлении оси
Y. Анализ по Фон
Мизесу.
|
|
Рис. 7. Динамическая визуализация
моделирования перелома диафиза бедренной кости в сложнонапряженном состоянии
при преднагрузке в проксимальной части бедренной кости в противоположном
направлении оси
Y. Анализ по
векторам главных напряжений.
|
|
Рис. 8. Визуализация внутренних
напряжений в проксимальном отделе при деформации кручения верхней части
бедренной кости при преднагрузке в проксимальной части бедренной кости в
противоположном направлении оси
Y.
|
Для
валидации математической модели мы сравнили морфологии переломов, возникающих в
аналогичных условиях на биоманекенах и в практических наблюдениях, в которых линия
перелома возникала в зоне действия травмирующего предмета с формированием
косопоперечной асимметричной линии перелома с типичными признаками растяжения
костной ткани на противоположной стороне от места воздействия. В отдельных
наблюдениях формировалась дополнительная винтообразная линия перелома на проксимальной
части кости, которая ближе к зоне вращения кости [3, 10] (Рис. 9-11).
|
Рис. 9. Изображение линии перелома
диафиза бедренной кости в месте локального ударного травмирующего
воздействия.
|
|
Рис. 10. Изображение поверхности
перелома диафиза бедренной кости в месте локального ударного травмирующего
воздействия с выраженной асимметричностью.
|
|
Рис. 11. Изображение винтообразной
линии перелома проксимальной части диафиза бедренной кости вне места
локального ударного травмирующего воздействия.
|
В
описанных натурных экспериментах свойства костной ткани не изучались,
отмечалось лишь отсутствие какой-либо патологии костной ткани. Также в
экспериментах на биоманекенах и практических наблюдениях не изучалось изменение
толщины компактной костной пластинки на протяжении диафиза бедренной кости в
зоне травмирующего воздействия.
При
анализе экспериментальных данных и практических наблюдений нами было
предположено, что формирование винтообразной линии перелома в проксимальной
части бедренной кости возникает при превышении напряжений костной ткани её
прочностных характеристик. Поэтому было решено провести дополнительный
эксперимент с приложением большей силы на проксимальную часть бедренной кости в
преднагрузке и удалением жесткой подложки без изменений остальных граничных условий.
В
эксперименте 3 при направлении преднагрузки в проксимальной части бедренной
кости величиной 300 Н по оси
Y
в результате моделирования ударного воздействия перпендикулярно к продольной
оси бедренной кости в условной средине диафиза бедренной кости сформировалась
косопоперечная асимметричная линия перелома с выкрашиванием элементов в линии
перелома и началом формирования линии перелома в зоне минимальной толщины
компактного костного вещества в травмируемой зоне проксимальнее места
приложения травмирующей силы бедренной кости на противоположной стороне от
нагружения и винтообразная линия перелома в проксимальной части диафиза
бедренной кости. Максимальные эквивалентные напряжения достигали 119,95
MPa
(Рис. 12-14).
|
Рис. 12. Визуализация конечного
этапа моделирования перелома диафиза бедренной кости в сложнонапряженном
состоянии при увеличенной до 300 Н преднагрузке в проксимальной части
бедренной кости по оси
Y.
Анализ по Фон Мизесу.
|
|
Рис. 13. Динамическая визуализация
моделирования перелома диафиза бедренной кости в сложнонапряженном состоянии
при увеличенной до 300 Н преднагрузке в проксимальной части бедренной кости
по оси
Y. Анализ по Фон Мизесу.
|
|
Рис. 14. Динамическая визуализация
моделирования перелома диафиза бедренной кости в сложнонапряженном состоянии
при увеличенной до 300 Н преднагрузке в проксимальной части бедренной кости
по оси
Y. Анализ по векторам
главных напряжений.
|
При
сравнении морфологических особенностей переломов и их моделировании в нашем
экспериментальном исследовании обнаружено соответствие локализации и характера
переломов. Тем не менее, вариабельность морфологических особенностей переломов
диафизов бедренных костей в нашем экспериментальном исследовании методом
математического моделирования имела закономерность в зависимости от:
•
направления
вектора действия травмирующей силы и величины преднарузки,
•
места
приложения травмирующей силы (локальный перелом возникает в зоне приложения
травмирующей силы с зоной разрыва на противоположной стороне, винтообразный
перелом возникает при значительной деформации вращения ближе к зоне вращения),
•
толщины
компактной пластины в зоне травмирующего воздействия на противоположной стороне
от места приложения травмирующей силы.
Конечно-элементный
анализ позволяет визуализировать и прогнозировать напряжения, возникающие в костной
ткани при комбинации ударного и вращающего воздействия (сложнонапряженное
состояние). Полученные при моделировании данные подтверждаются результатами
оригинальных натурных экспериментов и практических наблюдений.
С
помощью математического моделирования выявлена зависимость морфологических
характеристик переломов от места приложения травмирующей силы на передней
поверхности бедра, толщины компактной пластины, места силы приложения вращающих
нагрузок. При сочетании ударного и вращающего воздействия в локальном перелома
возникают комбинация признаков угловой деформации с асимметрией плоскости
перелома, характерной для винтообразного перелома. При увеличении вращения
бедренной кости ближе к месту вращения возникает дополнительная винтообразная
линия перелома.
Проведенные
эксперименты показывают возможность применения метода конечных элементов в
судебной медицине, что позволяет достоверно предсказывать процесс разрушения
биологических объектов при различных видах механического воздействия и
визуализировать его результаты на любом этапе эксперимента с большим
количеством отображенных данных. В будущем возможно использование этой
информации для решения обратной задачи – определения трасологических свойств
травмирующего орудия на основе морфологической картины разрушения. Это
подтверждает высокую эффективность конечно-элементного анализа в судебной
медицине.
1.
Streubel PN, Gardner MJ, Ricci WM.
Management of femur shaft fractures in obese patients.
Orthopedic Clinics of
North America.
2011 Jan;42(1):21-35. https://doi.org/10.1016/j.ocl.2010.07.004
2.
Кислов
М.А., Бахметьев В.И., Кильдюшов Е.М., Крупин К.Н. Математическое моделирование
перелома диафиза бедренной кости под острым углом / Судебно-медицинская
экспертиза. 2022. № 65(6). С. 37-41.
http://doi.org/10.17116/sudmed20226506137.
3.
Бахметьев
B.И., Кислов М.А., Степанян Н.А. Особенности разрушений костной ткани при
различных способах внешнего воздействия / Системный анализ и управление в
биомедицинских системах. 2006. № 5(1). С. 192-195.
4.
Официальная
страница продукта
Ansys
LS-DYNA.
Дата обращения 15.05.22.
https://www.ansys.com/products/structures/ansys-ls-dyna
5.
Kopperdahl DL, Keaveny TM. Yield strain
behavior of trabecular bone.
J Biomech. 1998;
31(7):601-608.
https://doi.org/10.1016/s0021-9290(98)00057-8
6.
Yosibash Z, Tal D, Trabelsi N. Predicting
the yield of the proximal femur using high-order finite-element analysis with
inhomogeneous orthotropic material properties. Philosophical Transaction of the
Royal Society a Mathematical, Physical and Engineering Sciences.
2010;368(1920):2707-2723.
https://doi.org/10.1098/rsta.2010.0074
7.
Keller TS, Mao Z, Spengler DM. Young's
modulus, bending strength, and tissue physical properties of human compact
bone. Journal of Orthopaedic Research. 1990;8(4):592-603.
https://doi.org/10.1002/jor.1100080416
8.
Chirchir H. Limited Trabecular Bone
Density Heterogeneity in the Human Skeleton. Anatomy Research International,
2016,
1(1):7.
https://doi.org/10.1155/2016/9295383
9.
Бахметьев
В.И. Определение вида и направления внешних воздействий по морфологии
разрушения длинной трубчатой кости. Судебная медицина. 2019. №5(S1). С. 155.
10.
Бахметьев
В.И., Кислов М.А. Определение вида внешнего воздействия на основе анализа морфологии
излома длинных трубчатых костей нижних конечностей Судебно-медицинская
экспертиза. 2008. № 51(6). С. 11-14.
Mathematical Modeling and Visualization of a Complex Stress State in Case of a Fracture of the Femoral Diaphysis
Authors: K.N. Krupin1,A, M.A. Kislov2,B, V.I. Bahmetev3,C, E.M. Kildyushov4,A
A Pirogov Russian Research Medical University, Moscow, Russia
B Sechenov First Moscow State Medical University (Sechenov University), Moscow, Russia
C Voronezh State Medical University named after N.N. Burdenko, Voronezh, Russia
1 ORCID: 0000-0001-6999-8524, krupin_kn@rsmu.ru
2 ORCID: 0000-0002-9303-7640, smedik@gmail.com
3 ORCID: 0000-0002-8770-1664, bahmetev@vsmaburdenko.ru
4 ORCID: 0000-0001-7571-0312, kem1967@bk.ru
Abstract
The purpose of this work is to establish the possibility of using the finite element analysis method to study complex stress states in case of a femur fracture with subsequent data visualization. Experimental data were obtained on a solid-state mathematical parametric model of the femur, created on the basis of computer tomogram data, and repeating studies on native biological objects. As a result of mathematical modeling, oblique transverse and helical fractures of the diaphysis of the femur with elements of helical deformation were studied. The application of finite element analysis made it possible to visualize and predict the stresses arising in bone tissue under the impact of a blunt solid object in a complex stress state and the morphological features of femoral shaft fractures under different torsional loading forces of the proximal part of the femur. The data on the mechanism and morphology of the femoral shaft fracture obtained during modeling are confirmed by the results of original full-scale experiments.
Keywords: femur, impact, complex stress state, finite element analysis, forensic medicine.
1.
Streubel PN, Gardner
MJ, Ricci WM. Management of femur shaft fractures in obese patients.
Orthopedic
Clinics of North America. 2011 Jan;42(1):21-35.
https://doi.org/10.1016/j.ocl.2010.07.004
2.
Kislov M.A., Bakhmetyev
V.I., Kildyushov E.M., Krupin K.N. Matematicheskoye modelirovaniye pereloma
diafiza bedrennoy kosti pod ostrym uglom. Sudebno-meditsinskaja ekspertiza.
2022;65(6):37-41. [In Russian].
http://doi.org/10.17116/sudmed20226506137.
3.
Bakhmetyev V.I., Kislov M.A., Stepanyan N.A.
Osobennosti razrushenij kostnoj tkani pri razlichnyh sposobah vneshnego
vozdejstvija / Sistemnyj analiz i upravlenie v biomedicinskih sistemah. 2006.
No. 5(1). P. 192-195 [In Russian].
4.
Oficial'naja stranica
produkta Ansys LS-DYNA. Data obrashhenija 15.05.22.
https://www.ansys.com/products/structures/ansys-ls-dyna
5.
Kopperdahl DL, Keaveny
TM. Yield strain behavior of trabecular bone.
J Biomech.
1998;
31(7):601-608.
https://doi.org/10.1016/s0021-9290(98)00057-8
6.
Yosibash Z, Tal D,
Trabelsi N. Predicting the yield of the proximal femur using high-order
finite-element analysis with inhomogeneous orthotropic material properties.
Philosophical Transaction of the Royal Society a Mathematical, Physical and
Engineering Sciences. 2010;368(1920):2707-2723.
https://doi.org/10.1098/rsta.2010.0074
7.
Keller TS, Mao Z,
Spengler DM. Young's modulus, bending strength, and tissue physical properties
of human compact bone. Journal of Orthopaedic Research. 1990;8(4):592-603.
https://doi.org/10.1002/jor.1100080416
8.
Chirchir H. Limited
Trabecular Bone Density Heterogeneity in the Human Skeleton. Anatomy Research
International, 2016,
1(1):7.
https://doi.org/10.1155/2016/9295383
9.
Bahmetyev V.I.
Opredelenie vida i napravlenija vneshnih vozdejstvij po morfologii razrushenija
dlinnoj trubchatoj kosti.
Sudebnaja
medicina.
2019;5(S1):155. [In Russian].
10.
Bakhmetyev V.I., Kislov M.A. Opredeleniye
vida vneshnego vozdeystviya na osnove analiza morfologii izloma dlinnykh
trubchatykh kostey nizhnikh konechnostey. Sudebno-meditsinskaya ekspertiza.
2008;51(6):11-14.
[In
Russian].