Full Text

Введение В настоящее время актуальной задачей выступает моделирование напряженно-деформированного состояния пористых многофазных сред, одной из которых является костная ткань. При моделировании поведения костной ткани [9, 25, 26, 30, 33] необходимо учитывать анизотропные свойства среды, определяемые особенностями различных видов костных органов в зависимости от их функционала, возраста, пола, среды обитания и т.п.; таким образом, данная задача остаётся актуальной. В случае моделирования адаптационных процессов [10-13, 20, 34], возникающих в костной ткани под воздействием внешних силовых воздействий, необходимо знать начальное состояние костного органа, а значит, распределение в теле анизотропных свойств. Для определения механических свойств пористого материала необходимо проведение натурных экспериментов, но данный метод является дорогостоящим и требует большого количества однотипных образцов; более того, его использование невозможно без разрушения самого костного органа. На сегодняшний день развитие компьютерной томографии позволяет получать информацию о структуре образца без разрушения, а распространённость компьютерных томографов в клинике расширяет область применения методики в задачах биомеханики. Развитие мощностей компьютерных томографов расширяет область применения метода на микроуровне (микробиология, кристаллография и т.п.) и макроуровне (медицина, машиностроение, нефтяная промышленность и т.п.). Вкупе с современными численными методами это позволяет решать задачи для структурно сложных систем и проводить численные эксперименты. Целью данной работы является создание методики анализа распределения ортотропных свойств костного органа по данным компьютерной томографии. Материалы и методы Компьютерная томография образца представляет собой цифровой прототип исследуемого объекта, где каждое число суть нормированный по шкале Хаунсфилда коэффициент рентгеновского ослабления. Определяя порог бинаризации, можно дифференцировать материал и воздух (поры). В полученном бинаризированном массиве представляется возможным определение замыкающих поверхностей, которые формируют объем объекта. При предлагаемом подходе предлагается дискретизировать исследуемый объем на репрезентативные элементы, для каждого из которых в последующем необходимо определить механические свойства в численных экспериментах. Получаемый набор механических констант для репрезентативных объемов подвергается статистическому анализу с целью восстановления аппроксимирующих функций. Из предположения, что анизотропия костной ткани определяется неравномерностью распределения пор по объему [15, 22, 32], в расчетах можно предположить, что сам материал костной ткани изотропен. В дальнейшем применяется метод представительных объемов, нашедший распространение в механике композиционных материалов и достаточно разработанный в настоящее время [1-6, 19, 27, 35]. Кроме того, он на основе данных компьютерной томографии применяется в геологии и биомеханике [28]. Существуют оценки адекватности получаемых результатов в сравнении с другими методами [16, 21, 29, 36] и для различных моделей компьютерных томографов [14]. Определить напряженно-деформированное состояние исследуемого объема можно с помощью метода граничных элементов [17] и метода конечных элементов [18, 24, 31, 36]. Последний используется в данной работе. В каждом представительном элементе можно построить регулярную конечно-элементную сетку по данным томографии: в этом случае физический размер элемента соответствует размеру вокселя томографии. Полученная конечно-элементная сетка трансформируется согласно бинаризированным данным компьютерной томографии. Для определения порога бинаризации в работе был использован метод Отцу. Затем для полученного конечно-элементного ансамбля проводятся численные эксперименты: три на осевое сжатие и три на чистый сдвиг (рис. 1). После решения задачи напряженно-деформированного состояния компоненты напряжения усредняются: (1) где σ - вектор напряжений; ε - вектор деформаций; j - индекс в нотации Фойгта. Таким образом, в результате численных экспериментов можно получить 36 линейных уравнений, в то время как неизвестных - 21. Для решения переопределенной системы применяется метод наименьших квадратов. а б Рис. 1. Схема нагружения в численных экспериментах: а - одноосное сжатие; б - сдвиг На следующем шаге необходимо определить оси ортотропии и величины механических характеристик, для чего компоненты тензора упругих констант выписываются относительно трех углов поворота [7], по которым минимизируется целевая функция: (2) Подробно методика описана в работе [8]. В результате расчетов для каждого репрезентативного элемента можно определить оси ортотропии и механические характеристики в этих осях. Полученные результаты представляют собой тензорное поле над объемом объекта. Для практических расчетов удобно использовать аппроксимации полученных механических свойств. Результаты и обсуждение Приведенная методика была реализована для определения характеристик диафиза бедренной кости крысы. Сканирование образца выполнено с применением микро-/нанофокусной системы рентгеновского контроля для компьютерной томографии и 2D-инспекции 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%). При анализе модуля Юнга по оси Y в поперечном направлении (см. рис. 3, в) - 1595,8 ± 6,68 МПа (0,41%), в продольном (см. рис. 3, г) - 1595,8 ± 6,93 МПа (0,43%). При анализе модуля Юнга по оси Z в поперечном направлении (см. рис. 3, д) величина составила 12254,3 ± 34,16 МПа (0,27%), в продольном (см. рис. 3, е) - 12254,3 ± 19,34 МПа (0,15%). а б в Рис. 2. Компьютерная томография костного органа (а), дискретизация на репрезентативные элементы (б) и полученные направления ортотропии (в) а б в г д е Рис. 3. Распределение средних значений модуля Юнга, МПа, по слоям в определенных направлениях: а - по оси Х в продольном направлении; б - по оси Х в продольном направлении; в - по оси Y в продольном направлении; г - по оси Y в поперечном направлении; д - по оси Z в продольном направлении; е - по оси Z в поперечном направлении Сравнивая среднее значение модулей Юнга по всем слоям в направлениях X и Y, получаем, что абсолютная разница между ними не превышает 0,098%, а значит, можно заключить, что образец обладает трансверсальной изотропией. Аналогично проделав описанные операции для модуля сдвига и коэффициента Пуассона, можно показать, что упругие константы не зависят от пространственной координаты. В результате расчетов было получено, что отклонение от среднего не превышает 0,33%. Таким образом, модуль сдвига в плоскости XY составил 24640 ± 41 МПа, а в двух других направлениях - 1726 ± 3 МПа. Осредненные значения упругих констант приведены в таблице. Осредненные значения упругих констант Упругая константа Индекс компоненты 11 22 33 12 23 13 E, МПа 1590 1590 12 205 - G, МПа - 24542 1718 1718 ν 0,87 0,021 0,021 Для наглядности на рис. 4 приведены поверхность тензора упругих констант и ее сечения плоскостями. а б в г Рис. 4. Тензорная поверхность (а) и ее сечения: б - в плоскости XY; в - в плоскости YZ; г - в плоскости XZ Заключение 1. В работе представлена методика для определения распределения упругих констант образца по данным томографии. В используемой методике изучаемый объект предлагается разбить на набор репрезентативных объемов, на основе данных томографии для каждого объема производится построение конечно-элементной модели и проводятся численные эксперименты. В результате можно получить распределение упругих констант в объекте. В дальнейшем для оценки получаемого распределения предлагается использовать методы статистического анализа. 2. В работе рассмотрен диафиз бедренной кости крысы, для описания которого использована данная методика. В результате расчетов было показано, что образец обладает трансверсальной изотропией, не зависящей от пространственной координаты; были определены средние величины упругих констант и их стандартные отклонения. 3. Результаты, полученные при применении описанной методики, могут быть улучшены путем увеличения количества численных экспериментов и регуляризаций задачи определения компонент тензора упругих констант. 4. Описанный подход может быть расширен не только на объекты из костной ткани, но в этом случае необходима некоторая априорная информация о структурных составляющих.

About the authors

N. V Kharin

O. V Gerasimov

P. V Bolshakov

A. A Khabibullin

A. O Fedyanin

M. E Baltin

T. V Baltina

O. A Sachenkov

References

Statistics

Views

Abstract - 117

PDF (Russian) - 60

Refbacks

• There are currently no refbacks.