автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное моделирование кровотока при наличии сосудистых имплантатов или патологий

На правах рукописи

Добросердова Татьяна Константиновна

Численное моделирование кровотока при наличии сосудистых имплантатов или патологий

13.18 - Математическое моделирование, численные методы и комплексы

программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Москва - 2013

1 7 ОКТ 2013

005535010

Работа выполнена в Федеральном государственном образовательном учреждении высшего профессионального образования "Московском государственном университете имени М.В.Ломоносова".

Научный руководитель:

доктор физико-математических наук, доцент Ольшанский Максим Александрович Официальные оппоненты: Нечепуренко Юрий Михайлович,

доктор физико-математических наук, доцент, Федеральное государственное бюджетное учреждение науки Институт вычислительной математики Российской академии наук, ведущий научный сотрудник Чупашн Александр Павлович,

доктор физико-математических наук, профессор, Федеральное государственное бюджетное учреждение науки Институт гидродинамики им. М.А. Лаврентьева Сибирского отделения Российской академии наук, заведующий лабораторией дифференциальных уравнений Ведущая организация:

Федеральное государственное образовательное учреждение высшего профессионального образования "Московский физико-технический институт (государственный университет)"

Защита состоится 31 октября 2013 г. в 14 часов на заседании диссертационного совета Д 002.045.01 при Федеральном государственном бюджетном учреждении науки Институте вычислительной математики Российской академии наук (ИВМ РАН), расположенном по адресу: 119333, г. Москва, ул. Губкина, 8.

С диссертацией можно ознакомиться в библиотеке ИВМ РАН. Автореферат разослан 28 сентября 2013 г.

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

диссертационного совета Д 002.045.01,

доктор физико-математических наук

Бочаров Г. А.

Общая характеристика работы

Актуальность темы исследования. По данным федеральной службы государственной статистики РФ за последние десять лет по причине болезней системы кровообращения (БСК) ежегодно происходит более 55% всех смертей в РФ. Эта причина лидирует во всем мире. В настоящее время значительно увеличивается количество стентированнй, ангиопластик, число применений современных технологии диагностики и лечения острой сосудистой патологии. В проекте государственной программы по разяитию здравоохранения в РФ до 2020 г. планируется дальнейшее усиление мер по оказанию специализированной медицинской помощи больным с острой сосудистой патологией. Данная тема привлекает внимание ученых и врачей во всем мире.

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

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

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

3

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

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

Исследования, вошедшие в диссертацию, были частично поддержаны грантами РФФИ 10-01-91055-НЦНИ_а, 11-01-00767-а, 11-01-00971-а, 12-01-00283-а, 12-01-33084-а, федеральной целевой программой "Научные и научно-педагогические кадры инновационной России".

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

Для достижения заданной цели требовалось разработать две стратегии. В первом случае для реализации поставленной задачи используется волоконная или пружинно-волоконная модель эластичной стенкн сосуда. С ее помощью можно воспроизвести зависимость трансмурального давления от площади поперечного сечения сосуда. Эта зависимость далее используется в одномерной модели глобального кровообращения. В рамках данного подхода необходимо было создать цельную технологию, позволяющую импортировать данные между моделями.

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

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

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

В работе предложены новые граничные условия на стыке областей разных размерностей при синтезе одномерной и трехмерной моделей течения жидкости. Эти условия обеспечили выполнение энергетического баланса для двумасштабной модели течения жидкости при использовании конвективной формы записи уравнений Навье-Стокса; выведена энергетическая оценка. Для вычислений предложена и численно исследована схема расщепления второго порядка точности. Подобные двухмасштабные модели могут использоваться для изучения кровотока в сети сосудов с произвольной геометрией как имплантатов, так и сосудов, их стенок, патологий. Никогда ранее не производился расчет кровотока многомасштабной моделью с сильно анизотропными включениями в трехмерной области, такими, как кава-фильтр. В диссертационной работе приведены результаты указанного численного эксперимента. Для решения уравнений Навье-Стокса в описанной двухмасштабной модели впервые используется метод бисопряженных градиентов с блочным псрсобуславлива-телем специального вида.

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

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

Разработанные модели являются этапом в создании "Виртуальной системы кровообращения". Проект реализуется в ИВМ РАН совместно с Василевским Ю.В., Симаковым С.С., Ивановым Ю.А., Крамаренко В.К. Данная технология, оснащенная удобным графическим интерфейсом на сенсорной панели с функциями мультитач, позволит делать гемодинамические расчеты для конкретного пациента (с использованием данных магнитно-резонансной или компьютерной томографии), учитывая влияние различных факторов, в частности, патологий н имплантатов. Разработка может использоваться для повышения точности прогнозирования хирургических операций. Положения, выносимые на защиту:

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

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

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

Степень достоверности и апробация результатов. Основные результаты диссертации докладывались на следующих семинарах и конференци-

б

ях: науный семинар института прикладной математики имени М.В. Келдыша РАН (ИПМ РАН. 2013 г.); семинар '''Mathematical modeling of natural disasters and technical hazards" (г.Сьон, Швейцария, 2013г.); международная конференция по математической теории управления и механике (г. Суздаль. Россия, 5-9 июля 2013 г.); научный круглый стол "Современные проблемы и инновационные перспективы моделирования кровообращения" (ФЦ сердца, крови и эндокринологии им. В.А. Алмазова, г. Санкт-Петербург, 24 июня 2013 г.); семинар кафедры вычислительной математики под руководством проф., д.ф.-м.н. A.C. Холодова (МФТИ, 2013 г.); научный семинар Института вычислительной математики РАН (ИВМ РАН. 2013 г.); семинар кафедры вычислительной математики под руководством проф., д.ф.-м.н. Г.М. Кобелькова (мех-мат МГУ, 2013 г.); научный круглый стол "Cardiovascular simulations: challenges and perspectives" (Университет Хьюстона, США, 29 апреля 2013 г.); день математического моделирования "Инновации в фармацевтике и медицине" (Москва, Россия. 14 ноября 2012 г.); VI Всероссийская конференция "Актуальные проблемы прикладной математики и механики", посвященная памяти академика А.Ф.Сидорова (Абрау-Дюрсо. Россия, 2012); 4th Workshop on Advanced Numerical Methods for Partial Differential Equation Analysis (Санкт-Петербург, Россия. 22 - 24 августа 2011); I (2010 г.), II (2011 г.), III (2011 г.), IV (2012 г.) конференции по математическим моделям и численным методам в биоматематике (ИВМ РАН, Москва, Россия); 52-я (2009 г.), 53-я (2011 г.), 55-я (2012 г.) научные конференции МФТИ "Современные проблемы фундаментальных и прикладных наук" (МФТИ, Россия, Московская область, г.Долгопрудный); Лобачевские чтения - 2009 (Казань, Россия).

Публикации. Материалы диссертации опубликованы в 12 печатных работах, из них 5 статей в журналах из Перечня ведущих рецензируемых научных журналов и изданий, рекомендованных ВАК. [1-5] и 7 тезисов докладов.

Личный вклад автора. Содержание диссертации и основные положения, выносимые на защит}-, отражают персональный вклад автора в опубликованные работы. Подготовка к публикации полученных результатов проводилась совместно с соавторами в работах [1-5]. Вклад соавторов равновелик. Диссертантом были реализованы составные модели течения крови в сети сосудов с патологиями или имплантатами на основе имеющихся моделей глобаль-

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

Структура и объем диссертации. Диссертация состоит из введения, трех глав, заключения и библиографии. Общий объем диссертации 102 страницы, включая 16 рисунков и 13 таблиц. Библиография включает 96 наименований на 12 страницах.

Содержание работы

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

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

В разделе 1.1 описывается используемая модель глобального кровообращения 1. В пупктс 1.1.1 приводится математическая постановка задачи. В основе модели лежат законы сохранения массы и импульса для вязкой несжи-

1 Холодов A.C., Симаков С.С. Численное исследование содержания кислорода в крови человека при низкочастотных воздействиях // Математическое моделирование. 2008. Т. 20(4). С. 87-102.

маемой жидкости:

ds d(Su)

дй t д(й2/2 + p/p) т+-al— = l>{t,x,s,u),

х € [0,1} — координата вдоль сосуда, длина которого I: í Е [0,Т] — время и интервале от начала расчетов до момента Т; S сечение сосуда; й — скорость, осредненная по сечению; р — трансмуральное давление; р — плотность крови; !.р, ф — заданные функции. В данной работе положим ср = 0, а ф будет задавать вязкое трение и определяться формулой:

ф = -16^шг/(5)/(&г2), (2)

2, Б > 1,

где S = SS~l, 77(5) „ „ , .

S — площадь поперечного сечения сосуда при р = 0; d — диаметр трубки; ц — коэффициент вязкости крови. Третьим уравнением модели является уравнение состояния, представляющее собой зависимость трансмурального давления от площади поперечного сечения. Для здоровых сосудов принимается эмпирическая функция

Р = pclf(S), (3)

Гехр(^-1)-1, П1 \\n(SS~i), S<S' [)

где со — скорость распространения магтых возмущений. Уравнение состояния является главной характеристикой эластичных свойств стенок сосудов в модели. Далее приводятся краевые условия для задачи (1) в точках стыковки сосудов. В разделе 1.1.2 обосновывается диссипативиость модели глобального кровообращения, выводится энергетическое равенство. В разделе 1.1.3 рассматриваются численные методы для гемодннамнчсскпх расчетов. Приводится двухслойная консервативная разностная схема для системы уравнений (1). Граничные условия дискретизируются с первым порядком точности, полученная нелинейная система рассчитывается методом Ныотопа. В разделе 1.1.4 рассматриваются допустимые дополнительные предположения и обобщения

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

В разделе 1.2 описывается один из таких методов. Суть его заключается в модификации уравнения состояния (3). В пункте 1.2.1 рассматриваются особенности этой функции, в частности, нелинейность и 5-образность ее графика. В пункте 1.2.2 приводится краткое описание имеющихся волоконной модели эластичной стенки сосуда и пружинно-волоконной модели стенки артерии с атеросклеротической бляшкой. Указанные модели позволяют рассчитывать реакцию стсики сосуда на деформацию. Полученная зависимость транс-мурального давления от площади поперечного сечения далее используется в модели глобального кровотока. В пункте 1.2.3 описано несколько численных экспериментов по моделированию течения крови в сети сосудов с атеросклерозом. Предполагается наличие как одной, так и нескольких атеросклероти-чсских бляшек разных размеров, сравниваются результаты расчетов в различных частях организма. Полученные данные подтверждают возможность проявления хорошо известных в медицине осложнений и симптомов данного заболевания.

Результаты первой главы опубликованы в работах [2-5].

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

В разделе 2.1 представлена модель течения жидкости в трехмерной области Г2зп, основанная на уравнениях Навье-Стокса. Здесь приводится математическая постановка задачи. Далее выводится энергетическая оценка для модели. Полученное равенство показывает, что при однородных граничных условиях наблюдается диссипация энергии. Уравнения Навье-Стокса днекре-тизируются по времени со вторым порядком точности. В зависимости от того, как аппроксимируется конвективный член, система дифференциальных уравнений оказывается линейной или нелинейной.

В разделе 2.2 описываются численные методы для решения полученной задачи. Если задача линеаризована, дальнейшая дискретизация проводится

А

1 out

VL

odown "in

ID

x=b

Рис. I. Схема стыковки П30 и областей.

методом конечных элементов. Полученная система уравнений решается методом бисопряженных градиентов с блочным треугольным переобуславливате-лем специального вида 2. Используемый переобуславливатель эффективен для достаточно большого диапазона параметра вязкости., а также его собственные значения могут быть ограничены константами, не зависящими от шага расчетной сетки. Если решаются нелинейные уравнения Навье-Стокса, используется метод Ньютона-Крылова. Все используемые численные методы реализованы с помощью программного пакета АшЗБ.

Раздел 2.3 посвящен сопряжению трехмерной модели течения жидкости с одномерной моделью глобальной циркуляции крови. Рассматривается трехмерная область Г2зо (Гои1; - граница вытекания) и присоединенная к ней одномерная П^™, граничащая с П3о в точке с1 (см. рис.. 1). В разделе 2.3.1 приводится обзор существующих граничных условий на стыке областей разных размерностей, требующих непрерывности каких-либо осредненных величин: нормальной компоненты скорости, осредненной по сечению, площади поперечного сечения, давления, нормальной компоненты тензора напряжений, входящей характеристики или потока. Наиболее популярной комбинацией являются требования непрерывности нормальной компоненты тензора напряжения и потока:

2 Olshanskii М.А.. Vassilevski Yu.V. Pressure Schur complement precotiditioners for the discrete Oseen problem// SIAM J.Sci.Comp. 2007. V.29. P.2686-2704

u • n Iis = Su.\x=d,

(6)

r,

где u — трехмерная скорость; р — давление в П3о; p,S,ü — давление, площадь поперечного сечения н скорость, осредненная по сечению, в П^р™; ¡у — вязкость жидкости; п — внешняя нормаль к поверхности. Однако, с условиями (5), (6) двухмасштабная модель не является диссипативной при однородных условиях на границе втекания в трехмерную область Г,п и на границе вытекания из одномерной области при х = е. Сохранение данного свойства желательно с точки зрения физического смысла, поскольку одномерная и трехмерная модели по отдельности являются дисспативными.

В работе Formaggia и др. предложено использовать условие (б) совместно с требованием непрерывности полного напряжения:

~ + + 2|u|2)n = {p+2ü2)\x=än "а Гои" (?) В этом случае двухмасштабная модель оказывается диссипативной, однако, условие (7) естественно для вихревой формы уравнений Навье-Стокса. Численный расчет таких уравнений затруднителен средствами существующих программных пакетов.

В диссертационной работе решаются уравнения Навье-Стокса в обычной конвективной форме, естественным граничным условием для которых является (5). Считая полное напряжение постоянным вдоль границы rout, как следует из (7), получим некорректное условие вытекания для простейшего течения Пуазейля. Для сохранения дпссипативностп двухмасштабной 1D-3D модели предлагается совместно с условием (5) использовать новое граничное условие, требующее непрерывности линейной комбинации потоков энергии и жидкости:

р | u-nds + l J |u|2(u • n)ds = {pSü + (8)

Гout Гoui

Если обозначить энергию одномерной модели для сосуда Clf™" через £w(t), а трехмерной модели для области Г230 через £з£>(0: верна следующая теорема об энергетическом балансе для двухмасштабной модели:

Теорема 1. Рассмотрим задачу течения жидкости в области Г2зп - Пщ™, основанную на уравнениях Навье-Стокса и системе уравнений (1), (3), (4), с граничными условиями (5), (8) на стыке областей разных размерностей. Пусть на остальных границах заданы однородные граничные условия: u;n =

12

0; ы|.г=е = 0. Тогда достаточно гладкое решение удовлетворяет следующему энергетическому равенству:

т т

£зо(Т) + £ю(Т) + и ||Уи||2с1Ь +

К^йЧх^. = £,о(0) + £1О(0) (9) о о п-д™

для любого Т > 0. Еаьи ф определена как (2), тогда К„(3) = 1б1/г)(3)3с1~2 > О, и совместная энергия модели монотонно убывает по вре.,иени для любого ненулевого решения:

Если к трехмерной области Пзо но границе Г;„ присоединить одномерный сосуд стыкующийся в точке 6 (см. рис. 1), возможным граничным условием может быть условие свободного стока в точке Ь н комбинации с требованием непрерывности нормальной компоненты скорости, осреднепной по сечению, или потока, или линейной комбинации потоков жидкости и энергии (аналогично условию (8)) на стыке областей Пзо и

Для расчета течения жидкости описанной двухмасштабной моделью в области (см. рис. 1) предложена схема расщепления, описанная в разделе 2.3.2 - 2.3.3. Обозначим через ¿¿".р™, ип и рп значения соответствующих параметров в момент времени Ь = £,г. Используя эти величины, вычислим йп+1,рп+1, 5П+1, и"+1 и рп+1 при £ = ¿„+1 (Д£ = £„+, - £„) в три этапа:

Шаг 1. Проинтегрируем (1) на интервале при £ 6 [£„, £п+1] с данной и(£п+1) в точке х = а и условием свободного стока в точке х = Ь.

Шаг 2. Из условия сопряжения моделей в точке а, вычислим и^ с использованием

{й,р,5} = {йп+1,Р™+\£'"+1}-Найдем р* и 5* как линейную экстраполяцию р\хы и •?|.-с=<г с временных слоев £п и £п_1 на слой Решим задачу Навье-Стокса в Пзо относительно ип+1,р"+1 с условием

Зип+1

на Г,

Шаг 3. Найдем 1£п+1|х^ из уравнения

(P'S^Tln+l + ^S'(йn+lf)\x=d = f | и"+1-пс1з + ^ | К+1|2(и"+1-п)с1з.

ГоиС Г'оиб

(10)

(Для нахождения также можно использовать условие непрерыв-

ности нормальной компоненты скорости, осредненной по сечению, или потока.) Теперь, используя й1111 для граничных условий в точке х = (1 и условие свободного стока в точке х = е, интегрируем (1) при £ 6 и находим йп+\рп+\, 5П+1 в

Данный алгоритм можно использовать с разными шагами по времени в одномерной и трехмерной областях.

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

В разделе 3.1 приведен первый эксперимент. Рассчитывается течение жидкости в области П^-Пзо-^т™ (см-

рис. 1), где трехмерная область Пзо является цилиндром, а аналитическое решение задачи известно. Сначала уравнения Навье-Стокса линеаризуются, дискретизируются методом конечных элементов, исследуется сходимость метода бисопряжениых градиентов с блочным пе-реобуславливаетелем специального вида. Установлено, что сходимость незначительно зависит от вязкости и улучшается при измельчении сетки. Надежность алгоритма при различных значениях параметра вязкости обнаружена только для достаточно малых значений шага по времени (в рассматриваемых случаях не больше 0.01). Наблюдается сходимость численного решения к аналитическому решению со вторым порядком точности.

Тот же эксперимент рассчитывается с использованием нелинейных уравнений Навье-Стокса. Метод Ньютона-Крылова, применяемый для их решения, оказывается надежным примерно в том же диапазоне параметров, имеет также второй порядок сходимости и такую же величину погрешности. Однако, время, затрачиваемое на расчет нелинейных уравнений, оказывается значительно большим, чем затрачиваемое на. расчет линейных, из-за большего количества линейных итераций на каждом шаге по времени. По этой причине в последующих расчетах уравнения Навье-Стокса линеаризуются.

14

В разделе 3.2 описана вторая серия экспериментов. Схема расщепления тестируется на задачах с известным аналитическим решением н использованием различных условий сопряжения решений на стыке областей разных размерностей. Во всех случаях установлен второй порядок сходимости.

В разделе 3.3 тестируется вычисление важных прикладных характеристик. таких как сила сопротивления и перепад давления, создаваемый жидкостью при обтекании препятствия. Также тестируется применение новых условий сопряжения 1D и 3D моделей (5), (8). Следуя работам Schäfer, Turek, Braack, в качестве трехмерной области рассматривается канал с прямоугольным сечением, в качестве препятствия — круговой цилиндр. Течение задается стационарным с числом Рейнольдса 20 или нестационарным с максимальным числом Рсйнольдса 100. Для указанных задач известен диапазон значений исследуемых характеристик. В описанных тестах не устанавливалось ограничений на условия вытекания, таким образом, мы использовали условия (5), (8), присоединив к трехмерному каналу одномерный сосуд.

Численные результаты хорошо приближают эталонные. Используемый в Пзо метод конечных элементов для задач с условиями сопряжения моделей (5), (8) на границе вытекания rout обеспечивает надежную и устойчивую аппроксимацию решения.

В разделе 3.4 описан расчет двухмасштабной моделью течения крови в вене с установленным кава-фильтром. В качестве трехмерной области брался участок нижней полой вены длиной 4.5с.м с эллиптическим сечением 1.6 х 2.4см. Кава-фильтр прикреплен на расстоянии 0.5см от границы втекания, его длина 2см, диаметр каждой из 12 ножек 0.5мм. Для такой области строилась адаптивная сетка. Кровь считается вязкой несжимаемой жидкостью с динамической вязкостью 0.0055 Пахе и плотностью 1 г/см3. Профиль скорости на границе втекания в П"р строился по данным доплерографии в нижней полой вене п был аппрокспмирован гладкой периодической функцией.

Результаты расчетов показали, что за кава-фильтром наблюдается замедление кровотока и обратные течения. Это объясняет повторное тромбообразо-вание на ножках капа-фильтра, являющееся побочным эффектом установки имплантата. Численное решение в этой области не осеспмметрично: возмуще-

ния, вызванные несимметричностью тетраэдральной сетки, достаточно значимы для появления неустойчивости в виде дорожки Кармана за препятствием.

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

Результаты третьей главы опубликованы в работе [1].

В Заключении кратко сформулированы результаты диссертационной работы.

Список публикаций

1. Dobroserdova Т. К., Olshanskii М. A. A finite element solver and energy stable coupling for 3D and ID fluid models // Computer Methods in Applied Mechanics and Engineering. 2013. V. 259. P. 166 - 176.

2. Vassilevski Yu., Simakov S., Salamatova V., Ivanov Yu., Dobroserdova T. Numerical issues of modelling blood flow in networks of vessels with pathologies / / R ussian Journal of Numerical Analysis and Mathematical Modelling. 2011. V. 26(6). P. 605-622.

3. Vassilevski Yu., Simakov S., Salamatova V., Ivanov Yu., Dobroserdova T. Vessel wall models for simulation of atherosclerotic vascular networks // Mathematical Modelling of Natural Phenomena. 2011. V. 6(7). P. 82-99.

4. Vassilevski Yu.. Simakov S., Salamatova V.; Ivanov Yu.; Dobroserdova T. Blood flow simulation in atherosclerotic vascular network using fiber-spring representation of diseased wall // Mathematical Modelling of Natural Phenomena. 2011. V. 6(5). P. 333-349.

5. Иванов Ю. А., Добросердова Т. К. Математическое моделирование влияния установки кава-фильтра на гемодинамику кровеносной системы / / Научно-технический вестник СПбГУ ИТМО. 2010. Т. 04(68). С. 94-98.

16

Подписано в печать 27.09.2013 Формат 60x88 1/16. Объем 1.0 п.л. Тираж 100 экз. Заказ № 1333 Отпечатано в ООО «Соцветие красок» 119991 г.Москва, Ленинские горы, д.1 Главное здание МГУ, к. А-102

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ИМЕНИ М.В. ЛОМОНОСОВА

На правах рукописи

04201362776

Добросердова Татьяна Константиновна

Численное моделирование кровотока при наличии сосудистых имплантатов или патологий

05.13.18 - Математическое моделирование, численные методы и комплексы

программ

ДИССЕРТАЦИЯ на соискание ученой степени кандидата физико-математических наук

Научный руководитель д. ф.-м. н.

Ольшанский Максим Александрович

Москва - 2013

Содержание

Введение ........................................................................4

Глава 1. Моделирование влияния патологий на кровоток посредством изменения упругой модели для стенок сосудов............20

1.1. Модель глобального кровообращения................................20

1.2. Учет патологий и имплантатов моделью глобального кровообращения ....................................................................35

1.3. Выводы..................................................................47

Глава 2. Сопряжение 10 модели глобального кровотока и ЗО модели течения жидкости в канале......................................49

2.1. Трехмерная модель течения несжимаемой вязкой ньютоновской жидкости................................................................49

2.2. Численное решение уравнений Навье-Стокса........................52

2.3. Сопряжение одномерной и трехмерной моделей течения жидкости для моделирования кровотока....................................58

2.4. Выводы..................................................................70

Глава 3. Численные эксперименты......................................71

3.1. Сравнение расчетов 10-30-1Б задачи с использованием линеаризованного и нелинейного уравнений Навье-Стокса..................71

3.2. Тестирование схемы расщепления с различными условиями сопряжения моделей......................................................76

3.3. Моделирование обтекания кругового цилиндра ....................79

3.4. Моделирование обтекания кава-фильтра............................83

3.5. Выводы..................................................................88

Заключение......................................................................90

Литература

Введение

Главной причиной смертности в мире являются сердечно-сосудистые заболевания. Атеросклероз — самое распространенное среди них. В результате болезни часто поражаются сразу несколько артерий, поэтому влияние и развитие патологического процесса необходимо рассматривать в сети сосудов. Для устранения стенозов производят стентирование артерий.

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

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

Современные представления о строении и функционировании сердечнососудистой системы были заложены революционными открытиями Вильяма Гарвея (1578-1657гг.) [2]. Он обосновал различие дыхания и пульса, подробно описал работу сердца, физиологию сосудов. Гарвеем были установлены путь движения крови, замкнутость и наличие двух кругов в системе кровообращения. Дальнейшее развитие в этой области обязано не только эмпирическим наблюдениям и догадкам врачей, физиологов и анатомов, но и быстро развивающимся точным наукам, предоставляющим информацию о свойствах газов, жидкостей, тканей, а также оборудование и методы исследований. В средние

века большинство ученых были специалистами сразу в нескольких областях, в частности, многие врачи были физиками. Именно поэтому некоторые физиологические явления пытались описать с точки зрения законов механики и математики. Так, например, фундаментальные уравнения гидродинамики — уравнения Эйлера и Бернулли, описывающие движение идеальной жидкости, впервые предназначались для описания движения крови по сети сосудов.

Множество открытий привело к усложнению наук и специализации ученых: медициной стали заниматься исключительно врачи, а детальным исследованием природы различных явлений — физики и инженеры. За последние десятилетия взаимодействие физики, математики и медицины приняло другую форму. На стыке этих наук возникла новая область исследований, а именно математическое моделирование физических процессов, в нашем случае кровообращения.

Модель аортальной компрессионной камеры (АКК), впервые предложенная Стивеном Хейлзом в 1733 г. [47], а дальше активно исследуемая и развиваемая Отто Франком [44], была одной из первых моделей для описания гемодинамики. В ее основу легла главная функция аорты и крупных артерий — трансформация дискретно поступающего в аорту сердечного выброса в непрерывный, несколько пульсирующий поток артериолярно-капиллярного русла. Подобная трансформация происходит благодаря упругим свойствам стенок. В системе кровообращения сердце можно рассматривать как генератор потока крови, АКК как источник давления, обеспечивающий непрерывный кровоток в капиллярном и, частично, в венозном русле. Движение крови в сосудах — сочетание работы этих двух генераторов. Внутренним сопротивлением АКК может рассматриваться входной импеданс артериальной системы, являющийся внешним сопротивлением для левого желудочка. Внешней же нагрузкой АКК является периферическое сопротивление артериальной системы, величина которого может изменяться в соответствии с текущими запросами на кровоснабжение тканей. Модуль артериальной эластичности характеризует взаимосвязь между

изменениями общего объема артериальной системы (Д V) и соответствующими изменениями давления (АР) заполняющей ее крови. В общем виде этот показатель может быть представлен соотношением:

п Ар

£ = дг

В ходе экспериментальных исследований (Соре, 1965), проводившихся на препаратах аорты и крупных артерий, была установлена высокая степень постоянства этого показателя в нормальных диапазонах изменения давления. На предположении относительного постоянства величины Е основывается теория АКК (Лайтфут, 1977), фундаментальное уравнение которой имеет следующий вид [14]:

н

где Р(Ь) - давление в АКК; Е - модуль артериальной эластичности; (¿(Ь) - входной кровоток артериальной системы; Я - периферическое сопротивление. С помощью этого уравнения можно получать различные зависимости от времени для давления Р(£) в артериальной системе.

В своих работах Уомерсли (1907—1958 гг.) рассматривал сосуд как однородную твердую трубку и использовал в ней линеаризованное уравнение Навье-Стокса [91]. Такой подход не позволял воспроизвести многие эффекты, например, отражение волн от областей бифуркации и т.п. Постепенно были приняты во внимание упругие свойства стенок сосудов [92, 95], с помощью Фурье-анализа объяснены нелинейные эффекты [93, 94], рассмотрены случаи с большой и маленькой вязкостью [58].

В настоящий момент используются различные подходы к моделированию кровообращения, в том числе стохастические методы [49], методы теории управления [26, 77], комбинированно фрактально-вейвлетный анализ [57] и другие. Два самых больших класса составляют электромеханические [34, 36, 43, 56, 64, 67, 83] и гидродинамические модели [16, 24, 27, 33, 35, 60, 63, 73, 79, 82].

В моделях первого класса сердечно-сосудистая система аналогична некоторой цепи переменного тока. Сердцу сопоставляется элемент питания с переменной электро-движущей силой, а сосудистому руслу - набор сопротивлений и емкостей. Полное сопротивление сосудов может быть разделено на несколько частей: сопротивление больших артерий, артериол, капилляров, вен. Сопротивление артериол самое значимое. Оно также называется периферическим. Крупные артерии заменяются описанной выше АКК, позволяющей учесть эластичные свойства этих сосудов.Теория электрических цепей в настоящее время хорошо исследована, и данный подход помогает в решении многих задач. Однако, при его использовании не достаточно хорошо воспроизводятся гемоди-намические параметры течения крови с учетом свойств сосудов и локальных особенностей. Кроме того, не гарантируется однозначное соответствие между цепью и деревом кровеносных сосудов. Для каждой цепи придется заново выводить описывающие ее уравнения, и чем сильнее детализация, тем сложнее эти уравнения. Такую модель сложно адаптировать под конкретного пациента, тем более учесть влияние имплантатов или патологий на гемодинамику. Поэтому в данной работе будет использоваться гидродинамический подход.

В рамках гидродинамического подхода кровь считается вязкой несжимаемой жидкостью, протекающей по сети эластичных трубок. Для описания движения жидкости используются уравнения гидродинамики в каждом сосуде. В зависимости от требуемой детализации и сложности исследуемого процесса используются модели разных размерностей. Одномерные модели, см. например [16, 24, 60, 73, 79], позволяют делать гемодинамические расчеты во всей сосудистой сети. Первым этапом их построения является создание дерева сосудов — одномерного графа. В этой структуре определяется, какие сосуды соединяются между собой, координаты областей стыковок, их степени, при необходимости геометрия соединений, и т.п. Следующим этапом является задание параметров модели. Для этого требуется указать длины, диаметры сосудов, в некоторых случаях скорости распространений пульсовых волн, сопротивления некоторых

областей и т.д. Эти данные, как правило, берут из медицинских источников, где они получены опытным путем. Когда определена область интегрирования задачи, для каждого сосуда выписываются уравнения гидродинамики относительно осредненных по сечению сосуда величин, а затем решения сшиваются с помощью граничных условий в областях стыковки. Изменение сосудистого дерева и параметров каждого отдельного сосуда в таких моделях не требует изменения постановки математической задачи. Одна из таких моделей подробно описывается в разделе 1.1. Там же выводится фундаментальное соотношение, отражающее баланс энергии для данной модели (пункт 1.1.2), и исследованы численные свойства (пункт 1.1.3).

Простота модели позволяет учитывать множество физиологических процессов и внешних воздействий, например, влияние различных сил, нервную регуляцию, ауторегуляцию, работу мышц и других органов, дыхание, крово-потери, перенос веществ и т.п (см. раздел 1.1.4). Некоторые авторы пытаются также брать во внимание влияние имплантатов в сосудистой системе. В работе [10] воспроизводится кровоток при наличии кава-фильтра с тромбом в вене. Препятсвие учитывается посредством уменьшения эффективного сечения сосуда, однако, при этом не берется во внимание упругое воздействие имплантата на стенку сосуда. В работе [72] исследуется осесимметричная задача о течении крови вдоль артерии со стентом. Эластичные свойства стенки, полагаемой тонкостенной мембраной из нелинейного вязкоупругого материала, вычисляются через уравнение равновесия для такой оболочки. Описываемый метод приемлем только для осесимметричных задач, где толщиной стенки можно пренебречь. Таким образом, существующие методы учета влияния имплантатов на гемодинамику в модели глобального кровообращения, не использующие дополнительных численных моделей, имеют ряд ограничений и допущений. Указанные подходы не являются универсальными: их область применения включает только частные случаи возможных послеоперационных ситуаций.

Для более точного описания влияния патологий и имплантатов на гемоди-

намику в данной работе предложено два подхода:

1. синтез модели глобального кровообращения с волоконной моделью эластичной стенки сосуда;

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

В главе 1 изложен первый метод, для которого ключевыми являются эластичные свойства стенки сосуда. В модели глобальной циркуляции крови, рассматриваемой в разделе 1.1, для их описания используется уравнение состояния. Оно представляет собой зависимость трансмурального давления от попречного сечения сосуда и более подробно обсуждается в пункте 1.2.1. При установке стента в артерию, кава-фильтра в вену, появлении атеросклеротической бляшки или другой патологии эластичные свойства стенки сосуда меняются. Новое уравнение состояния рассчитывается с помощью волоконной или пружинно-волоконной модели эластичной стенки сосуда (см. пункт 1.2.2) и далее используется моделью глобального кровообращения. Математическая постановка задачи включает только обыкновенные дифференциальные уравнения, что обеспечивает вычислительно эффективную технологию. Данный метод позволяет изучать гемодинамику во всем организме при наличии как одной, так и нескольких патологий или имплантатов, а также не требует от них осесимметричности. На его основе проведен ряд численных экспериментов по исследованию влияния атеросклероза на кровоток, описанных в разделе 1.2.3.

Второй рассматриваемый подход к моделированию гемодинамики при наличии патологий и имплантатов является многомасштабным. Многомасштабные модели стали очень популярны в последнее десятилетие. В них комбинируются модели разных размерностей, например, (Ю (модели, использующие электромеханические аналогии) и ЗБ [78], Ш и 2Б или Ш и ЗБ [66, 74, 84]. Двумерные и трехмерные гидродинамические модели основываются на уравнениях Навье-Стокса (см. раздел 2.1.1). Вычисления с их помощью оказывают-

ся довольно трудоемкими, поэтому они используются для описания кровотока только в небольших окрестностях: в областях стыковки между сосудами (метод декомпозиции области) [65], в месте бифуркации артерий [84], в интересующем сосуде [78] или группе сосудов [46] и т.п. Описываемые 2D и 3D модели можно усложнить, дабавив уравнения движения упругой стенки сосуда. В этом случае получается задача о взаимодействии жидкости и твердой структуры, далее FSI (Fluid Structure Interaction problem) [66, 73]. Такая задача лучше отражает реальные физические процессы и позволяет не только рассматривать кровоток внутри сосудов, но и учитывать распространение пульсовых волн по их стенкам. Это дополнение является следующим возможным этапом развития модели, описанной в разделе 2.1.1.

Совместное решение 3D задачи FSI и уравнений течения жидкости в одномерном сосуде, присоединенном к трехмерной области на границе вытекания, описывается в работе [41], а также в [90] с использованием метода конечных элементов. Аналогичная модель представлена в статье [84]: в качестве трехмерной области берется бифуркация сонных артерий, а одномерный граф сосудов строится для всей артериальной части сосудистой системы. Некоторые авторы рассматривают кровоток в окрестности патологии, например, стеноза [30]. Однако, при усложнении области, появлении анизотропных элементов, например, кава-фильтра, расчеты значительно усложняются и доступны не любыми методами. Вариационную постановку задачи для таких многомасштабных моделей можно найти в работе [30]. Данные методы особенно эффективны, когда важно знать не только особенности локального кровотока, но и гемодинамику в целом. Одной из основных трудностей данного подхода становится задание корректных граничных условий на стыке областей разных размерностей.

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

дивателем специального вида для решения линеаризованных уравнений Навье-Стокса и метод Ньютона-Крылова для решения нелинейных уравнений Навье-Стокса (пункт 2.2). В разделе 2.3.1 для двухмасштабной ШЗЭ модели течения крови приведен обзор используемых граничных условий на стыке областей разных размерностей. Кроме того, предложены новые условия, накладывающие непрерывность линейной комбинации потоков энергии и жидкости. Их использование гарантирует выполнение энергетического баланса, позволяет использовать метод расщепления 2-го порядка точности для решения задачи, а также не требует в трехмерной области выполнения условий и-п<0ии-п>0на границах втекания и вытекания соответственно, где и — трехмерная скорость, а п — внешняя нормаль к поверхности. Это принципиально для расчетов, где возникают обратные течения, например, при моделировании кровотока в нижней полой вене. В пунктах 2.3.2 - 2.3.3 диссертации предлагается схема расщепления для двухмасштабной модели.

В третьей гл