Содержание к диссертации
Введение
Глава 1 Современное состояние моделирования движения ньютоновских сред в сферических каналах 9
1.1 Основные математические модели движения ньютоновских вязких сред 9
1.2 Методы решения систем дифференциальных уравнений в частных производных 23
1.3 Аппроксимация и устойчивость конечно-разностных приближений модельных уравнений 28
1.4 Выводы. Цели и задачи исследования 33
Глава 2 Синтез математической модели движения ньютоновской сжимаемой жидкости в сферическом канале 35
2.1 Постановка задачи 35
2.2 Преобразование исходной системы уравнений к безразмерному виду 46
2.3 Численная схема решения системы уравнений движения вязкой сжимаемой жидкости и алгоритм решения 51
2.4 Устойчивость конечно-разностной схемы 63
Глава 3 Экспериментальное исследование основных гид родинамических параметров сферического канала 65
3.1 Численный анализ устойчивости математической модели 65
3.2 Результаты численного моделирования полей давления, скоростей и температуры 67
3.3 Установка для проведения физического эксперимента и его результаты 77
3.4 Сравнительный анализ вычислительного и физического экспериментов 84
Глава 4 Практическая реализация результатов исследования 87
4.1 Методика расчета пневмомеханического устройства со сферическим несущим элементом 87
4.2 Описание пакета предметно-ориентированного прикладного программного обеспечения 90
4.3 Технические решения 96
4.4 Расчет экономической эффективности внедрения разработанных пакетов программ и устройств 104
Основные выводы 106
Список использованных источников
- Методы решения систем дифференциальных уравнений в частных производных
- Преобразование исходной системы уравнений к безразмерному виду
- Результаты численного моделирования полей давления, скоростей и температуры
- Описание пакета предметно-ориентированного прикладного программного обеспечения
Введение к работе
Актуальность темы. Для проектирования новых и анализа существующих технических систем необходима разработка методов расчета явлений переноса, базирующихся на фундаментальных законах.
Достоверная информация о структуре полей скорости, давления и температуры при течении вязкой сжимаемой среды предопределяет рациональные режимы функционирования элементов в системах предметного назначения. Существующие подходы при синтезе многомерных математических моделей с распределенными параметрами, как правило, основывается на системе допущений, не позволяющей достаточно адекватно описывать реальную физическую картину. Например, пренебрежение фактором сжимаемости в дозвуковых течениях (число Маха М<1) и использование сжимаемой модели ньютоновской среды приводит к снижению прецизионности пневмомеханических элементов различных систем управления, неточному прогнозированию гидродинамической обстановки обтекания в узлах летательных аппаратов и т.д. Данное обстоятельство предполагает использование математических моделей, основывающихся на менее жестких допущениях. Задача осложняется тем, что разрешение такой ситуации приводит к необходимости решения нестационарных, сопряженных, существенно нелинейных уравнений Навье-Стокса и энергии в трехмерной постановке для областей сложной геометрии.
Классические аналитические методы неприменимы для решения задач в такой постановке, т.к. предназначены для анализа явлений, имеющих линейный характер, и поэтому на практике применяется численное интегрирование. Однако и численные методы, основывающиеся на технологии конечно-разностной аппроксимации уравнений и граничных условий, требуют
5 дополнительной адаптации для конкретных задач и минимизации использования ресурсов вычислительной техники для достижения необходимой точности.
В настоящее время подход, основывающийся на совместном применении численных методов и ЭВМ для решения разнообразных задач, трансформировался в новую современную технологию проведения теоретических исследований - вычислительный эксперимент. В связи с этим реализация эффективных численных методов и алгоритмов в виде компонентов проблемно-ориентированных программ для проведения вычислительных экспериментов актуальна.
В качестве предметной области, иллюстрирующей особенности использования разработанного инструментария, выбрана задача дозвукового течения вязкого газа в сферическом канале в элементах транспортных пневмосистем.
Диссертационная работа выполнялась на кафедре высшей математики Воронежской государственной технологической академии в соответствии с планом научно-исследовательских работ по теме: «Математическое обеспечение структурного и параметрического анализа технологических, технических и информационных систем» (№г.р. 01.20.0011235).
Цель работы. Разработка и реализация численного метода интегрирования уравнений движения вязкой сжимаемой среды в виде комплекса предметно-ориентированного программного обеспечения для проведения вычислительных экспериментов.
Для достижения указанной цели поставлены следующие задачи: проанализировать существующие математические модели и численные методы для проведения вычислительного эксперимента по определению основных гидродинамических характеристик потоков вязкой сжимаемой жидкости в сферических каналах; на основе уравнений Навье-Стокса синтезировать математическую модель осесимметричного напорного течения вязкой сжимаемой жидкости в сферическом канале и разработать численный метод интегрирования; провести вычислительный эксперимент по определению основных гидродинамических характеристик и сравнить полученные результаты с имеющимися данными и натурным экспериментом; разработать методику расчета и апробировать ее на системах пневмотранспорта.
Методы исследования. В работе использованы методы теории гидродинамических и теплообменных процессов, вычислительной математики и моделирования.
Научная новизна: математическая модель осесимметричного напорного движения вязкой сжимаемой среды в канале образованном концентрическими сферами, позволившая отказаться от допущений о несжимаемости среды в ламинарном режиме течения; явная конечно-разностная схема численного интегрирования нестационарной двумерной системы гиперболически-параболического типа с граничными условиями в виде Дирихле и фон Неймана, отличающаяся декомпозицией исходной задачи на последовательно решаемые несопряженные подзадачи гидродинамики и теплообмена; результаты вычислительных экспериментов по идентификации основных гидродинамических и тепловых характеристик течения вязкого сжимаемого газа в сферическом канале, их анализ и обобщение. Практическая значимость: алгоритм и программа расчета полей скорости, давления и температуры в сферическом канале при напорном движении вязкой сжимаемой среды; методика расчета для проектирования элементов систем со сферическими каналами;
7 - технические предложения по одновременному контролю массы и пространственного положения протекторного полотна в системах пневмотранспорта, защищенных патентами РФ №2194954 и №2183822.
Результаты работы апробированы на ОАО ПИТ "Воронежский шинный завод" и подтверждены соответствующей документацией.
Апробация работы. Основные положения диссертации докладывались и обсуждались в работе следующих научных форумов: «Математическое моделирование в технике и технологиях. ММТТ-12» (Великий Новгород, 1999); Восьмая научно-практическая конференция «Резиновая промышленность. Сырье, материалы, технология» (ГУЛ НИИШП, Москва, 2001); «Информационные технологии в управлении и моделировании» (Белгород, 2005); «Составляющие научно-технического прогресса» (Тамбов, 2005); «Качество науки - качество жизни» (Тамбов, 2005); «Глобальный научный потенциал» (Тамбов, 2005); «Достижения ученых XXI века» (Тамбов, 2005), а так нее отчетных научных конференциях Воронежской государственной технологической академии (Воронеж, 1999-2005).
Публикации. По теме диссертационной работы опубликовано 12 печатных работ, в том числе два патента Федерального института промышленной собственности.
Структура и объем работы. Диссертационная работа состоит из введения, четырех глав, основных выводов, списка литературы из 122 наименований. Материал диссертации изложен на 118 страницах, содержит 45 рисунков и 6 таблиц.
В первой главе рассмотрены основные математические модели движения ньютоновских вязких сред на основе обзора материалов отечественных и зарубежных литературных источников. Проанализированы методы решения систем дифференциальных уравнений в частных производных. Изучены вопросы, касающиеся устойчивости решения методом конечно-разностной
8 аппроксимации. Сформулированы задачи, решение которых способствует достижению поставленной цели исследования.
Во второй главе синтезирована математическая модель движения ньютоновской сжимаемой жидкости в сферическом канале. Разработана численная схема решения системы уравнений движения вязкой сжимаемой жидкости. Здесь же представлен расчетный алгоритм, реализующий расчетную схему. Выполнен анализ устойчивости предложенной конечно-разностной схемы.
Третья глава посвящена экспериментальным исследованиям. Численно проанализирована устойчивость программно реализованной математической модели. Подробно представлены результаты проведенных вычислительных экспериментов по определению гидро- и термодинамических параметров вязкой сжимаемой ньютоновской жидкости. Выполнена проверка адекватности математической модели на специально сконструированных экспериментальных установках. Проведен анализ правомерности допущений принятых при построении математической модели.
Практическая реализация результатов математического моделирования представлена в четвертой главе. Предложена методика и диаграмма по определению величины силы поддерживающей сферический элемент на слое сжимаемой среды. Представлено разработанное прикладное программное обеспечение для расчетов гидродинамических и теплофизических параметров вязкой сжимаемой жидкости. Описаны сконструированные пневмомеханические устройства со сферическими каналами. Показана экономическая эффективность внедрения разработанных пакетов программ и устройств.
В приложении представлены акт производственных испытаний разработанного оборудования, патенты РФ, дипломы выставок и научных конференций, свидетельствующие о практической ценности и полезности проведенных работ.
Методы решения систем дифференциальных уравнений в частных производных
Преимуществом данного модельного уравнения является наличие его аналитического решения. Бентон и Плацман получили 35 точных аналитических решений для данного одномерного уравнения [ПО]. Наличие точных решений уравнения (1.48) позволяет отрабатывать вычислительные алгоритмы и схемы.
Однако в случае двумерного нелинейного уравнения невозможно получить его аналитическое решение [91]. Поэтому для его решения необходимо применять численные методы.
Численные методы решения дифференциальных уравнений разделяются на релаксационные и маршевые. Первые применяются для решения внутренних задач, вторые для отыскания решений в незамкнутой области. Релаксационные методы применяются для решения внутренних задач, а маршевые для отыскания решений в незамкнутой области [4].
Получение конечно-разностных аналогов дифференциальных уравнений численно может быть реализован одним из четырех методов: 1. разложение выражения в ряд Тейлора; 2. полиноминальная аппроксимация; 3. интегральный метод; 4. метод контрольного объема.
Однако для аппроксимации модельных уравнений первые три метода могут являться безусловно несходящимися и неустойчивыми.
Достоинством метода контрольного объема заключается в получении некоторого осредненного значения, что положительно сказывается на устойчивости системы. Рассмотрение макроскопических физических зависимостей обеспечивает очевидное преимущество по сравнению с использованием математического аппарата непрерывных функция. Данный метод хорошо зарекомендовал себя для решения дифференциальных уравнений с течением невязких или сжимаемых сред [81].
Таким образом, в методе контрольного объема прослеживается фактическое выполнение физических законов, а не простое предположение о пределе приращения, стремящемся к нулю. Это и лежит в основе понятия консервативности конечно-разностного метода, то есть, обеспечивает выполнение определенных законов сохранения, справедливых для исходных дифференциальных уравнений. Такой подход позволяет избежать трудностей, связанных с аппроксимацией старших производных уравнений Навье-Стокса. Разделение вектора плотности на конвективную и вязкую составляющие дает возможность без труда использовать данную методику и для расчета движений идеальной среды [4, 81].
Консервативность предлагаемого метода позволяет по единому алгоритму осуществлять численное моделирование ряда сложных гидродинамических течений во всей области интегрирования, вплоть до достаточно больших чисел Рейнольдса.
Для рассматриваемого случая течения вязкого сжимаемого газа метод контрольного объема в переменных поля образует явную конечно-разностную схему, которая является условно-сходящейся и условно-монотонной. Схема в силу способа ее построения является консервативной по массе, по составляющей импульса и полной энергии.
Наиболее эффективным численным методом решения поставленной задачи является получение ее нестационарного решения. Интересующее нас стационарное решение получается как асимптотическое по времени [81].
Рассмотрим некоторые конкретные методы решения дифференциальных уравнений в частных производных и в частности модельного уравнения (1.48).
Наиболее простой и эффективной схемой решения уравнения Бюргерса является метод вперед по времени и центральных разностей по пространству (ВВЦП). В результате конечно-разностный аналог уравнения (1.48) принимает вид Л L + c-J« і± = ц і J— ±. (1.49) At 2Ax (Де)2 ;
Выражение (1.49) - явная одношаговая схема с погрешностью аппроксимации 0(At,(Ax)2). Для уравнения Бюргерса вводятся коэффициенты г = /jAt/(Ах)2, и v cAt/Ax. Используя эвристический анализ устойчивости [4] ограничение примет вид v2 2r. (1.50) Полезным параметром при решении модельного уравнения является сеточное число Рейнольдса, определяемое соотношением Re = cAx/v = v/r. (1.51) Этот безразмерный параметр характеризует отношение конвекции к диффузии. Условие устойчивости накладывает следующее ограничение на значение сеточного числа Рейнольдса 2v Reto 2/v. (1.52)
Важной характеристикой конечно-разностных схем является появление осцилляции (всплесков) решения. В работе [81] указывается, что критическим значением, при дискретном решении уравнений, является сеточное число Рейнольдса равное 2. Рост числа Re приводит к увеличению частоты разностного решения до частоты Найквиста, отвечающей минимально возможной длине волны Л = 2Ах, что для уравнения Бюргерса соответствует значению сеточного Re = 2. Поэтому на больших значениях числа Рейнольдса невозможно получить картину течения даже качественно за счет возникающей ошибки неразличимости. Однако выражение (1.52) показывает, что в действительности ограничение имеет другую форму записи. При значениях сеточного числа Рейнольдса немного больше 21 v возникают осцилляции, которые ведут к разрыву в решении.
Преобразование исходной системы уравнений к безразмерному виду
Получение численного решения системы уравнений с граничными условиями, проводилось заменой непрерывной области совокупностью изолированных точек. Данный метод, основанный на аппроксимации производных конечными разностями, позволяет перейти от системы дифференциальных уравнений в частных производных к системе алгебраических уравнений.
Из рис. 2.3 видно, что проекция рассматриваемой физической области решения может быть представлена в следующем виде (рис. 2.4).
В соответствии с принятой нормированной вычислительной областью решения модифицировались исходные безразмерные уравнения. Безразмерная радиальная координата вводится следующим образом: Л = - -, (2.77) П го где г0 г: - радиусы внутренней сферы и полусферической полости соответственно; г - текущий радиус. При этом значении безразмерного R радиуса находится в пределах от 0 до 1. Тогда r -(r]-r0) (R + a){r r0) (R + a)d, (2.78) r\-ro где а = rQ 1{гА - r0) - геометрический симплекс. Соответственно система модельных уравнений движения ньютоновской вязкой сжимаемой среды в сферическом канале (2.66)-(2.70) с граничными условиями (2.71)-(2.76) преобразована к нормированному виду, при этом значения координат изменяются в пределах 0 І? 1 и 0 1,
Стационарное решение отыскивалось из нестационарного решения задачи (2.79)-(2.89) маршевым методом. Интересующее нас стационарное решение получается как асимптотическое по времени [81]. В качестве маршевой координаты считалось время процесса.
В результате перехода от непрерывной области решения к ее дискретному аналогу, задача об отыскании величины давления P(R,,rj) сводится к определению вычислению значений P(i-AR,j-dkQ,n-h7j) в узлах сетки (рис. 2.6, 2.7) определяемых значениями величин i,j.
Для конечно-разностной аппроксимации значений давления в области решения использовалась трехточечной схемой по плоскости координат R и 0 (рис. 2.7), т.е. вычисление значения в точке i,j будем осуществлять по известным значениям в соседних точках і Таким образом, выбранная конечно-разностная схема аппроксимации исходной системы уравнений обеспечивает точность получения решения с первым порядком точности по времени и вторым порядком точности по пространственным координатам, или погрешностью в виде оГДт/ Дх Дл )2!.
Конечно-разностный аналог дифференциальных уравнений (2.79)-(2.83), с использованием трехточечной схемы с центральными разностями по координатам (2.90)-(2.92) и разностью вперед по времени (2.93) записывается
Вычисление значений на следующем временном слое осуществляется путем выражения соответствующее значение переменных из выражений (2.94), (2.95), (2.96) и (2.98)
Аппроксимация граничных условий второго рода выполняется «правыми» и «левыми» конечно-разностными аналогами, которые в общем случае для трехточечной схемы записываются в виде [14]: В соответствии с выражением (2.107) граничные условия второго рода
Система алгебраических уравнений (2.99)-(2.102) представляет собой конечно-разностный аналог исходной системы уравнений в частных производных, а граничные условия (2.103)-(2.106), (2.108)-(2.111) позволяют численно получить решение в узлах сетки дискретизированной физической области. Численную реализацию на ЭВМ осуществляется в следующей последовательности: а) вычисление значение плотности по выражению (2.101); б) расчет значений скоростей на основании уравнений (2.99) и (2.100); в) определение значений температуры по выражению (2.102).
В соответствии с данной последовательностью составлен алгоритм численного решения модельной системы. Блок схема разработанного алгоритма представлена на рис. 2.8.
Результаты численного моделирования полей давления, скоростей и температуры
Программная реализация синтезированной математической модели позволила провести ряд численных экспериментов для различных значений характерных геометрических и гидродинамических параметров. Типичные графические зависимости исследуемых параметров представлены на рис. 3.4-3.11.
Анализ полученных графических зависимостей показывает, что в силу малости зазора в сферическом канале давление во входном участке, по радиальной составляющей, для любого меридионального сечения успевает выровняться (рис. 3.4 а). Наблюдаются значительные потери давления во входной области потока связанные с разворотом течения и соответственно потерей кинетической энергии, в области по 0 от 0 до 0,1-Ю,4 (в зависимости от числа Маха). После этого участка давление стабилизируется, аналогично линейным каналам с постоянным сечением (рис. 3.4 а-3.11 а).
На рис. 3.9 видно, что в конце входного участка давление имеет дозвуковой скачок в узкой области по 0, с локальным минимумом и локальным максимумом в диапазоне 0,1-Ю,3, тем самым, увеличивая размеры входного участка. Полученный скачок имеет всего лишь одну колебательную амплитуду и быстро затухает после 0«О,4, с дальнейшей структурой изменения, как и для чисел Маха менее 0,2, такую ситуацию достаточно информативно отслеживает рис. 3.9 в.
Из анализа рисунков следует, неизотермическое течение в исследуемом диапазоне чисел Маха играет существенную роль только на входном участке, и поэтому при расчетах вдали от точки входа «0,3, можно приметь модель несжимаемого вязкой среды.
При числах Маха более 0,2 (рис. 3.10, 3.11) величина дозвукового скачка давления увеличивается, а сам он растягивается в меридиональном направлении. Область отрицательной радиальной скорости сужается (рис. 3.10 г) и появляется обратное течение в меридиональном направлении на поверхности внутренней сферы напротив точки входа.
Возникновение противотока на входе в канала изменяет поле температуры, которое становится существенно неоднородным, это иллюстрирует рис. 3.10 6. Учитывая, что при этом число Рейнольдса практически соответствует критическому значению перехода от ламинарного к турбулентному режиму течения, а следовательно границе применимости данной модели, тем не менее, в целом, получаются вполне корректные с точки зрения физических представлений результаты.
Проведенный вычислительный эксперимент с моделью для закритических значений чисел Рейнольдса (рис. 3.11) свидетельствует о правомерности применения данной модели в этой области. При этом изотермическое поле температур на входном участке приближается к адиабатному, по-видимому за счет более интенсивных вторичных течений на внутренней поверхности сферического канала напротив точки входа (рис. 3.11 в). Радиальная скорость становиться положительной во всей области за счет отражения потока от поверхности внутренней сферы, а область ее отрицательных значений смещается в область минимума дозвукового скачка.
При значениях чисел Рейнольдса более 3000 полученные гидро- и термодинамические характеристики не отвечают физическому смыслу для описания течений в данной области, поэтому необходимо применение математических моделей учитывающих возникающую турбулизацию.
Полученные результаты коррелируют с известными работами [4, 109]: критическое число Маха составляет порядка 0,4, что подтверждает, прежде всего, физическую адекватность применения полученной математической модели в указанном диапазоне чисел Маха. При этом учет сжимаемости особенно актуален в области чисел Рейнольдса от 600 и Маха от 0,08 до 0,4.
Это означает, что строгой границы по числу Маха отделяющего сжимаемое от несжимаемого течения вязкого газа отсутствует, и в каждом конкретном случае полученные результаты позволяют оценить ошибку применения модели для несжимаемой среды в зависимости от геометрических, гидродинамических и теплофизических масштабов.
Предлагаемая математическая модель течения вязкой сжимаемой среды в сферическом канале, позволила получить не только графики зависимостей для примитивных переменных, но и проанализировать влияние дифференциальных компонентов входящих в уравнения Навье-Стокса и сохранения энергии (рис. 3.12).
Градиенты давления представлены на рис. 3.12 (а — по угловой координате; б-в радиальном направлении), анализ показывает, что они имеют противоположное направление, и величина градиента по угловой координате на порядок больше, чем по радиальной.
Конвективные составляющие уравнений движения по координатам представлены нарис. 3.12 в и д, а вязкостные слагаемые соответственно г и е. Их рассмотрение позволяет сделать вывод, что основной перенос импульса течения происходит за счет конвективной составляющей, ее величина на 4 порядка больше вязкостной по абсолютному значению. Из рис. 3.12 (ж, з, и) следует, что во входном участке увеличение внутренней энергии происходит в основном за счет вязкого трения (диссипации механической энергии).
Описание пакета предметно-ориентированного прикладного программного обеспечения
На основании предложенного алгоритма разработано прикладное программное обеспечение, позволяющее численно определять гидро- и термодинамические характеристики потока вязкой сжимаемой среды в сферическом канале.
Программирование осуществлялось с использованием объектно-ориентированного языка Borland Delphi 7. Конечный программный продукт предназначается для расчетов на ПК под управлением ОС Windows. При разработке программного обеспечения использовано разделение расчета на потоки, что обеспечивает стабильность работы всей системы.
Расчет осуществляется по исходным данным, задаваемым в главном окне программы. Основными исходными данными для расчета являются количество узлов разбиения расчетной области, величина приращения по маршевой координате (времени), радиус сферы, высота ее подъема и расход подаваемой среды, причем последние три параметра задаются в размерном виде, а программа переводит их к безразмерному виду в процессе вычислений.
Контроль устойчивости расчетного процесса выполняется автоматически. Необходимые значения параметров устойчивости: сеточные числа Рейнольдса, число Куранта и максимально возможное приращение шага по времени проверяются на странице «Условия устойчивости», (общий вид которой представлен на рис. 4.4).
Результаты расчета, начиная с первого цикла, можно контролировать в режиме реального времени. Полученные данные для радиальной и меридиональной скорости, давления и температуры могут быть представлены в виде графиков (рис. 4.5), таблиц (рис. 4.6) или диаграмм.
Программное обеспечение предусматривает запись данных в файл, причем сохранение обеспечивается через заданное пользователем количество циклов расчета, разработанного формата ( .2Dg). Это позволяет просматривать развитие картины течения от начального состояния до стационарного режима. Просмотр осуществляется при помощи клавиш просмотра размещенных на панели инструментов окна графиков (рис. 4.5). Настройка «умолчаний» программы (рис. 4.9) предусматривает задание атмосферного и питающего давления, диаметра подводящего патрубка и условия окончания расчета, в качестве которых может использоваться достижение конкретного значения времени реального течения процесса или достижение стационарного режима. і Параметры Общие Сохранение Важные умолчания Методы решения Граничные условия ,,. ." . :,5. Па Па Важные умолчания Питающее давление; 125600 Атмосферное давление: 1013 Диаметр подводящего патрубка: 0,008 м Угол в = л/ р d Время остановки: 1 jjc Г" Стационарный режим через І80 -Ч циклов »та., ок Отмена Рис. 4.9 - Окно настройки параметров по умолчанию В программе предусмотрены следующие методы решения: - приближение тонкого слоя; - расчет в консервативных переменных; - метод Мак-Кормака. Выбор метода расчета выполняется установкой соответствующего включателя (рис. 4.10). Граничные условия и порядок точности аппроксимации граничных условий определяется в окне представленном на рис. 4.11.
Для обеспечения четвертого порядка точности при аппроксимации использована пятиточечная схема.
Задание плавно возрастающего источника (по экспоненциальному закону) позволяет избежать возникновения разрывов в решении при расчете в начальный период.
Наличие большого количества циклов для выхода на стационарный режим течения предопределяет системные требования к аппаратному обеспечению ПК. Для обеспечения «быстрого» выхода на стационарный режим (в течении 3-4 часов при сеточной области 21x21) необходимо использовать процессор с тактовой частотой не мене 1 ГГц.
