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

Statistics

Views

Abstract - 110

PDF (Russian) - 56

Refbacks

• There are currently no refbacks.