Математическое моделирование колебаний биологических тканей, насыщенных жидкостью

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Автореферат

диссертации на соискание ученой степени доктора физико-математических наук

Математическое моделирование колебаний биологических тканей, насыщенных жидкостью

Специальность: 05.13.18 - Математическое моделирование, численные методы и комплексы программ

Маслов Леонид Борисович

Санкт-Петербург - 2010

1. Общая характеристика работы

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

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

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

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

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

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

Цели исследования:

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

Разработать новые математические модели биологических тканей и органов опорно-двигательного аппарата человека на основе полученных уравнений пороупругости; провести оценку эффективных свойств тканей.

Выявить взаимосвязи между динамическими свойствами основных биомеханических структур опорно-двигательного аппарата нижней конечности человека и характеристиками движения внутритканевой жидкости; представить теоретическое обоснование методов вибрационной диагностики и поддержания состояния длинных трубчатых костей при остеопорозе.

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

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

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

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

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

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

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

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

7. Методами имитационного моделирования и вибрационной диагностики исследовать динамические характеристики вынужденных колебаний биомеханических структур опорно-двигательного аппарата голени человека.

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

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

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

Положения, выносимые на защиту:

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

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

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

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

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

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

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

Новизна исследования:

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

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

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

Математические модели твердых и мягких биологических тканей, образующих костные и мышечные структуры опорно-двигательного аппарата человека, как пороупругих сред, насыщенных жидкостью, имеют отличия от известных моделей в части более точного описания динамических свойств за счет учета инерционного взаимодействия фаз и вязкости, учета анизотропии эффективной среды в «u-p» переменных; точнее описывают резонансные режимы колебаний структур и микроперемещений внутритканевой жидкости.

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

Одномерные дифференциальные уравнения, описывающие продольные и поперечные колебания пороупругого стержня, выведены из построенных динамических уравнений связанной линейной задачи пороупругости в «u-p» переменных и более полно, чем известные аналоги, учитывают силовое, физическое и инерционное взаимодействие фаз. Выявленные в результате параметрического анализа полученных численно-аналитических решений характерные динамические свойства пороупругих систем и закономерности движения жидкости в порах в различных частотных диапазонах отсутствуют в литературе.

Разработанные вычислительные алгоритмы решения поставленных краевых задач пороупругости в «u-p» и «u-w» переменных реализуют трехмерный метод конечных элементов на основе общего подхода взвешенных невязок и отличаются тем, что рассматривают модель анизотропной среды, учитывают вязкость и перечисленные виды взаимодействия фаз.

Разработан проблемно-ориентированный комплекс MechanicsFE для расчета динамического напряженно-деформированного состояния пороупругой среды под действием внешних гармонических сил, позволяющий исследовать специфические задачи, не имеющие реализации в универсальных коммерческих программах. По сравнению с моделями пористых сред, включенными в другие программы, разработанный алгоритм учитывает инерционное, силовое и физическое взаимодействие фаз и собственную вязкость материала. Численная реализация основана на «u-p» формулировке задачи поровязкоупругости, что позволяет значительно уменьшить число неизвестных при расчете вынужденных колебаний. Результирующая система уравнений, записанная в комплексных переменных, решается с учетом разреженности и симметрии глобальной матрицы.

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

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

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

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

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

Исследование находится в русле приоритетного направления развития науки, технологии и техники РФ - живые системы, и вносит вклад в развитие критических биомедицинских технологий жизнеобеспечения человека, биоинженерии, технологий создания биосовместимых материалов.

Апробация результатов. Основные результаты работы были представлены на следующих всероссийских и международных научных форумах: Всерос. конференции, приуроченной к 20-лет. ИПРИМ РАН «Механика и наномеханика структурно-сложных и гетерогенных сред» (Москва, 2009); IV-X Всерос. конференции по биомеханике (Нижний Новгород, 1998-2008, Саратов, 2010); 3-й Всерос. НТК «Фундаментальные исследования в технических университетах» (С.-Петеpбуpг, 1999); XXX-XXXI Неделе науки СПбГТУ (С.-Петеpбуpг, 2001, 2002); III-IV МНТК «Физика и радиоэлектроника в медицине и биотехнологии» (Владимир, 1998, 2000); Межд. школе-семинаре «Математическое моделирование и биомеханика в современном университете» (Ростов-на-Дону, 2005, 2009); VI Съезде травматологов и ортопедов России (Нижний Новгород, 1997); Конгрессе травматологов-ортопедов России «Новые имплантаты и технологии в травматологии и ортопедии» (Ярославль, 1999); IX-XV МНТК «Бенардосовские чтения» (Иваново, 1999-2009); II Int. Conference on Tools for Mathematical Modelling (St.-Petersburg, 1999); III-VI Int. Workshop on Nondestructive Testing and Computer Simulations in Science and Engineering (St.-Petersburg, 1999-2002); 15th World Conference on Non-Destructive Testing (Roma, 2000); 13th Conference of European Society of Biomechanics (Wroslav, 2002); 12th Int. Conference on Mechanics in Medicine and Biology (Lemnos, 2002); XXIXиme Congrиs de la Sociйtй de Biomйcanique (Paris, 2004) и др.

Разработанные модели биологических тканей и результаты исследования динамических свойств биомеханических систем докладывались на научных семинарах кафедр общей физики, механики и процессов управления, биомеханики и валеологии (СПбГПУ), теоретической и прикладной механики (ИГЭУ), травматологии, ортопедии и ВПХ (ИГМА), Института проблем машиноведения РАН. Во время зарубежных стажировок на рабочих совещаниях в институте передовых технологий и центре мехатроники компании Самсунг Электроникс (Южная Корея), в лаборатории прикладной механики Национального центра научных исследований (LMARC CNRS, Франция) обсуждались особенности конечно-элементного моделирования динамики человека и антропоморфных механизмов, новые модели эндопротезов тазобедренного сустава.

Внедрение результатов исследования. Выявленные закономерности динамического поведения биомеханических структур использованы при разработке средств вибрационной диагностики ахиллова сухожилия и трехглавой мышцы голени, стабильности остеосинтеза большеберцовой кости, выполненных в соответствии с планами госбюджетных НИР ИГЭУ по единому заказ-наряду Минобрнауки РФ (1997-2010) на фундаментальные научные исследования в НИЛ «Биомеханика» (ИГЭУ-ИГМА-ИвГУ). Практические выводы и рекомендации из проведенных теоретических исследований отражены в отчетах: Разработка теоретических основ методов вибрационной диагностики, лечения и реабилитации состояния биомеханических объектов / Отчет по НИР, рук. Шапин В.И., РИ № 01.2001.03153; Диагностические методы идентификации, вибродиагностики и частотно-резонансной терапии биомеханических объектов / Отчет по НИР, рук. Шапин В.И., РИ № 02.2006.09975, и в рамках текущей НИР «Разработка вибрационных резонансных методов диагностики, идентификации и реабилитации биомеханических объектов», РИ № 1.2.06.

В конкурсе «Молодые ученые - малому предпринимательству», проводимому Фондом поддержки малого предпринимательства, в 2001 г. получена поддержка проекта «Компьютерное моделирование конструкции эндопротеза тазобедренного сустава». Методы и принципы компьютерного моделирования гетерогенных материалов и биомеханических структур использованы при разработке ножек эндопротезов тазобедренного сустава из композитных полимерных и пористых титановых материалов в рамках совместного проекта с группой исследователей лаборатории прикладной механики (LMARC CNRS, Франция). Закономерности движения человека, данные о динамических свойствах тканей и биомеханических структур голени человека, принципы конечно-элементного анализа динамики упругих сред применены при разработке опорно-двигательного механизма шагающего робота в центре мехатроники компании Самсунг Электроникс (US Patent 7303031; 2005/0173164).

Работа соискателя «Теория колебаний биомеханических структур как пористых насыщенных жидкостью сред» представлена в Сводный реестр Федерального агентства по образованию как наиболее значимый результат научных исследований, полученных работниками вуза (ИГЭУ) в 2009 году. В 2009 году на основе результатов многолетних исследований и достижений соискателя в ИГЭУ создан НОЦ моделирования и диагностики структурно-сложных сред и систем. Материалы диссертации и разработанное программное обеспечение внедрены в авторский учебный курс «Численные методы механики» (ИГЭУ), использованы при чтении лекций «Computational Biomechanics» в Высшей национальной школе механики и микротехники (ENSMM, Франция).

Разработанное программное обеспечение зарегистрировано в Российском агентстве по патентам и товарным знакам (Роспатент) и получено авторское свидетельство № 2009616004 «Программная система конечно-элементного анализа пороупругих механических конструкций MechanicsFE».

Публикации по теме исследования. По материалам диссертации опубликовано 63 печатные работы, в том числе монография, 20 статей в ведущих рецензируемых журналах согласно Перечню ВАК Минобрнауки РФ, 4 авторских свидетельства, 16 статей в сборниках и материалах конференций.

2. Содержание работы

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

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

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

, , (1)

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

В предположении упругой модели твердой фазы и модели идеальной сжимаемой жидкости, фазовые уравнения преобразуются в определяющие соотношения пороупругой среды в «u-p» переменных (перемещение упругого каркаса - давление поровой жидкости). Кинематической переменной жидкой фазы является вектор относительного перемещения или дивергенция этого вектора , имеющая физический смысл относительного изменения объемного содержания жидкости в порах. Для случая анизотропии упругих и гидростатических свойств уравнения примут вид:

, , (2)

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

В результате перехода к новым переменным возникают эффективные материальные константы, представленные тензором упругих модулей в дренированном состоянии и тензором коэффициентов Био:

, . (3)

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

, , (4)

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

Сила для схемы Рахматулина силового взаимодействия и совместного деформирования фаз складывается из равновесной и диссипативной составляющих:

, (5)

где В - тензор коэффициентов вязкого трения, выражаемый через тензор проницаемости и вязкость жидкости: , .

С учетом выражения тензора напряжений в идеальной сжимаемой жидкости и соотношения (5) уравнения (4) примут вид:

, , (6)

где - полная объемная сила; - полная плотность среды.

Динамические уравнения (6) содержат дополнительную переменную w, которая выражается с помощью комплексного преобразования Лапласа:

, (7)

где - приведенная комплексная гидравлическая проницаемость пористой среды.

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

, (8)

,

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

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

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

Рассмотрим модифицированные определяющие соотношения двухфазного материала в следующем виде:

, (9)

.

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

Формулы (9) учитывают собственные диссипативные свойства, обусловленные вязкостью материала фаз, в то время как взаимное влияние фаз друг на друга, осуществляемое на границе раздела материалов за счет сил вязкого трения, учитывается в уравнениях баланса количества движения. Показано, что при использовании модели вязкоупругого материала (9), могут быть составлены определяющие соотношения эффективной среды с шестью степенями свободы относительно векторов u и w с новыми материальными константами, описывающими среду в недренированном состоянии:

, (10)

,

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

Давление жидкости в порах входит во второе определяющее соотношение системы(10) и выражается через перемещения скелетона и флюида:

. (11)

В результате преобразований уравнений баланса количества движения сформулирована гиперболическая краевая задача пороупругости в «u-w» переменных, включающая в себя диссипативные слагаемые, пропорциональные вязкости жидкости и полной вязкости двухфазного материала:

(12)

Дифференциальные уравнения (12) являются основными расчетными соотношениями для решения динамической задачи пороупругости в «u-w» переменных. Согласно общей постановке задачи механики сплошной среды краевые условия представляют собой кинематические ограничения на степени свободы и заданные распределенные силы на частях поверхности:

, , (13)

, .

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

, ; (14)

,;

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

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

, (15)

,

где , - тензоры упругих модулей эффективной среды в дренированном состоянии и коэффициентов эффективных напряжений Био аналогично определениям (3) для случая простой пористости:

, . (16)

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

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

, (17)

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

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

, (18)

,

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

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

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

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

, (19)

в котором выражение тензорной функции имеет вид:

,

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

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

Система обыкновенных дифференциальных уравнений (19) относительно компонент дополнена начальным условием , согласно которому при нулевой пористости эффективные характеристики среды равны упругим модулям матрицы материала.

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

, (20)

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

(21)

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

Для расчета тензора коэффициентов эффективных напряжений и гидростатической константы используется принцип суперпозиции, согласно которому напряженно-деформированное состояние упругого тела, возникающее в результате действия комплекса нагрузок, может быть представлено как суперпозиция полей напряжений и деформаций от отдельных нагрузок. В случае представительного элемента пористого материала, насыщенного жидкостью, полный тензор эффективных упругих напряжений (2) складывается из чисто упругой части, соответствующей деформированию твердого дренированного каркаса , и гидростатической части, соответствующей давлению поровой жидкости . Последовательно рассматривая дренированное состояние среды, характеризуемое тензором деформаций твердой фазы и отсутствием давления жидкости в порах , и состояние II, которое возникает при отсутствии деформаций пористого каркаса , заполненного жидкостью с давлением , в конечном итоге получаем формулы:

, , (22)

где - модуль объемного сжатия жидкости.

Формулы (18)-(22) реализованы в программном модуле, предназначенном для расчета материальных характеристик тканей, требуемых для динамического пороупругого анализа. Хотя данные методы в целом известны в литературе, новизной и достоинством данной части работы является то, что впервые разработан универсальный математический алгоритм и проведен численный анализ эффективных характеристик твердых и мягких биологических тканей как пороупругих сред при произвольных значениях пористости в случае принятия модели трансверсально-изотропного тела.

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

Рис. 1. Графики эффективных модулей губчатого костного вещества в дренированном состоянии среды: а - модули Юнга и сдвига; б - коэффициенты Пуассона; в - коэффициенты эффективных напряжений; г - гидростатическая константа.

Таблица 1. Эффективные модули тканей, используемые в модели голени

Параметры среды	Тип костного вещества
	Компактное	Губчатое	Костномозговое
	0.05	0.1	0.2	0.65	0.8	0.995
, кг/м3	1494	1468	1416	1182	1021	1021
, ГПа	14.3	12.2	10.7	1.55	0.415	1.3610-4
, ГПа	22.3	21.1	13.1	3.73	0.513	1.3610-4
, ГПа	4.69	4.1	4.1	0.618	0.171	0.56610-4
, ГПа	6.98	6.26	4.52	0.97	0.190	0.56610-4
	0.424	0.408	0.301	0.253	0.217	0.20
	0.270	0.224	0.266	0.139	0.191	0.20
	0.241	0.416	0.504	0.942	0.974	1.00
	0.211	0.367	0.435	0.893	0.970	1.00
, ГПа	0.089	0.185	0.390	1.420	1.750	2.29
, м2	1.4710-20	1.4710-20	7.910-16	1.610-12	3.210-9	2.110-7
, м2	6.3610-13	6.3610-13	1.110-11	3.010-10	7.910-10	6.310-6
, Пас	110-3	110-3	37.510-3	37.510-3	37.510-3	37.510-3

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

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

 (23)

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

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

. (24)

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

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

Особую роль в системе (23) играют комплексные частотно-зависимые коэффициенты и , имеющие следующий вид:

,

(25)

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

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

В частности, показано, что область существования апериодического решения ограничена двумя кривыми в плоскости (,). Из диаграммы (рис. 2) следует, что в области изменения параметров существует особая точка, характеризуемая тем, что при любых значениях < 8 и < 27 имеет место колебательное движение, и, наоборот, в заштрихованной области колебания отсутствуют. При вариациях , происходит значительное изменение амплитудно-частотных характеристик системы (рис. 3). В высокочастотной области становится существенным влияние , , связанных с инерционным взаимодействием фаз, и амплитуды перемещений и давления сильно возрастают.

Рис. 2. Области существования гармонических колебаний: и - колебания отсутствуют; - в системе возможны затухающие колебания

Пятая глава диссертации посвящена разработке конечно-элементного алгоритма расчета колебаний пороупругого анизотропного тела.

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

а); ; б);

Рис. 3. Амплитудно-частотные характеристики (АЧХ) стержня при значениях параметров =1, =27 (а) и =1, =1 (б) и =10^-2, =10^-3 (---); =10^-1, =10^-3 (- - -); =10^-1, =10^-2 (····)

Тем самым, построены системы конечно-элементных уравнений динамики пороупругого тела в «u-w» переменных:

(26)

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

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

, (27)

где - вектор амплитудных значений узловых перемещений.

В случае смешанной «u-p» формулировки задачи пороупругости (8) дискретная система содержит изображения искомых функций и параметр преобразования Лапласа. Если нагрузки являются произвольными функциями от времени, то для нахождения динамического отклика системы на силовое воздействие требуется использовать обратное преобразование Лапласа и специальные численные алгоритмы. Однако если внешние нагрузки изменяются гармоническому закону, то отклик линейной системы также представляет собой гармонические во времени функции, а параметр преобразования Лапласа равен , что дает следующую систему матричных уравнений:

, (28)

,

где ,M- стандартные глобальные матрицы жесткости и массы; - дополнительная глобальная матрица массы; , - глобальные матрицы взаимного влияния распределения давления жидкости в порах и деформации упругого скелетона; D, - глобальные матрицы насыщения и проницаемости; U, P - глобальные векторы комплексных амплитудных значений перемещений и давлений в узлах конечно-элементной сетки.

Для разработки специализированного комплекса, позволяющего рассчитывать вынужденные пространственные колебания пористых упругих тел, насыщенных жидкостью или газом, были выбраны уравнения теории пороупругости в «u-p» переменных. Данный выбор обоснован меньшим числом скалярных переменных; прямым соответствием динамических уравнений относительно изображений задаче о гармонических колебаниях пороупругого тела; перспективной возможностью расширения данной формулировки на пористые материалы с двойной системой пор.

Однако, поскольку программа предназначена, в первую очередь, для анализа биологических тканей, проявляющих существенно диссипативные свойства, то для корректного численного решения в «u-p» формулировку включено слагаемое, учитывающее эффективную вязкость двухфазного материала. В результате первое уравнение системы (28) примет вид:

, (29)

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

.

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

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

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

а) б)

Рис. 4. Конечно-элементная модель голени человека: а - внешний вид; б - вид внутренних структурных частей голени: большеберцовой и малоберцовой кости, ахиллова сухожилия, трехглавой мышцы

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

...