BUILDING THE INHOMOGENEOUS FINITE ELEMENT MODEL BY THE DATA OF COMPUTED TOMOGRAPHY

Abstract


The aim of the work is to reveal a methodology for constructing a finite element model by tomography data. To evaluate the model, calculations of the femur were carried out. The relevance of this study is confirmed by the effect of the distribution of the mechanical properties of the bone on stress-strain state and the need for individualizing the approach to modeling. Numerical studies were performed using the finite element method in the Ansys software, computer tomography data processing was carried out in the Avizo software. The problem of a linear inhomogeneous elastic body was considered. Power functions of the optical density were used to determine the Young's modulus and the limiting voltage, in turn, the optical density was determined from linear relations depending on the Hounsfield numbers. For finite element model, the mechanical properties of the material were distributed for each element according to the tomography data. After solving the problem of the stress-strain state, at each node a factor of safety was determined adjusted for the properties of the material from the tomography data. Inhomogeneous and homogeneous models with average properties were built. Calculations for both models were performed. Numerical results clearly illustrate significant differences in the results of the stress-strain state of the inhomogeneous and homogeneous models of the organ. Inhomogeneous model allows us to evaluate the local strength of bone tissue taking into account individual characteristics.

Full Text

Введение Костная ткань обладает неравномерными механическими свойствами, что является следствием ее адаптационных свойств [3, 8, 19-21]. Снижение физической активности [24, 34], физиологические [11, 12] или возрастные особенности могут приводить к локальной потере прочности костных органов. В клинической практике это может обусловливать расхождение между реальным поведением органа и ожидаемым (моделируемым) результатом оперативного вмешательства. В ортопедии и травматологии множество исследований направлены на улучшение методик оперативного вмешательства [15, 17, 18], оптимизации используемых стержневых аппаратов [2, 4, 5, 7-10], но на практике, ввиду сложности прогнозирования поведения материала костной ткани [27, 30-32], многие решения о тактике операции приходится принимать хирургу на месте. В связи с этим актуальной задачей становится оценка механических свойств органа не в интегральном, но в распределенном смысле. Современные подходы к численному моделированию позволяют решать задачи для сложных областей [6, 13, 14, 16, 23, 25, 36], наряду с развитием компьютерной томографии, благодаря которой можно судить о внутреннем строении органа, задача индивидуализации при моделировании органов из костной ткани может предоставить полезную информацию в предоперационный период. Целью нашего исследования является построение конечно-элементной модели, позволяющей учитывать особенности распределения механических свойств (модуль Юнга и предел прочности) по данным компьютерной томографии, и определение ее напряженно-деформируемого состояния. Материалы и методы Постановка задачи Принятые допущения: при построении модели авторы предполагают, что материал изотропный и негомогенный. Более того, исходя из предположения наличия связи между физической плотностью и механическими характеристиками, будем использовать оптическую плотность для расчета модуля Юнга и предельных напряжений. Математическая модель Расчетная область включает в себя некоторое негомогенное тело. Обозначим его за V, а за ¶V - границу этого множества (свободную поверхность). Механическое поведение системы, занимающей область V в R3 с границей ¶V, в рамках линейной теории упругости описывается следующей системой уравнений: , , (1) , , (2) , , (3) где V = VȶV, - тензор напряжений; - тензор упругой деформации; - вектор перемещений; - тензор упругих свойств, - некоторая функция, определяющая негомогенность среды. Более того, предположим, что допустимые напряжения также зависят от пространственной координаты: , (4) Тогда для оценки напряженно-деформируемого тела введем коэффициент запаса, определяемый по формуле , (5) где - напряжения по Мизесу в заданной точке. В этом случае ключевую роль будут играть зоны, где введённый выше коэффициент (5) превышает единицу. Для построения моделей органов из костной ткани предлагается использовать данные компьютерной томографии, в этом случае для каждой точки можно сопоставить значение компьютерной томографии . Для удобства использования таких данных их принято нормировать (шкала Хаунсфилда): (6) где , - линейные коэффициенты ослабления для воды и воздуха соответственно. Из широкого спектра исследований [1, 3, 22, 26, 29, 33] известно, что существует связь между числами Хаунсфилда и оптической плотностью, упругими константами, предельными напряжениями. В общем виде эти зависимости можно записать так [26, 29, 33]: , (7) , (8) , (9) здесь коэффициенты a и b с соответствующими индексами определяются из эксперимента и могут варьироваться в зависимости от вида костной ткани или органа. Численное решение задачи Численная реализация проводилась на основе метода конечных элементов. Основная идея заключается в построении взаимосвязи между конечно-элементной сеткой и данными компьютерной томографии. Для этого необходимо на дискретизированном множестве, которое представляет собой набор данных узлов (их координаты) и элементов (номера узлов, образующие элемент), определить коэффициенты ослабления в узлах. (10) (11) где - конечно-элементная дискретизация объема ; - номер узла и его координаты; - номер элемента и образующие его номера узлов, - количество узлов и элементов соответственно. Далее проводится осреднение механических характеристик для каждого элемента, что позволяет, находясь в упругой постановке, использовать классические конечные элементы. Для определения коэффициентов запаса (5) создается массив допустимых напряжений, по размерности совпадающий с количеством узлов в сетке. Для обработки данных томографии был использован программный продукт Avizo, для конечно-элементного моделирования - ANSYS. Для конечно-элементного моделирования напряженно-деформированного состояния органа применялся четырехузловой конечный элемент с линейной аппроксимацией solid186. Для осреднения коэффициентов ослабления в элементе использовалось среднее арифметическое, что накладывает условия на форму конечного элемента (равномерность трех линейных размеров). Выражения (7)-(9) удобно выразить сразу через числа Хаунсфилда: (12) (13) Используемые значения коэффициентов в выражениях (12)-(13) приведены в табл. 1 [26, 29, 33, 28]. Таблица 1 Значения параметров Коэффициент Модуль Юнга E 46·103 2 Предельное напряжение 60 1,82 В случае конечно-элементной модели соотношения для определения коэффициента запаса (5) можно переписать в виде (14) Анализ напряженно-деформированного состояния органа в этом случае будем проводить по данным коэффициента запаса (14). Модельная задача Для расчетов была рассмотрена бедренная кость (см. рис. 1, а). Изучались следующие задачи: 1) одноосное сжатие, 2) изгиб и 3) совместный изгиб со сжатием. В табл. 2 приведены величины сил. Для оценки значимости предлагаемого подхода были рассмотрены две модели: первая негомогенная - напряженно-деформированное состояние определяется по описанной методике, вторая гомогенная - задача напряженно-деформированного состояния решается для некоторого осреднения механических свойств по всему объему: E = 32 ГПа, σ = 48 МПа. Таблица 2 Величины сил Расчетный случай F1, Н F2, Н 1 20 0 2 0 200 3 20 200 Образец бедренной кости подвергался сканированию на компьютерном томографе Vatech PaX-I 3D. Данные обрабатывались в программном комплексе Avizo. Порог экстерьера определялся по методу Отсу, методика обработки подробней изложена в работе [35]. Далее для множества точек, относящегося к органу, строилась поверхностная сетка и после этого объемная сетка. Для удовлетворения равномерности размеров элементов строилась довольно мелкая сетка (порядка 14·104 узлов и 7·105 элементов). Данные о сетке и значениях чисел Хаунсфилда в узлах экспортировались в файл, который загружался в программный продукт ANSYS. Там проводилось распределение осредненных (по элементу) механических характеристик, накладывались кинематические и статические граничные условия, производились решение а б D B A F F C в Рис. 1. Рассматриваемый образец: а - томография бедренной кости; б - диафизарный участок бедренной кости; в - расчетная схема и расчет коэффициента запаса (14). В расчетах рассматривался диафизарный участок бедренной кости (см. рис. 1, б). На рис. 1, в приведена схема граничных условий: к поперечному сечению AD прикладывалась нагружающая сила, на поперечном сечении BC фиксировались все перемещения узлов, величины сил приведены в табл. 2. Далее в программном комплексе ANSYS решалась задача напряженно-деформированного состояния и определялись коэффициенты запаса по формуле (12). Результаты и обсуждение Для построенной конечно-элементной сетки производилась ее загрузка в программный продукт ANSYS. Каждому элементу назначался параметр, определяющий значение модуля Юнга и предела прочности согласно (12), (13) (на рис. 2 приведена сетка с цветовой диаграммой, указывающей распределение модуля Юнга и предельных напряжений по объему). На рис. 2 можно видеть, что наибольшей жесткостью и прочностью обладает кортикальная кость (красные области на рисунке). Стоит прокомментировать синие вкрапления на внутренней поверхности кости, которым соответствуют малые модуль Юнга и предел прочности, это вызвано наличием пор, которые при сканировании ослабляют интегральные значения коэффициента ослабления в малой области. После получения решения анализировались поля распределений напряжений по Мизесу и коэффициента запаса, рассчитанного по соотношению (14). На ряду с негомогенной задачей рассматривалась связанная гомогенная задача (обладающая той же геометрией, теми же граничными условиями, но с осреднёнными по всему объему механическими характеристиками). σ 16 20 24 28 33 38 43 48 54 60 МПа E 0,80 2 4 7 10 14 19 24 30 36 ГПа Рис. 2. Распределение механических параметров по конечно-элементной модели: модуля Юнга и предельных напряжений Расчетный случай 1 Для случая сжатия негомогенного образца локализация максимальных напряжений по Мизесу происходила в середине кортикала кости на расстоянии примерно трети длины образца от заделки (см. рис. 3, а). Для гомогенного образца локализация максимальных напряжений по Мизесу происходила на внешней поверхности кортикала и незначительно перемещалась в продольном направлении к нагруженному торцу (см. рис. 4, а). Численные значения максимальных напряжений в обоих случаях отличались незначительно: 36 и 37 МПа для негомогенного и гомогенного образцов соответственно. Картины полей коэффициента запаса отличались существенно. а б в Рис. 3. Результаты для первого расчетного случая негомогенного образца: а - распределение напряжений по Мизесу, МПа; б, в - распределение коэффициента запаса а б в Рис. 4. Результаты для первого расчетного случая гомогенного образца: а - распределение напряжений по Мизесу, МПа; б, в - распределение коэффициента запаса В случае негомогенного образца коэффициент запаса превышал единицу на внешней и внутренней поверхностях кости, в области локализации максимальных напряжений по Мизесу, при этом эта область протяжённая: почти половина образца (рис. 3, б, в), наибольшее значение коэффициента запаса в этом случае составило 1,38. Для гомогенного образца наибольшее значение коэффициента запаса достигалось в той же области, где локализовались максимальные напряжения по Мизесу, наибольшее значение коэффициента запаса в этом случае составило 0,89 (см. рис. 4, б, в). То есть согласно негомогенной модели образец будет претерпевать разрушение, в то время как согласно гомогенной постановке критических значений напряжения не достигают. Разница в коэффициентах запаса составила 55%. Расчетный случай 2 Для случая изгиба негомогенного образца локализация максимальных напряжений по Мизесу происходила в середине кортикала кости на малом удалении от заделки (рис. 5, а). Для гомогенного образца локализация максимальных напряжений по Мизесу происходила на внешней поверхности кортикала и примерно на таком же расстоянии от заделки, как и в негомогенном случае (рис. 6, а). Численные значения максимальных напряжений в обоих случаях отличались незначительно: 24 и 23 МПа для негомогенного и гомогенного образцов соответственно. Различия картин распределения коэффициента запаса схожи с первым расчетным случаем. В случае негомогенного образца коэффициент запаса был близок к критическому значению (единице) на внешней и внутренней поверхностях кости, в области локализации максимальных напряжений по Мизесу, при этом эта область значительна (см. рис. 5, б, в), наибольшее значение коэффициента запаса в этом случае составило 0,95. Для гомогенного образца наибольшее значение коэффициента запаса достигалось в той же области, где и локализовались максимальные напряжения по Мизесу, наибольшее значение коэффициента запаса в этом случае составило 0,58 (см. рис. 6, б, в). То есть согласно негомогенной модели образец будет близок к разрушению, в то время как согласно гомогенной постановке критических значений напряжения не достигают с достаточным запасом. Разница в коэффициентах запаса составила 63%. Расчетный случай 3 Для совместного случая нагружения негомогенного образца локализация максимальных напряжений по Мизесу происходила в середине кортикала кости на расстоянии примерно трети длины образца от заделки (рис. 7, а). Для гомогенного образца локализация максимальных напряжений по Мизесу происходила на внешней поверхности кортикала и незначительно перемещалась в продольном направлении к нагруженному торцу (см. рис. 8, а). Численные значения максимальных напряжений в обоих случаях отличались незначительно: 57 и 53 МПа для негомогенного и гомогенного образцов соответственно. Картины полей коэффициента запаса отличались существенно. В случае негомогенного образца коэффициент запаса превышал единицу на внешней и внутренней поверхностях кости, в области локализации максимальных напряжений по Мизесу, при этом эта область протяжённая: почти половина образца и распространяется по окружному направлению (рис. 7, б, в), наибольшее значение коэффициента запаса в этом случае составило 2,33. а б в Рис. 5. Результаты для второго расчетного случая негомогенного образца: а - распределение напряжений по Мизесу, МПа; б, в - распределение коэффициента запаса а б в Рис. 6. Результаты для второго расчетного случая гомогенного образца: а - распределение напряжений по Мизесу, МПа; б, в - распределение коэффициента запаса а б в Рис. 7. Результаты для третьего расчетного случая негомогенного образца: а - распределение напряжений по Мизесу, МПа; б, в - распределение коэффициента запаса а б в Рис. 8. Результаты для третьего расчетного случая гомогенного образца: а - распределение напряжений по Мизесу, МПа; б, в - распределение коэффициента запаса Для гомогенного образца наибольшее значение коэффициента запаса достигалось в той же области, где и локализовались максимальные напряжения по Мизесу - на внешней поверхности кортикала, наибольшее значение коэффициента запаса в этом случае составило 1,28 (рис. 8, б, в). Разница в коэффициентах запаса - 82%. То есть при совместном нагружении образец будет разрушаться при много меньших силах, согласно расчету негомогенной модели, чем те, которые рассчитываются при гомогенной модели. Эти примеры иллюстрируют необходимость учета распределения механических свойств костного органа при моделировании напряженно-деформированного состояния. Качественным образом изменяются картины распределения коэффициентов запаса, значительно изменяются их количественные значения. Приведенная методика позволяют моделировать механическое поведение костных органов с учетом индивидуализации не только геометрии органа, но и его механических свойств. Выводы 1. В работе приведена методика для построения конечно-элементной негомогенной модели органа из костной ткани по данным компьютерной томографии. 2. Численные результаты модельных задачах иллюстрируют значимые различия в результатах напряженно-деформированного состояния органа и позволяют судить о локальной прочности костной ткани. 3. При моделировании принимались следующие допущения: негомогенный материал изотропный, существует связь между физической плотностью и механическими характеристиками, что позволяет использовать оптическую плотность для расчета модуля Юнга и предельных напряжений. 4. Применение описанной методики позволяет на основании анализа результатов коэффициентов запаса делать заключение о напряженно-деформированном состоянии объекта при различных внешних нагрузках. 5. Предлагаемая модель может быть расширена с учетом ортотропности материала.

About the authors

O. A Sachenkov

O. V Gerasimov

E. V Koroleva

D. A Mukhin

V. V Yaikova

I. F Akhtyamov

F. V Shakirova

D. A Korobeynikova

H. Ch Khan

References

  1. Акулич А.Ю., Акулич Ю.В., Денисов А.С. Экспериментальное определение разрушающих касательных напряжений трабекулярной костной ткани головки бедра человека // Российский журнал биомеханики. - 2010. - Т. 14, № 4. - С. 7-16.
  2. Акулич Ю.В., Акулич А.Ю., Денисов А.С. Влияние количества и размеров резьбовых фиксаторов на адаптационные изменения механических свойств губчатой костной ткани и усилие сжатия отломков после контролируемого остеосинтеза // Российский журнал биомеханики. - 2012. - Т. 16, № 2. - С. 21-29.
  3. Акулич Ю.В., Акулич А.Ю., Денисов А.С., Шайманов П.С., Шулятьев А.Ф. Уточнение индивидуальной зависимости модуля упругости трабекулярной костной ткани от объемного содержания матрикса // Российский журнал биомеханики. - 2014. - Т. 18, № 2. - С. 158-167.
  4. Акулич Ю.В., Подгаец Р.М., Скрябин В.Л., Сотин А.В. Исследование напряженно-деформированного состояния эндопротезированного тазобедренного сустава // Российский журнал биомеханики. - 2007. - Т. 11, № 4. - С. 9-35.
  5. Анисимов О.Г., Ахтямов И.Ф., Шигаев Е.С. Особенности стационарного этапа лечения переломов проксимального отдела бедренной кости. - Казань: ТаГраф, 2017.
  6. Андреев П.С., Коноплев Ю.Г., Саченков О.А., Хасанов Р.Ф., Яшина И.В. Математическое моделирование ротационной флексионной остеотомии // Научно-технический вестник Поволжья. - 2014. - № 5. - С. 18-21.
  7. Ахтямов И.Ф., Коваленко А.Н., Анисимов О.Г., Закиров Р.Х. Лечение остеонекроза головки бедра. - Скрипта, 2013. - 176 с.
  8. Ахтямов И.Ф., Шакирова Ф.В., Клюшкина Ю.А., Бакланова Д.А., Гатина Э.Б., Алиев Э.О. Анализ регенеративного процесса в области перелома большеберцовой кости // Травматология и ортопедия России. - 2016. - № 1 (79). - С. 100-108.
  9. Ахтямов И.Ф., Хаертдинов И.С., Шигаев Е.С., Коваленко А.Н., Гатина Э.Б. Лечение пострадавших с переломами проксимального отдела бедренной кости в условиях Больницы скорой медицинской помощи // Современное искусство медицины. - 2013. - № 1 (9). - С. 23-30.
  10. Ахтямов И.Ф., Шигаев Е.С., Анисимов О.Г. Особенности стационарного этапа лечения переломов проксимального отдела бедренной кости. - Казань: ТаГраф, 2017. - 222 с.
  11. Банецкий М.В. Биомеханическое обоснование использования вертлужного компонента при эндопротезировании тазобедренного сустава: дис. … канд. мед. наук. - М., 2008 - 94 с.
  12. Белецкий А.В., Ахтямов И.Ф., Богосьян А.Б., Герасименко М.А. Асептический некроз головки бедренной кости у детей. - Скрипта, 2010. - 255 с.
  13. Бережной Д.В., Сагдатуллин М.К., Султанов Л.У. Расчет взаимодействия деформируемых конструкций с учетом трения в зоне контакта на основе метода конечных элементов // Вестник Казанского технологического университета. - 2014. - Т. 17, № 14. - С. 478-481.
  14. Иванов Д.В., Барабаш А.П., Барабаш Ю.А. Интрамедуллярный стержень нового типа для остеосинтеза диафизарных переломов бедра // Российский журнал биомеханики - 2015. - Т. 19, № 1. - С. 52-64.
  15. Измайлова З.Т. Предоперационная диагностика модульной трансформации при чрескостном остеосинтезе бедренной кости // Российский журнал биомеханики. - 2009. - Т. 13, № 2. - C. 93-98.
  16. Саченков О.А., Хасанов Р.Ф., Андреев П.С., Коноплев Ю.Г. Численное исследование напряженно-деформированного состояния тазобедренного сустава при ротационной остеотомии проксимального участка бедренной кости // Российский журнал биомеханики. - 2016. - Т. 20, № 3. - C. 257-271.
  17. Климов О.В. Расчет и контроль биомеханической оси нижней конечности во фронтальной плоскости при ее коррекции по Илизарову // Российский журнал биомеханики - 2014. - Т.18, № 2. - С. 239-247.
  18. Менщикова Т.И., Долганова Т.И., Аранович А.М. Влияние силы мышц бедра и голени на опорные реакции стоп у больных ахондроплазией после коррекции роста // Российский журнал биомеханики. - 2014. - Т. 18, № 2. - С. 247-258.
  19. Тверье В.М., Никитин В.Н. Задача коррекции прикуса в зубочелюстной системе человека // Российский журнал биомеханики - 2015. - Т. 19, № 4. - С. 344-358.
  20. Тверье В.М., Няшин Ю.И., Никитин В.Н. Биомеханическая модель определения усилий мышц и связок в зубочелюстной системе человека // Российский журнал биомеханики - 2013. - Т. 17, № 2. - С. 8-20.
  21. Шилько С.В., Черноус Д.А., Бондаренко K.К. Метод определения in vivo вязкоупругих характеристик скелетных мышц // Российский журнал биомеханики - 2007. - Т. 11, № 1. - С. 45-54.
  22. Cuppone M., Seedhom B.B., Berry E., Ostell A.E. The longitudinal Young’s modulus of cortical bone in the midshaft of human femur and its correlation with ct scanning data // Calcified Tissue International. - 2004. - Vol. 74 (3). - P. 302-309.
  23. Davydov R.L., Sultanov L.U. Numerical algorithm for investigating large elasto-plastic deformations // journal of engineering physics and thermophysics. - 2015. - Vol. 88 (5). - P. 1280-1288. doi: 10.1007/s10891-015-1310-7
  24. Baltina T.V., Ahmetov N.F., Sachenkov O.A., Fedyanin A.O., Lavrov I.A. The influence of hindlimb unloading on bone and muscle tissues in rat model // BioNanoSci. - 2017. - № 7(1). - P. 67-69. doi: 10.1007/s12668-016-0288-8
  25. Galiullin R.R., Sachenkov O.A., Khasanov R.F. and Andreev P.S. Evalution of external fixation device stiffness for rotary osteotomy // International Journal of Applied Engineering Research. - 2015. - Vol. 10, № 24. - P. 44855-44860.
  26. Gupta S., Dan P. Bone geometry and mechanical properties of the human scapula using computed tomography data // Trends Biomater. Artif. Organs. - 2004. - Vol. 17(2). - P. 61-70.
  27. Heesakkers N., van Kempen R., Feith R., Hendriks J., Schreurs W. The long-term prognosis of Legg-Calvé-Perthes disease: a historical prospective study with a median follow-up of forty one years // International Orthopaedics. - 2015. - Vol. 39 (5). - P. 859-863.
  28. Isaksson H., Malkiewicz M., Nowak R., Helminen H.J., Jurvelin J.S. Rabbit cortical bone tissue increases its elastic stiffness but becomes less viscoelastic with age // Bone. - 2010. - Vol. 47(6). - P. 1030-1038.
  29. Kaneko T.S., Pejcic M.R., Tehranzadeh J., Keyak J.H. Relationships between material properties and CT scan data of cortical bone with and without metastatic lesions // Med. Eng. Phys. - 2003. - Vol. 25 (6). - P. 445-54.
  30. Nakamura N., Inaba Y., Machida J., Saito T. Rotational open-wedge osteotomy improves treatment outcomes for patients older than eight years with Legg-Calve-Perthes disease in the modified lateral pillar B/C border or C group // International Orthopaedics. - 2015. - Vol. 39 (7). - P. 1359-1364.
  31. Ozel B.D., Ozel D., Ozkan F., Halefoglu A.M. Diffusion-weighted magnetic resonance imaging of femoral head osteonecrosis in two groups of patients: Legg-Perthes-Calve and Avascular necrosis // Radiologia Medica. - 2016. - Vol. 121 (3). - P. 206-213.
  32. Pailhe R., Cavaignac E., Murgier J., de Gauzy J.S., Accadbled F. Triple osteotomy of the pelvis for Legg-Calve-Perthes disease: a mean fifteen year follow-up // International Orthopaedics. - 2016. - Vol. 40. - P. 115-122.
  33. Rho J.Y., Hobatho M.C., Ashman R.B. Relations of mechanical properties to density and CT numbers in human bone // Med. Eng. Phys. - 1995. - Vol. 17(5). - P. 347-55.
  34. Sachenkov O., Kharislamova L., Shamsutdinova N., Kirillova E. and Konoplev Yu. Evaluation of the bone tissue mechanical parameters after induced alimentary Cu-deficiency followed by supplementary injection of Cu nanoparticles in rats // IOP Conference Series: Materials Science and Engineering. - 2015. - Vol. 98. - P. 012-015. doi: 10.1088/1757-899X/98/1/012015
  35. Shigapova F.A., Baltina T.V., Konoplev Y.G., Sachenkov O.A. Methods for automatic processing and analysis of orthotropic biological structures by microscopy and computed// International Journal of Pharmacy & Technology. - 2016. - Vol. 8, № 3. - P. 14953-14964.
  36. Sagdatullin M.K., Berezhnoi D.V. Statement of the problem of numerical modelling of finite deformations // Applied Mathematical Sciences. - 2014. - Vol. 8, № 35. - P. 1731-1738. DOI: 12988/ams.2014.4283

Statistics

Views

Abstract - 53

PDF (Russian) - 56

Refbacks

  • There are currently no refbacks.

Copyright (c) 2022 Российский журнал биомеханики

This website uses cookies

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

About Cookies