Численный метод расчета нестационарных режимов гидравлических систем

Министерство образования и науки

Российской федерации

Федеральное агентство по образованию

Новосибирский государственный университет

Механико-математический факультет

Кафедра математического моделирования

Дипломная работа

МИТКИНОВА Владимира Евгеньевича

Численный метод расчета нестационарных режимов гидравлических систем

Научный руководитель

Д. ф.-м. наук, профессор

А.Ф. Воеводин

Новосибирск 2010

Содержание

Содержание………………………………………………………………………..2

Введение…………………………………………………………………………...3

Постановка задачи………………………………………………………………...4

Математическая модель узла и условия сопряжения…………………………..7

Численный метод………………………………………………………………….9

Вычисление коэффициентов……………………………………………………12

Исследование свойств схемы…………………………………………………...13

Вычислительный эксперимент………………………………………………….15

Выводы…………………………………………………………………………...19

Литература……………………………………………………………………….20

Введение

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

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

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

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

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

Постановка задачи

Уравнения, описывающие изотермические движение жидкости в трубах имеют вид [1,4]:

Здесь x-пространственная координата, t-время, ρ- плотность жидкости, g-ускорение силы тяжести, α-угол между осью канала и вектором силы тяжести. Последнее слагаемое во 2ом уравнение системы учитывает потери на трение жидкости о стенки канала, D-диаметр канала, λ-коэффициент трения.

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

где - скорость звука в жидкости.

Запишем систему в векторном виде:

Решая характеристическое уравнение:

находим собственные значения и характеристические направления нашей системы:

Заметим, что характеристики являются прямыми линиями, и их уравнения имеют вид: .

Также допускается запись системы (1), (2) в характеристической форме:

Определим инварианты Римана в виде:

Учитывая (3) и (4), запишем систему уравнений гидроудара через инварианты:

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

При и необходимо задавать по одному граничному условию, например:

Для расчетов удобно систему уравнений (1), (2) представить в безразмерном виде. Пусть , , , – безразмерные переменные, где – масштаб длины, - масштаб времени, - масштаб скорости, – масштаб давления. Тогда уравнения (1), (2) перепишутся в виде:

Выберем , , так, что , и обозначим – безразмерная скорость звука.

В результате получим систему уравнений (1), (2) в безразмерной форме:

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

Математическая модель узла и условия сопряжения

Место соединения труб с дополнительным отводом (подводом) жидкости будем описывать моделью с сосредоточенными параметрами (моделью узла). При этом будем считать, что в узле:

а) происходит полное мгновенное перемешивание жидкости;

б) нет местных сопротивлений.
Пусть уравнения (1), (2) заданы на отрезках (ребрах) графа, содержащего J вершин и K отрезков; - множество номеров отрезков, примыкающих к j-ому узлу; - множество номеров отрезков, из которых жидкость поступает в данный узел; - множество номеров отрезков, которые отводят жидкость из узла. Тогда для j-го узла можно записать следующие балансовые соотношения, вытекающие из закона сохранения массы [1,2]:

-скорость в концевом (начальном) сечении m-ой трубы.

-сосредоточенный подвод (отвод) жидкости в узле.

-площадь сечения m-ой трубы.

Рисунок 1

Для «стыковки» моделей трубы и узла сформулируем условия примыкания. Будем пренебрегать местными сопротивлениями вблизи узла по сравнению с гидравлическими сопротивлениями по длине труб. Тогда условия примыкания имеют вид:

.

Где - сосредоточенный параметр в узле (узловое давление), - давление в концевом (начальном) сечении m-й трубы.

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

Численный метод

Будем аппроксимировать систему (5), (6) схемой [2]:

Выражаем из (7) и подставляем в (8), получаем разностное уравнение:

Это уравнение можно записать в виде:

где:

.

Для решения уравнения (9) используем следующий способ дискретизации задачи [6].

Введем на отрезке [0,1] сетку с узлами в общем случае с неравномерным шагом .

Решение уравнения (9) будем искать например как решение задачи (преобразование Риккати)

Тогда для прогоночных коэффициентов и получим следующие уравнения:

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

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

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

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

Наряду с уравнением (13) на отрезке рассмотрим уравнение

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

и, кроме того,

Тогда с помощью уравнения (13) при и уравнения (18) при , полагая , , , , , на каждом отрезке будем иметь систему двух уравнений:

которые могут быть приведены к следующему виду

где

Приравнивая производные слева и справа, в узлах сетки , из уравнения (22) на отрезке и уравнения (23) на отрезке , получим разностные уравнения для неизвестных

Здесь коэффициенты , , , вычисляются по формулам (24)-(29), при этом:

Для замыкания системы (30) привлекаем граничные условия.

Вычисление коэффициентов

Как было отмечено ранее граничные условия для уравнений (15)-(17), (19)-(21) являются свободными параметрами и могут быть выбраны специальным образом с учетом свойств решений этих уравнений.

Выбирая, для уравнений (15), (19) граничные условия: , , для (16), (18) , и для (17), (21) , , в результате получаем:

Для вычисления интегралов в значениях коэффициентов , строим кубический сплайн Эрмита [7] для подынтегральной функции . При этом считается, что кроме значений функции на концах отрезков задаются значения производных и , которые вычислялись при помощи формул (22), (23).

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

Также для сравнения были проведены расчеты уравнений (31), (32) с помощью формулы трапеций.

Отметим, что коэффициенты системы (30) получаются в итоге постоянными величинами.

Исследования свойств схемы

Для исследования устойчивости метода прогонки и монотонности схемы используем принцип максимума, т.е. если в системе уравнений (30) , , и , то в матрице системы уравнений (30) присутствует диагональное преобладание, матрица является монотонной и выполняется оценка [6].

Рассмотрим коэффициенты системы уравнений (30):

Отсюда следует, что схема является монотонной и метод прогонки устойчив.

Устойчивость по времени нетрудно показать с помощью спектрального признака устойчивости [5].

Перепишем уравнение (9) в следующей форме:

и решение (33) будем искать в виде: Тогда подставляя это выражение в (33) получим следующее уравнение:

и следовательно,

Очевидно, что

Отсюда следует, что схема абсолютно устойчива.

Для исследования порядка аппроксимации системы (5), (6) схемой (7), (8) разложим точное решение системы по формуле Тейлора в окрестности точки , предполагая наличие непрерывных вторых производных:

Подставляя эти значения в (7), (8) и, используя, свойства и системы (5), (6) получим невязку:

Таким образом, учтя погрешность при вычислении коэффициентов (31), (32) , система (5), (6) имеет аппроксимацию .

Вычислительный эксперимент

Рассматривалась система трубопроводов, состоящая из 3 отрезков и 4 узлов. Длины всех отрезков одинаковы и равны 1, диаметры каналов одинаковы. Расчеты производились приведенным выше методом. Результаты расчетов сравнивались с точным решением системы (5), (6):

где выбиралась:

1) .

Начальные условия:.

Граничные условия: при

при

2) .

Начальные условия: .

Граничные условия: при

при

3) .

Начальные условия:.

Граничные условия: при

при

Таблицы погрешностей

Расчеты на момент времени .

1) Для первого случая:

2) Для второго случая:

3) Для третьего случая:

Замечание

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

2) С уменьшением шага наблюдается увеличение погрешности по , это связано с зависимостью показателя диагонального преобладания от шага .

Выводы

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

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

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

Литература

1. Воеводин А.Ф., Шугрин С.М. «Методы решения одномерных эволюционных систем». Новосибирск: Наука, 1993.

2. Воеводин А.Ф., Пономарев М.Ю. «Численный метод расчета нестационарных режимов гидравлических систем».

3. Годунов С.К. «Численное решение многомерных задач газовой динамики». Москва: Наука, 1976.

4. Меренков А.П. «Математическое моделирование и оптимизация систем тепло-, водо, нефте-, газоснабжения». Новосибирск: Наука, 1992.

5. Самарский А.А., Попов Ю.П. «Разностные схемы газовой динамики». Москва: Наука, 1975.

6. Воеводин А.Ф. «Метод факторизации для линейных и квазилинейных сингулярно возмущенных краевых задач для обыкновенных дифференциальных уравнений ». // Сиб. жур. вычисл. математики / РАН. Сиб. отд-ние.- Новосибирск, 2009.-Т. 12, №1.-с.1-15.

7. Волков Е.А. «Численные методы». Москва: Наука, 1987.