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

Многочисленными отечественными и зарубежными исследованиями доказано негативное влияние неблагополучной санитарно-гигиенической обстановки на состояние здоровья человека, в том числе на дыхательную систему [1 –3] . Существующие методы медицинской диагностики [4, 5] предоставляют возможность определять состояние пациента на момент обследования, но они не предназначены для получения прогнозных оценок, в том числе влияния вредных факторов. Прогнозные оценки изменения состояния здоровья населения в зависимости от различных воздействий окружающей среды, которые не всегда удается воспроизвести в натурных испытаниях (например, при выбросах загрязняющих веществ с концентрацией, существенно превышающей допустимую), позволяет сделать математическое моделирование. На протяжении нескольких лет авторами статьи ведется работа по созданию математической модели дыхательной системы человека, предназначенной в том числе для расчета объема поступления ингаляционным путем загрязняющих химических веществ из атмосферного воздуха в организм человека и последующего прогнозирования возникновения патологий органов дыхания, обусловленных негативным воздействием на них [6] .

При исследовании течения воздушной смеси в дыхательной системе широкое распространение получили уравнения газовой динамики [7, 8] , которые с приемлемой точностью описывают процессы в верхних и крупных нижних воздухоносных путях человека. Нижние воздухоносные пути представляют собой разветвленную сеть каналов. При увеличении числа генераций ветвления дыхательных путей возникают сложности как с построением реальной геометрии отдельных каналов и альвеол (в легких взрослого человека содержится около 600–700 млн альвеол, а также соединяющих их каналов [9] ), так и с резко возрастающими требованиями к вычислительным мощностям. Другим подходом к изучению движения воздуха в легких человека является их рассмотрение в виде

Статья опубликована в открытом доступе по лицензии CC BY 4.0

насыщенной пористой среды [6, 10 –13] ; движение воздуха через деформируемую пористую среду моделируется с привлечением теории фильтрации [14, 15] .

В разрабатываемой авторами математической модели биомеханика дыхания человека считается совокупностью взаимосвязанных процессов: 1) течения воздуха в недеформируемой воздухопроводящей зоне (воздухоносных путях) (оно подробно обсуждается в предыдущих работах авторов [6, 16, 17] ); 2) течения воздуха в упруго-деформируемой респираторной зоне (легких человека) [12, 13] . Последнее служит предметом настоящего исследования.

В предыдущей статье авторов [13] приведена математическая формулировка модели течения воздуха в легких человека, аппроксимированных двухфазной насыщенной пористой средой, одна из фаз которой — деформируемый скелет среды (легочная ткань), вторая фаза — газ, заполняющий поровое пространство; представлены разрешающие соотношения, необходимые для решения задачи. Легкие человека в процессе дыхания испытывают большие градиенты перемещений (компоненты которых по результатам предварительных расчетов достигают в некоторых подобластях, особенно в примыкающих к диафрагме и грудной клетке, 90% при спокойном дыхании), вследствие чего рассматриваемая задача относится к классу геометрически нелинейных связанных проблем. Получить ее аналитическое решение для двухфазной пористой среды, испытывающей большие перемещения, и фильтрации газа в ней не представляется возможным. В связи с этим предлагается прибегнуть к численным методам, содержащим пошаговые процедуры.

Численное решение нелинейной подзадачи деформирования двухфазной среды легких предлагается осуществлять методом конечных элементов (МКЭ) [18] , подзадачу фильтрации (относительное течения воздуха в деформируемой пористой среде легких) решать методом конечных (контрольных) объемов (МКО), который хорошо зарекомендовал себя при реализации задач механики жидкости и газа [19, 20] . Достоинством МКО является выполнение законов сохранения как на глобальном, таки локальном уровнях. Для каждого контрольного объема (КО) гарантируется соблюдение баланса полного потока, что является дополнительным плюсом при рассмотрении задач гидроаэродинамики, уравнений переноса. В предложенном алгоритме метод конечных объемов участвует только в расчете просачивания воздуха (обмена массой между соседними КЭ), входные данные (величины давления в элементах) получены с применением МКЭ. При этом в данной работе использованы симплекс-элементы; давление во всех элементах в каждый момент времени однородно. При решении уравнения Дарси (описывающего фильтрацию газа) с помощью МКЭ потребовался бы повышенный порядок аппроксимации (для определения градиента давления). Это повлекло бы увеличение объема вычислительных ресурсов. Решение же уравнения Дарси МКЭ с привлечением симплекс-элементов нуждается в создании дополнительной процедуры, которая и реализуется в представленной работе МКО.

Исследование процесса течения воздуха в дыхательной системе человека производится численно для ее индивидуальной трехмерной геометрии, восстановленной на основе данных компьютерной томографии [13, 16, 17, 21] .

• 2.    Восстановление трехмерной формы легких, задание закона движения стенок

Реконструкция реальной трехмерной формы легких человека выполнена исходя из данных компьютерной томографии, взятых из открытой базы данных I-ELCAP (International Early Lung Cancer Action Program) [22] . Обработка исходных двумерных изображений (томографических снимков в формате DICOM, размещенных в файлах с расширением dcm), процедура сегментации и восстановление трехмерной формы осуществлены посредством программного продукта с открытым доступом ITK-SNAP [23] . Полученная после применения процедуры сегментации трехмерная форма легких представляет собой поверхностную сетку, информация о которой хранится как список составляющих ее треугольных граней и нормалей к ним в файле формата stl. Сетка обработана в графическом редакторе Blender. В результате обработки построены поверхности бронхов (5 поколений дерева бронхов) и внешние поверхности легких, а также выполнено сглаживание этих поверхностей. На рисунке 1 представлена трехмерная форма легких до и после сглаживания.

При обработке в редакторе Blender из внутренней области легких (Рис. 1б ) исключены внутрилегочные участки бронхов (используемые при исследовании течения воздуха в воздухопроводящей зоне [17] ). На основе поверхностной сетки легких в комплексе ANSYS (модуль ICEM CFD) сгенерирована объемная конечно-элементная сетка, аппроксимирующая деформируемую насыщенную пористую среду легких для изучения фильтрации воздуха. Эта сетка представлена на рисунке 2. Здесь показаны области и границы, учитываемые в разрабатываемой модели легких человека. Символом Ω обозначена внутренность области легких, рассматриваемая как двухфазная

Рис. 1. Восстановленная трехмерная форма легких до обработки ( а ), после обработки сглаживанием ( б )

( a )

пористая среда; Γ _L — стенки легких; Γ _AW/_L — границы выходов из бронхов, являющиеся входами, через которые в легкие воздух поступает во время вдоха и выходит из легких во время выдоха. В рассматриваемой геометрии присутствует 31 вход в легкие из системы бронхов (14 входов в левом легком и 17 в правом).

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

Легкие находятся в грудной клетке (каждое в своем плевральном мешке). Они соприкасаются снизу с главной дыхательной мышцей — диафрагмой, разделяющей грудную и брюшную полости, по сторонам — со стенкой грудной клетки [9] . Движение воздуха в легких в процессе дыхания (вдох, выдох) обусловлено разностью давления на входах бронхов в легкие и изменяющимся давлением воздуха в самих легких вследствие изменения их объема. Изменение объема легких осуществляется за счет грудного дыхания (расширения/сжатия грудной клетки) и диафрагмального (брюшного) дыхания (опускания/подъема диафрагмы). При этом наибольшее изменение формы (максимальные перемещения) реализуется за счет движения диафрагмы и обеспечивает вентиляцию воздуха в 2/3 объема легких [24, 25] .

Наиболее детальное описание движения стенок легких получено с использованием 4D-компьютерной

Рис. 2. Конечно-элементная сетка легких для исследуемой области, заполненной двухфазной пористой средой (вид спереди)

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

томографии [26 –29] . Название происходит от добавления времени (как четвертого измерения) к традиционной трехмерной томографии (4D = 3D+ 1). К сожалению, авторам не удалось найти в открытом доступе данные об изменении 3-мерной геометрии легких в разные моменты дыхательного цикла. Стенки легких в процессе дыхания имеют сложный характер движения, для описания которого необходимо создание подмодели костно-мышечной системы и диафрагмы, окружающих дыхательные органы человека. В связи с этим в данной работе предлагается в первом приближении использовать упрощенные граничные условия, имеющиеся в обзорах результатов медицинских исследований, обеспечивающие интегральные характеристики — изменение общего объема воздуха в легких, смещение диафрагмы и грудной клетки.

В качестве начального состояния дыхательного цикла рассмотрим момент «конец выдоха–переход к вдоху». Для этого момента полагается, что скорости движения частиц воздуха равны нулю, перепад давления в легких и воздухоносных путях (на входах в легкие) отсутствует, давление воздуха в легких аналогично атмосферному). В процессе дыхания легкие расширяются в горизонтальной (0x 1 x 2 ) и вертикальной (фронтальной) (0x 1 x 3 ) плоскостях (Рис. 3а ). При спокойном дыхании расширение (экскурсия) грудной клетки — разность между показателями ее окружности на вдохе и выдохе, составляет 2-3 см [30] , на максимальном вдохе и выдохе—7 ^ 8.5см [30] , 6 ^ 8 см [31] . В зависимости от вертикальной координаты (x ₃ ) экскурсия грудной клетки (и, соответственно, легких) различается. В области вершины легкого она наименьшая, по мере приближении к основанию легких степень их расширения увеличивается. На уровне подмышечной впадины периметр грудной клетки изменяется в среднем на 3.52 см у мужчин и на 3.19 см у женщин; на уровне 3–4-го межреберного пространства, соответственно полу, — на 4.02 см и на 3.71 см; в области нижнего края грудины (на уровне мечевидного отростка) — на 4.52 см и на 4.12 см [32] . Средние показатели периметра грудной клетки зависят от пола, возраста, конституции человека: у мужчин 18 лет это 90.20 ± 5.24 см (среднее ± стандартное отклонение — M ± SD), 26-30 лет — 92.24 ± 4.94, у женщин 18 лет — 84.08 ± 4.38 см, 26-30 лет — 85.32 ± 4.72 [33] . Исходя из проведенных оценок степень расширения грудной клетки (и легких соответственно) при спокойном дыхании оценивается в 2–4.5% от исходного состояния.

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

При спокойном дыхании диафрагма, по разным источникам, смещается на 1 см [34] или 0.5 ^ 2 см [35] , при форсированном дыхании амплитуда ее движения может достигать 10 см [34] или 7 ^ 8 см [35] . В соответствии с данными МРТ из [36] , средняя амплитуда движения диафрагмы при спокойном дыхании составила 27.3 ± 10.2 мм (M ± SD). В работе [37] приведены результаты измерений, полученные с использованием цифровой рентгенографии, согласно которым амплитуда перемещения правого купола диафрагмы достигала 25.1 ± 6.8 мм, левого купола — 30.6 ± 8.2 мм.

В работах [10, 11, 38], посвященных численному моделированию процессов в легких, рассмотрен синусоидальный закон движения диафрагмы, в [10, 11] максимальное смещение в основании легких принималось равным 1.5 см. В [39] амплитуда смещения опорных точек, лежащих в основании правого легкого, составляла 12.7±6.1 мм (M±SD), в основании левого легкого — 16.5±4.5 мм; максимальное смещение достигало 24 мм.

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

При взятом в качестве начального отсчетного состояния дыхательного цикла моменте «конец выдоха–переход к вдоху» принимается периодический закон движения стенок; период дыхания при обычном спокойном дыхании равен t * =4 с [34] .

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

V s = (2n/t * )(5 | 1 | /2) 1 n sin(2nt/t * ), (1)

где V s — вектор их скорости, l—вектор длиной 1 1 | , соединяющий «центр» и точку поверхности легкого в начальный момент времени, l n — нормированный вектор, направленный вдоль l , δ — отношение максимального смещения (в направлении l n ) точки на поверхности легкого от значения | l | в начальном положении, δ | l | — максимальное смещение точки поверхности легкого от начальной точки в процессе дыхания, 2n/t * — циклическая частота.

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

V d = V s + V d . (2)

Здесь вектор V d описывает скорость неравномерного смещения по вертикальной координате, [м/с]. Вектор V d имеет только вертикальную ненулевую компоненту v * 3 (v * 1 = v * ₂ = 0); для каждого легкого она своя и определяется из соотношения:

* = (-(d ■ l(-Cn/t*)si„(2n(/(*), (3) I (max(d)) I где C — смещение точек вершины купола основания легкого без учета их смещения по вертикальной координате за счет первого слагаемого из правой части (2) за время t*/2 = 2 с (при расчетах в качестве C используются величины: CL для левого легкого, CR для правого), d — расстояние от точки в основании легкого до ближайшей точки контура, одновременно принадлежащего как границе основания, так и боковой стенке соответствующего легкого в начальный момент времени.

Максимальная удаленность от контура основания легкого (max(d)) имеет место в точке, являющейся вершиной купола диафрагмы (прилегающей к основанию легкого). Согласно соотношению (3) при d = max(d) скорость смещения точек во фронтальной плоскости максимальна, за 2 с вдоха вершина купола диафрагмы смещается на C . Вертикальная компонента скоростей точек контура основания легкого (d = 0) равна нулю; скорость смещения точек основания легкого определяется только за счет первого слагаемого в правой части соотношения (2) .

Результаты выполненного обзора медицинских исследований и осуществленных предварительных расчетов позволили идентифицировать параметры, входящие в закон движения стенок легких, используемый для задания граничных условий. Получены следующие значения: 5 = 0.0191, C L =0.013 м, C R = 0.012 м. При этих параметрах максимальное смещение по вертикальной координате основания левого легкого достигает 0.0155 м, правого легкого — 0.0146 м; объем газовой фазы в легких изменяется на 0.79 л (0.00079 м ³ ). Приведенные параметры соответствуют известным физиологическим данным для спокойного дыхания взрослого человека.

• 3.    Алгоритм решения связанной задачи

Связанная задача фильтрации воздуха в упруго-деформируемой пористой среде является нелинейной. Ее реализация осуществляется в виде пошаговой (по времени) процедуры при достаточно малых временных шагах.

Последовательно решаются подзадачи расчета параметров деформирования двухфазной пористой среды (с использованием МКЭ) и относительного течения воздуха в пористой среде (с использованием МКО).

На текущем ( n - м) шаге рассматривается отрезок времени [t n - i ,t n ]. Решение в скоростях отыскивается при времени t _{n -1} и граничных условиях (в скоростях) для этого же момента t _{n -1} при значениях переменных с конца предыдущего — (п — 1) - го, шага. Обсуждаемая задача сводится к двум подзадачам: подзадаче деформирования двухфазной среды в скоростях и подзадаче фильтрации (относительного течения воздуха в деформируемой пористой среде легких) с использованием результатов расчета давления, полученных при решении первой подзадачи. По определенным в момент t _{n -1} скоростям изменения искомых параметров путем интегрирования вычисляются значения параметров в момент t n ; для интегрирования применяется явная схема Эйлера. Процессы деформирования пористой среды и фильтрации воздуха в ней различаются по интенсивности и скорости распространения в пространстве, в связи с чем при решении подзадачи фильтрации полагается переменный (более мелкий) временной шаг; на произвольном n-м шаге решения подзадачи деформирования производится K _n >  1 решений подзадачи фильтрации.

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

Скорости частиц двухфазной среды находятся через скорости движения частиц каждой из фаз как среднемассовая скорость смеси. В силу того, что масса газа, содержащегося в объеме двухфазной среды легких, пренебрежимо мала по сравнению с массой их твердого каркаса, скорость двухфазной среды можно считать совпадающей со скоростью твердой фазы ( v = v s ). Скорость движения частиц воздуха в рассматриваемой точке можно выразить в виде суммы скоростей «переносного» движения (смещения частиц вместе с деформируемой твердой фазой во введенной лабораторной системе координат) и «относительного» движения, обусловленного просачиванием газа через твердый каркас.

Разрешающие соотношения МКЭ для подзадачи деформирован ия д вухфазной пористой среды легких человека, аппроксимированной сеткой из m тетраэдральных элементов (m =1 ,M), получены в [13] (см. формулу (29)). В этом соотношении суммирование по элементам осуществляется для компонент матриц с их глобальными номерами:

17_R^(m)q , _R(^m)r,,A, _R(^m)q । _R(^m)r,,A, _R^(m)q 3 ZB kke + B jia ^ f s B jie + B jia IB^B je

^{( m ) r 2} ' rmmqq, o^{(m)r 2} x o^(m)q 1 J p(m) L(m)e , B jja g ^^Y s ^B kke + B jia 3 ^ j B kke Jd^ (v q +

+

_R⁽m)r 1ft _R^(m)q B jia 2 S kj B ike

,⁽m)r 1^ Tm(^m)q_ _R^(m)r ¹ q Tm(^m)q

’ jia 2 S ^kj B kie    B jia 2 S ik B kje ⁺

1 d(m)r 1 q D(m)q D(m)r V D(m)q 1лo(m) СЛ1(т)в_|_ + Bjia ^ Sik Bjke ~Bjia ^kj Bike Jd^v

{ Mr

X X I ^P E- \ Го^(т)г /О     C D^(m)q  TD^(m)r ,Q      Q R^(m)ql z7 o(^m) l^(^m)e_|_

I \ o I LBjia QiplkSPj Bike   Bjia QpjlkSip Bike J d^    fvq m=b(m) PP'

M Γ

X                ( m ) r          ( m ) q      ( m ) q

⁺ L /   N ia   —T i B ^kkU +B lke n k n l T i

(m=1

(m)

Γ AW/L

\^L]mr==| E / Nm[T^rm/A.

J            m = l (m)>

Γ AW/L

В (4) приняты обозначения: индексы в виде греческих букв α и β — номера узлов, в виде латинских букв — координатные оси; Nim — функции формы (Na) (re) = 0 при a = в; Nim(re) = ^i; i,j = 1,3); Bjm) (r) = dNjem' (r)/dxi — производные функций формы; V^m)13 — вектор узловых скоростей; компоненты скоростей перемещений внутри элемента и компоненты градиента вектора скорости со скоростям перемещений в узлах связаны формулами i(m) = Nim4 (r)v(m)e, vj™) = dvjm / dxl = В^в4 (r^v^m)13; p/p — отношение плотностей двухфазной среды в текущей и отсчетной конфигурациях; Z — переменная, введенная в описание связи скорости изменения среднего напряжения в терминах тензора напряжений Кирхгоффа двухфазной среды (Σav) со скоростью изменения объема двухфазной среды (первым инвариантом тензора деформации скорости D) с учетом взаимодействия газа и легочной ткани: Sav = Ii(D)Z [12, 13,40]; j — второй параметр Ламе твердой фазы; is — объемная доля твердой фазы; Σij — компоненты взвешенного тензора напряжений Кирхгоффа двухфазной среды лл

i

л_

в актуальной конфигурации; S kj — компоненты девиатора тензора напряжений Кирхгоффа твердой фазы; Q iplk — компоненты тензора четвертого ранга Q , входящего в определение логарифмического спина [41, 42] ; П ^(m) — подобласть, занимаемая m-м конечным элементом двухфазной пористой среды легких; r AW/L — границы выходов из воздухоносных путей, которые одновременно служат входами для воздуха в легкие человека (Рис. 2) ; T i — компоненты вектора поверхностных сил (вектора напряжений); n k — компоненты вектора нормали.

Соотношение, являющееся основным в МКО для определения плотности газовой фазы в элементе легких — КО (его вывод представлен в [13] , см. выражение (32)), имеет вид:

(              a     (n - 1)    (n - 1) (n - 1) \

, (n-1) k(H) pfA(a) -pfC   eCA^X   (n-1)s(n-1)       At      ,  (n-1)™ pf -f  •     |e(n-1)|     |e

В качестве КО используются те же КЭ, которые применяются для решения подзадачи деформации пористой среды легких. На основе соотношения (5) определяется изменение плотности газовой фазы в КО (КЭ) за счет изменения их объема вследствие поступления/оттока массы (между соседними элементами или через границы (грани) Γ_AW/_L, являющиеся выходами из воздухоносных путей и одновременно входами в легкие), обусловленное градиентом давления газовой фазы.

Подзадачи деформирования двухфазной пористой среды легких и фильтрации воздуха в пористой среде взаимосвязаны между собой. Поля напряжений и перемещений, установленные из решения подзадачи деформирования, применяются далее для решения подзадачи фильтрации. Давление воздушной фазы зависит от изменения двух составляющих: объема легких в процессе деформирования и значения плотности воздушной фазы (производные по времени давления и плотности прямо пропорциональны). Переопределяемое при решении подзадачи фильтрации давление необходимо для решения подзадачи деформирования пористой среды на следующем временном шаге (входит в соотношение для расчета переменной Z, описывающей связь скорости изменения среднего напряжения двухфазной среды легких со скоростью изменения объема двухфазной среды, учитывающей взаимодействие газа и легочной ткани [13]). Таким образом, в процессе решения выполняется обмен данными между подзадачами.

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

Этап I. Постановка подзадачи деформирования двухфазной среды осуществлена в скоростной форме, наиболее подходящей для геометрически нелинейных задач [43], поэтому в момент времени tn-1 n-го шага решается подзадача деформирования двухфазной пористой среды в скоростях с использованием граничных условий в скоростях для момента времени tn-1 и значений переменных с конца (п — 1)-го шага по времени.

Определяются вектор узловых скоростей (v⁽ⁿ 1)), скорости изменения компонент тензора напряжения (S⁽   )),

тензора Альманси (A⁽   )) и давления газовой фазы за счет изменения объема (pfn 1)) (верхний индекс (п— 1)

обозначает их соответствие моменту времени t_n-i).

Этап II. Решается подзадача фильтрации (для объема V^(п 1)), находятся скорости изменения характеристик

газовой фазы (pf 1), pf 1)) на момент t_n-1и значения параметров (pfⁿ⁾

(п) ~(n)\                  , mf , pf ) на момент tn .

Этап III. Интегрируются скорости изменения искомых переменных, отвечающие моменту времени t_n-1,

и вычисляются их значения на момент времени tn для использования на следующем временном шаге (следует заметить, что перед началом первого шага по времени в цикле по времени в подзадаче деформирования задаются условия на момент времени t0, являющиеся начальными на 1-м временном шаге).

Начальный момент времени t₀. В качестве этого времени, как отмечено ранее, выбран момент перехода от выдоха к вдоху. При этом легочная ткань находится в естественном ненапряженном состоянии, давление газовой фазы в легких и на входе одинаковое и равно атмосферному: pf⁰⁾= patm; шаровая составляющая и компоненты девиаторной части (S) тензора напряжений твердой фазы являются нулевыми; взвешенный тензор напряжений Кирхгоффа двухфазной среды определяется шаровой частью, основной вклад в которую вносит давление газа, находящегося в легких (S⁽⁰⁾= Z^I); среднее напряжение взвешенного тензора напряжений Кирхгоффа есть среднее напряжение тензора Коши (Д® = j^V)); среднее напряжение двухфазной среды следует из соотношения, полученного в результате решения вспомогательной задачи о всестороннем сжатии/растяжении представительного объема двухфазной среды в форме шара, заполненного воздухом [40]:

■             . у \ ab(V—Vf) (\/VV  pv] — CVfPf

(AVf + BV И Л     J

Здесь a_av — среднее напряжение в терминах тензора напряжений Коши двухфазной среды; V, V — объем двухфазной среды в отсчетной и текущей конфигурациях; Vf — объем газовой фазы в текущей конфигурации; A = (2р+3а), B = 4р, C = (6р+3а), а, р — параметры Ламе твердой фазы.

Далее, через координаты узлов определяются функции формы [N]⁽⁰⁾и их производные [B]⁽⁰⁾, а также объем V⁽⁰⁾КЭ. Согласно известным экспериментальным данным объемная доля воздуха в легких варьируется от 81 до 87.4% [44–47]; величина пористости (доля газовой фазы) в начальный момент времени принимается равной 0.85 (Yf⁰⁾), объемная доля легочной ткани, соответственно, 0.15 (7⁽⁰⁾). На основе информации об объеме элемента и доле газовой фазы находится объем газовой фазы в отсчетной конфигурации, а также объем твердой фазы (легочной ткани), который в процессе дыхательного цикла меняется несущественно [48] и потому полагается постоянным на всех временных шагах.

По заданным параметрам упругости легочной ткани, пористости материала при атмосферном давлении газовой фазы вычисляется значение Z⁽⁰⁾. В начальном состоянии компоненты тензора деформации Альманси (A⁽⁰⁾) нулевые; собственные значения тензоров искажений, характеризующих отношение длин материальных отрезков в актуальной и отсчетной конфигурациях, направленных вдоль главных векторов левой и правой мер искажений pi и pi [43], равняются 1. Компоненты тензора спина, а также тензора четвертого ранга Q, отыскиваемые с использованием собственных значений тензоров искажений, в начальный момент времени являются нулевыми.

Рассмотрим детально реализацию алгоритма на произвольном n-м шаге по времени. Полагается, что все параметры (Z, S, S, Q, [B], объемы КЭ (двухфазной среды), объемы и объемные доли фаз) на момент начала данного шага t_n-1(конца предыдущего — (п — 1)-го шага) определены (на 1-м временном шаге они известны из начальных условий).

На Этапе I алгоритма решается подзадача деформирования двухфазной среды в скоростях. На основе разрешающего соотношения МКЭ (см. (4)) по значениям переменных (Z, S, S, Q, [B]) с конца предыдущего (п — 1)-го шага по времени (для момента времени t_n-1) и граничным условиям на начало текущего n-го шага (для момента времени t_n-1) формируется система линейных алгебраических уравнений (СЛАУ). В результате ее решения устанавливаются значения компонент вектора среднемассовой скорости двухфазной среды v^(n-1)в узлах КЭ сетки (этот вектор, в силу относительной малости удельного веса воздуха, с высокой степенью точности равен вектору скорости частиц твердой фазы; при этом вклад от скорости просачивания воздуха является малым). При реализации разрешающего соотношения (4) в части выполнения силовых граничных условий последний член левой части и правая часть, содержащие интегралы по поверхности, принимаются равными нулю.

Рис. 4. Схема алгоритма решения связанной задачи течения воздуха в упруго-деформируемой пористой среде

По вычисленным компонентам вектора скорости перемещений двухфазной среды v⁽ⁿ 1) и функций формы [B]⁽ⁿ 1) определяются компоненты градиента вектора скорости перемещений двухфазной среды Vv^(n-1), вычисляются компоненты тензора деформации скорости D^(n-1)= (vv⁽ⁿ 1) +v⁽ⁿ ^V)/2, его первый инвариант I1 (D^(n-1)) и компоненты девиатора d^(n-1)= D^(n-1)— (1/3)I1 ^D^(n-1)j I, компоненты тензора вихря W^(n-1)= (v^(n-1)V — V v^(n-1)j.2.

Материальная производная тензора напряжений Кирхгоффа с учетом переменных с конца предыдущего — (n— 1)-го, шага по времени (на момент времени t_n-1) рассчитывается на основе соотношения:

± (n 1) = I1(D(n-1))Z (n-1)I+2Md(n-1)7(n-1) ^^EnT1) •S(n-1) — S(n-1) ^nT1), где D — тензор деформации скорости двухфазной среды, I — единичный тензор, d — девиатор тензора деформации скорости, 7s — объемная доля твердой фазы, Qlog — логарифмический спин.

С использованием градиента вектора скоростей перемещений Vv^(n-1)и значения тензора деформации A^(n-1)с предыдущего шага по времени определяется скорость изменения тензора Альманси [43]:

A^(n-1)= 1(7v⁽n^-1)+Vv^T(n^-1)) — (Vv⁽n^-1)•A^(n-1)+ A^(n-1)•Vv^T(n^-1)).

Выше сказано, что давление воздуха зависит от двух составляющих: 1) от изменения объема (результата решения подзадачи деформирования); 2) от перетекания воздуха между КО (результата решения подзадачи фильтрации).

Скорость изменения давления газовой фазы за счет изменения объема (обозначаемая как pf 1)) становится известной на Этапе I алгоритма в результате вычисления первого инварианта тензора деформации скорости. Для установления связи pf 1) со скоростью изменения объема используются следующие рассуждения. В процессе движения по воздухоносным путям (включающим в себя носовую полость, глотку, гортань, трахею и 5 поколений дерева бронхов) вдыхаемый воздух нагревается/охлаждается до температуры тела. Полагается, что фильтрация воздуха в пористой среде легких происходит при постоянной температуре. На каждом малом промежутке времени массу газа, находящегося в альвеолах, можно считать постоянной; при изотермических условиях pf Vf = const (или pf pf1 = const). Соотношение для изменения давления газовой фазы в зависимости от изменения объема газовой фазы (или от изменения плотности), вытекающее из закона Бойля–Мариотта, записывается в виде pf /pf = —Vf Vf, или pf /pf = pf /Pf.

Первый инвариант тензора деформации скорости определяет скорость изменения объема двухфазной среды: I1 (D) = dV / dV. Изменение объема двухфазной среды обусловлено изменением объема газовой фазы (изменением объема твердой фазы можно пренебречь):

ЛЛ

VV

I1(D) = V « V.

Из соотношения (7) следует выражение:

Vf /Vf = I1(D) (V /Vf) = I1(D)/7f.

Тогда с использованием (6) и (8) скорость изменения давления газовой фазы за счет изменения объема (pf 1)) определяется как pf /pf = I1 (D)/Yf или pfn-1) = —pfn-1)I1(D(n-1))/7fn-1).

На Этапе II при реализации подзадачи фильтрации воздуха временной интервал от конца предыдущего/начала текущего среза по времени (момента tn-1) до конца текущего (n-го) шага (на момент tn), равный одному шагу по времени в решении подзадачи деформирования, разбивается на K временных шагов. Введем индекс k для обозначения номера шага по времени в подзадаче фильтрации (к =1 ,К). При k = 1 момент времени соответствует началу шага по времени tn-1 (в подзадаче деформирования), k = K — моменту времени tn; сумма временных шагов подзадачи фильтрации равна временному шагу подзадачи деформирования.

На начало шага по времени t_n-1 для решения подзадачи с использованием координат узлов (вычисленных при решении подзадачи деформирования двухфазной среды) определяются: координаты центров элементов; векторы, соединяющие центры соседних КО (e^^a))) и расстояние между этими центрами — | eCa(a) |; площади граней — S^n 1) (с привлечением операции векторного произведения) и вектор нормали к общей грани этих КО — п(П^-1). _D                                                               (n-1) I (n-1) I c(n-1) (n-1)

ВычисленныечерезкоординатыузловпараметрыecA(a), |^eCA(_a)|, S(_a) , n(a) полагаютсяпостояннымина всех K временных шагах решения подзадачи фильтрации.

Таким образом, на n-м временном шаге решения подзадачи деформирования (далее в приведенных соотношениях индекс n, обозначающий шаг по времени, опущен) на основе разрешающего соотношения МКО (5) с использованием начальных и граничных условий для газовой фазы определяется плотность газовой фазы на k-м шаге решения подзадачи фильтрации в КО (давление газовой фазы в КО получено из решения МКЭ задачи деформирования пористой среды):

.⁽k) X ^pf -z^ (a)

(—А-Dk(H) у f    Pf

(k-1)    (k-1)

pfA(a) — pfC

I e^(n-1) eCA(a)

(n-1) eCA(a) e^(n-1) eCA(a)

(n-1) (n-1)

·n S

(a)       (a)

Atk

• (^n-1)V (n-1)

+^p

(k-1) f

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

По найденной плотности газовой фазы находится масса газа в каждом КО:

(k)      (n-1)   (k)

mf = Vf    -Pf , и скорость изменения плотности газовой фазы

•(k-1) = fMkJl

^Pf       Atk

.

С использованием данных о скорости изменения плотности газовой фазы вычисляется скорость изменения давления в КО за счет фильтрации газа:

•(k-1) p_f

(k-1) ρ_f

Mf

Rθf

где pf — скорость изменения давления газовой фазы за счет фильтрации воздуха, R, Q f, Mf — универсальная газовая постоянная, температура и молярная масса газа соответственно.

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

Путем интегрирования скорости изменения давления газовой фазы за счет фильтрации в каждом КО пересчитывается значение давления газовой фазы на конец текущего шага решения подзадачи фильтрации:

t_k

(k)     (k-1)         •(k-1)

pf = ^pf ' / ^pf   d^

tk-1

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

На первом временном шаге решения подзадачи фильтрации (при k = 1) давление газовой фазы соответствует давлению конца предыдущего/начала текущего шага решения подзадачи деформирования (pf 1) = pf⁰= p⁽n 1)).

В результате выполнения Этапа II алгоритма (по окончании пошаговой процедуры для фильтрации (t_k= t_n)) становятся известными параметры газовой фазы на текущем — n-м, временном шаге (на момент tn).

На Этапе III устанавливаются значения параметров на момент времени tn (на конец n-го шага) для их использования на следующем шаге. По вычисленному вектору скорости перемещений рассчитываются за временной шаг приращения перемещений u⁽ⁿ⁾= v^(n—1)At и координат узлов x⁽ⁿ⁾= x^(n—1) + v^(n-1)4t. По переопреленным координатам находятся функции формы [N]⁽ⁿ⁾и их производные [B]⁽ⁿ⁾.

По известным координатам узлов пересчитываются объемы (V⁽ⁿ⁾, V^n)) и объемные доли фаз (fl⁽ⁿ⁾, Y⁽ⁿ⁾). Отношение плотностей двухфазной среды в текущей и отсчетной недеформированной конфигурации

О р / Р = V /V в каждом элементе отыскиваются как

P⁽ⁿ⁾/р⁽⁰⁾= V⁽⁰⁾/v⁽ⁿ⁾.

(n)

По пересчитанным объемам V_fпереопределяется плотность газовой фазы в каждом КО:

Pf = m?"/Vf

.

Далее скорость изменения давления газовой фазы за счет изменения объема (см. формулу (9)) интегрируется, и назначения давления в КЭ:

t_n

(n)    ~(n) . I -!-(n —1) 7^

Pf == p? +   p?  ’d..

tn-1

Путем интегрирования скорости изменения тензора деформации Альманси и суммирования его компонент с их приращениями (в базисе ЛСК), вычисленными на конец предыдущего (начало текущего) шага по времени, tn определяется тензор деформации Альманси на момент времени tn (A(n) = A(n—1) + ^ A(   )d£). Для тензора tn-1

деформации Альманси находятся его собственные значения Ai и соответствующие им главные векторы pi. По Ai на основе соотношения (Ai = (1 — 2Ai)^—1/2) вычисляются собственные значения правого тензора искажений Ai. С использованием главных векторов pi рассчитываются косинусы углов между ними и базисными векторами ЛСК ki . По вычисленным значениям λi и компонентам матрицы косинусов устанавливаются компоненты тензора четвертого ранга Q⁽ⁿ⁾, а по найденным D⁽ⁿ⁾, W⁽ⁿ⁾, Ai — компоненты тензора логарифмического спина Q^ [13].

Далее выполняется определение переменной Z⁽ⁿ⁾, учитывающей взаимодействие газа и легочной ткани. С помощью интегрирования вычисляется взвешенный тензор напряжений Кирхгоффа двухфазной среды (S⁽ⁿ⁾= t_n

s^(n—1) + у ±^(n—) d^) и его девиатор (S⁽ⁿ⁾). После этого осуществляется переход на следующий временной шаг. tn-1

Для реализации описанного выше алгоритма, основанного на методах КЭ и КО, авторами разработан программный комплекс (на языке программирования C++ с привлечением библиотеки Eigen). Будучи нелинейной и обладающей значительной вычислительной сложностью, задача решалась с помощью пошаговых процедур. Вычисления производились на многопроцессорных системах с применением технологии параллельных вычислений. Рабочий процессор имел 12 ядер с частотой 3.3 ГГц. Распараллеливание вычислений выполнялось с использованием директив стандарта OpenMP. По тестовым оценкам использование 12 ядер повысило производительность счета почти в 10 раз.

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

Для визуализации результатов вычислений на основе разработанного программного комплекса использовался открытый графический кроссплатформенный пакет ParaView [50].

• 4.    Результаты

В результате численного решения согласно предложенному алгоритму связанной задачи течения воздуха в упруго-деформируемой пористой среде легких получено пространственное распределение параметров легочной ткани и воздушной фазы в разные моменты одного дыхательного цикла (вдоха/выдоха), в том числе поля скоростей

Таблица 1. Значения параметров

Параметр (размерность) Значение Молярная масса воздуха (кг^моль-1) 28.9640-3 Температура газовой фазы (К) 309.75 Универсальная газовая постоянная (м2 кг2 •с-2 К-1 •моль-1) 8.3144 Динамическая вязкость воздуха (Па·с) 18.6^6 Плотность воздуха в легких в начальный момент времени, а также плотность воздуха, поступающего в легкие при температуре 30 °С (кг^м- ) 1.165 Проницаемость (м2) 1.5840-9 [11] Модуль Юнга легочной ткани (Па) 10 000 [38] Коэффициент Пуассона легочной ткани 0.4 [49] перемещений, перемещений, деформаций, напряжений, потоков воздуха. В связи с ограничением объема в данной статье приводятся лишь картины полей давления газовой фазы в легких и массовые потоки воздуха через границы ΓAW/L — выходы из системы бронхов (входы в легкие), через которые воздух поступает/выходит во время вдоха/выдоха.

На рисунке 5 представлены поля давления газовой фазы в разные моменты деформирования (с шагом 0.5 с). Длительность дыхательного цикла составляет 4 с. За время t от 0 до 2 с происходит вдох (Рис. 5а–д), в течение двух следующих секунд — выдох (Рис. 5д–и). В конце выдоха легкие и давление воздуха в них возвращаются в исходное состояние: рисунок 5 а соответствует моменту начала, рисунок 5 и — концу дыхательного цикла (t = 0 си t = 4 с).

В начальный момент времени t0 (Рис. 5а), в момент перехода от вдоха к выдоху (Рис. 5д) и по завершению выдоха (Рис. 5и) давление в легких равно атмосферному. Элементы сетки, расположенные в непосредственной близости к областям выхода бронхов/входов в легкие (отмечены кружками, см. Рис. 5в), имеют давление более близкое к атмосферному в течение всего дыхательного цикла (Рис. 5а–и).

Во время вдоха легкие расширяются, их объем увеличивается, а давление воздуха в них, соответственно, становится ниже. Области легких, прилегающие к диафрагме, подвергаются наибольшему растяжению (преимущественно по вертикальной координате), следовательно, здесь давление имеет наименьшее значение. Воздух из области высокого давления распространяется в область низкого давления и постепенно распространяется по всему объему легких; давление в легких выравнивается. В конце вдоха скорость деформирования замедляется (t = 2 с соответствует максимальному растяжению легких), и воздух распространяется по всему объему, давление воздуха в легких выравнивается с атмосферным (Рис. 5д).

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

В процессе решения задачи устанавливаются массовые потоки воздуха на выходах из системы бронхов/входах в легкие (на границах Γ_AW/_L), которые в дальнейшем используются в качестве граничных условий в подмодели течения воздуха в воздухоносных путях человека (Рис. 6). Всего в рассматриваемой геометрии присутствует 31 вход (14 входов в левом легком и 17 в правом). На рисунке 6а приведены графики изменения массовых потоков воздуха через каждый из входов в легкие в режиме спокойного дыхания. Массовые потоки представляются синусоидами с различной амплитудой. Отличия между потоками воздуха через сечения входов (Рис. 6а) обусловлены разными величинами градиентов давления и площади входов (последняя варьируется от 1.31·10^-6до 1.16·10^-5м²).

Из рисунка 6б, на котором показано изменение суммарного массового потока (в левое и правое легкое и суммарный через 31 вход в легкие) в процессе спокойного дыхания, следует, что массовые потоки воздуха в левом и правом легком примерно одинаковы. Средний суммарный расход воздуха в процессе вдоха (за 2 с) при спокойном дыхании для двух легких составляет 4.43·10^-4кг/с, что при плотности воздуха 1.165 кг·м^-3 соответствует 0.38 л/с. Следовательно, за 2 с вдоха в легкие человека поступает 0.76 л (0.00076 м³), что соответствует известным физиологическим данным.

Был рассмотрен еще один сценарий — учащенное дыхание с периодом t* = 2 с при неизменных остальных параметрах, используемых в граничных условиях задачи (1)–(3). Суммарный массовый поток для учащенного

(a)      f^V R                [.        (б)     t' R                I.         (в)

МММ

1*3 R                I.        (г)            »                1         (д)     H ^R                         (е)

МММ

(ж)

(з)

(и)

Рис. 5. Поле давления газовой фазы в легких человека (вид спереди; буквой R обозначено правое легкое, буквой L – левое) в разные моменты времени одного цикла дыхания (вдоха-выдоха) t, [c] : 0 (а); 0.5 (б); 1 (в); 1.5 (г); 2 (д); 2.5 (е); 3 (ж); 3.5 (з); 4 (и)

дыхания представлен пунктирной линией на рисунке 6б. При учащенном дыхании средний суммарный расход за время вдоха (за 1 с) составляет 8.82·10^-4кг/с, а в легкие поступает или 0.757 л (0.000757 м³) воздуха. Таким образом, предложенная математическая модель может быть использована для описания течения воздуха в легких при различных режимах дыхания.

• 5.    Заключение

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

Рис. 6. Массовый поток воздуха: через выходы из бронхов/входы в легкие при спокойном дыхании (а), суммарный массовый поток через входы легких: обоих (сплошная линия), левого (точечная линия), правого (штрихпунктирная линия), а также суммарный массовый поток через входы в обоих легких при учащенном дыхании (пунктирная линия) (б)

(a)

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

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

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

Исследования выполнены при финансовой поддержке Министерства науки и высшего образования Российской Федерации (проект № FSNM-2023-0003).