Метод оптимизации в гемодинамике Текст научной статьи по специальности «Математика»

УДК 533.70; 542.12; 532.5:532.135

Вестник СПбГУ. Сер. 1. 2012. Вып. 2

МЕТОД ОПТИМИЗАЦИИ В ГЕМОДИНАМИКЕ*

В. А. Цибаров1, Д. А. Юдинцева2

1. С.-Петербургский государственный университет, д-р физ.-мат. наук, профессор, [email protected]

2. С.-Петербургский государственный университет, аспирант, [email protected]

В проблеме гемодинамики живого организма можно выделить две основные задачи. Первая — это глобальная задача описания функционирования кровеносной системы. Для ее решения применяются [1] элементы теории управления. Это обусловлено тем, что в большой системе могут появляться дополнительные крупномасштабные процессы и соответствующие им связи, которые допускают огрубление описания, учитывающее только «следы более мелкомасштабных процессов».

Вторая задача заключается в описании гемодинамики участка кровеносного сосуда. Именно эта задача рассматривается в настоящей работе. Вторая задача, в свою очередь, также разделяется на две подзадачи: проблему реологии крови, являющейся гетерогенной средой (гидровзвесью), и вытекающую из нее вторую — проблему переноса лабораторных результатов на реальный сосуд (in vitro ^ in vivo).

Решение задачи о реологии крови осуществляется на основе стохастического подхода [2-5] в условиях слабо неравновесного течения. Замыкающие соотношения взяты из [3, 5, 6], но они дополнены возможностью учета напряжений сдвига. В настоящей работе при решении задачи о течении вязко-пластичной гетерогенной жидкости в сравнительно мелких сосудах, как и в [3, 5, 6], производится отображение участка кровеносного сосуда с квазиэллиптическим нормальным сечением, обладающего переменными и подвижными стенками, на круговой цилиндр единичного радиуса и единичной длины. Как и в [3, 6], течение крови в отображенной области описывается с помощью приближенных уравнений, имеющих минимальную интегральную квадратичную невязку внутри рассматриваемой области. Это позволяет «отщепить» задачу об отыскании гемодинамических профилей внутри отображенной области от задачи об отыскании сглаженных характеристик, зависящих от времени t и продольной координаты или от времени. В сочетании методов отображения, оптимизации и расщепления заключается главное отличие от традиционных приближенных методов решения задач о медленном течении вязкой жидкости. Еще одно отличие от традиционного подхода (см. [7-12]) заключается в удовлетворении решения системы уравнений гемодинамики (подобно [3, 6]) граничным условиям смешанного и скачкообразного типа, обычно называемым условиями скольжения, вместо условий прилипания, что позволяет учесть проскальзывание форменных элементов крови вдоль стенки здорового сосуда. Существен, как и в [5, 13], учет многочастичных столкновений форменных элементов, приводящий к более реальным значениям коэффициентов переноса. Кроме того, в настоящей работе формулируется условие возможности перехода in vitro ^ in vivo для гетерогенной среды в динамике.

* Работа выполнена при финансовой поддержке Министерства образования и науки РФ (проект 13.G25.31.0076). Доклад на Международной конференции по механике «Шестые Поляховские чтения» 31 января — 3 февраля 2012 г., Санкт-Петербург, Россия.

© В. А. Цибаров, Д. А. Юдинцева, 2012

1. Введение. В практических задачах часто применяется ньютоновская модель крови [1, 7, 8]. Однако уже в [9] отмечается недостаток такой модели, хотя и в этой монографии при обработке экспериментальных данных по эффективной вязкости крови практически не учитывается ее неньютоновский характер. Среди моделей, учитывающих неньютоновские свойства крови, получила распространение модель Кессона [10]. Применяется и модель вязко-пластичной среды Шведова—Бингама [11, 12]. В работах [6, 14] предлагается и используется модификация степенной вязкой модели Оствальда. Нашла применение [9, 11, 12] псевдоньютоновская модель крови с эффективной вязкостью, зависящей от эффективной (постоянной) скорости сдвига.

В настоящей публикации применяется модификация модели Шведова—Бин-гама—Оствальда, полученная кинетическими методами по результатам работ [3, 5, 15]. Используемая здесь кинетическая модель работы [5] удачно сочетает баланс простоты и строгости. Она позволяет воспользоваться известными результатами теории концентрированных газовзвесей и плотных газов [3], легко модифицируя их. Кроме того, в ней на уровне релаксационной модели учтен эффект многочастичных взаимодействий включений, что при больших концентрациях форменных элементов крови дает более реальные значения коэффициентов вязкости диспергированной фазы. Важным моментом является зависимость коэффициентов переноса диспергированной фазы, получаемых из стохастических законов сохранения, от внешних (гидродинамических) полей и от частоты соударения включений. Корреляция форменных элементов при соударении учитывается на уровне парной локально равновесной корреляционной функции х(с), зависящей от объемной доли (с) всех форменных элементов крови, и берется из [3].

Влияние объемной доли примеси на сдвиговую вязкость плазмы крови, рассматриваемой как ньютоновская среда, можно вычислить по эмпирической формуле из [16]:

и ( 2,5 с ] 2

^ = { 2 (1 — 1,35 с)} ' С<°'5' и

где и и и — сдвиговая вязкость гетерогенной и несущей среды, а с — объемная доля включений в гетерогенной среде. В этом случае удельная энергия Ер = 0р®р/4 хаотического движения включений и их наиболее вероятная скорость ар однозначно вычисляются по их объемной доле, полному давлению р и числу возбужденных степеней свободы ]р, поскольку полное давление и давление фаз (рр и pf) связаны соотношениями

1 2, / 1./1 (1 + 2кР \ 2££а/3 Р=Р1+Рр, РР = -ррафр = ср, фр = 1+4:сх ——--ар), ар = --.

2 р \1 + кр ; трар

Здесь кр € [0; 2/3] —безразмерный момент инерции включений [3], тр —средняя масса включения, еав = (0, 2595 ± 0, 0085) • 10-11 дин • см — глубина потенциальной ямы, полученная по экспериментальным данным из [12], е = 1 — с, а корреляционная функция х = {1 — 0,160 с2 — 0, 366 с3 + 16, 390 с4 + 0(с5)}/{е (1 — 1, 5 с)}, плотность крови р = pf + рр, pf = е^р = с7р, Yf и 7р —истинные (табличные) плотности фаз. Выражение для корреляционной функции применимо не только для монодисперсных гидро-и газовзвесей, но и для плотных простых газов.

Реальная кровь полидисперсна. В кинетической (стохастической) модели крови этот эффект, как ив [5], учитывается с помощью двух континуальных индивидуальных случайных переменных — объема и массы включения. В уравнениях гемодинамики он учитывается с помощью коэффициента хУ = 0, 96704 из [5].

Пусть i? и L — радиус и длина участка кровеносного сосуда, v — средняя скорость течения крови по сосуду, г/ср—частота пульса, a Sh = Ri/cp/v — число Струхаля. В соответствии с данными из [7, 9] у здоровых людей Sh ~ 10~3 + 10~2, кроме полых вен, у которых Sh ~ 10_1, а kz = R/L ~ 10~3 + 10~2. У части артериол может быть kz ~ 10_1. Приведенные оценки позволяют в отображенной области выделить ведущий оператор, соответствующий дифференцированию по радиальной координате. Это важно для реализации метода оптимизации.

2. Метод оптимизации. Как известно [7], реальный сосуд близок к квазиэллиптическому с полуосями a(z, t) и b(z, t), где z — продольная координата. Отобразим внутреннюю область участка сосуда на внутреннюю поверхность кругового цилиндра единичного радиус и единичной длины. Для этого осуществим переход от декартовой системе координат {x,y,z) к безразмерным декартовым координатам (ж, у, С), к размерным (r, в, z) и безразмерным £, ф, Z цилиндрическим координатам по формулам

ж = a(C,í) ж, y = b(C,t)y, z = L(t)(,

X = г cos в , у = г sin в , Z = Z , Г = \/X2 + у2 ,

ж = £ cos ф, у = £ sin ф, с = С • При такой замене справедливы следующие важные соотношения:

г = £Щ, Щ = у*'a2 cos2 ф + Ъ2 sin2 ф, в = arctg{(6/a) tg </>}, ф = arctg (у/ж), Vr = r, ve = гв, Vz = z = LZ + Lz, v^ = £, Уф = £ф, vc = Z, vx = a^ cos ф + ау^ cos ф — аУф sin ф, vy = b£ sin ф + bví sin ф + Ьуф cos ф, д Щ д д cos ф д sin ф д д sin ф д cos ф д

дв ab дф1 дх а д£ £а дф1 ду b д£ £b дф1

а также равенство

д ( dv\ д ( dv\ f cos2 ф sin2 ф\ д ( dv\

/sin2 ф cos2 ф\ 1 ( dv\

Точка над величиной означает дифференцирование по времени вдоль траектории.

Из последнего равенства следует отсутствие осесимметричного решения рассматриваемой задачи, кроме случаев ньютоновской жидкости или кругового цилиндра, на чем основаны лабораторные исследования. Но задачу можно оптимизировать, если такие операторы усреднить по всем углам ф G [0;2эт], что соответствует минимуму интегральной квадратичной невязки приближенного уравнения (или выражения) по указанной области [3]. Если ¡ = ¡(£,Z,t), то последнее выражение после усреднения примет вид

т, ^ 1 Í2 т, s и 11 д Г dv\ п2 2a2b2 В области (£, ф) боковая поверхность сосуда удовлетворяет уравнению £ — 1 = 0.

b = b(z,t).

Компоненты vr, vg и vz вектора скорости крови представимы [5, 6] в форме

b2 - a2 • /о^ £

Щ + Уф-—— sin (2</>) + —— (adcos2 ф + ЪЪsin2 ф)

2 V У

ab b a — a b ■

vg = — +£———sin (2</>), = wcL(t) + CL.

2

Для «отщепления» уравнения для продольной скорости течения крови член с градиентом давления усредняется по поперечному сечению кровеносного сосуда. При такой операции оптимизированное уравнение обладает минимумом интегральной квадратичной невязки по локальному сечению сосуда [3].

В практической медицине кровь, как правило, считают [7-9, 11] несжимаемой жидкостью, что неверно, вообще говоря, в силу динамических процессов в гетерогенной среде. Такой подход оправдан, если сжимаемость и нестационарный характер течения крови компенсируются тангенциальной компонентой ее вектора скорости. Указанное сужение класса рассматриваемых гемодинамических задач соответствует условию

= Bv = P"t -^°UtPUtfoutput R=Q {kz Sh) ; (2)

где F — площадь поперечного сечения, а V — объем участка сосуда. Словами «input» и «output» обозначены значения параметров на входе в участок сосуда и выходе из него. Условие (2) является следствием вида оператора дивергенции в переменных

(£,Ф,С):

1 д 1 дvф dvz V

Условие (2) позволяет решать оптимизированную в пространстве (£, ф, Z) задачу гемодинамики с плотностью цельной крови р, зависящей от времени как от параметра.

2.1. Оптимизированная система. В осесимметричном приближении в переменных (£, ф, Z) с учетом малости числа Струхаля и большого удлинения сосуда для степенной жидкости при наличии напряжений сдвига то = const (жидкости Шведова—Бингама—Оствальда) с точностью до членов порядка o (kz, Sh) она имеет вид

1 д . dvz , ,

Z дс ^ dZ

| = °. = KnSh,Kn*.), (4)

dp Кпд

Т = const, — = С\ = const, — = С*2 = 1 — С\ = const, (6)

рр

число Кнудсена Kn = Poov/ (pmR), Poo —характерное стандартное (табличное) значение вязкости крови, а pm и T — среднее ее давление (свое для каждого типа сосудов) и температура. Скорость v отнесена к v. Разность w скоростей фаз Vf и (плазмы крови и ее форменных элементов) вычисляется по формуле [5]

(7P-7f)(3g + 3)4 f tv = Vf — =-:---< р g к - Vp + V ■ П >, 7

Р 18pfipaRelpKcAf(3a + 2)plFy * У

r

где k — орт, направленный против сил тяжести, g — ускорение сил тяжести, П — тензор вязко-пластичных напряжений в цельной крови (в соответствии с моделью Шведова—Бингама—Оствальда). Через p в (4)—(7) обозначено превышение давления над гидростатическим, а — отношение вязкости внутри включения к сдвиговой вязкости ¡f несущей среды, dp — эффективный диаметр включения, вычисляемый по их среднему объему, Af = 1 +0,158Rep/3. Сдвиговая вязкость ¡f вычисляется по формуле (1), в которой ¡¡f заменяется на вязкость воды, а с — на объемную долю включений в плазме крови. Число Rep вычисляется по 7f, ¡¡f, dp и |w|. Параметр cpa = 1,187/lg (15, 385/с) при Rep < 0, 2 и cpa = /,Г0'9 при 0, 2 < Rep < 2000 (см. [16]). Коэффициент сферичности /с G (0; 1]. Показатель I = 0 при Rep < 0, 2 и I = \/l/fc — 1 при Rep е [0, 2; 2000], а

K = (1 - Л5)/До , До = (1 - Лс)4 (1 + 1,75 Ас + Л2) , Ас = с1/3.

Термическое уравнение состояния фазы включений имеет вид [5]

d ln pp

d ln Pp V Jp

2д

(8)

где 3Р —число степеней свободы у включения. Аппроксимация (8) политропой дает

, (9)

РрЪ PPab \PlabJ V JpA Cc / c=clab,a =0

поскольку политропа является интегралом движения. Индексом «lab» помечены данные лабораторного анализа. При оптимизации правой части уравнения (8) имеем

(ср. [5])

2 \ f ^ / ч Qab ^ V'pl^max)

??р = ( 1 + —j <¡ 1 +1п^р(Стах) +

1

'-max - Clab ^p(ciab) cmax clab

ln фр dc>.

Объемная доля С}аь определяется по гематокриту крови, стах < 2/3.

Реология крови различна внутри ядра, £ € [0;£о], и вне его (при £ € [£о;1]). Внутри ядра = то, скорость = Уо (неизвестна), = 0. По данным [12] напряжения сдвига то = 2,0 дин/см2 (у здоровых людей) и то = 2, 2 дин/см2 (у больных раком толстой кишки). Вне ядра (с учетом напряжений сдвига и антисимметричных напряжений)

— = — + Кп

V 8vq

/ Jabx 0,970

n = 0,185 .

д£ \PlabJ

n ^

dvz

pm pm

Время релаксации tp и вычисляются по формулам [6, 14]

\2

(10)

tp =

+

V^ií6(i + M Л+ Л 1 + 5^ 48V2 X I (6 + Шр) V V 1 + ^р

dp/Xv

+

768 1 + 2кр 2 2 257Г 1 + кр° Х

сарфр{1 + ер +е')

= 2, 5 х ■

¡f

m = 1 +

5а/7Г 1 f 6(1 + fcp)2

48л/2 X I (6 + 13Лр)

1 + 4 cx 0,4

1 + 5fcp 1 + kp

+

1728 1 + 2kp 22 + 25тг 1 + kp° X

Rep/ (XvX)

__rT _ PiaPdP

/ / -I _ _ / \ ? JA-fp -

cap^p(l + ep + e ) e pf

c

2

— a

p

c

2

a

p

При с < 0, 5 для приближенного определения рж/р^ можно воспользоваться экспериментальной формулой (1). Влияние внешнего (гидродинамического) поля на коэффициенты переноса описывается [3, 5, 6] безразмерными параметрами

5тг (1 + кр)2 А^/Йёр р{ Г „Л ,рЛ (л , п 1сотз2/ЗА За + 2

= Ьтги + EpJ fcf/rv-ep * - / _ л , Re2/3)

8а/2 (6 + 13/гр) с\ 7р V ^ V р ; а + 1

Параметр Ьр1 = Зя^®/ (^р®^, где к^Р и Яр —коэффициенты объемной и сдвиговой вязкости, вычисляемые по [3], характеризует антисимметричные напряжения в плазме.

2.2. Критерии отбора решений. На неизвестной границе £ = £0 должны выполняться условия = Уо и = то, а при £ =1 — условия

_ лД 4 - Лр(1 - 5kt) Qp ра

лр — -

Р 2 2(1 - kt)kp g Pp R

v

dv,

R d£

í=i

o

4B dap чп i í. vi t^ , ou п í.2

p_1 3(l + fcp)(l-fc2) + 2fcp(l-fct2)

?=1' 10(1 + M

a2(l) = ~ЛР~

(3 kp - 5)(1 - k?) , Я | Ap^ MNScP

+ 1 15 (1 + fcp) 10 ^ ~ J Л^ ' РГр- A- •

Число Прандтля Prp вычисляется по слабо неравновесному решению из [3], kt € [—1; —|—1] и kn € [— 1; 0] —коэффициенты восстановления касательной и нормальной составляющей импульса включения при не вполне упругом взаимодействии с границей, [p] — скачок величины p. К граничным условиям следует добавить условие задания средней скорости v и интегральное условие перехода in vitro in vivo:

2 i vz £d£ = 1, í í Uc(1,Z,t) dtdZ = 0, (11)

Jo J0 J0

где Uq(1, Z,t) —смещение вдоль сосуда. При записи (11) учтено, что dxdy = ab£d£d4i.

2.3. Решение оптимизированной системы. Решение системы (3)-(5), удовлетворяющее граничным условиям и первому из условий (11), внутри ядра имеет вид

С v С v R 2go(m + go) n + 2

t0 = £оКптс? 1 = £0Кптш = kz — — , В o =--—г-, то =——■,

2 dZ m(m +1 n +1

а решение вне ядра — вид

e-eo o j i-V°c(Í) f

m 5Cll + (m + 2)Ap/(l-eo) + Bo V U '

n

a

„ Со fm е0 л - (1-еоГ r tp41/(1+п)

m + 1 \ £ £2 / ш + 2 tp у р,

1 + + 7eff =

(m + 2)(l -»(°(1))й/(шД)

1 + (ш + 2) Лр/ (1 - £o)+Bo

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

Meff _ (m + 2)(í-v¡(í))vtp/R Ро

1 + (m + 2) Лр/ (1 - £о)+Вс

2(1 - to)1 П 1 + to + £on(m - 1)

n(m — 1) + 2 n(m — 1) + 1

(12)

При to = 0 формула (12) применялась в работах [6, 14].

Для сглаженного по углам ф G [0; 2п] значения Уф имеем Уф = £ Bv, уф (1) = Bv.

2.4• Нахождение сглаж.енных величин. Первый этап метода оптимизации привел к решению, содержащему неизвестный параметр у0(1) = vo(Z,t). Проинтегрируем уравнение количества движения цельной крови по поперечному сечению сосуда. В предположении непроницаемости его стенок получим совпадение нормальных проекций тензора напряжений в крови и в стенке сосуда. Предположим ее упругой. Эти предположения приводят к акустическому уравнению

где E — модуль Юнга, v — коэффициент Пуассона, U(1,Z,t) = U/L, Uq — продольное смещение жидкости, связанное с подвижностью реальных границ. Границы отображенной области неподвижны. Поэтому для уравнения (13) выдвигаются условия

U(1,С,0) = 0, U(1,0, t) = U(1,1,t) = 0 . (14)

К условиям (14) добавляется второе из условий (11).

3. Вычисленные параметры. В предположении слабого изменения со временем плотности крови и L из решения уравнения (13) при указанных выше условиях отбора получены главные члены выражений для U и у0(1) = dU/dt:

U(1,Z,t) = ^^ sin (2nnat) sin (2nnZ) , v°(1,C,t) = 2nna ^^ cos (2nnat) sin (2nnZ) .

В качестве примера течения вязко-пластичной жидкости Шведова—Бингама— Оствальда рассмотрено течение в артериоле. Для здоровой крови человека реологический коэффициент n = 0,088. В артериоле среднее давление pm = 0, 043 атм. За V примем допустимое значение 0,6 см/с. Пусть R = 0,03 см, 7ед = 1000 с-1,

= 3, 5 сП, L = 2, 6 см, kp =0,4 (как у шара). В этих условиях вычисления дают c =0, 57, vo = 2,14, to = 0, 71, tp = 0,038 с, наиболее вероятная скорость форменных элементов ap = 0,044 см/с. При kp = 2/3 (полые сферы) ap = 0, 045 см/с.

Метод оптимизации описывает поступательно-пульсирующее течение крови. Радиальная и тангенциальная скорости обусловлены подвижностью стенок сосуда и

х

переменностью площади его сечения по длине канала. В пределе очень протяженных сосудов (кг ^ 0) радиальная и тангенциальная скорости крови исчезают.

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

Литература

1. Лищук В. А. Математическая теория кровообращения. М.: Медицина, 1991. 256 с.

2. Цибаров В. А. Кинетика и гемодинамика крови. I // Вестн. С.-Петерб. ун-та. Сер. 1. 1997. Вып. 4 (№22). С. 101-107.

3. Цибаров В. А. Кинетический метод в теории газовзвесей. СПб.: Изд-во С.-Петерб. ун-та, 1997. 192 с.

4. Цибаров В. А., Цибарова Е. В. Сохастическая модель сложных сред // Аэродинамика / под ред. Р. Н. Мирошина. СПб.: НИИ Химии СПбГУ, 2002. С. 90-111.

5. Цибаров В. А. Стохастический метод в гемодинамике сосуда // Вестн. С.-Петерб. ун-та. Сер. 1. 2009. Вып. 1. С. 129-138.

6. Фомина О.Н., Цибаров В. А. Гемодинамика отрезка сосуда // Вторые Поляховские чтения: Избранные труды. СПб.: НИИ Химии СПбГУ, 2000. С. 179-189.

7. Бабский Е. Б., Зубков А. А., Косицкий Г. И. и др. Физиология человека. М.: Медицина, 1966. 656 с.

8. Парашин В. Б, Иткин Г. П. Биомеханика кровообращения. М.: Изд-во МГТУ им. Н. Э. Баумана, 2005. 224 с.

9. Левтов В. А., Регирер С. А., Шадрина Н. Х. Реология крови. М.: Медицина, 1982. 272 с.

10. Лосев Е. С., Нетребко Н. В. Моделирование реологичнского поведения крови в нестационарных сдвиговых течениях // Изв. АН. МЖГ. 1995. №6. С. 25-30.

11. Гидродинамика кровообращения / под ред С. А. Регирера. М.: Мир, 1971. 269 с.

12. Ганцев К. Ш. Гемореология у больных раком толстой кишки. — http://im.mtometeo.ru /УевШ1к/6/аг^с1е-2004-04-02. 6 с.

13. Долидович Н. Ю., Цибаров В. А. Медленные течения плотного газа // Пятые Поляховские чтения: Избранные труды. СПб.: Изд-во С.-Петерб. ун-та, 2008. С. 57.

14. Нагорный С. С, Цибаров В. А., Цой С. В. Вязкость крови как неньютоновской среды // Четвертые Поляховские чтения: Избранные труды. СПб.: НИИХ СПбГУ, 2009. С. 230-235.

15. Цибаров В. А. Сохастические законы сохранения в теории неньютоновских сред // Аэродинамика / под ред. Р. Н. Мирошина. СПб.: НИИ Химии СПбГУ, 2000. С. 93-119.

16. Разумов И. М. Пневмо- и гидротранспорт в химической промышленности. М.: Химия, 1979. 248 с.

Статья поступила в редакцию 26 декабря 2011 г.