Неявные разностные схемы для одномерных уравнений гемодинамики Текст научной статьи по специальности «Математика»

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ И

ПРОЦЕССЫ УПРАВЛЕНИЯ N. 3, 2023 Электронный журнал, рег. Эл. N ФС77-39410 от 15.04.2010 ISSN 1817-2172

http://diffjournal, spbu. ruf e-mail: jodiff@mail.ru

Прикладные задачи

Неявные разностные схемы для одномерных уравнений

гемодинамики

Г. В. Кривовичев, Н. В. Егоров

Факультет прикладной математики - процессов управления Санкт-Петербургского государственного университета

E-mail: g.krivovichev@spbu.ru, n.v.egorov@spbu.ru

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

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

по пространству. При проведении вычислительных экспериментов по моделированию кровотока в модельных сосудистых системах показано, что разработанные схемы позволяют производить расчеты за меньшее время, чем известные из литературы явные конечно-разностные и конечно-объемные схемы.

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

1 Введение

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

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

Конечно-объемные схемы особенно часто используются в задачах одномерной гемодинамики. При этом применяются различные подходы к реконструкции потока на границах конечных объемов и для построения монотонных схем. В [7, 8, 9] предложены консервативные схемы, точно воспроизводящие стационарные решения в случае нулевой скорости течения (т.н. dead-man equilibria). В [10] для моделирования течений в различных отделах аорты авторы предложили использовать конечную-объемную WEN О (weighted essentially non-oscillatory) схему. В [11] для проведения расчетов используются TVD (total variation diminishing) схемы. В [7, 8, 12] в рамках метода

конечных объемов предложено использовать схемы типа MUSCL (monotonie upwind scheme for conservation laws) и ENO (essentially non-oscillatory). В [12] эти схемы сравниваются друг с другом и со схемами, предложенными другими авторами.

Конечно-элементные схемы для уравнений одномерной гемодинамики, как правило, строятся на основе метода Галеркина [13] и Тэйлора — Галерки-на [14, 15, 16]. Схемы на основе последнего можно рассматривать как обобщения схемы Лакса — Вендроффа на случай конечно-элементных аппроксимаций по пространству.

Как отмечают авторы работы [12], конечно-разностные схемы для задач одномерной гемодинамики являются наиболее удобными для практической реализации. В задачах моделирования течения крови наиболее популярны схемы второго порядка аппроксимации, особенно схемы Лакса — Вендроффа и МакКормака [17, 18, 19, 20, 21, 22]. В [17] при решении конкретных задач показано, что такие схемы позволяют получать близкие по точности результаты. Из-за простоты в реализации и удобства распараллеливания такие схемы используются в том числе и в пакетах прикладных программ [23].

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

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

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

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

2 Основные уравнения

С позиций механики кровь описывается как вязкая несжимаемая жидкость [27], динамика которой описывается системой уравнений, состоящей из условия несжимаемости и уравнения движения. В результате осреднения по поперечному сечению сосуда эта модель сводится к одномерной модели [28, 29]:

дА дQ л дQ д ( д2\ АдР £,Л ^

Ж + Ж = 0 -Ж + дДат) + Адр = М

где £ — время, г — пространственная координата, А(£, г) — площадь поперечного сечения, Q(t, х) — объемный расход, а > 1 — коэффициент Буссинеска, р _ плотность, Р(£, г) — давление, f — вязкий член, вычисляемый таким образом:

, . 2пЯ .

/ = -Ттг\т=Е,

р

где Я(£, г) — радиус сосуда, Тгх — компонента тензора касательных напряжений. В большинстве работ по одномерным моделям кровь рассматривается как ньютоновская жидкость. Для такой модели вязкий член представляется

как [30]:

! (А,Я) = —кА,

где

а 2пр

к =--,

а — 1 р

где р = оспвЬ — коэффициент динамической вязкости.

А

и средней скорости и:

дА д(Аи) ди д ( и2 \ 1 дР _ и

~дЬ + дг = , ~дЬ + + ~рдг = — А' ^

Систему (2) можно получить из (1) в случае значения а =1, когда оно рассматривается только в левой части уравнения для А в (1). В работах [30, 31] показано, что значение этого коэффициента, рассматриваемое только в левой части уравнения, существенно не влияет на результаты расчетов течения крови в модельных сосудистых системах. В связи с этим на прак-

а=1

уравнения движения.

В случае, когда вязкие свойства крови не учитываются, (2) сводится к системе вида [4, 32, 33, 34, 35]:

дА д(Аи) ди д и2 1 дР

ж + -V = 0, ж + + = 0' (3)

Системы (2) и (3) замыкают уравнением состояния Р = Р(А). Как и во многих других работах, будем использовать уравнение, применяемое для моделирования течений в артериях [4]:

в(VÄ -JA), (4)

Р — РвхЬ = Р3, + (

Аа

где Рехг — внешнее давление, Аа и Ра — диастолические площадь и давление, в = т^л/лЕЬ, Е — модуль Юнга, Н — толщина стенки сосуда. Будем полагать, что для конкретного сосуда в системе эти параметры являются постоянными.

После подстановки (4) в (2) и (3), уравнения движения примут вид:

ди д (и2\ .... и

dU д (U2\

ж + Ш[т + iiA) = 0<

(6)

где

^(A) = х

д VA

dz '

X =

в

pAd"

3 Разностные схемы

Рассмотрим сетку по времени, построенную с шагом Д£ и узлами £п, п = 1, М, и сетку по пространству с узлами г^, г = 0,Ы + 1, построенную с шагом Н. Аппроксимируем осредненное условие несжимаемости

дА д (Аи)

+

= 0

д£ дг

и

как известная функция и*:

А^1 - Аг

At

+ L^(An+1) = 0,

где Ап = (Ао, АА^)т, Ап « А(ЬП) г), г = 0, N + 1: тор СА задается как:

(7)

а сеточный опера-

Г An+1U| - An+1U0*

i = 1,

1 2h

An+aU»-2 - 8An+11U*_ 1 + 8An+11U;+1 - An+^U*

4+1 ^ i+1

4+2 ui+2

12h

л n+1 T т* _ лп+1 т т* AN+1UN+1 aN-1uN-1

2h ,

i = 2, N - 1,

i = N,

где узлы го и гN+1 являются граничными.

Аппроксимируем (5) разностной схемой вида:

Uп+1 _ Uп

+ Lu (U n+1)+ ^(A*) = -K

U п+1

Д£ ..............."А" (8)

где сеточная функция ип вводится то аналогии с Ап, а оператор Си и ф(А*) задаются как:

f U2n+1U2* - U0n+1 u0*

i = 1,

Lu (Uin+1) = I

4h

Ur_+21Ui*_2 - 8Uin+L1Ui*_ 1 + 8Uin+1Ui*+1 - U+^U*

i+1 i+1

i+2 i+2

24h

ТТП+1т т* _ TTn+1T T*

UN+1UN+1 UN-1UN-1

4h ,

i = 2, N - 1,

i = N,

УД - у/Л* . 1

-тгт^-, i = 1,

2h ' '

m¡) = <

^Дй - z^Ati + 8\/Л+1 - vA+2

i-2 "V ^i-1

_ ^ 12h

\/an+1 - \JaN-1

i = 2,N — 1,

i = N.

2Н

при этом будем предполагать, что функция А* является известной. Аналогичным образом аппроксимируем уравнение (6): иП+1 _ иг

At

+ Lu (Un+1)+ ф(Л*) = 0. (9)

Процесс вычислений на одном шаге по времени [tn,tn+i] организуем следующим образом: на первом этапе по известному U* из (7) найдем значения An+\ а на втором этапе по известным A* и U* го (9) определим Un+1. Будем полагать, что U* = Un, A* = An+1.

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

Заметим, что в связи с тем, что A* и U* являются известными, схемы (7)-(9) являются линейными и их применение сводится к решению систем линейных алгебраических уравнений с пятидиагональными матрицами, для решения которых можно эффективно применять метод пятиточечной прогонки, который имеет сложность Q(N) [37].

4 Анализ устойчивости

Как и в большинстве работ, проведем исследование устойчивости только по начальным условиям, т.е. будем рассматривать случай неограниченного промежутка по пространству. Для анализа устойчивости воспользуемся методом гармоник (фон Неймана) [38, 39] в сочетании с принципом замороженных коэффициентов [40].

В схеме (7) «заморозим» значения U*, используя вместо них величину U = max |U*|. При введенных предположениях схема (7) примет следующий

г

вид:

дп+1 _ ли и

к к + ^(А—2 - 8АП+; + 8АП+; - лп+.]) = 0. (10)

Представим решение (10) в виде [38]:

АП = Ап(р)ехр(3рк), (И)

где р € [0,2п), ^2 = — 1, А(р) — спектральная функция. Проверим для схемы (10) выполнение спектрального критерия устойчивости |А(р)| < 1, Ур € [0, 2п). Подставляя (11) в (10), получим выражение для А:

1

Л =

Y

1 + Y2(exp(-2j^) - 8exp(-j^) + 8exp(j^) - exp(2j^))

где 7 = иАЬ/Н — число Куранта.

Упрощая выражение в знаменателе, получим:

а = - 1

1 + ^(8sin(^) - sin(2^)) 6

Из последнего выражения очевидно, что спектральный критерий устойчивости выполняется. Таким образом, схема (10) является безусловно устойчивой.

В схеме (8) произведем «заморозку» U*/2 и A*, рассматривая значения U = max |U*/2| и в = K max |1/A*|:

г г

иг1 +12 (un—2 - 8un_+i + 8un+11 - ий) + в Atur1 = щ. (12)

Представляя решение (12) аналогично (11), получим следующее представление для спектральной функции:

л =_ 1

1 + в At + — (8sin(p) — sin(2p)) 6

откуда следует выполнение спектрального признака устойчивости.

5 Вычислительные эксперименты

В настоящем разделе предложенные схемы расщепления будут применены к решению тестовых задач. Две задачи имеют аналитическое решение и при их численном решении покажем, что предложенные схемы сходятся с четвертым Н

расщепления с известными из литературы явными схемами.

5.1 Явные схемы

При проведении численных расчетов будем проводить сравнение схем (7)-(9) с явными схемами, которые строятся для системы уравнений гемодинамики, выписанной в дивергентной форме:

f+™

где U = (A, U)T.

Рассмотрим следующие схемы:

1. Схема Лакса — Вепдроффа. Эта схема на одном шаге по времени реализуется в два этапа: на первом производится расчет в полуцелых узлах

сетки:

где

им = 1 (U+1 + U) - At (F(U'+l) - F(U?)) + TB(u;,+1

U+1 = 1 (U+i + u?).

г+2 ~ 2К г+1

На втором этапе производится расчет значений в узлах с целыми индексами:

иг1 = и - Д (г (и+*) - Е (и^)) + ДШ(И?).

Схема имеет второй порядок сходимости по обеим переменным и является устойчивой при выполнении условия [17, 41]:

h

At <

гдес = J = J „-AA

max(|Ui| + ci)

p dA y '

2. Схема MUSCL (Monotonie upwind scheme for conservative laws). Эта схема основана на методе конечных объемов. При ее применении область разбивается на ячейки i, zi+1] с центрами в узлах zi. Средние значения решения по каждой ячейке вычисляются следующим образом [12]:

M

up = h У U(tn, z)dz.

Схема задается выражением:

ип+ = ип — А (гп+ 2 — ГП— 2) + Аьв(ип),

где Еп+ 1 и Еп_ 1 есть значения численного потока в граничных точках ячейки.

г+ 2 г 2

В рамках такой конечно-объемной схемы используется поток Русанова:

Fn+1 = - (Wun+ О + F fun, iЛ) — a(Un+1 + — Un+1 ),

2 2\ У г+ 2 -J V i+ 2 +JJ 2У г+ 2 + i+ 2

где a = max ^Лт 1 , Лт 1 , где Лт есть наибольшее по модулю собственное значение матрицы A систем (2)-(3), записанных в квазилинейной форме:

f + A(U, f = f

а Un+1 + и Un+ 1 есть предел справа Цм и предел слева U"- в узле zi +1. В

i+ 2 + i+ 2 2

2 2 a

посредством h/At, a U -5 U+ берутся как средние по соответствующим ячейкам, получим поток Лакса — Фридрихса и первой порядок аппроксимации по пространству. Но если использовать линейную реконструкцию потока, то будем иметь второй порядок [12].

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

hh

si - 2 + = si- 2Ds^ si+Ь = si + 2 г

где Dsi расчитывается таким образом:

si si -1 si+1 si h , h

Dsi = minmod

где

!min(x,y), x,y > 0, max(x,y), x,y < 0, 0, .

Условие устойчивости имеет следующий вид [7, 12]:

h

At

2 max(|Ui| + о)'

При проведении расчетов по явным схемам шаг по времени выбирался из условия:

At < ncFL min ( max (-, * . I I , (13)

s=1,S \i=1,N \ Cs,i + IUs,tU J где S есть число сосудов в системе.

5.2 Задачи на бесконечном промежутке

Рассмотрим модельные задачи, поставленные на всей вещественной прямой, для которых можно получить аналитические решения. Покажем, что численные решения, получаемые по предложенным схемам расщепления, сходятся с четвертым порядком по h.

Введем безразмерные переменные:

z ~ t Z A ~ U Z = ~T, t> = —, A = —, U = —, (14)

Lc Tc Ad Uc

где Lc и Tc есть характерные длина и время, Uc = Lc/Tc. В дальнейшем в этом разделе знак тильды у безразмерных переменных будем опускать.

5.2.1 Задача для невязкой модели

После подстановки (14) в (3), замкнутую посредством (4), получим систему вида:

dA д(AU) dU dU к dA

~dt + V = °, Ж + UdH + TAd-z =0, (15)

в

T~iTT/2i Lr* - -L_

ГДе K = 2 TAdpU2 ■

Предположим, что начальные условия представимы в виде:

A(0, z) = A0 + гф1(х) + £2ф2^) + ..., U(0, z) = eZi(z) + e2Qi{z) + ..., (16)

где £ есть малый безразмерный параметр, A0 > 0 есть постоянная величина, функции ф.1 (z), Zi(z) считаем ограниченными на всей прямой. Такое представление начальных условий отвечает идее возмущения стационарного решения

(A0, 0)

Будем искать решение задачи (15)—(16) в виде:

A(t, z) = Ao + £Ai + £2A2 + . . . , U(t, z) = £Ui + £2U2 + ....

Легко показать, что А0 = А0. В дальнейшем ограничимся только первым приближением по е. Для А1 и и1 получим систему:

дА1 дЩ ди к дА1

+ = 0, ^Г + ТА0~дх =0,

и1

д 2и1 2 д 2и1

= , (17)

dt2 dz2 '

1

где c = л/KAq .

0 А1

£

А1(Ь, г) = — Ао J дц-(т, г)((т + и)(г), 0

где /ш(г) определяется из начального условия. Для (17) получаются начальные условия следующего вида:

и1(0, г) = <,(,), *£(„,,) = — ^и). (18)

В случае С1(г) = С (г) ф1(г) = 0 получим решение (17)-(18) в виде суперпозиции бегущих волн:

и1(Ь,г ) = 2(С (г —сЬ) + ((г + сЬ)).

А1

А

А1(Ь,г) = АО (С (г — сЬ) — С (г + сЬ)).

При численном решении этой задачи вместо бесконечного промежутка будем рассматривать отрезок [—Ь, Ь] и временной промежуток [0,Т]. В граничных точках поставим условия совместности [42]:

(дИ . .тт.дИ'

ди + А±(И) ди

WU)^—+ A±(U)—]= 0, (19)

где 1± и А± есть левые собственные векторы и отвечающие им собственные значения матрицы А. В случае г = — Ь в (19) фигурируют А — < 0 и 1 —, в случае г = Ь используются А+ > 0 и 1+. Помимо (19) в граничных точках

зададим условия на значения и, вычисляемые посредством аналитического решения.

Производную по £, входящую в (19), аппроксимируем при £ = £п+1 с помощью левой разностной производной:

дИ ип+1 - Ип

диИ(£„+Ь*г) — ^Д-^, (20)

где г = 0 и г = N + 1. Производим е по г аппроксимируем конечными разностями второго порядка:

дИ \ -3Ип+1 + 4Ип+1 - Ип+1

дг ^ 2Н (о\)

дИ ч Ип+1 - 4Ип+1 + 3Ип+1 1 ;

\tn+i +l)

N-1 ^ N

дг '+1' " +1у 2Н

Нелинейные функции 1± (И) и Л±(И) аппроксимируем при £ = £п+1 следующим образом: 1±(Ип+1) — 1±(Ип), Л±(Ип+1) — Л±(Ип). Таким образом, в результате дискретизации условий совместности и с использованием условий для и получим выражения для значений А в граничных узлах в момент времени £ = £п+1.

При реализации неявных схем для расчета значений в приграничных узлах (% = 1 и г = требуется использовать значения в граничных узлах. Это приводит к необходимости решения глобальной системы уравнений относительно значений искомых функций. Для локальной реализации расчетов течений в конкретном сосуде (т.е. независимо от других сосудов), значения в граничных узлах должны быть известны. Локальная реализация важна ввиду возможности распараллеливания вычислений и характерна для явных схем (поскольку значения в граничных узлах для таких схем берутся с предыдущих слоев по времени).

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

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

тов, адекватные и превосходящие по производительности явные схемы результаты удалось получить просто с использованием значений с предыдущего слоя:

Ап+1 ^ дп Дп+1 г^ дп ип+1 ~ ип ип+1 ~ и и А0 ~ А0 , АН +1 ~ АМ +1, и0 ~ и0 , +1 ~ иМ +1.

Для проведения численных расчетов использовались следующие параметры: Ь = 10 Т = 1, ((г) = вт(3г) А0 = п, к = 0.9 и е = 0.01. Сетка по времени содержала 2000 узлов, сетки по пространству задавались с шагами Н = 1/(1.5)в, й = 1, 5. Норма погрешности вычислялась как сеточная Ь2 Ь = Т

Порядок сходимости оценивался посредством формулы Рунге для главного члена погрешности.

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

5.2.2 Задача для вязкой модели

В безразмерных переменных система (2) записывается следующим образом:

дА д(Аи) ди ди к дА и

ж + -V = 0 ж + иж + да = —еА, (22)

где е = КТС/АЛ.

Пусть начальные условия имеют следующий вид:

А(0, г) = А0 + еф1(г) + е2^(г) + ..., и(0, г) = и0 + е£(г) + е^(г) + ...,

где фг(г) и (г(г) обладают теми же свойствами, что в случае задачи для невязкой модели, а

А0

и

и0

считаются известными. Представим решение задачи

е

А(Ь, г) = А0 + еА1(Ь, г) + ..., и(Ь, г) = и + еих (Ь, г) + .... (23)

А0 = А0 и0 = и0

для нахождения первого приближения получим задачу вида:

^ + Сю+ Л,^ = 0, § + + и^ = — £, (24)

Ь г г Ь А0 г г А0

А1(0,г ) = ф1(г), и (0,г ) = (1 (г). (25)

1од2(И)

-0.5

Рис, 1: Графики зависимости логарифма Ь2-нормы погрешности от логарифма шага для задачи дня повязкой модели: 1 — схема расщепления; 2 — схема Лакса — Вендроффа; 3 — схема МиБСЬ; 4 — прямая с наклоном, равным двум; 5 — прямая с наклоном, равным четырем

Систему (24) можно записать в виде:

д и м д и

31 дг

(26)

где и = (Л\, и\)т, g = (0, — и0/А0)Т, а матрица М представляется в виде:

\Т

М =

Щ Ао

ио

Система (26) является строго гиперболической и ее матрицу можно представить в следующем виде М = ИЛЬ, где И и Ь есть матрицы правых и левых собственных векторов, а Л является диагональной матрицей из собственных значений. После умножения (26) слева на Ь получим систему относительно г) = Ьи^, г):

£+л 3;=* ^

где О = Начальное условие для w получается из начального условия и(0, г) = и0(г) следующим обр азом: w(0, г) = Ьи0(г).

Система (27) состоит из двух не зависящих друг от друга линейных уравнений переноса вида:

I+%=^

где О является константой. Решение начальной задачи с условием с(0,г) = с0(г) для этого уравнения имеет следующий вид:

с(Ь,г) = ОЬ + у(Ь,г),

где у(Ь, г) = с0(г — ЛЬ) является решением задачи с тем же начальным условием для уравнения

Зу Зу Ь г

После получения ^ решение и находится как и = Rw. Приме-

няя описанный подход, получим решение (24)—(25):

А1(Ь, г) = —:—1 (Ь22Ьпф1(г — Л1Ь) + ^2^2(1 (г — ЛхЬ)—

Ь11Ь22 — Ь12 Ь21 4

Ь12Ь

12^2101 (z - \2t) - L12L22Z1 (z - \2t)) ,

Ui(t, z) = - Uo t - ——1 fL2iLii0i (z - Ait) + L21 Li2Zi(z - Ait)-

Ao LiiL22 - Li2L2^

-LiiL2^i(z - A2t) - LiiL22Zi(z - A2t^ ,

где Lij представляют собой компоненты матрицы L.

Для получения численного решения на концах пространственного промежутка ставились условия совместности для случая вязкой модели:

TU + A±(U) Ж ,

в дополнение к которым ставились условия первого рода для скорости, вычисляемые по аналитическому решению при значениях z = ±L.

Численные расчеты производились для U0 = 1 при тех же значениях параметров, что и в случае невязкой модели на временной сетке из 5000 узлов и пространственных сетках с шагами h = 1/(1.25)s, s = 2, 7. На рис. 2 представлены графики зависимости логарифма нормы погрешности от логарифма шага. Как можно видеть, схема расщепления позволяет получать результаты с четвертым порядком, при этом значения погрешности в случае такой схемы меньше, чем для явных схем второго порядка.

1±(U^ — + А±(U) — - f) = 0, (29)

Ш

та

О

1од2(И)

-0.5

Рис, 2: Графики зависимости логарифма Ь2-нормы погрешности от логарифма шага для задачи дня вязкой модели: 1 — схема расщепления; 2 — схема Лакса — Вепдроффа; 3 — схема МиБСЬ; 4 — прямая с наклоном, равным двум; 5 — прямая с наклоном, равным четырем

5.3 Задача о течении в одиночном сосуде

Рассмотрим задачу о численном моделировании течения в одиночном сосуде, поставленную в [41]. В этой задаче рассматривается идеализированный сосуд длины Ь = 100 см, остальные параметры имеют следующие значения: ^ = 0.15 см, Е = 4• 106 дин/см2, р = 4• 10-2 дин/см2 •. Значение а берется равным 11/10, значение А^ равно п см2.

Начальные условия задавались как и(0, г) = 0 А(0,х) = Арасчеты производились на промежутке [0,0.5] с.

г=0

и (г, 0) = 1ехр(-104(£ - 0.05)2), п

в дополнение к которому ставилось условие совместности (29) для случая А_.

Л+

жающее условие, отвечающее характеристике, входящей в вычислительную область [43]:

1—(и)( — 0 = а (30)

В случае явных схем шаг ДЬ выбирался из (13) при псеь = 0.95 для схемы Лакса — Вендроффа и исрь = 0.45 для схемы МиЯСЬ. В случае схемы расщепления значение шага специально бралось больше, чем для неявных схем, но таким, чтобы не развивались численные неустойчивости, возникающие из-за наличия экстраполяции граничных значений на этапе предиктора. Расчеты производились па сетках с N = 500, 1000, 1500 и 2000 узлами. В табл. 1 представлено время расчетов и число шагов (в скобках) для рассмотренных схем. Как можно видеть, для случаев рассмотренных сеток наименьшее время расчетов и число шагов характерны для схемы расщепления.

Таблица 1. Время проведения расчетов в секундах и число шагов (в скобках) для разностных

схем в задаче о течении в одиночном сосуде

Схема N = 500 N = 1000 N = 1500 N = 2000

Лаке — Вендрофф 0.796 3.02 7.12 14.9

(441) (882) (1323) (1764)

MUSCL 3.21 11.5 26.1 44.8

(979) (1959) (2977) (4147)

Расщепление по физ. процессам 0.48 1.61 4.37 9.68

(270) (520) (770) (1030)

5.4 Задача о течении в сосуде с ветвлением

Рассмотрим систему, состоящую из основного сосуда (I) и двух дочерних сосудов (II) и (III) (рис. 3). Для проведения численных расчетов использовались параметры из работы [44]. Длина каждого сосуда бралась равной 20 см, начальный диаметр основного сосуда брался равным 1 см, диаметры дочерних сосудов выбирались равными 1/л/6 см. Плотность р бралась равной 1 г/см3. Значен ие Ad вычислялось как nR2, где R вычислялось через начальный диаметр. Параметр ß/Ad брался равным 32497 г/(см2с2) для основного сосуда и 79602 г/(см2с2) для дочерних сосудов.

Начальные скорости полагались равными нулю, начальные площади за-Ad

модели (3). На входе сосуда I значение U задавалось таким образом [44]:

U(t, 0) = Uo exp (-C(t - to)2) ,

I

Рис, 3: Схема модели сосуда с ветвлением

где U0 = 1 см/с, C = 5000 с-2, t0 = 0.05 с. В дополнение к этому условию задавалось условие совместности для характеристики, отвечающей А_. На выходах сосудов II и III задавались неотражающие условия и условия совместности.

В точке ветвления ставилось условие сохранения расхода:

Ui Ai = Un An + UIH Ain,

условие непрерывности давления:

Piii (Ai ii) + рЩ1 = Pj (Aj) + pU, Pji (Aji) + pUj = Pj (Aj) + pU,

и три условия совместности:

MUi) (^ + A+(Ui)^) = 0,

l_(Uii) (^ + A_(Uii)) = 0,

l_(Uiii) (ЩИ1 + A_(Uiii)d-Ull) = 0.

Производные, входящие в условия совместности, аппроксимировались посредством конечных разностей (20) (21). После их дискретизации для вы-

AU

алгераических уравнений. Для ее решения использовался метод Ньютона Рафсона. В качестве начального приближения использовались значения для момента t = tn. Такой выбор начального приближения обеспечивал быструю сходимость итерационного процесса 2-3 итерации при допустимой точности 10_6.

Расчеты производились на промежутке [0,0.4] с. В случае явных схем шаг по времени выбирался из условия устойчивости (13). Значение ucfl для

0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5

t, s t, s

Рис, 4: Графики Qu и Рц в средней точке дочернего сосуда, полученные для сетки из 1000 узлов: 1 — схема Лакса — Вепдроффа; 2 — схема MUSCL; 3 — схема расщепления

схемы MUSCL бралось равным 0.3, поскольку при больших значениях наблюдалось развитие численных неустойчивостей. В случае схемы Лакса Вендроффа бралось значение nCFL = 0.95. На рис. 4 представлены графики численного решения в центральной точке дочернего сосуда, полученного для сетки с 1000 узлами. Отметим, что графики давления аналогичны графикам, представленным на fig. 4 в [44]. В табл. 2 представлены время и число шагов для всех рассмотренных схем. Как можно видеть, при использовании схемы расщепления для невязкой модели расчеты можно производить за меньшее время, чем в случаях явных схем.

Таблица 2. Время проведения расчетов в секундах и число шах'ов (в скобках) для разностных

схем в задаче о течении в сосуде с ветвлением

Схема N = 500 N = 1000 N = 1500 N = 2000

Лаке Вендрофф 7.415 (1231) 28.44 (2463) 60.92 (3696) 115.6 (4929)

MUSCL 35.31 (3920) 138.4 (7847) 303.3 (11847) 544.3 (15921)

Расщепление по физ. процессам 2.171 (800) 11.31 (1500) 25.69 (2200) 68.37 (2900)

5.5 Задача о течении в модельной сосудистой системе

В работах [19, 20, 29] рассмотрена модельная сосудистая система (т.н. fractal-tree arterial system). Система представляет собой модель асимметричной артериальной системы, образующей дерево с элементарными ветвлениями. В каждой точке ветвления сосуд разделяется на два дочерних сосуда ¿и v, ра-

диусы которых вычисляются через радиус сосуда-родителя Rp так Rb = iRp и Rv = vRp. В настоящей работе, по аналогии с [29], использовались значения 1 = 0.9 и v = 0.6. Длина каждого сосуда вычислялась через радиус как L = 50R. Значения коэффициента ß/Ad определялись через значение радиуса по формуле:

= 4 (ki exp(k2R) + кз), Ad 3Jn

где к1 = 2 • 107 дин/см4, к2 = -22.53 см-1, к3 = 8.65 • 105 дин/см4. При численном моделировании предполагалось, что радиус входного сосуда равен 1 см. На входе в этот сосуд ставилось условие вида [29]:

и1(*, 0) = Пщ тах(^0,вт^п^т)) ,

где Тнн = 1 с и Ян = 50 см3/с. В дополнение к этому условию ставилось условие совместности вида (29). На выходе из каждого терминального сосуда ставилось аналогичное условие и неотражающее условие.

Рис. 5: Модель сосудистой системы из 7 сегментов

Для проведения численных расчетов бралось значение динамической вязкости, равное 0.035 с^дин/см2 [45]. Значение р бралось равным 1.06 г/см3 [41]. Коэффициент а брался равным 8/7 для того, чтобы воспроизвести приплюснутый профиль скорости, характерный для крови. При проведении расчетов рассматривалась модельная система из 7 сосудов, схема которой представлена па рис. 5. Расчеты производились па временном промежутке [0, 5]

с. При расчетах для схемы Лакса — Веидроффа использовалось значение nCFL = 0.95, для схемы MUSCL nCFL бралось равным 0.45. В табл. 3 представлены значения времени расчетов и число шагов за один период, равный 1 с для случаев трех разных сеток по пространству. Как можно видеть, для всех случаев наилучшие результаты характерны для случая схемы расщепления.

Таблица 3. Время проведения расчетов в секундах и число шагов (в скобках) за один период для разностных схем в задаче о течении в системе из семи сосудов

Схема N = 50 N = 100 N = 150

Лаке — Вендрофф 1.284 4.598 9.714

(796) (1609) (2422)

MUSCL 3.691 15.36 30.62

(1662) (3309) (4951)

Расщепление по физ. процессам 0.331 0.689 1.328

(470) (1100) (1900)

6 Заключение

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

Отметим следующие направления дальнейших исследований в рамках схем расщепления:

1. Проведение сравнения производительности схем при различных подходах к экстраполяции граничных условий, применяемых на этапе предиктора.

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

3. Разработка схем расщепления на основе других подходов к дискретизации по пространственной переменной. Большие перспективы в этой обла-

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

Список литературы

[1] Quarteroni A., Manzoni A., Vergara С. The cardiovascular system: mathematical modelling, numerical algorithms and clinical applications // Acta Numerica. 2017, vol. 26. P. 365-590.

[2] Audebert C., Bucur C., Bekheit M., Vibert E., Vignon-Clementel I., Gerbeau J. Kinetic scheme for arterial and venous blood flow, and application to partial hepatectomy modeling // Computer Methods in Applied Mechanics and Engineering, 2017, vol. 314. P. 102-125.

[3] Marchandise E., Willemet M., Lacroix V. A numerical hemodynamic tool for predictive vascular surgery // Medical Engineering and Physics, 2009, vol. 31. P. 131-144.

[4] Того E. F. Brain venous haemodynamics, neurological diseases and mathematical modelling: a review // Applied Mathematics and Computation, 2015, vol. 272. P. 542-579.

[5] Dobroserdova Т., Liang F., Panasenko G., Vassilevski Y. Multiscale models of blood flow in the compliant aortic bifurcation // Applied Mathematics Letters, 2019, vol. 93. P. 98-104.

[6] Xiao N., Alastruey J., Figueroa C. A systematic comparison between 1-D and 3-D hemodynamics in compliant arterial models // International Journal for Numerical Methods in Biomedical Engineering, 2014, vol. 30. P. 203-231.

[7] Delestre O., Lagree P.-Y. A 'well-balanced' finite volume scheme for blood flow simulation // International Journal for Numerical Methods in Fluids, 2013, vol. 72. P. 177-205.

[8] Delestre O., Ghigo A. R., Fullana J.-M., Lagree P.-Y. A shallow water with variable pressure model for blood flow simulation // Networks and Heterogeneous Media, 2016, vol. 11, no. 1. P. 69-87.

[9] Ghigo A. R., Delestre O., Fullana J.-M., Lagree P.-Y. Low-Shapiro hydrostatic reconstruction technique for blood flow simulation in large

arteries with varying geometrical and mechanical properties // Journal of Computational Physics, 2017, vol. 331. P. 108-136.

[10] Cavallini N., Caleffi V., Coscia V. Finite volume and WENO scheme in one-dimensional vascular system modelling // Computers and Mathematics with Applications, 2008, vol. 56, no. 9. P. 2382-2397.

[11] Huang P. G., Muller L. O. Simulation of one-dimensional blood flow in networks of human vessels using a novel TVD scheme // International Journal for Numerical Methods in Biomedical Engineering, 2015, vol. 31, no. 1. Art. no. e02701.

[12] Wang X., Fullana J.-M., Lagree P.-Y. Verification and comparison of four numerical schemes for a ID viscoelastic blood flow model // Computer Methods in Biomechanics and Biomedical Engineering, 2015, vol. 18. P. 1704 1725.

[13] Bessems D., Rutten M., van de Vosse F. A wave propagation model of blood flow in large vessels using an approximate velocity profile function // Journal of Fluid Mechanics, 2007, vol. 580. P. 145-168.

[14] Malossi A. C. I., Blanco P. J., Deparis S. A two-level time step technique for the partitioned solution of one-dimensional arterial networks // Computer Methods in Applied Mechanics and Engineering, 2012, vol. 237-240. P. 212— 226.

[15] Melicher V., Gajdosik V. A numerical solution of a one-dimensional blood flow model - moving grid approach // Journal of Computational and Applied Mathematics, 2008, vol. 215, no. 2. P. 512-520.

[16] Sherwin S. J., Formaggia L., Peiro J., Franke V. Computational modelling of ID blood flow with variable mechanical properties and its application to the simulation of wave propagation in the human arterial system // International Journal for Numerical Methods in Fluids, 2003, vol. 43. P. 673-700.

[17] Elad D., Katz D., Kimmel E., Einav S. Numerical schemes for unsteady fluid flow through collapsible tubes // Journal of Biomedical Engineering, 1991, vol. 13, no. 1. P. 10-18.

[18] Smith N. P., Pullan A. J., Hunter P. J. An anatomically based model of transient coronary blood flow in the heart // SIAM Journal on Applied Mathematics, 2002, vol. 62. P. 990-1018.

[19] Duanmu Z., Chen W., Gao H., Yang X., Luo X., Hill N. A. A one-dimensional hemodynamic model of the coronary arterial tree // Frontiers in Physiology, 2019, vol. 10. Art. no. 853.

[20] Olufsen M. S., Peskin C. S., Kim W. Y., Pedersen E. M., Nadim A., Larsen J. Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions // Annals of Biomedical Engineering, 2000, vol. 28. P. 1281-1299.

[21] Saito M., Ikenaga Y., Matsukawa M., Watanabe Y., Asada Т., Lagree P.Y. One-dimensional model for propagation of a pressure wave in a model of the human arterial network: comparison of theoretical and experimental results // Journal of Biomechanical Engineering, 2011, vol. 133. Art. no. 121005.

[22] Azer K., Peskin C. S. A one-dimensional model of blood flow in arteries with friction and convection based on the Womersley velocity profile // Cardiovascular Engineering, 2007, vol. 7. P. 51-73.

[23] Diem A. K., Bressloff N. M. VaMpy: A Python package to solve ID blood flow problems // Journal of Open Research Software, 2017, vol. 5. Art. no. 17.

[24] Huo Y., Kassab G. S. A hybrid one-dimensional/Womersley model of pulsatile blood flow in the entire coronary arterial tree // American Journal of Physiology-Heart and Circulatory Physiology, 2007, vol. 292, no. 6. P. H2623-H2633.

[25] Watanabe S. M., Blanco P. J., Feijoo R. A. Mathematical model of blood flow in an anatomically detailed arterial network of the arm // ESAIM: M2AN, 2013, vol. 47, no. 4. P. 961-985.

[26] Carson J., van Loon R. An implicit solver for ID arterial network models // International Journal for Numerical Methods in Biomedical Engineering, 2017, vol. 33, no. 7. Art. no. e2837.

[27] Kapo К., Педли Т., Шротер P., Сид У. Механика кровообращения. М.: Мир, 1981. 624 с.

[28] Formaggia L., Lamponi D., Quarteroni A. One-dimensional models for blood flow in arteries // Journal of Engineering Mathematics, 2003, vol. 47. P. 251276.

[29] Ghigo A. R., Lagree P.-Y., Fullana J.-M. A time-dependent non-Newtonian extension of a ID blood flow model // Journal of Non-Newtonian Fluid Mechanics, 2018, vol. 253. P. 36-49.

[30] Puelz C., Canic S., Riviere В., Rusin C. G. Comparison of reduced models for blood flow using Runge-Kutta discontinuous Galerkin methods // Applied Numerical Mathematics, 2017, vol. 115. P. 114-141.

[31] Krivovichev G. V. Computational analysis of one-dimensional models for simulation of blood flow in vascular networks // Journal of Computational Science, 2022, vol. 62. Art. no. 101705.

[32] Britton J., Xing Y. Well-balanced discontinuous Galerkin methods for the one-dimensional blood flow through arteries model with man-at-eternal-rest and living-man equilibria // Computers and Fluids, 2020, vol. 203. Art. no. 104493.

[33] Sheng W., Zheng Q., Zheng Y. The Riemann problem for a blood flow model in arteries // Communications in Computational Physics, 2020, vol. 27. P. 227-250.

[34] Того E. F., Siviglia A. Flow in collapsible tubes with discontinuous mechanical properties: Mathematical model and exact solutions / / Communications in Computational Physics, 2013, vol. 13. P. 361-385.

[35] Spiller С., Того E. F., Vazquez-Cendon M. E.,Contarino C. On the exact solution of the Riemann problem for blood flow in human veins, including collapse // Applied Mathematics and Computation, 2017, vol. 303. P. 178189.

[36] Яненко H. H. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967. 196 с.

[37] Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1978. 592 с.

[38] Richtmyer R. D., Morton К. W. Difference Methods for Initial-Value Problems. Florida, Krieger Publishing, 1967. 420 p.

[39] Самарский А. А., Гулин А. В. Устойчивость разностных схем. M.: Наука, 1973. 415 с.

[40] Бахвалов Н. С., Жидков Н. П.. Кобельков Г. М. Численные методы. М.: Лаборатория знаний, 2020. 636 с.

[41] Boileau E., Nithiarasu P., Blanco P. J., Muller L. О., Fossan F. E., Hellevik L.R., Donders W. P., Huberts W., Willemet M., Alastruey J. A benchmark study of numerical schemes for one-dimensional arterial blood flow modeling // International Journal for Numerical Methods in Biomedical Engineering, 2015, vol. 31. P. 1-33.

[42] Магомедов К. M., Холодов А. С. Сеточно-характеристические численные методы. М.: Юрайт, 2018. 314 с.

[43] Hedstrom G.W. Nonreflecting boundary conditions for nonlinear hyperbolic systems // Journal of Computational Physics, 1979, vol. 30, no. 2. P. 222-237.

[44] Xiu D., Sherwin S.J. Parametric uncertainty analysis of pulse wave propagation in a model of a human arterial network // Journal of Computational Physics, 2007, vol. 226, no. 2. P. 1385-1407.

[45] Razavi M.S., Shirani E. Development of a general methods for designing microvascular networks using distribution of wall shear stress // Journal of Biomechanics, 2013, vol. 46. P. 2303-2309.

Implicit finite-difference schemes for equations of one-dimensional

hemodynamics

G.V. Krivovichev, N.V. Egorov

Saint-Petersburg State University, Faculty of Applied Mathematics and Control

Processes

E-mail: g.krivovichev@spbu.ru, n.v.egorov@spbu.ru

Abstract. The paper is devoted to the construction and analysis of implicit finite-difference schemes for a system of one-dimensional equations of hemodynamics. The schemes are based on the use of finite differences, which approximate spatial derivative with the fourth order. The schemes are based on the splitting on physical processes. According to this approach, at one time step, two mechanical processes are considered: the deformation of the vessel filled with fluid and the fluid flow in the deformed vessel. This approach makes it possible to separately consider finite-difference schemes, which approximate governing equations. In the practical implementation of the proposed schemes, they are reduced to systems of linear algebraic equations with pentadiagonal matrices.

The stability analysis of constructed schemes is based on the von Neumann method and the principle of frozen coefficients. In the numerical solution of problems with known analytical solutions, it is demonstrated that the schemes lead to numerical solutions with a fourth-order convergence rate. In the computational experiments on simulation of blood flow in model vascular systems, it is demonstrated that the developed schemes make it possible to perform calculations in much less time than well-known explicit finite-difference and finite-volume schemes.

Keywords: one-dimensional models, hemodynamics, finite-difference schemes, splitting method, stability