Алгоритм численного анализа биологических тканей на основе модели двухфазной среды Текст научной статьи по специальности «Физика»

УДК 620.179.1

АЛГОРИТМ ЧИСЛЕННОГО АНАЛИЗА БИОЛОГИЧЕСКИХ ТКАНЕЙ НА ОСНОВЕ МОДЕЛИ ДВУХФАЗНОЙ СРЕДЫ

Маслов Л.Б. кандидат техн.наук

В статье рассмотрен алгоритм конечно-элементного анализа вынужденных колебаний модели биологической ткани в трехмерной постановке и представлены результаты численного анализа тестовых примеров. Ткань представлена как сплошная двухфазная среда, образованная упругим каркасом и жидкой фазой согласно теории пороупругости Бийо. Разработанный алгоритм расчета вынужденных колебаний тел, имеющих систему связанных пор, заполненных вязкой жидкостью или воздухом, может быть также применен для вибрационного анализа слоистых пористых структур или пенообразных тел и в геомеханике при анализе распространения упругих волн.

Введение

Известно, что твердые и мягкие ткани, формирующие основные элементы опорно-двигательной системы человека, такие как кости, мышцы и связки, могут при определенных ограничениях рассматриваться с позиций механики гетерогенных сред в виде сплошной двухфазной среды [1]. Предполагается, что в ткани могут быть выделены твердая фаза, представляющая собой упругий формообразующий скелет и воспринимающая основную силовую нагрузку, и жидкая фаза, образованная сложной по биологическому составу жидкостью с растворенными в ней питательными веществами. Математические модели, описывающие такую структуру, опираются на два фундаментальных подхода механики. В первом случае, основанном на теории гетерогенных сред, в явной форме записываются законы сохранения и уравнения движения для каждой фазы и всего тела в целом [2]. Во втором случае используется аппарат эффективных характеристик и ткань представляется как сплошная пористая среда, параметры которой выражаются в терминах эффективных напряжений и специфических физико-механических свойств материала [3, 4].

Задача математического моделирования вынужденных колебаний пороупругого тела актуальна для корректной научной интерпретации результатов измерений и теоретического обоснования методов резонансной вибрационной диагностики элементов опорно-двигательного аппарата человека, получивших активное развитие в настоящее время [5, 6]. Однако в известных библиографических источниках, как правило, рассматриваются вопросы распределения давления жидкости в порах при наличии внешнего гармонического воздействия, инерционными же эффектами, учитывающими возможное движение жидкости в порах, чаще всего пренебрегают [4]. Таким образом, в статье будет рассмотрен алгоритм приближенного расчета вынужденных колебаний пороупругого тела в трехмерной постановке и представлены результаты численного анализа.

Метод решения. Уравнения связанной динамической задачи пороупругости, сформулированные в теории Бийо [3, 4], могут быть записаны в следующем тензорном виде в объеме тела V:

^ , & 2» V-о + {V = р—у

£ = (Уи)5

о = 20(£ + 0Еу /(1 - 2у)) - арЕ

др „ 59

а—— V - к\р + а— =у Ы дг

(1)

где а и е - тензора эффективных напряжений и деформаций; 9 - объемная деформация; Е - единичный тензор; где V - дифференциальный оператор Гамильтона; и и ^ - вектора перемещений и заданных объемных сил; р - давление в порах; Y - плотность объемных источников; р, 6 и V - плотность, модуль сдвига и коэффициент Пуассона; к - коэффициент проницаемости, или гидравлическая проводимость; а - коэффициент эффективных напряжений Бийо; а - коэффициент, численно равный отношению количества жидкости, покинувшей или заполнившей поры, к заданному давлению при отсутствии объемной деформации упругого скелета [3]. В известной литературе коэффициент а не имеет собственного названия и, как правило, записывается через другие физико-механические константы.

Согласно [4], параметры к и а имеют следующий вид:

к =

к_

а = •

а 2(1 - 2у)(1 - 2у „)

20 (V н -V)

(2)

где к' - проницаемость - существенная константа, входящая в закон диффузии Дарси; р - вязкость жидкости в порах; Vu - коэффициент Пуассона при условии заполнения пор жидкостью.

Необходимо отметить важную особенность описания двухфазного континуума с помощью аппарата эффективных характеристик. Приведенные выше модули сдвига и коэффициент Пуассона, а также выражаемые через них модуль Юнга Е и объемный модуль упругости К имеют физический смысл упругих параметров, вычисляемых для представительного объема пористого материала, в котором поры не заполнены жидкостью, или осушены. В

англоязычной литературе материальные константы, определенные при таких условиях, обозначаются индексом «d», т.е. «drained». Для упрощения дальнейшего изложения индекс «d» будет опущен. В противоположность «осушенным» параметрам при данном подходе вводится система упругих модулей, вычисляемых для представительного объема пористого материала, в котором поры полностью заполнены жидкостью, не подвижной относительно твердого каркаса. Зарубежные авторы обозначают материальные константы, определенные при таких условиях, индексом «и», т.е. «undrained» [4, 7].

Граничные условия на поверхности тела S могут быть представлены в виде существенных условий, налагаемых на функции перемещений и давлений на части поверхности Si, и в виде естественных условий, налагаемых на поверхностные силы и потоки на части поверхности S2:

u = ~ на S', fs = n • о на S2 , S' + S'2 = S

~ dp (3)

p = p на S[, q = к— на S2 , S[ + S2 = S dn

Записанные выше уравнения представляют собой систему связанных дифференциальных уравнений относительно неизвестных физических полей

- перемещений, деформаций и напряжений формообразующего скелета и давления жидкости в порах.

Для приближенного решения представленной системы уравнений воспользуемся общим подходом теории численных методов - методом взвешенных невязок (МВН) [8]. Рассмотрим первое уравнение системы (1). Предполагая, что искомые функции являются приближенным решением задачи, условие минимизации функции невязки можно записать в виде

j (V- о + fv + fn ) • w dV = 0 (4)

V

где w - весовая функция; fin - вектор сил инерции; V

- объем тела.

Выражение (4) может быть преобразовано к так называемому ослабленному уравнению МВН, являющемуся исходным соотношением для многих численных методов, путем применения известной теоремы Остроградского-Гаусса о дивергенции, широко используемому в теории упругости [9]:

j (V-о) • w dV = j (V-(о • w) - о --Vw) dV=

j fs - w dS + j fv - w dV + j fin - w dV

S2

j n - (о - w) dS - j о - - Vw dV

(5)

5 V

где п - вектор внешней нормали к поверхности тела Б.

Накладывая ограничение на выбор весовой функции, а именно, учитывая, что весовая функция должна удовлетворять кинематическим (существенным) граничным условиям на части поверхности тела Б?, а также учитывая свойства симметрии тензора напряжений и соотношение Коши, связывающее тензор напряжений с полным вектором усилий на площадке (3), соотношение (4) преобразуается к виду

j о --(Vw)s dV

V

где Бг - часть поверхности тела, на которой заданы силовые (естественные) граничные условия.

Механическая интерпретация соотношения (5) достаточно очевидна. В самом деле, поскольку весовая функция может быть выбрана произвольно при выполнении определенных ограничений на ее вид, то в данном случае удобно принять за весовую функцию возможные перемещения точек сплошной среды w=5u. Тогда уравнение (5) будет представлять собой одну из форм принципа возможных перемещений: сумма работ внешних сил на возможных перемещениях равна вариации потенциальной энергии деформации тела, которая, в свою очередь, также может быть интерпретирована как внутренняя работа обобщенных сил (напряжений) на обобщенных перемещениях (деформациях). Напряжения при этом являются функциями перемещений и давления жидкости в порах согласно второму и третьему соотношениям системы (1).

Аналогичным образом преобразуется четвертое уравнение системы (1):

[ (а — -V- кУр + а —+ = 0 (6)

V & д1

где w - скалярная весовая функция.

Второй интеграл в соотношении (6) может быть вычислен аналогично предыдущему:

| V - (kVp)wdV = | {V - (кУрм) - кУр - Ум}с1У =

V V

= | п - (кУр -1 кУр - VwdV.

Выберем весовую функцию таким образом, чтобы она равнялась нулю на части поверхности Б?, где задано давление, и представим интеграл по поверхности в виде суммы двух интегралов. Принимая во внимание известную формулу вычисления производной по нормали к поверхности и учитывая естественное граничное условие на давление (3), окончательно получим второе исходное выражение для дальнейшей разработки алгоритма:

др

V V

Ш (7)

+ |а—wdV = | qwdS -^YwdV

V д V

Для численного решения уравнений (5) и (7), используя стандартный подход конечно-элементной дискретизации [8], введем аппроксимирующие функции для основных неизвестных задачи - перемещений точек твердой фазы и давления жидкости в произвольной точке конечного элемента:

j a—wdV + j kVp - VwdV

+

u( х, y, z, t) = N u (X, y, z )ue (t) = Nu (x, y, z) auU(t), (x, y, z) £ Ve p( X, y, z, t) = Np (x, y, z) pe (t) = = Np (x, y, z) aeP(t), (x, y, z) £ V&

Ne

Ne

(8)

где - матрицы, образованные локальными

интерполирующими функциями, определенными в области конечного элемента Vе; aце, apе - матрицы кинематических связей, состоящие из нулей и единиц; ие и pе - элементные вектора узловых переменных - перемещений и давлений в узлах конечного элемента, которые, в свою очередь, выражаются через глобальные вектора узловых переменных U и P. Последние содержат перемещения и давления в узлах глобальной конечно-элементной сетки.

Аналогичные интерполяционные соотношения можно ввести для функций, заданных в объеме тела, для объемных сил и внутренних источников, что соответствует так называемому изопараметри-ческому подходу в теории конечных элементов [8].

Функции, заданные на границе области, т.е. поверхностные силы и потоки, интерполируются по узловым значениям элементов, выходящих на поверхность области V, аналогичным образом.

Весовые функции w и № могут быть выбраны различным образом. Наиболее распространенным является вариант метода Бубнова-Галеркина, в котором весовые функции задаются в том же виде, что и приближенные решения:

х, у, г, г) = Nеи (х, у, г^е (г) = = N (х, у, г)а^(г), (х, у, г) е К w( х, у, г, г) = N р (х, у, г ) (г) =

= NP (x, y, z)aPW(t), (x, y, z) £ Vе

(9)

рр

Известное преимущество такого подхода состоит в образовании симметричных матриц, входящих в глобальную систему алгебраических уравнений.

В соответствии с общей идеей МКЭ интегралы по объему и поверхности тела представляются в виде суммы интегралов по конечным элементам, на которые разбита область и граница области. В результате уравнения (5) и (7) примут следующий вид:

Ne Ne

2 JwTfsdSe + 2 JwTfvdV

e=l a'e e=l Tre

Ne

+ 2 JwTfindVe = 2 J[(Vw)S]T GdVe

e=l V e

e=l V e

2 J awT —dVe + 2 J k (Vw)T (Vp)dV

e=1 Ve dt

+

e=1 Ve

Ne

+ 2 Ja wT d-dVe =

e=1 vU dt

Ne Ne

2 J wTqdSe -2 J wTydVe

e=1 r,''e S2

e=1 Vе

Уравнения (10) путем применения стандартных математических процедур метода конечных элементов приводятся к виду

мй + Кй - HP = + Р!В(1>,

йР + ОР + Сй = 0(1) + Г(0, (11)

где М, К - глобальные матрицы масс и жесткости; Н - глобальная матрица связи «давление - перемещения», или матрица влияния распределения давления жидкости в порах на движение твердой фазы; Fv(f), Fs(f) - глобальные узловые вектора, статически эквивалентные заданным объемным и поверхностным силам; й и О - глобальные матрицы насыщения и проницаемости; С - глобальная матрица связи «перемещения - давление», или матрица влияния перемещений твердого каркаса на распределение давления жидкости в порах; 0 и Г -глобальные узловые вектора, статически эквивалентные заданным объемным источникам и поверхностным потокам.

Глобальные матрицы масс и жесткости, а также глобальные вектора узловых сил имеют стандартный вид соответствующих матриц и векторов, получаемых в результате конечно-элементной дискретизации уравнений упругости [8], и поэтому здесь не приводятся. Остальные матрицы могут быть записаны следующим образом:

H = 2 (aU )T J (Be )T AeNepdVeap ,

e=1 Ve

Ne

D = 2 (ap )T J a (N p )T NepdVeaep , (12)

e=1 Ve

Ne

G = 2(ap)T J(Gp)TkGepdVeaep ,

e=1 Ve

где интегралы по элементам представляют собой соответствующие элементные матрицы; Ae - элементная матрица коэффициентов эффективных напряжений; Be - матрица градиентов; Gpe - матрица частных производных от интерполирующих функций, сходная с матрицей градиентов Be.

Можно показать, что матрицы H и C совпадают с точностью до операции транспонирования, что подчеркивает симметричность эффекта влияния движения фаз друг на друга в связанной задаче пороупругости.

При анализе динамического поведения биологических структур, необходимо учитывать вязкость живых тканей. В теории наследственной упругости и в работах, посвященных моделированию

e

мягких тканей, известно большое число определяющих соотношений вязкоупругого материала, задающих связь между скоростями деформаций и напряжений. Полученная выше система динамических уравнений (11) позволяет интерпретировать явление диссипации энергии в биологических тканях как результат вязкого взаимодействия жидкой и твердой фаз в пороупругом теле.

В самом деле, рассмотрим частный случай общей динамической задачи -вынужденные колебания тела под действием внешней гармонической силы.

Представим внешние силы, потоки и источники в гармоническом виде с помощью комплексной экспоненты :

Е(/) = , Е(0 = , Е,Е еЯе

О/) = , Г(0 = Гвш, О, ГеЯе где Ру и Рэ - глобальные вектора амплитудных значений объемных и поверхностных узловых сил; О и Г - глобальные вектора амплитудных значений потоков и источников; т - круговая частота колебаний; / - мнимая единица.

В этом случае частное решение, соответствующее установившимся вынужденным колебания, определяется в том же виде, что и правая часть дифференциальных уравнений:

и(0 = иёт/, Р(0 = Р ёт/, и, Р е С где и и Р - глобальные вектора комплексных амплитудных значений перемещений и давлений в узлах конечно-элементной сетки.

Подставляя выражения (13) и (14) в уравнения (11), получим линейную систему комплексных алгебраических уравнений относительно амплитудных значений перемещений и давлений в узлах:

(К - т2М)и - НР = К + Е

V * (15)

гт Си + (О + 7юБ)Р = О + Г

Поскольку матрицы Н и С совпадают с точностью до операции транспонирования, то глобальная матрица системы (15) может быть приведена к симметричному виду, а уравнения примут следующий вид:

(К -т2М)и - НР = е + Е

-Н Т и + (-Б + гО / т)Р = г (О + Г)/т

(16)

После решения полученной системы для фиксированных значений частоты внешней силы могут быть также вычислены амплитудные значения узловых скоростей и ускорений. При этом действительные амплитуды и фазы колебаний определены согласно известным формулам представления комплексных чисел в экспоненциальной форме.

Рассмотренный алгоритм реализован в авторском конечно-элементном комплексе МесИапю-бРЕ_РЕУ и позволяет рассчитывать установившиеся колебания пористых упругих тел, насыщенных вязкой жидкостью. Отметим, что в настоящее время авторам неизвестны описания доступных коммерческих программ, позволяющих проводить подобные расчеты.

Результаты

Для тестирования программы и верификации получаемых результатов необходимо иметь про-

стые аналитические решения. Как и в классической теории упругости, общая система уравнений поро-упругости (1) допускает частные случаи, описывающие продольные и поперечные колебания стержня:

-Еи'' + ри +ар' = аи'- кр'' + ар = у

Е^'7 + рБ\& + аР'' = ду = 11ууСв в

- а^'' - кР' + аР = у у = | уу Св

(17)

(18)

где и=и(х,1) - продольные перемещения; ч= -прогиб стержня; Р =| урс№ - эквивалентное давление в сечении; Б и Л - площадь и момент инерции поперечного сечения стержня; у - вертикальная координата.

В книге [10] приведены другой вывод этих уравнений и некоторые аналитические решения, полученные путем разложения неизвестных функций в ряд Фурье при определенных граничных условиях.

В качестве числового примера был рассмотрен стержень длиной 30 см, имеющий сечение в виде квадрата со стороной 1 см, нагруженный продольной или поперечной объемной силой интенсивности 1 Н/мм3 и 0.01 Н/мм3 соответственно, изменяющейся во времени по гармоническому закону. В первом случае стержень консольно закреплен на левом торце при свободном правом, а во втором граничные условия представляют собой шарнирное опирание при свободных для относительного движения жидкости торцах. Модули упругости и характеристики пористости соответствуют компактной ткани бедренной или берцовой кости [7]: модуль Юнга Е = 14.575 ГПа; модуль сдвига Э = 5.5 ГПа; плотность р = 1500 кг/м3; проницаемость к = 6.36-10-10 м2/Па-с; пористость ф = 0.04; объемный модуль упругости жидкости К = 2.3 ГПа. В качестве параметра, меняющегося в пределах от 0 до 1 [4, 7], примем коэффициент эффективных напряжений Бийо а.

Заметим, что применение формулы (2) в случае, когда коэффициент эффективных напряжений приближается к нулю, становится проблематичным. В самом деле, выражение коэффициента Пуассона в условиях насыщения Vu имеет следующий вид [4]:

3у + а£(1 - 2у)

,, - V / (19)

3 -аВ(1 - 2у)

а

Тогда при а = 0 коэффициент vu совпадает с коэффициентом Пуассона в условиях отсутствия жидкости в порах V, и, следовательно, множитель а, согласно формуле (2), становится неопределенным. Преодолеть возникшую проблему можно, воспользовавшись определением коэффициента Скемптона [4]:

аК г

В = --—-^-;--(20)

(а-ф(1 -а)) +фК

Г

где К - эффективный объемный модуль упругости пористого материала в условиях отсутствия жидкости в порах.

Подставляя выражение (20) в формулу (2) и выполняя ряд алгебраических преобразований, окончательно получим

3(1 -2v) -а2 +а(1 + ф)-ф

1 -

К

\\

К

а = ■

Г ))

(21)

20(1 + V)

В результате решения уравнений (17) и (18) получены аналитически точные амплитудно-частотные характеристики (АЧХ) продольных и поперечных колебаний одномерной модели стержня при рассмотренных выше граничных условиях. На рис.1, 2 показаны соответствующие резонансные кривые, отображающие зависимость амплитуды продольных перемещений и прогиба от частоты внешней гармонической силы в точках стержня, расположенных на правом свободном торце в первом случае или в среднем сечении во втором случае. На рис.3, 4 показаны резонансные кривые, соответствующие зависимости давления жидкости в порах от частоты продольной или поперечной возмущающей силы в точках стержня, расположенных на левом закрепленном торце стержня в первом случае или в среднем сечении стержня во втором случае. В качестве примера фазо-частотной характеристики (ФЧХ) приведена зависимость фазы поперечных колебаний точек среднего сечения стержня (рис.5) от частоты приложенной объемной нагрузки.

Аналитически рассчитанные кривые сравнивались с приближенными результатами, полученными на трехмерной конечно-элементной модели стержня, имеющего те же размеры и значения материальных констант и нагруженного объемной силой той же интенсивности. Модель стержня состоит из 10 конечных элементов, расположенных вдоль оси стержня. Каждый элемент представляет собой 20-узловой шестисторонний изопараметрический элемент, имеющий 4 степени свободы в узле: 3 компоненты вектора перемещений твердой фазы и значение давления жидкости. Общее число степеней свободы модели составляет 512, из которых 384 - кинематические переменные. Полученные в результате расчета амплитудно-частотные и фазо-частотные характеристики стержня в тех же сечениях, что и в случае аналитической модели стержня, представлены на рис. 1-5.

Обсуждение

Представленные графики (рис.1 - 5) амплитудно-частотных и фазо- частотных характеристик стержня, полученные аналитическим путем на основе одномерной модели и методом конечных элементов с использованием пространственной модели стержня в виде длинного параллелепипеда, демонстрируют хорошее соответствие точных и приближенных результатов. Имеющееся несоответствие вполне укладывается в рамки введенных допущений, накладываемых на соотношения между характерными размерами стержня, и допущений, связанных с характером распределения давления жидкости в порах по сечению стержня. Наибольшая ошибка при расчете резонансной частоты продольных и поперечных колебаний методом конечных элементов относительно аналитических данных составляет gf прод = 0.33 % при а = 0.4 в первом случае и gf попр = 0.8 % при а = 1 во втором случае.

Сравнительный анализ амплитудно-частотных характеристик стержня, совершающего

продольные или поперечные колебания (рис.1 - 5), позволяет сделать ряд выводов об упругих и дисси-пативных свойствах конструкции, построенной из двухфазного материала. Можно отметить два характерных свойства пороупругих систем, независимо от типа прикладываемой нагрузки.

Во-первых, резонансная частота стержня в обоих случаях смещается вверх при увеличении значения коэффициента эффективных напряжений Бийо а в заданных пределах от нуля до единицы. Как следует из уравнений (11), (17) и (18), при а = 0 система уравнений перестает быть связанной, распадаясь на два независимых уравнения, описывающих самостоятельное движение твердой фазы и распределение жидкости в порах, которое в данном случае равно нулю, если отсутствуют внешние потоки на поверхности $2 и внешнее давление равно

нулю на поверхности $1. Таким образом, коэффициент эффективных напряжений является основным параметром, отвечающим за эффект коуплинга в задаче пороупругости. Во-вторых, имеет место диссипация энергии в системе, выражаемая в том, что амплитуда колебаний становится ограниченной при совпадении частоты возбуждающей силы с собственной частотой колебаний стержня. Более того, при увеличении коэффициента эффективных напряжений от нуля до единицы эффект демпфирования колебаний усиливается. Интересно отметить, что возможно практически полное исчезновение резонансного пика при больших значениях этого параметра, что наблюдается при а > 0.8 в случае поперечных колебаний стержня (рис. 2). Еще более необычные резонансные свойства проявляет система, совершающая изгибные колебания, при очень больших значениях коэффициента а, а именно, при его максимальном значении, равном единице, амплитудно-частотная характеристика имеет сильно сглаженный резонансный пик на частоте 125 Гц, которая значительно меньше собственной частоты колебаний идеально упругого стержня, а также резонансных частот пороупругого стержня при значениях коэффициента эффективных напряжений а = 0.4 и а = 0.8.

Из общих теорем аналитической механики известно, что о наличии диссипации в системе также свидетельствуют фазо-частотные кривые. В качестве примера в статье приведены фазо-частотные характеристики поперечных колебаний стержня (рис. 5). При нулевом значении коэффициента эффективных напряжений Бийо сдвиг фазы колебаний на 180° происходит мгновенно при переходе системы через резонанс. В случае же его увеличения до единицы наблюдается значительное сглаживание ступенчатой функции.

Поведение давления в порах в целом повторяет характер зависимости продольных и поперечных перемещений твердого каркаса стержня от частоты внешнего воздействия (рис.3, 4). Логично предположить, что основной вклад в движение жидкой фазы вносят упругие перемещения твердого каркаса, который физически вовлекает в свое движение жидкость, содержащуюся в порах. Имеющее место относительное движение жидкости, видимо, незначительно, что подтверждает некоторое смещение резонансного пика в сторону более высоких частот по сравнению с АЧХ колебаний точек твердого каркаса пористого материала. Данный эффект стано-

вится более заметным при увеличении коэффициента эффективных напряжений.

Интересно отметить существенное различие в чувствительности двух систем к значению коэффициента Бийо, что проявляется в значительном смещении резонансной частоты при рассмотрении продольных колебаний стержня. Максимальное повышение резонансной частоты составляет в этом случае 122.3 % при а = 1, в то время как аналогичный параметр для поперечных колебаний составляет лишь 6.05 % при а = 0.8. Заметим, однако, что влияние коэффициента эффективных напряжений на изгибные колебания оказывается более значительным в смысле повышения внутреннего трения в

системе и, как следствие, увеличения диссипации выше критического значения, что приводит к значительному снижению амплитуды колебаний. Причиной данного эффекта может быть как различное соотношение между направлениями прикладываемой нагрузки (вдоль оси стержня при продольных колебаниях и перпендикулярно - при поперечных) и преимущественным направлением движения жидкости в порах (вдоль продольной оси стержня в обоих случаях), так и значительное различие в собственных частотах идеально-упругого стержня на растяжение и изгиб.

АЧХ - Растяжение консольного стержня (nu=0.325)

alpha-0.0 (ANLT ) --alpha=0.4 (ANLT) -------alpha=0.E (ANLT ) ----alpha=1.0 (ANLT)

ж alpha=0.0 (FE3D ) о alpha=0.4 (FE3D ) + alpha=0.! (FE3D ) - alpha=1.0 (FE3D )

Рис. 1. Амплитудно-частотные характеристики продольных колебаний консольного стержня, рассчитанные в средней точке правого свободного торца стержня аналитически (Д1\ИТ) и численно (РЕЗР).

АЧХ - Изгиб опертого стержня (пи=0.325)

частота, Гц

а1рИа-0.0 (АЫ1_Т) --а1рИа=0.4 (АИИ) -------а1р И а=0 Ё (АМТ) ---- а1р11а= =1.0 (АИИ)

Ж а1рИа=0.0 (РЕЗО ) о а1рИа=0.4 (РЕЗО ) а1рИа=0.£ (РЕЗО ) а1р11а= =1.0 (РЕЗО )

Рис. 2. Амплитудно-частотные характеристики поперечных колебаний шарнирно опертого стержня, рассчитанные в среднем

сечении стержня аналитически (ДЫЬТ) и численно (РЕ3Р).

АЧХ - Растяжение консольного стержня (пи=0.325)

частота, Гц

-а1рИа=0.0 (АМН)--а1рИа=0.4 (А[\11_Т) -------а1рИа=0.8 (АМТ) ----а1рИа=1.0 (АЮ)

ж а!рИа=0.0 (РЕЗО) о а1рИа=0.4 (РЕЗО) + а1рИа=0.8 (РЕЗО) - а1рИа=1.0 (РЕЗО )

Рис. 3. Резонансные кривые, соответствующие зависимости давления жидкости в порах от частоты продольной возмущающей силы, рассчитанные в средней точке левого закрепленного торца стержня аналитически (ДЫЬТ) и численно (РЕЗР).

АЧХ - Изгиб консольного стержня (пи=0.325)

200 т —— —— —— —— —— ——

120 140 160 180 200 220 240 260 280

частота, Гц

а1рИа-0.0 (АМ1_Т) а1рИа-0.4 (АГЛТ) ------ - а1р11а=0.£ (АМ1_Т) ----а1рИа=1.0 (АШ")

Ж а!рИа=0.0 (РЕЗР ) о а!рИа=0.4 (РЕЗО ) а1р(1а=0.£ (РЕЗР ) а1рИа=1.0 (РЕЗР )

Рис. 4. Резонансные кривые, соответствующие зависимости давления жидкости в порах от частоты поперечной возмущающей силы, рассчитанные в среднем сечении стержня аналитически (АГ\ИТ) и численно (РЕЗО).

ФЧХ - Изгиб консольного стержня (пи=0.325)

а1рИа-0.0 (АЫ1_Т) --а1рИа=0.4 (АЫЬТ) -------а1р|1а-0.Е (АМТ) ---- а1рИа= =1.0 (А1\11_Т)

Ж а!рИа=0.0 (РЕЗО ) о а1р11а=0.4 (РЕЗР ) а1рИа=0.£ (РЕЗО ) а1рИа= =1.0 (РЕЗО )

Рис. 5. Фазо-частотные характеристики поперечных колебаний шарнирно опертого стержня, рассчитанные в среднем сечении

стержня аналитически (АЫЬТ) и численно (РЕ3й).

Заключение

Полученные результаты численного моделирования имеют хорошее совпадение с аналитическими данными, что подтверждает работоспособность предложенного конечно-элементного алгоритма и позволяет использовать программный комплекс Мес11атс8РЕ_РЕУ для расчета сложных биомеханических систем, насыщенных жидкостью, таких, как длинные трубчатые кости опорно-двигательного аппарата, связки, хрящи, сухожилия и скелетные мышцы, для математического описания которых аппарат теории двухфазных сред представляется наиболее естественным. Разработанный подход может быть применен для вибрационного анализа слоистых пористых структур и пенообразных тел, широко используемых в системах звукоизоляции, а также в геомеханике при анализе распространения упругих волн в насыщенных породах и грунтах.

Список литературы

1. Регирер С.А., Штейн А.А. Методы механики сплошной среды в применении к задачам роста и развития биологических тканей // Современные проблемы биомеханики. 1985. № 2. С.5-37.

2. Нигматулин Р.И. Динамика многофазных сред. М.: Наука, 1987. 464 с.

3. Biot M.A. General theory of three-dimensional consolidation // J. Appl. Physics. 1941. V. 12. No 2. P.155-164.

4. Cowin S.C. Bone poroelasticity // J. Biomech. 1999. V. 32. No 3. P.217-238.

5. Roberts S.G., Steele C.R. Efficacy of monitoring long-bone fracture healing by measurement of either bone stiffness or resonant frequency: numerical simulation // J. Orthop. Res. 2000. V. 18. №. 5. P.691-697.

6. Georgiou A.P., Cunningham J.L. Accurate diagnosis of hip prosthesis loosening using a vibrational technique // Clin. Biomech. 2001. V. 16. № 4. P.315-323.

7. Smit T.H., Huyghe J.M., Cowin S.C. Estimation of the poroelastic parameters of cortical bone // J. Biomech. 2002. V. 35. P.829-835.

8. Zienkiewicz O.C., Taylor Robert L., Taylor R.L. Finite Element Method: Volume 1, The Basis. Butterworth-Heinemann, 2000. 712 p.

9. Лурье А.И. Теория упругости. М.: Наука, 1970. 932 c.

10. Cederbaum G., Li L.P., Schulgasser K. Poroelastic structures. Elsevier Science Ltd., 2000. 158 p.