УДК 531/534:57+612.7
Российский
Журнал
Биомеханики
www. biomech. ас. ru
МЕТОД ОПИСАНИЯ ТЕЧЕНИЯ И ОПРЕДЕЛЕНИЯ РЕОЛОГИЧЕСКИХ КОНСТАНТ ВЯЗКОПЛАСТИЧЕСКИХ БИОМАТЕРИАЛОВ. ЧАСТЬ 1
© 2002 г. С.Л. Гавриленко*, Р.А. Васин**, С.В. Шилько*
*Институт механики металлополимерных систем им. В.А. Белого НАН Беларуси, ул. Кирова, 32а, 246050, Гомель, Беларусь, e-mail: Shilko_mpri@mail.ru
**Институт механики Московского государственного университета, Воробьевы горы, 119899, Москва, Россия, e-mail: Vasin@imec.msu.ru
Аннотация. Получены определяющие уравнения и выполнена идентификация моделей биотканей, демонстрирующих вязкопластические свойства при наличии существенной нелинейности скоростной чувствительности. В качестве примера рассматривается определение реологических констант крови в условиях течения Куэтта и Пуазейля. Развивается приближенный метод расчета характеристик течения на основе предположения об «эквивалентной вязкости». Дана оценка устойчивости алгоритма при наличии погрешности исходных данных.
Ключевые слова: кровообращение, вязкопластичность, неньютоновская жидкость, эквивалентная вязкость, течение Куэтта и Пуазейля.
Введение
К вязкопластическим биологическим средам, демонстрирующим одновременно вязкие свойства и порог текучести, относятся как многокомпонентные жидкости, подобные крови [1,2], так и относительно твердые ткани внутренних органов. Так, по мнению авторов, вязкопластические модели могут привлекаться и для феноменологического описания формоизменения (движения) тканей и органов в условиях структурных изменений адаптационно-компенсаторного характера при силовом воздействии. К ним следует отнести костную ткань [3] и пародонт зубочелюстной системы [4], выполняющие опорную и соединительную функции. Для указанных биотканей имеет место ярко выраженный пороговый характер адаптивной реакции. Так, например, «вязкие» перемещения зуба при ортопедической коррекции, вообще говоря обусловленные процессами резорбции и остеогенеза, инициируются при определенном уровне механических напряжений. В целях биомеханического описания важно получить определяющие соотношения подобных объектов и на их основе определить соответствующие реологические константы.
Следует отметить, что изучением подобных сред занимались многие исследователи [5-11]. В большинстве случаев [8-11] полученные результаты относятся либо к модели нелинейной вязкой жидкости, имеющей нулевой предел текучести, либо описывают среду Шведова-Бингама с единичным коэффициентом скоростной чувствительности. Однако одновременный учет порогового напряжения и нелинейной вязкости биоматериалов в соответствии с механической моделью, показанной на рис. 1, очевидно, имеет принципиальное значение, что и является предметом настоящего исследования.
Рис. 1. Механическая модель нелинейной вязкопластической среды.
Анализ экспериментальных данных по реологии крови
Экспериментально показано, что кровь является существенно нелинейной вязкопластической средой [1,2]. Механическая нелинейность крови обусловлена ее сложным составом, в особенности, наличием фазы гематокрита, состоящей из
эритроцитов (красных кровяных телец) в плазме. Пластичность крови обусловлена вращением, дисковой формой, деформацией и ориентацией эритроцитов, образованием и распадом эритроцитных агрегатов в виде «монетных столбиков» в процессе течения. Неньютоновский характер течения крови, выявляемый в ротационных вискозиметрах, также объясняют агрегацией и морфологией эритроцитов. Важно, что реологические свойства крови существенно различаются при их измерении вискозиметрами куэттовского и пуазейлевского типов по сравнению с in vivo [1].
Сказанное свидетельствует о том, что течение крови неудовлетворительно описывается зависимостями, предложенными для ньютоновской жидкости с
постоянной вязкостью р = const, например, формулой Пуазейля для течения в капилляре
8 = ^, (1)
где Q - объем кровотока; r, l - радиус и длина капилляра; Ар - перепад давления; р -вязкость. В связи с этим рассмотрим более точные вязкопластические модели и возможность аппроксимации последних для описания течения биологических сред.
Зависимость напряжений сдвига т от скорости сдвига у для стационарных течений учитывает формула Кессона
Л=л/^0+kVt, (2)
где То - предельное напряжение сдвига; к - кессонова вязкость.
Метод «эффективной вязкости»
В работе [6] рассмотрен метод исследования пространственных течений вязкопластических сред с единичным коэффициентом скоростной чувствительности, основанный на введении «эффективной вязкости». Показано, что данный метод в первом приближении может быть применен для описания течений Куэтта и Пуазейля.
В настоящей работе показано, что для существенно нелинейных сред данный метод имеет невысокую точность и его следует применять с большой осторожностью.
Полные уравнения Г. Генки, описывающие пространственные течения вязкопластических сред [10], не удается решить из-за существенно нелинейной зависимости тензора напряжений от тензора скоростей деформации. Часто делают некоторые предположения относительно поля скоростей и применяют полуобратный метод, как это было проделано в работе [7], или используют замену уравнений течения вязкопластических сред уравнениями течения вязких жидкостей (уравнениями Навье-Стокса). В этом случае возникает проблема выбора коэффициента вязкости среды.
Целью данной работы было выяснение, применим ли данный метод к существенно нелинейным средам, и идентификация констант вязкопластической среды т0, п, К в условиях течения Пуазейля при известном расходе и перепаде давления.
Описание вязкопластической модели биоматериалов
Для общего теоретического описания биоматериалов предлагается следующий путь перехода от уравнений Генки к уравнениям Навье-Стокса для нелинейных вязкопластических сред.
Имеем следующие определяющие соотношения [9]:
^ =2 -%;
т = т 0 + К
(н Лт (3)
I при Т>Т0 (Н = 0 при Т<Т0),
180 )
где - девиатор тензора напряжений, ^ - девиатор тензора скоростей деформации, К - пластическая вязкость, Т0 - предельное напряжение сдвига, т - параметр скоростной чувствительности.
Здесь предлагается проводить линеаризацию следующим образом. Комплекс
(н Лт-1
Т0 + КI — I заменяется параметром це, который предлагается называть
V80 )
«эквивалентной вязкостью» рассматриваемой вязкопластической среды:
Тп ит-1
= -0 + К -Г (4)
С физической точки зрения эта замена означает, что факторы предельного напряжения сдвига и вязкости вязкопластической среды заменяются некоторой эквивалентной вязкостью. Следовательно, под эквивалентной будем понимать вязкость ньютоновской жидкости, оказывающей такое же сопротивление своему относительному перемещению, что и данная вязкопластическая среда.
Предлагается определять эквивалентную вязкость, как характеристику данного течения, следующим образом.
Будем считать, что определяемая из (4) величина це не зависит от переменных интегрирования. В этом случае система уравнений движения вязкопластической среды будет описывать течение вязкой ньютоновской жидкости с коэффициентом динамической вязкости, равным ^е. Определяем поле скоростей и давлений данного течения, а также интенсивность скоростей деформаций. Подставляя найденное значение интенсивности скоростей деформации в (4), получим в общем случае трансцендентное уравнение относительно ^е. Решая его, найдем эквивалентную
вязкость как функцию времени и пространственных координат. Используя далее различные способы осреднения, найдем цє как некоторую характерную величину процесса, не зависящую от переменных интегрирования.
Точное решение задачи о движении нелинейной вязкопластической среды
Ниже приводится точное решение задачи о движении нелинейной вязкопластической среды в круглом цилиндрическом капилляре при стационарном перепаде давления.
Для нашего материала имеем определяющие соотношения (3). Уравнения движения среды в трубе (в цилиндрической системе координат г, р, z) имеют вид:
5а гг +1 дст гф ^ до„
дr
да
r дф
дr
да
фг + і дафф |
дz
да
а - а
rz rr фф
|-----------—+ pFr =p ar;
дф
дг+2 афr +p ^=p aф;
(5)
і да zrn
zr I 1 Zip ^ zr
дr
r
дф
|
|
дz
- + p Fz =p az ,
где о у - компоненты тензора напряжений,
F, F., Fz )
вектор массовых сил.
Задачу будем решать полуобратным методом. Относительно поля скоростей сделаем следующие предложения: Уг = 0; Уф = 0; У2 = У2 (г). При данных
предположениях тензор скоростей деформации имеет следующий вид:
0 ; 8 гф 0 ; 8 фф 0 ; 8 фz 0 ; 8 zz 0 ; 8 rz о Vz (r) .
(б)
Из определяющих соотношений следует запись тензора напряжений:
-
а; =
p 0 аг
o - p o
а rz 0 - p
л
rz
д p
С учетом (6) имеем следующие уравнения движения = 0
и
dа„
dr
а
dp
dp
Ap
1 Ap
+ -р = 0. Предположим также, что -j- = const = ^-^аzr = -^^рr
dz
dz
l
2l
„ Ц і Ар
интенсивность тензора напряжений т = Ьу&у = 2 /”г, интенсивность девиатора
тензора скоростей деформации Н = ^2£у ^ = |У.,'(г).
Тогда из определяющих соотношений для интенсивностей:
Ap
Т
r = тп + K
}’-(г )'
8 0
1 j ap j
2 l
Г -Ir
Kn
(7)
где n = —. m
Поскольку У2 (г) убывает при г а, то У2 (г) < 0 и перед скобками должен стоять минус. Далее считаем, что в0 = 1-с-1, поэтому в формулах эту величину опускаем.
а
m
n
8
V; (в )= 0.
Из условия прилипания рассматриваемой среды к стенкам капилляра имеем
„я .1 2 / г -Т0
С 1 2/ Г1 Ы А
ёг= С1 “к7ЕЖ/ГГ) 2“Гг-Т0
п +1
У
V(а) = 0 ^ С, =4
2/
1 \^Р\
2 Т
Кп |Ар|(п +1)
В итоге получим следующую формулу:
^ \ \ \ \п +1 С
У г )=-! 2/ Г1М ' '
а - тп
2/
кп |Ар|(п +1)
2/
а -тг
1 1АР|
2 /
п +1
Г - Тг
(8)
(9)
(10)
Находим радиус области зон застоя из условия т(г) ) = Т0. Тогда
1 1АР|
2 /
2т0/
|Ар|
что совпадает с радиусом зон застоя для линейных
вязкопластических сред [6].
Вычислим расход среды Q в произвольном сечении капилляра:
а Г0
Q =| 2жгУ ; (г )ёг + 12жгУ ; (г0 )ёг ,
(11)
Г0
0
где У ; (г0 ) - скорость движения среды в зоне застоя.
После подстановки скоростей в (11) и интегрирования получим:
Л 1а„1 У+1
Q =
2п/
1 АР
2 1а -Т0
V У
Кп |Ар| (п + 1)(п + 2)(п + 3)
При п = 1 получим известную формулу Букингэма [8]:
Кп + 1)(п + 2)а2 + 2г<2 + 2аг0 (п +1)| . (12)
Q =
п/
1 Ар
к| Ар | з ■ 41
после преобразования которой получим:
^2~Та-Т0) (ба2 + 2г02 + 4аг0 ),
= жё4 Ар I. - 4 г0 +1Г г0 ^ 16К ■ 8 / | 3 а 3 !а
(13)
(14)
где ё = 2а - диаметр капилляра.
4
Пример расчета параметров течения
Рассмотрим решение методом эффективной вязкости задачи о течении пластической среды через сечение капилляра при стационарном перепаде давления.
В соответствии с вышеуказанным порядком определения эквивалентной вязкости находим вначале решение соответствующей задачи для вязкой жидкости с динамическим коэффициентом вязкости Г1е. Воспользовавшись, например, решением, приведенным в [8], запишем:
У; =
1 ёр 1 ёУ; 1 ёр
4ле ё; ёг 2ле ё;
(15)
г
Здесь
ёр
ёг
- постоянный перепад давления, Я - радиус сосуда, г - текущая
радиальная координата (ось г направлена вдоль центральной оси). Находим интенсивность скоростей деформации
Н =
ёУ 1 ёр
ёг 2Ле ёг
г.
Среднее значение интенсивности скоростей деформации по радиусу области вязкопластического течения вычисляется по формуле:
Я
< Н > = 1
Я - г
| Нёг
1 ёр
4Ле ёг
(Я + го).
(16)
Подставляя в (4) среднее значение Н из (16), получим уравнение для определения эквивалентной вязкости. Уравнение имеет следующий вид:
Ле =
|Д ёр
(4Ле )т-‘ ёг
т-1
(Я + г0 )
т-1
+
4тоЛе
ёр
ёг
(Я * г0 )
(17)
Отсюда следует, что в области вязкопластического течения г0 < г < Я
выражение для скорости У2 (г) имеет вид
V (г )=|
(Я + г>)- 4г>
1/т
■(я2 - г2
).
(18)
а =
п
(4м)п
>2 + г2
)|,
(19)
^ ) (Я + г0)
В интервале г0 < г < 0 ,5Я приближенная оценка, вычисленная методом эффективной вязкости, отличается от точной не более чем на 30% для любых нелинейно-вязкопластических материалов. Для объемного расхода имеем формулу
'Ц}п (Я - г0 п{ 2(Я!
что совпадает с приближенной формулой, приведенной в работе [6].
Чтобы сравнить, насколько величина расхода, определяемая по приближенной формуле отличается от расхода, определяемого точной формулой, оценим относительную погрешность:
|Л0| Г 2 '1 (1 * к 2 )п *1 Хп * 2)(п * з)
2 V--------ч----2----1----ч---1. (20)
^ (п * 1 )(п * 2) * 2к * 2к(п * 1)
При п = 1 величина расхода (19) совпадает с формулой, приведенной в [6], и максимальная погрешность достигает 12,1% при к = 0,41, что дает достаточно надежную расчетную оценку. Однако с увеличением п погрешность растет и при п = 4 достигает 60%, т.е. для существенно нелинейных вязкопластических сред приближенная зависимость (19) расхода оказывается неприменимой.
Г
0
Выводы
Метод «эффективной вязкости» может быть применен в задаче определения момента в условиях течения Куэтта. В особенности, если вязкая составляющая существенно больше пластической, в первом приближении вязкопластическое сопротивление может быть описано в терминах «эффективной вязкости». При
сравнимых величинах ^Hm /sо и То применение метода не гарантирует нужной точности.
Для течения Пуазейля метод “эффективной вязкости” может быть применен при n, близком к единице, и при достаточно большом перепаде давления. Для существенно нелинейных материалов, когда значение п близко к 0,25, метод неприменим.
При идентификации нелинейного вязкопластического материала обнаружено, что при малых То (То < 0,3 К) процесс нахождения То неустойчив. При То > 0,5 К погрешность не превышает 10%, остальные константы (n, K) вычисляются со значительно более высокой точностью.
Разработанная модель учитывает реологические свойства крови и может быть использована для уточненного описания гемодинамических процессов.
Литература
1. Регирер С.А. Биофизические основы движения крови по сосудам (Глава 14). В книге: Руководство по кардиологии. Т. 1. Структура и функция сердечно-сосудистой системы в норме и при патологии. Под ред. Чазова Е.И. М.: Медицина, 1982.
2. Nyashin Y.I., Nyashin M. Y., Shabrykina N.S. Models of microcirculation and extravascular fluid exchange // Russian Journal of Biomechanics. 2002. V. 6 №: 1. P. 62-77.
3. Akulich Yu.V., Podgaets R.M. Simulation of the bone tissue adaptation process under the proportional
change of the load // Russian Journal of Biomechanics. 2002. V. 6, № 1. Р. 25-33.
4. Руководство по ортопедической стоматологии. Под ред. В.Н. Копейкина. М.: Медицина, 1993.
5. Ильюшин А.А. К вопросу о вязкопластическом течении материала. Труды конф. по пластичности. Деформации. М.: Изд-во АН СССР: 5-18, 1936.
6. Гноевой А.В., Климов Д.М., Чесноков В.М. Об одном методе исследования пространственных течений вязкопластичных сред // Известия РАН, Механика твердого тела. 1993. № 3. С. 150-158.
7. Гавриленко С.Л., Шилько С.В., Васин Р.А. Идентификация определяющих соотношений вязкопластического материала в условиях течения Куэтта // Прикладная механика и техническая физика. 2002. Т. 43. № 3. С. 117-124.
8. Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1973.
9. Огибалов М.М., Мирзоджанзаде А.Х. Нестационарные движения вязкопластичных сред. М.: Изд-во МГУ, 1997.
10. Ишлинский А.Ю. Уравнения деформирования не вполне упругих и вязкопластичных тел // Известия АН СССР. Отделение технических наук. 1945. № 1-2. С. 34-45.
11. Shulman Z.P., Mansurov V.A., Makhaniok A.A. Simulation of non-steady flow of linear viscoplastic media // Inz. chem. i process. 2000. V. 21. № 1. P. 219-224.
A METHOD FOR DETERMINING FLOW AND RHEOLOGICAL CONSTANTS OF VISCOPLASTIC BIOMATERIALS. PART 1
S.L. Gavrilenko (Gomel, Belarus), R.A. Vasin (Moscow, Russia),
S.V. Shilko (Gomel, Belarus)
A model of biotissues is proposed to display viscoplastic properties in the presence of strong nonlinearity of velocity sensitivity and corresponding basic relations are derived. By way of example, determination of rheological constants of blood are considered under the Couette-Poiseuille flow conditions. An approximate calculation method of flow characteristics is described based on the hypothesis of equivalent viscosity. The algorithm stability is estimated in view of the initial data errors.
Keywords: blood circulation, viscoplasticity, non-Newtonian liquid, equivalent viscosity, Couette-Poiseuille’s flow.
Получено 13 сентября 2002 года