Full Text

Введение В настоящее время одной из актуальных проблем клинической практики выступает оценка напряженно-деформированного состояния органов при проведении хирургического вмешательства. Первоначально на основе параметрической цифровой модели органа производится локальная оценка прочности кости после установки имплантата, в том числе при различных возможных патологиях [1-3]. Необходимость такой оценки обеспечивается в том числе физиологическими процессами, возникающими в костной ткани под действием силовых факторов. Происходит так называемая перестройка костной ткани, что влечет за собой изменения структурных свойств исследуемых органов. Этот эффект неоднократно был замечен в рамках экспериментов in vivo на животных [4-6]. Существуют подходы к моделированию таких процессов. Разрабатываемые модели позволяют определить архитектуру костной ткани после перестройки на основе численных расчетов [7-9]. Необходимость учета структурных свойств органов продиктована особенностями локальной прочности. Один из подходов оценки механических свойств костной ткани основан на гомогенизации исследуемой области с применением различных методов. Одним из распространенных подходов является метод представительного объемного элемента [10, 11]. Такой способ позволяет учитывать в том числе и нелинейность материала [12]. Применение данных компьютерной томографии объекта для построения представительного объемного элемента расширило возможности ее использования [12, 13]. При этом стоит отметить трудозатратность такого метода. К другому подходу гомогенизации можно отнести тензор структуры, определяющий основные направления распределения свойств материала [14, 15]. Существующие физические соотношения позволяют использовать тензор структуры в расчетах [16]. В первую очередь получили большое распространение вопросы гомогенизации области расчета путем введения эффективных свойств [17-19]. В данном случае некоторой нетривиальной геометрии ставится в соответствие упрощенная область с усредненными параметрами, позволяющая производить расчет с приближенными характеристиками объекта. К другому способу моделирования можно отнести построение негомогенной модели на основе вычисления локальных свойств для каждого представительного объема в отдельности [20-22]. Данный подход обладает существенным недостатком в виде ресурсоемкости в силу необходимости проведения колоссального количества расчетов. Можно отметить обособленный подход к моделированию гетерогенных сред, в рамках которого формирование конечно-элементного ансамбля осуществляется на основе данных компьютерной томографии. Учет неоднородного распределения свойств при таком подходе осуществляется при формировании расчетных матриц для метода конечных объемов [23]. Отдельно были рассмотрены задачи построения расчетной сетки по данным компьютерной томографии [24-25]. Обобщение вышеуказанных методик было произведено в работе [26] для персонализированных расчетов по данным медицинских изображений. Предлагаемый подход учета данных компьютерной томографии в трехмерной постановке был реализован на основе метода конечных элементов [27]. В этой же работе приведена оценка влияния плотности заполнения вокселями конечно-элементного ансамбля на результаты расчетов. Цель данного исследования - разработка конечного элемента для формирования расчетного ансамбля и моделирования механического поведения с учетом данных двумерных медицинских изображений. 1. Материалы и методы Исследование вводит в рассмотрение методологию моделирования механического поведения системы, принадлежащей некоторому пространству сплошного материала Ω и соответствующему ему дискретному пространству Ω'. Исследуемый объект площади S лежит в Ω из R2. Данные компьютерной томографии, соответствующие объекту, определяются пространством Ω' и состоят из микрообластей со значениями рентгеновской плотности, которые могут быть визуализированы с использованием цветового спектра. Данные объекта дискретной области S' из Ω' могут быть получены путем разделения площади S на большое число микроподобластей (пикселей), размер которых зависит от разрешающей способности измерительного устройства. Подразумевается, что объект области S состоит из материала и пор. Первым предположением выступает наличие связи между значениями в пикселях изображения (цифрового прототипа) и упругими свойствами материала. Вторым предположением полагается, что материал сплошной и изотропный, а анизотропию на макроуровне определяет наличие и распределение пор. На практике к таким изображениям могут быть отнесены изображения, получаемые на основе рентгеновской томографии. Так, данные компьютерной томографии могут быть использованы при определении локальных механических параметров. В этом случае необходимо понимать взаимосвязь между механическими параметрами и значениями, полученными из шкалы Хаунсфилда. На рис. 1 представлено пространство сплошного материала Ω и связанное с ним дискретное пространство Ω'. а б Рис. 1. Схематическое представление исходного объекта: а - пространство сплошного материала Ω; б - дискретное пространство соответствующих данных с изображений Ω' Fig. 1. A schematic representation of the original object: (a) solid material space Ω; (b) discrete space Ω' of the corresponding data from images Определим основные возможности предложенного подхода в случае плоских деформаций. Рассмотрим принцип реализации конечного элемента неоднородной среды. За основу был взят известный метод построения четырехузлового плоского изопараметрического конечного элемента с линейной аппроксимацией геометрии и поля перемещений. Использовались следующие функции для аппроксимации радиус-вектора r и вектора перемещений ϑ: (1) (2) где x и y - вектор координат в глобальной декартовой системе координат; u и v - вектор перемещений в глобальной декартовой системе координат; n - локальный номер узла конечного элемента; Nn(ξ, η) = (1 + ξnξ)∙(1 + + ηnη)/4 - функции формы, ξn и ηn - локальные координаты узлов элемента, ξ и η - координаты в локальной системе координат. Данные соотношения могут быть переписаны в матричной форме: (3) где [N] - матрица аппроксимирующих функций; ϑe - вектор узловых перемещений. Деформация среды может быть описана с помощью компонент линейных и сдвиговых деформаций, представленных в виде приведенного вектора ε и выраженных через перемещения с помощью известных соотношений вида (4) которые также могут быть представлены в матричном виде: (5) Напряженное состояние определяется тензором напряжений через компоненты главных напряжений σxx, σyy, σzz, т.е. нормальных для случая, когда тангенциальные напряжения σxy = σxz = σzy равны нулю. Здесь под осями x-y понимаются главные оси тензора напряжений. Тензор напряжений, подобно тензору упругих деформаций, может быть записан в виде приведенного вектора σ. Закон Гука, связывающий приведенные векторы напряжений и деформаций, может быть представлен в следующем виде: (6) Для вычисления локальной матрицы жесткости использовались известные соотношения вида [28] (7) где Ve - объем конечного элемента в пространстве конечно-элементной сетки Ωe; [D(r)] - тензор упругости; [B(r)] - матрица, связывающая приведенный вектор упругих деформаций и вектор узловых перемещений: (8) С применением гипотезы об анизотропии, возникающей вследствие наличия пор, матрица жесткости конечного элемента относительно дискретного пространства Ω' может быть вычислена с помощью следующего соотношения: (9) где ω(r) - функция Ω→Ω' и [D] - тензор упругости изотропного материала, который может быть выражен через модуль упругости Юнга и коэффициент Пуассона. В следующих уравнениях будет опущено интегрирование в направлении оси z. В случае перехода к локальным координатам уравнение (9) может быть переписано в виде (10) где |J(ξ, η)| - определитель якобиана преобразования координат; he - размер элемента в направлении оси z. Численное интегрирование локальной матрицы жесткости проводилось с использованием метода средних прямоугольников. В качестве конечного элемента неоднородной среды был выбран выпуклый четырехугольник с линейными ребрами. Точки интегрирования для такого элемента представляют собой локальные координаты геометрических центров пикселей входного изображения, таким образом: (11) где ξi и ηi - локальные координаты точек интегрирования; ∆ξ и ∆η - шаг интегрирования в направлении двух локальных координатных осей; ω(ξi, ηi) - веса квадратурной формулы, определяемые значениями с изображений в текущей точке интегрирования; I и J - количество квадратурных точек в направлении двух локальных координатных осей в элементе. Для определения принадлежности пикселя области конечного элемента было использовано следующее неравенство: (12) где rA - радиус-вектор точки A (центр пикселя); rI - радиус-вектор точки I; ni - вектор нормали к грани i (i от 1 до 4); δ - параметр точности. На рис. 2 представлены численные данные (а) и геометрия конечного элемента (б), наложенная на дискретное пространство компьютерной томографии (I, J, K, L - узлы конечного элемента). Если пиксель полностью принадлежит области конечного элемента, параметр точности δ определяется нулевым значением. В процессе интегрирования возникают некоторые трудности вследствие внутренней дискретности пространства данных со входных изображений. а б Рис. 2. Схематический вид исследуемой области: а - бинаризованная область S' дискретного пространства компьютерной томографии Ω'; б - пространство конечно-элементной сетки Ωe, наложенное на дискретное пространство Ω' (зеленым цветом отмечены пиксели, полностью принадлежащие внутренней области конечного элемента, оранжевым цветом обозначены пиксели, принадлежащие области относительно их геометрического центра, синим цветом определены пиксели, касающиеся границы конечного элемента некоторой окрестностью геометрического центра, черным цветом - пиксели, лежащие вне интегрируемой области) Fig. 2. A schematic view of the investigated area: (a) binarized area S' of the discrete space Ω' of the computed tomography; (b) the space of the finite element mesh Ωe superimposed on the discrete space Ω' (the pixels that completely belong to the inner region of the finite element are marked in green, the pixels belonging to the area relative to their geometric center are marked in orange, the pixels touching the boundary of the finite element by some neighborhood of the geometric center are defined in blue; in black one can see pixels lying outside the integrated area) Существует несколько способов определения вхождения пикселя в область конечного элемента - в этом случае рассматривается некоторая окрестность точки интегрирования. На рис. 2, б представлены варианты со следующими условиями принадлежности области элемента: 1) геометрический центр пикселя лежит внутри области (зеленый и оранжевый цвета); 2) некоторая окрестность пикселя имеет ненулевое пересечение с областью элемента (оранжевый и синий цвета). В силу предположения о разрывности значений напряжений в узлах элементов необходимо провести анализ напряжений, основанный на локальном усреднении напряженного состоянии по объему каждого конечного элемента. Вектор усредненных напряжений может быть найден с помощью выражения [28, 29] (13) где V - объем конечного элемента в пространстве сплошного материала Ω; V' - объем в дискретном пространстве компьютерной томографии Ω'. Для сегментации исследуемого объекта на изображении использовались методы, основанные на оптической плотности [30, 31] с модификациями, описанными в [32, 33]. Для автоматизации построения расчетной сетки использовался подход фильтрации сетки, описанный в [34]. Оценка влияния плотности заполнения конечного элемента пикселями на точность результата подробнее представлена в работе [27]. Весовая функция может быть определена в зависимости от типа входных двумерных изображений. Так, для данных томографии существуют формулы для пересчета оптической плотности и механических характеристик (модуль Юнга, предел прочности) [17, 20]. 3. Результаты и обсуждение 3.1. Тестовые задачи Для оценки сходимости предложенного подхода были рассмотрены тестовые задачи на основе образцов прямоугольной геометрии. Алгоритм был реализован в системе компьютерной алгебры Wolfram Mathematica. Сравнение результатов было проведено в соответствии с данными, полученными в универсальном программном комплексе конечно-элементного анализа ANSYS. Таблица 1 Геометрические и физические параметры модели Table 1 Geometric and physical parameters of the model a, мм b, мм h, мм E, ГПа ν 10,0 10,0 0,1 20 0,3 В табл. 1 приведены геометрические и физические параметры тестовых задач: a и b - размеры геометрии тестового образца вдоль горизонтальной (Ox) и вертикальной (Oy) осей соответственно, h - толщина в направлении нормали к плоскости образца (Oz), E - модуль упругости Юнга, ν - коэффициент Пуассона. Для сравнения полученных результатов были рассмотрены перемещения в направлении оси Oy в узлах нагруженной линии. 3.2. Сплошной образец В первой тестовой задаче рассматривался образец с квадратной геометрией, состоящий из сплошного изотропного материала. Размер массива данных компьютерной томографии составлял 100×100 пикселей. Прикладываемая к верхней грани равномерно распределенная нагрузка равнялась 3000 Н. Нижняя грань фиксировалась в перемещениях в направлении двух координатных осей. Для сравнения результатов было произведено моделирование в системе ANSYS образцов с двумя типами сеток, состоящих из 1 и 10 000 конечных элементов соответственно. Результаты сравнивались с данными модели, построенной на основе предложенного подхода с использованием 1, 4 и 9 конечных элементов. Таблица 2 Результаты численного расчета тестовой задачи: средние перемещения в нагруженных узлах вдоль вертикальной оси Table 2 Results of numerical calculations of the test problem: average displacements in loaded nodes along the vertical axis ANSYS, 1 КЭ ANSYS, 10 000 КЭ 1 КЭ 4 КЭ ∆y, мм -1,4553 -1,48171 -1,4553 -1,4714 Исходя из предположения о непрерывности материала, весовые функции принимали значения, равные единице. Количество пикселей, принадлежащих области одного конечного элемента, варьировалось в зависимости от его размера. Численное интегрирование локальной матрицы жесткости осуществлялось на основе стандартного метода интегрирования средними прямоугольниками. В связи с этим ошибка результатов была связана с ошибкой метода интегрирования и способа усреднения значений по верхнему краю (табл. 2). а б б г Рис. 3. Деформированное состояние сплошного образца: а - геометрия модели с граничными условиями; б, в, г - 1, 4 и 9 КЭ соответственно Fig. 3. The deformed state of a solid specimen: (a) a model geometry with boundary conditions; (b), (c) and (d) are 1, 4 and 9 FEs, respectively На рис. 3 представлено деформированное состояние образца в зависимости от числа конечных элементов. Рис. 4. Поле узловых перемещений вдоль вертикальной оси (10 000 КЭ, ANSYS) Fig. 4. The field of nodal displacements along the vertical axis (10 000 FEs, ANSYS) Рис. 4 отображает деформированное состояние вдоль вертикальной оси, вычисленное в системе ANSYS. Свойства материала определяют симметрию в распределении полученных значений поля деформаций. 3.3. Образец с треугольным вырезом В следующем эксперименте рассматривался образец 100×100 пикселей с треугольным вырезом, определяющим остаточную рамку шириной 0,1 мм. В этом случае весовая функция принимала значения двух видов (0 и 1), в зависимости от локальных координат текущей точки интегрирования. Параметры модели аналогичны постановке первой тестовой задачи, прикладываемая нагрузка соответствовала 300 Н. а б в г Рис. 5. Образец с треугольным вырезом: а - геометрия выреза; б, в, г - деформированное состояние модели из 1, 4 и 9 КЭ соответственно Fig. 5. A sample with a triangular cut: (a) a cut geometry; (b), (c) and (d) is the deformed state of the model with 1, 4 and 9 FEs, respectively Полость, образованная надрезом, влияет на податливость образца к внешней нагрузке, что определяется повышенной деформацией геометрии рассматриваемой области (рис. 5, б-г). Следует отметить, что при небольшом количестве конечных элементов аппроксимация геометрии образца в деформированном состоянии является угловатой (рис. 6). В связи с этим было предложено рассматривать модели, построенные по данным компьютерной томографии с использованием 100 и 400 конечных элементов соответственно. В этом случае модуль Юнга был равен 200 ГПа, а приложенная нагрузка - 20 Н. Для автоматизации построения такой расчетной сетки использовался подход фильтрации сетки [34]. Рис. 6. Поле узловых перемещений вдоль вертикальной оси (10 000 КЭ, ANSYS) Fig.6. The field of nodal displacements along the vertical axis (10 000 FEs, ANSYS) При построении численных моделей учитывался тот факт, что конечные элементы, целиком принадлежащие областям, в которых нет материала (значения в пикселях равны нулю), определяются нулевой локальной матрицей жесткости. При отсутствии связи в узлах с элементами, соответствующими ненулевой матрице жесткости, это приводит к вырождению глобальной матрицы жесткости. Таким образом, подобные «пустые» элементы заранее исключались из сетки (рис. 7, а, б). а б Рис. 7. Деформированное состояние модели: а - 100 КЭ, б - 400 КЭ Fig. 7. The deformed state of the model: (a) 100 FE, (b) 400 FE 3.4. Образец с неравномерной плотностью В следующей тестовой задаче рассматривалось моделирование области с неравномерным распределением плотности материала. Законы распределения неоднородности взяты из работы [26]. Граничные условия задачи и механические параметры аналогичны предыдущим постановкам тестовых задач. На рис. 8, а приведено изображение, определяемое функцией а на рис. 9, а - функцией где i и j - индексы пикселя в матрице изображения. а б в г Рис. 8. Первый образец с неравномерной плотностью: а - первый тип распределения плотности; б, в, г - деформированное состояние модели с 1, 4 и 9 КЭ соответственно Fig. 8. The first sample with an uneven density: (a) the first type of density distribution; (b), (c) and (d) the deformed state of the model with 1, 4 and 9 FEs, respectively а б в г Рис. 9. Второй образец с неравномерной плотностью: а - второй тип распределения плотности; б, в, г - деформированное состояние модели с 1, 4 и 9 КЭ соответственно Fig. 9. The second sample with the uneven density: (a) the second type of density distribution; (b), (c) and (d) the deformed state of the model with 1, 4 and 9 FEs, respectively Исходя из характера деформации, можно сделать вывод, что модель, построенная по данным со входных изображений, способна описывать поведение анизотропного материала с достаточной точностью путем учета неоднородности при интегрировании матрицы жесткости каждого конечного элемента. Характер деформации геометрии моделей с анизотропными свойствами (рис. 8, б-г, 9, б-г) существенно отличается от результатов, полученных с использованием моделей, основанных на предположении об изотропии материала (см. рис. 4, б-г). При этом снижается требование к количеству элементов, а анизотропность среды учитывается за счет интегрирования локальных матриц жесткости. 3.5. Модельная задача В постановке модельной задачи использовалось изображение среза дистального отдела бедренной кости крысы, полученное с использованием данных компьютерной томографии. Разрешение исходного изображения составляло 996×795 пикселей (рис. 10, а). Для оценки влияния разрешающей способности сканирующего устройства на качество получаемых результатов исходное изображение было сжато с потерей качества. Таким образом, проводились расчеты также на аналогичных изображениях, но с уменьшенным разрешением - 332×265 и 166×132 пикселей соответственно. Таблица 3 Геометрические и механические параметры образца модельной задачи Table 3 Geometric and mechanical parameters of the model problem sample a, мм b, мм h, мм E, ГПа ν 8,79 11,01 0,1 12 0,3 Механические и геометрические параметры испытания представлены в табл. 3: a и b - размеры геометрии образца модельной задачи вдоль горизонтальной и вертикальной осей Ox и Oy соответственно, h - толщина в направлении ортогональной к плоскости образца оси Oz, E - модуль упругости Юнга, ν - коэффициент Пуассона. Сравнение полученных результатов осуществлялось на основе значений перемещений в направлении вертикальной оси Oy в узлах верхней линии (рис. 11 красным цветом обозначена точка, для которой производилось сравнение полученных величин). Геометрия образца была аппроксимирована 70 конечными элементами (7 элементов по горизонтальной оси и 10 - по вертикальной, рис. 10, б). Рассматриваемый образец подвергался воздействию одноосной продольной сжимающей вертикальной нагрузки, равномерно распределённой по верхним узлам и равной 100 Н. Величина нагрузки была выбрана на основе проводимых экспериментальных исследований. Узлы, принадлежащие нижней линии, фиксировались в перемещениях вдоль направления двух координатных осей (Ox и Oy). Для предотвращения вырождения глобальной матрицы жёсткости из сетки заранее исключались конечные элементы, соответствующие нулевой плотности, определяющей отсутствие материала (рис. 10, в). а б в Рис. 10. Дистальная часть бедренной кости крысы: а - компьютерная томография среза; б - конечно-элементная сетка, наложенная на изображение; в - сетка с граничными условиями после удаления элементов с нулевой плотностью материала Fig. 10. The distal part of the rat femur: (a) the computed tomography of the section; (b) the finite element mesh superimposed on the image; (c) the mesh with boundary conditions after removing elements with zero material density Полученная картина перемещений отражает асимметрию в геометрии испытуемого образца. Характер наблюдаемой деформации также объясняется неравномерной плотностью материала. Следует отметить, что разница в результатах, полученных по изображениям с различным разрешением, незначительна (рис. 11, а-в). а б в Рис. 11. Деформированное состояние модели: а - 996×795; б - 332×265; в - 166×132 пикселей Fig. 11. The deformed state of the model: (a) 996×795, (b) 332×265, (c) 166×132 pixels Незначительное различие результатов свидетельствует о возможности уменьшения количества точек интегрирования вследствие ослабления разрешения компьютерной томографии (табл. 4). Подобная процедура может уменьшить объем вычислений без значительного увеличения погрешности результатов. Таблица 4 Перемещения красной точки вдоль горизонтальной и вертикальной осей Table 4 Displacements of the red point along the horizontal and vertical axes 996×795 332×265 166×132 ∆x, мм 0,126686 0,126512 0,136352 ∆y, мм -0,060812 -0,060795 -0,066001 Заключение В статье представлен метод построения плоского конечного элемента, способного учитывать неравномерное распределение свойств материала на основе изображений исследуемой области. Такой подход позволяет производить построение сетки конечных элементов на основе связи между значениями в пикселях и упругими свойствами материала. В статье представлено решение показательных тестовых задач, позволяющих оценить обоснованность предложенной методики. Алгоритмы численных расчетов реализованы в собственном программном обеспечении. Представленные результаты отображают данные численного эксперимента в модельной постановке: данные компьютерной томографии среза бедренной кости крысы использовались в качестве входного изображения исходной области. Полученные значения сравнивались с результатами решения соответствующих задач в системе ANSYS. Представленные деформации рассматриваемых областей показали адекватность получаемых результатов, а также указали на возможность значительного сокращения времени расчета путем разумного уменьшения разрешения исходного изображения. В дальнейшем предполагается развитие данной методики с использованием других типов конечных элементов и улучшенных методов интегрирования конечных элементов.

About the authors

O. V Vorobiev

Kazan (Volga region) Federal University

E. V Semenova

Kazan (Volga region) Federal University

D. A Mukhin

Kazan (Volga region) Federal University

E. O Statsenko

Kazan (Volga region) Federal University

T. V Baltina

Kazan (Volga region) Federal University

O. V Gerasimov

Kazan (Volga region) Federal University

References

Statistics

Views

Abstract - 333

PDF (Russian) - 205