Построение, численное моделирование и анализ одномерной модели гемодинамики

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

Рубрика Программирование, компьютеры и кибернетика
Вид курсовая работа
Язык русский
Дата добавления 24.09.2012
Размер файла 3,9 M

Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже

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

38

Размещено на http://www.allbest.ru/

Содержание

  • Введение
  • 1. Основная часть
  • 1.1 Вывод модели
  • 1.2 Анализ математических свойств модели
  • 1.3 Граничные условия
  • 1.4 Реализация модели
  • 1.4.1 Сведение к обыкновенному дифференциальному уравнению
  • 1.4.2 Линеаризация граничных условий
  • 1.4.3 Метод ортогональной прогонки
  • 1.4.4 Теория возмущений (замечание о точности решения)
  • 1.4.5 Генерация матриц граничных условий для метода ортогональной прогонки для системы из 55 артерий
  • 1.4.6 Описание работы программы
  • 1.4.7 Описание структуры программы
  • 1.5 Тестовые расчёты
  • 1.5.1 Моделирование гемодинамики в норме
  • 1.5.2 Моделирование физической нагрузки и тахикардии
  • 1.5.3 Применение метода итераций для оптимизации тестовых расчётов
  • Заключение
  • Список литературы
  • Приложение

Введение

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

Артериальное давление (АД) является важным параметром, характеризующим работу кровеносной системы. Нормальные численные значения АД лежат в пределах 120-135 мм рт ст для систолического и 75-85 мм рт. Ст. для диастолического в режиме бодрствования. Стойкое превышение этих значений характеризуется как артериальная гипертония (АГ), которая является одним из самых распространённых заболеваний ССС, и без должного и вовремя начатого лечения может привести к поражению различных органов человека. АГ занимает первое место по смертности в большинстве развитых стран мира, поэтому разработка мер и средств по ее предупреждению, лечению и диагностике является одной из важных задач современной науки.

Моделирование движения крови по организму человека является важной задачей. Созданием моделей процессов, происходящих в ССС, люди начали заниматься достаточно давно. Началом формального описания физиологии кровообращения можно считать работу Уильяма Гарвея "Анатомическое исследование о движении сердца и крови у животных (Exercitatio anatomica de motu cordis et sanguinis in animalibus), выпущенной в 1628 году. Он заключил, что кровь в организме человека движется по кругу. В этой книге Гарвей точно описал работу сердца и выделил малый и большой круги кровообращения. [6] Наряду с этим, Гарвей доказал, что сердце ритмически бьется до тех пор, пока в организме теплится жизнь, причем после каждого сокращения сердца наступает короткий перерыв в его работе, во время которого этот важный орган отдыхает (таким образом он описал систолу и диастолу).

Математические модели стали появляться гораздо позже. Одномерное моделирование артериальной системы человека было предложено Эйлером в 1775 году, который вывел частные дифференциальные уравнения, выражающие закон сохранения массы и закон сохранения импульса для идеальной (невязкой) жидкости. Для того, чтобы замкнуть систему уравнений, он предложил два возможных, но экспериментально нереалистичных, уравнения состояния, описывающих поведение эластичной стенки в зависимости от изменений давления в сосуде. Вероятно, он не определил природу потока, как волновую, и не мог найти решения своих уравнений, испытывая, как он сам сказал, "непреодолимые трудности". Волновая природа артериального течения была впервые научно описана Юнгом, который вывел волновую скорость, используя аргументы, базирующиеся на интуиции и аналогии с ньютоновской теорией скорости звука в воздухе [3]. Уравнения, выведенные Эйлером в 1775 году - это система нелинейных дифференциальных уравнений в частных производных, аналогичная волновому уравнению для мелкой воды в гидродинамике или уравнению, описывающему одномерное невязкое течение в газовой динамике.

Следует отметить, что физиологические параметры артериальной системы человека таковы, что большинство характеристик потока сохраняются при линеаризации системы уравнений. Это, главным образом, было использовано в методе Вумерслея (Womersley, 1957), который линеаризовал двумерные уравнения потока в прямых, круглых эластичных трубках и получил волновое решение с помощью преобразования Фурье.

Развивая математический аппарат, предложенный Эйлером, Навье (1822) и Стокс (1845) вывели дифференциальные уравнения, описывающие движения вязкой несжимаемой жидкости (уравнения Навье-Стокса). В 1842 году физик и врач Пуазейль связал расход жидкости с давлением. В начале XX в Франк предложил идею, что система кровообращения аналогична электрической цепи, что положило начало новому способу представления ССС, используемому до настоящего времени.

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

Используя предположение, что плотность крови постоянна, физические процессы, происходящие в артериальной системе человека, можно описать трёхмерными нестационарными уравнениями Навье-Стокса, представляющими собой систему дифференциальных уравнений второго порядка в частных производных, совместно с уравнениями динамики эластичных оболочек сосудов, предложенными Пуазейлем [4]. Аналитическое решение этих уравнений найдено лишь для некоторых частных случаев. Применение такой сложной многомерной модели требует нахождения численных решений задач со свободной границей для уравнения Навье-Стокса в сложных областях. Это приводит к большим вычислительным затратам, поэтому многомерные модели на практике для глобального описания артериальной системы не применяются, а используются только для детального описания кровотока в интересующих исследователя локальных областях. Более того, часто нет необходимости в детальном представлении процесса гемодинамики для дальнейшего анализа медиками-специалистами, поэтому даже упрощённые модели могут обеспечить достаточным объемом информации при разумных вычислительных затратах. Таким образом, в настоящее время для глобального описания кровотока используют различные приближения исходной трёхмерной модели. Модели кровообращения, описывающие гемодинамику системы "в целом", являются важнейшим инструментом для исследования взаимодействия гемодинамических параметров, изучения и диагностики повреждений ССС, оценки последствий тех или иных методов вмешательства в её функционирование и т.д.

В настоящее время в моделировании гемодинамики используются несколько моделей, отличающихся степенью детальности описания системы. Наиболее простыми (и характеризующимися "грубым" приближением реальной трехмерной модели ССС) являются модели с сосредоточенными параметрами (или "нульмерные" (0D) модели). Главной их особенностью является сопоставление кровеносной системе электрических цепей. При таком подходе аналогами давления и объемного кровотока являются напряжение и электрический ток, а аналогами сосудов - сопротивления отдельных сегментов электрической цепи.0D модели являются упрощением одномерных (1D) моделей гемодинамики. Одномерные модели являются гидравлическим приближением исходной математической 3D модели, и получаются её усреднением по "поперечному" направлению. Третья группа упрощённых моделей - это гибридные модели. Они представляют собой стыковку, например, 3D и 0D моделей, 3D и 1D моделей и т.д.

В настоящее время наиболее интересной для анализа представляется 1D модель в силу предоставляемой ею возможности проследить за изменением параметров на всей протяженности сосуда по сравнению с 0D моделью, а также ввиду большей простоты моделирования при достаточной детальности описания системы по сравнению с 2D и 3D моделями. Для решения уравнений 1D модели могут быть применены различные вычислительные методы, характеризующиеся различной точностью, устойчивостью к погрешностям, скоростью сходимости итераций и простотой использования.

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

1. Основная часть

1.1 Вывод модели

Одномерная модель описывает течение крови в артериях и ее взаимодействие с подвижными стенками. Пусть t - время, а (x,y,z) - декартовы координаты. Будем считать, что идеализацией артерии является эластичная трубка в виде прямого кругового цилиндра. Рассмотрим кусок этой трубки длины L=const от z=0 до z=L (Рис.1.1).

Рис.1.1 Упрощённая геометрия кровеносного сосуда.

Через D (t) мы будем обозначать пространственную область, которая является указанной трубкой, заполненной кровью. Эта область меняется со временем под воздействием пульсирующей жидкости (крови). Осевое сечение {z=const} цилиндра D (t) будем обозначать через S=S (t,z). Описанная упрощенная геометрия артерии естественным образом приводит нас к использованию цилиндрических координат (r, и, z).

Основные уравнения выводятся при следующих упрощающих предположениях:

1) Осевая симметрия. Все величины не зависят от угловой координаты и. Как следствие, каждое осевое сечение S остается круговым в течении всего времени движения стенок сосуда. Радиус трубки R является функцией времени t и осевой координаты z.

2) Радиальные смещения. Стенка сосуда перемещается только вдоль радиального направления, т.е. в каждой точке поверхности трубки вектор смещения

з = з,

где з = R ? смещение относительно исходного радиуса .

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

4) Постоянство давления на каждом осевом сечении. Предполагается, что давление p постоянно на каждом сечении, т.е. оно зависит только от z и t.

5) Отсутствие массовых сил. Мы пренебрегаем массовыми силами. (Xотя, например, учет силы гравитации не создает больших трудностей и связан только лишь с добавкой к давлению слагаемого вида gh. Немного сложнее учесть изменение гравитационных сил, возникающего при вставании человека из горизонтального положения.)

6) Преобладание осевой компоненты скорости. Считается, что компоненты скорости, ортогональные оси z, являются пренебрежимо малыми по сравнению с осевой компонентой скорости . Предполагается, что

(1.1)

где - средняя скорость на каждом осевом сечении, а v: [0, 1] > R - профиль скорости. (Тот факт, что профиль скорости явно не зависит от времени и пространственных переменных, вообще говоря, противоречит экспериментальным наблюдениям, а также численным результатам для полной трехмерной модели. Однако, мы можем считать, что v является профилем скорости для усредненной конфигурации течения.)

Через A=A (t,z) мы будем обозначать площадь осевого сечения S (t,z). Она находится по формуле

. (1.2)

Средняя скорость

(1.3)

где

(1.4)

средний объёмный кровоток. Из (1.1) следует, что

(1.5)

Введём в рассмотрение коэффициент

(1.6)

Который называют коэффициентом Кориолиса. Используя неравенство Коши-Буняковского, из (1.1) мы получаем, что б ? 1. В общем случае этот коэффициент является величиной, зависящей от времени и пространственных переменных, но в рамках нашей модели б = const, что является следствием упрощающего предположения (1.1).

В качестве возможного профиля скорости можно взять параболический профиль что соответствует известному решению Пуазейля, характерному для стационарных течений в круговых трубках. В этом случае б = 4/3. Однако, известно, что для кровотока в артериях профиль скорости в среднем достаточно "плоский”. Как правило, он описывается функцией где обычно г = 9 (величина г = 2 соответствует параболическому профилю). Тогда б = (г + 2) / (г + 1) = 1.1 Более того, выбор б = 1, соответствующий "абсолютно плоскому” профилю скорости (полагая формально г = ?, имеем v (y) = 1 для и v (1) = 0), значительно упрощает последующий анализ модели. Заметим, что в большинстве работ (см., например, [7]) по умолчанию рассматриваются модели с б = 1.

Одномерная модель выводится с помощью интегрирования уравнений Навье-Стокса по произвольному осевому сечению S, т.е. усреднением этих уравнений по "поперечному" направлению.

Одномерная модель представляется системой из двух дифференциальных уравнений с частными производными:

(1.7)

где A, Q - неизвестные,

б - коэффициент Кориолиса,

с - плотность,

- коэффициент трения, который зависит от выбранного профиля скорости, т.е. от выбора функции v в (1.1).

Система, рассматривается при для интервала времени (0,T).

Для того чтобы замкнуть систему надо задать давление. Для задания закона давления необходимо учесть движение стенки кровеносного сосуда (в нашем случае артерии). Для описания упругих напряжений стенки сосуда, на практике чаще всего используется линейный закон (закон Гука)

,

где , E - модуль Юнга, - толщина стенки сосуда.

Тогда с учётом (1.2) получаем

где , .

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

(1.8)

И система (1.7), (1.8) является системой уравнений в частных производных первого порядка.

Рассматривается случай б = 1, (для параболического профиля)

н - постоянный коэффициент вязкости.

Случай б = 1 соответствует профилю скорости с г = ?. Для этого профиля скорости коэффициент трения . Заметим, что > ? при г > ?. Поэтому, уравнения (1.7), (1.8) с = 8рн являются некоторым огрублением общей одномерной модели гемодинамики. С другой стороны, модель (1.7), (1.8) с = 8рн часто вполне адекватно отражает реальные физические процессы, протекающие в артериальной системе человека. Поэтому и мы численные эксперименты проводим именно для этой модели, т.е. для уравнений (1.7) с б = 1 и = 8рн (и используем закон давления (1.8)).

1.2 Анализ математических свойств модели

Для простоты будем предполагать, что и (т.е. , и модуль Юнга E = const). С учетом (1.8) система уравнений (1.7) записывается в квазилинейном виде

(2.1)

для вектора неизвестных, где

, .

При естественном требовании A>0 система (2.1) является одномерной квазилинейной гиперболической системой уравнений. Матрица В имеет два вещественных и различных собственных значения

, (2.2)

где

(). Причём при стандартных физиологических условиях типичные значения скорости крови u = Q/A и механических характеристик стенок сосуда таковы, что , т.е. и . (2.3)

Таким образом, гиперболическая система (2.1) на интервале (0,L) имеет по одной уходящей характеристике на границах z=0 и z=L и поэтому она требует по одному граничному условию на каждой границе.

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

, (2.4)

где

.

Обсудим теперь кратко вопрос о приведении системы (2.1) к каноническому виду. Заметим, что гиперболическая система двух уравнений всегда может быть приведена к каноническому виду локально, т.е. в достаточно малой окрестности произвольной точки U. Однако, для частного случая б = 1 существуют также и глобальные инварианты Римана. Пусть (,) и (,) пары левых и правых собственных векторов матрицы B. Введем в рассмотрение матрицы

, , ,

где, не нарушая общности, LR = I. Тогда матрица B может быть представлена в виде

B = LЛR,

а система (2.1) эквивалентно переписана так:

.

Инвариантами Римана (для нашей системы) являются, как известно, величины и , удовлетворяющие соотношениям

,

Полагая , приводим систему (2.1) к диагональному виду

Нахождение инвариантов Римана и подробно обсуждается в [2]. Мы же здесь приведем только результат (для случая б = 1):

(2.5)

Заметим, что знание инвариантов Римана важно, например, для правильной постановки граничных условий при z = 0 и z = L. Так помимо того, что в силу (2.3) на каждой границе для системы (2.1) должно быть поставлено по одному граничному условию, необходимо, чтобы эти граничные условия разрешались для "уходящих” инвариантов Римана, т. е для инвариантов Римана, соответствующих уходящим характеристикам. Вопрос о граничных условиях обсуждается в следующем параграфе.

1.3 Граничные условия

Рассмотрим вначале случай одной артерии. Как уже отмечалось, для артерии длины L на левой (при z = 0) и на правой (при z = L) границах необходимо поставить по одному граничному условию. С математической точки зрения абсолютно неважно какими будут эти условия, а важно лишь то, чтобы они могли быть приведены к виду, когда "уходящие” инварианты Римана выражаются через "приходящие”. Так как в нашем случае в силу (2.3)"уходящим” инвариантом Римана на левой границе является , а на правой границе - , необходимо, чтобы граничные условия приводились к виду

, , (3.1)

где и - некоторые функции.

С другой стороны, граничные условия должны быть разумными с физической точки зрения.

На левой границе таким разумным граничным условием является задание давления:

, (3.2)

где - некоторый входной профиль давления, (в нормальном случае мы рассматривали кусок синусоиды) моделирующий давление на выходе из сердца (Рис.3.1.). С учётом (1.8), (2.5) можно проверить, что граничное условие (3.2) является хорошим и с математической точки зрения, поскольку оно может быть приведено к виду первого условия из (3.1).

Рис.3.1 График давление на выходе из сердца.

Что касается правой границы, то мы поставили там условие постоянства давлений,

=9.32*104Дин/см2=70 мм рт ст, (3.3)

которое также является хорошим с математической точки зрения.

С физической точки зрения, является вполне логичным задавать

в качестве правого граничного условия - условие постоянства площади осевого сечения A, на концах последних артерий

.

с физической точки зрения, является вполне логичным задавать

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

Кроме этого, рассматривался вариант задания, в качестве левого граничного условия, в виде зависимости потока от времени

,

но оказалось, что это некорректно ввиду условия 3.1.

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

Граничные условия в точках ветвления артерий.

Рис.3.2 Одномерная модель бифуркации артерии (метод разбиения области).

Артериальная система характеризуется наличием ветвлений. Течение крови в местах ветвления артерий является существенно трехмерным. Однако, используя метод разбиения области, можно описать его в рамках одномерной модели. На Рис.3.2 изображена схема одномерной бифуркации. Мы упрощаем реальную геометрическую структуру бифуркации, полагая, что ветвление происходит в одной точке, и пренебрегая эффектами, связанными с углами ветвления.

Для того, чтобы решить три системы уравнений в областях (основная ветвь), и соответственно, нам снова необходимо найти интерфейсные условия в точке z = Г (= Г±). Так как системы являются гиперболическими и имеют по одной уходящей характеристике на границе {z = Г, }, требуется поставить три интерфейсных условия. Для уравнения неразрывности (первое в системе (1.7)) получаем непрерывности объемного кровотока в точке z = Г (условие сохранения массы крови при переходе через точку бифуркации), т.е.

при z = Г. (3.4).

Что касается интерфейсного условия для уравнения движения (второго в системе (1.7)), то в литературе часто используется условие непрерывности давления p.

В работах Formaggia и Quarteroni [8] было показано, что для случая б = 1 интерфейсным условием, которое для задачи в подобластях , и гарантирует те же условия "численной" устойчивости, что и для исходной задачи в области D, является условие непрерывности на разрыве полного давления

. (3.5)

Таким образом, второе условие - непрерывность полного давления, в точке ветвления:

при z=Г. (3.6)

В итоге мы получаем три требуемых интерфейсных условия (3.4) и (3.6).

Аналогичным образом можно получить одномерную модель гемодинамики для всей артериальной системы человека, состоящей из 55 крупных артерий тела человека (Рис 3.3), список которых был использован в работе Lamponi (Табл.3.1) [1].

Таблица 3.1 Список и характеристики 55 крупных артерий тела человека.

#

Название артерии

Длина

(см)

Площадь сечения (см2)

Beta (в)

(кг/с2)

Коэффициент отражения (Rt)

1

Ascending Aorta

4.0

5.983

388

-

2

Aortic Arch I

2.0

5.147

348

-

3

Brachiocephalic

3.4

1.219

932

-

4

R. Subclavian I

3.4

0.562

1692

-

5

R. Carotid

17.7

0.432

2064

-

6

R. Vertebral

14.8

0.123

10360

0.302

7

R. Subclavian II

42.2

0.510

1864

-

8

R. radial

23.5

0.106

11464

0.273

9

R. Ulnar I

6.7

0.145

8984

-

10

R. Interosseous

7.9

0.031

51576

0.319

11

R. Ulnar II

17.1

0.133

9784

0.298

12

R. internal Carotid

17.6

0.121

10576

0.261

13

R. external Carotid

17.7

0.121

9868

0.26

14

Aortic Arch II

3.9

3.142

520

-

15

L. Carotid

20.8

0.430

2076

-

16

L. internal Carotid

17.6

0.121

10576

0.261

17

L. external Carotid

17.7

0.121

9868

0.264

18

Thoracic Aorta I

5.2

3.142

496

-

19

L. Subclavian I

3.4

0.562

1664

-

20

Vertebral

14.8

0.123

10360

0.302

21

L. Subclavian II

42.2

0.510

1864

-

22

L. Radial

23.5

0.106

11464

0.274

23

L. Ulnar I

6.7

0.145

8984

-

24

L. Interosseous

7.9

0.031

51576

0.319

25

L. Ulnar II

17.1

0.133

9784

0.298

26

Intercostals

8.0

0.196

3540

0.209

27

Thoracic Aorta II

10.4

3.017

468

-

28

Abdominal I

5.3

1.911

668

-

29

Celiac I

2.0

0.478

1900

-

30

Celiac II

1.0

0.126

7220

-

31

Hepatic

6.6

0.152

4568

0.308

32

Gastric

7.1

0.102

6268

0.307

33

Splenic

6.3

0.238

3224

0.31

34

Superior Mesenteric

5.9

0.430

2276

0.311

35

Abdominal II

1.0

1.247

908

-

36

L. Renal

3.2

0.332

2264

0.287

37

Abdominal III

1.0

1.021

1112

-

38

R. Renal

3.2

0.159

4724

0.287

39

Abdominal IV

10.6

0.697

1524

-

40

Inferior Mesenteric

5.0

0.080

7580

0.306

41

Abdominal V

1.0

0.578

1596

-

42

R.common Iliac

5.9

0.328

2596

-

43

L.common Iliac

5.8

0.328

2596

-

44

L. external iliac

14.4

0.252

5972

-

45

L. internal Iliac

5.0

0.181

12536

0.308

46

L. Femoral

44.3

0.139

10236

-

47

L. deep Femoral

12.6

0.126

10608

0.295

48

L. posterior Tibial

32.1

0.110

23232

0.241

49

L. anterior Tibial

34.3

0.060

36972

0.239

50

R. external Iliac

14.5

0.252

5972

-

51

R. internal Iliac

5.1

0.181

12536

0.308

52

R. Femoral

44.4

0.139

10236

-

53

R. deep Femoral

12.7

0.126

10608

0.296

54

L. posterior Tibial

32.3

0.110

23232

0.241

55

R. anterior Tibial

34.4

0.060

36972

0.239

Рис.3.3 Схема дерева из 55 крупных артерий человека.

Для этого на входе первой артерии (аорты) ставится граничное условие вида (3.2), моделирующее работу сердца, на выходе конечных артерий граничные условия вида (3.3), а в точках ветвления артерий мы требуем выполнения интерфейсных условий (3.4) и (3.6). При этом в каждой отдельной артерии течение крови описывается системой гемодинамики (1.7). Считая, что исходный радиус артерий не меняется с изменением их длины, мы ставим начальные данные для площади осевого сечения A для каждой артерии:

для (3.7)

где - площади сечения отдельных артерий взяты из таблицы параметров артериальной системы человека (Табл.3.1).

Что касается начального условия для объемного кровотока Q, то на практике в численных расчётах артериальной системы обычно используют условие .

Значения таких физических параметров как длина сосудов и коэффициент также были взяты из работы Lamponi [1].

1.4 Реализация модели

1.4.1 Сведение к обыкновенному дифференциальному уравнению

Рассмотрим уравнение (2.1), к которому свелась система гемодинамики

(4.1)

где вектор неизвестных, и при б = 1

, .

Требуется найти распределение U=U (t,z) вдоль оси z в любой момент времени t, при заданных начальных и граничных условиях.

Дискретизация по времени. Рассмотрим пространственно-временную систему координат {t,z}. В полуполосе построим по оси t полосы с шагом по времени

k=0,1,… Пусть Тогда

И вместо (4.1) получаем

(4.2)

Перепишем последнее выражение в виде

(4.3)

Таким образом, у нас получилось дифференциальное уравнение

(4.4)

где матрица G имеет вид

(4.5)

а вектор f есть

. (4.6)

Здесь в элементы G и f входят известные величины , которые берутся из предыдущего временного слоя.

1.4.2 Линеаризация граничных условий

Граничные условия (3.1), (3.3), (3.6) не являются линейными. Для линеаризации произведём следующее преобразование

(4.7)

Для линеаризации (3.6) также используются ис предыдущего шага по времени:

(4.8)

Обозначим

, (4.9)

, (4.10)

, (4.11)

тогда

. (4.12)

После этого, граничные условия (3.2), (3.3) для обыкновенного дифференциального уравнения (4.4), с учётом (4.9) - (4.12) формулируются в виде

, (4.13)

. (4.14)

Граничное условие (3.6) (в общем случае ) будет выглядеть так

. (4.15)

Решение. Для численного решения дифференциального уравнения (4.4) с граничными условиями (4.13), (4.14), (3.3), (4.15) используем метод ортогональной прогонки.

1.4.3 Метод ортогональной прогонки

Метод ортогональной прогонки предназначен для численного решения краевой задачи

, , , (4.16)

где А - матрица размера nЧn, f - вектор-функция размера n, L,R - прямоугольные матрицы размера kЧn и pЧn, их ранги соответственно равны k и p. Элементы матрицы A (x) и вектор-функции предполагаются непрерывными на отрезке. Метод основан на представлении решения задачи (4.16) через решения серий задач Коши. Метод ортогональной прогонки зарекомендовал себя на практике как эффективное и надежное средство численного решения краевой задачи. Он обладает повышенной устойчивостью к погрешностям округлений при реализации метода на компьютере.

Подробное описание алгоритма метода ортогональной прогонки представлено в работах [9] и [10].

1.4.4 Теория возмущений (замечание о точности решения)

Пусть краевая задача

, , , (4.17)

правильно поставлена. (Здесь A (x) - матрица размера n Ч n, L, R - прямоугольные матрицы размера k Ч n и p Ч n соответственно, причем k + p =n, вектор-функция f (x) - размера n). Это означает, что краевая задача (4.17) имеет единственное решение для любых , ш и любой непрерывной на отрезке [0, d] вектор-функции f (x). Тогда решение задачи представляется c помощью матриц Грина G (x, s), , (см. [11])

, (4.18)

где - матрица-функция размера n Ч k, решение задачи

матрица-функция размера n Ч p, решение задачи

где - нулевая матрица размера l Ч m.

Будем говорить, что задача (4.17) хорошо обусловлена, если для любых x и s из интервала верны оценки

(4.19)

Из (4.18) следует, что для решения хорошо обусловленной задачи (4.17) справедлива оценка

. (4.20)

Это позволяет утверждать, что решение задачи (4.17) устойчиво по отношению к возмущениям. Предположим, что мы одновременно рассматриваем две "близкие" краевые задачи. Первая задача состоит в отыскании на отрезке 0 ? x ? d решения u (x) системы

,

удовлетворяющего граничным условиям

, .

Во второй задаче требуется на том же отрезке найти решение системы

(4.21)

для которого

,

.

Матрицы , , входящие в граничные условия, предполагаются имеющими то же число строк и столбцов, что и L, R соответственно. Утверждение, что первая и вторая задачи "близкие”, можно понимать как утверждение, что имеют место неравенства

(4.22)

где

- достаточно маленькое число.

Всюду f (x), - непрерывные вектор-функции. Оказывается, что из разрешимости первой задачи при не слишком большом е вытекает разрешимость второй.

Теорема 1. Пусть u (x) - решение правильно поставленной задачи (4.17), для функций Грина, которой справедливы оценки (4.19).

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

Тогда при

е < min{1/ [2K (2 + d)], 1/ (2 + d) } (4.23)

существует и единственно - решение задачи (4.21) и для него верна оценка

, (4.24)

причем µ - число обусловленности задачи равно

µ = K (2 + d) (1 + 2KF),

где

.

(Доказательство - в [10].)

Оценка означает, что решение u (x) краевой задачи непрерывно зависит от коэффициентов A, L, R и правых частей , ш, f (x). В оценку непрерывности входит длина d отрезка, на котором ищется решение, и оценка K норм матриц Грина.

1.4.5 Генерация матриц граничных условий для метода ортогональной прогонки для системы из 55 артерий

Граничные условия выглядят так

Lu=l, Ru=r

В случае 55 артерий L, R - матрицы размера (55Ч110), вектор неизвестных

(110*1), вектора и соответственно (55*1).

Артериальное дерево укладывается определённым образом как показано на рисунке 4.1

Рис.4.1 Укладка артериального дерева.

Рассмотрим, как выписываются матрицы L, R и вектора и .

Рассмотрим k-ый шаг по времени.

Граничное условие для первой артерии выражается формулой (4.13), таким образом первая строчка матрицы L будет выглядеть так:

Первый элемент вектора

.

Граничные условия для всех конечных артерий выражаются формулой (4.14), и если m - конечная артерия, то соответствующая ей строчка в матрице будет

И соответствующий ей элемент вектора

.

Граничные условия для бифуркации артерий, выражаются формулами (3.3), (4.15) или в общем случае, когда i артерия разветвляется на j и l

Соответствующие этим условиям строчки матрицы L

,

соответствующие элементы вектора

.

Аналогичным образом составляется матрица R и вектор .

1.4.6 Описание работы программы

Для работы с данной моделью был создан программный код на языке Java для интеграции ее с базой данных BMOND (http://bmond. biouml.org) посредством пакета программ BioUML (http: // www. biouml. org), что позволило создать матрицы для всех 55 сосудов, необходимые для расчетов. BioUML позволяет создавать формализованное описание комплексных систем и их моделирование и обеспечивает удобный пользовательский интерфейс для ввода модели.

Для работы с данной моделью, было создано приложение-плагин "Hemodynamics" для BioUML и реализован пользовательский интерфейс, а также был создан новый тип диаграмм "Arterial tree" (Рис.4.1 и Рис.4.2.). С помощью этого типа диаграмм можно создавать формальное представление ССС как системы графов, используя сосуды как ребра графа и точки их ветвления (node) как вершины графа. Первой вершиной является сердце, от которого отходит первый сосуд (аорта) и дальше строится всё артериальное дерево. Сосудистое дерево является бинарным, то есть от одного сосуда могут ответвляться только два сосуда.

После создания схемы сосудистого древа были заданы дополнительные параметры, необходимые для дальнейших вычислений. Входными данными для программы являются начальный и конечный моменты времени T (tStart, tFinal), шаг разбиения по времени (tDelta), число разбиений одного сосуда по длине (Number of segmentation) и число разбиений сосуда для интегрирования (Segmentation for integrating), задаваемые в закладке Parameters (Рис.4.2.). В закладке Table of vessels задаётся длина сосуда, площадь его сечения (в начальный момент времени) и коэффициент (Рис.4.1.).

Выходными данными программы являются два массива матриц - A [i], Q [i], где сохранено значение площади сечения сосудов и потока в каждой точке разбиения по времени и по пространству, которые используются BioUML для построения графиков зависимости давления в заданных точках заданных сосудов от времени.

Рис.4.1 Пользовательский интерфейс пакета программ BioUML для плагина "Hemodynamics" и типа диаграмм "Arterial tree". Активна закладка Table of vessels.

Рис.4.2 Пользовательский интерфейс пакета программ BioUML для плагина "Hemodynamics" и типа диаграмм "Arterial tree". Активна закладка Parameters.

1.4.7 Описание структуры программы

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

вершины: сердце, точка разветвления сосудов, контрольная точка (в которой можно наблюдать изменение давления);

рёбра - сосуды.

Свойства сосудов: длина, радиус, коэффициент .

Свойства точки ветвления: номер сосуда, который разветвляется; номера двух дочерних сосудов.

Система BioUML обеспечивает пользовательский интерфейс для графического представления редактирования структуры графа. Для удобства редактирования параметров вершин графа (параметры сосудов) в виде таблицы создана специальная вкладка.

Таблица 4.1 UML диаграмма, описывающая структуру классов в пакете "Hemodynamics”.

Сплошные стрелки - включение одного объекта в другой. Пунктирные стрелки - реализация интерфейса. Красные - использование одним объектом другого.

Таблица 4.1 UML диаграмма, описывающая структуру классов в пакете "Hemodynamics”. (Продолжение.)

Описание работы алгоритма

При создании графа сосудистого древа в BioUML, в программе автоматически создаётся объект класса ArterialTreeDiagramType, затем с помощью метода convert класса DiagramToArterialTreeConvertor из диаграммы генерируется объект (atm) класса ArterialBinaryTreeModel.

После того, как мы запускаем моделирование (нажимаем на кнопку "Start Model"), запускается метод solve из класса HemodynamicsModelSolver с параметром atm.

В методе solve реализована ортогональная прогонка, которая в том числе использует метод QR-разложения матрицы (QRDecomposition).

Для генерации матриц граничных условий L и R используется класс ArterialBinaryTreeMatrixResolver, в котором есть функция resolveMatrix, использующая метод traverseTree, реализующий интерфейс TreeTraversal. С Помощью TreeTraversal реализуется обход дерева.

Чтобы сгенерировать программный код для Matlab, в тесте создаётся объект класса ArterialBinaryTreeModel, для него запускается метод getMatlabPresentation класса ArterialBinaryTreeMatlabGenerator, после запуска теста программа автоматически генерирует m-файл, который в свою очередь является готовой программой на Matlab'е.

Вторая часть таблицы 4.1 демонстрирует реализацию пользовательского интерфейса. Класс HemodynamicsViewPart содержит ArterialBinaryTreePane - панель "Hemodynamics” в BioUML. На этой панели есть две вкладки ”Parameters" и ”Table of Vessels”, которые реализованы в классах ParameterListTab и TableOfVesselsTab соответственно.

1.5 Тестовые расчёты

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

1.5.1 Моделирование гемодинамики в норме

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

Рис.5.1 Графики давления крови в восходящей ветви аорты (красным), дуге аорты (синим) и плечеголовной артерии (зеленым).

Рис.5.2 Графики давления крови в восходящей ветви аорты (желтым), правой сонной артерии (красным) и ее внешней (зеленым) и внутренней (синим) ветвях.

Рис.5.3 Графики давления крови в восходящей ветви аорты (зеленым), в левой (красным) и правой (синим) почечных артериях.

Рис.5.4 Графики давления крови в восходящей ветви аорты (желтым), бедренной (красным) левой задней голенной (синим) и левой передней голенной (зеленым) артериях.

1.5.2 Моделирование физической нагрузки и тахикардии

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

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

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

Далее приведены графики давления в некоторых артериях в зависимости от времени (Рис 5.6-5.9), которые были получены с помощью одномерной модели, в которой на вход был подан профиль давления, приведённый на Рис.5.5.

Рис.5.6 Графики давления крови в восходящей ветви аорты (красным), дуге аорты (синим) и плечеголовной артерии (зеленым) при физической нагрузке (начинается с 30 сек.).

Рис.5.7 Графики давления крови в восходящей ветви аорты (желтым), правой сонной артерии (красным) и ее внешней (зеленым) и внутренней (синим) ветвях.

Рис.5.8 Графики давления крови в восходящей ветви аорты (зеленым), в левой (красным) и правой (синим) почечных артериях при физической нагрузке (начинается с 30 сек.).

Рис.5.9 Графики давления крови в восходящей ветви аорты (желтым), бедренной (красным) левой задней голенной (синим) и левой передней голенной (зеленым) артериях при физической нагрузке (начинается с 30 сек.).

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

Затем были проведены эксперименты по моделированию тахикардии. Входной профиль давления был взят из результата модели Солодянникова, и подставлен на вход в одномерную модель гемодинамики. Результаты эксперимента приведены в дипломной работе Семисалова Б.

1.5.3 Применение метода итераций для оптимизации тестовых расчётов

В начальный момент времени мы задаём нереалистичные условия (Q (0) =0, A (z) =const, , d - длина сосуда). Поэтому выход на действительное решение происходит после определённого промежутка времени. Если применить метод добавления итераций, то выход происходит быстрее.

Рассмотрим уравнение (4.1)

(5.1)

,

где . Тогда уравнение (5.1) будет следующим

. (5.2)

В пункте 4.1 мы использовали уравнение

, (5.3)

которое является приближением уравнения (5.2).

В модель вместо уравнения (5.3) можно ввести метод итераций

1 шаг: вычисляем также, как и в пункте (4.1) по формуле (5.3)

. (5.4)

2 шаг: вычисляем по формуле:

. (5.5)

n шаг: вычисляем по формуле:

. (5.6)

Для (5.5) и (5.6) получается дифференциальное уравнение, аналогичное (4.4)

, (5.7)

где матрица G имеет вид

(4.5)

а вектор f есть

. (4.6)

Рис.5.6 Графики давления в сосудах при шаге разбиения по времени =0.02, одна итерация.

Рис.5.7 Графики давления в сосудах при шаге разбиения по времени =0.02, две итерации.

Заключение

В ходе выполненной работы были получены результаты:

· Созданы компьютерные программы на Java и на MatLab, реализующие одномерную модель гемодинамики с помощью метода прямых и метода ортогональной прогонки, позволяющие проводить моделирование и исследование артериальной системы человека.

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

· Данная программа в виде плагина интегрирована в пакет программ Biouml, созданный в Институте Системной Биологии, предназначенный для формального описания биологических систем, создания моделей и их изучения и способный работать с различными базами данных.

· Проведены эксперименты с созданной моделью, в ходе которых рассматривались различные состояния сердечно-сосудистой системы:

o нормальное состояние в покое;

o при физической нагрузке;

o в случае патологического состояния - тахикардии.

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

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

гемодинамика одномерная модель дифференциальная

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

[1] Lamponi D. One dimensional and multiscale models for blood flow circulation. // Lausanne, EPFL, 2004.

[2] Formaggia L., Veneziani A. Reduced and multiscale models for the human cardiovascular system. // Von Karman, 2003. Lecture notes of the 7th VKI lecture series, "Biological Fluid Dynamics”

[3] Sherwin S., Franke V., Peiro J., Parker K. One-dimensional modeling of a vascular network in space-time variables. // Journal of Engineering Mathematics, Volume 47, Numbers 3-4, December 2003, pp.217-250 (34).

[4] Педли Т. Гидродинамика крупных кровеносных сосудов. // М: Мир, 1983.

[5] Павловский Ю.Н., Регирер С.А., Скобелева И.М. Гидродинамика крови // Итоги науки. М.: ВИНИТИ, 1970. Сер. Гидромеханика, 1968.

[6] Quarteroni A. Modeling the Cardiovascular System - A Mathematical Adventure: Part I.

// Siam News, Volume 34, Number 5.2001

[7] Абакумов М.В., Гаврилюк К.В., Есикова Н.Б., Кошелев В.Б., Лукшин А.В., Мухин С.И., Соснин Н.В., Тишкин В.Ф., Фаворский А.П. Математическая модель гемодинамики сердечно-сосудистой системы. // Дифференциальные уравнения. 1997, 33 (7), с.892-898.

[8] Quarteroni A., Formaggia L. Mathematical modelling and numerical simulation of the cardiovascular system. // In: Ciarlet P. G. (ed.). Handbook of numerical analyis, V.12, special volume: Ayache N. (ed.).computational Models for the Human Body. Amsterdam, Elsevier Science & Technology, 2003,P.3-129.

[9] Бибердорф Э.А., Попова Н.И. Глава 3: Формальное описание и моделирование сердечно-сосудистой системы. // А.М. Блохин, Л.Н. Иванова и др. Система кровообращения и артериальная гипертония: биофизические и генетико-физиологические механизмы, математическое и компьютерное моделирование. (Монография). Издательство СОРАН, Новосибирск, 2008. (В печати).

[10] Кузнецов С.В. Развитие метода ортогональной прогонки.

Диссертация. // Институт математики СО РАН, Новосибирск: 1988.

[11] Неймарк М.А. Линейные дифференциальные операторы. // М: Наука, 1969.

Приложение

1. Пример сгенерированного Matlab кода

function [Energy] =hemodynamicModelSimulation ()

Nfun=2; %***********number of unknown functions*****************

M=5; %*************segmentation for integrating*******************

%******************physical parameters**************************

rho=1;

h_=0.05;

nu=0.035;

E=3*10^6;

N=10; %***********number of segmentation************************

tau=0.05; %*********time step***********************************

t=0; %***********start time************************************

T=50; %***********end time************************************

Nvet=55; %********number of vessels*****************************

%******************vessels length******************************

length (1) =4.0;

length (2) =2.0;

length (54) =32.3;

length (55) =34.4;

%******************vessels initial area***************************

A0 (1) =5.983;

A0 (2) =5.147;

A0 (54) =0.11;

A0 (55) =0.06;

%*************sign*******************************************

znak= [1; - 1; 1; - 1; 1; 1; - 1; 1; 1; - 1; - 1; 1; 1; - 1; - 1; 1; - 1; 1; - 1; - 1; 1; - 1; 1; - 1; 1; - 1; - 1; 1; - 1; - 1; 1; - 1; - 1; 1; - 1; 1; 1; - 1; - 1; 1; 1; 1; - 1; - 1; - 1; 1; - 1; - 1; 1; 1; - 1; - 1; 1; - 1; - 1];

%******************derived physical parameters*******************

gamma (1) =388.0;

gamma (2) =348.0;

gamma (54) =23232.0;

gamma (55) =36972.0;

Kr=8*pi*nu;

%******************segmentation*******************************

for j=1: Nvet,

h (j) =length (j) /N;

end;

%******************dimension*********************************

n_a=Nfun*Nvet;

%******************solution***********************************

V=zeros ([Nfun, (N+1),Nvet]);

for i=1: N+1,for j=1: Nvet,

V (1, i,j) =A0 (j);

end;

end;

%*************time iteration***********************************

Kmax=round (T/tau);

Energy=zeros (Kmax);

for pr=1: 4,for k=1: Kmax,

t=t+tau;

Ordinata (k) =t;

%*************marching **************************************

%*************left boundary condition***************************

for j=1: Nvet,

An (j) =V (1,1,j);

Qn (j) =V (2,1,j);

a (j) =gamma (j) / (A0 (j) * (sqrt (An (j)) +sqrt (A0 (j))));

b (j) = (rho*Qn (j)) / (2*An (j) ^2);

c (j) =gamma (j) / (sqrt (An (j)) +sqrt (A0 (j)));


Подобные документы

  • Разностная схема решения уравнения теплопроводности. Численное решение уравнения теплопроводности в табличном процессоре Microsoft Ехсеl и в пакете математических расчётов MathCAD. Расчёт методом прогонки. Изменение пространственной координаты.

    дипломная работа [248,4 K], добавлен 15.03.2014

  • Разработка программы на языке С++ для решения дифференциального уравнения Лапласа в прямоугольной области методом сеток. Численное решение задачи Дирихле для уравнения Лапласа, построение сетки и итерационного процесса. Листинг и результат программы.

    курсовая работа [307,5 K], добавлен 30.04.2012

  • Решение задачи Коши для дифференциального уравнения методом Рунге-Кутта и Адамса с автоматическим выбором шага и заданным шагом. Интерполирование табличной функции. Численное решение системы линейных алгебраических уравнений методами простой итерации.

    методичка [35,8 K], добавлен 15.03.2009

  • Численное решение задачи Коши для обыкновенного дифференциального уравнения первого и второго порядка методом Эйлера и Рунге-Кутты и краевой задачи для ОДУ второго порядка с применением пакета MathCad, электронной таблицы Excel и программы Visual Basic.

    курсовая работа [476,2 K], добавлен 14.02.2016

  • Математическое описание алгоритмов схемы и операций для уравнения Лапласа. Изучение разностной схемы "крест" для нахождения численного решения эллиптического уравнения, задача Дирихле. Использование указателей в среде Matlab для решений методом Гаусса.

    дипломная работа [859,3 K], добавлен 23.10.2014

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

    курсовая работа [1,7 M], добавлен 09.06.2012

  • Разработка прикладного программного обеспечения для решения расчетных задач для компьютера. Численное интегрирование - вычисление значения определённого интеграла. Проектирование алгоритма численного метода. Тестирование работоспособности программы.

    курсовая работа [1,1 M], добавлен 03.08.2011

  • Расчет тепловой схемы с применением методов математического моделирования. Разработка алгоритма реализации модели. Составление программы для ПЭВМ, ее отладка и тестирование. Проведение численного исследования и параметрическая оптимизация системы.

    курсовая работа [2,8 M], добавлен 01.03.2013

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

    курсовая работа [67,9 K], добавлен 07.12.2009

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

    курсовая работа [79,2 K], добавлен 25.06.2011

Работы в архивах красиво оформлены согласно требованиям ВУЗов и содержат рисунки, диаграммы, формулы и т.д.
PPT, PPTX и PDF-файлы представлены только в архивах.
Рекомендуем скачать работу.