Full Text

Введение В настоящее время доказано негативное влияние неблагополучной санитарно-эпидемиологической обстановки на состояние здоровья человека [1-3], в том числе на органы дыхания [4, 5]. Существующие методы медицинской диагностики [6, 7] позволяют оценивать состояние здоровья пациента на момент измерения, но не дают возможности прогнозировать изменение этого состояния на сколь-нибудь длительный период и количественно оценивать влияние комплекса факторов среды обитания на это изменение. Для прогнозирования развития нарушений в человеческом организме, обусловленных действием факторов среды обитания, научным коллективом, к которому относятся авторы данной статьи, разрабатывается многоуровневая математическая модель для описания эволюции функциональных нарушений человеческого организма [8]. Данная модель включает подмодели «мезоуровня» пищеварительной [9], иммунной и эндокринной систем [10]. Необходимость учитывать доставку химических веществ из атмосферного воздуха в организм человека ингаляционным способом, а также прогнозировать развитие заболеваний болезней органов дыхания обусловливает необходимость разработки углубленной подмодели «мезоуровня» дыхательной системы человека [11], которую в силу ее сложности можно рассматривать как отдельную модель. На текущий момент в рамках разработки математической модели дыхательной системы авторами решена задача о нестационарном течении воздуха (многокомпонентной смеси газов с твердыми частицами) и оседании твердых частиц различных фракций в крупных нижних воздухоносных путях (трахея и первые 3 генерации бронхов) [11], рассмотрено течение воздуха при наличии обструктивных нарушений бронхов. Воздухоносные пути органов дыхания человека по существу являются биологическими каналами сложной формы; в этой связи для описания механических аспектов течения воздуха в дыхательных путях широко используются модели газовой динамики в трехмерной постановке (например, для верхних дыхательных путей - в работах [12, 13], в нижних - в [14, 15]). В современных работах наблюдается тенденция перехода к персонализированной медицине, в том числе рассматривается реальная геометрия дыхательных путей, полученных с использованием компьютерной томографии [16, 17]. С точки зрения моделирования процесса дыхания в запыленных условиях и для моделирования переноса лекарственных препаратов ингаляционным способом представляют интерес работы [18, 19]. С ростом вычислительных мощностей и развитием вычислительных методов увеличивается степень детализации респираторного тракта; так, в [14] моделируется 17 генераций дыхательных путей, рассматривается 1400 элементов каналов. При моделировании движения воздуха по воздухоносным путям с использованием средств вычислительной аэромеханики возникают существенные сложности при описании геометрии всей иерархии дыхательных путей вплоть до альвеол (в легких взрослого человека содержится около 600-700 млн альвеол) и соединяющих их каналов. При этом детальное математическое описание процесса дыхания требует чрезвычайно больших вычислительных ресурсов. С целью преодоления указанных трудностей представляется возможным для описания течения газовой среды в мелких дыхательных путях и альвеолах легких использовать модель пористой среды с привлечением соотношений теории фильтрации [20]. Подобный подход был предложен также в работах [21-23], однако в цитируемых работах рассматривалось лишь течение газа через недеформируемую пористую среду с неизменными свойствами под действием заданного градиента давления, совместное решение задач деформирования легочной ткани и фильтрации не рассматривалось. Иначе говоря, пористая среда рассматривалась как «замороженная», недеформируемая, проводящая для газовой фазы среда. В то же время именно задача фильтрации воздуха под действием градиентов давления, создаваемых деформируемой пористой средой человека, представляется приемлемой моделью для описания движения газа в легких; при этом она оказывается существенно более сложной, чем модель фильтрации газа через «пассивную» среду. В частности, в предлагаемой модели необходимо учитывать геометрическую нелинейность среды, обусловленную большими градиентами перемещений, взаимодействие газа и легочной ткани. В рамках поставленной цели - разработка подмодели «мезоуровня» дыхательной системы человека - в настоящей работе основное внимание уделено формулировке необходимой для достижения поставленной цели конститутивной модели для описания движения воздуха в пористой среде легких. Концептуальная постановка задачи Дыхательная система человека представляет собой сложную биомеханическую систему, обеспечивающую доставку воздуха из атмосферы в легкие и газообмен между воздухом в легких и кровью. Структурными элементами дыхательной системы являются верхние воздухоносные пути (носовая полость, носоглотка, ротоглотка) и нижние дыхательные пути (трахея, система бронхов, легкие). Нижние дыхательные пути представляют собой древовидную структуру: трахея разделяется на два главных бронха, входящие в левое и правое легкие; в легких человека бронхи продолжают ветвление на более мелкие дыхательные пути; ветвление заканчивается альвеолами, имеющими сферическую форму (рис. 1, 2). Стенки альвеол, окутанные плотной сетью кровеносных капилляров, являются мембраной, через которую осуществляется обмен газами (О2 и СО2) между воздухом в легких и кровью посредством диффузии за счет градиента концентраций. Диаметр одной альвеолы составляет 0,25-0,3 мм [24, 25], в одном легком содержится порядка 300 млн альвеол, общая площадь их газообменной поверхности составляет около 80 м2 [26], что позволяет обеспечивать эффективный газообмен. Движение воздуха из атмосферы в легкие и обратно происходит за счет разницы между атмосферным давлением и давлением в легких. Изменение давления в легких обусловлено изменением их объема. У среднестатистического взрослого индивида изменение объема легких в процессе спокойного дыхания составляет порядка 15 %, в процессе глубокого - порядка 55 %. К изменению объема легких приводит изменение объема грудной клетки, которое осуществляется благодаря работе костно-мышечной системы. Моделирование последней представляет собой весьма сложную самостоятельную проблему, выходящую за рамки предлагаемой работы. Полная математическая модель дыхательной системы человека в силу своей сложности может рассматриваться совокупностью подмоделей, взаимосвязанных между собой и описывающих движение воздуха в воздухоносных путях (верхних, крупных нижних), деформируемых легких человека, газообмен между воздухом в легких и кровью. В стенках верхних и крупных нижних воздухоносных путей содержится хрящевая ткань, которая не позволяет им деформироваться во время дыхания. Подмодель движения воздуха в воздухоносных путях была разработана на предыдущих этапах исследования и представлена в статьях [11, 20], течение воздуха в ней описывается уравнениями газовой динамики. В данной статье особое внимание уделяется подмодели движения воздуха в деформируемых легких человека. Легкие человека, образованные мелкими дыхательными путями и альвеолами с содержащимся в них воздухом, будем рассматривать сплошной деформируемой насыщенной пористой средой, заключенной во внутреннюю камеру с изменяющимся объемом (подвижными стенками). Насыщенная газом пористая среда представлена двухфазной сплошной средой, одна из фаз которой - деформируемый скелет (матрикс, каркас) среды, описываемый моделью деформируемого твердого тела, вторая фаза - газ, заполняющий поровое пространство. Выделим представительный объем двухфазной среды, содержащий не менее 103 альвеол; введение данного понятия необходимо для описания легких как пористого континуума и движения газа с использованием теории фильтрации. Для описания относительного движения газовой фазы в пористой среде легких используется закон фильтрации Дарси [28, 29]. Легкие человека неоднородны, поэтому в общем случае состав и свойства различных участков легких отличаются друг от друга. Проницаемость легочной ткани в общем случае анизотропна и описывается тензором второго ранга. Рис. 1. Дихотомическое деление нижних воздухоносных путей согласно морфометрической модели Вейбеля [25, 24] Fig. 1. Dichotomous division of the lower airways according to the Weibel morphometric model [25, 24] Рис. 2. Увеличенный фрагмент терминальных бронхиол, заканчивающихся альвеолами [27] Fig. 2. Enlarged fragment of terminal bronchioles ending in alveoli Скорость частиц двухфазной среды определяется через скорости движения частиц каждой из фаз как среднемассовая скорость смеси [30, 31]. В силу того что масса газа, содержащегося в легких, пренебрежимо мала по сравнению с массой частиц твердого каркаса, скорость двухфазной среды можно считать совпадающей со скоростью твердой фазы ( ). Абсолютную скорость движения газа (за счет переноса вместе с деформируемой средой) можно считать совпадающей со скоростью частиц твердой фазы (скелета) и со скоростью двухфазной среды, относительная скорость движения газа (относительно твердой фазы) в пористой среде легких, возникающая за счет градиента давления, определяется с использованием методов теории фильтрации. Воздух в общем случае представляет собой многокомпонентную смесь газов, на данном этапе работы будем рассматривать его как гомогенный газ. В полости носа воздух нагревается/охлаждается до температуры тела, при исследовании фильтрационного течения воздуха в пористой среде легких изменением температуры воздуха пренебрегается. Вследствие геометрической нелинейности и неопределенности конфигурации в каждый момент деформирования для исследования предпочтительно использовать скоростную постановку [32] краевой задачи теории упругости. В качестве определяющего соотношения (ОС) в этом случае может быть использован гипоупругий закон (или эквивалентный ему гиперупругий закон) [33, 34]; в качестве меры деформации используется логарифмическая мера деформации [35]: , . (1) где - взвешенный тензор напряжений Кирхгоффа двухфазной среды, состоящей из воздуха и твердой фазы, связанный с тензором напряжений Коши ( ) через отношение плотностей двухфазной среды в отсчетной ( ) и текущей ( ) конфигурациях ( [36]); П - 4-валентный тензор упругих характеристик (принимается гипотеза об изотропии свойств); - логарифмический спин; Н - тензор деформации Генки, определенный в актуальной конфигурации; верхний индекс - обозначение логарифмической коротационной производной; [37]; - тензор скорости деформации (в нелинейной механике его обычно называют «тензор деформации скорости» - «деформатор» над полем скорости перемещений [32]); все характеристики в ОС относятся к двухфазной среде. Следует отметить, что гипоупругий закон (1) для изотропного материала удовлетворяет требованиям независимости от выбора системы отсчета [38], отсутствия диссипации энергии и замкнутости траектории напряжений при циклическом упругом деформировании [33]. Газ в легких (альвеолах и каналах) и легочная ткань взаимодействуют друг с другом; изменение объема легких вызывает изменение давления газа, находящегося в легких, которое, в свою очередь, приводит к изменению напряжений в легочной ткани. Давление газа действует во все стороны одинаково, не меняя форму материала (девиаторную часть тензора деформаций), а лишь меняя его объем (шаровую часть). В то же время сопротивление газа изменению формы (сдвиговым деформациям) пренебрежимо мало по сравнению со сдвиговой жесткостью твердого тела, что позволяет не учитывать вклад газовой среды в девиатор тензора напряжений двухфазной среды. В силу вышесказанного тензор напряжений для двухфазной среды определяется суммой девиатора напряжений в твердой фазе, умноженного на долю твердой фазы, и шаровой составляющей для двухфазной среды легких: , (2) где - скорость изменения среднего напряжения Кирхгоффа двухфазной среды; - единичный тензор; - девиатор тензора напряжений Кирхгоффа в твердой фазе; - объемная доля твердой фазы. В силу малости изменения плотности твердой фазы девиатор тензора напряжений Кирхгоффа равен девиатору тензора напряжений Коши. Скорость изменения девиатора напряжений твердой фазы определяется согласно соотношению , (3) где - девиаторная составляющая тензора скорости деформации; - параметр Ламе (модуль сдвига твердой фазы). В разрабатываемой математической модели предполагается использование граничных условий кинематического типа. Области входа крупных воздухоносных путей в легкие полагаются неподвижными (относительно тела человека). В начальный момент времени легочная ткань находится в естественном ненапряженном состоянии (момент перехода от выдоха к вдоху), когда давление в легких равно давлению на входе и равно атмосферному. В силу существенной нелинейности задачи для ее решения необходимо разработать и использовать пошаговую процедуру с применением достаточно малых (с точки зрения достижения требуемой точности) шагов по времени. На каждом малом промежутке времени массу газа, находящегося в альвеоле, можно считать постоянной; при изотермических условиях произведение давления и объема газовой фазы постоянно. Математическая постановка задачи Постановка краевой задачи описания легких как пористой среды содержит следующие соотношения: - уравнение равновесия в скоростях: , ; (4) - соотношение для скорости изменения взвешенного тензора напряжений Кирхгоффа: , , (5) - соотношения для тензора скорости деформации и его девиатора: , ; (6) - соотношение для логарифмической коротационной производной логарифмической меры деформации: ; (7) - соотношение для тензора вихря: ; (8) - соотношение для компонент логарифмического спина в базисе главных векторов меры Н: (9) - соотношение для описания относительного движения газовой фазы, полученное из закона Дарси: , , ; (10) - соотношение для изменения давления газовой фазы в зависимости от изменения объема газовой фазы, полученное из закона Бойля - Мариотта: , , ; (11) - начальные условия: (12) , , , ; (13) - граничные условия: , , , (14) , , , (15) где - замкнутая область двухфазной насыщенной пористой среды легких человека ( ); ΩЛ - внутренность области легких; ГЛ - граница области легких; ГДП/Л - границы выходов из крупных дыхательных путей, являющихся одновременно входами в легкие человека; - оператор Гамильтона, определенный в актуальной конфигурации; σ - тензор напряжений Коши для двухфазной среды; v - скорость частиц двухфазной среды; - скорость газовой фазы; - тензор проницаемости (второго ранга) пористой среды, м2; - динамическая вязкость, Па∙с; - тензор логарифмических деформаций твердой фазы (мера деформации Генки), определенный в терминах актуальной конфигурации; - давление газовой фазы, Па; - объем газовой фазы; - тензор вихря; - относительные главные удлинения. Основной трудностью при формулировке конкретных определяющих соотношений (5) для рассматриваемой двухфазной среды является построение связи скорости изменения среднего напряжения и изменения объема двухфазного континуума. При этом газ (воздух) находится большей частью в альвеолах, соединенных друг с другом и с воздухоносными путями системой весьма тонких в сечении каналов, что затрудняет перетекание воздуха из внутренних областей в бронхи и обратно при изменениях давления газовой среды. Для описания этого тела с весьма сложной геометрией предложено использовать модель пористой среды с пористостью, определяемой отношением суммарного объема альвеол и каналов к величине рассматриваемого представительного объема. Действительно, очевидна обозначенная выше проблема, что твердый несплошной деформируемый каркас при отсутствии в порах какой-либо среды ведет себя совершенно иначе, чем пористая среда, содержащая в себе газ с затрудненным истечением его из рассматриваемого объема. В связи с указанной сложностью осуществлена постановка и решение вспомогательной задачи. Предполагая, что на малом шаге по времени (напомним, что решение задачи предполагается реализовать с использованием пошаговой процедуры) изменением количества газа в рассматриваемом объеме можно пренебречь, в качестве расчетной схемы для определения указанной связи гидростатических характеристик используется представительный объем в форме замкнутой сферы, заполненной воздухом. Масса оболочки равна массе твердотельного каркаса, объем внутренности сферы в отсчетной конфигурации равен суммарному объему альвеол и каналов в естественной конфигурации представительного объема. К сожалению, авторам не удалось получить в аналитическом виде решения геометрически нелинейной задачи о всестороннем сжатии тонкостенного полого шара, заполненного воздухом. Численное решение не позволяет сформулировать искомую связь ~I1(D) в виде простого алгебраического уравнения, необходимого для создания эффективной вычислительной процедуры. В связи с этим было получено аналитическое решение соответствующей геометрически линейной задачи. На основе этого решения записано соотношение для связи скорости изменения среднего напряжения и скорости изменения объема I1(D) (I1(D) - первый инвариант тензора D) двухфазной среды легких, учитывающее взаимодействие между легочной тканью и воздухом, содержащимся в легких [39]. Далее с помощью численного решения соответствующей геометрически нелинейной задачи проведена оценка возможности использования аналитического решения задачи в геометрически линейной постановке для умеренно больших изменений объема, характерных для процесса дыхания. Полученное соотношение для скорости изменения среднего напряжения в терминах тензора напряжений Кирхгоффа двухфазной среды имеет вид (16) где , , ; , - параметры Ламе твердой фазы; - объемная доля газовой фазы; - объем двухфазной среды в актуальной и отчетной конфигурациях. Численное решение вспомогательной задачи для определения скорости изменения среднего напряжения Рассматривается однородный изотропный упругий тонкостенный полый шар, заполненный воздухом, подвергаемый равномерному внешнему сжатию (или расширению). Объем газа, содержащегося в шаре, соответствует объему газообразной фазы в представительном объеме легких в рассматриваемый момент процесса; объем стенки - объему твердой фазы, общий объем - объему двухфазной среды. Сжатие шара происходит при постоянной температуре газа (процесс изотермический). Для описания нелинейного поведения однородного изотропного материала используется гиперупругий неогуковский закон. Рассматривается геометрически нелинейная одномерная (с учетом сферической симметрии) задача, в качестве меры деформации используется мера Генки. На внешней стенке шара заданы кинематические граничные условия (всестороннее сжатие/растяжение); на внутренней стенке шара - силовые (давление газа). Среднее напряжение (гидростатическое напряжение) для двухфазной среды определим как радиальную компоненту напряжения на внешней поверхности шара ( ). Весь процесс нагружения был разбит на K временных шагов. Для определения отсчетной естественной конфигурации на внутренней и внешней поверхностях оболочки предварительно задавалось давление, равное атмосферному. На каждом временном шаге использовались соответствующие граничные условия: изменение давления (как функция от изменения внутреннего объема) на внутренней границе и приращение перемещений на внешней. Численное решение геометрически нелинейной задачи было получено с использованием программного продукта ANSYS Mechanical APDL (ANSYS Parametric Design Language). В результате численного решения была получена радиальная компонента тензора напряжений Коши в узлах сетки на внешней поверхности шара в процессе пошагового нагружения (для различных моментов нагружения). Аналогичные граничные условия были использованы в полученном на предыдущих этапах работы аналитическом решении вспомогательной задачи (в соотношении для определения радиальной компоненты тензора напряжений Коши на внешней поверхности шара) [39]. Сравнение аналитического и численного решений (радиальных компонент напряжений) выполнялось на каждом временном шаге для нескольких сценариев всестороннего сжатия и расширения сферы при реальных параметрах легочной ткани (модуль Юнга взят равным 60 кПа, коэффициент Пуассона - 0,45) [40]. Относительная погрешность численного и аналитического решений на каждом временном шаге определялась по соотношению (17) где - k-й временной шаг; - величина радиального напряжения, полученная аналитически; - величина радиального напряжения, полученная численно; - относительная погрешность, %. Определены также средняя и максимальная относительные погрешности по соотношениям: , (18) , (19) где - средняя относительная погрешность; - максимальная относительная погрешность; K - количество временных шагов. Результаты Сравнение аналитического и численного решений было выполнено для двух сценариев. Первый сценарий подразумевал уменьшение внутреннего объема за счет всестороннего сжатия на 52,95 %; второй сценарий соответствовал увеличению объема за счет всестороннего расширения на 17,61 %. Результаты сопоставления решений приведены в таблице. Для сравнения в процессе спокойного вдоха объем легких человека изменяется примерно на 15 %, во время максимально глубокого вдоха - около 55 %. Результаты сравнения аналитического и численного решения Results of comparison of analytical and numerical solutions Сценарий Изменение внутреннего объема, % Средняя относительная погрешность по временным шагам, % Максимальная относительная погрешность, % 1 -52,95 1,46 3,22 2 17,61 0,09 0,24 В первом сценарии средняя погрешность между аналитическим и численным решением составила 1,46 %, во втором сценарии - 0,09 % (см. таблицу). На рис. 3 приведена зависимость радиального напряжения на внешней поверхности сферы от изменения объема для линейной и нелинейной постановок для первого сценария, а также изменение максимальной (по модулю) разности численного и аналитического решения. Рис. 3. График изменения значений радиальной компоненты тензора напряжения Коши на внешней поверхности сферы в зависимости изменения объема для линейной и нелинейной постановок для первого сценария и изменения погрешности численного и аналитического решений Fig. 3. Graph of the change in the values of the radial component of the Cauchy stress tensor on the outer surface of the sphere depending on the volume change for linear and nonlinear statements for the first scenario and changes in the error of numerical and analytical solutions Различие численного и аналитического решения в процессе изменения объема (в результате всестороннего сжатия сферы на 52,95 %) постепенно увеличивается и достигает максимального значения 3,22 %. Во втором сценарии максимальная погрешность составила 0,24 %. Основываясь на приведенных результатах, можно констатировать, что полученное аналитическое решение находится в удовлетворительном соответствии с решением аналогичной задачи в нелинейной постановке как при спокойном, так и при глубоком дыхании. Соотношение (16) дополняет постановку краевой задачи описания легких как пористой среды (соотношения (4)-(15)) и может быть использовано для описания связи скорости изменения среднего напряжения и изменения объема двухфазного континуума для умеренно больших деформаций. Заключение и перспективы Таким образом, представлена конститутивная модель насыщенной пористой среды с учетом взаимодействия воздуха и твердого деформируемого каркаса, необходимая для описания течения воздуха в легких человека. Представлена концептуальная и математическая постановки задачи описания процесса дыхания. С помощью численного решения вспомогательной задачи о связи среднего напряжения и изменения объема при всестороннем растяжении-сжатии тонкостенной сферы в геометрически нелинейной постановке подтверждена гипотеза о приемлемости аналитического решения, полученного в геометрически линейной постановке, для умеренно больших деформаций. Реализация предложенной модели в связи с ее существенной нелинейностью предполагает использование пошаговой процедуры. При разработке модели остаются сложные моменты, решение которых еще предстоит получить (в частности, установление компонент тензора проницаемости для различных участков пористой среды легких с учетом анизотропии и нелинейности (зависимости от изменения плотности пористой среды)). Данная задача является предметом дальнейших исследований. Алгоритм решения задачи с использованием пошаговой процедуры, а также результаты по течению газов в легких человека будут представлены в следующих публикациях авторов.

About the authors

P. V Trusov

Federal Scientific Center for Medical and Preventive Health Risk Management Technologies for the Population; Perm National Research Polytechnic University

N. V Zaitseva

Federal Scientific Center for Medical and Preventive Health Risk Management Technologies for the Population

M. Yu Tsinker

Federal Scientific Center for Medical and Preventive Health Risk Management Technologies for the Population

References

Statistics

Views

Abstract - 385

PDF (Russian) - 269