CC BY
49
УДК 531/534:57
Я. А. Гатаулин, математик, А. Д. Юхнев, науч. сотрудник,
ФГБОУ ВПО «Санкт-Петербургский государственный политехнический университет» М. А. Попов, мл. науч. сотрудник,
Д. И. Курапеев, канд. мед. наук, зав. НИЛ биопротезирования и кардиопротекции, ФГБУ «ФЦСКЭ им. В. А. Алмазова» Минздрава России
Численное моделирование кровотока в обшей сонной артерии с S-образной извитостью
Ключевые слова: вычислительная гидродинамика, общая сонная артерия, S-образная извитость, сдвиговые напряжения, индекс колебаний сдвиговых напряжений, GIMIAS, ANSYS CFX.
Keywords: сomputational fluid dynamics, common carotid artery, S-shaped tortuosity, shear stresses, oscillatory shear index, GIMIAS, ANSYS CFX.
По данным компьютерной томографии разработана математическая модель кровотока в общей сонной артерии c S-образной извитостью. Для описания течения в общей сонной артерии с S-образной извитостью использована модель вязкой несжимаемой жидкости на основе уравнений Навье—Стокса. На входе в общую сонную артерию задается физиологическая кривая расхода, полученная методом фазово-кон-трастной магнитно-резонансной ангиографии. Рассматривается влияние плоского и параболического профилей скорости на входе в артерию на скорость течения и сдвиговые напряжения на стенке. Определены зоны с низкими значениями сдвиговых напряжений на стенке и высоким индексом колебаний сдвиговых напряжений, предрасположенные к развитию атеросклероза.
Введение
На сегодняшний день заболевания, связанные с нарушением мозгового кровообращения, вышли на ведущее место среди причин смерти и инвалиди-зации в развитых странах. По данным Всемирной организации здравоохранения с 2000 по 2010 год от церебрального инсульта умерли более 5 млн человек, а из 15 млн выживших более 80 % остались инвалидами. В Российской Федерации ежегодно регистрируется 0,4—0,5 млн новых случаев инсульта. Таким образом, проблема лечения пациентов, перенесших острое нарушение мозгового кровообращения, не только влечет за собой огромные финансовые потери, но и является острой социальной проблемой. Для ее решения в последнее время
активно используются методы математического моделирования физиологических процессов.
Разработка аппарата математического моделирования кровотока предполагает:
• построение замкнутой механико-математической модели процесса, описывающей поведение биологической среды на основе системы уравнений в частных производных механики сплошных сред;
• разработку замыкающих систему уравнений механики сплошных сред реологических соотношений, описывающих поведение той или иной среды (для гидродинамики это уравнения состояния, для механики деформируемого твердого тела — соотношения между компонентами тензоров напряжений и деформаций);
• корректную математическую постановку задачи: представление замкнутой системы уравнений механики сплошных сред, постановку необходимых для ее решения начальных и граничных условий, условий на контактных границах;
• разработку или реализацию вычислительных методов, адаптированных к специфике решения конкретной задачи;
• разработку алгоритма численного решения задачи и его программную реализацию;
• численное решение задачи и визуализацию полученных результатов [1].
Существует достаточно много численных и экспериментальных работ [2], посвященных изучению течения в бифуркации сонной артерии. Однако течения в общей сонной артерии с патологической извитостью изучены недостаточно.
Среди трудов отечественных исследователей следует отметить работу О. Е. Павловой и др. [3], в которой проведено численное исследование движения крови в анатомически реальной здоровой и патоло-
гически извитой бифуркации сонной артерии человека. Физиологические значения средней скорости кровотока в сечении бифуркации сонной артерии были получены с использованием неинвазивной технологии. 3D-модель была построена на основе данных компьютерной томографии. Предполагалось, что материал стенки сосудов линейно-упругий и изотропный. Сравнительный анализ полученных числовых данных здоровой и патологически извитой сонных артерий показал значительное влияние изгиба артерии на характер кровотока. Максимум давления смещается к выпуклой поверхности извитости, а максимум скорости смещается в противоположную сторону. Это приводит к образованию существенной асимметрии потока крови в поперечном сечении извитости. За счет указанного явления происходит заметное уменьшение сдвиговых напряжений на выпуклой стороне стенки. Было установлено, что наличие извитости дополнительно создает благоприятные условия для развития атеросклероза.
В численном исследовании влияния степени кривизны общей сонной артерии на скошенность профиля скорости [4] по данным магнитно-резонансной томографии построена реальная геометрия общей сонной артерии, на основе которой был создан набор идеализированных моделей сонной артерии, отличающихся радиусом кривизны. Расчеты показали, что модели общей сонной артерии с умеренной кривизной не демонстрируют развитого или осесимметрич-ного профиля скорости. Авторы делают вывод, что развитый профиль в общей сонной артерии — это «скорее исключение, чем правило» [4].
Цели нашего исследования — провести расчет кровотока в общей сонной артерии с S-образной извитостью с помощью численного моделирования и выявить зоны, предрасположенные к развитию атеросклероза.
Медицинские аспекты
поражения брахиоцефальных сосудов
По виду инсульты делятся на геморрагические (кровоизлияние) и ишемические. Поражения второго типа преобладают и составляют около 80 % от всех случаев острого нарушения мозгового кровообращения. Основными причинами ишемических инсультов являются атеросклеротическое поражение артерий и, реже, их патологические извитости. По современным представлениям, причина патологической извитости врожденная либо развивается при гипертонической болезни. Особенностью клинического течения патологических извитостей является их длительная бессимптомность, первым проявлением может стать острое нарушение мозгового кровообращения. Патологические извитости классифицируют по виду пространственного строения (рис. 1): изгибы (^ и S-образные), перегиб (кин-
а)
б)
в)
Рис. 1
Виды патологических извитостей сонных артерий: а — петля (койлинг); б — двойная петля; в — перегиб (кинкинг)
кинг) и петля (койлинг). Многолетние наблюдения показывают, что наиболее опасными видами патологических извитостей являются перегиб и петля.
Патофизиологические механизмы влияния сдвиговых напряжений на сосудистую стенку
В результате многочисленных исследований было замечено, что атеросклеротический процесс носит локальный характер, предпочтительно располагаясь в местах со сложной геометрией (бифуркациями, изгибами). Особенностью действия сдвиговых напряжений является то, что в процесс вовлекается только эндотелий и не оказывается воздействие на остальные слои сосудистой стенки, тогда как давление и растяжение влияют на все слои практически равномерно по всей толщине сосудистой стенки. В первых исследованиях предполагали, что высокие значения сдвиговых напряжений, вызывая повреждение эндотелия, способствуют развитию атеросклеротического процесса [5]. В настоящее время принята теория о проатерогенном значении низких величин сдвиговых напряжений [6]. В экспериментальных работах по изучению характера кровотока в коронарных артериях у пациентов, умерших от ишемической болезни сердца, с помощью визуализации течений полистироловыми микросферами [7] показано, что атеросклеротические бляшки локализуются на внешней части стенки бифуркации сосуда в местах со сложными вихревыми потоками, где имеют место срыв потока и, как следствие, низкие значения сдвиговых напряжений.
Сдвиговые напряжения оказывают деформирующее воздействие на эндотелиальные клетки, выстилающие внутреннюю стенку сосуда. Эндоте-
лиальные клетки артерий принимают участие в процессах регуляции тонуса кровеносных сосудов, работе свертывающей системы, выполняют барьерную функцию. Находясь на границе двух сред, сосудистый эндотелий опосредует действие биологически активных веществ крови, обеспечивая механочув-ствительность стенок сосудов, присущую как крупным, так и мелким артериям.
Эндотелиальные клетки, в том числе их ядра, ориентируются в направлении движения потока крови, и в эндотелии меняется скорость синтеза некоторых веществ. Изменения сдвиговых напряжений воспринимаются эндотелиальными клетками в виде деформаций их стенки. При этом эндотелий способен синтезировать и выделять высокоактивное сосудорасширяющее вещество — эндотелиальный релаксирующий фактор (ЭРФ), представляющий собой оксид азота. Кроме непосредственного действия на компоненты сосудистой стенки, окись азота влияет на активность форменных элементов крови, эффективно ингибирует как агрегацию, так и адгезию тромбоцитов и лейкоцитов к эндотелию сосудов [8].
Считается, что механорецепция опосредована через модуляцию селективных для К+- и Са+2-каналов при растяжении клеточной мембраны. Кроме того, предполагается, что в этом процессе участвует особая субпопуляция К+-каналов эндотелия (К+Са+2), чувствительных к осцилляциям потока [9].
Изучение влияния сдвиговых напряжений на атерогенез было затруднено в связи с отсутствием методик их объективной оценки. Начиная с 1990-х годов, развитие вычислительной техники и математического аппарата открыло возможность применения методов компьютерного моделирования.
Интересной особенностью изучения данной темы является то, что она требует тесного взаимодействия врачей со специалистами технического профиля, изучающими вопросы гидродинамики. В настоящее время считается, что физиологические значения сдвиговых напряжений для артерий находятся в диапазоне от 1 до 7 Па, значения меньше 1 Па способствуют развитию атеросклероза, больше 7 Па — могут повреждать эндотелий вплоть до денудации сосудистой стенки [5].
Для прогнозирования атеросклеротического процесса общепринятым стало вычисление таких характеристик потока, как полные и осевые сдвиговые напряжения, осредненные по времени сдвиговые напряжения и индекс колебаний сдвиговых напряжений.
Построение 3D-модели сосуда по данным компьютерной томографии
В качестве источника геометрической модели использованы данные обследования больного с извитостью общей сонной артерии, которому выпол-
нена компьютерная томография брахиоцефальных артерий с применением контраста по стандартной методике. Результаты исследования представляют собой файл, содержащий набор из 655 срезов (слай-сов) с шагом 0,75 мм. Размер каждого слайса — 512 х 512 пикселей. Обработка слайсов и ЭБ-модель выполнены с использованием методов обратного инжиниринга, трехмерная модель анатомической структуры создана с помощью оцифровки физического прототипа и с использованием данных обследования. Для выполнения этих задач задействованы приложения Slicer3D, GIMIAS, VolView, ParaView — мультиплатформенные программные продукты с открытым исходным кодом.
Чтение данных из файла формата DICOM выполнено программой 3DSlicer 4.0.1. На первом этапе применяется метод сегментации — разбиение изображения на непохожие по интенсивности серого цвета области. Предполагается, что области соответствуют реальным объектам или их частям, а границы областей — границам объектов. Методы сегментации можно разделить на два класса: автоматические, не требующие взаимодействия с пользователем, и интерактивные, использующие пользовательский ввод непосредственно в процессе работы. Затем данные экспортируются в приложение GIMIAS, позволяющее выполнять построение и редактирование сетки, содержащей до нескольких сот тысяч элементов. На заключительном этапе данные геометрии с помощью программы ParaView сохраняются в формате STL и могут быть использованы для задания граничных условий (рис. I, см. обложку).
Построение математической модели кровотока
Для описания течения в общей сонной артерии с S-образной извитостью использована модель вязкой несжимаемой жидкости на основе уравнений Навье—Стокса
V • V = о,
+ Р(т^ ^ V = ^ + + (^?)*),
где V = i + j + k — оператор набла; V — ох ду дz
вектор скорости жидкости; p — давление жидкости; ( )* — (сопряженный тензор). Для крови выбрана ньютоновская модель жидкости, с плотностью р = = 1000 кг/м3 и коэффициентом динамической вязкости т = 0,004 Па • с. Течение ламинарное, пульсирующее. Для стенки артерии выбрана жесткая модель. Вопрос влияния упругости стенки предполагается рассмотреть в дальнейших исследованиях. На рис. II (см. обложку) изображена вычислительная сетка, которая построена в программе-генераторе сеток 1СЕМ CFD у.12 и содержит порядка 300 000 элементов.
0,025
0,02
0,015
ч ©
0,01
0,005
t / T
Рис. 2
Зависимость расхода жидкости от времени на входе в общую сонную артерию [10]
Кривая расхода на входе в левую общую сонную артерию (рис. 2), измеренная методом фазово-контраст-ной магнитно-резонансной ангиографии, взята из работы [10]. Длительность цикла T составляет 0,7 с.
В работе рассмотрены два случая задания входного профиля скорости — плоский профиль и параболический профиль. Цель — показать чувствительность течения к входному граничному условию.
Расчеты выполнены в программном комплексе ANSYS CFX v.12, который является мощным современным средством для решения задач гидрогазодинамики и теплообмена по методу контрольных объемов. Конвективные слагаемые аппроксимированы по высокоточному методу High Resolution. Метод решения уравнений — в конечно-объемной формулировке, в которой связь скорости и давления осуществляется по схеме Rhie-Chow. Шаг по времени равен 0,01 с, количество итераций на одном шаге по времени — 30. Полученные решения являются численно независимыми от сетки, от шага по времени, от числа итераций на каждом шаге по времени.
Анализ влияния входного профиля скорости на структуру течения
Влияние входного граничного условия (плоский или параболический профили скорости) проанализируем на течении во втором изгибе общей сонной артерии (рис. 3). В выбранном для анализа сечении в фазу уменьшения расхода (Ь / T = 0,43) формируется самый значительный по длине артерии отрыв потока с низкими значениями сдвиговых напряжений на стенке.
Отличия между значениями скорости и сдвиговых напряжений на стенке при разных входных профилях скорости в фазе замедления потока (рис. 4, 5) более значительны по сравнению с фазой ускорения потока. Особенно это касается поля
Рис. 3
Отрыв потока в фазу замедления течения
(Ь / Т = 0,43) на втором изгибе общей сонной артерии
Рис. 4 Профили скорости в поперечном сечении второго изгиба общей сонной артерии в момент времени t/T = 0,43 в расчетах с разными профилями скорости на входе: а — плоский; б — параболический
скорости, когда в фазе замедления потока профили скорости имеют качественные различия. Плоский профиль скорости создает большие значения максимальной скорости и максимальных сдвиговых напряжений на стенке, чем параболический. Сред-
И
40 80 120 160 200 240 280 320 360
j, °
Рис. 5
Величина сдвиговых напряжений на стенке в поперечном сечении второго изгиба общей сонной артерии в момент времени Ь / Т = 0,43:
1 — плоский вид входного профиля скорости;
2 — параболический вид входного профиля скорости
0
4
3
0
ние по контуру сдвиговые напряжения в расчетах с плоским входным профилем скорости на всем протяжении цикла превышают на 5—8 % напряжения, полученные в расчетах с параболическим профилем скорости. Таким образом, расчеты с плоским профилем скорости дают оценку снизу по величине зон минимальных сдвиговых напряжений и тем самым указывают наиболее вероятные зоны атеросклеро-тических осложнений. В дальнейшем для анализа используются расчеты с плоским профилем скорости на входе.
Гемодинамический анализ течения
Гемодинамический анализ течения в общей сонной артерии призван показать зоны, опасные с точки зрения развития атеросклероза. Для этого рассчитываются локальные и осредненные распределения сдвиговых напряжений на стенке сонной артерии, индекс колебаний сдвиговых напряжений и показываются области с малыми значениями сдвиговых напряжений на стенке и большими значениями индекса колебаний сдвиговых напряжений, которые, как известно, являются факторами развития атеросклероза.
Осредненные сдвиговые напряжения AWSS (time-averaged wall shear stress) определяются по формуле
T
AWSS = 1Г It \dt,
Т J ' w
T J w 0
где |tw| — величина мгновенного вектора сдвиговых напряжений на стенке tw ; T — длительность сердечного цикла.
Индекс колебаний сдвиговых напряжений OSI (oscillatory shear index) определяется по следующей формуле
f
OSI =
J twdt
1 -
Л
0
Мгновенные (рис. 6) и осредненные по времени (рис. III, a, см. обложку) распределения сдвиговых напряжений меньше по величине в местах отрыва потока, который возникает в поворотах сосуда, составляют порядка 0,5 Па. Пониженные, опасные с точки зрения развития атеросклероза сдвиговые напряжения на стенке наблюдаются в фазе замедления потока и винтовым образом ориентируются вдоль общей сонной артерии. Осредненное поле сдвиговых напряжений на стенке (рис. III, a, см. обложку) показывает, что наиболее опасные места для развития атеросклероза — зоны перед и за первым поворотом, а также за вторым поворотом общей сонной артерии. В них величина сдвиговых напряжений на стенке не достигает физиологического
3,5
3,0
2,5
2,0
CS И 1,5
а р 1,0
0,5
0
0,5
1,0
/
£ / Vv * \
V—т^ГС 1J Ii
f>2
\/1 /
1 » \ / / 1
/
40 80 120 160 200 240 280 320 360
j, °
Рис. 6 Полные (1) и осевые(2) сдвиговые напряжения на стенке в поперечном сечении второго изгиба общей сонной артерии в момент времени í / Т = 0,43:
1 — полные WSS; 2 — осевые
значения 1 Па. Все эти области характеризуются тем, что в их окрестности происходит отрыв потока вследствие действия центробежной силы, которая появляется тогда, когда частицы жидкости движутся по искривленным траекториям. Центробежная сила отклоняет поток вдоль радиуса кривизны от внутренней стенки к внешней. Таким образом, на внутренней стенке поворота скорость становится меньше, чем на внешней. При небольшом радиусе кривизны поворота и большой скорости потока за поворотом формируется зона возвратного течения.
Сдвиговые напряжения на стенке растут с увеличением расхода (рис. 7). Величина сдвиговых напряжений совпадает с их осевой составляющей в фазу ускорения потока и сильнее всего отличается в фазу замедления потока, когда возникает локальное вих-реобразование, приводящее к появлению поперечных сдвиговых напряжений. В конце пульсации
6
CS
И
0,6
t / Т
Рис. 7
Характерное изменение величины сдвиговых напряжений по времени цикла в зонах низких напряжений
0
5
4
3
2
1
0
0,5
0,4
0,3
м о
0,2
0,1
• •!
•7
.Crt i * SV * » • • • • •
— tO INsS • * * •• •
2
AWSS, Па
Рис. 8
Соотношение индекса колебаний сдвиговых напряжений и осредненных по времени сдвиговых напряжений на стенке общей сонной артерии
появляется возвратное течение, что вызывает сдвиговые напряжения, направленные против основного потока (отрицательные осевые сдвиговые напряжения на стенке на рис. 6). В зонах возвратного течения (рис. III, а, см. обложку) большую часть цикла (примерно 80 %) сохраняются малые значения сдвиговых напряжений на стенке (рис. 7).
Индекс колебаний сдвиговых напряжений также используется как показатель для оценки риска развития атеросклероза. Наиболее опасны его большие значения. Индекс колебаний сдвиговых напряжений максимален (порядка 0,3) в местах отрыва потока (рис. III, б, см. обложку), что характеризует большую частоту смены направления вектора сдвиговых напряжений в течение цикла.
Соотношение малых сдвиговых напряжений и больших значений индекса колебаний сдвиговых напряжений, встречающееся в литературе [11], было также построено для рассматриваемого сосуда (рис. 8). Видно, что только малые значения сдвиговых напряжений характеризуются частой сменой направления, а большие значения сдвиговых напряжений чаще всего постоянны по направлению.
Заключение
Разработан алгоритм построения геометрии кровеносных сосудов на основе метода сегментации по данным компьютерной томографии, с помощью которого построена 3Б-модель общей сонной артерии с Б-образной извитостью. Проведен расчет течения в общей сонной артерии с использованием модели вязкой несжимаемой жидкости на основе трехмерных нестационарных уравнений Навье—Стокса в программном комплексе АЫБУБ СЕХ.
Рассмотрено влияние плоского и параболического профилей скорости на входе в артерию на скорость течения и сдвиговые напряжения на стенке. Влияние более значительно в фазу замедления потока, нежели в фазе ускорения. Особенно это касается поля скорости, когда в фазе замедления потока профили скорости даже качественно не совпадают. Плоский профиль скорости создает большие значения максимальных скорости и сдвиговых напряжений на стенке сосуда, чем параболический. Средние по контуру сдвиговые напряжения в расчетах с разными входными профилями скорости в фазе увеличения расхода отличаются не более чем на 10 %. Поскольку расчеты с плоским профилем скорости на входе дают более высокие значения сдвиговых напряжений, то благодаря им получена оценка снизу по величине зон минимальных сдвиговых напряжений. Таким образом, указаны наиболее вероятные зоны атеросклеротических осложнений.
Осредненное поле сдвиговых напряжений на стенке показывает, что наиболее опасные места для развития атеросклероза — зоны с отрывом потока, которые наблюдаются с внутренней стороны изгибов сосудов. В них величина сдвиговых напряжений на стенке мала и не достигает нижней границы физиологического значения — 1 Па. Отрыв приводит к тому, что в фазу замедления потока формируется небольшое возвратное течение. Индекс колебаний сдвиговых напряжений достигает максимального значения 0,3 в местах отрыва потока.
Литература
1. Петров И. Б. Математическое моделирование в медицине и биологии на основе моделей механики сплошных сред // Труды МФТИ. 2009. Т. 1, № 1. С. 5-16.
2. Steinman D. A., Thomas J. B., Ladak H. M. et al. Reconstruction of carotid bifurcation hemodynamics and wall thickness using computational fluid dynamics and MRI // Magnetic Resonance in Medicine. 2002. Vol. 47. P. 149-159.
3. Павлова О. Е., Иванов Д. В., Грамакова А. А. и др. Гемодинамика и механическое поведение бифуркации сонной артерии с патологической извитостью // Изв. Саратовск. ун-та. Нов. сер. 2010. Т. 10, вып. 2. С. 66-73.
4. Manbachi A., Hoi Y., Wasserman B. A. et al. On the shape of the common carotid artery with implications for blood velocity profiles // Physiol. Meas. 2011. Vol. 32. P. 1885-1897.
5. Fry D. Acute vascular endothelial changes associated with increased blood velocity gradients // Circulat. Res. 1968. Vol. 22. P. 165-197.
6. Caro C., Fitz-Gerald J., Schroter R. Atheroma and arterial wall shear: observation, correlation, and proposal of a shear dependent mass transfer mechanism for atherogenesis // Proc. R. Soc. Lond. (Biol). 1971. Vol. 177. P. 109-159.
7. Asakura T., Karino T. Flow patterns and spatial distribution of atherosclerotic lesions in human coronary arteries // Circulat. Res. 1990. Vol. 66. P. 1045-1066.
8. Leitinger N., Oguogho A., Rodrigues M. et al. The effect of NO/ERF and monocytes/macrophages on LDL-oxidation // J. Physiol. Pharmacol. 1995. Vol. 46, N 4. P. 385-408.
9. Hutcheson I., Griffith T. Heterogeneous populations of K+ channels mediate ERF release to flow but not agonists in rabbit aorta // Am. J. Physiol. 1994. Vol. 266. P. 590-596.
0
1
3
4
Иллюстрации к статье Я. А. Гатаулина, А. Д. Юхнева, М. А. Попова, Д. И. Курапеева
«ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КРОВОТОКА В ОБЩЕЙ СОННОЙ АРТЕРИИ С Б-ОБРАЗНОЙ ИЗВИТОСТЬЮ»
С. 27-33.
Рис. I. КТ-изображение дуги аорты и брахиоцефальных сосудов
Рис. II. Вычислительная сетка общей сонной артерии с Б-образной извитостью
Рис. III. Осредненные по времени сдвиговые напряжения (а) и индекс колебаний сдвиговых напряжений (б) на стенке общей сонной артерии с S-образной извитостью (стрелками показаны зоны, предрасположенные к развитию атеросклероза)
