Full Text

Введение В настоящее время актуальной задачей является моделирование напряженно-деформированного состояния пористых или многофазных сред. Применение аппарата механики сплошных сред для моделирования пористых или многофазных сред позволит использовать имеющиеся классические решения для анализа большого количества задач. Но остается открытым вопрос приведения таких сред к сплошным. В настоящее время в связи с развитием неразрушающих методов контроля можно проводить анализ структуры различных гетерогенных материалов. Актуальной задачей является определение связи между структурными и механическими свойствами материала [1-8]. Эта задача остро стоит в клинической медицине [9-11], в частности учет изменения костной ткани под внешним силовым воздействием [12-14]. Так, информация о качестве костной ткани может сыграть решающую роль при планировании лечения, кроме того, такие данные позволят улучшить качество биомеханического моделирования суставов и органов опорно-двигательной системы [9-11]. Методы построения репрезентативной модели для анизотропных сред, в частности, в задачах биомеханики были рассмотрены ранее [15-17], в том числе и с привлечением новых математических объектов, например, таких как тензор структуры [18]. На основе обобщения рассмотренных подходов [15-17] в настоящей работе предлагается методика построения репрезентативного элемента пористого материала на основе данных компьютерной томографии (КТ). 1. Материалы и методы Построение репрезентативной модели В описываемом подходе предполагается, что сам материал изотропен, а анизотропность репрезентативного элемента достигается наличием пор. Размеры и распределение пор по объему количественно характеризуют анизотропные свойства [18-20]. Как уже упоминалось, информация о структуре материала получается на основе КТ. В этом случае входными данными для реконструкции геометрии исследуемого объекта является трехмерный массив данных коэффициента ослабления рентгеновского излучения [21]. Физический размер элемента этого массива соответствует размеру вокселя, а размер всего массива - размерности рабочего поля КТ. На практике данные коэффициента ослабления (коэффициента Хаусфилда) нормируются [21]. Для дифференциации материала от воздуха определяется порог бинаризации. Тогда массив данных можно трансформировать в бинарный, где 1 означает наличие материала, а 0 - его отсутствие (пора). Для оптимального определения порога бинаризации в работе был использован метод Отцу [22]. Для построения конечно-элементной модели строится гексагональная равномерная сетка, в которой размер каждого элемента определяется размерами вокселя. Область сетки соответствует области данных КТ. Для построения репрезентативной модели сетка модифицируется согласно бинарным данным, в этом случае удаляются элементы, соответствующие порам [23, 24]. Схематично этот процесс проиллюстрирован на рис. 1. Рис. 1. Схематичное изображение построения конечно-элементной модели Fig. 1. Schematic representation of constructing the finite element model Анизотропная репрезентативная модель Тензор упругих констант обладает базовой симметрией вида Для упрощения представления о тензоре 4-го ранга в дальнейшем будем использовать нотацию Фойгта. Благодаря этому можно записать тензор в виде симметричной матрицы 6х6. В этом случае закон Гука представим в следующем векторно-матричном виде: В общем случае тензор упругих констант содержит 21 независимую компоненту. Для определения этих констант проводились численные эксперименты. На основе метода конечных элементов для репрезентативной модели проводилось 6 испытаний: 3 испытания на одноосное сжатие (рис. 2, а) и 3 на чистый сдвиг (рис. 2, б) [18, 25-29]. В практических расчетах актуальным вопросом является выбор размера репрезентативного элемента [30], конечно-элементной сетки, используемой для элемента. В данной работе этот вопрос не рассматривается. В численных экспериментах производилось кинематическое нагружение образца. В этом случае на одной грани (пассивной) накладывались условия отсутствия перемещений (грань ABCD на рис. 2), а на противоположной (активной) - перемещения, соответствующие заданным деформациям (грань A’B’C’D’ на рис. 2). После решения задачи напряженно-деформированного состояния (НДС) на пассивной грани определяются средние по всем компонентам напряжения. Деформации, соответствующие перемещениям, должны быть такими, чтобы материал находился в упругой зоне. а б Рис. 2. Схема нагружения в численных экспериментах: а - одноосное сжатие; б - сдвиг Fig. 2. Loading scheme in numerical experiments: a - is a uniaxial compression; b - is shift Для каждого из расчетных случаев в численном эксперименте получается 6 осредненных по пассивной грани напряжений. Таким образом, для всех расчетных случаев получается 36 уравнений, в то время как неизвестных - 21. В нотации Фойгта это может быть представлено выражениями, приведенными ниже, в которых индекс j отвечает за номер численного эксперимента (в описываемой методике j от 1 до 6): Для решения переопределенной системы был применен метод наименьших квадратов (МНК). Таким образом, можно получить осредненные анизотропные свойства репрезентативной модели. Ортотропная модель подразумевает наличие трех ортогональных плоскостей симметрии. В этом случае компоненты тензора упругих констант описываются девятью независимыми значениями. Соответственно, для заданного анизотропного тензора упругих констант можно сформулировать две задачи: определение направлений осей ортотропии и определение независимых компонент. Ортотропная репрезентативная модель Для определения осей ортотропии рассмотрим преобразование компонент тензора при повороте [31]: , где Из условия обращения в ноль соответствующих компонент тензора в осях ортотропии можно определить целевую функцию [2, 16]: Эта функция в осях ортотропии должна обращаться в ноль, а значит, можно сформулировать задачу минимизации целевой функции по искомым углам [2, 15]: В осях ортотропии тензор упругих констант имеет вид На практике такая идеализация реализуется редко ввиду погрешностей, связанных с наличием шумов, бинаризацией и дискретизацией данных КТ, а также погрешностей при решении переопределённой задачи. Поэтому после поворота тензора на углы ортотропии производится принудительный пересчет компонент из условия обнуления соответствующих компонент, при этом пересчет ненулевых компонент производится из условия равенства напряжений на одном и том же поле деформаций. Таким образом, можно получить осредненные ортотропные свойства репрезентативной модели и направление осей ортотропии. Модельная задача Рассмотрим применение описанной методики на примере пористого образца, изготовленного из пластика «Ласил Каст 12». Пористая структура изготавливалась путем активной химической реакции во время отливки, рассмотрена область 20×20×60 мм. Образец подвергался сканированию (рис. 3) на компьютерном томографе Vatech PaX-I 3D [32], размер вокселя составил 0,2 мм. В этом случае возможные микропоры, обладающие меньшим размером, чем размер вокселя, влияют на интегральную величину коэффициента ослабления в сканируемой точке. При необходимости методику бинаризации можно дополнить учетом этих свойств, но в данной работе это не рассматривается. Для бинаризации полученных данных был использован метод Отцу [22]. Полученные данные дискретизировались для конечно-элементного анализа согласно описанной выше методике. Рис. 3. Данные компьютерной томографии Fig. 3. Computed tomography data В расчетах для изотропного материала были использованы следующие значения механических констант: модуль Юнга - 2 МПа, коэффициент Пуассона - 0,3. Образец был разделен на три репрезентативных элемента (рис. 4, а). Для каждого образца были рассчитаны анизотропные и ортотропные (с определением осей ортотропии) тензоры упругих констант по описанной выше методике. Расчеты выполнялись в программном комплексе ANSYS. Для расчетов использовался восьмиузловой конечный элемент с линейной аппроксимацией solid185. В табл. 1 приведены значения деформаций, используемых в численных экспериментах. Таблица 1 Величины деформаций, используемых в численных экспериментах Table 1 Strain values used in numerical experiments Заданная деформация Номер эксперимента 1 2 3 4 5 6 1/30 0 0 0 0 0 0 1/30 0 0 0 0 0 0 1/30 0 0 0 0 0 0 1/60 0 0 0 0 0 0 1/60 0 0 0 0 0 0 1/60 Описанная методика была применена для изучаемого образца. Рассмотрим подробнее компоненты тензора упругих констант на примере одного из репрезентативных элементов. Для краткости в дальнейшем будем опускать размерность компонент, понимая ее как МПа, а также приводить только половину значений компонент ввиду симметричности тензора. 2. Результаты и обсуждение После интегрирования полученных в результате решения задачи НДС и решения переопределенной системы уравнений можно получить компоненты тензора упругих констант: После минимизации целевой функции для данного тензора были определены углы поворота: . Поворот тензора на эти углы приводит тензор упругих констант к следующему виду: Как уже было отмечено выше, ввиду накопления погрешностей тензор в осях ортотропии отличается от теоретического вида. Принудительное представление тензора приводит к виду Для каждого репрезентативного элемента были определены такие матрицы анизотропии и ортотропии. Для оценки качества такого осреднения был проведён проверочный расчет. Для этого образец с нижнего торца жестко фиксировался, а к верхнему торцу прикладывалась нагрузка 400 Н (см. рис. 4, а). В расчетах рассматривалось три модели: первая модель - прямое моделирование, то есть построение конечно-элементного ансамбля на основе бинаризированных данных КТ всего образца; вторая состояла из репрезентативных элементов с анизотропными свойствами; третья состояла из репрезентативных элементов с ортотропными свойствами, определенными в найденных осях ортотропии. Таблица 2 Усредненные значения расчетных значений интенсивности напряжений (МПа) для репрезентативной области Table 2 Results of averaged stress intensities (MPa) for the representative region Тип модели Кол-во конечных эл. 1-й элемент 2-й элемент 3-й элемент Изотропная 3000000 1,23 1,20 1,30 Анизотропная 3 1,04 1,22 1,72 Ортотропная 3 23,01 21,54 18,75 24 1,78 1,40 2,10 81 1,43 1,03 0,99 192 1,35 1,01 0,99 375 1,28 1,00 0,99 375000 1,01 1,01 1,02 Для количественной оценки репрезентативных моделей использовалась интенсивность напряжений, которую стоит понимать как инвариантную величину получаемых тензоров напряжений. На рис. 4, б приведены результаты интенсивности напряжений для первой модели. Стоит отметить большое количество концентратов напряжений, связанных с численной реализацией. На рис. 4, в приведены аналогичные расчеты для второй модели. Можно отметить, что максимальные значения напряжений локализуются в тех же областях, где и концентраты напряжений для первой модели. Для третьей модели возникали большие концентраты напряжений, поэтому элементы подвергались сгущению сетки. На рис. 4, г приведены результаты интенсивности напряжений для третьей модели, в которой каждый элемент был сгущен на 53 элементов. а б в г Рис. 4. Расчетная модель (а), интенсивность напряжений (МПа): б - прямая модель; в - анизотропная модель; г - ортотропная модель Fig. 4. Calculation model (a), stress intensity (MPa): b - is a direct model; с - is an anisotropic model; d - is an orthotropic model Для количественной оценки результатов производилось осреднение интенсивности напряжений в областях репрезентативных элементов. В табл. 2 приведены результаты такого осреднения, а для ортотропной модели отдельно приведены результаты для сгущенной сетки. Можно отметить наличие сходимости при сгущении сетки для ортотропной модели. Стоит отметить достаточно большое расхождение результатов для третьего репрезентативного элемента. Это может быть объяснено увеличением пористости в образце (порядка 50 %), что наводит на мысль о несостоятельности выбранных размеров репрезентативного элемента для этой области. Выводы 1) Приведена методика для построения репрезентативной анизотропной и ортотропной моделей по данным компьютерной томографии. 2) Приведены результаты применения описанной методики. Оценка результатов производилась в терминах инварианта (интенсивность напряжений). 3) При моделировании принимались следующие допущения: сам материал предполагается изотропным, анизотропность возникает ввиду наличия пор, и количественно анизотропия характеризуется размерами и законом распределения пор. 4) Применение описанной методики позволяет на основании данных томографии строить осреднённые модели сплошной среды. Но стоит отметить необходимость оценки размера репрезентативного элемента по отношению к размерам пор для повышения качества результатов, что планируется сделать в дальнейших исследованиях.

About the authors

N V Kharin

Kazan Federal University

O V Vorobyev

Kazan Federal University

D V Berezhnoi

Kazan Federal University

O A Sachenkov

Kazan Federal University; Kazan National Research Technical University named after A.N. Tupolev

References

Statistics

Views

Abstract - 250

PDF (Russian) - 103