CC BY
48
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ВИБРАЦИОННЫХ ПОТОКОВ ЖИДКОСТИ В СИСТЕМЕ ПОР БОЛЬШЕБЕРЦОВОЙ КОСТИ
Л. Б. Маслов1, Д. Г. Арсеньев2, А. В. Зинковский3
1. Ивановский государственный энергетический университет,
канд. техн. наук, доцент, leonid.maslov@tipm.ispu.ru, leonid-maslov@mail.ru
2. С.-Петербургский государственный политехнический университет, д-р техн. наук, профессор, imop@imop.spbstu.ru
3. С.-Петербургский государственный политехнический университет, д-р биол. наук, профессор, bioval@imop.spbstu.ru
1. Введение
Закон функциональной адаптации костной ткани к меняющимся внешним силовым условиям, известный как закон Вульфа [1], является основополагающей гипотезой различных теорий костной перестройки. Ключевым нераскрытым фактором к пониманию механобиологических процессов костной перестройки остается вопрос, как специфические клетки кости, ответственные за процессы резорбции и формирования новой ткани, воспринимают механический стимул и реагируют на него, чтобы привести физиологическое состояние костного вещества в соответствие изменившимся внешним условиям. Поскольку основная масса костных клеток располагается на стенках канальцев и поверхностях лагун, входящих в систему микропор кости, предполагается, что возмущения, вносимые внешней механической нагрузкой в установившееся движение жидкости в транспортной системе кости, могут обеспечивать передачу управляющих сигналов между клетками кости в процессе ее функциональной адаптации [2].
Определенный вклад в решение указанных проблем могут внести методы математического моделирования, основанные на законах механики многофазных сред. Авторы статьи [3] используют коммерческий программный комплекс ABAQUS (ABAQUS Inc.) и линейную изотропную модель пороупругого материала [4] для расчета квазистати-ческого напряженно-деформированного состояния призматического образца компактной ткани под действием гармонической растягивающей силы и изгибающего момента. Показано, что распределение давления в однородной лакунарно-каналикулярной сети в поперечном сечении образца существенно зависит от частоты воздействия и коэффициента просачивания. Задача расчета напряженного состояния модели остеона под действием нагрузки, типичной для длинных трубчатых костей, рассмотрена в работе [5]. Численными расчетами доказано, что нормальная физиологическая нагрузка может приводить к повышению давления жидкости в лакунарно-каналикулярной системе пор остеона в несколько десятков раз большему, чем вызванное изменением кровяного давления. Пороупругая модель ячейки компактной ткани, образованной Гаверсовым каналом и окружающим матриксом, и гармонические колебания призматического образца кортикальной кости рассмотрены в работе [6]. Показано, что амплитуды сдвиговых напряжений и давления жидкости в Гаверсовых каналах возрастают с увеличением частоты физиологически реалистичной нагрузки и могут достигать уровня нескольких паскалей, которого, как показано в других работах, достаточно для генера-
© Л. Б. Маслов, Д. Г. Арсеньев, А. В. Зинковский, 2009
ции морфологического отклика клеток кости. В рассмотренных работах используются модели простых геометрических форм, представляющие собой только один вид костного вещества — компактного или губчатого. Однако экспериментально было показано, что потоки жидкости в компактном веществе кортикального слоя длинных трубчатых костей обусловлены не только деформацией твердого матрикса, но и избыточным ин-трамедуллярным давлением, возникающим в результате нагружения кости [7]. Таким образом, необходимо рассматривать более сложные многофазные модели, включающие в себя основные типы костного вещества, формирующие кость как биомеханическую структуру.
2. Математическая модель костной ткани
Согласно теории пороупругости Био предполагается, что в пористом материале могут быть выделены твердая фаза, представляющая собой упругий формообразующий скелет и воспринимающая основную силовую нагрузку, и жидкая фаза, полностью или частично заполняющая поры.
Определяющие соотношения двухфазного материала, сформулированные Био на основе феноменологического подхода механики сплошных сред [8], в случае анизотропии упругих и гидравлических свойств могут быть записаны следующим образом:
ст(1)(и, и) = С ••е(и) + де(И), (1а)
а(2)(и, И) ••Е = д ••е(и)+ Яе(И), (1б)
где и, И — вектора перемещений твердой и жидкой фаз; е —тензор малых упругих деформаций твердой фазы; е — объемная деформация жидкой фазы; <г(1) ,а(2) —тензоры условных напряжений в твердой и жидкой фазах; С — тензор четвертого ранга эффективных упругих модулей твердой фазы; д — тензор коэффициентов, определяющих взаимное влияние деформаций одной из фаз на возникающие напряжения в другой фазе; Е — единичный тензор второго ранга; Я — гидростатическая константа, которую можно интерпретировать как эффективный модуль объемного сжатия жидкой фазы.
Уравнения (1а, б) связывают между собой эффективные напряжения в твердом
каркасе и жидкости с деформациями твердой и в жидкой фаз с помощью эффектив-
ных модулей. В соответствии с механическим смыслом объемной упругой деформации объемная деформация жидкой фазы равна изменению объема жидкости в порах по отношению к первоначальному объему, т. е. характеризует изменение доли жидкости в представительном элементарном объеме гетерогенного материала. Давление в порах р с точностью до пористости и в соответствии с известным правилом знаков теории упругости является условным гидростатическим давлением, возникающем в жидкой фазе:
-РР = + ^22) + (2)
где р —давление в порах, ^ — пористость материала.
Вывод уравнений динамики пороупругого эффективного континуума при использовании смешанной и - Р формулировки приводит к следующим определяющим соотношениям пороупругого континуума:
<г(и,р) = С • •е(и) - Ар = Сте(и) - Ар, (За)
С(и,р) = А • -е(и) + ц)2Я-1р, (Зб)
где <г(и,р) = С4 • -е(и) —тензор напряжений в точках твердой фазы, вызываемый только деформациями упругой матрицы материала; С — тензор четвертого ранга упругих модулей эффективной пороупругой среды в дренированном состоянии; А — тензор Био, учитывающий вклад напряжений жидкой фазы в выражение полного тензора напряжений двухфазной среды; £ — скалярная переменная, численно равная относительному изменению содержания жидкости в порах.
Тензор коэффициентов эффективных напряжений Био А и тензор упругих модулей среды в дренированном состоянии С выражаются через исходные материальные константы двухфазной среды С^, д^ и Я следующим образом:
Математическая модель биологических структур как пороупругого материала, насыщенного жидкостью, и численный алгоритм реализованы в авторском конечноэлементном комплексе МесЬапісвЕЕ [9]. Система конечно-элементных уравнений, описывающих установившиеся гармонические колебания эффективной среды на основании определяющих уравнений пороупругого тела (3а) и (3б), имеет следующий вид:
где M, K — стандартные глобальные матрицы массы и жесткости; L — дополнительная глобальная матрица массы; Hi, H2 —глобальные матрицы взаимного влияния распределения давления жидкости в порах и деформации упругого скелетона; D и G — глобальные матрицы насыщения и проницаемости; U, P — глобальные векторы комплексных амплитудных значений перемещений и давлений в узлах конечно-элементной сетки; F — глобальный узловой вектор объемных и поверхностных сил.
3. Расчет вынужденных колебаний большеберцовой кости
Рассмотренная модель пороупругой среды была использована для расчета вынужденных колебаний и давления в порах изолированного образца большеберцовой кости голени человека под действием поперечной силы (Fx = Fy = 1 H), изменяющейся по гармоническому закону. Модель большеберцовой кости среднестатистического человека была разработана на основе фотографических снимков поперечных срезов, представленных в компьютерной базе данных проекта Visible Human Национального Института Здоровья США [10]. Граничные условия, наложенные на перемещения упругого скелетона, моделируют одноосное шарнирное соединение большеберцовой кости в коленном и голеностопном суставах. Внешняя поверхность кости представляет собой непроницаемую стенку, что означает отсутствие заданных потоков жидкости через поверхность. Значения материальных констант тканей, полученные расчетным путем согласно алгоритму вычисления эффективных модулей композитных материалов [11], представлены в таблице 1.
Согласно вибрационному анализу имеют место основные изгибные формы колебаний большеберцовой кости преимущественно в сагиттальной (YZ) и фронтальной (XZ) плоскостях (рис. 1, 2). Резонансные частоты, соответствующие изгибным модам, равны 201 Гц и 279 Гц для первой моды в плоскости (YZ) и (XZ), и 713 Гц и 831 Гц —
(4)
Cd = C - R-1QQ.
(5)
(K - w2(M - iwL)) U - (H1 + iwH2) P = F, -(H1 + iwH2)T U + (-D + iw-1 G) P = 0,
(б)
(7)
Таблица 1. Расчетные значения эффективных пороупругих характеристик структурных элементов большеберцовой кости
Струк- турный элемент кости Пори- стость Эффективные характеристики пороупругой среды
Плот- ность Модули упругости среды в дренированном состоянии Тензор Био Гидр. конст.
1-І V — у — С-у г С-ху г — г/у г 1'г у Ахх — Ауу я
% кг/м15 МПа МПа МПа МПа МПа
компактная кость 10 1469 11200 21100 6190 4050 0,224 0,331 0,567 0,494 169
губчатая кость 80 1104 434 1510 333 175 0,104 0,241 0,984 0,961 1790
костный мозг 96 1021 26,7 26,7 11,1 11,1 0,202 0,202 0,999 0,999 2200
Рис. 1. Амплитудно-частотные характеристики модели большеберцовой кости, рассчитанные с помощью авторского программного комплекса МесЬапісвРЕ в трех узлах среднего поперечного сечения вдоль оси У, относящихся к различным видам
тканей:------АЧХ в точке компактного вещества;------— АЧХ в точке губчатого
вещества;.....— АЧХ в точке костного мозга, заполняющего внутренний канал.
(д)
(е)
Рис. 2. Резонансные формы колебаний модели кости и распределение амплитуд давления (кПа) по поверхности, рассчитанные с помощью авторского программного комплекса МесЬашсвРЕ: (а), (б) — первые изгибные формы в сагиттальной (201 Гц) и фронтальной (279 Гц) плоскостях; (г), (д) —вторые изгибные формы в сагиттальной (713 Гц) и фронтальной (831 Гц) плоскостях; (в), (е)—смешанные изгибные формы с кручением вокруг продольной оси (582 Гц) и продольной деформацией (952 Гц) соответственно.
(д) (е)
Рис. 3. Вектора амплитудных значений потоков жидкости (мм/с) в порах компактной ткани большеберцовой кости, рассчитанные с помощью авторского программного комплекса МесЬашсвРЕ на резонансных формах колебаний: (а), (б) —первые изгибные формы в сагиттальной (201 Гц) и фронтальной (279 Гц) плоскостях; (г), (д) —вторые изгибные формы в сагиттальной (713 Гц) и фронтальной (831 Гц) плоскостях; (в), (е) — смешанные изгибные формы с кручением вокруг продольной оси (582 Гц) и продольной деформацией (952 Гц) соответственно.
для второй моды. Кроме того, в исследуемый диапазон попали две смешанные формы колебаний на частотах 582 и 956 Гц. Резонансной частоте 582 Гц соответствует изгиб с кручением относительно продольной оси, а частоте 956 Гц — изгиб с продольной деформацией.
Модальный расчет кости, проведенный в классической упругой постановке с помощью коммерческого конечно-элементного комплекса ANSYS (ANSYS Inc.), показал наличие собственных частот и форм колебаний, в целом соответствующих найденным резонансным частотам и формам колебаний кости в пороупругой постановке. Однако постановка при использовании тензора упругих модулей в дренированном состоянии дает заниженные значения собственных частот по сравнению с разработанной моделью двухфазной среды на 1%, 1.6, 14.6, 3.2, 1.7 и 11.7% соответственно. Упругая постановка при использовании значений модулей в недренированном состоянии дает завышенные значения собственных частот на 0.8%, 3.1, 6.7, 2.0, 5.1 и -0.1% соответственно. Таким образом, за исключением последней продольно-изгибной моды колебаний пороупругий алгоритм дает значения частот, лежащих в диапазоне между двумя упругими постановками.
Интересно отметить, что первые формы колебаний кости имеют большие амплитуды колебаний в поперечной плоскости, что может быть использовано при разработке методов резонансной вибрационной диагностики повреждений опорно-двигательного аппарата человека, т. к. это облегчает процесс измерения колебаний. Однако для детальной разработки алгоритма диагностики необходимо также рассмотреть чувствительность резонансных частот к механобиологическому состоянию кости в зоне перелома. С точки зрения вопроса восстановления костной ткани в зоне перелома, более важным представляется характер распределения давления (рис. 2) и вибрационных потоков жидкости в порах в компактной костной ткани (рис. 3). Это связано с тем, что именно благодаря движению внутритканевой жидкости обеспечиваются как доставка питательных веществ и строительных микроэлементов клеткам, так и регуляция работы всей механо-чувствительной системы костных тканей. Графики амплитудных значений вибрационных потоков показывают, что наибольший эффект оказывают вторые изгибные формы колебаний, поскольку потоки в компактной ткани достигают существенных значений при намного меньших перемещениях упругого скелетона на соответствующих резонансных формах колебаний. Отметим также, что старшие, и особенно продольные, формы колебаний не только повышают «коэффициент полезного действия» вибрационной стимуляции, но также снижают возможные болевые ощущения, которые могли бы быть вызваны поперечными смещениями значительной амплитуды.
4. Заключение
Полученные результаты позволяют использовать разработанный подход для расчета сложных биомеханических систем, насыщенных жидкостью, таких как длинные трубчатые кости опорно-двигательного аппарата, связки, хрящи, сухожилия и скелетные мышцы, для математического описания которых аппарат теории двухфазных сред представляется наиболее естественным. Идентификация механического аспекта, регулирующего процесс остеосинтеза, может дать теоретическую базу для разработки биомеханических методов и устройств для лечения остеопороза, интенсификации процесса восстановления кости после перелома, стимуляции врастания костной ткани в материал имплантанта.
1. Wolff J. Das Gesetz der Transformation der Knochen. Berlin: A. Hirchwild (1892). Translated as: The Law of Bone Remodeling / Edited by P. Maquet and R. Furlong. Berlin: Springer-Verlag, 1986.
2. Регирер С. А., Щадрина Н. Х. Движение крови и интерстициальной жидкости в костной ткани: обзор // Изв. РАН: Механика жидкости и газов. 1999. №5. С. 4-28.
3. Manfredini P., Cocchetti G., Maier G., Redaelli A., Montevecchi F. M. Poroelastic finite element analysis of bone specimen under cyclic loading // J. Biomech. 1999. Vol. 32. N2. P. 135144.
4. Rice J. R., Cleary M. P. Some basic stress-diffusion solutions for fluid saturated elastic porous media with compressible constituents // Rev. Geophys. Space Phys. 1976. Vol. 14. P. 227-241.
5. Zhang D., Weinbaum S., Cowin S. C. On the calculation of bone pore water pressure due to mechanical loading // Int. J. Solids Structures. 1998. Vol. 35. N 34-35. P. 4981-4997.
6. Swan C. C., Lakes R. S., Brand R. A., Stewart K. J. Micromechanically based poroelastic modeling of fluid flow in Haversian bone // J. Biomech. Engng. 2003. Vol. 125. N 1. P. 25-37.
7. Qin Y.X., Lin W., Rubin C. T. The pathway of bone fluid flow as defined by in vivo in-tramedullary pressure and streaming potential measurements // Annals Biomed. Engng. 2002. Vol. 30. P. 693-702.
8. Biot M. A. Theory of propagation of elastic waves in a fluid-saturated porous solid, part I: low frequency range // J. Acoust. Soc. Am. 1956. Vol. 28. N 2. P. 168-178.
9. Маслов Л. Б. Применение теории Бийо к исследованию вынужденных колебаний пористых структур // Механика композиционных материалов и конструкций. 2005. Т. 11. №2. С. 276-297.
10. Ackerman M. J. The Visible Human Project // Proceedings of the IEEE, March 1998. Vol. 86. №3. P. 504-511.
11. Арсеньев Д. Г., Зинковский А. В., Маслов Л. Б. Эффективные упругие характеристики анизотропной модели пористого биологического материала, насыщенного жидкостью // Научно-технические ведомости Санкт-Петербургского государственного технического университета. 2008. №3. Т. 59. С. 230-236.
Статья поступила в редакцию 19 марта 2009 г.
