Construction of a Representative Model Based on Computed Tomography

Abstract


The simulation of the stress-strain state of porous or multiphase media is an important task nowadays. The application of the mathematical model of continuum mechanics to such media will make it possible to extend the scope of the problems to be solved. The development of non-destructive methods of control, such as computed tomography, allows obtaining data on the structure of various heterogeneous materials. This task is especially important in the areas of clinical medicine and biology. The paper presents the method aimed at determining mechanical properties of a representative element using computed tomography. Based on the finite element method for a given region, a finite-element ensemble is constructed using the scanning data on a computer tomograph of a real sample. For the obtained sample, numerical experiments are performed in the kinematic formulation, after which the problem of the stress-strain state is solved. The stresses obtained as a result of the calculations are averaged and used to determine the components of the elastic constant tensor. Thus, the anisotropic properties of the representative element are determined. To determine the orthotropic properties of the representative element, a target function is introduced, the arguments of which are unknown directions of orthotropy. These unknown directions are determined from the condition of minimizing the objective function. The transformation of the rotation to an anisotropic matrix of elastic constants makes it possible to determine the components of the elastic constant tensor in the orthotropic axes. As an illustration of the technique, calculations of a porous sample are given in the paper, and the obtained results are evaluated. For the quantitative comparison, the invariant of the stress tensor is used. The obtained results illustrate not only a sufficient accuracy of describing the medium in terms of continuity, but also a discrepancy in the results in the case of large porosity.

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

  1. Gross T., Pahr D.H., Zysset P.K. Morphology-elasticity relationships using decreasing fabric information of human trabecular bone from three major anatomical locations // Biomech. Model Mechanobiol. - 2013. - Vol. 12. - P. 793-800.
  2. Schwen L.O., Wolfram U., Rumpf M. Determining effective elasticity parameters of microstructure materials // 15th Workshop on the Finite Element Method in Biomedical Engineering, Biomechanics and Related Fields. - July 2008. - P. 41-26.
  3. Cowin S.C. Continuum Mechanics of Anisotropic Materials. - New York: Springer-Verlag, 2013. - 425 p.
  4. Standardized nomenclature, symbols, and units for bone histomorphometry: a 2012 update of the report of the ASBMR Histomorphometry Nomenclature Committee / D.W. Dempster, J.E. Compston, M.K. Drezner [et al.] // J. Bone Miner Res. - 2013. - Vol. 28. - P. 1-16.
  5. Trabecular bone score (TBS): available knowledge, clinical relevance, and future prospects / V. Bousson, C. Bergot, B. Sutter, P. Levitz, B. Cortet // Osteoporos Int. - 2012. - Vol. 23. - No. 5. - P. 1489-1501.
  6. The predictive value of trabecular bone score (TBS) on whole lumbar vertebrae mechanics: an ex vivo study / J.P. Roux, J. Wegrzyn, S. Boutroy, M.L. Bouxsein, D. Hans, R. Chapurlat // Osteoporos Int. - 2013. - Vol. 24. - No. 9. - P. 2455-2460.
  7. Added value of trabecular bone score over bone mineral density for identification of vertebral fractures in patients with areal bone mineral density in the non-osteoporotic range / K. Nassar, S. Paternotte, S. Kolta, J. Fechtenbaum, C. Roux, K. Briot // Osteoporos Int. - 2014. - Vol. 25. - No. 1. - P. 243-249.
  8. Spine bone texture assessed by trabecular bone score (TBS) predicts osteoporotic fractures in men: The Manitoba Bone Density Program / W.D. Leslie, B. Aubry-Rozier, L.M. Lix [et al.] // Bone. - 2014. - Vol. 67. - P. 10-14.
  9. Коноплев Ю.Г., Митряйкин В.И., Саченков О.А. Применение математического моделирования при планировании операции по эндопротезированию тазобедренного сустава // Учен. зап. Казан. ун-та. Сер.: Физико-математические науки. - 2011. - Т. 153, № 4. - С. 76-83.
  10. Численное исследование влияния степени недопокрытия вертлужного компонента на несущую способность эндопротеза тазобедренного сустава / Ю.Г. Коноплев, А.В. Мазуренко, О.А. Саченков, Р.М. Тихилов // Рос. журн. биомеханики. - 2015. - Т. 19, № 4. - С. 330-343.
  11. Численное исследование напряженно-деформированного состояния тазобедренного сустава при ротационной остеотомии проксимального участка бедренной кости / О.А. Саченков, Р.Ф. Хасанов, П.С. Андреев, Ю.Г. Коноплев // Рос. журн. биомеханики. - 2016. - Т. 20, № 3. - С. 257-271.
  12. Экспериментальное определение тензора структуры трабекулярной костной ткани / А.А. Киченко, В.М. Тверье, Ю.И. Няшин, А.А. Заборских // Рос. журн. биомеханики. - 2011. - Т. 15, № 4. - С. 78-93.
  13. Киченко А.А., Тверье В.М., Няшин Ю.И. Математическое описание поведения губчатой костной ткани под нагрузкой // Математическое моделирование в естественных науках. - 2013. - № 1. - С. 84-85.
  14. Дядюкина А.Д., Киченко А.А. Математическое моделирование трабекулярной костной ткани // Математическое моделирование в естественных науках. - 2016. - Т. 1. - С. 627-630.
  15. Direct Mechanics assessment of elastic symmetries and properties of trabecular bone architecture / R. Huiskes, B. Van Rietbergen, A. Odgaard, J. Kabel // J. Biomech. - 1996. - Vol. 29(12). - P. 1653-1657.
  16. Bone Volume Fraction and Fabric Anisotropy Are Better Determinants of Trabecular Bone Stiffness Than Other Morphological Variables / G. Maquer, S.N. Musy, J. Wandel, T. Gross, P.K. Zysset // Journal of Bone and Mineral Research. - 2015. - Vol. 30. - No. 6. - P. 1000-1008. doi: 10.1002/jbmr.2437
  17. Guidelines for assessment of bone microstructure in rodents using micro-computed tomography / M.L. Bouxsein, S.K. Boyd, B.A. Christiansen, R.E. Guldberg, K.J. Jepsen, R. Müller // J. Bone Miner Res. - 2010. - Vol. 25. - P. 1468-1486.
  18. Zysset P.K., Curnier A. An alternative model for anisotropic elasticity based on fabric tensors // Mechanics of Materials. - 1995. - Vol. 21. - Iss. 4. - P. 243-250.
  19. Zysset P.K. A review of morphology-elasticity relationships in human trabecular bone: theories and experiments // J. of Biomechanics. - 2003. - Vol. 36. - P. 1469-1485.
  20. Gross T., Pahr D.H., Peyrin F., Zysset P.K. Mineral heterogeneity has a minor influence on the apparent elastic properties of human cancellous bone: a SRmCT-based finite element study // Computer Methods in Biomechanics and Biomedical Engineering. - 2012. - Vol. 15. - No. 11. - P. 1137-1144.
  21. Fosbinder R., Orth D. Essentials of Radiologic Science. - Lippincott Williams & Wilkins, 2011. - 392 p.
  22. Otsu N. A threshold selection method from gray-level histograms // IEEE Transactions on Systems, Man, and Cybernetics. - 1979. - Vol. 9. - Iss. 1. - P. 62-66.
  23. Kazembakhshi S, Luo Y. Constructing anisotropic finite element model of bone from computed tomography (CT) // Biomed Mater. Eng. - 2014. - Vol. 24. - No. 6. - P. 2619-26.
  24. Comparison of mixed and kinematic uniform boundary conditions in homogenized elasticity of femoral trabecular bone using micro finite element analyses / J. Panyasantisuk, D.H. Pahr, T. Gross, P.K. Zysset // J. Biomech Eng. - 2015. - Vol. 137. - No. 1. - 011002. doi: 10.1115/1.4028968
  25. Pahr D.H., Zysset P.K. A comparison of enhanced continuum FE with micro FE models of human vertebral bodies // Journal of Biomechanics. - 2009. - Vol. 42. - P. 455-462.
  26. Dependence of mechanical properties of trabecular bone on plate-rod microstructure determined by individual trabecula segmentation (ITS) / B. Zhou, X.S. Liu, J. Wang, X.L. Lu, A.J. Fields, X.E. Guo // J. Biomech. - 2014. - Vol. 47. - No. 3. - P. 702-708.
  27. HR-pQCT-based homogenised finite element models provide quantitative predictions of experimental vertebral body stiffness and strength with the same accuracy as mFE models / D.H. Pahr, E. Dall’Ara, P. Varga, P.K. Zysset // Comput Methods Biomech Biomed Engin. - 2012. - Vol. 15. - No. 7. - P. 711-720.
  28. Assessment of transverse isotropy in clinical-level CT images of trabecular bone using the gradient structure tensor / D. Larsson, B. Luisier, M.E. Kersh [et al.] // Ann Biomed Eng. - 2014. - Vol. 42. - No. 5. - P. 950-959.
  29. Hazrati Marangalou J., Ito J., van Rietbergen K. A novel approach to estimate trabecular bone anisotropy from stress tensors. // Biomech. Model Mechanobiol. - 2015. - Vol. 14. - No. 1. - P. 39-48. doi: 10.1007/s10237-014-0584-6.
  30. Mirkhalaf S.M., Andrade Pires F.M., Ricardo Simoes. Determination of the size of the Representative Volume Element (RVE) for the simulation of heterogeneous polymers at finite strains // Finite Elements in Analysis and Design. - 2016. - Vol. 119. - No. 15. - P. 30-44. DOI.org/10.1016/j.finel.2016.05.004
  31. Лехницкий С.Г. Теория упругости анизотропного тела. - М.: Наука, 1977. - 416 с.
  32. Gerasimov O., Koroleva E., Sachenkov O. Experimental study of evaluation of mechanical parameters of heterogeneous porous structure // IOP Conference Series: Materials Science and Engineering. - 2017. - Vol. 208. - 012013. doi: 10.1088/1757-899X/208/1/012013

Statistics

Views

Abstract - 80

PDF (Russian) - 38

Cited-By


PlumX


Copyright (c) 2018 Kharin N.V., Vorobyev O.V., Berezhnoi D.V., Sachenkov O.A.

Creative Commons License
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.

This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies