Full Text

Введение В современных научно-исследовательских и опытно-конструкторских работах все активнее развивается практика математического моделирования объектов методом конечных элементов (МКЭ), который позволяет решать широкий спектр задач и, в частности, эффективно проводить анализ напряженно-деформированного состояния (НДС) твердых тел. Развитие представленного математического моделирования направлено на получение конечно-элементной (КЭ) модели деформируемого твердого тела (ДТТ), максимально соответствующей реальному объекту, что напрямую зависит от задания в модели геометрических параметров и механических характеристик. Определение этих характеристик в КЭ-модели ДТТ в современной практике предполагает допущения, вызванные, например, их осреднением, что в итоге оказывает влияние на точность результатов использования (анализа) математической модели. То есть традиционно механические характеристики ДТТ определяются на основе результатов испытания стандартных образцов, что и приводит к использованию в КЭ-модели их осредненных значений. Тогда как в пределах детали и даже образца наблюдается существенная неоднородность (рис. 1) механических характеристик материала. Представленные упрощения в модели реальных изделий (ДТТ) могут привести к неточной оценке их надежности и долговечности. а б в Рис. 1. Примеры материалов с неоднородной структурой: а - костная ткань [1]; б - металл; в - композиционный материал Fig. 1. Examples of materials with heterogeneous structure:а - bone; b - metal; c - composite material 1. Материалы и методы исследования Для решения представленной проблемы, построения и анализа НДС КЭ-моделей реальных ДТТ, в настоящей работе предлагается использовать технологию сканирования [2], которая может быть построена на любом физическом принципе (рентгеновский, акустический, радиационный и др.) с получением в результате цифрового изображения в некотором наборе сечений изучаемого ДТТ. В настоящее время, например, интенсивно развивается и используется сканирование компьютерным томографом (КТ). Однако на пути реализации этого подхода, определения реального изменения механических характеристик материала ДТТ и его индивидуальной геометрии стоит проблема разработки комплекса математических методов и алгоритмов компьютерной реализации, позволяющих использовать технологию КТ-сканирования. Рис. 2. Общая блок-схема математического моделирования ДТТ с неоднородными механическими характеристиками материала и геометрией на основе сканирования Fig. 2. Flowchart of mathematical modeling of deformable solid body with heterogeneous mechanical properties of the material and geometry based on scanning Общая блок-схема алгоритма построения КЭ-моделей реальных ДТТ с использованием технологии КТ-сканирования представлена на рис. 2. На начальном этапе (см. рис. 2, блок 1) проводится предварительная оценка неоднородности структуры и геометрии изучаемого ДТТ по параметрам: степень неоднородности структуры и механических характеристик; особенность и степень нерегулярности геометрии, наличие контактирующих тел, внешние условия, обзор данных натурных испытаний и др. В зависимости от результатов предварительной оценки изучаемого ДТТ проводятся выбор типа КТ (рис. 3) и формирование настроек (рис. 2, блок 2) процесса сканирования [3, 4, 5]. Важными настройками для получения необходимого результата являются только некоторые. К ним относятся: шаг сканирования, основное и дополнительные направления сканирования, размер и тип формата полученных данных сканирования, которые формируются в виде растровых изображений сечений ДТТ (см. рис. 2, блок 3). Каждое из этих изображений сохраняется в памяти компьютера в виде матрицы индексов цвета (см. рис. 2, блок 4). На рис. 4 представлено определение системы координат пакета растровых изображений и способ выделения области изучаемого объекта на изображении. Для примера выбрано изображение сечения бедренной кости человека, в котором индексы цвета изменяются от 0 (абсолютно черный цвет) до 255 (абсолютно белый цвет). Геометрия сечения определяется, как внешним контуром, так и контурами внутренней структуры материала в нем (см. рис. 2, блок 5). Для реализации их построения используются два этапа. На первом из них, проводится выделение областей с ярко выраженными изменениями механических характеристик. Оно осуществляется с использованием математических функций «0-1» и «0-1-2». Его целью является построение предварительных контуров в сечении (внешней геометрии и внутренней структуры) и определение временного центра тяжести в них (рис. 5). Линия контура определяется координатами пикселей на краях областей. На втором этапе относительно контуров, полученных на первом этапе, проводится уточнение геометрии внешнего контура и внутренних контуров структуры материала в сечении. Оно осуществляется с использованием метода исследования градиентов изменения индексов цвета пикселей [6, 7, 8]. а б Рис. 3. Медицинский (а) и промышленный (б) томографы и пакет растровых изображений, полученных при сканировании ДТТ Fig. 3. Medical (a) and industrial (b) tomography and a package of raster images obtained by scanning deformable solids а б Рис. 4. Система координат растрового изображения (а) и выделение в нем области исследования (б) Fig. 4. The coordinate system of the raster image (a) and the field of investigation in it (b) а б Рис. 5. Предварительное построение контура геометрии и внутренней структуры ДТТ с помощью специальных математических функций «0-1» (а) и «0-1-2» (б) Fig. 5. Preliminary construction of the geometry contour of deformable solids and its internal structure with special mathematical functions "0-1" (a) and "0-1-2" (b) Рис. 6. Уточнение контура геометрии ДТТ и его внутренней структуры на основе исследования градиента изменения индекса цвета Fig. 6. The contour refinement of deformable solids geometry and its internal structure based on the color index change gradient Рис. 7. Исследование изменения градиента по лучам, перпендикулярным предварительной геометрии Fig. 7. Investigation of the gradient variation along the rays perpendicular to the preliminary geometry а б в Рис. 8. Каркасная модель ДТТ бедренной кости (а) и верхней части зуба человека (б); корректировка центров тяжести сечений (в) Fig. 8. Frame model of a deformable solid body (bones (a) and upper part of a human tooth (b)); correction of the gravity centers of the sections (c) Контур сечения может иметь как симметричный, так и произвольный вид. Исследование градиента изменения индексов цвета пикселей производится на отрезках, перпендикулярных предварительному контуру сечения (рис. 6, 7). Более подобное описание представленного метода определения контуров геометрии показано в работах [6, 9, 10] и реализовано посредством специально разработанной программы [11, 12] обработки и идентификации растровых изображений сканирования КТ ДТТ. Полученные контуры всех сечений ДТТ объединяются для создания каркасной модели (рис. 8, а, б) относительно направления КТ-сканирования. Расстояние между сечениями равно шагу сканирования. Центр тяжести сечений (рис. 8, в) чаще всего не является вертикальной осью каркасной модели. При необходимости, например при изучении чистого сжатия или растяжения, проводится корректировка положения центров тяжести сечений по осям используемой системы координат. Дальнейшая работа общего алгоритма (см. рис. 2, блок 6) связана с тем, что каркасная модель ДТТ используется для генерации сетки его КЭ-модели [9, 13]. В зависимости от степени регулярности геометрии ДТТ используются два метода генерации сетки: 1. Регулярность геометрии при переходе от каркасной к КЭ-модели позволяет использовать КЭ типа гексаэдр (HEX) напрямую. 2. В случае нерегулярности геометрии переход от каркасной к КЭ-модели осуществляется с дополнительным построением объёмной модели твердотельной геометрии ДТТ, и чаще всего используется КЭ типа тетраэдр (TET). Процесс генерации сетки КЭ-модели проводится в программном комплексе пре- и постпроцессора MSC Patran [13, 14, 15]. Для определения механических характеристик материала в узлах и конечных элементах (см. рис. 2, блок 7) используются подходы, предложенные в работах [1, 6, 9, 16], с помощью специальной программы [11, 12]. Последовательно выполняются следующие действия: - вычисляется весовой коэффициент для определения модуля упругости относительно значения индекса цвета пикселей; - совмещаются узлы сетки КЭ-модели с координатами пикселей; - определяются механические характеристики материала узлов и элементов КЭ-модели; - осуществляется преобразование структуры данных в виде файла *.pcl для возможности использования в КЭ-модели [16, 17, 18, 19]. Предполагается, что зависимость между индексами цветов пикселей и механическими характеристиками является линейной, что впоследствии подтверждается через оценку сходимости результатов расчета с результатами эксперимента. Таким образом, значение весового коэффициента модуля упругости определяется по формуле , (1) где - среднее значение модуля упругости, полученное на образцах стандартных испытаний; - среднее значение индекса цвета относительно всех сечений. Наиболее точным является определение среднего значения индексов цветов на основе математического ожидания для всех сечений по формуле , (2) где - значение индексов цветов; - вероятность появления . Совмещение узлов сетки КЭ-модели с координатами пикселей (рис. 9) и определение индексов цвета в узлах КЭ-модели ДТТ проводятся на основе интерполяции узла относительно сечений каркасной модели (рис. 10). В целом различаются два случая: 1. Если узлы КЭ-модели лежат в плоскости какого-либо сечения каркасной модели, то узлу присваивается индекс цвета, полученный посредством математического ожидания относительно соседних с узлом пикселей (рис. 9), согласно формуле , (3) где - значение индексов цветов ближайших к узлу пикселей; , - расстояния до центров тяжести, ближайших к узлу 4 пикселей. Рис. 9. Определение значения индекса цвета в узле КЭ-модели Fig. 9. Determining the color index value in the FE model node 2. Если узлы КЭ-модели лежат между сечениями (рис. 10, а), то дополнительно, производится процедура линейной интерполяции относительно сечений, между которыми находится рассматриваемый узел (рис. 10, б). На рис. 10 представлен самый общий случай интерполяции по прямой наклонной к продольной оси каркасной модели сканируемого ДТТ. То есть завершающая формула [9] определения индекса цвета nP в узле Р на основе интерполяции имеет вид . (4) В определении модуля упругости каждого конечного элемента КЭ-модели относительно номеров индексов цвета, полученных в узлах КЭ-модели, используется среднее арифметическое значение индексов цвета его узлов: , (5) где - весовой коэффициент модуля упругости, определяемый по формуле (1); - индекс цвета i-го узла; - число узлов в конечном элементе. Более подробно данный метод представлен в работах [1, 9, 16], его алгоритмизация реализована в виде специальной программы [11, 12]. Результатом выполнения всех предыдущих этапов (рис. 2, блоки 1-7) является КЭ-модель (рис. 11, 12), в которой представлено поле изменения модуля упругости как основной механической характеристики ДТТ и индивидуальность его геометрии. Дальнейшее построение КЭ-модели (см. рис. 2, блок 8) связано с заданием в ней граничных условий закрепления, внешней нагрузки и условий контактного взаимодействия ДТТ. После чего производится анализ КЭ-модели (см. рис. 2, блок 9). а б Рис. 10. Интерполяция по прямой наклонной к оси Z Fig. 10. Interpolation in a straight line inclined to the Z axis а б в Рис. 11. Распределение неоднородности модуля упругости: а - растровое изображение сканирования; б - изменение индексов цвета в сечении; в - изменение модуля упругости в сечении в КЭ-модели Fig. 11. Inhomogeneity distribution of the elasticity modulus: a - raster scan image; b - change of color indices in the section; c - change in modulus of elasticity in cross section in the FE model На рис. 11 представлены: растровое изображение сканирования в сечении (рис. 11, а); изменение значений (неравномерность) индексов цвета (рис. 11, б) и полученное на их основе поле изменения модуля упругости в сечении КЭ-модели (рис. 11, в). а б Рис. 12. Распределение величины модуля упругости в объемных КЭ-моделях: а - КЭ типа тетраэдр; б - КЭ типа гексаэдр Fig. 12. The distribution of the modulus of elasticity in FE models: a - model with tetrahedron; b - model with hexahedron На рис. 12 представлено изменение модуля упругости в объемной КЭ-модели ДТТ в виде фрагмента бедренной кости человека с использованием конечных элементов типа тетраэдр (рис. 12, а) и гексаэдр (рис. 12, б). Анализируя результаты, представленные на рис. 11 и 12, необходимо отметить, что область сечения имеет условно две зоны: силовую, где E изменяется от 1,7E4 до 2,1E4 МПа, и несиловую, где E изменяется от 0,8E4 до 1,7E4 МПа. На рис. 13 и 14 представлено доказательство достоверности результатов разработанного математического моделирования относительно данных натурного эксперимента [20, 21, 22, 23, 24]. На рис. 13, а представлен принцип вырезания образцов из фрагмента реальной бедренной кости человека, используемых для исследования их механических характеристик, главным образом модуля упругости. Образцы вырезаются по главным осям: вертикальной (осевой) - 1, окружной (тангенциальной) - 2 и радиальной - 3 с размерами 3,5×3,5×10 мм. Вдоль этих осей и проводится исследование. Также представлен образец для уточнения механических характеристик в осевом направлении (на рис. 13, а, слева). а б Рис. 13. Ориентация образцов по осям: 1 - вертикальной (осевой); 2 - окружной; 3 - радиальной Fig. 13. Orientation of samples: 1 - axial; 2 - circumferential; 3 - radial На рис. 13, б представлены аналогичные по геометрическим размерам виртуальные образцы, вырезанные из КЭ-модели. Результаты исследования представлены на рис. 14 в виде доказательства точности и сходимости (рис. 14, а) численного решения МКЭ [9, 25] относительно данных натурного эксперимента. Результаты показывают (рис. 14, б), что, для достижения необходимой точности численного решения МКЭ требуется КЭ-сетка плотностью 3 и более конечных элемента на 1 мм. Также результат показывает, что свойство неоднородности механических характеристик материала ДТТ может быть отражено набором конечных элементов, каждый из которых со своим модулем упругости и изотропной структурой материала. Число конечных элементов 2 3 1 Рис. 14. График сходимости нормального напряжения в центре образцов, в направлениях: 1 - осевом; 2 - окружном; 3 - радиальном Fig. 14. Convergence graph of normal stress in the center of the sample in: 1 - axial; 2 - circumferential; 3 - radial 2. Построение и анализ НДС КЭ-моделей реальных ДТТ В качестве примера исследования НДС реальных ДТТ используются следующие объекты: зуб человека в челюсти и его зуб в челюсти с композитной пломбой. Представленный выбор объектов исследования не принципиален, но обусловлен двумя основными обстоятельствами: высокой степенью неоднородности механических характеристик материала кости с индивидуальной её геометрией, а также высоким уровнем развития технологии и качества сканирования КТ в медицине и технике. Также важным преимуществом представленного объекта исследования является то, что в практике изучения механических характеристик костной ткани человека существует достаточно большое количество результатов натурных испытаний [20, 21, 26, 27, 28, 29] (экспериментальных опытов). На рис. 15, в представлен вид сечения зуба человека в челюсти с его внутренней структурой, включающей области: эмаль, дентин, пористый канал и др [30]. На растровом изображении сканирования эти области обозначены (рис. 15, а, б) и имеют свои диапазоны изменения индексов цвета: эмаль (191-255), дентин (70-190), пористое волокно каналов (16-69). Для построения индивидуальной геометрии зуба (рис. 16) используется каркасная модель, которая состоит из 60 горизонтальных сечений (образец 1) и 31 вертикального сечения (образец 2). Это усложнение объясняется тем, что ДТТ в виде зуба человека в челюсти имеет высокий уровень нерегулярности геометрических параметров. На основе каркасной модели зуба строится трехмерная геометрическая модель, которая используется для генерации объемной сетки КЭ-модели с применением конечных элементов типа тетраэдр. а б в Рис. 15. Структура механических характеристик ДТТ в виде зуба человека: а - исходное растровое изображение КТ-сканирования; б - результат его обработки; в - эскиз структуры зуба человека Fig. 15. Structure of mechanical characteristics of deformable solids in the form of a human tooth: a - the initial raster image of the CT scan; b - the result of its processing; c - sketch of the human tooth structure а б Рис. 16. Каркасная модель зубов человека: а - образец 1; б - образец 2; построены на основе обработки растровых изображений их сканирования КТ Fig. 16. A frame model of human teeth, based on the processing of raster images of their scanning CT а б Рис. 17. Схема контактного взаимодействия КЭ-модели зуба: а - без продукта; б - с продуктом Fig. 17. Scheme of contact interaction of tooth FE model: a - without product; b - with product На рис. 17 представлена расчетная схема анализа НДС зуба человека в челюсти с граничными условиями и внешней нагрузкой. Кроме этих параметров, на схеме дополнительно представлены условия для решения задачи контактного взаимодействия зубов «без продукта» (рис. 17, а) и «с продуктом» (рис. 17, б). На рис. 18 представлено поле изменения свойств материала в модели в виде распределения величины модуля упругости в его объемной КЭ-модели. Это распределение полностью соответствует общим представлениям о структуре свойств материала зуба. Ниже, на рис. 19 и 20, представлены результаты анализа НДС КЭ-моделей зуба человека в челюсти в виде распределения эквивалентных напряжений (показаны в некотором срединном сечении). Эти результаты представлены в сравнении: на рис. 19, б и рис. 20, б приведены напряжения при использовании моделирования структуры механических характеристик, полученных на основе сканирования; на рис. 19, в и рис. 20, в приведены напряжения при использовании осреднения механических характеристик в КЭ-модели. Рис. 18. Трехмерная КЭ-модель зуба с распределением механических характеристик материала (модуля упругости) Fig. 18. Three-dimensional of FE model with the distribution of mechanical characteristics of the material (elastic modulus) а б в Рис. 19. Результат анализа НДС - эквивалентное напряжение КЭ-модели зуба (а) c учетом реального изменения механических характеристик материала (б) и без учета этого изменения (усредненного значения) (в) Fig. 19. The result of stress-strain state analysis - equivalent stress of FE model tooth (a) based on actual changes in the material and mechanical characteristics (b) and without this change (averaged value) (c) а б в Рис. 20. Результат анализа НДС - эквивалентное напряжение КЭ-модели зуба с продуктом (а) c учетом реального изменения механических характеристик материала (б) и без учета этого изменения (усредненного значения) (в) Fig. 20. The result of stress-strain state analysis - equivalent stress of FE model tooth with product (a) based on actual changes in the material and mechanical characteristics (b) and without this change (averaged value) (c) Результаты анализа НДС моделей (рис. 19, б; рис. 20, б) с учетом моделирования неоднородности механических характеристик материала показывают, что максимумы напряжений наблюдаются в точках контакта на вершинах зубных бугров. Большую часть нагрузки воспринимают эмали и часть дентина зуба. Пульпа, которая располагается во внутренней части зуба, практически не воспринимает нагрузку. Картина распределения напряжений соответствует реальной работе зуба, когда нагрузка воспринимается эмалью и дентином. Без учета моделирования неоднородности значений механических характеристик материала (рис. 19, в; рис. 20, в) с использованием их осреднения анализ НДС ДТТ показывает некорректную картину распределения напряжений. В частности, качественная картина совпадает, например, в виде концентрации напряжений в местах контакта между зубами. Однако количественная картина величин напряжений в области зубов не характеризует их работу как некоторую систему деформируемого тела (эмаль - дентин - пористая ткань). 3. Построение и анализ НДС КЭ-модели зубов человека с композитной пломбой в контакте Пломба в стоматологии используется для заполнения в зубе полости, возникшей в результате его повреждения. Она способствует восстановлению его нормальной функции. С технической (инженерной) точки зрения пломба является одним из составляющих элементов сборной конструкции и может изготавливаться из материалов широкого класса: амальгамы, пластмассы, керамики, металла и др. Наибольший интерес на современном этапе развития этих материалов представляют пломбы, изготовленные из композитов на полимерной армированной основе, не уступающих по механическим характеристикам и прочности эмали естественного зуба [31, 32, 33, 34]. а б в Рис. 21. Результат анализа НДС КЭ-модели зуба: а - трехмерная модель зуба с композитной пломбой; б - с повреждением, в - с пломбой Fig. 21. Result of the stress-strain state analysis FE model tooth: a - three-dimensional tooth model with a composite seal; b - with damage, c - with a filling Композитный материал отличает способность противостоять значительным вертикальным и горизонтальным нагрузкам, при этом его прочность на сжатие колеблется от 220 до 450 МПа. В представленном расчете для пломбы используется материал Grandio с модулем упругости 13910 МПа. Для моделирования повреждения (например, в виде кариеса) используется геометрия пломбы (рис. 21, а), для которой задается материал с очень низким модулем упругости, что практически исключает этот материал из работы деформируемого тела. Результаты анализа представлены на рис. 21, б, в. В первом случае (рис. 21, б), рассматривается вариант верхнего зуба с повреждением эмали. Картина напряжений показывает, что часть зуба с повреждением не воспринимает нагрузку. После замены поврежденной части зуба пломбой из композитного материала результат анализа (рис. 21, в) показывает, что распределение напряжений в зубе практически соответствует распределению напряжений в здоровом зубе (см. рис. 19, б), что свидетельствует о восстановлении функциональности зуба. Заключение Представленный комплекс математических методов моделирования интерпретации растровых изображений сканирования КТ позволяет проводить построение КЭ-моделей для любых сложных структур механических характеристик материала реальных ДТТ и их индивидуальной геометрии. Результаты решения контактной задачи взаимодействий ДТТ также подтверждают точность и надежность применения представленной технологии КЭ-моделирования на основе сканирования реальных ДТТ с учетом неоднородности механических характеристик материала и индивидуальности геометрии. Разработанная технология математического моделирования может быть применена для любых физических принципов сканирования, таких как рентгеновский, ультразвуковой, лазерный и др., и для любых типов материалов, когда полученная в результате сканирования информация сформирована в виде цифрового (растрового) изображения.

About the authors

A A Pykhalov

Irkutsk State University of Railway Transport; Irkutsk National Research Technical University

V L Duong

Irkutsk National Research Technical University

V G Tolstikov

Irkutsk National Research Technical University

References

Statistics

Views

Abstract - 274

PDF (Russian) - 85