Одномерная модель описывает течение крови в артериях и ее взаимодействие с подвижными стенками. Пусть 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 Анализ математических свойств модели

Рис.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 крупных артерий тела человека.

Рис.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) используем метод ортогональной прогонки.

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

, , , (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 Тестовые расчёты