CC BY
26
Russian Journal of Biomechanics
www.biomech.ac.ru
НЕИЗОТЕРМИЧЕСКОЕ СТОХАСТИЧЕСКОЕ СМАЗЫВАНИЕ ТАЗОБЕДРЕННОГО СУСТАВА ЧЕЛОВЕКА С РАЗЛИЧНЫМИ ЧАСТОТАМИ И АМПЛИТУДАМИ
К.Х. Вежхольский
Base Technique Department, Maritime University of Gdynia, Morskastr. 83, PL-81-225, Gdynia, Poland, e-mail: wierzch@am.gdynia.pl
Кафедра основ техники, Морской Университет, Гдыня, Польша
Аннотация. Данная статья содержит аналитические и численные расчеты давления и нагрузки при гидродинамической смазке зазора в тазобедренном суставе человека в условиях неизотермического нестационарного периодического движения. Рассмотрена динамика колебаний с различными частотами и амплитудами головки бедра, вертлужной впадины и синовиальной жидкости. Принята во внимание шероховатость суставных поверхностей. Измерения поверхностей хрящей выполнены с помощью лазерного датчика. Предполагается, что сферическая головка бедра в тазобедренном суставе человека совершает движения в двух направлениях - окружном и меридиональном. Течение синовиальной жидкости описывается уравнениями сохранения импульса и уравнением непрерывности. Численные расчеты выполнены в пакете Mathcad 11 Professional Program с использованием метода конечных разностей. Этот метод обеспечивает устойчивость численного решения уравнений в частных производных и дает правильные значения давления и нагрузки в тазобедренном суставе человека.
Ключевые слова: тазобедренный сустав человека, смазывание
тазобедренного сустава, периодическое движение, различные частоты и амплитуды.
1. Введение
Тазобедренный сустав в синовиальной оболочке представляет собой настолько сложную систему, что даже современные технологии не могут полностью повторить ее. В этой системе используется относительное скольжение гладких сферических поверхностей для обеспечения вращения конечности при приложении к ней различно направленных и зависящих от времени нагрузок.
В данной работе предполагается, что сферическая головка бедра совершает вращательные движения вокруг двух осей с различными угловыми скоростями, см. рис. 1. Символом ф обозначена окружная координата, символом r - координата в
© Вежхольский К.Х., 2005
направлении высоты зазора, символ & обозначает меридиан костной поверхности вращения, или долготную координату. В работах [1, 5, 9] обсуждаются
гидродинамические и механические параметры тазобедренного сустава человека. Как правило, там используются как аналитические, так и численные методы.
В работах [3, 6, 13, 14, 18, 24] исследуются силы трения в различных суставах человека для различной геометрии соприкасающихся костных поверхностей при изменяющейся высоте суставного зазора и несимметричном течении синовиальной жидкости. В вышеупомянутых работах [1, 5, 7, 8, 9, 12, 19, 20, 21, 23] не рассматривалось периодическое движение с различными частотами головки бедра и вертлужной впадины тазобедренного сустава человека. Также не принималось во внимание вязкоупругое смазывание при случайных условиях. Новизна данной работы состоит в вычислении нагрузки в тазобедренном суставе человека для шероховатых поверхностей кости и хряща, которые смазываются благодаря периодическому течению с различными частотами и амплитудами.
При несимметричном течении синовиальной жидкости три компоненты уф, уг, У& скорости жидкости по направлениям ф, г, & зависят от переменных ф, г, & и времени Давление р и вязкость синовиальной жидкости зависят от ф, &, Высота в зазора есть функция переменных ф, & и ¿. Полученное в данной работе решение задачи о смазке дает давление и нагрузку для неньютоновой синовиальной жидкости во втором порядке точности. Предполагается, что колебания с различными частотами и амплитудами в окружном и меридиональном направлениях происходят в синовиальной жидкости и на поверхностях головки бедра и вертлужной впадины. Учитываются случайные изменения высоты зазора между костями.
2. Основные уравнения
В тазобедренном суставе человека синовиальная жидкость течет между головкой бедра и вертлужной впадиной. Задача о смазке в суставном зазоре будет
решена с помощью уравнения сохранения импульса, уравнения непрерывности и уравнения энергии [20, 21]:
Кроме того, общее определяющее соотношение Ривлина и Эриксена во втором порядке точности записывается в следующем виде [2, 15, 16, 17]:
где 8 - тензор напряжений, р - давление, I - единичный тензор, Аі и А2 - первые два тензора Ривлина - Эриксена и а, Р - три материальные константы - динамическая вязкость). Через Т обозначена температура, к - теплопроводность синовиальной жидкости, еу - теплоемкость синовиальной жидкости. Через А1 и А2 обозначены симметричные матрицы, определенные в [15, 16, 17]:
движется относительно вертлужной впадины в окружном направлении ф и меридиональном направлении д. Поверхности вертлужной впадины и головки бедра колеблются в направлениях ф и д с различными амплитудами и частотами. Кроме того, вертлужная впадина колеблется в направлении высоты зазора. Такое движение суставных поверхностей вызывает течение синовиальной жидкости в суставном зазоре. Угловую скорость вращения головки бедра в окружном направлении обозначим через ш1, в меридиональном направлении - через ш3. В суставном зазоре имеется несимметричное нестационарное течение синовиальной жидкости, при этом жидкость проявляет вязкоупругие и нестационарные свойства. Центробежными силами пренебрегаем. Введем следующие обозначения: и = &1Я - окружная скорость, р = р0 -постоянная плотность синовиальной жидкости, ^ = ^0Л1 - переменная вязкость, ^0 -характерное значение динамической вязкости, в - зависящая от времени высота суставного зазора, Я - радиус головки бедра, 10 - характерное значение времени; течение жидкости предполагается изотермическим.
Предполагаем, что произведение чисел Дебора (ОвЬвгаИ) Бе = Рш1/^0 и Струхала 81х=ЯШ0 , то есть Бе 81х, и произведение числа Рейнольдса Яе=рив/^0, безразмерного зазора у и числа Струхала, то есть Ке81ху, являются величинами одного порядка и что Бе 81х >> Бе [14]. Учтем эти оценки в основном уравнении течения тонкого слоя в сферическом суставном зазоре [14]. Пренебрегая величиной у^в/Я « 10_ и центробежными силами, получим в вышеупомянутых предположениях следующие уравнения в сферических координатах ф, г, д [12, 19, 20, 21, 22]:
Div S = р dv/dt, div v =0,
р d(cvT)/dt = div(^ grad T) = div(vS) - v Div S.
(1)
S = - p I+^A1+a(A1)2+PA2,
(2)
Ai = L+LT, A2 = grad a +(grad a)T+2LTL, a = L v +, L = grad v, LT = (grad v)T, (3)
где v - скорость (м/с), t - время (с), a - ускорение (м/с2). Сферическая головка бедра
dt
(4)
dr ’
(5)
дУф ^т? • (?1дк д - + К БІЙ \ — I—- Н----
I К ) дг д$
дф
= 0.
д г
к-
дТ_ д г
Гду„ у (^
V дг )
+
Vд г )
= 0,
(7)
(8)
где 0 <ф < 2л 01, 0 < 01 <1, Ьт = лЯ/8 < Э < лЯ/2 = Ъц, 0 < г < в, в - высота зазора. Члены с множителем Р описывают влияние вязкоупругих свойств синовиальной жидкости на эффекты смазки. Конвективными членами пренебрегаем. В левых частях (4) и (6) оставлены только производные по времени от компонент скорости. В уравнении энергии учтены диссипативные члены.
д
3. Метод решения задачи о периодическом смазывании
Так как уравнения (4) - (7) линейны, то можно разделить стационарное и нестационарное течения жидкости. Учтем касательные и вертикальные ускорения точек суставной поверхности. Запишем компоненты скорости жидкости и давление в виде сходящихся рядов [3, 11]:
ад
V = уг(0)(ф,г, Э)+ 2 ^}(ф,г, Э)ехр(1кШо*), I = Ф, &, (9)
к=1
ад
V = у|:0)(ф, г, Э) + 2 у(гк )(ф, г, Э)ехр(1к ш0* ), (10)
к=1
ад
Т = Т (0)(ф, г, Э) + £ Т(к) (ф, г, Э) ехр(1к ш0* ), (11)
к=1
ад
р = р (0)(ф, Э) + 2 Р(к )(Ф, Э)ехР('кш0*). (12)
к=1
Через ш0 обозначена частота колебаний (с-1); это частота периодических возмущений
нестационарного течения жидкости. Через I = л/-Г обозначена мнимая единица. Неизвестные функции с верхним индексом (0) есть компоненты скорости и давление для стационарного течения не вязкоупругой жидкости. Неизвестные функции с
верхним индексом (к) для к = 1, 2, 3,..., являются поправками, обусловленными
нестационарными вязкоупругими свойствами жидкости. Высота зазора имеет следующий вид [4, 10]:
ад
вш = в(0) +в + 5 = в(0)(ф, Э) + ^в(к )(ф, Э)ехр(1кшг*) + 5, (13)
к=1
где 5 (£) - случайная добавка, возникающая из-за колебаний и неровностей поверхности, отсчитываемых от среднего уровня, £, - случайная величина, в -нестационарное возмущение высоты зазора, вызываемое нестационарными условиями работы сустава, в(к) - не зависящие от времени коэффициенты в возмущениях высоты зазора, в(0) - начальное значение высоты зазора, в0 - общая высота зазора,шг - частота колебаний в направлении высоты зазора. Детали метода решения см. в Приложении 1.
головка бедра
3=tcR/8
4. Граничные условия
Расположение области смазывания Q: 0 < ф < л, лЯ/8 < 9<лЯ/2 показано на
рис. 2. Вертлужная впадина движется в окружном ф и меридиональном $
направлениях. Кроме того, учитываем зависящее от времени касательное ускорение точек поверхности головки бедра. Касательные (в направлениях ф и 9) компоненты скорости (U) на поверхности головки бедра и компоненты скорости (V) на поверхности вертлужной впадины, а также температура имеют следующий вид:
ад ад
U i = Uo + ZUik exp(ik®¿), Vi = Zvj exp0^),
k=i k=i (14)
ад v '
T = To при r = 0, T1 = f (o) +Z Tk exP (ik&iJ) при r =e,
k=1
U V f(k)
Uik = -рт, Vik = J0 = c0ns1, Tk ^A=®10R sin (9 / R), UаД-®зоR sin ф ,
где Uik - не зависящие от времени амплитуды касательных колебаний на поверхности головки бедра, i = ф, 9, к = 1, 2, 3,...., Vik, Р - не зависящие от времени амплитуды
колебаний касательной скорости и температуры на поверхности вертлужной впадины,
i = ф, 9; к = 1, 2, 3,.... Черезf0 обозначена не зависящая от времени часть температуры вертлужной впадины. Через q9V, шфМ, o9v, о9и обозначены различные частоты колебаний вертлужной впадины (нижний индекс V) и головки бедра (нижний индекс u) в направлениях ф и 9. Эти частоты отличаются от частоты ш0, соответствующей синовиальной жидкости. Окружная скорость на поверхности сферической головки бедра не зависит от ф и достигает максимального значения QiR на экваторе при 9 = Rtc/2, а нулевого значения - на полюсе при 9 = 0 или 9 = Rtc. Меридиональная скорость на поверхности сферической головки бедра не зависит от 9 и достигает максимального значения q3R при ф = л/2, а нулевого значения - при ф = 0 и ф = л. Следовательно, не зависящие от времени части окружной и меридиональной скоростей для головки бедра при стационарном движении имеют вид:
Ut0 =q1Rsin(9/R), U90 =q1Rsin^), 0 <9<Rл,0 <ф<л. (15)
Обозначим через ш1 и ш1о угловую скорость и поправку к ней для вращения головки бедра в ф - направлении; через ш3 и ш30 - угловую скорость и поправку к ней для вращения головки бедра в 9-направлении. Внутри суставного зазора на скорости (9, 10) накладывается граничное условие (14) 1. Приравниваем первые стационарные члены правой части (9), представляющие собой компоненты скорости жидкости, стационарным членам правой части (14)1. Компоненты скорости жидкости на поверхности неподвижной вертлужной впадины равны нулю. На поверхности головки бедра равна нулю только вертикальная компонента скорости жидкости, так как учитываются вязкие свойства жидкости. Компоненты скорости жидкости на поверхности головки бедра в направлениях ф и 9 равны соответствующим компонентам скорости головки бедра. Граничные условия на стационарную часть компонент скорости жидкости для i = ф, 9; к = 1, 2, 3,... на поверхности головки бедра при r = 0 и на поверхности вертлужной впадины при r = s имеют следующий вид:
уф0)(ф,r = о,9) = U*, (9), у£0)(ф,r = 0,9) = U90 (ф), (16)
уг(0)(ф, r = 0,9) = 0, T (0)(ф, r = 0,9) = T0, (17)
.(0)(ф,r = s,9) = 0, при i = ф, r,9; T(0)(ф,r = s,9) = /(0). (18)
Теперь наложим нестационарную часть граничных условий (14)1 на
нестационарную часть компонент скорости (9, 10, 11). Однако сопоставить нестационарные члены в рядах (14)1 и (9, 10, 11) невозможно, так как частоты в обеих частях уравнений различны. Поэтому ряды в правых частях (14)1, (14)2, описывающие колебания с известными амплитудами ифк, U9k, ¥фк, V9k, Тк и известными частотами шфи Q9u, Qv Q9v, преобразуются (см. Приложение 2) в ряды с некоторыми новыми, неизвестными, амплитудами ифк ,U^k ,V*k , V*, T* и известной частотой ш0, отвечающей синовиальной жидкости. В Приложении 2 также приведено вычисление амплитуд и* ,vk [2, 9]. В результате вычислений новые неизвестные амплитуды
Uфк,U9*k,V*k,V*k,Wk*,Т* получаются как функции известных амплитуд U9k, U9k, Vфk, V9k,
Wk, Tk, известных, экспериментально измеренных, частот Шфм, 09u, 0ф., Q9v, Qr и известной частоты ш0, отвечающей колебаниям синовиальной жидкости [2]. Например, через шфМ, ю9и обозначены частоты колебаний головки бедра в ф и 9 направлениях. Через Q(pv, Q9v обозначены частоты колебаний вертлужной впадины в ф и 9 направлениях. Теперь можно приравнять нестационарные члены в правых частях (9, 10) нестационарным членам в правых частях (14)1, (14)2. Граничные условия на нестационарные поправки к температуре и к компонентам скорости i = ф, 9; к= 1, 2, 3,... имеют на поверхности головки бедра при r = 0 и на поверхности вертлужной впадины при r = s следующий вид [3, 11]:
v;'(ф, r = 0,9) = и;(и,к,), v9k>(ф, r = 0,9) = U9k(U9k,09,), (19)
vfV r = 0,9) = 0, r = 0,S) = 0, (20)
»¡к'(ф, r = s, 9) = V^Vk, и,.), при i = ф, 9, T(k) (ф, r = s, 9) = T*. (21)
Высота суставного зазора со временем меняется. Вертикальная компонента скорости жидкости на поверхности вертлужной впадины равна производной по времени от высоты зазора. Следовательно, используя (13), получаем:
да Яр
Z vf) (Ф, r = р, 9) exp (¡;ш ot) = -рО- =
к= (22)
да да да
= Zp( k )ikrar exp (\k&rt) = i Z Wk exp (ik ort) = iZ W* exp (ikra0t).
к=1 k=1 k=1
Приравнивая амплитуды в обеих частях выписанных выше рядов, находим:
vf)(Ф,r =р,9) = ¡w;(Wk,и,), Wk = e(k)ки, (23)
для к = 1, 2, 3,... , р(к) = В^^/к5, где 0 < Be <1. Безразмерный коэффициент Be отвечает за не зависящие от времени возмущения высоты зазора. Не зависящая от времени средняя высота зазора с возмущениями имеет вид [3]:
у(к)(ф, r = р, 9) = iw; (Wk, и,), Wk = р( к )к и, (23)
где
р(0) (ф, 9 / R) = А рх cos ф sin 9 / R + А рy sin ф sin 9 / R -Арz cos 9 / R - R +
+[(A рx cos ф sin 9 / R + А рy sin ф sin 9 / R -Арz cos 9 / R)2 + (25)
+(R + рmin )(R + 2D + рmin )]
и ts - средний временной период возмущений суставного зазора, ^e - вещественная часть комплексного числа. Центр сферической головки бедра расположен в точке 0(0, 0, 0), а центр сферического хряща - в точке 0i( х-Арх, у-Ару, z+АрД Эксцентриситет обозначен через D (см. рис. 3). Вывод соотношения (25) приведен в Приложении 2 (см. уравнения A2.20 - A2.22).
5. Стохастические уравнения для распределения давления и сил трения
Нахождение компонент v(0\ vf) скоростей жидкости для граничных условий (16), (18)1, (19) и для постоянного значения динамической вязкости в направлении высоты зазора приведено в Приложениях 1 и 3. Граничные условия (16), (18) 1, (19), (21)1 добавляются к системе уравнений (A1.1) - (A1.3), описывающих стационарное течение, а граничные условия (19), (21)1 - к системе уравнений (A1.6) - (A1.8) для поправок (см. Приложение 1). Тогда получаем искомые компоненты скорости стационарного течения v^0), v9(0) и поправки v^k), v9(k) для к = 1, 2, 3,..., вызванные вязкоупругими свойствами и нестационарностью течения. Результаты вычислений приведены в Приложении 3.
Подставим функции v^0), v9(0) (см. (A3.1), (A3.2) в Приложении 3) в уравнение непрерывности для стационарного течения (A1.4) и подставим найденные функции v^k), v9(k) (k = 1, 2, 3,..., см. (A3.3), (A3.4)) в уравнение непрерывности для поправок (A1.9). Тогда после интегрирования по r получим компоненту скорости v/0) в направлении высоты зазора и получим поправки v/k) (к = 1, 2, 3,...) к этой компоненте (см. Приложение 4). Учитывая граничные условия (17)1, (18)1, (20)1, (23) для v/0) (см. (A4.1)) и поправок vr(k) (к = 1, 2, 3,.) к ней (см. (A4.2)), получим модифицированное уравнение Рейнольдса для давления р(0) и поправок к нему p(k) (к = 1, 2, 3,...). Детали вычислений приведены в Приложении 5.
Рис. 3. Изменение высоты зазора со временем, центр сферической головки бедра и вертлужной впадины, значения скорости на границе
Модифицированное уравнение Рейнольдса для стационарного течения синовиальной жидкости без вязкоупругих свойств определяет искомое давление р^ и имеет вид:
1 Г 9^ д I в3 др(0)
--f'OQf'f' I - ----<-------------
R [ R ,ідф{лоЛі дФ
+ R — <
dp
д9 І Л0Л1 д9
(0) . 9 1 sin — > =
R
= 6ш1 R sin І — + 6q3R2 sin (ф) — ^ R) дф д9
• Г —
в Sin l —
IR)
(26)
в Q: 0 < ф < л, tcR/8 <9 < tcR/2.
Далее рассмотрим модифицированные уравнения Рейнольдса для поправок p(k) (k = 1, 2, 3,...) к давлению, вызванных периодическим течением жидкости с вязкоупругими свойствами. Эти уравнения выведены в Приложении 5 (см. (A5.6)). Умножая (A5.6) на exp (ik&0t) и приравнивая вещественные части, получим следующую последовательность модифицированных уравнений Рейнольдса:
3
в
В3 дР
ек
Лк, дФ
+Я— дЭ
в3 дРек • (ЭЛ --------—81П! _
Лк, дЭ IЯ У
Э
= -12Щ*ЯбіпI — Iбіп(кш0/) +
+ 6 [и; («з)+V
дв ,, ч кш0р0 д і 3
— соб (кш0/) + ——(в3 дф V 0 ' 12 дф
(в3ХЛк )
-12
Я
т_* дв _( . Э\дв
V*—+я! біп— IV* —
фк
дф
Я
дЭ
• СОБ
(^)+ 6[и*к («1 )+^к ]Я
_д_
дЭ
• (Э
В біп I —
IЯ ,
СОБ (кш0/) +
кт0р0 д
12 дЭ
в3 X, біп I—
Лк
Я
(27)
для к = 1, 2, 3,..., 0 < ф < 2 к 91, 0 < 91 < 1, Ът = кЯ/8 < Э< кЯ/2 = Ъ^, 0 < г < в и
Рек = ОДР(к) ЄХР(Ік®0/)], Хлк =
БІп(кш0Ґ ) СОБ(кш0Ґ )
Лк
- +
Лк,
ЛкЪ
— =
Лк,
'О
_ Л0Л1
чЛк у
Лк
= лт
ЛкЪ
Г\}
чЛк у
к ш0Р 1 |2~
Лк
(28)
, Лк =Л0Л1 + '^0^ ІЛк Г =(Л0Л1 )2 +(к®0р)2.
Уравнение (27) определяет поправки Рск к давлению, вызванные нестационарным течением синовиальной жидкости с вязкоупругими свойствами. Новые амплитуды имеют вид (см. Приложение 2):
‘‘"’Л». - (29)
и * = и*А V V * = ^ V ш * =
^ ік к5 іик’ ік к5 іхк* к к
где I = ф, Э и к = 1, 2, 3,... Функции Бык, Б^, $гк зависят от частот шфу, шфМ, о^, »»«, »г (см. Приложение 2). В частном случае ш^ = »¡м = ш0 и шг ^ ш0 находим 8ик = 1, Sivk = 1, следовательно, ик = ик, Ук = Ук и 8гк ^1, где 0 < Ве << 1 (см. Приложение 2, (А2.18)). Уравнения (27) при к = 1, 2, 3,... дают окончательные модифицированные стохастические уравнения Рейнольдса для поправок [3]:
1
Я біп
__5[ Ев3) дЕ(Рек) 1+ Я Ев3) дЕ(Рек) 8ІП Э1 =
Эдф| Лк, дф і дЭ| Лк, дЭ Яі Я
= -12т ЯВ, ібіпЭ|в»>^ 5ІП(^> -12ІУ дЕ<В> . Э V дЕ(в)'
. Э
I—
Я
}гк
к 4
СОБ (кш0/) х-----+ 6
к
т10ЯI б1п ЯЯ 1 + Кр0
ф0 дф + ЯIБ1П Я/а0 дЭ
дЕ(в) СОБ(кт/) | Ю0Р0 А[Е(в3)х
\ I— Л
(30)
+6Я [Ш30 Я (ЯП ф)+К>0
_д_
дЭ
Е(в)б1п— Я
дф к5 12к4 дф
СОБ (кш0/)
к5
|
т0р0 д
12к 4 дЭ
Лк
Э
|
Е(в )ХЛк81П-
Я
где 0 < ф < 2к91, 0 < 91 < 1, Ьт = кЯ/8 < Э < кЯ/2 = Ъц, 0 < г < в. При тг = ш0 получим Сгк = 1, где 0 < В <<1. Уравнение (30) определяет средние значения для поправок к давлению Е(Рек). Среднее значение суммарного давления и оператор осреднения имеют следующий вид [3]:
2
1
х
а=л/2
нагрузка
давление
Рис. 4. Расположение сил трения в окружном и меридиональном направлениях
E (p) = E[ p '”>] + £ E[ Pck ],
k=1
E(*> =J (*) yf (8)d8, f (8).
35 32c7 0
(c2 - 82) при - c < 8 < +c, при 18 > 9,
(31)
где / - плотность распределения случайной шероховатости хряща. Через е обозначена половина диапазона изменения случайной толщины. Функция / обращается в нуль при е = ± 3а, где а - стандартное отклонение. После вычислений находим:
+да +да -і
Е(в) =|є/(5)^8 = Г,6'0', Е(є3) = |є3/(5)аб = (гу0>)3(1+3р2), 0= - (32)
—да —да 5
где Г5 = 1+Гу. На рис. 4 показано распределение давления р(ф, Э) на сферической поверхности головки бедра (ф и Э - сферические координаты). Полная нагрузка на головку бедра получается следующим образом:
С,., = // Е[р(ф, Э)]ап(ф, Э), (33)
п(ф,э)
где Я - интеграл по поверхности головки бедра. Элемент площади в интеграле имеет
п
вид:
d Q =
Я r Я r
^*0 ^ ^*0 X
Яф ЯЭ
d ф d Э при r0= ix+jy + kz ,
(34>
r=R
где
x = r cos ф sin($/R), y= r sin ф sin($/R), z = r cos($/R) 0 < r< R. (35)
Через r0 обозначен радиус-вектор точки поверхности головки бедра (см. рис. 4); 0 < ф < 2лть 0 < ci< 1, tcR/8 < $< tcR/2. Подставляя (35) в (34), получим:
dQ =R si n(Э/R>dфd(Э/R>.
(36)
На рис. 4 показано расположение сил трения в направлениях ф и д. Окружная и меридиональная компоненты силы имеют вид:
Л
□(Ф,а)'
Л,
¿О
д r
dSX = jj I ^
s 0(ф,&) \
dO
д r
d 0(ф, $).
(37)
У r=s
6. Экспериментальные измерения
Плотность распределения высоты суставного зазора находилась путем экспериментальных измерений шероховатости образцов хряща и соответствующих стандартных отклонений. Измерения выполнялись с помощью лазерного микродатчика, установленного на аппаратуре Rank Taylor Hopson Talyscan 150 и управляемого посредством компьютерных программ Talyscan Expert и Microsoft Exel. Измерения проводились для 29 образцов; получены параметры St, Sz, Sa шероховатости поверхности (в мкм). Здесь: St - разность максимальной высоты и максимальной глубины на поверхности головки бедра, Sz - среднее арифметическое пяти максимальных высот и пяти максимальных глубин, Sa - стандартное отклонение для плотности распределения шероховатости поверхности хряща. Измеренные значения St находятся в интервале от 9,79 мкм до 24,7 мкм. Измеренные значения Sz находятся в интервале от 8,52 мкм до 14,7 мкм. Окончательно вычисленные значения Sa находятся в интервале от 0,78 мкм до 1,96 мкм. На рис. 5 показана карта поверхности нормального суставного хряща с высотой неровностей до 24,7 мкм и сечения 1 и 2 этой поверхности. Образец имел 2 мм в длину и 2 мм в ширину.
7. Численные расчеты
Распределение давления р(0) и поправки к нему р(1), р‘\ р‘3),... определялись в области смазывания О. Суммарное давление равно атмосферному рл на границе области О, показанной на рис. 3 и на основе медицинской информации, определяемой неравенствами О: 0 < ф < л, 'кRJ8 <а3 = д <лЯ/2. Это сечение сферического углубления. Численные расчеты проводились для области О, R = 0,0265 м, ш1 = 1,40 с-1, ш3 = -0,45 с-1, ш0 = 500 с-1, Ш10 = 0,100 с, Ш30 = 0,025 с, Дбх = 2,5 ? м, Дбу = 0,5 ? м, Де2 = 2,0 ? м, п0= 0,15 Па, р0=1000 кг/м3. Минимальное значение высоты зазора равнялось вт;п = 4,8 дм, максимальное значение высоты зазора равнялось втяУ = 10,50 дм. Учитывалась измеренная шероховатость поверхностей. Значения давления вычислялись в следующие моменты времени: t = 0 с, t = 0.3л/ш0 с, t = л/ш0 с, t = 1.7л/ш0 с, t = 2л/ш0 с,.... Вначале были проведены численные расчеты без учета вязкоупругих свойств синовиальной жидкости, то есть при р0 = 0,00000 Па2 и р0 = 0,00020 Па2.
На вертлужной впадине принимались следующие амплитуды колебаний касательной скорости: ¥ф0 = 0,001 м/с, ¥д0 = -0,0002 м/с. Коэффициент возмущений высоты зазора принимался равным Ве = +0,0002. Принимались во внимание случайные эффекты для 0 < р < 1.
Численное нахождение зависящего от времени распределения давления проводилось с помощью уравнения (30), а высоты зазора с помощью (24, 25). Результаты показаны на рис. 6, 7, 8, 9. Временной период возмущений составлял t = 2л/ш0.
цм
24,7 цм
2 мм
длина образца
длина образца
Рис. 5. Измеренная поверхность нормального хряща в тазобедренном суставе
и два ее сечения
На рис. 6 показано изменение со временем распределения давления для р = 0 и Р0 = 0. Стохастическое описание шероховатости хряща не принималось во внимание, так как р = 0. Вязкоупругие свойства синовиальной жидкости не учитывались, так как Р0 = 0,0000 Па2.
Для моментов времени / = 0, / = 0,3к/ш0, / = гс/ш0, / = 1,7к/ш0, / = 2к/ш0 получены следующие значения нагрузки: 1292,7 Н, 1655,0 Н, 1306,5 Н, 936,7 Н, 1292,7 Н.
На рис. 7 показано изменение со временем распределения давления для р = 0 и Р0 = 0,00020 Па2. Стохастическое описание шероховатости хряща не принималось во внимание, так как р = 0. Вязкоупругие свойства синовиальной жидкости учитывались, так как Р0 = 0,00020 Па2.
Для моментов времени / = 0, / = 0,3к/ш0, / = к/ш0, / = 1,7к/ш0, / = 2к/ш0 получены следующие значения нагрузки: 1289,0 Н, 1839,1 Н, 1309,2 Н, 749,5 Н, 1289,0 Н.
Рис. 6. Распределение давления в зазоре тазобедренного сустава, вызванное вращением сферической головки бедра одновременно в окружном ф-направлении и меридиональном &-направлении. Случайные эффекты и вязкоупругость синовиальной жидкости не учитываются, так как р = 0 и р0 = 0. Вычисления проведены для ненулевых значений угловых скоростей ю1 = 1,40 с-1 ,©3 = - 0,45 с-1, ненулевых значений возмущений угловых скоростей ©10, ©30 при нестационарном течении и для частоты возмущений высоты зазора ©0 = 500 с-1, Я = 0,0265 [м], ц = 0,15 [Пас], ©0 = 500 [1/с], ©1 = 1,4 [1/с], ©з = 0,45 [1/с],
©10 = 0,1 [1/с], ©зо = 0,025 [1/с], р0 = 0, р = 0, смазываемая поверхность = 20,38 [см2]
На рис. 8 показано изменение со временем распределения давления для р =1/3 и р0= 0. Стохастическое описание шероховатости хряща принималось во внимание, так как р = 1/3. Вязкоупругие свойства синовиальной жидкости не учитывались, так как р0 = 0,0000 Па2.
Для моментов времени t = 0, t = 0,3л/ш0, t = л/ш0, t = 1,7лУю0, t =2л/ю0 получены следующие значения нагрузки: 1258,1 Н, 1612,6 Н, 1270,8 Н, 909,5 Н, 1258,1 Н.
На рис. 9 показано изменение со временем распределения давления для р = 1/3 и р0 = 0,00020 Пас . Стохастическое описание шероховатости хряща принималось во внимание, так как р = 1/3. Вязкоупругие свойства синовиальной жидкости учитывались, так как р0 = 0,00020 Па .
Для моментов времени t = 0, t = 0,3я/ю0, t = л/ю0, t = 1,7я/ю0, t = 2я/ю0 получены следующие значения нагрузки: 1254,7 Н, 1792,8 Н, 1273,4 Н, 726,4 Н, 1254,7 Н.
Легко видеть, что распределения давления и нагрузка, представленные на рис. 6,
7, 8, 9 для моментов времени t = 0 с и t = 2л/ш0 с, то есть для t = 2кл/ю0 с и к = 0, 1, 2,..., имеют одинаковые значения. Распределения давления и нагрузка в моменты времени t = (2к - 1)л/ю0 с для к = 1, 2, 3, 4,... также имеют одинаковые значения.
=1,798-106 Псі
т-1,5 [МПа]
-0,75 [МПа]
=2,443-10° На Сс=1839,1 Н
г 1,5 [МПа]
.0,75 [МПа]
*=я/©о с
р^Ы^-Ю6 1Ла С,=1308,2 Н
¿=1,771/о0 с Ршах=1,176-106 Па См-=749,5 Н
Рис. 7. Распределение давления в зазоре тазобедренного сустава, вызванное вращением сферической головки бедра одновременно в окружном ф-направлении и меридиональном 0-направлении. Случайные эффекты не учитываются, так как р = 0; вязкоупругость синовиальной жидкости учитывается, так как р0 = 0,00020 Па2. Вычисления проведены для ненулевых значений угловых скоростей о1 = 1,40 с-1 , о3 = -0,45 с-1, ненулевых значений возмущений угловых скоростей о10, о30 при нестационарном течении и для частоты возмущений высоты зазора о0 = 500 с-1. Я = 0,0265 [м], р = 0,15 [Пас], о0 = 500 [1/с],
©і = 1,4 [1/с], ю3 = 0,45 [1/с], о10 = 0,1 [1/с], о30 = 0,025 [1/с], р0 = 0,00020 [Пас2], р =0,
смазываемая поверхность = 20,38 [см2]
Верхняя левая иллюстрация на рис. 6-9 показывает распределение давления в начале и в конце периода возмущений движения сустава. Иллюстрации, отвечающие ї = л/ш0 на рис. 6-9, показывают распределение давления в середине периода возмущений движения сустава. Затем распределения давления возвращаются к показанному на первой иллюстрации.
Влияние случайных изменений высоты зазора и вязкоупругих свойств синовиальной жидкости на нагрузку показано на рис. 10 и обсуждается ниже.
На рис. 10 показано изменение нагрузки со временем в течение периода возмущающих эффектов (частота о0 = 500 с-1) для четырех вариантов условий, отвечающих четырем случаям на рис. 6, 7, 8, 9: (1) р = 0, р0 = 0; (2) р = 0, р0 = 0,0002 Па2; (3) р = 1/3, р0 = 0; (4) р = 1/3, р0 = 0,0002 Па2. Если р = 0, то случайные эффекты не учитываются. Если р0 = 0, то синовиальная жидкость является ньютоновской.
На рис. 11 (а) показано изменение нагрузки со временем в течение периода возмущений при трех различных частотах возмущений: о0 = 100 с-1, о0 = 500 с-1, о0 =
1000 с-1 для ньютоновской (р0 = 0,00000 Па2) и вязкоупругой (р0 = 0,00020 Па2) синовиальной жидкости.
На рис. 11 (Ь) показано изменение нагрузки со временем в течение периода возмущений (ш0 = 500 с-1) для следующих значений параметра р0, описывающего вязкоупругие свойства синовиальной жидкости: р0 = 0,00001 Па2, р0 = 0,00005 Па2, р0 = 0,00010 Па2, Р0 = 0,00015 Па2, р0 = 0,00020 Па2, р0 = 0,00025 Па2, Р0 = 0,00030 Па2, Р0 = 0,00035 Па2, р0 = 0,00040 Па2.
Размеры неровностей на поверхности хряща лежат в интервале от 0,2 цм до 0,9 цм. Поэтому период 2я/ш0 возмущений смазывающего течения синовиальной жидкости может достигать больших значения, и частоты возмущений в этом случае очень малы.
Г 1.5 [МПа]
0.75 [МПа]
t=0 , /=2т1/со0 с Ртах=1,777-106 Па Ctot=1258,1 Н t=0,377/CD0 С Ртах=2Д93-106 Па С,0^1612,6 Н
-1,5 [МПа]
-0,75 [МПа]
^-0
t= 71/С00 С
ртах=1,786-106 Па Сс=1270,8 Н
1,7^7./ Ю 0 С Ртах=1,369-106 Па С,0,=909,5 Н
Рис. 8. Распределение давления в зазоре тазобедренного сустава, вызванное вращением сферической головки бедра одновременно в окружном ф-направлении и меридиональном &-направлении. Случайные эффекты учитываются, так как р=1/3; вязкоупругость синовиальной жидкости не учитывается, так как р0 = 0,00000 Па2. Вычисления проведены для ненулевых значений угловых скоростей «1=1,40 с-1 , «3= -0,45 с-1, ненулевых значений возмущений угловых скоростей й10, «30 при нестационарном течении и для частоты возмущений высоты зазора «0=500 с-1. Я = 0,0265 [м], ц = 0,15 [Пас], «0 = 500 [1/с], «1 = 1,4 [1/с], «3 = 0,45 [1/с], й10 = 0,1 [1/с], «30 = 0,025 [1/с], р0= 0, р=1/3, смазываемая
поверхность = 20,38 [См2]
1,5 [МПа]
-0,75 [М Па]
І-0
ї=0, ї=27і/(»о с ртаХ=1,775-106Па С,0,=1254,7 Н
¿=0,37 <Ю0 с ртаХ=2,407-106Па С,„,=1792,8 Н
1,5 [МПа]
-0,75 [МПа]
■*-0
ї=7і/С00 с
:^1,7^7'8-106 Па С„=1273,4 Н
¿=1,777/ю0 с Ртах=1,159-106 Па С„=726,4 Н
Рис. 9. Распределение давления в зазоре тазобедренного сустава, вызванное вращением сферической головки бедра одновременно в окружном ф-направлении и меридиональном &-направлении. Случайные эффекты учитываются, так как р = 1/3; вязкоупругость синовиальной жидкости учитывается, так как ро=0,00020 Па2. Вычисления проведены для ненулевых значений угловых скоростей «1 = 1,40 с-1 , «3 = - 0,45 с-1, ненулевых значений возмущений угловых скоростей й10, «30 при нестационарном течении и для частоты возмущений высоты зазора «0 = 500 с-1. Я = 0,0265 [м], ц = 0,15 [Пас], «0 = 500 [1/с], «1 = 1,4 [1/с], «3 = 0,45 [1/с], «10 = 0,1 [1/с], «30 = 0,025 [1/с], р0 = 0,00020 [Пас2],
р = 1/3, смазываемая поверхность = 20,38 [см2]
2000
К
о
ст
1800
1600
1400
1 -р=0, Р0=0 3 -~р=1/3,Р0=0 2 ^р=0, р0=0,00020 Па2 4 р=1/3, Р0=0,00020 Па2
ю0=500 1/с
&
к
1200
1000
800
600
7/юо
2я/ю0
0 0,001 0,002 0,003 0,004 0,005 0,006 0,007 0,008 0,009 0,010 0,011 0,012 0,01 3 х, с
Рис. 10. Зависимость нагрузки от времени за период для четырех предположений, отвечающих следующим случаям: (1) р = 0, р0 = 0; (2) р = 0, р0 = 0,00020 Па2; (3) р=1/3,
р0 = 0; (4) р= 1/3, р0= 0,00020 Па2
X
О
и
а
із
д
2500
2000
1500
1000
500
27/ю0
0,06
ї, с
Рис. 11 (а). Зависимость нагрузки от времени за период возмущений для двух констант псевдовязкости р0 = 0,00000 Па2 и р0 = 0,00020 Па2 и трех различных частот возмущений
о
Я
к
2500
2000
1500
1000
500
0
ю0= 100 с , ю0= 500 с , ю0= 1000 с
ро=0,00001
р0=0,00005
в0=0.00010
р0=0,00015 •Р0=0,00020 Ро=0,00025 ю0=500 1/с
-ж-Ро=0,00030 о=0,00035 Р0=0,00040
0 0,001 0,002 0,003 0,004 0,005 0,006 0.007 0,008 0,009 0.010 0,011 0,012 0,013
X, С
Рис. 11 (Ь). Зависимость нагрузки от времени за период возмущений при частоте возмущений ю0 = 500 с-1 для различных констант псевдовязкости
0,00001 Па2 < р0 < 0,00040 Па2
0
Кривые 1 и 4, приведенные на рис. 11 (а), показывают, что для сравнительно малых значений частоты возмущений ю0 влияние вязкоупругих свойств синовиальной жидкости на суставную нагрузку пренебрежимо мало. Если же частота ю0 принимает значения от 500 с-1 до 1000 с-1, то вязкоупругие свойства синовиальной жидкости (Р0 = 0,00020 Па2) дают увеличение нагрузки в диапазоне от 15% до 60% по сравнению с нагрузкой при р0 = 0,0000 Па2 (ньютоновская жидкость).
8. Заключение
Полученное уравнение для распределения давления при нестационарных периодических условиях, с учетом случайно меняющегося зазора, вязкоупругих свойств синовиальной жидкости и различных частот и амплитуд колебаний головки бедра и вертлужной впадины приводит, в частном случае, к известной форме уравнения Рейнольдса для стационарного течения.
Из численных расчетов следует, что вязкоупругость синовиальной жидкости приводит к увеличению суставной нагрузки на 15%, а в некоторых случаях на 60%. Влияние вязкоупругости синовиальной жидкости на нагрузку очевидно связано с периодичностью возмущений в суставе. Если частота возмущений больше 100 с- , то вязкоупругость синовиальной жидкости существенно влияет на суставную нагрузку. При частотах, меньших 100 с-1, влиянием вязкоупругости на суставную нагрузку можно пренебречь.
Из численных расчетов следует, что стохастическое описание неровностей костных поверхностей и толщины пленки синовиальной жидкости изменяет суставную нагрузку на 11%. Из этих расчетов также следует, что учет смазывания при одновременном вращении головки бедра в двух направлениях позволяет установить точную картину распределения давления в суставе.
Понимание механизма смазывания сустава синовиальной жидкостью с вязкоупругими свойствами весьма важно. Например, такое понимание может существенно улучшить диагностику заболеваний суставов.
Приложение 1
Подставим бесконечные ряды (9) - (12) в систему уравнений (4) - (7) и приравняем члены с одинаковыми верхними индексами в скобках и одинаковыми степенями экспоненциальных функций.
Приравнивая члены с верхним индексом нуль, получим следующие уравнения стационарного движения ньютоновой жидкости:
Л 1 г^ф(0) я Г аУ°Л
0 = —— СОБеС | — |—--ъц0 —
я
Я) дф
дт
Лг
дт
0 =
др
(0)
дт
др(0) д ґ
о =-----------+ Л0— Лі
дЗ дт , 1
д”ф0) Н я • н д
—— + Я БІДІ — I—— Н----------------------
I Я ) дт дЗ
дф
д_ д т
дуЗ0)^ дт 'З'
я 8ІПІЯI уЗ0)
= о,
к-
дТ(0) д т
\
Ґ
ду!0) У Гду(0) V
+ '
дт
V дт )
= 0,
(А1.1)
(А1.2)
(А1.3)
(А14)
(А1.5)
для 0 < ф < 2л01, 0 < 01 <1, Ьт = тЯ/8 < 3 <тЯ/2 = 0 < г < в.
Система (А1.1) - (А1.5) определяет искомые давление р(0), температуру Т-0) и
(0) (0) (0) п
компоненты скорости уф , Уг , у3 по направлениям ф, г, 3.
Приравнивая члены с верхним индексом к = 1, 2,..., получим следующие уравнения движения для поправок к-го порядка, учитывающих нестационарность течения и вязкоупругие свойства жидкости:
Г яук) >
1
(к)
ikш0р0уф) =-----cosec I — ,
0 0 ф R | R J
а +— dr
Лк‘
dr
(A1.6)
0 =
dp
(к)
dr
dp(k) d
----------1----
d3 dr
Лк
) H R ■ i»W) h d
—— + R sin I — I—— H--
dф v R J dr d3
dvf) ^ dr
З
R sin IR J v3k)
= 0,
(A1.7)
(A1.8)
(A19)
для к = 1, 2, 3,..., 0 < ф < 2лть 0 < с1 < 1, Ьт = лЯ/8 < З <лЯ/2 = 0 < т < в; лк =
ЛоЛ1+і^ЮоР.
Система (А1.5) - (А1.8) определяет искомые поправки р(к) к давлению и поправки V,
(k) v (к)
ф 5 5
у3(к) к компонентам скорости. Кроме того, получаем уравнение энергии для определения поправок Тк) к температуре. Верхние индексы к = 1, 2, 3,... обозначают порядок поправки.
Приложение 2
Проведем преобразование бесконечного ряда периодических функций с известными амплитудами Aimk и известными частотами oim в бесконечный ряд с новыми неизвестными амплитудами A*mk и известной частотой ш0 (колебаний в синовиальной жидкости). Это преобразование проводится следующим образом:
ад ад
Z A*mk exp (ik^) = Z Aimk exP (^Ш^) , (A2.1)
k—1 k=1
где в частном случае i = ф, m = u,v; k =1, 2, 3,... и
AUk - Uk, Aik - С Auk - Uik, 4* - . (A2.2)
Уравнение (A2.1) описывает следующие случаи:
ад ад ад ад
ZUi exPOw)=ZU# exP(ik^), exPOw)=2^ exP(ik^),
k—1 k—1 k—1 k—1 / а ч
ад ад ад ад
ZU3k exP (ik^) = ZU3k exP (ik^) , ZF»k exP (ik^) = ZF3k exP (^Ш^) .
k—1 k—1 k—1 k—1
Раскладывая экспоненциальную функцию в обеих частях (A2.1), получим для левой части (A2.1):
^1(1 Н ^ ^ - і ^ Н..............................) + ■■.. Н
2 6
\3
+4т2(1 + 2І®0* - 2(Ю0^)2 - І ( 30 ) +.............) + .... + (А2.4)
,<* /•, 9(ш0ґ)2 . 9(ш0ґ)3 .
+ЛШ3(1 + 3^ - і^^~ +.) +....
Для правой части (А2.1) получаем:
(®гт*)2 : КтО'
4т1(1 + і®гт* +.......) +.... +
26
+4т2(1 + 2іЮгт* - 2(®гт*)2 - і ^ ^ +....) + .... + (А2.5)
+4т3(1 + 3І®гт* - - і ) +....) + ....
Приравнивая обе части (А2.4) и (А2.5), получим систему линейных
алгебраических уравнений:
4т1 Н 4т2 Н 4т3 = 4т1 Н 4т2 Н 4т3,
4т1 + 2 Ат2 + 34т3 = '°гт (4т1 + 2 Ат2 + 34т3Х (А2.6)
1 А* + 2 А* +9 А* = О2(° А + 2 А +9 А )
2 іт1 т2 “ 4 Лгт3 т V ^ Лг' т1 “ ^т2 “ 4 Лг т3 Л
где
Огт = ~, (А2.7)
®0
для і = ф, З; т = и, V.
Решая систему (А2.6), найдем искомые амплитуды:
12 3
4т1 = ^ 4т 1 (2 - 5гт )(3 + 2Огт ) + ^ Ат2 (1 - 5гт )(3 + 4Огт ) + ^ 4т3 (1 - 5гт )(2 + 3Огт X
1 1 3
Ат2 = ^ 4т1 (1 - 0гт )(3 - 4Огт ) + ^ Ат2 (1 - 2Огт )(3 - Цт ) + ^ 4т3 (1 - 0гт )(1 - 6Огт X (А2. 8)
* 2 4 1 2
4т3 = - ^ 4т1 (1 - ^т )(2 - 0іт ) - ^ Ат2 (1 - 20гт )(1 - 0іт ) - ^ 4т3 (4 - 1 Цт + 90гт ),
где і = ф, З; т = и, V. Например, символ А*гтк для к = 1, 2, 3,... описывает следующие амплитуды:
А* = ТІ * А* = ТІ* А* = ТІ*
^фи1 ^ф1’ ^фи 2 ^ф2’ ^фи 3 ф3 ’......................
А* = И * Аа = И * А* = И *
У1Зи1 З1 ’ Зи2 З2 ’ Зи3 (~/З3’-"- (А2 9)
Т7-* Т7-* Л* Т7-* ^ " '
А = V А = V А = V
^хфу1 у ф^ ^фу2 у ф!^ фу3 ^ф3’.
А* = V * А* = V * А* = V *
У1Зу1 З1 ’ Зу2 уЗ^ Зу3 ^ З3 ’--------
Из (11) и (А2.2) для і = ф, З; т = и, V следует:
и.
А = и =
Лфи1 ^ф1 15
и
фД
и
фД
А„„, = ий1 =
3Д
15
А = и =
фи 2 ф2 25
А = и = и
Л3и 2 32 25
и
фД
3Д
А = и =
фи 3 ф3 35
А = и = и
3и 3 ^33 ~5
3Д
(А2.10)
V V V
А = V = ^0 А = V = ^0 А = V = ^0
фу1 *ф1 15 5 фу2 *ф2 25 5 фу3 *ф3 35 5...
А = V = А = V = --3о А = V = ^
3у1 '31 15 ’ 3у2 32 2^ 3у3 ^ 33 35
Например, амплитуды колебаний в окружном ф и меридиональном 3-направлениях имеют следующий вид:
1 2 3
ип = "5 иг1 (2 - Я ги )(3 + 2Яги ) + иг2(1 - Яги )(3 + 4Яги ) + ^ ^'3 (1 - Я ги )(2 + 3Яги ),
1 1 3
и* = 5и (1 - Яги )(3 - 4Яги) + 5и2(1 - 2^ )(3 - 8Яги) + 5 игЪ (1 - Яги )(1 - 6яы ), (A2.11)
2 4 1
иг3 = - ^ иг1 (1 - Яги )(2 - Яги ) - ^ иг 2(1 - 2Яги )(1 - Яги ) - ^ ^'3 (4 - 1 8Яги + 9Яги X
1 2 3
^г 1 = ^ ^г 1(2 - Ягу )(3 + 2ягу ) + ^ ^г 2(1 - Ягу ) (3 + 4Ягу ) + ^ ^3(1 - Ягу )(2 + 3Ягу X
1 1 3
^г2 = 5 ^г 1(1 - Ягу )(3 - 4ягу ) + ^ Vг 2(1 - 2Ягу )(3 - 8Ягу ) + ^ Vг3(1 - Я гу )(1 - 6Ягу X
2 4 1
V 3 =- 5 V 1(1 - Ягу )(2 - Ягу ) - 5 V 2(1 - 2Ягу )(1 - ) - 5 V 3(4 - Ш* + 9^ ),
где г = ф, 3. Используя (А2.10)-(А2.11), найдем:
и V
и* = гД С V* = ф0 С
к ~ к5 гк ~ к 5
жк ’
(А2.12)
для г = ф, 3 и к = 1, 2, 3,., где:
Ст1/15 =
Ст2/25 = 5
Ст3,/3! = 5
2 3
(2 - Ягт )(3 + 2я гт ) + ^5 (1 - Ягт )(3 + 4Ягт ) + 3Т (1 - Ягт )(2 + 3Ягт )
1 3
(1 - Ягт )(3 - 4Ягт ) + (1 - 2Я гт )(3 - 8Ягт ) + 3? (1 - Я гт )(1 - 6Ягт )
4
1
-2(1 - Ягт )(2 - Ягт ) - ТТ (1 - 2я гт )(1 - Ягт ) - ТТ (4 - 1 8Ягт + 9Ягт )
2
35
(А2.13)
где т = и, V и г = ф, 3. Для шгт = ш0 получим Ягт = 1.Следовательно, 8гтк= 1 для т = и, у и
* *
игк = игк, Vгk = Vгk. Из медицинских экспериментов известно, что 0,95 <Ягтк <1.05. Поэтому следующие неравенства выполнены:
С
С
С.,
+0,964 <+1,034, + 0,016 <^т52 <+0,048, -0,006 <^т3 <+0,021,... (А2.13)
1
2
35
Применяя аналогичный алгоритм вычислений для температуры (14)2 и для выражения (22), получим решение, аналогичное (А2.8). Искомые амплитуды для компоненты скорости в направлении высоты зазора имеют вид:
где
1 2 3
W = 5 Wi(2 - s )(3 + 2sr) + 5 W2(1 - S )(3 + 4^.) + 5 W3(l - S )(2 + 3^.),
1 1 3
W = 5 W1(1 - s )(3 - 4^.) + 5 W2(1 - 2яг )(3 - 8^) + 5 Wi(1 - s )(1 - 6^),
2 4 1
W3 =- 5 W1(1 - s )(2 - s) - 5 W2(1 - 2^. )(1 - s) - 5 W3(4 - 18яг + 9s2),
(A2.14)
где
®г
Sr = —
ш0
Учитывая, что s(k) = Bes(0)/k5, находим из (23):
(A2.15)
W = в(1)ш = Bs(0V,, W = 2s(2)® = 2B s(0V / 25 = -V B в(0)ш
(0)
r e r’ 2
r er
2
4 e r ’
W = 3s(3V = 3Д,(0V/35 = -1-Bs(0)Qr,...., Wk = -1£es(0y.
34 r k k e
r er
(A2.16)
Тогда из (A2.12) - (A2.14) следует, что
,(0)j
,(0)j
14
24
k4
v14=5
2
3
^=5
* 5/з4=5
(2 - sr )(3 + 2sr ) + 24 (1 - Sr )(3 + 4sr ) + 3^ (1 - sr )(2 + 3sr )
1 3
(1 - sr )(3 - 4sr ) + 24(1 - 2sr )(3 - 8sr ) + 3^(1 - sr )(1 - 6sr )
4
1
-2(1 - s)(2 - s) - — (1 - 2Sr)(1 - s) --(4 - 18^ + 9sr2)
r 34
(A2.18)
Для or = ш0 находим sr = 1 и, следовательно, Srk = 1. В этом случае Wk = Wk и
Tk*= Tk.
Теперь выведем выражение для высоты гладкого зазора s(0). Соотношения между прямоугольными (х, y, z) и сферическими (ф, r, 3 ) координатами (см. рис. 3) имеют вид:
х = r 8Іп(31)ео8ф, y = r sin^^si^, z = r cos(31), 0 <r<R, 3^3/R. (A2.19)
Графически положения центра сферической головки бедра 0(0, 0, 0) и центра сферической вертлужной впадины 01(x-As1, y-As2, z+As3) показаны на рис. 3. Уравнение поверхности сферической вертлужной впадины в окрестности центра 01(x-As1, y-As2, z+Аєз) можно записать в виде:
(x-As02+ (y-As2)2+ (z+A83)2=(R+D+8 min)2, D = [ (AS1)2+(AS2)2+(AS3)2]0,5. (A2.20)
Подставляя (A2.19) в (A2.20), получим:
(r cosф sin31-As1)2+(r sinф sin31-As2)2+(r cos31+As3)2=(R+D+smin)2. (A2.21)
Высота зазора
в(0)(ф,30 = г-Я. (А2.22)
Найдем г из (А2.21) и подставим эту величину (А2.22). Тогда получим окончательное выражение для высоты зазора в виде (22).
Приложение 3
Предположим, что динамическая вязкость синовиальной жидкости не меняется в направлении высоты зазора. Интегрируя (А1.1) и (А1.3) дважды по переменной г1 и учитывая граничные условия (16), (18), получим компоненты скорости синовиальной жидкости в окружном и меридиональном направлениях в виде:
1 1 f-Э^ф(0)
v„ —-----COSeC I — I--------s s (i — s) + и ф0(
ф 2ВДR IR) а„ ' ' ""
v<0) — — ^--|-cosecj ^ \^— s'2s(1 — s) + U„,(1 — s), (A3.1)
1 аР(0)
v»0) — —~--------s-s(1 — s) + UЭ0(1 — sX (A3.-)
2Л0Л1 03
где s = r/s, 0 < r <s, 0 < ф < 2л01, 0 < 01 < 1, tcR/8 < 3 < tcR/2.
Интегрируя (A1.6), (A1.8) дважды по переменной r1 и учитывая граничные условия (19), получим поправки к вышеуказанным компонентам скорости, вызванные нестационарностью течения и вязкоупругими свойствами жидкости [3]:
V») — п Y + и * sinh [(s — r) A ] + у silnh (rA>) (A3 3)
Vi ikYk Uik • л / j \ ik • i / л \, (A3.3)
sinh (sAk) sinh (sAk)
Yk -[1—exp (rAk )]—[1—exp (sAk)] sinh (sAk), (A34)
п ^п ^_L_^ A ^ /ikffl0p0 (A35)
„k ko0p0R sin(3 / R) 0ф ’ * kш0р0 03 k v Лk '
где к = 1, 2, 3,..., при I = ф, 3, 0 < ф < 2л01з 0 < 01 < 1, лЯ/8 < 3 < лЯ/2, 0 < г <в.
Приложение 4
Интегрируем уравнение непрерывности (А1.4) по переменной г1 с учетом граничных условий (17). Компонента скорости синовиальной жидкости в направлении высоты зазора принимает тогда следующий вид:
v<0> — — Icosec j 3lf ^r — cosec j f A
r r j r ) J0 а„ l r \ J0 аз
v"’sin U ,
dr. (A4.1)
Интегрируем уравнение непрерывности (А1.9) по переменной г1 с учетом граничных условий (20) 1. Тогда поправки к вышеупомянутой компоненте скорости, вызванные нестационарностью течения и вязкоупругими свойствами жидкости, принимают вид:
к> 1 (-&yf0V í-&yf5
vy =-----cosec\ — И —^—dr - cosecl — І —
' R \ R) J0 Єф l R) J0 es
vs>sin (I.
dr,
(A4.2>
где к = 1, 2, 3,..., 0 < ф < 2л0ь 0 < 01 < 1, лЯ/8 < 3 < лЯ/2, 0 < г < в.
Учитывая граничное условие (18) для г = г на компоненту (А4.1), найдем:
Є av(0> є
+ R J
o аФ 0
A
as
v|0> sin l 1
R
dr = 0.
Учитывая граничное условие (23) на компоненту (A4.2), найдем:
^ + R J
О аФ 0
А
as
vSk> sin l 1
R
dr = - i Wk* R sin (1
0 < ф < 2л91з 0 < 91 < 1, tcR/8 < S < tcR/2, 0 < r < є.
Приложение 5
(A4.3)
(A4.4)
Подставим компоненты скорости (A3.1), (A3.2) в (A4.3). Тогда получим модифицированное уравнение Рейнольдса (26). Используя правило дифференцирования интеграла с переменными пределами и граничное условие (21)ь найдем из (A4.4):
¿jk*+R^k)sin (Rdr=-¡w*r sin Ш+V* | + rvrsin (Ш ■ (A51)
0 ^ 0
Подставив в (A5.1) поправки к компонентам скорости (A3.3) для i =ф, S, найдем:
А
аф
(k>
J Ykdr
+ R —
as
ар
(k>
's'
sin \ —
as l r ,
o
J Wkdr
+
+
Aj kra0p0 аф
i
sinh [(j- r > Ak ]
sinh [j Ak ]
dr+n J -
sinh (rAk) dr
sinh (J Ak)
+
T Jsinh[<s-r>Ak]d. + V. J
sk I . ! г TI sk I
o sinh [j Ak ]
sinh (rAk)
+ V a,,\----p—^dr
0 sinh [j Ak ]
(A5.2>
= -kQ0p0W;R sin \ І+ ^0
s y kШ0Р0
i
* ає * ає . ( s
_ фк аф+ sk а!sin l r ,
для к = 1, 2, 3,..., 0 < ф < 2л01, 0 < 01 < 1, лЯ/8 < 3 < лЯ/2, 0 < г < в.
Для дальнейшего упрощения (А5.2) необходимо вычислить следующие интегралы:
k ®0p0
U* J
sinh [(j- r > Ak ] sinh [j Ak ]
dr =
k ®0p0
i
U* J
ÁJ-r) Ak
- e
-(j-r)Ak
Є Ak -Є Ak
e - e k
dr =
k^°p° U* tanh í^-
J-12 j3 Ak + 0 (j4 )
(A5.3>
, при i = ф, s.
ЄГАк _ Є~ГАк
QZAk _ Є_ЄAk
dr =
V *
= s +—tanh i^Ak-
A I 2 .
кш0р0 S sinh rAk
i s
кШ0Р0 _ О (s4)
(A5.4)
12 Лк
кш0р0 Г erAk _ e_Ak
!М_0 Гsinh rAk dr = v* 0P0 Г
i * sinh s Ak lk i *
ЄЄ Ak _ e S Ak
dr =
1
= _ 2 ікю0р0^к
s_ 12s3 A>2 + O (s4)
(A5.5)
для І = ф, 3.
Подставляя (А5.3), (А5.4), (А5.5) в (А5.2), получим модифицированное уравнение Рейнольдса для поправок к давлению:
s3 др{к) Лк дФ
+ R — dS
3
(к)
S . f S | др1
— sinI — I-------------
Лк V R J дд
S
= 12i WJR sin I — I _ 12
R
тґ* дs * . f S 1 дs
+ sin I r JaS
Г* 1
' фк J
фк дф дs 1 д f s3 ^
+
дф 12 дф
i д
чЛк у
ik®oрo
+
[ • fSY 1 д [s3 'sY
s sin I — I sin — I
_ IR J_ 12 дS _Лк vR J_
+6 U; (s)+v*
+6 R [идк (ф)+VSL J
для к = 1, 2, 3,..., 0 < ф < 2л0і, 0 < 9і < 1, тсЯ/8 < S <kRJ2, 0 < r < s.
Благодарность
(A5.6)
^00р0
Данный исследовательский проект был поддержан Ассоциацией по передаче знаний им. Марии Кюри в рамках Шестой Программы Европейского Сообщества; номер контракта - МТКО-СТ-517226.
Список литературы
1. Dowson, D. Bio-Tribology of Natural and Replacement Synovial Joints / D. Dowson, V.C. Mow, A. Ratcliffe, S.L-Y. Woo // Biomechanics of Diarthrodial Joint, New York, Berlin, Londyn, Paris, Tokyo, Hong Kong. - 1990. - Vol. 2. - P. 305 - 345.
2. Fung, Y.C. Bioviscoelastic Solids / Y.C. Fung // Springer Verlag, Biomechanics, Mechanical Properties of Living Tissues. - 1993. - P. 242 - 230.
3. Kaliski, S. Drgania i fale w cialach stalych / S. Kaliski // Warszawa, PWN. - 1966 (in Polish).
4. Lin, J.-R. Surfaces Roughness Effect on the Dynamic Stiffness and Damping Characteristics of Compensated Hydrostatic Thrust Bearing / J.- R. Lin // Int. J. Machine Tools Manufact. - 2000. - Vol. 40. -P. 1671 - 1689.
5. Maurel, W. Biomechanical Models for Soft Tissue Simulation / W. Maurel, Y. Wu, D. Thalmann // Springer Verlag, Berlin Heidelberg. - 1998.
6. Mow, V.C. Basic Orthopedic Biomechanics / V.C. Mow, G.A. Atesian // Philadelphia, Raven Publishers. -1997.
7. Mow, V.C. Fluid transport and mechanical properties of articular cartilage / V.C. Mow, M.H. Holmes, W.M. Lai // Journal of Biomechanics. - 1984. - Vol. 17. - 337-394.
8. Mow, V.C. Biomechanics of Diarthrodial Joints / V.C. Mow, A. Ratcliffe, S. Woo // Berlin, Heidelberg, New York, Springer Verlag. - 1990.
9. Mow, V.C. Friction, Lubrication and Wear of Diarthrodial Joints / V.C. Mow, L.J. Soslowsky // Raven Press. Basic Orthopedic Biomechanics. - 1991. - P. 254 - 291.
10. Mow, V.C. Cell Mechanics and Cellular Engineering / V.C.Mow, F. Guilak // Berlin, Heidelberg, New York, Springer Verlag. - 1994.
11. Osinski, Z. Teoria drgan / Z. Osinski // Warszawa, PWN. - 1997 (in Polish).
12. Ping, Huang. Study on thin Film Lubrication with Second - Order Fluid / Ping Huang, Zhi-Heng Li, Meng
Yong-Gang, Wen Shi-Zhu // ASME, Journal of Tribology. - 2002. - Vol. 124. - P. 547 - 552.
13. Sobczyk, M. Statystyka / M. Sobczyk // Warszawa, PWN. - 1996 (inPolish).
14. Steinhagen, J. The pathophysiology of cartilage diseases / J. Steinhagen, B. Kurz, O. Niggemeyer, J. Bruns
// Ortopedia Traumatologia Rehabilitacja. - 2001. - Vol. 3, No. 2. - P. 163 - 168.
15. Teipel, I. The Impulsive Motion of a Flat Plate in a Viscoelastic Fluid / I. Teipel // Springer Verlag, Acta Mechanica. - 1981. - Vol. 39. - P. 277 - 279.
16. Truesdell, C.A. First Course in Rational Continuum Mechanics / C.A. Truesdell // John Hopkins University, Baltimore. - 1972.
17. Truesdell, C. Hypo-elasticity / C. Truesdell // J. Rational Mechanics and Analysis. - 1955. - Vol. 4. - P.83
- 133.
18. Ungethum, M. Tribologie in Medizin / M. Ungethum, W. Winkler-Gniewek // Tribologie
Schmierungstechnik. - 1990. - Vol. 5 - P. 268 - 277.
19. Wierzcholski K. Oil Velocity and Pressure Distribution in Short Journal Bearing Under Rivlin Ericksen
lubrication. 1998, SAMS, System Analysis Modelling and Simulations OPA Overseas Publishers. Assoc.
N.V., Vol. 32. - P. 205 - 228.
20. Wierzcholski, K. The Method of Solutions for Hydrodynamic Lubrication by Synovial Fluid Flow in Human Joint Gap / K. Wierzcholski // Control and Cybernetics. - 2002. - Vol. 31, No. 1. - P. 91 - 116.
21. Wierzcholski, K. Tragfahigkeiten fur Nichtstationare Schmierung von Menschlichen Huftgelenken / K. Wierzcholski // Expert Verlag. Tribologie und Schmierungstechnik. - 2005. - Vol. 1 (193). - P. 5 - 14.
22. Wierzcholski, K. Comparison Between Impulsive and Periodic non-Newtonian Lubrication of Human Hip Joint / K. Wierzcholski // Polish Academy of Sciences. Engineering Transactions. - 2005. - Vol. 53, No.1.
- P. 55 - 100.
23. Wu, J.Z. Artificial Joint Mechanics with Biphasic Cartilage Layer Under Dynamic Loading / J.Z Wu, W. Herzog, M. Epstein // Journal of Biomechanics Engineering. - 1998. - Vol. 120. - P. 77 - 84.
24. Xiao, H. Hypo-Elasticity Model Based Upon the Logarithmic Stress Rate / H. Xiao, O. T. Bruhns, A. Meyers // J. Elast. - 1997. - Vol. 47. - P. 51 - 68.
NON-ISOTHERMAL STOCHASTIC LUBRICATION OF HUMAN HIP JOINT IN PERIODIC MOTION WITH VARIOUS FREQUENCIES AND AMPLITUDES
K.Ch. Wierzcholski (Gdynia, Poland)
This paper presents a concept of analytical and numerical calculations of pressure and capacity for hydrodynamic lubrication of human hip joint gap in non-isothermal unsteady periodic motion. Here dynamic vibrations with different frequencies and amplitudes in the bonehead, acetabulum and synovial fluid are presented. The roughness of the joint surfaces is taken into account. The measurements of the cartilage surfaces are performed by the laser sensor. We assume that the spherical bonehead in human hip joint moves in two directions namely in circumferential and meridional ones. Synovial fluid flow is described by the equations of conservation of momentum and the continuity equation. Numerical calculations are performed in Mathcad 11 Professional Program, by using the method of finite differences.
This method assures stability of numerical solutions of partial differential equations and gives proper values of pressure and capacity forces occurring in human hip joints.
Key words: human hip joint, hip joint lubrication, periodic motion, various frequencies and amplitudes.
Получено 25 октября 2005
