Содержание к диссертации
Введение
1. Математическое моделирование гемодинамики одиночного сосуда
1.1 Математическая модель
1.2 Численный метод
1.3 Тестовые расчеты
Задача с правыми частями специального вида задача о распаде произвольного разрыва
Моделирование сердечного выброса. линейная модель .
Моделирование сердечного выброса. нелинейная модель
2. Многомасштабное моделирование
2.1 Нульмерная модель системы кровообращения в целом
2.2 Сопряжение одномерной модели сосуда с нульмерной моделью кровеносной системы
2.3 Сопряжение двумерной модели аорты с нульмерной моделью кровеносной системы
3. Моделирование течений крови в артериальном дереве. исследование эффектов, возникающих при пережатии сосуда манжетой
3.1 Моделирование гемодинамики в артериальном дереве в невозмущенных условиях
3.2 Моделирование пережатия плечевой артерии в артериальном дереве
Заключение
Список литературы
- Численный метод
- Моделирование сердечного выброса. линейная модель
- Сопряжение одномерной модели сосуда с нульмерной моделью кровеносной системы
- Моделирование пережатия плечевой артерии в артериальном дереве
Введение к работе
На сегодняшний день, благодаря накоплению специальных знаний и опыта, математическое моделирование стало мощным инструментом анализа процессов, протекающих в системе кровообращения. В особой мере это касается изучения комплексного воздействия разнообразных факторов на характеристики системы. Возможности прямого измерения при этом, как правило, ограничены, в то время как построенная с использованием доступных экспериментальных данных математическая модель позволяет обеспечить подробную детализацию и оценить взаимное влияние различных параметров друг на друга, а также на функционирование системы в целом. Важнейшим условием эффективности вычислительного эксперимента является адекватность математической модели протекающим в системе кровообращения физическим процессам.
Укрупнешю иерархия математических моделей гемодинамики, используемых в современной вычислительной практике, может быть представлена следующим образом.
1) Многомерные модели.
Течение крови в системе кровообращения в общем случае описывается трехмерными нестационарными уравнениями для вязкой в общем случае неньютоновской жидкости совместно с уравнениями динамики эластичных оболочек сосудов [1, 2]. Это связано с необходимостью учета реальных свойств крови, пространственной геометрии сосудов, влияния вязкости, взаимного влияния гидродинамики сосудов и их деформации. В камерах сердца и крупных кровеносных сосудах, где характерные размеры микроструктур значительно меньше характерных масштабов течения, кровь обычно рассматривается как однородная ньютоновская среда. Исключение могут составлять лишь области с малыми сдвиговыми напряжениями, в которых реализуются условия, благоприятные для агрегации микроэлементов. Задача моделирования течения крови в значительной мере облегчается тем обстоятельством, что практически во всех отделах кровеносной системы наблюдается ламинарный режим течения. Исключение составляет лишь часть аорты, где во время выброса крови из сердца наблюдается турбулентное течение, которое, однако, быстро переходит в ламинарное. Таким образом, одной из основных проблем построения вычислительного алгоритма является необходимость решения уравнений Навье-Стокса в областях с подвижными криволинейными границами. Этому вопросу посвящена, в частности, работа [3].
Другой проблемой являются чрезвычайно высокие требования, предъявляемые к вычислительным ресурсам [4, 5]. Положение несколько облегчается, если течение рассматривается как двумерное.
Однако вследствие необходимости учета импульсного нестационарного характера функционирования системы кровообращения вычислительные затраты и в этом случае остаются достаточно высокими, даже при использовании эффективных численных алгоритмов. В связи с вышесказанным многомерные модели не применяются для описания системы кровообращения в целом, а используются для получения детальной картины течения в характерных локальных зонах. Постановка граничных условий при этом должна отражать взаимное влияние рассматриваемого участка и других частей кровеносной системы. В этой связи все большее распространение получают многомасштабные модели, о которых будет сказано ниже.
2) Одномерные модели.
Одномерные модели предполагают работу с осредненными по поперечному сечению потока гемодинамическими параметрами. Поскольку течение в кровеносном сосуде направлено, главным образом, вдоль его оси, допущение об одномерности течения часто используется при исследовании волновых процессов в кровеносной системе. В этом случае движение крови по сосуду описывается гиперболической системой нелинейных дифференциальных уравнений в частных производных, выражающих законы сохранения массы и импульса. Вязкие эффекты описываются входящим в правую часть уравнения сохранения импульса выражением для силы трения. С целью учета взаимодействия течения и стенки сосуда система, как правило, замыкается полуэмпирической зависимостью площади поперечного сосуда от давления в потоке. Другой подход заключается в совместном решении уравнений гемодинамики и уравнений динамики оболочки сосуда. Проблемы динамического поведения кровеносных сосудов, как деформируемых оболочек рассмотрены в работе [6].
Использование одномерной модели для анализа течения в каком-либо отдельно взятом сосуде сопряжено с проблемой постановки граничных условий. В особой мере это касается правой, выходной, границы. Одним из возможных вариантов решения этой проблемы, также как и для многомерных моделей, является многомасштабный подход. Другой, более точный, но требующий более высоких вычислительных затрат, путь заключается в решении задачи для всей кровеносной системы в целом при использовании на каждом из участков одномерной модели с характерными для него параметрами. Несмотря на сложность решения задачи в такой постановке, в настоящее время имеются успешные реализации подобного моделирования. Наиболее полное комплексное решение проблемы предложено в работах [7, 8]. Здесь система кровообращения описана графом, ребра которого соответствуют отдельным кровеносным сосудам кровеносной системы или жгутам функционально однородных мелких сосудов, а вершинам графа приписаны функциональные свойства участков ветвления сосудов, мышечных тканей или отдельных органов. В качестве краевых условий для такой модели используются параметры сердечного выброса и разрежение, создаваемое сердцем в венозной части. Условия на границах сосудов получаются автоматически в процессе решения. При этом в узлах ветвления сосудов используются соотношения, выражающие закон сохранения массы и постоянства интеграла Бернулли. Аналогичная модель, но с меньшей степенью детализации, представлена в [9]. Модель охватывает 55 артерий, в каждой из которых течение описывается одномерной моделью. В отличие от [8] здесь не рассматривается венозная часть, а в выходном сечении артериального дерева задаются неотражающие условия.
При реализации рассмотренных комплексных моделей возникают два сорта проблем. Первый - обусловлен повышенной чувствительностью модели к заданию начальных данных, параметров сосудов и узлов сопряжения [8], второй - связан с необходимостью достижения определенной точности решения при разумных затратах вычислительных ресурсов. Последним обстоятельством продиктовано стремление как к использованию распараллеливания вычислений [9], так и к повышению эффективности вычислительных алгоритмов решения уравнений гемодинамики для одиночного сосуда. Именно этому вопросу уделяется основное внимание в первой главе диссертации, где рассматриваются различные аспекты применения к решению уравнений гемодинамики TVD-монотонизированных схем второго порядка точности по пространству и времени.
3) нульмерные (дискретные модели).
Данный класс моделей применяется для описания системы кровообращения в целом, которая представляется последовательностью характерных участков с соответствующим набором переменных по времени параметров течения (давление, расход крови и т.д.) [напр., 10]. Появление подобных моделей восходит к работам, посвященным аналоговому моделированию системы кровообращения с использованием электрических цепей [напр., И, 12, 13]. Поскольку модели данного класса не имеют пространственного разрешения, они обычно называются нульмерными моделями. Математически такая модель представляет собой систему обыкновенных дифференциальных уравнений, в результате решения задачи Коши для которой, определяются временные зависимости осредненных по пространству гемодипамических параметров на характерных участках кровеносной системы. Несмотря на то, что нульмерные модели не дают подробной картины течения, интерес к ним не угасает, что в значительной мере обусловлено их активным использованием в многомасштабпом моделировании.
4) Многомасштабные модели
Как уже отмечалось выше, одной из принципиальных проблем, возникающих при исследовании гемодинамики с использованием многомерных и одномерных моделей, является постановка граничных условий.
С одной стороны, необходимость учета взаимного влияния различных частей кровеносной системы зачастую приводит к невозможности изолированного рассмотрения какого-либо участка, а с другой - применение комплексных моделей, охватывающих систему в целом, связано с высокими, а в случае многомерных моделей - непомерными затратами вычислительных ресурсов. В этой связи разумным компромиссом является сопряжение моделей различной размерности. Схематично такой подход проиллюстрирован на рис.1, заимствованном из работы [14].
Здесь одномерная модель гемодинамики аорты и трехмерная модель, описывающая течение крови в зоне бифуркации (разветвления) сонной артерии, сопряжены с нульмерной моделью кровеносной системы в целом. Это дает возможность автоматически получить корректные краевые условия на границах аорты и зоны бифуркации и в тоже время исследовать ответную реакцию различных участков кровеносной системы на процессы, происходящие в этих областях. Отметим также работу [16], в которой проведено сопряжение одномерной и трехмерной моделей сосуда, что позволило существенно уменьшить нефизические эффекты в трехмерной модели, связанные с отражениями пульсовых волн.
Однако, в ряде случаев, когда необходимо получение детальной картины течения в определенном сосуде, местоположение которого на упрощенной нульмерной модели невозможно определить, а также, когда необходима подробная информация о гемодинамических параметрах в кровеносной системе, использование многомасштабной модели является неэффективным. В частности, с такими проблемами приходится сталкиваться при моделировании пережатия плечевой артерии манжетой в условиях измерения давления аускультативным или осциллометрическим методами. С целью детального анализа подобных процессов в настоящей работе применяется модификация модели артериального дерева, предложенной в работе [9]. Такой подход позволяет не только получить временные распределения гемодинамических параметров, но и детально изучить пространственную картину течения в условиях пережатия сосуда.
Немаловажным аспектом является выбор высокоточного, экономичного алгоритма решения задачи гемодинамики.
В связи с тем, что уравнения гемодинамики относятся к классу гиперболических уравнений, приведем краткий обзор численных методов решения задач рассматриваемого класса.
Исторически первыми были созданы методы первого порядка точности. Обзор методов первого порядка точности может быть найден в работах [17-20]. Среди таких методов выделим метод, предложенный Годуновым С.К., основанный на использовании решения задачи о распаде произвольного разрыва. В дальнейшем на базе этого метода был развит класс численных схем типа Годунова, в которых применяются как точное, так и приближенные решения задачи Римана.
Схемы первого порядка точности обладают сильными диссипативными свойствами, что приводит к «размазыванию» решения в областях с резким изменением расчетных параметров и, как следствие, вынуждают проводить расчеты на очень мелких сетках. Поэтому дальнейшее развитие методов сквозного счета было направлено на повышение порядка точности. В этой связи отметим известные схемы второго порядка точности - Лакса-Вендроффа, Мак-Кормака, схема «чехарда» и другие.
Среди неявных схем можно отметить такие как схема Бима-Уорминга, Кранка-Николсона. Несмотря на то что, как правило, неявные методы решения гиперболических уравнений безусловно устойчивы, тем не менее их применение не обеспечивает преимущество относительно явных схем. При использовании неявной схемы на каждом временном слое получается система связанных между собой уравнений. Таким образом, любое возмущение (например, ошибка округления) влияет на решение во всех узлах на следующем временном слое.
Затем были созданы схемы третьего порядка точности [21, 22]. Также были предложены схемы четвертого и более высоких порядков точности [23]. Подробный обзор развития явных разностных схем приведен в работе [24]. Сравнительные исследования свойств схем могут быть найдены в работах [25, 26].
Увеличение порядка точности разностной схемы приводит к тому, что область, на которой «размазана» зона с резким изменением расчетных параметров уменьшается. Однако вместе с тем повышение порядка приводит к возникновению ложных осцилляции решения в таких областях. Как известно из теоремы Годунова С.К., в линейном случае монотонность можно обеспечить только в разностных схемах первого порядка точности. Поэтому на пути дальнейшего развития методов основной проблемой при построении разностных схем стало выполнение следующего условия: численная схема не должна обладать сильными диссипативными свойствами и должна сходиться к физически корректному решению. То есть одновременно с повышением порядка точности назрела необходимость создания механизма подавления ложных осцилляции, обусловленных дисперсионными свойствами схем.
Одним из способов их устранения являются различные искусственные добавки, такие, например, как, искусственная вязкость [27].
Обзоры такого рода подходов приведены в [28, 29, 30]. Обобщение на многомерный случай рассмотрено в обзоре [31].
К недостаткам этой группы методов можно отнести, то, что ее применение зачастую приводит к существенным изменениям в решении [32, 33] . Кроме того, такого рода методы требуют тщательного подбора коэффициентов.
Другой способ, позволяющий избежать ложных осцилляции решения заключается в применении схем переменного порядка точности. В областях, где имеет место резкое изменение расчетных параметров могут быть использованы схемы первого порядка точности, там же где имеет место гладкое решение возможно использование схем более высокого порядка точности.
Развитие этого подхода, названного гибридным, восходит к работам отечественных авторов [34, 35]. Как правило, под ним понимают применение нелинейных схем, которые могут локально менять свои свойства, например, порядок аппроксимации, в зависимости от структуры решения. В частности, гибридная схема может осуществлять сквозной счет схемой второго порядка точности в областях, где имеет место гладкое решение и переходить на схему первого порядка точности в тех областях, где решение содержит большие перепады расчетных параметров. Такой подход позволяет избежать нефизических осцилляции решения, свойственных для традиционных схем второго порядка точности и вместе с тем избежать «размазывания» решения, характерного для методов первого порядка.
Простейшая гибридная схема представляет из себя комбинацию двух схем: kSl + (1 - k)S2. Здесь Si — это схема первого порядка точности, . — схема второго порядка, к — коэффициент гибридности, где 0 к 1.
Такая схема была пременена в [34] для линейного уравнения переноса.
Переключение между Si и S2 происходило на основе анализа отношения второй конечной разности решения к ее первой разности.
Затем было предложено обобщение гибридности на случай системы уравнений [36]. Применялась комбинация схем Лакса-Фридрихса [18] первого порядка точности и схемы Лакса-Вендроффа второго порядка точности [37, 38].
Были созданы и другие гибридные схемы на переменном разностном шаблоне, которые в зависимости от характера течения используют либо центральные, либо ориентированные разности [39,40,41,42].
Гибридным неявным разностным схемам с учетом характеристических направлений посвящены работы [43, 44]. Гибридизация явных схем с учетом характеристических направлений обсуждается в работах [45, 46,47,48].
В работах [49, 50, 51] предложен гибридный метод, позволяющий повысить порядок точности при помощи процедуры коррекции потоков. Основная идея такого рода методов состоит в следующем. На первом шаге (предиктор) применяется схема, в которую включена искусственная вязкость. На втором шаге, называемом антидиффузионным, решение частично уточняется. При этом антидиффузионные поправки ограничиваются таким образом, чтобы они не вносили новых экстремумов в решение, а также не увеличивали (не уменьшали) значения локальных максимумов (минимумов), которые имели место на первом шаге вычислений. Такие условия фактически эквивалентны условиям неувеличения полной вариации численного решения (TVD, total variation diminishing). Хотя этот метод показал на практике свою эффективность, он, тем не менее, не имеет строгого теоретического обоснования. Позднее были развиты схемы, которые опираются на TVD-условие более явно. Дальнейшее развитие монотонизированных схем связано с работами Хартена [52]. Он предложил использовать специальную кусочно-линейную (кусочно полиномиальную) реконструкцию сеточных функций, при помощи которой достигается выполнение TVD-условия для численного решения.
Наклоны расчетных параметров внутри дискретных ячеек расчетной сетки для выполнения TVD-свойства ограничиваются специальными функциями -ограничителями (limiters). Обзор ограничителей приведен в работах [53, 54].
Большинство современных схем переменного порядка точности используют кусочно-полиномиальную реконструкцию сеточных функций, удовлетворяющую TVD-условию.
Дальнейшее развитие TVD-подхода привело к созданию ряда модификаций, схем UNO [55], ENO [57], WENO [58] и др. Основная идея схем UNO состоит в том, что количество локальных экстремумов N[q] в сеточной функции q не возрастает со временем: N[qk+]] N[qk]. Схемы такого рода допускают возрастание полной вариации сеточной функции в пределах порядка точности схемы. Основная идея ENO и WENO-схем состоит в том, что реконструкция функции осуществляется на переменном шаблоне. Шаблон выбирается из соображений локальной гладкости решения и отсутствия ложных осцилляции в окрестности разрывов. ENO - схемы используют только один (оптимальный в некотором смысле) шаблон для реконструкции функции, в то время как WENO-схемы используют выпуклую комбинацию всех возможных шаблонов с весовыми коэффициентами. Коэффициенты определяются исходя из локальной гладкости решения, полученного с использованием данного шаблона. Появление нового класса методов привело к значительному улучшению качества численных решений по сравнению с классическими разностными методами фиксированного порядка точности. Схемы, основанные на TVD свойстве, а также их модификации относятся к классу схем высокого разрешения.
Однако, использование ENO, WENO-схем а также ряда других модификаций повышенного порядка точности сопряжено с определенными сложностями. Реализация такого рода схем требует значительных вычислительных затрат, что обусловлено перебором шаблонов. Кроме того, поскольку условие невозрастания вариации решения строго не соблюдается, в решении зачастую возникают нефизические эффекты. В связи с вышесказанным в качестве базового подхода к построению вычислительных алгоритмов в настоящей работе применялись TVD-монотонизированные схемы повышенного порядка точности.
Завершая обзор математических моделей гемодинамики и численных методов решения гиперболических систем, сформулируем основные цели работы:
разработка алгоритмов для численного моделирования гемодинамических процессов в системе крупных артерий;
компьютерная реализация многоуровневой системы математических моделей гемодинамики, разработка комплекса программ;
численное моделирование эффектов, возникающих в результате пережатия сосуда манжетой при измерении давления аускультативным и осциллометрическим методами. В первой главе рассмотрены вопросы математического моделирования течений крови в одиночном сосуде. Для рассматриваемой задачи применялись TVD
монотонизированные схемы, обеспечивающие второй порядок и позволяющие избежать ложных осцилляции решения.
В первой главе приведены результаты численного моделирования ряда характерных тестовых задач. Рассмотрены задача с правыми частями специального вида, задача о распаде произвольного разрыва и задача о моделировании сердечного выброса. Исследована сходимость и рассмотрены вопросы выбора оптимальных сеточных параметров.
Во второй главе рассматриваются вопросы многомасштабного моделирования гемодинамики. Для постановки адекватных граничных условий необходимо учитывать взаимосвязь исследуемого сосуда с сопряженными участками. Один из возможных путей решения этой проблемы - сопряжение уравнений гемодинамики одиночного сосуда с уравнениями, описывающими кровеносную систему в целом.
Рассматривается нульмерная модель сердечно-сосудистой системы, включающая шесть основных участков, каждый из которых характеризуется двумя параметрами - артериальным давлением и расходом крови. Применен подход, позволяющий с помощью нульмерной модели получать граничные условия для одномерной модели сосуда. Этот подход заключается в том, что одномерная модель сосуда «врезается» в нульмерную модель сердечнососудистой системы на определенном участке. Кроме того, рассмотрена задача сопряжения двумерной модели аорты с нульмерной моделью кровеносной системы. Такой подход позволил получить нестационарные профили скорости течения крови и провести сравнительный анализ результатов полученных для одномерной и двумерной модели сосуда.
В третьей главе применяется полномасштабный подход к моделированию гемодинамики. Рассматривается модель, включающая 55 основных крупных артерий человека, в каждой из которых течение описывается одномерной моделью. Получены решения для невозмущенных условий, проведен сравнительный анализ результатов.
С использованием модификации модели гемодинамики, основанной на переменной ригидности сосудистой стенки, осуществлено численное моделирование процессов, протекающих в сердечно-сосудистой системе при измерении артериального давления аускультативным или осциллометрическим методами. Анализируется детальная картина течения. Получена характерная для осциллометрических методов колоколообразная зависимость амплитуды колебаний площади поперечного сечения плечевой артерии от пережимающего давления. На основе полномасштабной модели артериального дерева показано, что возникновение тонов Короткова связано с нелинейным характером поведения упругих свойств оболочки сосуда. Научная новизна Разработаны эффективные вычислительные алгоритмы моделирования гемодинамики в системе крупных артерий. Предложена модификация модели гемодинамики на случай пережатия сосуда манжетой, основанная на переменной ригидности сосудистой стенки. Исследована детальная структура течения при пережатии. На основе полномасштабной модели артериального дерева показано, что возникновение тонов Короткова при измерении давления аускультативным методом связано с нелинейным характером поведения упругих свойств оболочки сосуда. Получена колоколообразная зависимость амплитуды колебаний площади поперечного сечения плечевой артерии от пережимающего давления, характерная для осциллометрических методов измерения артериального давления.
Апробация работы
Материалы диссертации докладывались и обсуждались на семинаре кафедры «Вычислительная математика и программирование» под руководством чл.-корр. РАН, профессора Пирумова У.Г., на XIX международном семинаре по струйным, отрывным и нестационарным течениям (Санкт-Петербург, 2002), на IV, V, VI международных конференциях по неравновесным процессам в соплах и струях (Санкт-Петербург, 2002, Самара, 2004, Санкт-Петербург, 2006), на XII, XIV международных конференциях по вычислительной механике и современным прикладным программным системам (ВМСППС) (Владимир, 2003, Алушта, 2005), на семинаре отдела новых методов диагностики НИИ кардиологии им. А.Л. Мясникова, на семинаре кафедры «Вычислительные методы» факультета вычислительной математики и кибернетики МГУ. Публикации
По материалам настоящей диссертационной работы опубликованы тезисы IV, V международных конференциях по неравновесным процессам в соплах и струях, XII и XIV международных конференциях по вычислительной механике и современным прикладным программным системам [60-64], также 3 статьи [65, 66, 67]
Численный метод
Как отмечалось во введении, стремление к описанию с помощью одномерной модели всей кровеносной системы стимулирует поиск высокоточных алгоритмов для решения задачи гемодинамики одиночного сосуда. При различных режимах функционирования в сосудах возможно появление зон со значительным продольным изменением гемодинамических параметров. Как следствие, применение схем первого порядка по пространству является малоэффективным в силу сильной схемной диссипации. В то же время традиционные схемы второго порядка обладают собственными дисперсионными свойствами, проявляющимися в виде ложных осцилляции численного решения в областях больших градиентов. Учитывая, что при определенных условиях в сосуде возможно возникновение высокочастотных колебаний физической природы, исследование которых представляет самостоятельный интерес, эффекты схемной диссипации и схемной дисперсии численного решения являются крайне нежелательными и требуют особого внимания при разработке вычислительных алгоритмов. Так авторами работы [8] использовалась схема второго порядка с искусственными вязкостью и дисперсией, в результате чего им удалось добиться подавления отмеченных негативных эффектов.
В настоящей диссертационной работе для решения системы уравнений гемодинамики (1), (3) используется TVD-подход (Total Variation Diminishing). Как отмечалось во введении, этот подход связан с построением схем, которые уменьшают или сохраняют полную вариацию функции, не допуская тем самым появления ложных осцилляции. Основная идея TVD-подхода состоит в том, что расчет ведется всюду со вторым порядком точности, кроме зон с резким изменением параметров, где схема автоматически переключается на первый порядок точности. Этот переход обеспечивается с помощью специальных функций - лимитеров (ограничителей). Вопросы построения TVD-схем подробно освещены, например, в [75, 76]. Следуя этим работам, рассмотрим общую схему построения вычислительного алгоритма применительно к задаче гемодинамики (10).
В результате интегрирования (10) по /-ой ячейке ( є [ ,.,/2 ,+1/2] ) разностная схема может быть представлена в виде: ".А Ах = -F1+F1+ \fdx at и- - J (13) где qt=— \qdx - среднее по ячейке; индекс /±— соответствует значениям параметров на границах ячейки.
Потоки на границах ячейки определяются из следующего соотношения: ( 2 ( F\+F\ /+- /+ V 2 2 w i+— /+ V 2 и (14) где R\A\R-\\A\ = diag(\Zl\) / = 1,2. Здесь все матрицы вычисляются при q \ ( І+- 1+ \ 2 1J индексы L, к означают, что соответствующие параметры вычисляются в левой и правой ячейках соответственно. Компоненты вектора q на границах ячеек вычисляются с помощью разложений в ряд Тейлора. При этом с целью монотонизации схемы для производных, входящих в разложения, применяются ограничители (limiter): (?Дд=(? ),+-Л , дд„ \дх ), limiter faL-faT fa),-fa ), (15) fa) =faL ,+i + /+i / limiter faL-faL"! fa),+,-fa)/ J к = 1,2 (16) где а? Л fa),-fa L ftc Ax
Данная схема относится к классу схем типа Куранта-Изаксона-Риса [см, напр, 75]. Схемы такого рода основаны на приближенном решении задачи Римана для локально линеаризованной системы уравнений. В работе используются различные ограничители: minmod, van Leer, МС, superbee, МГУ. Обозначим в = . , ,— -. Тогда выражения для fa)/ -fa)/-. ограничителей имеют следующий вид [74,75,76]: minmod : limiter(0) = van Leer: limiter(#) = МС: 0, 6 0 0, O 0 1; 1, 0 \ в + \в\ \ + \в\ = max ( 0, mm , 2, Superbee: limiter(#) = [0, 6 0 max(min(l, 26 ),min(2, в)) МГУ: limiter(#) = max 0, mini 20, - + -, 21 . Для аппроксимации уравнения (13) применяется метод Рунге-Кутта второго порядка. Обозначим Ах, ,- - \fdx г г , Ч "2 J (17) Тогда алгоритм перехода на новый временной слой можно записать в виде: я]=Ъ-Ця"), (18) q)"+l)=0.5(q:+q ,-L(q t)) (19) Такая модификация метода Рунге-Кутта является SSP-методом (Strong Stability Preserving High - order Time Discretization) [77]. Проведенные в настоящей диссертационной работе расчеты показали, что применительно к задачам рассматриваемого класса применение схем с повышенным порядком аппроксимации по времени дает существенный выигрыш по сравнению с методами первого порядка.
Моделирование сердечного выброса. линейная модель
To есть на характеристиках выполняются следующие соотношения: и-щ± \ -ds = Q (32)
Соотношение (32) выполняется в общем случае, для произвольной зависимости p(s) (из физического смысла ясно, что функция p(s) должна быть непрерывно-дифференцируемой, с положительной областью определения и монотонно возрастающей). Будем искать решение в виде q(t, х) - q( ) = q х-х , t . Здесь хс - — - центр веера волн разрежения, в данном случае середина сосуда. dx х-хс dt t-Є =%=и±с (33)
Подставив (33) в (32) получим соотношение, связывающее автомодельную переменную и площадь поперечного сечения сосуда s внутри веера: +с± \ -ds = О (34)
Обобщенно решение задачи о распаде разрыва, используя (29) и (32), можно записать в виде: F(s) = u1-ux+fx(s,sl) + f\s,s2) = 0 (35) \ -ds; s s (36) s-s / = T(p-p \ s s p s + s
Поскольку F(s) монотонно возрастающая функция, то ее исследование позволяет априорно определить конфигурацию, возникающую после распада разрыва: F(s ) 0 - влево веер, F(s ) 0 - влево ударная волна, F(s2) 0 - вправо веер, F(s2) 0 - вправо ударная волна. Решив уравнение F(s)-0, получим значения sup. Подставив их в (28), для ударной волны получим скорость распространения ударной волны, а для веера волн разрежения в (32) получим скорость на крайней характеристике веера. Параметры внутри веера разрежения могут быть рассчитаны из (34) как функции переменной %. В частном случае, для зависимости (24) можно записать (34) в виде: + с±2(с-са)=0 (37) Откуда следует: , (38) Для зависимости (24) можно записать (36) в виде: №-A s s f = s-s i(s + s ) s s (39) 2y Необходимо отметить, что в случае линейной зависимости p(s) можно аналитически получить соотношения, связывающие параметры до и после ударной волны. Обозначим М = W-ua (40)
Несложно получить следующие зависимости: / Л sa 4 М2 \ ! + ,!+ (41) s 2\ \ М1 (42) М1 и = иа+-Маа ( ч 3-J1 + (43)
Также можно получить выражения для числа М отраженной ударной волны в зависимости от параметров за падающей ударной волной: Мов = W-u 1 и М„в=--±\\-г + 16 т 4 4\с2 (44)
Здесь знак «+» соответствует случаю отражения УВ от левой стенки, «-» - от правой. Отражение характеристик от границы расчетной области для веера волн разрежения схематически изображено на рис. 4 (пунктирными линями обозначено продолжение характеристик вне расчетной области). Для областей 1, 2 решением являются начальные условия (26), для областей 3, 4 решение получено выше. Получим решение в областях 5, 6, 7 в предположении, что характеристики в области 5 являются прямыми линиями. Рассмотрим характеристику , = их - с1.
Сопряжение одномерной модели сосуда с нульмерной моделью кровеносной системы
Для моделирования течений вязкой несжимаемой жидкости в осесимметричном сосуде используются уравнения Навье-Стокса в размерном виде. Уравнения количества движения и уравнение неразрывности в цилиндрической системе координат имеют вид: ди ди ди 1 дР и + 11 + V— = + — dt дх дг р дх р + — (д2и д2и 1 ди} 2 а.. 2 дх дг г дг (70) dv dv dv \дР fi(d\ д\ \dv v) — + М— + V— = + — —- + —- + : dt дх дг p дг р\дх дг г дг г (71) ди v dv . —+-+—=0 дх г дг (72)
Здесь и, v - осевая и радиальная составляющая вектора скорости соответственно, Р - давление.
Расчеты проводятся для аорты длиной L = 4 см и переменной площадью поперечного сечения s{p), задаваемой соотношением (24). Уравнения (70)-(76) замыкаются следующими краевыми условиями: ди дх L = 0, ди дх = 0,и =v,(4- дг = 0, (73) дх = 0, dv_ дх :0,vr=vr(4vfl=0, (74) где Ух{х) осевая составляющая скорости движения границы, a vr(jc) -радиальная составляющая скорости движения границы.
С целью осуществления «врезки» двумерной модели аорты в нульмерную модель кровеносной системы проводилось сопряженное решение уравнении (70)-(72) и системы ОДУ (68). Для численного решения уравнений Навье-Стокса в области с подвижной границей использовалась программа, разработанная в работе [81].
В результате расчетов были получены нестационарные поля скорости и давления, а также геометрическая форма и положение подвижной степки сосуда. На рис. 21 приведены профили продольной скорости в серединном сечении сосуда в различные моменты времени. Можно отметить, что в периоды времени, соответствующие диастолической фазе, возникают пристенные возвратные течения, которые затем развиваются в обратные течения по всему поперечному сечению сосуда. Заметим, что полученные профили сильно отличаются от параболических, соответствующих стационарному течению в круглой трубе. В этой связи использование квазистационарного подхода (см. (2)) может рассматриваться лишь в качестве первого приближения к реалиям столь сложного процесса.
На рис.22 приведены временные зависимости средней скорости в срединном сечении сосуда, полученные на различных сетках. На графике видна хорошая согласованность решений, полученных по двумерной и одномерной моделям. Зависимости давления на входном и выходном сечении аорты представлены на рис. 23.
Таким образом, во второй главе осуществлена реализация многомасштабного подхода к численному моделированию гемодинамики. Проведено сопряжение одномерной модели аорты с нульмерной моделью кровеносной системы. Полученные результаты достаточно хорошо воспроизводят качественную картину эволюции давления. Кроме того, рассмотрена задача сопряжения двумерной модели аорты с нульмерной моделью кровеносной системы. Такой подход позволил получить нестационарные профили скорости течения крови и провести сравнительный анализ результатов, полученных для одномерной и двумерной модели сосуда.
Рассмотрим модель сердечно-сосудистой системы человека, состоящую из 55 крупных артерий (рис.24) [9]. Каждый из сосудов, входящих в модель, может быть представлен одномерной моделью гемодинамики. В этом случае течение в кровеносном сосуде (рис.2) описывается системой нелинейных гиперболических уравнений (1)-(3).
В модели артериального дерева [9], используемой в диссертационной работе, в качестве зависимости (2) принята функция вида:
Здесь И0 и s0 - толщина стенки сосуда и площадь поперечного сечения сосуда в невозмущенном состоянии при р = р0, р0- невозмущенное давление, ра1 внешнее давление, Е - модуль Юнга, v - коэффициент Пуассона 0 = -, что соответствует несжимаемой стенке). Константы /?, .s0 - являются характеристиками сосуда, их значения для каждой артерии приведены в таблице 1[9].
Граничные условия для системы (1)-(3) моделируют взаимное влияние сосудов и работу сердца. Все сосуды, входящие в рассматриваемую модель можно классифицировать следующим образом:
1. Восходящая аорта. На левой границе задан поток крови, выбрасываемый сердцем в кровеносную систему. На правой границе аорта разделяется на 2 «дочерних» сосуда (бифуркация);
2. Сосуды, выходящие из «родительского» сосуда и разделяющиеся на правой границе на два «дочерних» (бифуркация);
3. Терминальные артерии - сосуды, не имеющие «дочерних» сосудов.
На бифуркациях в качестве граничных условий используются равенства расходов и давлений в "родительском" и "дочерних" сосудах: U.S. =U,S, +U-.S, . . 2 2 зз (?6) Рі=Рі=Рз Здесь индекс 1 соответствует «родительскому сосуду», 2 и 3 - «дочерним» сосудам. Для терминальных артерий на правой границе могут быть реализованы различные виды граничных условий. В настоящей работе используются как неотражающие граничные условия [см. напр., 75]:
Моделирование пережатия плечевой артерии в артериальном дереве
Поскольку, как правило, манжета для измерения давления накладывается на левую руку, исследовалось пережатие левой плечевой артерии. Был проведен следующий вычислительный эксперимент. Сначала расчет выполнялся в условиях отсутствия манжеты, то есть рех1 = 0. Затем давление в манжете нагнеталось до 10 мм рт.ст. со скоростью 5 мм рт.стУс, после чего оно поддерживалось постоянным в течение времени, необходимого для получения установившегося решения (как показывают расчеты, для этого достаточно 3 секунд). Далее процесс нагнетания давления пережатия повторялся, т.е. давление нагнеталось до 20 мм рт.ст. с той же скоростью, поддерживалось постоянным и т.д, вплоть до достижения давления пережатия 120 мм рт.ст. При этом для остальных артерий рех, оставалось равным нулю в ходе всего процесса.
Общая схема вычислительного эксперимента проиллюстрирована на рис.32, где показаны изменение со временем давления пережатия (штриховая линия, правая шкала) и пульсаций среднеинтегралыюй площади поперечного сечения сосуда (сплошная кривая, левая шкала).
На рис. 33-38 представлены временные распределения гемодинамических параметров в центре левой плечевой артерии в невозмущенных условиях pai = 0 мм рт.ст. и при ее пережатии ра1 = 80 мм рт.ст. Можно отметить увеличение скорости течения крови в пережатом сосуде, а также уменьшение расхода крови, связанное с уменьшением площади поперечного сечения сосуда вследствие его пережатия. Кроме того, на всех распределениях гемодинамических параметров, полученных в условиях пережатия присутствует небольшой фазовый сдвиг по сравнению с аналогами для невозмущенных условий. На рис.37 представлены распределения расхода крови для артерии, являюшейся «родительской» по отношению к левой плечевой артерии - левой подключичной. На графике видно изменение амплитуды колебаний. На рис. 38 приведено распределение расхода крови для «дочерней» по отношению к исследуемой, левой лучевой артерии. Можно отметить наличие небольшого фазового сдвига, однако амплитуда колебаний практически не изменилась.
На рис. 39 и 40 представлены пространственные картины течения, наблюдаемые на периоде в невозмущенных условиях и при давлении пережатия ра1 = 80 мм рт.ст. В момент времени t = 0 (здесь и далее момент времени t = 0 с соответствует аналогичному моменту на рис. 25.) процесс находится в диастолической стадии. При / = 0.17 с возмущения, связанные с распространением пульсовой волны, достигают левой плечевой артерии. В последующие моменты времени можно видеть распространение возмущений, связанных с прохождением пульсовой полны через артерию (/ = 0.2с, / = 0.25с). В момент времени / = 0.3 с отраженная от бифуркации волна движется в обратном направлении. При / = 0.6 с стадия распространения возмущений в сосуде заканчивается. Поскольку решение имеет периодический характер с периодом Г = 0.8 с, начиная с момента времени t = T, картина течения повторяется.
Сравнивая кривые на рис. 39 и 40, можно отметить возрастание крутизны профилей при пережатии артерии (рис.40) по сравнению с невозмущенным состоянием (рис.39). Это обстоятельство свидетельствует об актуальности использования монотонизированных схем повышенного порядка точности для решения задач рассматриваемого класса.
Представляет интерес сопоставление результатов численного моделирования с эффектами, наблюдаемыми (и используемыми) в ходе измерения давления аускультативным и осциллометрическим методами. Они хорошо известны и подробно описаны в литературе (см. напр. [87]).
В основе наиболее точного, аускультативного, метода лежит появление слышимых тонов (тонов Короткова) при пережатии сосуда манжетой. С целью исследования этого феномена в настоящей диссертационной работе был проведен спектральный анализ колебаний площади поперечного сечения на правой границе левой плечевой артерии. Это позволило получить коэффициенты изменения амплитуды колебаний для различных частот в зависимости от давления пережатия. Результаты представлены на рис.41. Здесь все амплитуды отнесены к соответствующим значениям для невозмущенного сосуда. Кривая получена для частоты 1.25 Гц, кривая 2 - для частоты 12.5 Гц, 3-25 Гц. Особый интерес представляет поведение кривой 3, соответствующей частоте, лежащей в слышимом человеческим ухом диапазоне. Видно значительное повышение амплитуды при достижении давления пережатия, равного диастолическому, и резкое падение амплитуды при значениях пережимающего давления, превышающих величину систолического давления. Это подтверждает известную гипотезу о том, что появление тонов Короткова связано с нелинейным характером поведения упругих свойств оболочки сосуда. Отметим, что ранее подобный вывод был сделан в работе [88] на основе численных результатов, полученных с использованием нульмерной модели гемодинамики.
Другая группа измерительных методов (осциллометрические методы) основана на обработке осциллограмм, регистрирующих колебания давления в манжете при пережатии артерии. Поскольку они обусловлены колебаниями стенки сосуда, в данной диссертационной работе был проведен анализ изменения амплитуды колебаний средней по длине плечевой артерии площади поперечного сечения при изменении пережимающего давления рех1.
Полученные результаты представлены па рис. 42. На графике можно видеть существенное усиление колебаний, когда давление пережатия близко к диастолическому, дальнейшее нарастание амплитуды и последующее ослабление колебаний по мере приближения ра1 к систолическому давлению.
