CC BY
5
Герасимов О.В., Рахматулин Р.Р., Балтина Т.В., Саченков О.А. Определение эффективных механических параметров на основе данных цифрового прототипа. Российский журнал биомеханики, 2023, Т. 27, № 3, С. 53-66 DOI: 10.15593/RZhBiomech/2023.3.04
РОССИИСКИИ ЖУРНАЛ БИОМЕХАНИКИ № 3,2023
RUSSIAN JOURNAL OF BIOMECHANICS
https ://ered.pstu. ru/index.php/rjb
Научная статья
DOI: 10.15593/RZhBiomech/2023.3.04 УДК: 531/534: [57+61]
ОПРЕДЕЛЕНИЕ МЕХАНИЧЕСКИХ СВОЙСТВ КОСТНОЙ ТКАНИ ЧИСЛЕННО-ЦИФРОВЫМ МЕТОДОМ НА ОСНОВЕ ДАННЫХ КОМПЬЮТЕРНОЙ ТОМОГРАФИИ
О.В. Герасимов, Р.Р. Рахматулин, Т.В. Балтина, О.А. Саченков
Казанский федеральный университет, Казань, Россия
О СТАТЬЕ
Получена: 10 июля 2023 Одобрена: 30 августа 2023 Принята к публикации: 31 августа 2023
Ключевые слова:
негомогенные среды, методы неразрушающего контроля, численное моделирование, компьютерная томография, пористые структуры, костная ткань
АННОТАЦИЯ
Моделирование элементов негомогенных сред является актуальным направлением при изучении композитных материалов и объектов биомеханики. Необходимость проведения исследований определяется недостаточной разработанностью методик описания механических параметров с учётом анизотропии материалов. В этом случае одним из направлений выступает применение цифрового прототипирования и методов численного моделирования. В области биомеханики применимость подходов к численному моделированию на основе данных с изображений коррелирует с клинической практикой, что обусловливается особенностями методов неразрушающего контроля и позволяет получать результаты, оказывающие влияние на качество проводимого лечения.
В работе рассматривается метод восстановления механических свойств костной ткани на основе вычислительных и натурных экспериментов. Моделирование базировалось на методе конечных элементов и предполагало применение данных цифрового прототипа, полученного проведением компьютерной томографии. Восстановление расчётной области выполнялось фильтрацией исходной сетки относительно доли содержания материала. Представленный метод позволяет восстанавливать значения модуля упругости и касательных модулей ткани на основе линейности тензора упругих констант относительно модуля Юнга.
В ходе работы было проведено моделирование костных органов свиней. Содержание животных и экспериментальные процедуры осуществлялись с соблюдением биоэтических норм. Физические испытания и численное моделирование выполнялось в условиях действия изгиба. Валидация полученного распределения значений напряжённо-деформированного состояния выполнялась по данным натурных экспериментов. Представленное решение показало слабую зависимость результатов (менее 1 %) относительно степени дискретизации расчётной области, что позволило проводить вычисления с применением меньшего количества конечных элементов, избегая применения ресурсоёмких методов восстановления геометрии образца. Полученные в рамках работы результаты соответствуют литературным данным и определяют истинные значения упругих констант костной ткани, применение которых в расчётах позволяет проводить оценку прочности элементов негомогенной структуры в условиях действия внешних нагрузок с учётом пространственного распределения материала.
©ПНИПУ
© Герасимов Олег Владимирович - научный сотрудник, e-mail: valeriv.karasikov@qmail.com ID: 0000-0002-8297-8437
© Рахматулин Рамиль Русланович - студент бакалавриата, e-mail: rrahmatulin852@qmail.com ¡D: 0009-0004-4984-9344
© Балтина Татьяна Валерьевна - доцент, к.б.н., доцент, e-mail: tvbaltina@gmail.com : 0000-0003-3798-7665 © Саченков Оскар Александрович - доцент, к.ф.-м.н., заведующий кафедрой, e-mail: 4works@bk.ru ¡D: 0000-0002-8554-2938
Эта статья доступна в соответствии с условиями лицензии Creative Commons Attribution-NonCommercial 4.0 International License (CC BY-NC 4.0)
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License (CC BY-NC 4.0)
Введение
Одним из наиболее актуальных направлений в современной персонализированной медицине выступает исследование механических свойств элементов скелета человека, подверженных влиянию внешних нагрузок или находящихся в состоянии реабилитации после различного рода травм или прямого хирургического вмешательства. Предметом для изучения в этом случае оказываются костные органы, материал которых представляет собой пористую негомогенную среду. Разработка методик описания подобного рода структур предполагает необходимость учёта уникальных свойств анизотропии ткани. Более того, каждая кость обладает определёнными геометрическими особенностями, накладывающими ограничения на способ восстановления расчётной области.
Изучение биологических материалов, таких как костная ткань, требует разработки подходов, относящихся к методам неразрушающего контроля, так как предполагает дальнейшее применение в клинической практике. Результаты, получаемые на основе таких подходов, должны быть эффективны в определении подходящего курса лечения пациента ещё на этапе диагностирования. В этом случае область проводимых исследований устанавливается
физиологическими особенностями органов: костная ткань представляет собой живой биологический материал, обладающий функцией реадаптации -перестройки трабекулярной структуры согласно возникающим внешним усилиям [3, 37]. Таким образом, применение методов описания механических свойств, в том числе оценка прочности и расчёт напряжённо-деформированного состояния, может быть направлено на изучение влияния факторов снижения двигательной активности человека (гипокинезия) [39], условий отсутствия массовых сил (гипогравитация) [17], а также проведение хирургических операций, предполагающих структурные изменения костного органа (эндопротезирование, имплантирование) [1, 6].
Среди подходов к исследованию механических свойств материала наиболее распространённым выступает совместное применение численных алгоритмов и данных цифрового прототипа [26, 34, 30, 14], построение которого может быть проведено на основе методов восстановления изображений исходной области. Данные цифрового прототипа позволяют интерпретировать усреднённые значения в элементах дискретизации снимка (пикселях или вокселях для двухмерного и трёхмерного случаев, соответственно) в качестве механических параметров среды. Существуют работы, посвящённые экспериментальному изучению связи между данными с изображений и свойствами материала [32]. Одним из наиболее распространённых способов получения снимка объекта применительно к
объектам персонализированной медицины выступает проведение рентгеновской компьютерной томографии [25, 35, 10]. В трёхмерном случае полученное таким образом изображение представляет собой целочисленный массив данных, каждое значение которого соответствует усреднённой проницаемости среды в текущем микроэлементе объёма. Как правило, данные значения проходят процедуру нормировки согласно количественной шкале рентгеновской плотности Хаунсфилда и могут быть пересчитаны в механические параметры среды [13].
В общем случае применение данных цифрового прототипа в численных расчётах предполагает построение методики учёта негомогенных свойств материала на основе их пространственного распределения. В этом случае одним из подходов выступает метод гомогенизации исходной среды [7]. Подход основан на сведении некоторых начальных анизотропных свойств к ортотропии на основе одного из методов усреднения. В первую очередь к таким методам можно отнести построение распределения значений средней длины перехвата [4, 29] с последующей аппроксимацией одним из известных способов. Данная процедура позволяет провести вычисление компонент тензора структуры, приближённо описывающего направление осей ортотропии. Существующие физические соотношения позволяют использовать его в расчётах [28]. Недостатком такого подхода выступают как необходимость определения констант на основе проведения натурных испытаний, так и предварительное вычисление тензора структуры для каждого представительного элемента. Более того, степень дискретизации исходной области оказывает значительное влияние на сходимость полученных результатов. Иным способом гомогенизации исходной области выступает метод представительного объёмного элемента, который требует проведения серии вычислительных экспериментов для определения констант тензора упругих постоянных. Существенным недостатком данного подхода выступает необходимость в моделировании каждой элементарной единицы изображения отдельным конечным элементом, на основе набора которых выполняется построение представительного элемента для гомогенизации. Аналогично способу усреднения свойств методом построения тензора структуры, подход оказывается в значительной степени зависимым от параметров расчётной сетки. Таким образом, на текущее состояние отсутствует достаточная разработанность методов решения задач, направленных на описание механических свойств элементов негомогенных сред.
Отдельно можно выделить проблему восстановления расчётной области на основе данных цифрового прототипа. Существуют методы, в основе
которых лежит учёт структурного распределения материала в пространстве [38, 18, 27]. Одни из них предполагают сглаженную методику восстановления геометрии образца [15], другие - приведение области к ступенчатому виду [33]. Методы отличаются классом расчётных моделей, которые могут быть применены к построенным на их основе сеткам. Преимущественно к данным компьютерной томографии в большинстве случаев вычисления предполагают применение ортогональных сеток, соответствующих пространству координат изображения. В этом случае устанавливается относительно простая связь между массивом значений и координатами узлов численной модели. В основе таких подходов, основанных на применении данных цифрового прототипа, выполняется предварительная обработка исходного изображения, позволяющая выделять как внешнюю геометрию образца, так и его внутреннюю структурную составляющую,
характеризуемую наличием нескольких сред.
Целью работы выступает реализация метода для определения упругих констант элементов негомогенных пористых сред на основе данных компьютерной томографии. За основу был взят метод прямого учёта структурной анизотропии материала при построении численной модели на основе цифрового прототипа объекта [36, 21, 22]. Предложенный в данной работе алгоритм позволяет определять значения механических параметров на основе постпроцессорных вычислений, полученных по результатам вычислительного и натурного экспериментов.
Материалы и методы
Построение численной модели
Численное моделирование проводилось на основе метода конечных элементов. Построение сетки выполнялось с использованием трёхмерного восьми-узлового изопараметрического конечного элемента сплошной среды с линейной аппроксимацией. Функции интерполяции перемещений имели следующий вид:
{§} = ]
И)
(1)
где {9} - вектор координат точки в локальной системе (£, п, С).
Функции формы АХ{9}) представимы как
Nс) = 1 (1+£• ^)(1л,-)(1+С- С). (2)
Предполагается, что сам материал матрикса сплошной и изотропный, а структурная анизотропия возникает в силу его неравномерного распределения в
пространстве [28]. Таким образом, компоненты матрицы упругости в каждом микрообъёме постоянны, и их значения определяются материалом соответствующего вокселя. Данная зависимость компонент матрицы упругости от вектора пространственных координат {г} может быть выражена следующим образом:
|_ В(|г})]=[В]• Ц|г}}, (3)
где [П] - матрица упругих постоянных для упругого материала среды, а ю({г}) - функция весовых коэффициентов, значения которой соответствуют обработанным данным с изображений.
Применение данной функции по отношению к пористым средам предполагает выделение двух фаз материала: упругого (костная ткань) и фонового (вещество в порах). Таким образом, принимая матрицу упругости в качестве функции координат пространства, на основе известных соотношений для определения локальной матрицы жёсткости [5] возможно произвести учёт анизотропных свойств материала по данным о распределении механических свойств ткани при построении численной модели:
[ ^ Н. [5 ({0})]г •
•[ В][ В ({0})]| 3 ({0})Ю( Г ({9})) &,
(4)
где Vе - объём конечного элемента; [5({9})] - матрица дифференцирования, связывающая деформации и перемещения; |./({9})| - якобиан преобразования координат.
Так как исходные данные компьютерной томографии представляют собой регулярную ортогональную сетку вокселей, то численное интегрирование может быть проведено при помощи метода средних прямоугольников, где в качестве узловых точек принимались геометрические центры микрообъёмов. Более подробное описание основных соотношений метода представлено в работах [22, 20, 19].
Постпроцессорный анализ
На основе полученных из численного решения значений компонент узловых перемещений может быть проведён расчёт напряжённо-деформированного состояния модели. В этом случае аналогичным образом может быть введена аппроксимирующая величина, соответствующая компонентам напряжений или деформаций:
У
({0}) = 2> ({0})
У,
(5)
где - узловые значения компонент напряжений или деформаций.
Локальное усреднение данной величины может быть выполнено на основе интегрирования по дискретному
объёму конечного элемента Vе:
^} = v{f ({е})}dv.
(6)
Оценка достоверности полученного распределения напряжений и деформаций может быть получена на основе локальной ошибки энергии по напряжениям:
1
He =-
2l Ve ^ [ D (И)Г ^ dV , (7)
где {До} - вектор локально усреднённых значений ошибки напряжений в узлах конечного элемента [22].
Данная величина может быть нормирована по энергии деформации
ие = [ {а}" {ё} ^ (8)
согласно следующей формуле:
He =
He
Ue + He
100% .
(9)
Оценка полученных из численных расчётов значений нормированной ошибки энергии позволяет определять области, соответствующие допустимой погрешности вычислений, что может быть применено для проведения постпроцессорного анализа результатов моделирования: применение распределения значений ошибки энергии позволяет выполнить построение эффективного алгоритма определения зон костного органа, подверженных разрушению под воздействием внешних нагружающих факторов [22, 12].
Определение эффективных свойств
Предположение о линейности и изотропности моделируемого материала позволяет установить линейную зависимость между компонентами матрицы упругих постоянных [Б] и модулем упругости Юнга Е:
[ D (E )]= *[Dy
E .
(10)
Более того, интеграл для вычисления локальной матрицы жёсткости также предполагает линейность относительно матрицы упругих постоянных:
K ([D])
■[ D].
(11)
Принимая во внимание особенности алгоритма сборки глобальной матрицы жёсткости [Кт] на основе компонент соответствующих локальных матриц жёсткости [Ке;] (где / - номер конечного элемента) из [5], получаем
Km ([D])
[D].
(12)
Таким образом, с учётом соотношений (10)-(12) можно сделать вывод о наличии линейной зависимости
между результатами, полученными на основе данных вычислительных и натурных экспериментов. В этом случае построение диаграмм нагружения образца позволяет восстанавливать соотношение между значениями возникающих во время испытания усилий
при соответствующих перемещениях:
—
-¡Г, (13)
к = -
где Fact - значение усилия, полученное из натурного эксперимента, а Feq - величина эквивалентной силы, соответствующей численному нагружению образца.
Таким образом, пропорциональность между истинным модулем упругости Юнга Ef и начальным Einit, устанавливаемым на этапе построения численной модели, может быть выражена в виде
Eeff = к ■ Einit .
(14)
Коэффициент пропорциональности k может быть определён на основе угла наклона линейных участков диаграмм нагружения образца:
tan (p¿)
к =
tan (а)
(15)
где р,- - угол наклона /-го линейного участка на диаграмме, соответствующей натурному эксперименту, а а - угол наклона графика, полученного на основе результатов численного моделирования.
Необходимо отметить, что, несмотря на рассмотрение модели деформирования образца согласно линейному закону, в силу наличия неоднородного распределения материала в объёме возникает его локальное разрушение, которое обусловливает возможность выделения нескольких линейных участков на основе диаграммы физического разрушения образца:
Рг : Ф- (16)
В этом случае предполагается восстановление эффективных значений как модуля упругости (первый линейный участок), так и касательных модулей, соответствующих упрочнению (остальные диапазоны линейной зависимости).
Величина реакции на эквивалентную силу может быть вычислена на основе результатов численного моделирования согласно следующей формуле:
Feq = J ^n H{n} dA ,
(17)
где {п} - вектор, устанавливающий направление воздействия индентора; Ап — проекция нагруженной площади на перпендикулярную к вектору {п} плоскость, а [о] - тензор напряжений, усреднённый по поверхности кинематически нагруженных в области диафиза конечных элементов.
Применение описанного алгоритма по отношению к элементам пористой костной структуры позволяет восстанавливать эффективные значения модуля упругости Юнга, соответствующие истинным механическим параметрам твёрдой костной ткани.
Проведение испытаний
Испытания проводились на образцах трубчатых костей конечностей вьетнамских свиней-самцов весом 15-20 кг. Животных помещали в индивидуальные клетки в стандартных лабораторных условиях с неограниченным доступом к пище и воде и 12-часовым циклом день/ночь. Животных включали в эксперименты после периода адаптации продолжительностью не менее 7 дней. В экспериментах использовались кости животных после контузионной травмы позвоночника. Животных умерщвляли на 8-й неделе после травмы спинного мозга. Эвтаназия проводилась под глубокой анестезией путём передозировки ингаляционного наркозного средства (изофлурана). Извлекались кости передних и задних конечностей.
Протокол экспериментов, включающий анестезию, хирургическое вмешательство, послеоперационный уход, тестирование и эвтаназию, был одобрен Комитетом по уходу за животными Казанского государственного медицинского университета (протокол № 5 от 20 мая 2020 года). Все экспериментальные процедуры проводились в соответствии со стандартами, чтобы свести к минимуму страдания животных и размер экспериментальных групп.
Численное моделирование проводилось на основе данных, полученных путём проведения рентгеновской компьютерной томографии. Сканирование образцов выполнялось с применением микро-/нанофокусной системы рентгеновского контроля для компьютерной томографии и 20-инспекции Phoenix V\tome\X S240 в лаборатории рентгеновской компьютерной томографии Института геологии и нефтегазовых технологий Казанского (Приволжского) федерального университета. Система оснащена двумя рентгеновскими трубками: микрофокусной с максимальным ускоряющим
напряжением 240 кВ мощностью 320 Вт и нанофокусной с максимальным ускоряющим напряжением 180 кВ мощностью 15 Вт. Для первичной обработки данных и создания объёмной (воксельной) модели образца на базе рентгеновских снимков (проекций) использовалось программное обеспечение datos\x reconstruction. Зафиксированный в держателе образец помещался на вращающийся столик камеры компьютерного томографа на оптимальном расстоянии от источника рентгеновского излучения. Съёмка проводилась при ускоряющем напряжении 90-100 кВ и токе 140-150 мА.
Вычислительные и натурные эксперименты проводились на образце предплечевой кости. Объём данных компьютерной томографии - 752*752*752 вокселя, физический размер которых составлял 0,2*0,2*0,2 мм. Таким образом, габариты области сканирования - 150,4*150,4*150,4 мм. На рис. 1 представлена визуализация обработанных данных компьютерной томографии образца. Бинаризация исходного массива значений проводилась на основе метода Оцу, позволяющего провести разделение упругого костного материала и вещества в порах.
К исходному образцу прикладывались кинематические граничные условия, соответствующие испытанию на изгиб (см. рис. 2). Дистальные области подвергались жёсткому защемлению в зонах контакта с
поверхностью костного соприкосновения A, B, C, и D):
органа:
{и (X, у, z)} = 0: Vx, у, z е Sd
Sdof = SA ^>SB ^>SC ^>SD ■
границы
(18)
К верхней области диафиза прикладывалось перемещение и0 в поперечном направлении: участок соприкосновения ЕО) с помощью индентора заданной геометрии (рис. 3, синий круг: чёрным маркером отмечена середина контактного взаимодействия).
{и (X, у, z)J = {0, —и0, 0}Г : Vx, у, z е S
disp ,
Sdisp = S
(19)
FG ■
Рис. 1. Обработанное изображение предплечевой кости
Рис. 2. Схема кинематического нагружения образца костного органа
Рис. 3. Проведение натурного эксперимента на трёхточечный изгиб
На свободной поверхности костного органа статические граничные условия соответствовали нулевым:
{q} = 0: Vx,у,z е Sfree,
< \ (20)
Sfree ~SV disp ) -
где dV — внешняя поверхность геометрии образца.
Физические испытания проводились на разрывной машине УТС 110М-100. Проведение натурного эксперимента предполагало увеличение
прикладываемого усилия до разрушения образца в поперечном направлении (перелом костного органа с образованием трещины). Хрупкое разрушение материала, как правило, соответствовало области приложения кинематического нагружения (рис. 3: контакт в синем круге).
Восстановление расчётных сеток проводилось на основе различной степени дискретизации исходного изображения. Моделирование предполагало применение 20, 40 и 80 конечных элементов на стороне области (рис. 4). Фильтрация исходной регулярной ортогональной сетки проводилась путём удаления конечных элементов с относительным содержанием костной ткани менее 5 %:
V0 = ,1 е[1; N ], (21)
где V0 — восстановленная расчётная область; V, — подобъём, соответствующий /-ому конечному элементу, и N — их количество.
Область каждого конечного элемента может быть определена в следующем виде:
VI : х е[х,;х, + а],у е[у,;у + Ъ], 2 е[2;2 +с],
(22)
где {х,, у,, х,} соответствуют координатам узла с наименьшей нормой радиуса-вектора в глобальной системе координат, а {а, Ь, с} - геометрические размеры конечных элементов регулярной ортогональной сетки в направлении трёх декартовых координатных осей.
Граничные условия прикладывались к полученным таким образом расчётным сеткам в соответствии с характером нагружения образца в физическом эксперименте. Дальнейшие результаты будут представлены для ансамбля, построенного на основе сетки с 80 конечными элементами на стороне области.
Исходный модуль упругости Юнга костной ткани принимался равным 1 ГПа, коэффициент Пуассона - 0.3. Восстановление эффективных механических параметров
Рис. 4. Ортогональные конечно-элементные сетки различной степени дискретизации. Количество конечных элементов на стороне: а - 20; б - 40; в - 80
материала на основе представленного алгоритма предполагало дальнейший пересчёт начально установленного значения модуля упругости и, соответственно, касательных модулей, исходя из результатов, представленных в виде диаграмм нагружения образца. В силу линейной зависимости результатов вычислений относительно прикладываемого перемещения рассматривалось кинематическое нагружение, соответствующее 1 мм. Значения, полученные для перемещений на следующих этапах нагружения образца определялись в соответствии с диапазоном, не превосходящим разрушение костного органа во время проведения натурного эксперимента.
Общее время вычислений для каждой расчётной сетки составило около 15 минут.
Оценка напряжённо-деформированного состояния
На этапе постпроцессорного анализа вычислялось напряжённо-деформированное состояние в узлах сетки. На основе полученных результатов, с учётом неравномерного распределения материала, выполнялось усреднение значений в объёме каждого конечного элемента по локальным данным компьютерной томографии. Результаты вычислений приводились к главным компонентам тензора напряжений, на основе
а
б
в
которых выполнялось вычисление значений по Мизесу. Полученное таким образом распределение позволяло учитывать относительное содержание костной ткани в объёме. Достоверность усреднённого напряжённо -деформированного состояния оценивалось на основе вычисления нормированной ошибки энергии по напряжениям. В этом случае исследовалось влияние доли содержания упругого материала в объёме конечного элемента в совокупности с прикладываемыми граничными условиями. Анализ результатов осуществлялся на основе двух типов сеток: расчётной, соответствующей регулярной ортогональной конечно-элементной сетке, и постпроцессорной, представляющей собой интерполяцию полученных результатов на геометрию исходного образца, восстановленную по данным о распределении плотности материала. Определение областей, наиболее подверженных разрушению в условиях действия внешних нагрузок, проводилось на основе выбора подобъёмов расчётной дискретизации области, соответствующих наименьшей ошибке результатов вычислений: конечные элементы с напряжённым состоянием, превышающим предел
прочности костного материала, отмечались как критические.
Результаты
На рис. 5, а представлены результаты постпроцессорного анализа для сетки с применением 80 конечных элементов на стороне для исходной дискретизации области: получено распределение значений нормированной ошибки энергии по напряжениям. Красным цветом отмечены области, соответствующие наибольшим значениям (в процентах). Следует отметить большие величины ошибки (более 80 %) в областях, граничащих с поверхностью костного органа, а также в зонах кинематического нагружения модели (на дистальных участках). В областях, интересующих нас с точки зрения расчётов (верхняя граница диафиза), нормированная ошибка не превышала 50 %. В других зонах значения оказывалась не более 20 %.
Рис. 5, б отображает распределение напряжений по Мизесу. Полученные значения соответствуют прикладываемому перемещению величиной 1 мм.
б
Рис. 6. Результаты вычислений, перенесенные на цифровой двойник: а - первая главная компонента тензора
напряжений, МПа; б - первое главное направление
Наибольшие значения напряжений по Мизесу отмечаются в области приложения кинематического нагружения (зона контактного взаимодействия поверхности костного органа с индентором в области диафиза) и соответствуют 83 МПа, что также характеризуется величиной нормированной ошибки энергии в данной области, находящейся в диапазоне 2040 % (рис. 5, а). Значения напряжений по Мизесу в области жёсткого защемления не превышали 40 МПа и соответствовали нормированной ошибке энергии величиной 20-55 %.
Определение критических областей, подверженных локальному разрушению образца, осуществлялось на основе вычисления главных компонент тензора напряжений. В этом случае анализу подвергались первая и третья главные компоненты (рис. 6, а и рис. 7, а) с их соответствующими направлениями (см. рис. 6, б и рис. 7, б). Наибольшие абсолютные значения первой главной компоненты тензора напряжений соответствовали области ограничения в перемещениях на левом дистальном участке образца и составляли 40 МПа (рис. 6, а), второй главной компоненты - области
кинематического нагружения на верхней части диафиза костного органа и равнялись 78 МПа (рис. 7, а).
Наибольшие отрицательные значения для обеих главных компонент соответствовали области контакта с индентором. Полученные направления главных компонент (рис. 6, б и рис. 7, б) определяли линии тока образующихся напряжений и соответствовали характеру возникновения поперечной эквивалентной силы реакции.
На рис. 8 представлены результаты проведения вычислительного и натурного экспериментов. Данные, соответствующие физическому испытанию, позволяют выделить два линейных участка, на основе которых могут быть восстановлены истинные значения механических параметров материала (14), (15) - модуля упругости и касательного модуля соответственно. Значения коэффициентов детерминации R2 для аппроксимаций линейных участков приведены на рис. 8. Диаграмма численного моделирования, отмеченная оранжевым цветом, получена путём экстраполяции решения по значениям вычисленных усилий для перемещения величиной 0 и 1 мм соответственно.
а
б
Рис. 7. Результаты вычислений, перенесенные на цифровой двойник: а - вторая главная компонента тензора
напряжений, МПа; б - второе главное направление.
Интерполяция выделенных линейных участков, соответствующих диаграмме нагружения физического образца, проводилась при значениях коэффициентов детерминации, близких к единице (см. рис. 8: тонкие жёлтая и зелёная линии, Я2 > 0,998). Таким образом, на основе полученных диапазонов линейной зависимости были вычислены коэффициенты пропорциональности для определения эффективных значений параметров материала (13)-(15). Вариация соответствующих значений модуля упругости и касательного модуля в зависимости от степени дискретизации расчётной области представлена в таблице.
Величина стандартного отклонения для значений, как модуля упругости, так и касательного модуля составила приближённо 0,0125 ГПа, соответствующий коэффициент вариации (стандартное отклонение, отнесённое к среднему) - 0,706 %.
Обсуждение
Достоверность значений, полученных при определении напряжённо-деформированного состояния, оценивалась на основе вычисления нормированной ошибки энергии. В этом случае наибольший интерес
представляли конечные элементы с максимальными усреднёнными по объёму значениями напряжений по Мизесу. Корректная интерпретация и соответствующий анализ полученных данных предполагали следующее понимание процедуры локального усреднения на основе данных компьютерной томографии объекта: область интегрирования, соответствующая объёму конечного элемента, рассматривалась в качестве подобъёма с неизвестным пространственным распределением усредняемой величины, но с определённым средним значением. Таким образом, в случае исследования напряжений по Мизесу возникновение ложно высоких или низких значений объяснялось неравномерным объёмным распределением материала, что приводило к необходимости определения ошибки вычислений для каждого конечного элемента сетки.
Представленные результаты показали, что наибольшие значения нормированной ошибки энергии по напряжениям соответствуют конечным элементам на границе, что возникает в силу уменьшения жёсткости, определяемой компонентами локальной матрицы жёсткости. В интересующих с точки зрения расчётов областях ошибка энергии не превышала 50 %, что позво-
а
1.0 1.5 2.0 Перемещение, мм Рис. 8. Диаграмма нагружения образца
ляет получать достоверную оценку результатов, полученных на этапе постпроцессорного анализа.
Сопоставление результатов численного
моделирования с данными натурного эксперимента на основе вычисления реакции на эквивалентную силу (рис. 8) позволило провести расчёт значений эквивалентного модуля упругости и касательного модуля для костного материала с учётом выделения двух соответствующих линейных участков. Определение таких диапазонов линейной зависимости перемещения от усилия может быть произведено в автоматическом порядке [9, 8]. Образование второго линейного участка диаграммы нагружения физического образца объясняется локальным разрушением костного материала: в этом случае происходит нарушение целостности внешнего кортикального слоя костного органа и образование трещин на поверхности приложения кинематических граничных условий в области диафиза [11, 31]. Данный эффект достоверно описывается результатами численного моделирования, определяющими возникновение максимальных (по модулю) значений напряжений по Мизесу в зоне контактного взаимодействия с индентором (рис. 5, б).
Следует отметить, что описание разрушения образца в зонах с напряжениями, превышающими предельные значения [16, 23], может быть проведено с учётом применения в расчётах модуля упругости Юнга, близкого к истинному значению для соответствующего материала. Таким образом, определение параметров материала выступало следующим этапом проведения расчётов. В этом случае полученные значения для реакций на эквивалентные силы позволяют восстанавливать эффективный модуль упругости и касательный модуль, величина которых оказывается близкой к истинному значению для «чистой» костной ткани (без пор) [24], так как предполагает
интегрирование с учётом вокселей, соответствующих упругому материалу. Результаты расчётов показали незначительную зависимость эффективных значений относительно степени дискретизации исходной сетки, что объясняется особенностями построения численной модели на основе данных цифрового прототипа и выступает в качестве следствия, определяющего преимущество данного подхода.
Значения эффективного модуля упругости Юнга, полученные на основе сопоставления результатов вычислительных и натурных экспериментов, коррелируют с данными, представленными в литературе [24]. Применение полученного таким образом модуля упругости в расчётах позволяет восстанавливать напряжённо-деформированное состояние, близкое к действительному. В этом случае подобъёмы с усреднёнными напряжениями, превышающими предел прочности для костной ткани, определяют критические области, наиболее подверженные разрушению в условиях действия приложенных нагрузок [12].
Характер разрушения образца зависит от типа костного органа и от его геометрических и механических характеристик. Описание возникающих концентраторов напряжений может быть проведено на основе какой-либо главной компоненты тензора напряжений. В работе [22] было показано, что наибольшие значения первой главной компоненты соответствуют области образования трещины, наблюдаемой в натурном эксперименте. Такой характер разрушения проявляется при испытании трубчатых костей животных [2]. В данной работе испытание образца со провождалось локальным разрушением кортикального слоя в области
приложения индентора. При такой модели образования трещин эффективным способом оказывается перестроение конечно-элементной сетки в области, близкой к образованию смятия костной ткани, и
Значения эффективных механических параметров
Механический параметр Кол-во конечных элементов
20 40 80
Модуль упругости, ГПа 1,77 ± 0,01 1,78± 0,01 1,75± 0,01
Касательный модуль, ГПа 1,27± 0,01 1,29± 0,01 1,26± 0,01
соответствующий пересчёт матрицы жёсткости конечно -элементного ансамбля с учётом деформированного состояния. Данная процедура позволяет учитывать нелинейность деформирования образца, однако в значительной степени усложняет аппарат взвешенного интегрирования по данным компьютерной томографии области.
Погрешность расчётных значений, полученных на основе вычислений с применением данных компьютерной томографии, варьируется в пределах 1 % относительно степени дискретизации расчётных сеток и может быть обусловлена как неточностью приложения кинематических граничных условий, соответствующих характеру нагружения физического образца в натурном эксперименте, так и влиянием сходимости численного метода с точки зрения способа интегрирования локальной матрицы жёсткости каждого конечного элемента и в силу образованного количества степеней свободы модели. Таким образом, предложенный метод прямого учёта данных о пространственном распределении плотности материала позволяет восстанавливать эффективные механические свойства материала на основе меньшей степени дискретизации исходного изображения, что определяет оптимальность времени вычислений относительно размера получаемой
Список литературы
1. Ананьева А.Ш., Бараева Л.М., Быков И.М., Веревкина Ю.В., Курзанов А.Н. Моделирование повреждений костных структур в экспериментах на животных // Инновационная медицина Кубани. - 2021. -№ 1. - С. 47-55. DOI: 10.35401/2500-0268-2021-21-1-4755.
2. Ахметзянова А.И., Шарафутдинова К.Р., Сабирова Д.Э., Балтин М.Э., Герасимов О.В.,Балтина Т.В., Саченков О.А. Оценка влияния тяжести травмы спинного мозга на механические свойства костей задних конечностей опытных крыс // Российский журнал биомеханики. -2022. - № 4. - С. 45-55. DOI: 10.15593/RzhBiomeh/2022.4.04.
3. Киченко А.А., Тверье В.М., Няшин Ю.И., Осипенко М.А., Лохов В.А. Постановка начально-краевой задачи о перестройке трабекулярной костной ткани // Российский журнал биомеханики. - 2012. - Т. 16, № 4 (58). - С. 36-52.
системы уравнений и исключает необходимость в применении сложных алгоритмов восстановления геометрии образцов.
Заключение
В работе предложен метод восстановления значений модуля упругости и касательных модулей на основе совместного применения результатов вычислительных и натурных экспериментов. Методика предполагает применение численного моделирования методом конечных элементов, основанного на взвешенном интегрировании матриц жёсткости конечных элементов сетки по данным компьютерной томографии образцов. Полученное таким образом решение в перемещениях позволяет восстанавливать напряжённо-
деформированное состояние, достоверность которого может быть определена распределением значений нормированной ошибки энергии по напряжениям. Валидация результатов численного моделирования по данным компьютерной томографии выполнялась на основе проведения физических испытаний.
Представленные результаты показали применимость метода взвешенного интегрирования для расчёта элементов пористой структуры. Использованный в работе способ восстановления расчётной области позволил проводить вычисления на сетках с небольшим количеством элементов. Погрешность результатов вычислений определялось неточностью приложения граничных условий, соответствующих натурному эксперименту. Полученные в рамках предложенного подхода эффективные значения механических параметров соответствуют истинным значениям для костной ткани и могут быть использованы для определения критических областей костных органов, подверженных разрушению в условиях действия внешних нагрузок.
4. Киченко А.А., Тверье В.М., Няшин Ю.И., Заборских А.А. Экспериментальное определение тензора структуры трабекулярной костной ткани // Российский журнал биомеханики. - 2011. - Т.15, № 4. - С.78-93.
5. Крылов О.В. Метод конечных элементов и его применение в инженерных расчётах // Учеб. пособие для вузов. - М.: Радио и связь. - 2002. - 104 с. - ISBN 5-256-01627-X.
6. Маслов Л.Б., Дмитрюк А.Ю., Жмайло М.А., Коваленко А.Н. Исследование прочности эндопротеза тазобедренного сустава из полимерного материала // Российский журнал биомеханики. - 2022. -№ 4. - С. 19-33. DOI: 10.15593/RzhBiomeh/2022.4.02.
7. Саченков О.А., Герасимов О.В., Королёва Е.В., Мухин Д.А., Яикова В.В., Ахтямов И.Ф., Шакирова Ф.В., Коробейникова Д.А, Чжи К.К. Построение неоднородной конечно-элементной модели по данным компьютерной
томографии // Россиийский журнал биомеханики. -2018. - Т. 22, № 3. - С. 290-302. DOI: 10.15593/RJBiomeh/2018.3.05.
8. Саченков О.А., Яикова В.В., Харин Н.В. Программа для автоматического определения механических свойств костной ткани по экспериментальным данным. Свидетельство о регистрации программы для ЭВМ, № 2020615249, 19.05.2020. Заявка № 2020613959 от 03.04.2020.
9. Яикова В.В., Герасимов О.В., Харин Н.В, Балтина Т.В., Саченков О.А. Автоматическое определение механических свойств костной ткани по экспериментальным данным // Информационные технологии и нанотехнологии. - 2021. - Т. 3. - С. 034353.
10. Bagirov, A., Suvarly, P., Ogaryov, E., Yeltsin, A., Mininkov, D., Tagizade, A. Multislice computed tomography in the complex assessment of deformities of long tubular bones of the lower extremities: prospective cohort study // N.N. Priorov Journal of Traumatology and Orthopedics. - 2023. - Vol. 29. - P. 269-277. DOI: 10.17816/vto111559.
11. Bakhmet'ev, V. Destruction of a long tubular bone in combined mechanical and thermal exposures // Sudebno-Meditsinskaia Ekspertiza. - 1990. - Vol. 33. - P. 19-20.
12. Bolshakov P.V., Sachenkov O.A. Destruction simulation for the inhomogeneous body by finite elementmethod using computed tomography data, Russian Journal of Biomechanics. - 2020. - Vol. 24, No. 2. - P. 248-258. DOI: 10.15593/RzhBiomeh/2020.2.12.
13. Choi, S., Park, Y.-K., Kim, J.H., Moon, H., Kwon, W.-K., Ham, C. Clinical Importance of Hounsfield Unit in Computed Tomography of Sub-Axial Cervical Vertebral Body // Asian Journal of Pain. - 2022. - Vol. 8, No. 2. DOI: 10.35353/ajp.2022.00059.
14. Correia Marques, F., Boaretti, D., Walle, M., Scheuren, A., Schulte, F., Müller, R. Mechanostat parameters estimated from time-lapsed in vivo micro-computed tomography data of mechanically driven bone adaptation are logarithmically dependent on loading frequency // Frontiers in Bioengineering and Biotechnology. - 2023. - Vol. 11. DOI: 10.3389/fbioe.2023.1140673.
15. Cremers, D., Rousson, M., Deriche, R. A review of statistical approaches to level set segmentation: Integrating color, texture, motion and shape // International Journal OJ Computer Vision. - 2006. - Vol. 72. - P. 195-215.
16. Crenshaw T.D., Peo Jr.E.R., Lewis A.J., Moser B.D. Bone Strength as a Trait for Assesing Mineralization in Swine: A Critical Review of Techniques Involved // J. Ani Sci. -1981. - Vol. 53, No. 3. - P. 827-835. DOI: 10.2527/JAS1981.533827X.
17. Fomina, E., Lysova, N., Chernova, M., Khustnudinova, D., Kozlovskaya, I. Comparative Analisys of Efficacy of Countermeasure Provided by Different Modes of Locomotor Training in Space Flight // Fiziologiia Cheloveka. - 2016. -Vol. 42. - P. 84-91.
18. Frey, P.J. Generation and adaptation of computational surface meshes from discrete anatomical data // Int. J. Numer. Meth. Engng. - 2004. - Vol. 60. - P. 1049-1074. DOI: 10.1002/nme.992.
19. Gerasimov O., Kharin N., Statsenko E., Mukhin D., Berezhnoi D., Sachenkov O. Patient-Specific Bone Organ Modeling Using CT Based FEM // Lecture Notes in Computational Science and Engineering. - 2022. -Vol. 141. - P. 125-139. DOI: 10.1007/978-3-030-878092 10.
20. Gerasimov O., Sharafutdinova K., Rakhmatullin R., Baltina T., Baltin M., Fedianin A. Application of a digital prototype for CT-based bone strength analysis // VIII International Conference on Information Technology and Nanotechnology. - 2022. - P. 1-6. DOI: 10.1109/ITNT55410.2022.9848693.
21. Gerasimov, O., 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. - P. 87-97. DOI: 10.15593/RJBiomech/2019.1.10.
22. Gerasimov, O., Kharin, N., Fedyanin, A., Bolshakov, P., Baltin, M., Statsenko, E., Fadeev, F., Islamov, R., Baltina, T., Sachenkov, O.A. Bone Stress-Strain State Evaluation Using CT Based FEM // Frontiers in Mechanical Engineering. -2021. - Vol. 7. DOI: 10.3389/fmech.2021.688474.
23. Imai K. Computed Tomography-Based Finite Element Analysis to Assess Fracture Risk and Osteoporosis Treatment // Wjgem. - 2015. - Vol. 5, No. 3. - P. 182-187. DOI: 10.5493/wjem.v5.i3.182.
24. Kieser, D., Kanade, S., Waddell, J., Kieser, J., Theis, J.-C., Swain, M. The deer femur - A morphological and biomechanical animal model of the human femur // BioMedical Materials and Engineering. - 2014. - Vol. 24. -P. 1693-703. DOI: 10.3233/BME-140981.
25. Li, W., Jian, Y., Zhou, X., Wang, H. In situ tensile damage characterization of C/C composites through X-ray computed tomography and digital volume correlation // Ceramics International. - 2022. DOI: 10.1016/j.ceramint.2022.11.231.
26. Luan, S., Chen, E., John, J., Gaitanaros, S. A data-driven framework for structure-property correlation in ordered and disordered cellular metamaterials // Arxiv - 2023. DOI: 10.48550/arXiv.2304.04809.
27. Madi K., Forest S., Boussuge M., Gailliégue S., Lataste E., Buffiére J.-Y., Bernard D., Jeulin D. Finite element simulations of the deformation of fused-cast refractories based on X-ray computed tomography // Computational Materials Science. - 2007. - Vol. 39, No. 1. - P. 224-229. DOI: 10.1016/j.commatsci.2006.01.033.
28. Maquer, G., Musy, S., Wandel, J., Gross, T., Zysset, P. Bone Volume Fraction and Fabric Anisotropy Are Better Determinants of Trabecular Bone Stiffness Than Other Morphological Variables // Journal of Bone and Mineral Research. - 2014. - Vol. 30. DOI: 10.1002/jbmr.2437.
29. Mathieu Simon, Michael Indermaur, Denis Schenk, Seyedmahdi Hosseinitabatabaei, Bettina M. Willie, Philippe Zysset. Fabric-elasticity relationships of tibial trabecular bone are similar in osteogenesis imperfecta and healthy individuals // Bone. - 2022. - Vol. 155. - P. 116282. DOI: 10.1016/j.bone.2021.116282.
30. Mazur, K., Krawczuk, M., Dqbrowski, L. A new finite element with variable Young's modulus // International Journal for Numerical Methods in Biomedical engineering. - 2023. DOI: 10.1002/cnm.3712.
31. Peña Fernández, M., Schwiedrzik, J., Bürki, A., Peyrin, F., Michler, J., Zysset, P., Wolfram, U. In situ synchrotron radiation ^CT indentation of cortical bone: Anisotropic crack propagation, local deformation, and fracture // Acta Biomaterialia. - 2023. - DOI: 10.1016/j.actbio.2023.04.038.
32. Rho J.Y., Hobatho M.C., Ashman R.B. Relations of Mechanical Properties to Density and CT Numbers in Human Bone. // Med. Eng. Phys. - 1995. - Vol. 17, No. 5. -P. 347-355. DOI: 10.1016/1350-4533(95)97314-f.
33. Rietbergen B., Weinans H., Huiskes R., Odgaard A. A new method to determine trabecular bone elastic properties and
loading using micromechanical finite-element models // Journal of Biomechanics. - 1995. - Vol. 28, No. 1. -P. 69-81. DOI: 10.1016/0021-9290(95)80008-5. 37.
34. Sanchez-Molina, D., Garcia-Vilana, S., Saez, L., Lluma, J. A strain rate dependent model with decreasing Young's 38. Modulus for cortical human bone // Biomedical Physics & Engineering Express. - 2023. DOI: 10.1088/2057-1976/acd458.
35. Shi, D., Zhang, B., Liu, C., Wang, L., Yang, X., Luo, R. In-situ study on compressive behaviors of different types of 3D SiC/SiC composites using X-ray computed tomography and 39. digital image correlation // Journal of Materials Research and Technology. - 2023. - Vol. 22. DOI: 10.1016/j.jmrt.2022.12.178.
36. Vorobiev O.V., Semenova E.V., Mukhin D.A., et al. The Image-Based Finite Element Evaluation of the Deformed
Финансирование. Исследование выполнено при финансовой поддержке, выделяемой Казанскому федеральному университету
по государственному заданию в сфере научной деятельности, проект № FZSM-2023-0009.
Конфликт интересов. Авторы заявляют об отсутствии конфликта интересов.
DETERMINATION OF THE BONE TISSUE MECHANICAL PROPERTIES BY A NUMERICAL-DIGITAL METHOD USING CT DATA
O.V. Gerasimov, R.R. Rakhmatulin, T.V. Baltina, O.A. Sachenkov
Kazan (Volqa region) Federal University, Kazan, Russia
ABSTRACT
Modeling of elements of inhomogeneous media is an actual direction in the study of composite materials and biomechanics objects. The need for research is determined by the insufficient development of methods for describing mechanical parameters, taking into account the anisotropy of materials. In this case, one of the directions is the use of digital prototyping and numerical modeling methods. In the field of biomechanics, the applicability of approaches to numerical modeling based on image data correlates with cases of clinical practice, which is due to the peculiarities of non-destructive testing methods and allows obtaining results that affect the quality of treatment.
The paper considers a method for restoring effective mechanical properties of a material based on computational and field experiments. The simulation was based on the finite element method and assumed the use of data from a digital prototype obtained by computed tomography. The restoration of the calculated area was performed by filtering the initial grid relative to the proportion of the content of the material. The presented method makes it possible to restore the effective values of the elastic modulus and tangential modules based on the relationship established within the framework of the approach between the mechanical properties of the material and the loading parameters of the model sample.
In the course of the work, the bone organs of pigs were modeled. Animal husbandry and experimental procedures were carried out in compliance with bioethical norms. Physical tests and numerical modeling were carried out under the conditions of bending action. Validation of the obtained distribution of values of the stress-strain state was performed according to the data of full-scale experiments. The presented solution showed a weak dependence of the results (less than 1%) relative to the degree of discretization of the computational domain, which made it possible to perform calculations using fewer finite elements, avoiding the use of resource-intensive methods to restore the geometry of the sample. The results obtained in the framework of the work correspond to the literature data and determine the true values, the use of which in calculations makes it possible to assess the strength of the elements of an inhomogeneous structure under the action of external loads, taking into account the spatial distribution of the material.
©PNRPU
State // PNRPU Mechanics Bulletin. - 2021. - No. 2. -P. 44-54. DOI: 10.15593/perm.mech/2021.2.05. Yogev, O., Antonsson, E.K. Growth and development of inhomogeneous structures // Proceedings of ICED. - 2007. Young, P., Beresford-West, T., Coward, S., Notarberardino, B., Walker, B., Abdul-Aziz, A. An efficient approach to converting three-dimensional image data into highly accurate computational model // Philosophical Transactions. Series A, Mathematical, physical, and engineering sciences. - 2008. -Vol. 366. - P. 3155-3173. DOI: 10.1098/rsta.2008.0090. Zorbas, Y.G., Deogenov, V.A., Tsiamis, C., Yerullis, A.L. Bone mineral density during prolonged hypokinesia and rehydration in healthy subjects // Trace Elements and Electrolytes. - 2008. - Vol. 25. - P. 100-110. DOI: 10.5414/TEP25100.
ARTICLE INFO
Received: 10 July 2023
Approved: 30 August 2023
Accepted for publication: 31 August 2023
Key words:
inhomogeneous media, non-destructive testing methods, numerical modeling, computed tomography, porous structures
