МОДЕЛИРОВАНИЕ БИОЛОГИЧЕСКИХ ПРОЦЕССОВ И СИСТЕМ Презентационные материалы к курсу лекций для студентов направления «Биотехнические системы и технологии» Профессор каф. «Приборостроение», ДГТУ Литвин А.В..

Лекция 1 « Методологическая основа моделирования » Объектом (лат. objectum – предмет) является все то, на что направлена человеческая деятельность. Объекты характеризуются процессами, которые протекают в них или в которых они участвуют.

Лекция 1 « Методологическая основа моделирования » БИООБЪЕ́КТ 1. Живой организм как объект какого-либо исследования. 2. Живой организм как объект наблюдения или изучения

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

Лекция 1 « Методологическая основа моделирования » Для получения информации в биологических системах применяют различные методы: оптические, электрометрические, методы хемилюминесценции, рентгеноструктурный анализ, ЯМР- и ЭПР-спектроскопию, лазерную спектроскопию, метод меченых атомов и др.

Лекция 1 « Методологическая основа моделирования » Особенности человека как объекта исследования: невозможность получения полного описания физиологических процессов, происходящих в живом организме; использование в целях диагностики отдельных параметров или их совокупности; информация о состоянии биологического объекта представляется как качественными величинами, так и процессами (уровень артериального давления, температура тела, ЭКГ, ЭЭГ и др.); в практической медицине чаще всего применяется одномоментное точечное измерение того или иного параметра; многие изучаемые физиологические параметры человека не имеют конкретного, выраженного в физических величинах определения.

Лекция 1 « Методологическая основа моделирования » Процесс (от лат. processus – продвижение) – последовательная смена состояний объекта во времени. Природа объекта может быть произвольной: материальный (природный или искусственный) или идеальный (понятие, теория и т.п.), соответственно объект порождает материальный или идеальный процесс лат. состояний

Лекция 1 « Методологическая основа моделирования » Биологический процесс – это процесс, происходящий в живом организме. Биологические процессы состоят из ряда биохимических реакций или других событий, которые проводят к изменениям биообъектов.

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

Лекция 1 « Методологическая основа моделирования » Аналогией называют суждение о каком- либо частном сходстве объектов, причем такое сходство может быть существенным или несущественным. Существенность сходства (различия) зависит от уровня абстрагирования и в общем случае определяется конечной целью исследования.

Лекция 1 « Методологическая основа моделирования » Гипотезы и аналогии должны обладать наглядностью и сводиться к удобным для исследователя логическим схемам, которые называются моделями.

Лекция 1 « Методологическая основа моделирования » Модель (лат. modulus – мера) – это объект, заменяющий оригинальный объект с целью изучения некоторых свойств оригинала. Моделированием называется замещение одного объекта другим с целью получения информации о свойствах объекта- оригинала с помощью объекта-модели путем проведения экспериментов с моделью

Лекция 1 «Методологическая основа моделирования» ТЕОРИЯ ПОДОБИЯ - учение об условиях подобия физических явлений. Опирается на учение о размерности физических величин.. Планирование эксперимента (ПЭ) комплекс мероприятий, направленных на эффективную постановку опытов. Основная цель ПЭ достижение максимальной точности измерений при минимальном количестве опытов и сохранении статистической достоверности результатов

Лекция 1 « Методологическая основа моделирования » Классификация видов моделирования

Лекция 1 « Методологическая основа моделирования » Классификация математических моделей При синтезе и анализе систем, объектов и процессов используют различные виды математических моделей в зависимости от уровня иерархии, степени декомпозиции системы, аспекта, стадии и этапов синтеза и анализа. На любом уровне иерархии объект представляют в виде некоторой системы, состоящей из элементов. В этой связи различают математические модели элементов и систем

Лекция 1 « Методологическая основа моделирования » В общем случае уравнения математической модели связывают физические величины, которые характеризуют состояние объекта. Такими величинами являются: скорости и силы – в механических системах; расходы и давления – в гидравлических и пневматических системах; температуры и тепловые потоки – в тепловых системах; токи и напряжения – в электрических системах

Лекция 1 « Методологическая основа моделирования » Величины, характеризующие состояние технического объекта в процессе его функционирования, называют фазовыми переменными (фазовыми координатами). Вектор фазовых переменных задает точку в пространстве, называемом фазовым пространством. Фазовое пространство, в отличие от геометрического, многомерное. Его размерность определяется количеством используемых фазовых координат. Обычно в уравнениях математической модели фигурируют не все фазовые переменные, а только часть из них, достаточная для однозначной идентификации состояния объекта. Такие фазовые переменные называют базисными координатами..

Лекция 1 « Методологическая основа моделирования ». К математическим моделям предъявляются требования адекватности, экономичности и универсальности. Модель считается адекватной, если отражает исследуемые свойства с приемлемой точностью. Точность оценивается степенью совпадения предсказанных в процессе вычислительного эксперимента на модели значений выходных параметров с истинными их значениями.

Лекция 1 « Методологическая основа моделирования » Погрешность модели ε по всей совокупности т учитываемых выходных параметров оценивается одной из норм вектора εM =(ε1, ε2,…, εm) или где εj – относительная погрешность модели по j- му выходному параметру.

Лекция 1 « Методологическая основа моделирования » ε j = (ŷ j -y j )/y j ŷ j, где ŷ j значение j-го выходного параметра, полученное в результате вычислительного эксперимента; y j – значение того же параметра, принятое за истинное значение.

Лекция 1 « Методологическая основа моделирования » Классификация математических моделей.

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

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

Лекция 1 « Методологическая основа моделирования » На макроуровне объект рассматривают как динамическую систему с сосредоточенными параметрами. Математические модели макроуровня представляют собой системы обыкновенных дифференциальных уравнений. Эти модели используют при определении параметров объекта и его функциональных элементов.

Лекция 1 « Методологическая основа моделирования » На микроуровне объект представляется как сплошная среда с распределенными параметрами. Для описания процессов функционирования таких объектов используют дифференциальные уравнения в частных производных

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

Лекция 1 « Методологическая основа моделирования » Структурные модели имеют форму таблиц, матриц и графов. Наиболее перспективно применение древовидных графов типа И-ИЛИ-дерева. Такие модели наиболее широко используют на мета уровне при выборе принципа функционирования.

Лекция 1 « Методологическая основа моделирования » Морфологическое И/ИЛИ-дерево Морфологическое И/ИЛИ - дерево является одним из способов представления модели морфологического множества уровня идентификации. Морфологическое дерево имеет два типа вершин: ИЛИ - вершины, представляющие собой классификационные признаки и И - вершины, представляющие агрегацию подсистем. ИЛИ - вершины обычно изображают не закрашенными (белым цветом), а И - вершины изображают черными.

Лекция 1 « Методологическая основа моделирования » Морфологическое И/ИЛИ -дерево

Лекция 1 « Методологическая основа моделирования » Для анализа иерархической структуры могут применять теорию графов. Это позволяет перейти от графической модели к математической, в которой описание ведется по уравнениям, аналогичным законам Кирхгофа в электротехнике или уравнениям гидравлики.теорию графов модели

Лекция 1 « Методологическая основа моделирования » Иерархическая структура часто изображается в виде дерева, то есть графа без замкнутых маршрутов, с расположением вершин по определенным уровням, например, как показано на рис.1. Вершина верхнего уровня (на рисунке 0) называется корнем.

Лекция 1 « Методологическая основа моделирования » Граф, представленный на рис.1, соответствует И- дереву: вершины, которые расположены на одинаковых уровнях, являются обязательными элементами вышерасположенных систем. Так, для вершины 0.1 обязательные элементы 1.1, 1.2, а для вершины , 3.2 и 3.3. Например, сердечно-сосудистая система состоит из сердца, И аорты, И артерий, И артериол, И капилляров.И- дереву

Лекция 1 « Методологическая основа моделирования » Наряду с И-деревом используют ИЛИ- дерево, в котором на одинаковых уровнях располагаются вершины возможных элементов структур, их варианты. Например, прибор для измерения давления может иметь ртутный капилляр, ИЛИ трубку Бурдона, ИЛИ мембрану.

Лекция 1 « Методологическая основа моделирования » Часто применяют И-ИЛИ- дерево, которое соединяет уровни с обязательными элементами структуры с уровнями вариантов всех или части этих элементов (рис.2). Сочетание И- и ИЛИ- уровней может быть произвольным и не обязательно они должны чередоваться.

Лекция 2. «Применение теории графов в моделировании» Первая работа по графам была опубликована Леонардом Эйлером в 1736 г. Она содержала решение задачи о кенигсбергских мостах (рис), условие которой следующее: можно ли совершить прогулку таким образом, чтобы, выйдя из любого места, вернуться в него, пройдя один раз по каждому мосту? Эйлер обозначил части суши точками а, b, с, d – вершинами графа, а мосты, которые связывают части города ребрами, соединяющими соответствующие вершины. В результате получился граф, Эйлер дал отрицательный ответ на поставленный вопрос и доказал, что подобный маршрут имеется только для такого графа, каждая из вершин которого связана с четным числом ребер.

Лекция 2. «Применение теории графов в моделировании» Многие задачи моделирования сводятся к рассмотрению совокупности объектов, свойства которых описываются связями между ними. Интерес могут представлять различные связи и отношения между событиями, состояниями, функциональными органами биологической системы и т.п. В подобных задачах удобно изображать рассматриваемые объекты точками, называемыми вершинами, а связи между ними – линиями, называемыми ребрами. Множество вершин V, связи между которыми определены множеством ребер E, называют графом и обозначают G = (V, Е)

Лекция 2. «Применение теории графов в моделировании» Ориентированные графы Часто связи между объектами характеризуются определенной ориентацией. Например, в соединительных проводах электрической цепи, отношения между людьми могут определяться подчиненностью. Ориентированные связи характеризуют переход системы из одного состояния в другое. Для указания направления связи между вершинами графа соответствующее ребро отмечается стрелкой. Ориентированное таким образом ребро называют дугой, а граф с ориентированными ребрами – ориентированным графом или орграфом

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

Лекция 2. «Применение теории графов в моделировании» Типы конечных графов Если множество вершин графа, конечно, то он называется конечным графом. Конечный граф G = (V, E), содержащий р вершин и q ребер, называется (p, q)- графом. Пусть V ={v1, v2, …, vp} и Е ={e1, e2, …, eq} соответственно множества вершин и ребер (р, q)-графа. Каждое ребро еk E соединяет пару вершин vi, vj V, которые называются граничными. Ребро, граничными вершинами которого является одна и та же вершина, называется петлей. Ребра с одинаковыми граничными вершинами являются параллельными и называются кратными. Граф может содержать и изолированные вершины.

Лекция 2. «Применение теории графов в моделировании» Число ребер, связанных с вершиной vi (петля учитывается дважды), называют степенью вершины и обозначают через δ(vi) или deg (vi). Степень изолированной вершины равна нулю: δ(v4) = 0. Вершина степени единицы называется концевой или висячей вершиной. В любом графе сумма степеней всех вершин равна удвоенному числу ребер, а число вершин нечетной степени всегда четно.

Лекция 2. «Применение теории графов в моделировании» В орграфе различают положительные δ+(vi) и отрицательные δ-(vi)степени вершин. Суммы положительных и отрицательных степеней всех вершин орграфа равны между собой и равны также числу всех дуг. Крайние три вершины графа на рисунке имеют δ + (v)=2 и δ - (v)=1, центральная вершина – δ - (v)=3. Суммы положительных и отрицательных степеней всех вершин орграфа равны между собой и равны также числу всех дуг

Лекция 2. «Применение теории графов в моделировании» Граф, степени всех вершин которого одинаковы и равны r, называется однородным (регулярным) r-й степени. Полный граф с п вершинами всегда однородный степени п– 1, а пустой граф – однородный степени 0. Граф третьей степени называют кубическим. Он всегда имеет четное число вершин.

Лекция 2. «Применение теории графов в моделировании» Граф без петель и кратных ребер называют простым или обыкновенным (а). Граф без петель, но с кратными ребрами называют мультиграфом. Наиболее общий случай графа, когда допускаются петли и кратные ребра, называют псевдографом (б).. Если граф не имеет ребер (Е = 0), то все его вершины изолированы (V0), и он называется пустым или нуль-графом (в).

Лекция 2. «Применение теории графов в моделировании» Две вершины vi, vj V графа G = (V, Е) называются смежными, если они являются граничными вершинами ребра. Отношение смежности на множестве вершин графа можно определить, представив каждое ребро как пару смежных вершин, т.е. e k = (v i, v j ),i, j = 1, 2,..., q.

Лекция 2. «Применение теории графов в моделировании» Для неориентированных графов такие пары неупорядочены - e k = (v i, v j )=(v i v j ), а для орграфов – упорядочены, v i и v j означают соответственно начальную и конечную вершины дуги e k. Петля при вершине v i в обоих случаях представляется неупорядоченной парой (v i, v i ).

Лекция 2. «Применение теории графов в моделировании» Отношение смежности можно представить матрицей смежности. Строки и столбцы матрицы смежности соответствуют вершинам графа. Матрица смежности неориентированного графа симметрична, а орграфа – в общем случае несимметрична. В столбцах и строках, соответствующих изолированным вершинам, все элементы равны нулю. Элементы матрицы простого графа равны 0 или 1, а элементы главной диагонали матрицы нулевые

Лекция 2. «Применение теории графов в моделировании» Например, для графа, приведенного ниже имеем следующую матрицу смежности. В крайних столбце и строке указаны степени вершин abcddeg a02013 b20215 c02013 d

Лекция 2. «Применение теории графов в моделировании» Инцидентность элементов графов Если вершина vi, является концом ребра ek, то они инцидентны друг другу. Для орграфов различают положительную инцидентность, если дуга исходит из вершины, и отрицательную инцидентность, если дуга заходит в вершину.

Лекция 2. «Применение теории графов в моделировании» Рассматривая инцидентность вершин и ребер (р, q)-графа, можно представить его матрицей инцидентности размером р×q, строки которой соответствуют вершинам, а столбы – ребрам. e1e2e3e4e5e6e7 a b c d

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

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

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

Лекция 2. «Применение теории графов в моделировании» Графы Понтрягина–Куратовского являются непланарными. Полный пятиугольник а) представляет собой простой неплоский граф с минимальным числом вершин. Двудольный граф является моделью известной задачи о трех домах и трех колодцах.

Лекция 2. «Применение теории графов в моделировании» Строго доказывается, что граф является плоским только тогда, когда он не содержит подграфа, гомеоморфного одному из графов Понтрягина– Куратовского.

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

Лекция 2. «Применение теории графов в моделировании» Части графа Граф G = (V, Е) является частью графа G = (V, E), если Часть графа, которая наряду с некоторым подмножеством ребер графа содержит все вершины графа, называется суграфом.

Лекция 2. «Применение теории графов в моделировании» Деревья и лес Связные ациклические графы называются деревьями. Дерево на множестве р вершин всегда содержит q=p-1 ребер. При удалении хотя бы одного ребра дерево распадается на компоненты, каждая из которой представляет собой дерево или изолированную вершину.

Лекция 2. «Применение теории графов в моделировании» Маршрут длины т определяется как последовательность т ребер графа не обязательно различных, таких, что вершины двух соседних ребер совпадают. Например, (е 1, e3, е 4, е 6). Замкнутый маршрут приводит в начальную вершину, например (е 1, e4, е 3), он проходит через вершины (а, b, d).

Лекция 2. «Применение теории графов в моделировании» Маршрут, все ребра которого различны, называется цепью, а маршрут, для которого различны все вершины, называется простой цепью. Замкнутая цепь называется циклом, а простая цепь – простым циклом. Так, на графе (е 1, e4, е 7) - простая цепь, а (е 1, e4, е 3) – простой цикл

Лекция 2. «Применение теории графов в моделировании» Цикл, который содержит все ребра графа, называется эйлеровым циклом (задача о кенигсбергских мостах), а граф, в котором имеется такой цикл, называется эйлеровым графом

Лекция 3. Полюсные графы В зависимости от числа полюсов различают двухполюсные и многополюсные компоненты, которые называют соответственно двухполюсниками и многополюсниками. На риунке показано соединение двух трехполюсников (A, B), четырехполюсника (С) и трех двухполюсников (D, E, F). :

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

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

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

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

Лекция 3. Полюсные графы Если уравнение двухполюсника представимо в явном виде относительно поперечной переменной η= f y (ξ), то переменную η можно рассматривать как реакцию на воздействие ξ. Аналогично, если уравнение двухполюсника представлено в виде ξ = f z (η), то а переменную ξ можно рассматривать как реакцию на воздействие η. Двухполюсники, допускающие описание относительно обеих переменных, называются взаимоопределенными

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

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

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

Лекция 3. Полюсные графы Идеальные электрические двухполюсники: а – резистор; б – конденсатор; в – катушка индуктивности; г – источник напряжения; д – источник тока

Лекция 3. Полюсные графы Емкость – компонент, накапливающий электрическую энергию. Заряд q(t) емкости связан с напряжением u c (t) соотношением q(t) = Cu c (f), где С - емкость. Ток i c (t), протекающий через емкость, выражается как производная заряда по времени, следовательно: где S = 1/C называется инверсной емкостью.

Индуктивность – компонент, накапливающий магнитную энергию. Магнитный поток ψ(t) линейной индуктивности пропорционален протекающему в ней току i L (t) Напряжение uL(t) на индуктивности : где L – индуктивность; Г= 1/L - инверсная индуктивностью

Лекция 3. Полюсные графы Источники энергии в электрических цепях представляются идеальными двухполюсниками двух типов: Источник напряжения – это двухполюсник напряжение в котором определяется функцией времени e(t) и не зависит от протекающего по нему тока, т.е. u E (t) = e(t). Источник тока – это двухполюсник, ток в котором определяется функцией времени j(t) и не зависит от приложенного напряжения, т.е. i J (t)= j(t).

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

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

Лекция 3. Полюсные графы Уравнения связей выражаются законами Кирхгофа, представляющими условие непрерывности для токов и условие равновесия для напряжений в любой момент времени t: алгебраическая сумма токов для любой вершины равна нулю Σi(t)=0 алгебраическая сумма напряжений в любом контуре равна нулю Σu(t)=0

Лекция 3. Полюсные графы Механические поступательные системы Идеальные пассивные двухполюсники механических систем – это механическое сопротивление, масса и упругость. Перемещение x(t) и скорость v(t) являются продольными переменными, а сила f(t) – поперечной переменной.

Лекция 3. Полюсные графы Идеальные механические двухполюсники: а – сопротивление; б – масса; в – упругость; г – источник скорости; д – источник силы

Лекция 3. Полюсные графы Механическое сопротивление - компонент, который отражает превращение механической энергии в тепло. В простейшем случае это превращение происходит в результате вязкого трения, сила которого f B (t) пропорциональна относительной скорости v(t) трущихся тел. 1/B - механическое поступательное сопротивление, B – податливость

Лекция 3. Полюсные графы Масса– компонент, накапливающий кинетическую энергию и обладающий механической инерцией. Зависимость между силой инерции f M (t) и перемещением x M {t) или скоростью v(f) массы M относительно выбранной точки отсчета выражается уравнениями:

Лекция 3. Полюсные графы Упругость – компонент, накапливающий потенциальную энергию. Этот двухполюсник можно представить как пружину. Предполагается, что пружина не обладает массой и сила f K (t) реакции пропорциональна относительному перемещению x K (t) ее концов: К – жесткость пружины; 1/К - гибкость.

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

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

Лекция 3. Полюсные графы Уравнения связей механической поступательной системы выражают условие равновесия сил и условие непрерывности для скоростей (или перемещений): алгебраическая сумма сил для любой вершины равна нулю (принцип Даламбера): Σf(t)=0; алгебраическая сумма скоростей (перемещений) в любом контуре равна нулю: Σv(t) = 0

Лекция 3. Полюсные графы Механические вращательные системы Соотношения для механических вращательных систем аналогичны соотношениям для поступательных систем. Перемещению x(t) соответствует угол поворота φ(t), линейной скорости v(t) – угловая скорость ω(t), силе f(t) – вращающий момент μ(t). Для обозначения компонент механических вращательных систем можно использовать те же символы, что и для поступательных систем.

Лекция 3. Полюсные графы Вращательное сопротивление характеризует рассеивание механической энергии в тепло за счет вязкого трения, где B – крутильное сопротивление; 1/B – инверсное сопротивление.

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

Лекция 3. Полюсные графы Вращательная упругость – компонент, накапливающий потенциальную энергию вращательного движения, где К – крутильная жесткость; 1/К– гибкость.

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

Лекция 3. Полюсные графы Пример построения схемы и графа механической вращательной системы. Узлы графа соответствуют вращающимся массам, а направления ребер принимаются в соответствии с выбранным положительным направлением отсчета угла поворота. Параметры J1 и J2 означают моменты инерции роторов, В1 и В2 – вязкое трение в опорах, a K- жесткость вала.

Лекция 3. Полюсные графы Пневматические и гидравлические системы Движение газа или жидкости в ограниченной среде характеризуется зависимостью между давлением p(t) и потоком g(t), который выражается как количество молекул, проходящих в единицу времени. При этом поток рассматривается как поперечная величина, а давление (разность давлений) – как продольная величина.

Лекция 3. Полюсные графы Сопротивление – двухполюсник, учитывающий рассеивание энергии за счет вязкого трения. Его уравнение может быть представлено в одной из двух форм. G и R называются гидравлическими проводимостью и сопротивлением G=R -1, R=G -1. Параметр p R (t) - это разность давлений на концах двухполюсника при потоке g R (t).

Лекция 3. Полюсные графы Идеальные пневматические (гидравлические) двухполюсники: а – сопротивление; б – инертность; в – упругость; г – источник давления; д – источник потока

Лекция 3. Полюсные графы Инертность является двухполюсником, который характеризует противодействие изменению потока газа или жидкости в среде и описывается соотношениями. где L – пневматическая (гидравлическая) инертность. Величина p L (t) представляет собой разность давлений на концах этого двухполюсника при потоке g L (t).

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

Лекция 3. Полюсные графы Источники энергии в пневматических системах представляются идеальными двухполюсниками двух типов: источником давления и источником потока, которые определяются соответственно задающими давлением p(t) и потоком g(t).

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

Лекция 3. Полюсные графы Физическая система Аналогии между компонентами Поперечные Продольные Электрическая Ток i(t)Заряд q(t) Напряжение u(t) Потокосцепление ψ(t) Механическая поступательная Сила Импульс силы Скорость Перемещение Механическая вращательная Момент Импульс момента Угловая скорость Угол поворота Гидравлическая Объемный поток Объем жидкости Разность давлений Импульс давления Тепловая Теплоотдача Количество тепла Разность температур

Лекция 3. Полюсные графы Физическая система Аналогии между компонентами Электрическая Сопротивление RЕмкость СИндуктивность L Механическая поступательная Инверсное сопротивление Масса Упругость Механическая вращательная Инверсное сопротивление Момент инерции Вращательная упругость Гидравлическая Гидравлическое сопротивление Гидравлическая емкость Гидравлическая индуктивность Тепловая проводимость Теплоемкость

При построении математической модели физической системы исходные данные должны содержать сведения о структуре системы и свойствах входящих в нее компонент Исходные данные Топологические уравнения Компонентные уравнения Преобразование уравнений Однородная система координат Неоднородная система координат Сокращенная система координат Уравнения сечений Уравнения контуров Уравнения переменных состояния

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

Лекция 3. Полюсные графы В результате преобразования топологических и компонентных уравнений получаем систему уравнений, которую можно представить в матричной форме следующим образом:WX = QF. Матрица W и матрица Q, элементы которых выражаются через параметры компонент и интегро-дифференциальные операторы, определяют систему уравнений относительно вектора переменных X. Вектор F содержит заданные функции, характеризующие независимые источники.

Лекция 3. Полюсные графы Уравнение WX = QF в дифференциальной форме может быть преобразовано в уравнение переменных состояния, которое для линейной системы имеет вид где х – вектор переменных состояния и v – задающий вектор, связанные между собой матрицами A и B.

Лекция 4 Экспериментально- статистическое моделирование Объектами исследования - это реальные физические объекты (схемы, устройства, установки и т.д.), а также процессы, протекающие в них. Исследуемый объект (процесс) может быть представлен "черным ящиком«. Хi – входные управляемые параметры; Zi – входные контролируемые, но неуправляемые воздействия; Wi – неуправляемые и неконтролируемые переменные, Yi – выходные параметры,.

Лекция 4. Экспериментально- статистическое моделирование Входные параметры X и Z называют факторами, выходные параметры Y – откликами. Факторы могут быть количественными и качественными. Качественные факторы можно сделать количественными, построив условную порядковую шкалу. Наиболее простой пример такого кодирования - оценка качества по пятибалльной системе.

Лекция 4. Экспериментально- статистическое моделирование Уровень фактора – это его фиксированное значение. Нижний уровень фактора – это минимальное значение фактора, Хmin. Верхний уровень фактора – это максимальное значение фактора, Хmах. Основной уровень фактора – это его среднеарифметическое значение.

Лекция 4. Экспериментально- статистическое моделирование Интервал варьирования фактора – абсолютная величина равная половине диапазона изменения фактора: Факторное пространство –это пространство контролируемых переменных. Функция отклика – математическая зависимость отклика от рассматриваемых факторов. Поверхность отклика - пространство выходных переменных, т.е. геометрическое представление функции отклика.

Лекция 4. Экспериментально- статистическое моделирование Эксперимент - система операций, воздействий и (или) наблюдений, направленных на получение информации об объекте при исследованиях. Целью эксперимента является выявление свойств объектов, проверка справедливости гипотез для дальнейшего изучения объекта. Активный эксперимент - эксперимент, при котором предполагается вмешательство в процесс и возможность выбора в каждом опыте тех уровней факторов Х, которые представляют интерес. Пассивный эксперимент – это эксперимент, при котором исследователь не может воздействовать на изучаемый объект. Задача пассивного эксперимента – регистрация изменений факторов Z и откликов Y.

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

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

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

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

Лекция 4. Экспериментально- статистическое моделирование Полиномиальные модели Если вид функции отклика неизвестен, то она заменяется полиномом. Наибольшее распространение получили полиномиальные модели вида: для однофакторного эксперимента: полином первой степени полином второй степени полином m-й степени

Лекция 4. Экспериментально- статистическое моделирование для двухфакторного эксперимента (число факторов n=2) полином первой степени полином второй степени для многофакторного эксперимента (число факторов n) - полином второй степени

Лекция 4. Экспериментально- статистическое моделирование К факторам предъявляются следующие требования: Область определения факторов должна быть дискретной и ограниченной, каждый фактор может иметь три значения. Интервал варьирования, обычно в исследованиях интервал варьирования берется как 10…50 % от основного уровня фактора. Факторы должны быть управляемыми и однозначно определяемыми. Для каждого фактора должна быть возможность установки и поддержания требуемого уровня в течение всего опыта и возможность изменения уровня фактора в соответствии с планом эксперимента.

Лекция 4. Экспериментально- статистическое моделирование Факторы должны иметь непосредственное воздействие на исследуемый объект (процесс) и не являтся функцией других факторов. Факторы не должны быть коррелированны между собой. Между факторами недолжно существовать линейной связи. Факторы должны быть совместимы, т.е. при любом их сочетании, исследуемый объект должен сохранять работоспособность. Целесообразно располагать факторы в порядке убывания или возрастания их ожидаемого влияния на выходной параметр - f>fmin, где f –степень ожидаемого влияния на отклик.

Лекция 4. Экспериментально- статистическое моделирование Для упрощения расчетов коэффициентов уравнения регрессии, а следовательно, и построения модели исследуемого объекта (процесса), производят операцию кодирования переменных (переход от физической величины к ее изображению в относительных единицах) по формуле, где – x i кодированное текущее значение переменной (фактора); - X i текущее значение переменной в натуральных величинах; X i 0 - значение переменной на основном уровне в натуральных величинах; - h i интервал варьирования

Лекция 4. Экспериментально- статистическое моделирование Экспериментальные планы, по которым проводятся двухуровневые эксперименты, называются планами типа 2n, где n - количество факторов. А эксперимент, в котором реализуются все возможные комбинации уровней всех факторов, называется полным факторным экспериментом. Для двухуровневого эксперимента число всех возможных сочетаний факторов N=2n. Для полного факторного эксперимента принято обозначение ПФЭ 2 n. ПФЭ 2 n применяется для линейных моделей.

Экспериментально-статистическое моделирование Матрица спектра плана однофакторного эксперимента предусматривает поочередное варьирование каждого из факторов, в то время как все остальные факторы стабилизированы на основном уровне. Предположим, что проводится однофакторный эксперимент, при котором диапазоны изменения факторов: температуры -20…..+20 С (фактор 1); влажности 50…80% (фактор 2); давления 1….1,2 МПа (фактор 3). Номер опыта X1X1 X2X2 X3X3 1min00 2max00 30min0 40max0 000min 000max

Экспериментально-статистическое моделирование Значения факторов в натуральных величинах будут равны Номер опыта X1X1 X2X2 X3X о С50%1 Мп 2 20 о С50%1 Мп 3-20 о С80%1 Мп 420 о С80%1 Мп 5-20 о С50%1.2 Мп 620 о С50%1.2 Мп 7-20 о С80%1.2 Мп 820 о С80%1.2 Мп

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

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

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

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

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

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

Моделирование и анализ динамических систем Нормальная форма Коши где U – вектор фазовых переменных; F – вектор функция; t – время; начальные условия U0=U| t=0.

Моделирование и анализ динамических систем Интегрирование системы ОДУ в форме Коши заключается в определении значений U(t) на интервале времени 0 – tкон при заданных начальных условиях U0. При решении этой задачи на интервале интегрирования 0 – tкон выделяется конечное значение точек, в которых определяется значение U. Интервал между соседними точками называется шагом численного интегрирования и равен hi=(ti+1 – ti).

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

Моделирование и анализ динамических систем Алгоритмы решения систем дифференциальных уравнений Задача Коши заключается в решении систем обыкновенных дифференциальных уравнений первого порядка, представляемых в виде: где j=1÷n – номер каждой зависимой переменной уj, х – независимая переменная

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

Моделирование и анализ динамических систем Метод Эйлера – Коши – простейший метод первого порядка для численного интегрирования дифференциальных уравнений. Он реализуется следующей рекуррентной формулой: где h – шаг интегрирования (приращение переменной х). Этот метод обладает большой погрешностью и имеет систематическое накопление ошибок. Погрешность метода R~(h2), т.е. пропорциональна h2.

Моделирование и анализ динамических систем Метод трапеций – одна из модификаций метода Эйлера второго порядка. Он реализуется применением на каждом шаге формулы где Kj1= hFj(xi,Yji), Kj2= hFj(xi+h,Yji+Kj1). Этот метод дает погрешность R~ (h3) и относится к общим методам Рунге – Кутта

Моделирование и анализ динамических систем Метод Рунге – Кутта четвертого порядка является наиболее распространенным методом решения систем ОДУ при шаге h = const. Его достоинством является высокая точность – погрешность R~ (h5) - и меньшая склонность к возникновению неустойчивости решения. Алгоритм реализации метода Рунге – Кутта заключается в циклических вычислениях Yj(i+1) на каждом i+1-м шаге по следующим формулам

Моделирование и анализ динамических систем Алгоритм реализации метода Рунге – Кутта

Моделирование и анализ динамических систем Автоматическое изменение шага в ходе решения систем дифференциальных уравнений необходимо, если решение требуется получить с заданной точностью. При высокой точности (погрешность ε = 1*10-3) и при решении в виде кривых с сильно различающейся крутизной автоматическое изменение шага обеспечивает уменьшение общего числа шагов в несколько раз, резко уменьшает вероятность возникновения числовой неустойчивости

Моделирование и анализ динамических систем Метод Рунге – Кутта с автоматическим изменением шага заключается в том, что после вычисления Yj(i+1) с шагом h все вычисления проводятся повторно с шагом h/2. Полученный результат Y*j(i+1) сравнивается с Yj(i+1). Если, вычисления продолжают с шагом h, в противном случае шаг уменьшают. Если это неравенство слишком сильное, шаг, напротив, увеличивают.

Моделирование и анализ динамических систем 1. Задаем число уравнений n, погрешность ε=E, начальный шаг интегрирования h и начальные значения x0, y10,y20, …, yn0. 2. С помощью пяти циклов с управляющей переменной j=1, 2,..., n вычисляем коэффициенты, значение Yj(i+1)и погрешность: 3. Проверяем выполнение условий. Если первое условие не выполняется, делим шаг h на 2 и повторяем вычисления с п.2, восстановив начальные значения Yj(i+1). Если это условие выполняется и выполняется второе условие, значения xi+1=xi+h и Yj(i+1) выводятся на печать. Если второе условие не выполняется, шаг h увеличивается вдвое и вычисления повторяются с п.2.

Моделирование и анализ динамических систем Инерционные методы измерения параметров линейного движения твердого тела основаны на измерении силы инерции F(t), которая пропорциональна массе m и ускорению тела a. На рисунке показана схема механической системы прибора, состоящей из груза 1 массой m, движущегося поступательно по направляющим 2, и пружины 3

Моделирование и анализ динамических систем В соответствии со вторым законом Ньютона получим дифференциальное уравнение, где f - коэффициент вязкого трения, H*c/m; k- жесткость пружины, H/m; y(t) - перемещение тела, м; F(T) - внешняя сила, приложенная к телу, H.

Моделирование и анализ динамических систем Коэффициент вязкого трения определяется по формуле f= S/, где - коэффициент динамической вязкости масла, Па с; S – площадь поверхности контакта, м 2, - толщина слоя масла, м.

Жесткость цилиндрической пружины определяется по формуле где G - модуль упругости второго рода, для сталей G =0.8*1011 Па, i - рабочее число витков пружины, i = iсв+0.5; iсв - число свободных витков пружины; r - средний радиус витка пружины, м. Моделирование и анализ динамических систем

Преобразуем уравнение (94), разделив каждый его член на жесткость пружины k, получим: где m/k=T2, T – постоянная времени системы, с; f/k= 2 T, - коэффициент затухания системы; F(t)/k=y0 - начальное перемещение тела под действием внешней силы. Полученное дифференциальное уравнение является уравнением второго порядка, и для численного решения его необходимо преобразовать в систему дифференциальных уравнений первого порядка.

Моделирование и анализ динамических систем Запишем уравнение (96) в явном виде обозначив, получим систему ОДУ

Моделирование и анализ динамических систем Текст программы mech.m: global T dzeta; % Глобальные переменные T=0.05; dzet=[ ]; % Коэффициенты демпфирования y0=[5 0]; % Начальные условия t0=0.0; tfin=40*T; tspan=[t0 tfin]; % Временной интервал моделирования for i=1:3; dzeta=dzet(i); tspan, y0); % Решение ОДУ figure(1); plot(t,y(:,1),'k') % Графики переходных процессов hold on end grid; title('Transient response of the spring-mass-damper system') xlabel('Time, c.'); ylabel('Mass translations amplitude, mm') hold off

Моделирование и анализ динамических систем Система ОДУ записывается в виде отдельного m-файла, к которому осуществляется обращение из головной программы: % m- файл (odefun.m), описывающей систему ОДУ function ypr=odefun(t,y); global T dzeta % Система ОДУ ypr=[y(2); ((1e-3/T^2)-(y(2).*(2*dzeta/T))-y(1)./T^2))];

Моделирование и анализ динамических систем Графики переходных процессов

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

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

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

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

Моделирование и анализ статических процессов и систем Метод Гаусса, или метод последовательного исключения неизвестных, основан на приведении матрицы коэффициентов aij к треугольному виду. При этом алгоритм решения системы следующий [12]. 1. Двумя циклами с управляющими переменными i=1, 2,..., п и j = 1, 2,..., п вводим коэффициенты аij и bi, образующие массивы A (i, j)и B(i). 2. В трех циклах выполняем прямой ход исключения переменных путем преобразования коэффициентов по формулам: где i=1, 2,..., п - 1; j=i+l, i + 2, …, п и k = i+1, i + 2, …, п. В конце этих преобразований получаем хn = bn/аnn

Моделирование и анализ статических процессов и систем 3. Выполняем обратный ход, в котором определяются корни xn- 1, хn-2,..., x2, х 1, проводя вычисления по формулам: где i = п - 1, п – 2, …, 2, 1, j = i+1, i + 2,..., п. В результате формируется массив X (i) неизвестных xn-1, хn-2,..., x2, х Выводим значения корней из массива Х(i).

Моделирование и анализ статических процессов и систем Для системы (n = 3) 4x1+0.24x2 - 0,08x3 = 8; 0,09 х 1,+ 3x2 - 0,15x3 = 9; 0,04x1 - 0,08x2 + 4x3, = 20, используя систему MatLab, найдем значения корней: % Матрица коэффициентов A=[ ; ; ]; B=[8 9 20]'; % Вектор правых частей

Моделирование и анализ статических процессов и систем % Определение корней % Обратным делением матрицы A на вектор B X=A\B X = % Умножение инверсной матрицы A на вектор B X=inv(A)*B X =

Моделирование и анализ статических процессов и систем Если коэффициенты аii близки к 0, может наступить аварийный останов ЭВМ либо из-за погрешностей округления точность сильно ухудшается. Это является недостатком простейшей реализации метода Гаусса. Метод Гаусса с выбором главного элемента заключается в том, что при прямом ходе производится выбор наибольшего по модулю (главного) элемента и перестановка строк или столбцов. Последнее исключает деление на 0, если матрица А (i, j) содержит нулевые элементы. Для ПЭВМ, ведущих вычисления с числами с плавающей запятой, достаточен выбор А(i, i) 0.

Модели квазиодномерной гемодинамики Основные направления математического моделирования гемодинамики М.В. Абакумов, А.Я. Буничева, В.Б. Кошелев, С.И.Мухин, Н.В. Соснин, А.П. Фаворский, А.Б. Хруленко Факультет вычислительной математики и кибернетики МГУ. Опубликовал dodo.inm.ras.ru

Основные направления математического моделирования гемодинамики Моделирование течения крови в отдельном сосуде 2 D и 3D модели течения в крупных сосудах (уравнения Навье-Стокса) (моделирование упругости стенки сосуда, турбулентность течения, многокомпонентность крови (не ньютоновская жидкость), взаимодействие со стенками сосуда, области бифуркации сосудов, моделирование тромбообразования, стенозов, аневризмов и т.д.) 2D и 3D модели течения крови в мелких сосудах (капиллярах)с учетом реологии. Моделирование течения крови в сердце (2D и 3D модели) Моделирование течения крови в сети сосудов (дерево сосудов, замкнутая система) для исследования общих закономерностей течения крови. - На основе балансных соотношений - На аналогиях с «электрической цепью» - Квазиодномерное приближение: - соответствует типу сосудистой сети - дает возможность описать систему кровообращения в целом - позволяет отслеживать параметры течения крови вдоль сосуда - позволяет учесть особенности каждого сосуда - является основой для построения разномасштабных моделей -предоставляет возможность расчета переноса веществ кровью -предоставляет возможность использовать различные модели сосудов и органов -использует доступные физиологические данные -обладает хорошей точностью -предъявляет разумные требования к вычислительным мощностям

Комплексная нелокальная математическая модель сердечно-сосудистой системы Базовая модель описание течения крови в сосуде квазиодномерными уравнениями гемодинамики + нелокальные граничные условия + условия сопряжения в точках бифуркации Базовая модель описание течения крови в сосуде квазиодномерными уравнениями гемодинамики + нелокальные граничные условия + условия сопряжения в точках бифуркации 1. Создание математической модели течения крови в замкнутой системе сосудов (графе сосудов) произвольной топологии 2. Разработка эффективных моделей различных органов, сопряженных с работой сердечно-сосудистой системы (в том числе – точечной модели сердца ) 3. Создание эффективных однородных методов описания графа сосудов и численного решения глобальной математической модели S p Модели сосудов Ki dn ey Renal outflow q kid Модели почки Q t Модели сердца Однородная консервативная неявная разностная схема для системы уравнений на графе

Комплексная нелокальная математическая модель сердечно-сосудистой системы 6. Проведение вычислительных экспериментов в интересах фундаментальной и практической физиологии и медицины 5. Построение и анализ точных решений системы уравнений гемодинамики на графе 4. Создание интерактивного программного комплекса со средствами подготовки и обработки данных. Моделирование церебральной гемодинамики Моделирование влияния гравитационных нагрузок на сердечно-сосудистую систему g Моделирование почечной регуляции давления Моделирование влияния физических нагрузок... Моделирование переноса веществ кровью по графу сосудов с учетом процессов сорбции-десорбции

Адекватность свойств модели законам функционирования сердечно- сосудистой системы. Q t + + Kidne y Renal outflow q kid Программный комплекс CVSS позволяет с необходимой точностью рассчитывать гидродинамическую картину течения крови на графе, который физиологически адекватен сердечно-сосудистой системе человека, воспроизводить основные характеристики работы кровеносной системы Формализация системы кровообращения в граф Использование различных моделей элементов системы кровообращения

Обозначения Обозначения x u(t,x) L D(t,x), S=D 2 /4 Сосуд Локальная координата

x- координата вдоль оси сосуда, t - время, S(x,t) - площадь сечения сосуда, u(x,t) - скорость движения крови (направлена вдоль оси сосуда), p(x,t) - давление в крови, =const - плотность крови, L- длина сосуда, F(x,t) объемная плотность силы. Одиночный сосуд рассматриваем как трубку кругового сечения, протяженную по сравнению со своими поперечными размерами. Под эластичностью стенок понимается возможность изменения сечения сосуда под действием давления. X2X2 X1X1 S(x 2 ) S(x 1 ) X В основу описания движения крови в сосуде положены законы сохранения массы и импульса (количества движения). Интегральный вид уравнений S=S(p) – Эмпирическая зависимость площади сечения сосуда от давления. Система уравнений гемодинамики в квизиодномерном приближении Модель гемодинамики в одном сосуде

Последовательность математических моделей Модели кровеносных сосудов Сосуд Сосуд Жесткая трубка Эластичная трубка S=constS=const S const S=S(p) S=S(p,u,Q, … ) Диаметр сосуда может быть постоянным или переменным и может зависеть от любого числа физических и физиологических параметров, таких как давление, коэффициент эластичности, гравитация и т.д. Эту зависимость мы будем называть уравнение состояния. Кроме того, будем считать стенки сосуда тонкими. Q(S,p,u)=const

S p p min p max S min S max P S P min P max S min S max Характерный экспериментальный вид S(p) Простейшее приближение Нелинейное приближение Уравнение состояния

Типичное «уравнение состояния» сосуда S p p min p max S min S max В области нормальных значений давления зависимость площади поперечного сечения сосудов от давления близка к линейной. Простейшая математическая модель уравнения состояния S p p min p max S min S max Зависимость площади поперечного сечения от давления для большинства сосудов может моделироваться этим уравнением как в нормальном состоянии, так даже и в случае некоторых патологий. Заметим, что параметры S min, S max, P min, P max могут быть функциями времени t и пространственной переменной x. Например, стеноз, вызванный атеросклерозом, может быть промоделирован таким образом.

Пример специального «уравнения состояния» p0p0p0p0 pQ Q const p min p max p1p1p1p1 Эффект Остроумова-Балисса в церебральных артериях Этот эффект можно моделировать уравнением состояния следующего вида Необходимо иметь в виду, что различные формы уравнений состояния могут порождать специфические математические проблемы.

Сосуды - ребра графа Элементы модели 1. Система сосудов (вся сердечно- сосудистая система или ее часть) - граф сосудов. 2. Вершины графа : - области бифуркаций сосудов; - мышечные ткани; - отдельные органы. 2. Вершины графа : - области бифуркаций сосудов; - мышечные ткани; - отдельные органы. Почка, печень, кишечник, селезенка, легкие, … Ткань, мышцы вершины Области бифуркации сосудов моделируются законом сохранения потока вещества и условием непрерывности давления или интеграла Бернулли. Области фильтрации крови через ткани моделируются законом сохранения потока вещества и законом фильтрации Дарси.Области бифуркации сосудов моделируются законом сохранения потока вещества и условием непрерывности давления или интеграла Бернулли. Области фильтрации крови через ткани моделируются законом сохранения потока вещества и законом фильтрации Дарси. Каждый отдельный орган описывается специальной моделью, которая в простейшем случае является точечной.Каждый отдельный орган описывается специальной моделью, которая в простейшем случае является точечной. Модели вершин могут быть любого уровня сложности.Модели вершин могут быть любого уровня сложности.

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

Элементы модели 3. Сердце описывается двух или четырех камерной моделью. 3. Сердце описывается двух или четырех камерной моделью. легкие Большой круг кровообращения желудочек предсердие «Двухкамерная» модель сердца состоит из двух элементов- предсердия и желудочка и работает как насос. В течение систолы кровь из желудочка поступает в аорту. Этот процесс регулируется набором параметров: ударный объем, текущий объем предсердия и желудочка, давление в аорте и т.д. В течение диастолы желудочек наполняется. «Четырех камерная» модель объединяет две «двухкамерные» модели со своими параметрами. Пример «двухкамерной» модели сердца QVQV QAQA,, p mm Hg p top p bot t sys t dias t, сек

Модели сердца Неконсервативная модель График функции сердечного выброса приближает экспериментальную кривую потока или давления Q 0 t S t D t – продолжительность систолы – продолжительность диастолы Консервативная модель (V K – текущий объем крови в желудочке) Желудочек V K QAQA предсердие QVQV Регуляция по величине сердечного выброса по конечнодиастолическому объему V KD : V surg = K f V KD Регуляция по времени систолы и диастолы : V surg =const – ударный выброс – объемы желудочка в конце диастолы и систолы заданная параметрическая функция выброса вычисляемый текущий объем желудочка Поток крови Q V, поступающий в предсердие, определяется сердечным выбросом и потоком крови во всей системе. Такая модель позволяет сохранять объем циркулирующей крови, исследовать механизмы регуляции.

vv Математическая модель на графе ОБОЗНАЧЕНИЯ S(t,x) – площадь поперечного сечения u(t, x) - скорость потока крови p(t,x) - давление t - время x - локальная пространственная координата - плотность крови ( = const). F T – сила трения F T – сила тяжести 2. Каждой вершине графа, соответствующей области бифуркации сосудов, сопоставлено уравнение неразрывности и условие непрерывности интеграла Бернулли. 3. Каждой вершине графа, соответствующей тканям, сопоставлено уравнение сохранения вещества и уравнение фильтрации Дарси. 1. Каждому ребру графа сопоставлена система уравнений гемодинамики i,j –номера всех ребер, соединенных с каждой вершиной бифуркации, z i =±1 K d - коэффициент фильтрации

s = s ( p ) Свойства уравнений гемодинамики ( ГД ) Система уравнений гемодинамики ( ГД ) имеет гиперболический тип при условии, что для уравнения состояния выполнено неравенство dS/dp>0. В этом случае имеется две характеристики, два соотношения на характеристиках, скорость распространения малых возмущений. Соотношения на характеристиках Характеристики Скорость малых возмущений Так как давление в кровеносной системе мало отклоняется от своего среднего значения, в ряде случаев поведение системы удовлетворительно описывается линеаризованными уравнениями гемодинамики ( ЛГД ).

4. В граничных вершинах графа должны быть заданы краевые условия. Это, например, самосогласованная модель сердца или некоторые законы изменения потоков или давления. Граничные вершины Вход ( Q(t) ) Выход( p(t) ) Пример графа сосудов головного мозга 5. При моделировании процессов переноса растворенных в крови веществ в вершинах графа формулируются дополнительные уравнения, описывающие эти процессы.

Перенос вещества с массовой концентрацией C с учетом диффузии (D=const- коэффициент диффузии) по сосуду с переменным сечением описывается дифференциальным уравнением число Маха М

Уравнения гемодинамики на ребре i Уравнения сопряжения и граничные условия Неявная разностная схема, i=1,2,3,4 - весовые коэффициенты ( на каждом ребре) Неявная аппроксимация уравнений сопряжения и граничных условий a S,a u – коэффициенты искусственной вязкости i, I – коэффициенты искусственной дисперсии, например 2 = Численный алгоритм и программный комплекс Разностная схема апробирована на точных аналитических решениях

Линеаризация по Ньютону Линеаризованное разностное уравнение неразрывности в каждом внутреннем узле j дискретной сетки на каждом ребре графа Линеаризованное разностное уравнение движения

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

Структура программного комплекса CVSS (CardioVascular Simulating System) CVSS Pre-processor Расчет начальных данных База данных Начальные данные Контроль данных Метод 1 Метод 2 Решение в линейном приближении Блок контроля текущих результатов расчета Post-processor Solver Графический редактор Расчет переноса

Численные методы и алгоритмы 3. Разработана упрощенная конечно - разностная схема для решения уравнений гемодинамики. 1. Разработан специальный формат описания произвольного графа сосудов. 4. Полная нелинейная система разностных уравнений решается с использованием итерационных методов (метод Ньютона). 5. Линеаризованная система разностных уравнений решается с использованием прямых методов. 2. На каждом ребре графа использована однородная консервативная разностная схема второго порядка аппроксимации.

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

Общее решение ЛГД уравнений на любом i-ом ребре графа представляет собой суперпозицию двух произвольных волн, распространяющихся в противоположных направлениях (|u|

Коэффициент прохождения волны скорости через вершину графа из ребра j в ребро i Бегущие волны определяются следующей формулой: Они численно равны отношению амплитуды волн скорости и давления до и после их взаимодействия с вершинами графа сосудов. - «транспортные коэффициенты». Коэффициент отражения волны скорости от вершины графа на ребре i ji i

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

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

Граф артерий мозга Артерии головного мозга Моделирование церебрального кровообращения Построение графа для системы сосудов - первый этап моделирования. Пример графа сосудов головного мозга, включающий в себя артерии до уровня третьей бифуркации. Построение графа для системы сосудов - первый этап моделирования. Пример графа сосудов головного мозга, включающий в себя артерии до уровня третьей бифуркации. Ткань мозга Коллатерали Сердце Руки «Точечная» модель остальной части сосудистой системы Виллизиев круг Представленный граф включает в себя сердце, дугу аорты, позвоночные артерии, сонные артерии, Виллизиев круг, схематично представлены руки, артерии P1,P2,P3, A1,A2,A3,M1,M2,M3, коллатерали и некоторые другие артерии. Венозный возврат представлен схематично. Влияние всей остальной части кровеносной системы описывается «точечной» моделью.

Вычислительный эксперимент Пациент П. имеет 70% стеноз правой внутренней каротидной артерии и 90% стеноз левой каротидной артерии. Параметры сосудов его головного мозга (длины, диаметры, коэффициенты эластичности и т.д.) и сердца были установлены в ходе клинических исследований. В ходе операционного вмешательства предполагалась окклюзия (пережатие) ряда артерий (точки 2-9 на рисунке). Вопрос: Что произойдет с распределением потоков крови по различным отделам головного мозга при окклюзии? Пациент П. имеет 70% стеноз правой внутренней каротидной артерии и 90% стеноз левой каротидной артерии. Параметры сосудов его головного мозга (длины, диаметры, коэффициенты эластичности и т.д.) и сердца были установлены в ходе клинических исследований. В ходе операционного вмешательства предполагалась окклюзия (пережатие) ряда артерий (точки 2-9 на рисунке). Вопрос: Что произойдет с распределением потоков крови по различным отделам головного мозга при окклюзии? Результаты моделирования Допустимым является уменьшение притока крови не более чем на 20% First column is the flow without occlusions. Возможно ли это? В случае пациента П. компенсация происходит за счет коллатерального кровообращения. Потоки крови в коллатералях Распределение потоков крови сильно меняется после окклюзии Точки окклюзии

Вычислительный эксперимент Анализ изменений объема крови в венозной и артериальной части мозга для оценки ликвород инамики головного мозга по модели Келли. Изменение объема крови в сосудах выше Веллизиева круга в расчетах составляет 1,8 мл за один период сокращения сердца, что согласуется с экспериментальными данными в соответствии с которыми между головным и спинным мозгом за период сокращения сердца циркулирует около 1 мл ликвора. Граф церебральных сосудов + элементы большого круга кровообращения : - двухкамерное сердце - дуга аорты - артерии, вены, ткани рук - точечные сопротивления и обобщенные сосуды с соответствующими объемами и резистивными свойствами Модель замкнута Взаимовлияние давления в аорте и в головном мозге 11 22

Моделирование большого круга кровообращения сердце почки мозг ноги кишечник Граф большого круга Топология графа. Параметры сосудов (артерий, вен, резистивных сосудов, емкостных вен и т.д.) и сердца. Коэффициенты фильтрации (закон Дарси) для каждого органа. Сердце представлено само- согласованной двухкамерной моделью. aurical Ventrical V K QAQA QVQV Q 0 t S t D t Ударный объем сердца 85 ml, t s =0.3 s, t d =0.5 s. Ударный объем сердца 85 ml, t s =0.3 s, t d =0.5 s.

Квазипериодический режим в большом круге кровообращения Ударный объем сердца 85 мл, продолжительность систолыts=0.3 с, диастолы - td=0.5 с. Течение крови (после периодов работы сердца) в сердечно- сосудистой системе выходит на режим, при котором величины максимального и минимального давления в сосудах практически не меняются в течение длительного времени (сутки). Давление в аорте Давление в артерии верхней конечности Давление в вене верхней конечности Давление в резистивном сосуде брыжеечной артерии (в начале и конце сосуда)

Квазипериодический режим в большом круге кровообращения сосуда на графе Название сосудаОСК, численный расчет, мл/с ОСК, мед. данные * мл/с 2Чревный ствол 24,321,2 ± 2,5 56,57Подключичная артерия 3,61,15 - 5,26 50,49Позвоночная артерия 21,55 ± 0,55 Объемная скорость кровотока. (*Ультразвуковая доплеровская диагностика сосудистых заболеваний. Под ред. Ю.М. Никитина, А.И. Труханова, 1998) Построенная модель применена для исследования влияния таких факторов, как: параметры резистивных сосудов; продолжительность систолы; величина ударного объема сердца; величина коэффициента вязкости крови; физическая нагрузка; регуляторная функция почки; на артериальное давление и на кровоток во всей системе.

Изменение режима работы сердца Изменение тонуса сосудов Изменение наполнения тканей Изменение давления в системе и аорте Моделирование барорецепторной нейрогенной регуляции Барорецепторы

Повышение артериального давления Увеличение импульсации от барорецепторов Реакция ЦНС Снижение тонуса сосудов. Увеличение заполненности тканей кровью. Уменьшение частоты сокращений сердца. Механизм нейрорегуляции настроен на поддержание характерного давления в сосуде с барорецепторами, которые реагируют на отклонение,, где – среднее давление в сосуде. Понижение давления: обратные процессы Моделирование барорецепторной нейрогенной регуляции Принципиальная схема нейрогенной регуляции Модель изменения тонуса сосудов Повышение давления приводит к увеличению площади сечения и к уменьшению жесткости

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

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

Моделирование почечной регуляции давления Renal pressure Volume of CVS Простейшая модель почечной регуляции: если среднее почечное давление становится больше, чем некоторое p * поч, то начинается экспоненциально нарастающий (по сравнению с нормальным q * поч ) сброс жидкости Renal normal outflow Renal outflow Kidney Renal outflow q kid, мм.рт.ст А Выведение и потребление жидкости q, см 3 /с q * поч Выведение жидкости почкой выведение поступление

- Амплитуда пульсовой волны в i-ом сосуде в «норме» - Амплитуда для пораженной заболеванием сосудистой системы 1 - H=0%, 2 - H=25%, 3 - H=50% 4 - H=75%, 5 - H=90% 1 - R=1, 2 - R=0.5, 3 - R= R=0.01 Получены качественные и количественные зависимости симптоматики заболевания от степени поражения сосудистой системы. Моделирование неспецифического аортоартериита

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

Моделирование гравитационного воздействия g g Исследуется влияние гравитации на гемодинамику человека. Моделирование на полном графе (большой круг кровообращения + церебральное кровообращение) позволяет установить изменения в гемодинамике, например, в случае быстро меняющейся гравитации. Объем кровотока в мозге сильно падает под действием нарастающей гравитации. Меняется распределение потоков крови по отделам головного мозга.

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

Пространственный граф системы кровообращения человека в силовом поле