CC BY
23
DOI: 10.15593/RZhBiomeh/2019.3.11 УДК 531/534: [57+61], 539.3
МЕТОДИКА ОПРЕДЕЛЕНИЯ ОРТОТРОПНЫХ СВОЙСТВ КОСТНОГО ОРГАНА ПО ДАННЫМ КОМПЬЮТЕРНОЙ ТОМОГРАФИИ
Н.В. Харин1, О.В. Герасимов1, П.В. Большаков1, А.А. Хабибуллин1, А.О. Федянин2, М.Э. Балтин2, Т.В. Балтина2, О.А. Саченков1,3
1 Кафедра теоретической механики института математики и механики им Н.И. Лобачевского Казанского (Приволжского) федерального университета, Казань, 420008, ул. Кремлевская, 18, e-mail: 4works@bk.ru
2 Кафедра физиологии человека и животных института фундаментальной медицины и биологии Казанского (Приволжского) федерального университета, Казань, 420008, ул. Карла Маркса, 74, e-mail: tvbaltina@gmail.com
3 Кафедра основ конструирования Казанского национального исследовательского технического университета им. А.Н. Туполева, Казань, 420111, ул. Карла Маркса, 10
Аннотация. Актуальной задачей современной биомеханики костной ткани является анализ ее упругих свойств. Развитие компьютерной томографии вместе с численными методами позволяет использовать численные эксперименты для определения механических характеристик биологических объектов. В работе представлена методика, основанная на гипотезе, что анизотропия костной ткани вызвана неравномерностью распределения пор по объему. В предлагаемом подходе исследуемый орган предлагается разбить на набор репрезентативных объёмов, для каждого из них проводятся численные эксперименты, которые позволяют определить компоненты тензора упругих характеристик, направление осей ортотропии. Описанная методика применена для определения характеристик диафиза бедренной костной крысы. Показано, что получаемые упругие константы можно считать постоянными в органе, а тело - трансверсально изотропным.
Ключевые слова: костная ткань, ортотропия, численные эксперименты, репрезентативный объем.
Введение
В настоящее время актуальной задачей выступает моделирование напряженно-деформированного состояния пористых многофазных сред, одной из которых является костная ткань. При моделировании поведения костной ткани [9, 25, 26, 30, 33] необходимо учитывать анизотропные свойства среды, определяемые особенностями различных видов костных органов в зависимости от их функционала, возраста, пола, среды обитания и т.п.; таким образом, данная задача остаётся актуальной. В случае моделирования адаптационных процессов [10-13, 20, 34], возникающих в костной ткани под воздействием внешних силовых воздействий, необходимо знать начальное
© Харин Н.В., Герасимов О.В., Большаков П.В., Хабибуллин А.А., Федянин А.О., Балтин М.Э., Балтина Т.В., Саченков О.А., 2019
Харин Никита Вячеславович, лаборант-исследователь кафедры теоретической механики, Казань Герасимов Олег Владимирович, лаборант-исследователь кафедры теоретической механики, Казань Большаков Павел Владиславович, лаборант-исследователь кафедры теоретической механики, Казань Хабибуллин Азат Айдарович, аспирант кафедры теоретической механики, Казань Федянин Артур Олегович, аспирант кафедры физиологии человека и животных, Казань Балтин Максим Эдуардович, аспирант кафедры физиологии человека и животных, Казань Балтина Татьяна Валерьевна, к.б.н., доцент кафедры физиологии человека и животных, Казань Саченков Оскар Александрович, к.ф.-м.н., научный сотрудник кафедры теоретической механики, Казань
состояние костного органа, а значит, распределение в теле анизотропных свойств. Для определения механических свойств пористого материала необходимо проведение натурных экспериментов, но данный метод является дорогостоящим и требует большого количества однотипных образцов; более того, его использование невозможно без разрушения самого костного органа.
На сегодняшний день развитие компьютерной томографии позволяет получать информацию о структуре образца без разрушения, а распространённость компьютерных томографов в клинике расширяет область применения методики в задачах биомеханики. Развитие мощностей компьютерных томографов расширяет область применения метода на микроуровне (микробиология, кристаллография и т.п.) и макроуровне (медицина, машиностроение, нефтяная промышленность и т.п.). Вкупе с современными численными методами это позволяет решать задачи для структурно сложных систем и проводить численные эксперименты.
Целью данной работы является создание методики анализа распределения ортотропных свойств костного органа по данным компьютерной томографии.
Материалы и методы
Компьютерная томография образца представляет собой цифровой прототип исследуемого объекта, где каждое число суть нормированный по шкале Хаунсфилда коэффициент рентгеновского ослабления. Определяя порог бинаризации, можно дифференцировать материал и воздух (поры). В полученном бинаризированном массиве представляется возможным определение замыкающих поверхностей, которые формируют объем объекта. При предлагаемом подходе предлагается дискретизировать исследуемый объем на репрезентативные элементы, для каждого из которых в последующем необходимо определить механические свойства в численных экспериментах. Получаемый набор механических констант для репрезентативных объемов подвергается статистическому анализу с целью восстановления аппроксимирующих функций.
Из предположения, что анизотропия костной ткани определяется неравномерностью распределения пор по объему [15, 22, 32], в расчетах можно предположить, что сам материал костной ткани изотропен. В дальнейшем применяется метод представительных объемов, нашедший распространение в механике композиционных материалов и достаточно разработанный в настоящее время [1-6, 19, 27, 35]. Кроме того, он на основе данных компьютерной томографии применяется в геологии и биомеханике [28]. Существуют оценки адекватности получаемых результатов в сравнении с другими методами [16, 21, 29, 36] и для различных моделей компьютерных томографов [14]. Определить напряженно-деформированное состояние исследуемого объема можно с помощью метода граничных элементов [17] и метода конечных элементов [18, 24, 31, 36]. Последний используется в данной работе.
В каждом представительном элементе можно построить регулярную конечно-элементную сетку по данным томографии: в этом случае физический размер элемента соответствует размеру вокселя томографии. Полученная конечно-элементная сетка трансформируется согласно бинаризированным данным компьютерной томографии. Для определения порога бинаризации в работе был использован метод Отцу. Затем для полученного конечно-элементного ансамбля проводятся численные эксперименты: три на осевое сжатие и три на чистый сдвиг (рис. 1).
После решения задачи напряженно-деформированного состояния компоненты напряжения усредняются:
1 г
-¡д}сМ = С-Щ, (1)
где а - вектор напряжений; 8 - вектор деформаций; ] - индекс в нотации Фойгта.
Таким образом, в результате численных экспериментов можно получить 36 линейных уравнений, в то время как неизвестных - 21. Для решения переопределенной системы применяется метод наименьших квадратов.
а б
Рис. 1. Схема нагружения в численных экспериментах: а - одноосное сжатие; б - сдвиг
На следующем шаге необходимо определить оси ортотропии и величины механических характеристик, для чего компоненты тензора упругих констант выписываются относительно трех углов поворота [7], по которым минимизируется целевая функция:
36 -63 ~ 6 ~
Е (С) + 2Е Е (С) +4 Е (С)
у _ '=1 1=4_1=4 ] =1_г,] =4,]Фг__^2)
Е (с; )2+4 Е (с; )2 '
г, 1=1 1,1=4
Подробно методика описана в работе [8]. В результате расчетов для каждого репрезентативного элемента можно определить оси ортотропии и механические характеристики в этих осях. Полученные результаты представляют собой тензорное поле над объемом объекта. Для практических расчетов удобно использовать аппроксимации полученных механических свойств.
Результаты и обсуждение
Приведенная методика была реализована для определения характеристик диафиза бедренной кости крысы. Сканирование образца выполнено с применением микро-/нанофокусной системы рентгеновского контроля для компьютерной томографии и 2.0-инспекции Phoenix V\tome\X S240 в лаборатории рентгеновской компьютерной томографии Института геологии и нефтегазовых технологий Казанского (Приволжского) федерального университета. Система оснащёна двумя рентгеновским трубками: микрофокусной с максимальным ускоряющим напряжением 240 кВ, мощностью 320 Вт и нанофокусной с максимальным ускоряющим напряжением 180 кВ, мощностью 15 Вт. Для первичной обработки данных и создания объёмной (воксельной) модели образца на базе рентгеновских снимков (проекций) используется программное обеспечение datos\x reconstruction. Зафиксированный
в держателе образец помещается на вращающийся столик камеры рентгеновского компьютерного томографа на оптимальном расстоянии от источника рентгеновского излучения. Для полученных образцов была использована нанофокусная трубка. Съёмка проводилась при ускоряющем напряжении 90-100 кВ и токе 140-150 мА.
Диафиз был разделен в поперечном сечении на 10 репрезентативных элементов и на 22 репрезентативных элемента в продольном направлении (всего 220 репрезентативных элементов). На рис. 2, а приведены данные томографии, на рис. 2, б - разбиение на репрезентативные элементы. Каждый репрезентативный элемент представляет собой куб с размером 42*42*42 вокселя, размер вокселя -9,73 мкм.
Материал костной ткани задавался следующими характеристиками: модуль Юнга, равный 20 ГПа, и коэффициент Пуассона, равный 0,3 [23, 32]. После задания свойств материала проводились численные эксперименты по описанной методике [8]: определялись оси ортотропии и константы тензора упругости. На рис. 2, в представлены полученные направления осей ортотропии: для наглядности срединная поверхность образца аппроксимирована полигонами. Одно направление ортотропии совпадало с продольным направлением, два других располагались в поперечной плоскости.
Для анализа характера распределения упругих констант в пространстве весь набор данных был сгруппирован по продольным и поперечным направлениям. После этого в каждой из групп оценивалось распределение механических характеристик. Для обеих групп было получено нормальное распределение, а значит, для описания состоятельны средняя величина и стандартное отклонение. При анализе модуля Юнга по оси Х в поперечном направлении (рис. 3, а) величина составила 1590,42 ± 6,29 МПа (0,39%), в продольном (см. рис. 3, б) - 1597,37 ± 6,62 МПа (0,41%). При анализе модуля Юнга по оси У в поперечном направлении (см. рис. 3, в) -1595,8 ± 6,68 МПа (0,41%), в продольном (см. рис. 3, г) - 1595,8 ± 6,93 МПа (0,43%). При анализе модуля Юнга по оси 2 в поперечном направлении (см. рис. 3, д) величина составила 12254,3 ± 34,16 МПа (0,27%), в продольном (см. рис. 3, е) -12254,3 ± 19,34 МПа (0,15%).
а б в
Рис. 2. Компьютерная томография костного органа (а), дискретизация на репрезентативные элементы (б) и полученные направления ортотропии (в)
д е
Рис. 3. Распределение средних значений модуля Юнга, МПа, по слоям в определенных направлениях: а - по оси Х в продольном направлении; б - по оси Х в продольном направлении; в - по оси У в продольном направлении; г - по оси У в поперечном направлении; д - по оси 2 в продольном направлении; е - по оси 2 в поперечном направлении
б
а
Сравнивая среднее значение модулей Юнга по всем слоям в направлениях X и У, получаем, что абсолютная разница между ними не превышает 0,098%, а значит, можно заключить, что образец обладает трансверсальной изотропией. Аналогично проделав описанные операции для модуля сдвига и коэффициента Пуассона, можно показать, что упругие константы не зависят от пространственной координаты. В результате расчетов было получено, что отклонение от среднего не превышает 0,33%. Таким образом, модуль сдвига в плоскости ХУ составил 24640 ± 41 МПа, а в двух других направлениях - 1726 ± 3 МПа. Осредненные значения упругих констант приведены в таблице.
Осредненные значения упругих констант
Упругая константа Индекс компоненты
11 22 33 12 23 13
E, МПа 1590 1590 12 205 -
G, МПа - 24542 1718 1718
V 0,87 0,021 0,021
Для наглядности на рис. 4 приведены поверхность тензора упругих констант и ее сечения плоскостями.
в г
Рис. 4. Тензорная поверхность (а) и ее сечения: б - в плоскостиXY; в - в плоскости YZ; г - в плоскости XZ
Заключение
1. В работе представлена методика для определения распределения упругих констант образца по данным томографии. В используемой методике изучаемый объект предлагается разбить на набор репрезентативных объемов, на основе данных томографии для каждого объема производится построение конечно-элементной модели
и проводятся численные эксперименты. В результате можно получить распределение упругих констант в объекте. В дальнейшем для оценки получаемого распределения предлагается использовать методы статистического анализа.
2. В работе рассмотрен диафиз бедренной кости крысы, для описания которого использована данная методика. В результате расчетов было показано, что образец обладает трансверсальной изотропией, не зависящей от пространственной координаты; были определены средние величины упругих констант и их стандартные отклонения.
3. Результаты, полученные при применении описанной методики, могут быть улучшены путем увеличения количества численных экспериментов и регуляризаций задачи определения компонент тензора упругих констант.
4. Описанный подход может быть расширен не только на объекты из костной ткани, но в этом случае необходима некоторая априорная информация о структурных составляю щих.
Благодарности
Работа выполнена при финансовой поддержке Российского научного фонда, проект № 18-75-10027.
Список литературы
1. Каюмов Р.А. Расширенная задача идентификации механических характеристик материалов по результатам испытаний конструкций из них // Известия РАН. Механика твердого тела. - 2004. -№ 2. - С. 94-103.
2. Каюмов Р.А. Связанная задача расчета механических характеристик материалов и конструкций из них // Известия РАН. Механика твердого тела. - 1999. - № 6. - С. 118-127.
3. Каюмов Р.А., Нежданов Р.О., Тазюков Б.Ф. Определение характеристик волокнистых композитных материалов методами идентификации. - Казань: Изд-во Казанск. гос. ун-та, 2005. - 258 с.
4. Каюмов Р.А., Страхов Д.Е., Шакирзянов Ф.Р., Гимранов Л.Р., Мангушева А.Р. Идентификация жесткостных характеристик композита // Вестник Технологического университета. - 2016. - Т. 19, № 24. - С. 109-112.
5. Каюмов Р.А., Луканкин С.А., Паймушин В.Н., Холмогоров С.А. Идентификация механических характеристик армированных волокнами композитов // Ученые записки Казанского университета. Серия: Физико-математические науки. - 2015. - Т. 157, № 4. - С. 112-132.
6. Левандовский А.Н., Мельников Б.Е., Шамкин А.А. Моделирование пористого материала методом конечных элементов // Строительство уникальных зданий и сооружений. - 2017. - № 2 (53). -С. 61-77.
7. Лехницкий С.Г. Теория упругости анизотропного тела. - Изд. 2-е, перераб. и доп. - М.: Наука, 1977. - 416 с.
8. Харин Н.В., Воробьев О.В., Бережной Д.В., Саченков О.А. Методика построения репрезентативной модели по данным компьютерной томографии // Вестник Пермского национального исследовательского политехнического университета. Механика. - 2018. - № 3. - С. 95-102.
9. Саченков О.А., Герасимов О.В., Королева Е.В., Мухин Д.А., Яикова В.В., Ахтямов И.Ф., Шакирова Ф.В., Коробейникова Д.А., Хань Х.Ч. Построение негомогенной конечно-элементной модели по данным компьютерной томографии // Российский журнал биомеханики. - 2018. - Т. 22, № 3. - С. 332-344.
10. Чикова Т.Н., Киченко А.А., Тверье В.М., Няшин Ю.И. Биомеханическое моделирование трабекулярной костной ткани в состоянии равновесия // Российский журнал биомеханики. - 2018. -Т. 22, № 3. - С. 282-291.
11. Чикова Т.Н., Киченко А.А., Тверье В.М., Няшин Ю.И. Моделирование перестройки трабекулярной костной ткани в ветви нижней челюсти человека // Российский журнал биомеханики. - 2018. - Т. 22, № 3. - С. 292-300.
12. Baltina T., Ahmetov N., Sachenkov O., Fedyanin A., Lavrov, I. The Influence of Hindlimb Unloading on Bone and Muscle Tissues in Rat Model // BioNanoScience - 2017. - Vol. 7, iss. 1. - P. 67-69.
13. 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, iss. 3. - P. 864-867.
14. 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.
15. 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.
16. Faisal T.F., Awedalkarim A., Chevalier S., Jouini M.S., Sassi M. Direct scale comparison of numerical linear elastic moduli with acoustic experiments for carbonate rock X-ray CT scanned at multi-resolutions // Journal of Petroleum Science and Engineering. - 2017. - Vol. 152. - P. 653-663. DOI: 10.1016/j.petrol.2017.01.025
17. Fernandes G.R., Ohland G.A., Vieira J.F. A boundary element formulation to perform elastic analysis of heterogeneous microstructures // Engineering Analysis with Boundary Elements. - 2018. - Vol. 87. -P. 47-65. DOI: 10.1016/j.enganabound.2017.11.006
18. Gagliardi D., Sansalone V., Desceliers C., Naili S. Estimation of the effective bone-elasticity tensor based on ^CT imaging by a stochastic model. A multi-method validation // European Journal of Mechanics, A/Solids. - 2018. - Vol. 69. - P. 147-167. DOI: 10.1016/j.euromechsol.2017.10.004
19. Gerasimov O., Kharin N., Vorob'yev O., Semenova E., Sachenkov O. Determination of the mechanical properties distribution of the sample by tomography data // IOP Conf. Series: Journal of Physics: Conf. Series. - 2019. - Vol. 1158. - 022046. DOI: 10.1088/1742-6596/1158/2/022046
20. Gerasimov O., Yaikova V., Baltina T., Baltin M., Fedyanin A., Zamaliev R., Sachenkov O. Modeling the change in the stiffness parameters of bone tissue under the influence of external loads // IOP Conf. Series: Journal of Physics: Conf. Series. - 2019. - Vol. 1158. - 022045. DOI:10.1088/1742-6596/1158/2/022045
21. Greatrex F., Montefiori E., Grupp T., Kozak J., Mazza C. Reliability of an integrated ultrasound and stereophotogrammetric system for lower limb anatomical characterisation // Applied Bionics and Biomechanics. - 2017. - Vol. 2017. - Article ID 4370649. - 8 p. DOI: 10.1155/2017/4370649
22. Gross T., Pahr D.H., Peyrin F., Zysset P.K. Mineral heterogeneity has a minor influence on the apparent elastic properties of human cancellous bone: a SRmCT-based finite element study // Computer Methods in Biomechanics and Biomedical Engineering. - 2012. - Vol. 15, № 11. - P. 1137-1144.
23. Gross T., Pahr D.H., Zysset P.K. Morphology-elasticity relationships using decreasing fabric information of human trabecular bone from three major anatomical locations // Biomech. Model Mechanobiol. - 2013. -Vol. 12. - P. 793-800.
24. Hage I.S., Hamade R.F. Micro FEM simulations of single-cutting-edge sawing of cortical bone // ASME International Mechanical Engineering Congress and Exposition, Proceedings (IMECE). - 2013. - Vol. 3 A. DOI: 10.1115/IMECE2013-62589
25. 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
26. Jeanneret R., Arson C., Vennat E. Homogenization of dentin elastic properties based on microstructure characterization, statistical back-analysis, and FEM simulation // Poromechanics. - 2017. - P. 1339-1346. DOI: 10.1061/9780784480779.166
27. Kharin N., Vorob'yev O., Bolshakov P., Sachenkov O. Determination of the orthotopic parameters of a representative sample by computed tomography // IOP Conf. Series: Journal of Physics: Conf. Series. -2019. - Vol. 1158. - 032012. DOI:10.1088/1742-6596/1158/3/032012
28. Knowles N.K., Ip K., Ferreira L.M. The effect of material heterogeneity, element type, and down-sampling on trabecular stiffness in micro finite element models // Annals of Biomedical Engineering. - 2019. -Vol. 47, iss. 2. - P. 615-623. DOI: 10.1007/s10439-018-02152-6
29. Liu B., Lu S., Wu Y., Xing Y. Three dimensional micro/macro-mechanical analysis of the interfaces of composites by a differential quadrature hierarchical finite element method // Composite Structures. -2017. - Vol. 176. - P. 654-663. DOI: 10.1016/j.compstruct.2017.05.068
30. Marcian P., Florian Z., Horackova L., Kaiser J., Borak 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, SSP. - P. 362-365.
31. Marcian P., Wolff J., Horackova L., Kaiser J., Zikmund T., Borak L. Micro finite element analysis of dental implants under different loading conditions // Computers in Biology and Medicine. - 2018. - Vol. 96, № 1. - P. 157-165.
32. Maquer G., Musy S.N., Wandel J., Gross T., Zysset P.K. Bone volume fraction and fabric anisotropy are better determinants of trabecular bone stiffness than other morphological variables // Journal of Bone and Mineral Research. - 2015. - Vol. 30, № 6. - P. 1000-1008. DOI: 10.1002/jbmr.2437
33. Ridwan-Pramana A., Marcian P., Borak L., Narra N., Forouzanfar T., Wolff J. Finite element analysis of 6 large PMMA skull reconstructions: A multi-criteria evaluation approach // PLoS ONE. - 2017. - Vol. 12, iss. 6. - P. e0179325.
34. Shigapova F.A., Gerasimov O.V., Sachenkov O.A. The evolution of the trabecular bone at a constant combined loading // International Journal of Pharmacy and Technology. - 2016. - Vol. 8, № 4. -P.24261-24271.
35. Sokolowski D., Kaminski M. Homogenization of carbon/polymer composites with anisotropic distribution of particles and stochastic interface defects // Acta Mechanica. - 2018. - Vol. 229, iss. 9. - P. 3727-3765. DOI: 10.1007/s00707-018-2174-7
36. Wang K., Sun W. An updated Lagrangian LBM-DEM-FEM coupling model for dual-permeability fissured porous media with embedded discontinuities // Computer Methods in Applied Mechanics and Engineering. - 2019. - Vol. 344. - P. 276-305. DOI: 10.1016/j.cma.2018.09.034
TECHNIQUE FOR DETERMINING THE ORTHOTROPIC PROPERTIES OF THE BONE ORGAN ACCORDING TO COMPUTER TOMOGRAPHY
N.V. Kharin, O.V. Gerasimov, P.V. Bolshakov, A.A. Khabibullin, A.O. Fedyanin, M.E. Baltin, T.V. Baltina, O.A. Sachenkov (Kazan, Russia)
The actual task of modern biomechanics of bone tissue is to determine the elastic properties. Nowadays computed tomography coupled with computer technology allows the use of numerical experiments to determine the mechanical properties of biological objects. This paper presents such a technique. The technique is based on the hypothesis that the anisotropy of bone tissue is determined by the uneven distribution of pore volume. In the proposed approach, the organ under study is proposed to be divided into a set of representative volumes. A set of numerical experiments is carried out for each volume, which makes it possible to determine the components of the elastic tensor. Then, the axes of orthotropy and elastic constants are determined. The described technique is applied for diaphysis of the femoral bone rat. It is shown that the resulting elastic constants remain constant in the organ, and the body is transversally isotropic.
Key words: bone tissue, orthotropy, numerical experiments, representative volume.
Получено 9 апреля 2019
