DESTRUCTION SIMULATION FOR THE INHOMOGENEOUS BODY BY FINITE ELEMENT METHOD USING COMPUTED TOMOGRAPHY DATA

Abstract


The work is devoted to modelling the bone organ destruction using the finite-element method and taking into account the CT data. In order to evaluate the model, the simulation results for the femur are presented. Such an approach allows to individualize the simulation of the organ and determines the relevance of the study. Numerical studies using the finite element method were performed, an eight-node finite element with linear approximation was used. A linear formulation for an inhomogeneous elastic body was considered. Young’s modulus and ultimate stress were determined using power functions. The dependence of mechanical properties on the optical density was assumed, which, in turn, was determined by linear relationships depending on the Hounsfield numbers. For discretization by the finite element method, the distribution of the mechanical properties of the material for each element was performed according to tomography data. After solving, the stress-strain state problem safety factor was calculated in each node. Based on the calculated safety factor, a decision is made either to thicken the mesh or to remove the finite element. The paper presents the results of calculations using the method for two problems: an inhomogeneous sample and specimen with averaged mechanical properties. The numerical results in these problems illustrate the significant differences in the results of the stress-strain state of the organ and allow us to judge the nature of the destruction of the organ.

Full Text

Введение При оценке механических свойств костной ткани необходимо учитывать индивидуальные свойства и их распределение в изучаемых органах. Так, на свойства костной ткани оказывают большое влияние физиологические процессы [8, 9] и снижение физической активности [2, 9, 17, 18], которые могут внести существенный вклад в структуру костной ткани, что влечет за собой изменение механических свойств. В клинической практике для достижения хорошего качества лечения необходимо иметь информацию о состоянии костной ткани, в том числе такие данные могут определить и тактику лечения. В настоящее время в современной ортопедии и травматологии существует множество исследований, направленных на улучшение методик оперативного вмешательства [11-13], а также оптимизации используемых стержневых аппаратов и эндопротезов [1, 3-7, 20]. Несмотря на то, что существующие решения позволяют улучшить качество результатов оперативного вмешательства, в конечном счете они все равно зависят от полноты знаний о локальных механических свойствах органа. В связи с этим актуальной задачей становится задача об определении распределения механических свойств органа. Современные подходы к численному моделированию в сочетании с методами неразрушающего контроля позволили расширить класс решаемых задач [10, 26, 27, 34-36], более того, нашли развитие и численные методы, напрямую учитывающие структурные данные [29, 30] и даже позволяющие предсказать поведение органа в последующем [19, 24, 25]. Преобладающая доля работ, посвященных оценке механических свойств костной ткани, направлена на оценку упругих констант среды [15, 16, 21, 22, 28, 32, 33]. В настоящей работе предлагается подход к оценке прочности всего органа. Целью данного исследования является построение конечно-элементной модели, позволяющей c учетом распределения механических свойств (модуля упругости Юнга и предела прочности) по данным компьютерной томографии описать механизм разрушения органа под действием внешнего воздействия. Материалы и методы Постановка задачи Принятые допущения При построении модели мы предполагаем, что материал костной ткани является изотропным, линейным и негомогенным, а анизотропия возникает ввиду распределения пор по материалу. Более того, из предположения наличия связи между физической плотностью и механическими характеристиками использовалась оптическая плотность для расчета модуля упругости Юнга и предельных напряжений. Математическая модель Расчетная область включает в себя некоторое негомогенное тело. Его объем - V, а ¶V - граница этого множества (свободная поверхность). Механическое поведение системы, занимающей область V в R3 с границей ¶V, в рамках линейной теории упругости описывается следующей системой уравнений: , (1) , (2) , (3) где V°= VȶV; - вектор перемещений; - тензор напряжений; - тензор упругой деформации; - тензор упругих свойств; - некоторая функция, определяющая негомогенность среды. Используется подход, описанный ранее в работах [13, 15, 16], который предполагает, что допустимые напряжения зависят от пространственной координаты и их величина может быть определена по оптической плотности. Для удобства вводится коэффициент запаса прочности (4) где - напряжения по Мизесу в заданной точке, - допустимые напряжения в заданной точке. Для определения упругих констант и допустимых напряжений используются известные соотношения [23, 31], построенные на числах Хаунсфилда. Величины констант, входящих в данные соотношения, могут быть найдены из эксперимента. В общем виде эти зависимости можно выразить через числа Хаунсфилда и выписать в следующем виде: (5) (6) Здесь коэффициенты a и b с соответствующими индексами определяются из эксперимента и могут варьироваться в зависимости от вида костной ткани или органа. Идея подхода основывается на вычитании подобласти Vcrit расчетного объема V, в которой коэффициент запаса прочности более единицы. Численное решение задачи Численная реализация проводилась на основе метода конечных элементов. Подход для построения негомогенной модели подробно описан в прошлых работах [13, 15, 16]. Отметим лишь ключевые моменты метода. Дано дискретизированное множество конечных элементов, которое представляет собой данные о конечно-элементной сетке, а также величины коэффициента рентгеновского ослабления в узлах: (7) (8) где - конечно-элементная дискретизация объема - номер узла и его координаты; - номер элемента и образующие его номера узлов; - количество узлов и элементов соответственно. Для каждого элемента определяются осредненные механические характеристики по данным компьютерной томографии. Величины допустимых напряжений хранятся в отдельном массиве, который адресно связан с номерами элементов. Подход к моделированию разрушения основан на итерационном алгоритме, представленном ниже. Для обработки данных томографии был использован программный продукт Avizo, для конечно-элементного моделирования - Ansys. Для конечно-элементного моделирования напряженно-деформированного состояния органа применялся четырехузловой конечный элемент с линейной аппроксимацией Solid186. Используемые значения коэффициентов в выражениях (5)-(6) приведены в таблице. Значения параметров Параметр Модуль Юнга 46·103 2,00 Допустимые напряжения 60 1,82 Алгоритм Входные данные: сетка конечных элементов, граничные условия. Выходные данные: сетка конечных элементов, решение в узлах. 1: Решение задачи НДС 2: N = 1; Iter =1 3: Для элемента под номером N определяем: Минимальный размер элемента L, напряжения по Мизесу, осреднённые по элементам допустимых напряжений crL и расчетный коэффициента запаса kN, 4: Если kN ≥ 1, то Если L > crL, то Сохранение номера элемента в список rMesh Иначе Сохранение номера элемента в список kMesh 5: Если N ≤ nElem, то N = N + 1 переход на пункт 3 6: Сгущение всех элементов из списка rMesh 7: Удаление всех элементов из списка kMesh 8: Пересчет количества элементов в переменной nElem 9: Eсли Iter ≤ MaxIter, то Iter = Iter + 1 переход на пункт 1 Для конечно-элементной модели соотношения для определения коэффициента запаса прочности (5) можно переписать в виде (9) Анализ НДС органа в этом случае будем проводить по данным коэффициента запаса прочности (9). Модельная задача Для расчетов была рассмотрена бедренная кость кролика (рис. 1, а). Образец испытывался на изгиб. В этом случае перемещения узлов на одном торце фиксировались, а к узлам другого торца последовательно прикладывались перемещения. Перемещения распределялись по экспоненциальному закону в диапазоне (0, 0,5] мм со сгущением к правому концу интервала. Для оценки значимости предлагаемого подхода были рассмотрены две модели. Первая включала в себя негомогенность: напряженно-деформированное состояние определялось по описанной методике. Вторая основывалась на гомогенности - задача напряженно-деформированного состояния, где предполагаются следующие осреднённые по всему объему механические свойства: E = 32 ГПа, [σ]= 48 МПа. Образец бедренной кости подвергался сканированию на компьютерном томографе Vatech PaX-I 3D. Данные обрабатывались в программном комплексе Avizo. Для удовлетворения условия равномерности размеров конечных элементов была использована мелкая сетка (порядка 14·104 узлов и 7·105 элементов). Методика обработки данных компьютерной томографии и распределения механических свойств по органу подробно изложена в работе [13]. а б Рис. 1. Рассматриваемый образец: а - томография бедренной кости с выделенной расчетной областью; б - расчетная схема Далее проводились расчеты согласно описанному выше алгоритму, для оценки получаемых результатов на фиксированном торце определялось суммарное реактивное усилие. Результаты и обсуждение Для построенной конечно-элементной сетки производился импорт геометрии из программного комплекса Avizo в программный продукт Ansys. Каждому элементу назначался параметр, определяющий значение модуля Юнга и предел прочности согласно (12), (13) (на рис. 2 приведена сетка с цветовой диаграммой, указывающая распределение модуля Юнга и предельных напряжений по объему). Далее в Ansys с использованием макросов был реализован описанный в методах алгоритм. На каждом этапе производилась выгрузка номеров элементов, которые должны быть удалены, выгрузка реактивной силы в заделке, а также картины распределения напряжений по Мизесу. [σ] 16 20 24 28 33 38 43 48 54 60 МПа E 0,80 2 4 7 10 14 19 24 30 36 ГПа Рис. 2. Распределение механических параметров: модуля Юнга и предельных напряжений по конечно-элементной модели В расчетах было получено, что явное трещинообразование в случае как негомогенной, так и гомогенной модели реализуется при различных перемещениях. Для удобства оценки результатов использовалась реактивная сила в заделке. Так, для гомогенного случая сначала наблюдалось «выкашивание» материала (рис. 3, а) на внешней поверхности кости (при этом нагрузка составляла порядка 630 Н). Стоит отметить, что рост трещины, возникающей ближе к зоне заделки на определённом этапе расчета, останавливается. Вторая трещина, находящаяся подальше от заделки, развивается вплоть до разрушения, при этом наблюдается плавное снижение нагрузки до 550 Н (рис. 3, в). После отметки нагрузки 550 Н следует ее резкое падение, что характеризует разрушение костного органа. 0 13 26 40 54 67 80 94 107 121 а 0 42 84 126 169 211 253 295 338 380 б 0 51 103 155 206 258 310 361 413 464 в Рис. 3. Поле распределения напряжений по Мизесу (МПа): а - начальный этап; б - образование поперечной трещины; в - поперечная трещина в развитии Для негомогенного случая, начальная величина силы на тех же перемещениях составила порядка 820 Н (рис. 4, а), явного «выкрашивания» материала не наблюдалось. Локализация наибольших напряжений для негомогенного случая возникает в срединной поверхности кортикала (между внешней и внутренней поверхностями кости) рядом с областью заделки бедренной кости (рис. 4, б). Аналогично гомогенному случаю возникают две поперечные трещины, однако при нагрузке порядка 820 Н. Развитие трещин отличается от гомогенного случая. Рост трещины, возникающей ближе к зоне заделки, развивается вплоть до разрушения, при этом наблюдается плавное снижение нагрузки до 802 Н (рис. 4, в). После отметки величины нагрузки в 802 Н следует ее спад, т.е. наступает разрушение. 0 18 37 56 74 93 112 130 149 168 а 0 44 88 132 176 220 264 308 352 396 б 0 47 94 141 188 235 282 330 377 424 в Рис. 4. Поле распределений напряжений по Мизесу (МПа): а - начальный этап; б - образование поперечной трещины; в - поперечная трещина в развитии Рост второй трещины на определённом этапе расчета останавливается. Показано, что зона роста трещин и их характер роста для разных случаев отличаются. Для гомогенной модели разрушение начинается с внешней поверхности кости и сопровождается выкрашиванием (рис. 5, а), а для негомогенного случая разрушение начинается изнутри (рис. 5, б). На основе численных экспериментов были построены диаграммы в осях усилие-прогиб (рис. 5, в), величина усилия в этом случае оценивалась на основе данных, полученных в заделке. Ранее было отмечено, что характер падения нагрузки в обоих случаях отличается, как и нагрузка трещинообразования. Для гомогенной модели характерна меньшая нагрузка трещинообразования - 630 Н, при этом полное разрушение кости наблюдается при большем прогибе 0,51 мм и усилии 550 Н. Для негомогенного случая - при перемещении порядка 0,3 мм и усилии 810 Н. а б в Рис. 5. Трещина в гомогенном (а) и в негомогенном (б) образцах; диаграмма усилие-перемещение (в) Выводы 1. В работе приведена методика для моделирования разрушения органа негомогенной структуры на основе метода конечных элементов с учетом данных компьютерной томографии. 2. Численные результаты в приведенных задачах иллюстрируют значимые различия напряженно-деформированного состояния органа в гомогенной и негомогенной постановках и позволяют судить о локальной прочности костной ткани. 3. При моделировании принимались следующие допущения: а) негомогенный материал является изотропным; б) существует связь между физической плотностью и механическими характеристиками, что позволяет использовать оптическую плотность для расчета модуля упругости Юнга и предельных напряжений. 4. Применение описанной методики позволяет сделать вывод о характере разрушения органа при различных внешних нагрузках, с учетом распределения свойств материала. 5. Предлагаемая модель может быть расширена с учетом ортотропности материала.

About the authors

P. V Bolshakov

O. A Sachenkov

References

  1. Акулич Ю.В., Акулич А.Ю., Денисов А.С. Влияние количества и размеров резьбовых фиксаторов на адаптационные изменения механических свойств губчатой костной ткани и усилие сжатия отломков после контролируемого остеосинтеза // Российский журнал биомеханики. - 2012. - Т. 16, № 2. - С. 21-29.
  2. Акулич Ю.В., Акулич А.Ю., Денисов А.С., Шайманов П.С., Шулятьев А.Ф. Уточнение индивидуальной зависимости модуля упругости трабекулярной костной ткани от объемного содержания матрикса // Российский журнал биомеханики. - 2014. - Т. 18, № 2. - С. 158-167.
  3. Акулич Ю.В., Подгаец Р.М., Скрябин В.Л., Сотин А.В. Исследование напряженно-деформированного состояния эндопротезированного тазобедренного сустава // Российский журнал биомеханики. - 2007. - Т. 11, № 4. - С. 9-35.
  4. Анисимов О.Г., Ахтямов И.Ф., Шигаев Е.С. Особенности стационарного этапа лечения переломов проксимального отдела бедренной кости // ТаГраф. - 2017.
  5. Ахтямов И.Ф., Коваленко А.Н., Анисимов О.Г., Закиров Р.Х. Лечение остеонекроза головки бедра // Скрипта. - 2012.
  6. Ахтямов И.Ф, Хаертдинов И.С., Шигаев Е.С., Коваленко А.Н., Гатина Э.Б. Лечение пострадавших с переломами проксимального отдела бедренной кости в условиях Больницы скорой медицинской помощи // Современное искусство медицины. - 2013. - № 1 (9). - С.23-30.
  7. Ахтямов И.Ф., Шакирова Ф.В., Клюшкина Ю.А., Бакланова Д.А., Гатина Э.Б., Алиев Э.О. Анализ регенеративного процесса в области перелома большеберцовой кости // Травматология и ортопедия России. - 2016. - № 1 (79). - С.100-108.
  8. Банецкий М.В. Биомеханическое обоснование использования вертлужного компонента при эндопротезировании тазобедренного сустава: дис. … канд. мед. наук. - М., 2008. - 94 с.
  9. Белецкий А.В., Ахтямов И.Ф., Богосьян А.Б., Герасименко М.А. Асептический некроз головки бедренной кости у детей // Скрипта. - 2010.
  10. Иванов Д.В., Барабаш А.П., Барабаш Ю.А. Интрамедуллярный стержень нового типа для остеосинтеза диафизарных переломов бедра // Российский журнал биомеханики - 2012. - Т. 19, № 1. - С. 52-64.
  11. Измайлова З.Т. Предоперационная диагностика модульной трансформации при чрескостном остеосинтезе бедренной кости // Российский журнал биомеханики. - 2009. - Т. 13, № 2. - C. 93-98.
  12. Климов О.В. Расчет и контроль биомеханической оси нижней конечности во фронтальной плоскости при ее коррекции по Илизарову // Российский журнал биомеханики - 2014. - Т.18, № 2. - С. 239-247.
  13. Саченков О.А., Герасимов О.В., Королева Е.В., Мухин Д.А., Яикова В.В., Ахтямов И.Ф., Шакирова Ф.В., Коробейникова Д.А., Хань Х.Ч. Построение негомогенной конечно-элементной модели по данным компьютерной томографии // Российский журнал биомеханики. - 2018. - Т. 22, № 3. - С. 332-344.
  14. Саченков О.А., Хасанов Р.Ф., Андреев П.С., Коноплев Ю.Г. Численное исследование напряженно-деформированного состояния тазобедренного сустава при ротационной остеотомии проксимального участка бедренной кости // Российский журнал биомеханики. - 2016. - Т. 20, № 3. - C. 257-271.
  15. Харин Н.В., Воробьев О.В., Бережной Д.В., Саченков О.А. Методика построения репрезентативной модели по данным компьютерной томографии // Вестник Пермского национального исследовательского политехнического университета. Механика. - 2018. - № 3. - С. 95-102.
  16. Харин Н.В., Герасимов О.В., Большаков П.В., Хабибуллин А.А., Федянин А.О., Балтин М.Э., Балтина Т.В., Саченков О.А. Методика определения ортотропных свойств костного органа по данным компьютерной томографии // Российский журнал биомеханики. - 2019. - Т. 23, № 2. - С. 460-468.
  17. Baltina T.V., Ahmetov N.F., Sachenkov O.A., Fedyanin A.O., Lavrov I.A. The influence of hindlimb unloading on bone and muscle tissues in rat model // BioNanoSci. - 2017. - Vol. 7, no. 1. - P. 67-69. doi: 10.1007/s12668-016-0288-8
  18. Baltina T., Sachenkov O., Gerasimov O., Baltin M., Fedyanin A., Lavrov I. The influence of hindlimb unloading on the bone tissue’s structure // BioNanoScience. - 2018. - Vol. 8, no. 3. - P. 864-867. doi: 10.1007/s12668-018-0551-2
  19. Bobrovskaia A.S., Gavriushin S.S., Mitronin A.V. Evaluation of adhesive bond strength of dental fiber posts by “torque-out” test // Advances in Intelligent Systems and Computing. - 2020. - Vol. 902. - P. 261-270.
  20. Bolshakov P., Raginov I., Egorov V., Kashapova R., Kashapov R., Baltina T., Sachenkov O. Design and optimization lattice endoprosthesis for long bones: manufacturing and clinical experiment // Materials. - 2020. - Vol. 13. - P. 1185. doi: 10.3390/ma13051185
  21. Carniel T.A., Klahr B., Fancello E.A. A multiscale numerical approach for the finite strains analysis of materials reinforced with helical fibers // Mechanics of Materials. - 2018. - Vol. 126. - P. 75-85.
  22. Carniel T.A., Klahr B., Fancello E.A. On multiscale boundary conditions in the computational homogenization of an RVE of tendon fascicles // Journal of the Mechanical Behavior of Biomedical Materials. - 2019. - Vol. 91. - P. 131-138.
  23. Cuppone M., Seedhom B.B., Berry E., Ostell A.E. The Longitudinal Young’s Modulus of Cortical Bone in the Midshaft of Human Femur and its Correlation with CT Scanning Data // Calcified Tissue International. - 2004. -Vol. 74, no. 3. - P. 302-309.
  24. Demishkevich E., Gavriushin S. Simulation of the long-term orthodontic tooth movement using finite-element analysis of viscoelastic model for periodontal ligament // Series on Biomechanics. - 2018. - Vol. 32, no. 4. - P. 56-62.
  25. Demishkevich E.B., Gavriushin S.S. A viscoelastic model of the long-term orthodontic tooth movement // Advances in Intelligent Systems and Computing. - 2020. - Vol. 902. - P. 315-322.
  26. Eggermont F., Derikx L.C., Free J., van Leeuwen R., van der Linden Y.M., Verdonschot N., Tanck E. Effect of different CT scanners and settings on femoral failure loads calculated by finite element models // Journal of Orthopaedic Research. - 2018. - Vol. 36, iss. 8. - P. 2288-2295.
  27. Eggermont F., Derikx L.C., Verdonschot N., Van Der Geest I.C.M., De Jong M.A.A., Snyers A., Van der Linden Y.M., Tanck E. Can patient-specific finite element models better predict fractures in metastatic bone disease than experienced clinicians // Bone and Joint Research. - 2018. - Vol. 7, iss. 6. - P. 430-439.
  28. Gabidullin M.G., Kayumov R.A., Rakhimov R.Z., Temlyakov A.V. Inter-relation between structures and heat-transfer properties of porous ceramic materials // Stroitel'nye Materialy. - 2005. - Vol. 9. - P. 62-66.
  29. Gerasimov O.V., Berezhnoi D.V., Bolshakov P.V., Statsenko E.O., Sachenkov O.A. Mechanical model of a heterogeneous continuum based on numerical-digital algorithm processing computer tomography data // Russian Journal of Biomechanics. - 2019. - Vol. 23, no. 1. - P. 87-97. doi: 10.15593/RJBiomech/2019.1.10
  30. Giovannelli L., Rodenas J.J., Navarro-Jimenez J.M., Tur M. Direct medical image-based Finite Element modelling for patient-specific simulation of future implants // Finite Elements in Analysis and Design. - 2017. - Vol. 136. - P. 37-57.
  31. Gupta S., Dan P. Bone geometry and mechanical properties of the human scapula using computed tomography data // Trends Biomater. Artif. Organs. - 2004. - Vol. 17, no. 2. - P. 61-70.
  32. Hettich G., Schierjott R.A., Ramm H., Graichen H., Jansson V.b, Rudert M., Traina F., Grupp T.M. Method for quantitative assessment of acetabular bone defects // Journal of Orthopaedic Research. - 2018. doi: 10.1002/jor.24165.
  33. Kayumov R.A., Muhamedova I.Z., Tazyukov B.F., Shakirzjanov F.R. Parameter determination of hereditary models of deformation of composite materials based on identification method // Journal of Physics: Conference Series. - 2018. - Vol. 973, no. 1. - P. 012006.
  34. Marcián P., Florian Z., Horáčková L., Kaiser J., Borák L. Microstructural finite-element analysis of influence of bone density and histomorphometric parameters on mechanical behavior of mandibular cancellous bone structure // Solid State Phenomena. - 2017. - Vol. 258. - P. 362-365.
  35. Marcián P., Wolff J., Horáčková L., Kaiser J., Zikmund T., Borák L. Micro finite element analysis of dental implants under different loading conditions // Computers in Biology and Medicine. - 2018. - Vol. 96. - P. 157-165.
  36. Ridwan-Pramana A., Marcian P., Borak L., Narra N., Forouzanfar T., Wolff J. Finite element analysis of large PMMA skull reconstructions: a multi-criteria evaluation approach // PLoS ONE. - 2017. - Vol. 12. - e0179325. doi: 10.1371/journal.pone.0179325.

Statistics

Views

Abstract - 67

PDF (Russian) - 36

Refbacks

  • There are currently no refbacks.

Copyright (c) 2022 Российский журнал биомеханики

This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies