CC BY
70
УДК 519.63
МОДЕЛИРОВАНИЕ ДВИЖЕНИЯ ВЯЗКОЙ НЕОДНОРОДНОЙ ЖИДКОСТИ В КРУПНЫХ
КРОВЕНОСНЫХ СОСУДАХ
Д. А. Долгов, Ю. Н. Захаров
MODELLING OF VISCOUS INHOMOGENEOUS FLUID FLOW IN LARGE BLOOD VESSELS
D. A. Dolgov, Yu. N. Zakharov
Работа выполнена при поддержке проектной части Государственного задания № 1.630.2014/К.
В работе рассматривается математическая модель, описывающая течение неоднородной несжимаемой жидкости с переменной вязкостью в канале с гибкими стенками, а также метод ее численного решения. Приведены результаты моделирования аневризмы стенки кровеносного сосуда, а также движения примеси в нем.
In this paper we propose a mathematical model for describing the viscous inhomogeneous fluid flow in a canal with flexible walls. We present the results of modeling of a blood vessel aneurysm, and the flow of admixture inside the vessel.
Ключевые слова: вязкая неоднородная жидкость, аневризма сосуда, метод погруженной границы.
Keywords: viscous inhomogeneous fluid, aneurysm of blood vessel, immersed boundary method.
Введение
Трудно переоценить важность медицинских исследований кровеносной системы человека, т. к. знания в этой области чрезвычайно практичны и значимы. В последнее время интерес к математическому моделированию кровеносной системы существенно возрос (см. [16] и цитируемую здесь литературу). Существует множество исследований по данной тематике, в результате которых выделилось два основных подхода к решению подобных задач.
Первый подход связан с использованием конечно-элементных методов для моделирования течения крови в сосудах [3; 15; 17]. Он позволяет хорошо учитывать сложную геометрию задачи, однако необходимость учитывать взаимодействие жидкости и гибких стенок требует постоянного перестраивания расчетной сетки, чтобы удовлетворять меняющейся геометрии исследуемого объекта. Это приводит к существенным затратам работы компьютерных программ расчётов.
Второй подход, который и рассматривается в данной работе, связан с методом погруженной границы [4; 10; 12; 14]. Метод также может применяться в задачах со сложной геометрией, но при этом не требует модификации сетки.
В связи с необходимостью моделирования более сложных задач, появляются развития метода погруженной границы. В [5] предложена формулировка этого метода для случая трехмерного течения двух не перемешивающихся (разделенных гибким препятствием) компонент жидкости разной вязкости и плотности. В [8; 9] показано применение метода для решения двумерных задач о течении двухкомпонентной жидкости.
В данной работе мы предлагаем описывать движение крови в упругих крупных кровеносных сосудах как трехмерное нестационарное течение вязкой несжимаемой жидкости с переменной плотностью и вязкостью [6; 7; 11]. Таким образом, целью работы является построение математической модели и метода решения задачи о возникновении и развитии аневризмы стенок кровеносных сосудов с учетом неоднород-
ной структуры крови, а также о движении примеси (форменных элементов) внутри сосуда.
Постановка задачи
Рассмотрим нестационарную задачу о течении крови внутри сосуда. Кровь состоит из плазмы и взвешенных в ней форменных элементов. Стенки сосуда являются гибкими и изменяют свою форму в зависимости от течения крови. Сосуд может быть подвержен аневризме, т. е. "выпячиванию" стенок вследствие их растяжения.
Исходя из этого, будем моделировать кровь как вязкую, несжимаемую неоднородную двухкомпо-нентную жидкость с переменной вязкостью, стенки сосуда - как непроницаемую для жидкости поверхность цилиндрической формы, обладающую определенной жесткостью. Под воздействием давления жидкости сосуд деформируется. Аневризму, возникающую на стенках, будем моделировать с помощью локального уменьшения ее жесткости.
Так как источником движения крови в сосудах является давление, создаваемое сокращением сердца, то задачу о ее движении опишем следующей нестационарной системой дифференциальных уравнений Навье-Стокса [7]:
nu 1 ^
— + (и -V)u =--Vp + Va + f
dt p
V-и = 0
с начальными и краевыми условиями:
(1) (2)
= 0
с начальными условиями:
c(x,0) = c0(x),x eQ
и краевыми условиями на границе втекания: c( Х, t ) 1г2 = cs ( Х, t )-
суд в равновесное состояние. В силу того, что в данной работе рассматриваются небольшие деформации стенок сосуда, мы используем простую формулу для нахождения силы деформации в зависимости от смещения относительно исходного положения:
и (x ,0) = uo, и г = ив, и 1г2, г3 = о (3)
P iГ2 = Pin , P iГ2 = Pont (4)
где x = ( x, y, z ) eQ , и = (и, V, w) - вектор
скорости, U в - скорость, с которой двигаются стенки сосуда при деформации, p = p(x, t) - плотность, p = p(x, t) - давление,
a = ju(Vu +(Vu) ) - вязкий тензор напряжений,
j = j(x, t) - вязкость жидкости, f = f (x, t) -вектор массовых сил, который в дальнейшем используется для определения формы сосуда.
Область Q - сосуд с границей
Г =Г1 иГ2 иГ3,
где Г1 - стенка кровеносного сосуда, Г2 и Г3 - области втекания и вытекания (рис. 1).
Плотность p и вязкость jl определяются следующими соотношениями [7]:
j = c(j -j) + jj, (5)
P = c(P2 -Pl) + (6)
где Pi, j - плотность и вязкость жидкости (плазмы), а P2, j - плотность и вязкость примеси (форменных элементов), С - концентрация примеси. Концентрация С = c(x, t) , С e [0,1] примеси определяется как решение уравнения:
дс _
_ + и-Vc = 0 dt
F = k||X - X„||,
(10)
(7)
(8) (9)
где X, Х0 - функции, которые описывают поверхность сосуда в момент времени ? и в начальный момент времени, к - коэффициент жесткости.
Для более общего случая могут быть использованы формулы, приведенные в работах [4; 13].
Метод решения
В соответствии с методом погруженной границы [12], будем рассчитывать течение жидкости в параллелепипеде О , включающем в себя О , на границах которого задано условие прилипания. Для расчета течения жидкости будем использовать прямоугольную равномерную разнесенную сетку О И с шагами
И^, И^, Игк , а для определения деформации поверхности сосуда Г1 введем дополнительную область Г с лагранжевой системой координат, соотнесенной со стенками сосуда. На границе Г построим сетку Г и , узлы которой соответствуют точкам на Г . Алгоритм решения состоит из нескольких шагов: на сетке О И решаем задачу (1) - (4); затем решаем уравнение конвекции (7), т. е. определяем концентрацию примеси в области решения и пересчитываем значение плотности и вязкости. Далее определяем деформацию сосуда под воздействием жидкости, т. е. находим смещение узлов сетки ГИ относительно предыдущего момента времени, и с помощью уравнения (10) вычисляем значение сил, противодействующих деформации. После этого находим новое распределение массовых сил в уравнении движения жидкости.
Поставленная дифференциальная задача (1) - (9) решается методом конечных разностей. Для решения (1) - (4) будем использовать схемы расщепления по физическим факторам [1]:
= -(ип - V)u* - - Va + fn At p
Здесь С 0, - заданные функции.
Отсутствие задания одной компоненты вектора скорости на участках втекания-вытекания является одной из проблем при численном решении задач подобного типа. Она решается с помощью использования исходных уравнений (1) - (4) на границах Г2 и
Г3 для вычисления недостающих компонент вектора
скорости (подробнее см. [7]).
Для моделирования движения стенок сосуда необходимо определять силы, которые возвращают со-
pApn+i - (Vp- pn+i) =
п+1
P2Vu At
и
и
1
At
=--Apn
P
(11)
(12)
(13)
Численная реализация схемы состоит из 3-х этапов. Сначала по известным значениям скорости с
предыдущего временного слоя находится промежу-
*
точное поле и . Для этого уравнение (11) решается методом стабилизирующей поправки [2]. Затем, путем численного решения (12) с использованием мето-
да бисопряженных градиентов, определяется новое поле давления. И на последнем этапе восстанавливается окончательное поле вектора скорости по явным формулам (13).
После нахождения параметров течения жидкости необходимо вычислить новые значения плотности и вязкости. Для этого, используя полученные значения компонент скорости, делается шаг по времени для уравнения конвекции (7), и производится пересчет значений плотности и вязкости по формулам (5), (6).
Далее нам необходимо определять деформацию стенок сосуда под воздействием жидкости, а также распределение массовых сил в уравнении движения жидкости исходя из возникшей деформации. В работах [4; 12] показано, что для этого мы можем использовать следующие уравнения:
ах
ÔX
-(q, t ) = F u ( x, t) -S( x - X (q, t ))dxdydx (14)
мени. Точками обозначена форма сосуда в данный момент времени.
На рис. 2 показано изменение формы стенок сосуда под воздействием жидкости за счет локального снижения жесткости стенки. В верхней половине центральной части сосуда коэффициент жесткости имеет
значение k = 0.4-103. Плотность и вязкость были постоянны.
îïîîssîîïï
f (x, t) = J F (q, t) 8(x - X (q, t))dqdrds (15)
г
где q = (q, r, s) £ Г - точка поверхности сосуда,
X = X (q, t) - функция, описывающая поверхность сосуда в момент времени t, F = F (q, t) - сила сопротивления деформации в данной точке сосуда, u (x, t) - вектор скорости течения,
f (x, t) - вектор массовых сил,
8 - дельта-функция Дирака.
Используя уравнения (14), (15), которые численно интегрируются с помощью какой-либо квадратурной формулы, и уравнение (10), мы можем рассчитать деформацию, которой подвергаются стенки сосуда при данном давлении жидкости, и то, как повлияют возникшие силы сопротивления на течение.
Результаты
Приведем результаты некоторых методических расчетов с различными коэффициентами жесткости для случаев постоянной и переменной вязкости и плотности, цель которых заключается в демонстрации работоспособности описанного метода и возможности получать с его помощью картины деформации стенок сосуда и распространения примеси. Во всех расчетах был использован сосуд, начальная форма которого
являлась круговым цилиндром с длинной l = 1, радиусом r = 0.11 и жесткостью стенок k = 4 • 103 . Давление на входе и на выходе pin = 2 , pout = 1.
На рис. 2 - 4 приведены расчёты движения жидкости и стенки сосуда, где область Q имела параметры 1.0 X 0.5 X 0.5, шаги по пространственной
сетке ^h hx = hy = hz = 0.01, шаг по времени At = 0.01, n указывает количество шагов по вре-
а)
ï S 8 ï î Î %
б)
Рис. 2. Изменение формы сосуда под действием течения.
Внутри сосуда отображены треки некоторых
частиц рх = рг = 1, ¡Л1 = jU2 = 1 • 10 2 а) n = 100, б) n = 400, в) n = 1200
Как видно из рис. 2, в области, где стенка имеет меньшую жесткость, сосуд деформируется под воздействием давления жидкости.
Второй пример (рис. 3) отображает распространение примеси в жидкости при образовании «аневризмы». На Г2 задается постоянный приток примеси с концентрацией cs = 0.5.
Рис. 3. Изменение формы сосуда и распространение в нем примеси. На входе задается постоянный
приток примеси рх = 1, р2 = 2 , ¡Лх = 1 -10 2,
р2 = 2 -102, с, |Г2 = 0.5 а) п = 0, б) п = 400, в) п = 1200
Как следует из рис. 3, «аневризма» небольшого размера мало влияет на распространение примеси.
Третий пример показывает распространение уже имеющейся внутри сосуда примеси. В центре сосуда
имеется сужение до радиуса г2 = 0.09 . Центральная
часть заполнена вязкой примесью с концентрацией
с0(X, у, £) = 1 в области (X, у, £) еОс и
с0 (X, y, z) — 0 в остальной части
Рис. 4. Размытие примеси в деформированном сосуде. Область Q c в центре сосуда имеет повышенную плотность и вязкость рх — 1,
р2 — 2 , р1 — 1 • 10, р2 — 2 • 10~2 а) n = 0, б) n = 50, в) n = 150
Этот пример показывает, область в центре сосуда с повышенной вязкостью и плотностью размывается потоком жидкости меньшей вязкости и плотности.
Заключение
Построенная модель течения крови с переменной плотностью и вязкостью позволяет получать картины изменения стенок сосуда и распространения примеси под воздействием течения неоднородной жидкости.
сосуда.
Литература
1. Белоцерковский О. М. Численное моделирование в механике сплошных сред. М.: Наука, 1984. 520 с.
2. Яненко Н. Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск.: Наука, 1967. 197 с.
3. Black M. M., Howard I. C., Huang X., Patterson E. A. A three-dimensional analysis of a bioprosthetic heart valve // J. Biomech. 1991. 24(9). P. 793 - 801.
4. Boyce E. G. Immersed boundary model of aortic heart valve dynamics with physiological driving and loading conditions // International Journal for Numerical Methods in Biomedical Engineering. 2011. 1 - 29.
5. Fai T. G., Boyce E. G., Mori Y., Peskin C. S. Immersed boundary method for variable viscosity and variable density problems using fast constant-coefficient linear solvers I: Numerical method and results // SIAM Journal on Scientific Computing. 2013. 35(5). B1132 - B1161.
6. Geidarov N. A., Zakharov Y. N., Shokin Yi. I. Solution of the problem of viscous fluid flow with a given pressure differential // Russian Journal of Numerical Analysis and Mathematical Modeling. 2011. V. 26. № 1. P. 39 - 48.
7. Gummel E. E., Milosevic H., Ragulin V. V., Zakharov Y. N., Zimin A. I. Motion of viscous inhomogeneous incompressible fluid of variable viscosity // Zbornik radova konferencije MIT 2013. Beograd. 2014. 760 p. (Proceedings of International Conference "Mathematical and Informational Technologies MIT-2013" Врнячка Баня, Сербия, Будва, Черногория, 5.09.2013 - 14.09.2013). Kosovska Mitrovica. 2014. P. 267 - 274.
8. Jian D., Robert D. G., Aaron L. F. An immersed boundary method for two-fluid mixtures // Journal of Computational Physics. 2014. Volume 262. P. 231 - 243.
9. Lee P., Boyce E. G., Peskin C. S. The immersed boundary method for advection-electrodiffusion with implicit timestepping and local mesh refinement // Journal of computational physics. 2010. Volume 229. P. 5208 - 5227.
10. Ma X., Gao H., Boyce E. G., Berry C., Luo X. Image-based fluid-structure interaction model of the human mitral valve // Computers & Fluids. 2013. 71. P. 417 - 425.
11. Milosevic H., Gaydarov N. A., Zakharov Y. N. Model of incompressible viscous fluid flow driven by pressure difference in a given channel // International Journal of Heat and Mass Transfer. 2013. July. Vol. 62. P. 242 - 246.
12. Peskin C. S. Numerical Analysis of Blood Flow in the Heart // JC. 1977. 25. P. 220 - 252.
13. Peskin C. S. The immersed boundary method // Acta Numerica. 2002. 11. P. 479 - 517
14. Pilhwa L., Boyce E. G., Peskin C. S., The immersed boundary method for advection-electrodiffusion with implicit timestepping and local mesh refinement // Comput Phys. 2010. July 1. 229(13).
15. Taylor C. A., Hughes T. J. R., Zarins C. K. Finite Element Modeling of Blood Flow in Arteries // Computer Methods in Applied Mechanics and Engineering. 1998. Vol. 158. P. 155 - 196.
16. Yoganathan A. P., He Z. M., Jones S. C. Fluid mechanics of heart valves // Annu. Rev. Biomed Eng. 2004. Vol. 6. P. 331 - 362.
17. Zhang Y, Bajaj C. Finite element meshing for cardiac analysis // ICES Technical Report. 2004. Р. 4 - 26.
Информация об авторах:
Долгов Дмитрий Андреевич - аспирант кафедры вычислительной математики КемГУ.
Dmitry A. Dolgov - post-graduate student at the Department of Computational Mathematics, Kemerovo State University.
(Научный руководитель - Ю. Н. Захаров).
Захаров Юрий Николаевич - доктор физико-математических наук, профессор, заведующий кафедрой вычислительной математики КемГУ, [email protected].
Yury N. Zakharov - Doctor of Physics and Mathematics, Professor at the Department of Computational Mathematics, Kemerovo State University.
Статья поступила в редколлегию 12.03.2015 г.
