Full Text

Введение В настоящее время наиболее перспективным направлением в моделировании поведения гетерогенных сред является использование данных компьютерной томографии [15, 16]. Этот подход позволяет получать информацию о структуре сложных анизотропных, мелкозернистых или пористых материалов, а также иных тел, характеризующихся неоднородностью. Эта задача особенно актуальна в ортопедической клинической практике. Так, дефекты костной ткани, связанные с всевозможными патологиями, могут оказать большое влияние на качество проводимого лечения, потому на этапе диагностики необходимо получить максимально возможную информацию о прочности и жесткости костной ткани [6, 12, 18, 20, 21]. Существует несколько подходов к построению моделей элементов пористых сред, деформирующихся под действием внешних нагрузок. В первую очередь к ним можно отнести аппроксимацию распределения неоднородности методом средних пересечений линий [1, 2, 10]: в этом случае формулируются соотношения, связывающие компоненты тензора упругих констант и тензора структуры [1, 2, 7, 9, 10], характеризующего осреднённое направление пор. Другой подход - сведение анизотропии материала к ортотропии путём определения констант из численных экспериментов [5, 8, 13, 14, 17, 22, 23, 24]. В данной работе предлагается подход, основанный на учёте особенностей структуры пористого материала, выявленных по данным компьютерной томографии при построении численной модели. Проведение компьютерной томографии элемента пористой среды предполагает создание его цифрового прототипа, представляющего собой трехмерный целочисленный массив определенной структуры [15, 16]. Элементами этого массива являются единицы или нули, которые характеризуют наличие или отсутствие вещества в определенном микроэлементе объема и отражают рентгеновскую плотность согласно шкале Хаунсфилда. Таким образом, цифровой прототип представляет структуру элемента пористой среды в виде совокупности элементарных микрообъемов (кубиков), про каждый из которых известно - содержит он вещество (кости) или нет. По этим данным на основе какого-либо приближенного метода можно построить дискретную механическую модель элемента пористой среды. Наиболее удобным с алгоритмической точки зрения в этом случае методом дискретизации (а в дальнейшем и расчета) будет метод конечных элементов. Наивысшая точность расчета будет достигнута в случае моделирования каждого микрообъема трехмерным конечным элементом сплошной среды [15, 25, 26, 27]. Но в этом случае затраты ресурсов ЭВМ на создание дискретной модели, постпроцессорную обработку результатов и на этапе процессорных вычислений превысят все разумные пределы. Поэтому представляется целесообразным увеличить размеры конечных элементов, а каждый определенный в цифровом прототипе микрообъем считать окрестностью точки интегрирования при численном вычислении матрицы жесткости конечного элемента [19, 28]. Остается неясным вопрос определения числа квадратур в каждом из направлений при интегрировании локальной матрицы жесткости. Если число точек интегрирования по каждой координате мало, то точность решения задач на основе такого элемента может быть невысокой, так как в самом простом случае при интегрировании придется применять метод прямоугольников. Однако использование конечных элементов с очень большим числом квадратур хоть и повышает точность интегрирования в пределах элемента, но может увеличить жесткость расчетной области из-за возможно малого количества элементов. Поэтому целью работы представляется реализация методики статического расчета элементов трехмерных пористых объектов на основе трехмерного изопараметрического линейного конечного элемента сплошной среды, построенного на основе выявленного по данным компьютерной томографии исследуемой области ее цифрового прототипа. Материалы и методы Создание цифрового прототипа Для некоторого элемента расчетной области с помощью методов компьютерной томографии определяется соответствующий ему набор данных, представляющих собой трехмерную структуру. Интенсивность окрашивания микрообъемов структуры, слагающих расчетную область, отражает их рентгеновскую плотность согласно шкале Хаунсфилда. Процесс создания цифрового прототипа исследуемого объекта предполагает его разделение на большое число виртуальных микрокубиков (вокселей) размером с координатами центра вокселя , где линейные размеры обычно бывают равны друг другу и определяются разрешающей способностью компьютерного томографа. Значения данных компьютерной томографии, соответствующие вокселю, бинаризируются по заданному порогу, определяемому, например, методом Отсу, отделяя плотную костную структуру от вещества в порах: если значение компьютерной томографии в текущем вокселе больше порогового значения, то этот воксель характеризует костную ткань, в противном случае - пору (рис. 1). Рис. 1. Значения данных компьютерной томографии после бинаризации для одного слоя в локальной области Численное интегрирование матрицы жесткости с учетом данных компьютерной томографии Рассмотрим методику построения конечного элемента пористой сплошной среды на основе широко известной методики [3, 4, 11, 29] построения восьмиузлового трехмерного изопараметрического конечного элемента сплошной среды с линейной аппроксимацией геометрии и поля перемещений в локальных координатах. В рамках изопараметрической концепции для аппроксимации геометрии (радиуса-вектора точки и исходных перемещений (вектора перемещений точки ) используется одинаковая система функций: , (1) , (2) где - известные линейные функции формы, - локальные координаты узлов элемента, - проекции вектора перемещений на орты глобальной декартовой системы координат . Соотношения (2) могут быть переписаны в матричной форме , (3) где - матрица аппроксимирующих функций, - вектор узловых перемещений. Деформацию среды описывают компоненты линейных деформаций и деформации сдвига , которые представляются в виде приведенного вектора и выражаются через перемещения (2) известными соотношениями Коши , (4) которые, в свою очередь, также могут быть записаны в матричной форме . (5) Напряженное состояние представляет тензор напряжений Коши в виде компонент нормальных и касательных напряжений, который, как и тензор малых деформаций, можно записать в форме приведенного вектора . Закон Гука, связывающий приведенные векторы напряжений и деформаций, представляется в форме , (6) где - матрица упругих постоянных кусочно-однородного тела, которую можно представить в виде [5] . (7) Здесь - тензор упругости для изотропного материала, - постоянные Ламе, определяемые через модуль упругости Юнга и коэффициент Пуассона в виде ; (8) - некая весовая функция, определяемая точкой в пространстве по данным компьютерной томографии. Для вычисления матрицы жесткости конечного элемента используются известные соотношения [4, 29] вида , (9) где - матрица, связывающая приведенный вектор малых деформаций и вектор узловых перемещений [3, 4, 29]: . (10) При переходе к локальным координатам соотношение (9) можно переписать в виде (11) где - детерминант матрицы Якоби преобразования координат. В отличие от построения локальной матрицы жесткости изотропной сплошной среды, когда для вычисления интеграла (11) используется метод Гаусса, для пористых сред приходится использовать метод центральных прямоугольников. Конечный элемент пористой сплошной среды, как и в случае изотропии, представляет собой выпуклый шестигранник с линейчатыми четырехузловыми боковыми поверхностями. Квадратурными точками в таком элементе являются координаты центров вокселей из цифрового прототипа модели, только их, как и размеры вокселей, необходимо пересчитать в локальных координатах для каждого элемента. Тогда (12) где - локальные координаты точек интегрирования; - величина шага в трёх направлениях в локальных координатах; - веса квадратурной формулы, определяемые значением данных компьютерной томографии в точке интегрирования; - число квадратурных точек по каждой локальной координате в элементе. В ходе работы была реализована программа на языке программирования С++ для расчётов с использованием восьми узлового конечного элемента с тремя степенями свободы в каждом узле: перемещения вдоль соответствующих декартовых координатных осей. Решение задач осуществлялось на компьютере следующей комплектации: процессор - AMD Ryzen 7 1700 с 8 физическими и 16 логическими ядрами с частотой 3,7 ГГц, оперативная память - G.Skill Aegis 16GB DDR4 16GISB K2 3000 C16 с частотой 2933 МГц, материнская плата - MSI B350M MORTAR. Ввиду больших трудозатрат при интегрировании локальной матрицы жесткости в программе было реализовано распараллеливание алгоритма с использованием технологии OpenMP, в рамках которого расчет локальной матрицы жесткости для каждого конечного элемента производился на разных потоках. Полученные результаты визуализировались в программе ParaView. Результаты и обсуждение Для оценки влияния плотности заполнения конечного элемента вокселями были решены тестовые задачи, результаты которых сравнивались с результатами решения в пакете Ansys 14.5. Тестовые задачи Для оценки влияния плотности заполнения конечного элемента вокселями были решены тестовые задачи для одного конечного элемента в условиях различного нагружения. Рассмотрим полученные результаты для нагружения сжимающей силой: к верхней грани прикладывалась нагрузка Н, равномерно распределённая по четырём узлам; нижние узлы фиксировались в перемещениях по трём направлениям, модуль Юнга составлял 2 ГПа, коэффициент Пуассона - , размер куба - мм. Варьируя величину вокселей и устанавливая в каждом из них значение из данных компьютерной томографии, равное единице, произвели серию численных экспериментов для сплошного материала. В табл. 1 отображены количество вокселей и их соответствующие размеры, используемые при расчётах. На рис. 2 приведены результаты, полученные численным интегрированием методом центральных прямоугольников для сплошного изотропного материала с использованием разного количества вокселей внутри конечного элемента. Таким образом, при увеличении количества вокселей решение стремится к результатам, полученным в программном комплексе Ansys. Исходя из зависимости времени выполнения программы от количества вокселей внутри конечного элемента, можно сделать вывод о том, что для наиболее эффективного решения (оптимального по точности и времени выполнения) предпочтительно использовать порядка 50 вокселей на стороне куба. Результаты для других случаев нагружения конечного элемента аналогичны. С целью оценки влияния густоты сетки при постоянной плотности компьютерной томографии были решены тестовые задачи для одноосного сжатия и изгиба. Задача ставилась для геометрии в виде параллелепипеда с одним характерным размером. В этом случае на одном торце фиксировались все перемещения в узлах, на противоположном - прикладывалась распределённая по узлам нагрузка (см. рис. 3). Рассматривался образец со следующими параметрами данных компьютерной томографии: количество вокселей в направлении осей Ox, Oy, и Oz - соответственно, размер вокселя - мм. Свойства материала использовались из предыдущей тестовой задачи. Параметры конечно-элементной сетки приведены в табл. 2. Исследовались отдельно два вида нагружения: влияние сжимающей нагрузки вдоль оси Oy и действие изгибающей силы в направлении оси Oz (решение для силы, действующей вдоль оси Ox, аналогично в силу симметрии). На рис. 4 отображены результаты, полученные для трёх видов конечно-элементной сетки. Следует отметить, что относительная погрешность решения несущественно зависит от количества вокселей внутри одного конечного элемента, что позволяет в допустимых пределах изменять размер и количество конечных элементов с целью лучшей аппроксимации геометрии. Таблица 1 Количество вокселей на ребре куба и их соответствующие размеры Nvox 10 16 20 25 32 40 50 64 80 100 125 160 200 - Svox, мм 2 1,25 1 0,8 0,625 0,5 0,4 0,3125 0,25 0,2 0,16 0,125 0,1 - 0,07 0,06 0,05 0,04 0,03 0,02 0,01 0 0 50 100 150 200 250 а б Рис. 2. Результаты численных расчётов: а - относительная погрешность значений перемещений для узлов с приложенной нагрузкой; б - время выполнения программы, с; ось абсцисс - количество вокселей по одной координате A A-A b h P2 P2 P1 A x z L y z Рис. 3. Расчётная схема: P1, P2 - прикладываемая сила; b, L и h - параметры геометрии в направлении осей Ox, Oy и Oz соответственно Таблица 2 Параметры конечно-элементной сетки Количество вокселей в одном конечном элементе 1 5 1 2 10 2 3 20 4 40 35 30 25 20 15 10 5 0 1003 253 503 503 1003 0 50 100 150 200 250 300 350 Рис. 4. Относительная погрешность значений перемещений для узлов с приложенной нагрузкой: красная линия - результаты для модели в условиях действия сжимающей нагрузки, синяя линия - под действием изгибающей силы; ось абсцисс - количество конечных элементов Модельная задача В модельной задаче использовались данные компьютерной томографии для диафиза бедренной кости крысы (рис. 5, а). Сканирование выполнялось с применением микро- и нанофокусной системы рентгеновского контроля для компьютерной томографии и 2D-инспекции Phoenix V|tome|X S240 в лаборатории рентгеновской компьютерной томографии Института геологии и нефтегазовых технологий Казанского (Приволжского) федерального университета. Система оснащена двумя рентгеновскими трубками: микрофокусной с максимальным ускоряющем напряжением 240 кВ мощностью 320 Вт и нанофокусной с максимальным ускоряющем напряжением 180 кВ мощностью 15 Вт. Для первичной обработки данных и создания объёмной (воксельной) модели образца на базе рентгеновских снимков (проекций) использовалось программное обеспечение datos|x reconstruction. Зафиксированный в держателе образец помещался на вращающийся столик камеры рентгеновского компьютерного томографа на оптимальном расстоянии от источника рентгеновского излучения. Съёмка проводилась при ускоряющем напряжении 90-100 кВ и токе 140-150 мА. Размер исследуемой области - мм, количество вокселей в направлении соответствующих координатных осей Ox, Oy и Oz - , размер вокселя - мкм. а б в Рис. 5. Исследуемый образец диафиза бедренной кости крысы: а - представление данных компьютерной томографии; б - сетка 1; в - сетка 3 В табл. 3. отображены параметры расчетных конечно-элементных сеток, а на рис. 5, б и в приведены конечные элементы для сетки 1 и 3 соответственно. Узлы нижнего торца фиксировались в перемещениях по трём направлениям, на верхнем торце прикладывалась равномерно распределённая по узлам продольная сжимающая сила, равнодействующая которой равна Н. Модуль упругости Юнга принимался равным 30 000 МПа, коэффициент Пуассона - 0,3. На рис. 6, а и б приведены распределения осевых перемещений, полученных с использованием двух типов конечно-элементных сеток. На сетке с меньшим количеством конечных элементов (первая сетка) значения перемещений усредняются одним конечным элементом по толщине, при дискретизации по толщине двумя конечными элементами (для второй и третьей сеток) порядок аппроксимации распределения перемещений на торце повышается. На рис. 6, в показаны доверительные интервалы средних перемещений узлов верхнего торца объекта для конечно-элементных сеток различной густоты, увеличение стандартного отклонения объясняется увеличением количества узлов на исследуемом торце и, как следствие, изменением формы его деформированного состояния. Время выполнения расчетов составило: сетка 1 - ч, сетка 2 - ч, сетка 3 - ч. Видно, что затраченное время реализации алгоритма решения задачи на ЭВМ для различных сеток сопоставимо. Это объясняется тем, что формируемые матрицы жесткости конечного элемента для разных сеток вычисляются в одних и тех же квадратурных точках. Однако эта ресурсозатратная операция для конкретной расчетной области проводится только один раз, а при дальнейших расчетах сформированная ранее матрица жесткости просто считывается с носителя информации. Таблица 3 Параметры разбиения на конечно-элементную сетку Сетка 1 Сетка 2 Сетка 3 По толщине 1 2 2 По высоте 4 10 10 В окружном направлении 11 30 60 -0,043 -0,035 -0,03 -0,025 -0,02 -0,015 -0,01 0 -0,16 -0,12 -0,1 -0,08 -0,06 -0,04 0 а б -0,14 -0,12 -0,1 -0,08 -0,06 -0,04 -0,02 0 44 600 1200 в Рис. 6. Распределение осевых перемещений, мм, в рассматриваемом образце в условиях действия равномерно распределённой по узлам сжимающей нагрузки: а - сетка 1; б - сетка 3, в - среднее продольное перемещение верхнего торца со стандартным отклонением при различном разбиении (ось абсцисс - количество конечных элементов) Заключение 1. В работе представлен один из возможных численных подходов к описанию процессов деформирования гетерогенных сред под действием внешних нагрузок, основанный на моделировании структуры расчетной области с учетом данных ее компьютерной томографии. 2. Для оценки влияния данных компьютерной томографии на сходимость предложенной численной методики были решены тестовые задачи. Полученные результаты применялись для определения подходящих параметров (числа вокселей в конечном элементе) модельной задачи. 3. Полученные численные результаты отображают влияние точности аппроксимации геометрии образца, а также иллюстрируют зависимость поля перемещений и напряженно-деформированного состояния от структуры материала. 4. Данный подход позволяет исследовать поведение пористого материала под действием внешних нагрузок на основе оптической плотности материала. Исследования могут быть расширены с использованием другого вида обработки данных компьютерной томографии, путем усовершенствования формулы численного интегрирования локальной матрицы жесткости, а также с применением иных типов конечных элементов. Вопрос построения полей напряжений и деформаций остался вне рамок исследований данной работы, поскольку прямое построение напряженно-деформированного состояния в рамках предложенной модели метода конечных элементов проводится легко и однозначно, но возникает вопрос достоверности полученных результатов. Вполне возможно, что наилучшую точность можно будет получить путем осреднения напряженно-деформированного состояния, вычисленного в вокселях, но этот дискуссионный вопрос требует дополнительных экспериментальных исследований.

About the authors

O. V Gerasimov

D. V Berezhnoi

P. V Bolshakov

E. O Statsenko

O. A Sachenkov

References

Statistics

Views

Abstract - 75

PDF (Russian) - 64

Refbacks

• There are currently no refbacks.