Применение метода контрольного объема для моделирования нестационарных процессов в трубопроводных системах Текст научной статьи по специальности «Строительство и архитектура»

УДК 536.2(075)

DOI 10.18698/2309-3684-2016-4-3446

Применение метода контрольного объема для моделирования нестационарных процессов в трубопроводных системах

© В.Ю. Волков1, Л. А. Голибродо1, И.Г. Зорина2, О.В. Кудрявцев1, А. А. Крутиков1, А.П. Скибин1

1АО ОКБ «ГИДРОПРЕСС», Подольск, 142103, Россия, 2МГТУ им. Н.Э. Баумана, 105005, Россия

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

Ключевые слова: метод контрольного объема, нестационарные процессы, глобальный градиентный алгоритм, метод расчета.

Введение. В общем виде полная математическая модель движения жидкости состоит из системы уравнений, описывающей нестационарное турбулентное движение вязкой жидкости [1, 2]: уравнения движения

где р — плотность; и, Uj — проекция вектора скорости и на соответствующие оси; I — время; х^, х^ — оси декартовой системы координат (/, j = 1, 2, 3); Ту — тензоры напряжений; р — статическое

давление; S — источник объемных сил.

Однако для промышленности для расчета гидравлических трубопроводных систем требуется быстро получить интегральные результаты (расходов и гидравлических напоров), поэтому математическая модель сводится к решению двух законов Кирхгофа: сохранения мас-

уравнение неразрывности

сового баланса и потерь давления между потребителями. Таким образом, задача по определению интегральных характеристик для потерь давления (напора) между узлами i и j и описывается выражениями [2, 3]:

Ht - Hj = hj = Q + mQ22; - Di = 0,

i

где H — напор в узле, м; h — потери напора, м; i, j — индексы узлов гидравлической связи; r — коэффициент сопротивления; Q — объемный расход, м3/с; n — показатель степени, зависящий от режима течения; m — коэффициент местного сопротивления; Di —

расход в i-м узле (если величина положительная, то он определяет сток массы), м3/с.

Стационарные задачи потокораспределения сетевой гидравлики достаточно хорошо изучены. Так, первые достижения по решению контурных задач трубопроводных систем относятся еще к трудам Харди Кросса [4]. Наиболее универсальным методом решения задачи потокораспределения является глобальный градиентный алгоритм (Global Gradient Algorithm, GGA), предложенный Тодини и Пилати в 1987 г. [5]. Данный метод наиболее прост и эффективен с точки зрения реализации и вычислительной производительности [6-9]. На этом алгоритме основано наиболее известное отраслевое программное обеспечение (ПО) для гидравлических расчетов EPANET [10]. В 2009 г. GGA был расширен Густолизи возможностью эквивалентирования, что позволило сократить временные затраты на расчет [11]. Однако все упомянутые алгоритмы априори являются стационарными.

Во многих случаях для анализа изменения работы системы во время ее эксплуатации требуется решение нестационарных задач, связанных с наполнением и опорожнением гидравлических емкостей. Однако, как было показано в 2011 г. в статье [12], при решении с помощью стандартных программ и методов можно получить осцилляции давления и расходов даже для очень упрощенной сети, состоящей из трех узлов и трех связей (рис. 1). Пример такого расчета с помощью стандартного гидравлического ПО EPANET для шага по времени 5 мин представлен на рис. 2, на котором точками показаны значения в узлах и связях, линиями — численное решение.

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

Н, м Q, кг/с

Рис. 2. Результаты расчета напора H и массового расхода Q с помощью ПО EPANET (—, —, —) и точное решение (Exact, , , )

В настоящей статье представлено расширение разработанного ранее метода контрольного объема (МКО) для расчета стационарных процессов потокораспределения в гидравлических сетях [13]. Данное расширение МКО разработано для медленно протекающих нестационарных процессов в гидравлических сетях аналогично методам, изложенным в статьях [12-15], и не предназначено для расчета быстро протекающих местных процессов, таких как гидроудар [16].

Метод контрольного объема для расчета стационарных процессов. Для расчета гидравлической сети применен метод контрольного объема. Построение гидравлической сети осуществляется таким образом, чтобы центру каждого контрольного объема (КО) соответствовал узел гидравлической сети, а грани контрольных объемов — серединам гидравлических связей (трубопроводов), как представлено на рис. 3.

При построении дискретного аналога проводится линеаризация уравнения сохранения движения [13]. При записи дискретного аналога совершается переход от конвективных потоков Г к среднеинте-гральной скорости в связи.

Ж

АХ

Ы2

8х

Р Ж

Ре

Р) 1 5« '|----'е ' •- Е

Рис. 3. Контрольные объемы для уравнения движения (а) и уравнения неразрывности (б)

Таким образом, дискретные аналоги уравнения движения можно представить

2 АХ дР

и* = —

Р и

\£ дх

= —й/Р Р

РЖ

АХ

(1)

где — коэффициент дискретного аналога уравнения движения, вычисляемый по формуле

2Ах

=' I Iе'

Р Ы £

(2)

Здесь АХ — длина гидравлической связи, м; р — плотность движущейся среды, кг/м3; и* — скорость в гидравлической связи, м/с; £ — коэффициент гидравлического сопротивления; Р — давление, Па; Р, Ж, Е, Б — индексы, относящиеся к центру КО; е, s — индексы, относящиеся к грани КО.

Дискретный аналог уравнения неразрывности запишем в виде

арРр = о^Рж ^ ®еРЕ + Ь,

(3)

где коэффициенты дискретного аналога уравнения неразрывности определяются по формулам

Ое = МА, =рфЛ, =рф&, Ь = £ ар = аж + аЕ ^ (4)

5хе

5х

Б

5х„

Здесь Б — площадь поперечного сечения гидравлической связи, м ; Ь — источниковый член, кг/с.

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

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

Шаг по времени задается искусственно и определяется как

* *

А/ = / - / , где /, / — индексы, относящиеся к текущему и предыдущему моментам времени.

В качестве гидравлических емкостей рассмотрим баки постоянного сечения. В этом случае объем жидкости в баке V = БИ. Тогда изменение объемного расхода q в узлах за интервал At описывается

как Atq = Б (И (/) - И (/*)).

Таким образом, изменение объемного расхода

q = Б(И(/) - И(/*)) q = А/ '

Переходя к массовым расходам Q, имеем Q = qp. Перепишем выражение, заменив напор И на давление Р:

= Б (Р(/) - Р(/*)) Е А/ "

Следовательно, источниковый член дискретного аналога уравнения неразрывности имеет вид

Q = ^Р-А-Р* Е А/ Е А/

или при стандартной записи дискретного аналога

Q = аРР - аРР*.

Для решения системы уравнений (1)-(4) применена следующая модифицированная итерационная процедура:

1) вводятся предполагаемые поля скорости и давления, задается шаг по времени;

2) рассчитываются значения с1 для всех гидравлических связей по формулам (2);

3) определяются коэффициенты дискретного аналога поля давления (4) и поля давления и градиентов давления, а также массовые потоки через грани КО;

4) рассчитываются значения скорости для каждой гидравлической связи с использованием коэффициентов дискретного аналога уравнения движения и полученного градиента поля давления;

5) возврат к п. 2 до тех пор, пока не будет достигнута сходимость, ее критерием является достижение невязки узловых массовых балансов, значения которого меньше заданного;

6) переход на следующий шаг по времени, пока не закончится расчет.

Для разработанного метода контрольного объема не требуется выделение гидравлических контуров, он сводится к решению единого поля давления сразу для всей расчетной области. При этом, как было показано в работе [13], этот метод позволяет решать задачи по-токораспределения размерностью более 10 млн узлов и связей. Разработанный метод стабильнее, чем EPANET при решении задач большой размерности (более 2000 узлов) [13].

Тестирование метода КО было проведено на ряде задач с различными шагами по времени. В качестве наиболее характерных рассмотрены задачи расчета гидравлических сетей из статей [12, 15].

Простая система. Система представляет собой два соединенных между собой бака постоянного сечения, наполненных водой (см. рис. 1). Начальный уровень воды (напор) в баках 1 и 2 равен 20 и 30 метров соответственно. Баки посредством двух трубопроводов (трубы 2 и 3 на рис. 2) соединены с резервуаром 3, в котором поддерживается постоянное нулевое давление. Диаметр баков равен 3,56 м. Перепад давления на участках трубопроводов рассчитывается с помощью уравнения Хазена — Вильямса [6]. Для всех трубопроводов коэффициент шероховатости Хазена — Вильямса принимается равным 130. Длина всех трубопроводов составляет 100 м, диаметр трубы 1 равен 200 мм, а труб 2 и 3 — 100 мм. В начальный момент времени расход по трубопроводам отсутствует. В нулевой момент времени вода из баков начинает перетекать по трубопроводам до тех пор, пока баки полностью не опустеют.

На рис. 4 представлены результаты расчета данной системы на различных шагах по времени (5, 15, 60 и 120 мин соответственно). Сравнение результатов, полученных с помощью метода контрольного объема, проводилось с данными статьи [15], где для расчета подобных систем применен обобщенный глобальный градиентный алгоритм (Generalized Global Gradient Algorithm, G-GGA).

Полученные результаты согласуются с результатами статьи [15]. Погрешность не превышает 0,35 % относительно G-GGA. При этом максимальные абсолютные отклонения расходов не превышают 0,5 кг/с, а напоров — 0,1 м.

Погрешность расчета определялась по формуле:

где Е — значение переменной в 1-й момент времени; индексы МКО, О-ООЛ — значения, полученные с помощью МКО и данных статьи [15].

Н, м

-40 -80 -120

0 0,5 1,0 1,5 2,0 2,5 3,0 3,5 (,ч 0 0,5 1,0 1,5 2,0 2,5 3,0 3,5 (,ч

а

Н, м

0 0,5 1,0 1,5 2,0 2,5 3,0 3,5 Н, м

б, КГ/с

40

20 -

0

-20 {

-40 /

60

0 0,5 1,0 1,5 2,0 2,5 3,0 3,5

—*-*-*-* -20

0 1 2 3 4 5 67 ?, ч

01 2 3 4 5 6 7 ч

Н, м

Q, кг/с

0 1 234 567 1,ч

01 2 3 4 5 6 1 1,ч

Рис. 4. Результаты расчета с использованием различных шагов по времени:

, , — МКО; , , — в-ввЛ (в скобках указаны номера узлов/связей):

а — Д/ = 5 мин; б — Д/ = 15 мин; в — Д/ = 60 мин; г — Д/ = 120 мин

Полученные результаты расчета также согласуются с данными, представленными в статье [12]. При этом в гидравлической сети не возникает осциляций давления для всех шагов по времени, в отличие от решения с помощью ЕРАКЕТ.

Расчет сложной системы. Для проверки численной устойчивости МКО была рассмотрена гидравлическая сеть Апулии, представленная на рис. 5. Постановка задачи аналогична данным статьи [15]. Перепад давления на участках трубопроводов рассчитывается с помощью уравнения Дарси — Вейсбаха. Каждый из 23 узлов гидравлической сети представляется в виде бака с постоянной площадью поперечного сечения 10 м2. Для всех расчетов шаг по времени был равен 10 мин.

Были рассмотрены три вычислительных эксперимента переходных процессов в гидравлической сети:

• во всех баках в нулевой момент времени задан постоянный нулевой напор (все баки пусты и наполняются во время переходного режима);

• во всех баках в нулевой момент времени задан постоянный напор 72,8 м (все баки наполнены одинаково, и во время переходного режима происходит их опустошение);

• в начальный момент в баках 1-6, 16-19, 21 и 23 задан нулевой напор, в других — постоянный напор 72,8 м.

Для всех рассматриваемых случаев в узле 24 (резервуар) поддерживается постоянный напор, равный 36,4 м.

Таким образом, рассмотрены задачи наполнения и опустошения баков. Результаты расчета (изменение массовых расходов и узловых давлений) представлены на рис. 6. Полученные с помощью метода контрольного объема результаты согласуются с данными статьи [15]. При этом в гидравлической сети не возникает осциляций давления. Можно заключить, что разработанный метод устойчив при решении переходных процессов для задач наполнения и опорожнения гидравлических емкостей независимо от размерности задачи.

б, кг/с

8 10 (,ч

Рис. 6. Результаты гидравлического расчета сети Апулии:

а — эксперимент 1; б — эксперимент 2; в — эксперимент 3

Заключение. Представленное расширение разработанного метода контрольного объема для расчета нестационарных процессов потоко-распределения в гидравлических сетях разработано для медленно протекающих процессов в гидравлических сетях и не предназначено

для расчета быстро протекающих местных явлений, таких как гидроудар [12].

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

ЛИТЕРАТУРА

[1] Себиси Т., Брэдшоу П. Конвективный теплообмен. Физические основы и вычислительные методы. Москва, Мир, 1987, 592 с.

[2] Фомин А.А., Фомина Л.Н. О стационарном решении задачи течения несжимаемой вязкой жидкости при больших числах Рейнольдса. Математическое моделирование и численные методы, 2015, № 4 (8), с. 92-102.

[3] Меренков А.П., Хасилев С.Ю. Теория гидравлических цепей. Москва, Наука, 1985, 279 с.

[4] Cross H. Analysis of flow in networks of conduits or conductors. Engineering Experiment Station, University of Illinois, 1936, 38 p.

[5] Todini E., Pilati, S. A gradient method for the analysis of pipe networks. International Conference on Computer Applications for Water Supply and Distribution, UK, Leicester Polytechnic, 1987, 20 p.

[6] Zheng Y.W., Rong H.W., Walski T.M., Yang S.Y., .Bowdler D., Baggett C.C. Efficient pressure dependent demand model for large water distribution system analysis. 8th An. Intern. Symp. on Water Distribution System Analysis, USA, Cincinnati, Ohio, 2006, pp. 1-15.

[7] Creaco E., Franchini M. Comparison on Newton-Raphson global and loop algorithms for water distribution network resolution. Journal of Hydraulic Engineering,, 2013, pp. 313-320.

[8] Todini E. On the convergence properties of the different pipe network algorithms. 8th Annual Water Distribution Systems Analysis Symposium. USA, Cincinnati, Ohio, 2006, pp. 1-16.

[9] Todini E., Rossman L.A. Unified Framework for Deriving Simultaneous Equation Algorithms for Water Distribution Networks. Journal of Hydraulic Engineering, 2013, vol. 139, no. 5, pp. 511-526.

[10] Rossman L.A. EPANET 2, User's Manual. Water Supply and Water Resources Division National Risk Management Research Laboratory Cincinnati, OH 45268, 2000, 200 p.

[11] Giustolisi O., Laucelli D., Berardi L., Savic D.A. Computationally Efficient Modeling Method for Large Water Network Analysis. Journal of Hydraulic Engineering, 2012, vol. 134, no. 4, pp. 313-326.

[12] Todini E. Extending the Global Gradient Algorithm to Unsteady Flow Extended Period Simulations of Water Distribution Systems. Journal of Hydroinfor-matics, 2011, vol. 13, no. 3, pp. 167-180.

[13] Patankar S.V. Numerical heat transfer and fluid flow. Taylor and Francis, 1981, 196 p.

[14] Belova O., Skibin A., Volkov V. Control-volume Method for Extralarge Network Hydraulic Analysis. Journal of Hydroinformatics, 2014, no. 70, pp. 123-131.

[15] Giustolisi O., Berardi L., Laucelli D. Generalizing WDN Simulation Models to Variable Tank Levels. Journal of Hydroinformatics, 2012, no. 14, pp. 562-573.

[16] Wood D.J. Waterhammer Analysis — Essential and Easy (and Efficient). Journal of Environmental Engineering, 2005, vol. 131, no. 8, pp. 1123-1131.

Статья поступила в редакцию 21.11.2016

Ссылку на эту статью просим оформлять следующим образом: Волков В.Ю., Голибродо Л.А., Зорина И.Г., Кудрявцев О.В., Крутиков А.А., Скибин А.П. Применение метода контрольного объема для моделирования нестационарных процессов в трубопроводных системах. Математическое моделирование и численные методы, 2016, № 4 (12), с. 34-46.

Волков Василий Юрьевич окончил МГТУ им. Н.Э. Баумана. Канд. техн. наук, инженер-конструктор 1-й категории АО ОКБ «ГИДРОПРЕСС». Автор более 25 публикаций. e-mail: vasya-volkov@yandex.ru

Голибродо Лука Антонович окончил МГТУ им. Н.Э. Баумана. Инженер-конструктор 1-й категории АО ОКБ «ГИДРОПРЕСС». Автор более 20 публикаций. e-mail: luka.golibrodo@mail.com

Зорина Ирина Григорьевна окончила МГТУ им. Н.Э. Баумана. Канд. техн. наук, доцент кафедры «Высшая математика» МГТУ им. Баумана. Автор более 20 публикаций. e-mail: zorina.ig@gmail.ru

Кудрявцев Олег Вячеславович окончил МГТУ им. Н.Э. Баумана. Инженер-конструктор 1-й категории АО ОКБ «ГИДРОПРЕСС». Автор более 20 публикаций. e-mail: ovkudr@yandex.ru

Крутиков Алексей Александрович окончил МГТУ им. Н.Э. Баумана. Канд. техн. наук, инженер-конструктор 1-й категории АО ОКБ «ГИДРОПРЕСС». Автор более 35 публикаций. e-mail: krutikov_aa@mail.ru

Скибин Александр Петрович окончил МГТУ им. Н.Э. Баумана. Канд. техн. наук, начальник группы АО ОКБ «ГИДРОПРЕСС». Автор более 80 публикаций. e-mail: askibin@yandex.ru

Applying the control-volume method to extended period simulations in pipe network hydraulics

© V.Yu. Volkov1, L.A. Golibrodo1, I.G. Zorina2, O.V. Kudryavtsev1, A.A. Krutikov1, A.P. Skibin1

JSC OKB GIROPRESS, Podolsk, 142103, Russia 2Bauman Moscow State Technical University, Moscow, 105005, Russia

For modeling piping systems we made a transition from the mass balance equations, based on 1m and 2m Kirchhoff laws, to the mathematical description of a hydraulic network using the continuity equation discretization. For this purpose we applied a control-volume method. This paper introduces an extension of the developed control-volume method for extended period simulations in hydraulic networks. This extension is devel-

oped for slow time-varying conditions in the hydraulic networks and is not intended to calculate rapidly occurring local phenomena such as waterhammer. The control-volume method was successfully applied to test tasks.

Keywords: control-volume method, global gradient algorithm, extended period simulations, transient.

REFERENCES

[1] Sebisi T., Bredshou P. Konvektivnyy teploobmen. Fizicheskie osnovy i vychis-litelnye metody [Convective heat transfer. Physical fundamentals and computational methods], Moscow, Mir Publ., 1987, 592 p.

[2] Fomin A.A., Fomina L.N. Matematicheskoe modelirovanie i chislennye metody — Mathematical Modeling and Computational Methods, 2015, no. 4 (8), pp. 92-102.

[3] Merenkov A.P., Hasilev S.Yu. Teoriya gidravlicheskikh tsepey [Theory of hydraulic networks]. Moscow, Nauka Publ., 1985, 279 p.

[4] Cross H. Analysis of flow in networks of conduits or conductors. Engineering Experiment Station, University of Illinois, 1936, 38 p.

[5] Todini E., Pilati S. A gradient method for the analysis of pipe networks. International Conference on Computer Applications for Water Supply and Distribution, Leicester Polytechnic, UK, 1987, 20 p.

[6] ZhengY.W., Rong H.W., Walski T.M., Yang S.Y., .Bowdler D., Baggett C.C. Efficient pressure dependent demand model for large water distribution system analysis. 8th Annual Intern. Symp. on Water Distribution System Analysis, USA, Cincinnati, Ohio, 2006, pp. 1-15.

[7] Creaco E., Franchini M. Comparison on Newton-Raphson global and loop algorithms for water distribution network resolution. Journal of Hydraulic Engineering, 2013, pp. 313-320.

[8] Todini E. On the convergence properties of the different pipe network algorithms. 8th Annual Water Distribution Systems Analysis Symposium. USA, Cincinnati, Ohio, 2006, pp. 1-16.

[9] Todini E., Rossman L.A. Unified Framework for Deriving Simultaneous Equation Algorithms for Water Distribution Networks. Journal of Hydraulic Engineering, 2013, vol. 139, no. 5, pp. 511-526.

[10] Rossman L.A. EPANET 2, User's Manual. Water Supply and Water Resources Division National Risk Management Research Laboratory Cincinnati, OH 45268, 2000, 200 p.

[11] Giustolisi O., Laucelli D., Berardi L., Savic D.A. Journal of Hydraulic Engineering, 2012, vol. 134, no. 4, pp. 313-326.

[12] Todini E. Journal of Hydroinformatics, 2011, vol. 13, no. 3, pp. 167-180.

[13] Patankar S.V. Numerical heat transfer and fluid flow. Taylor and Francis, 1981, 196 p.

[14] Belova O., Skibin A., Volkov V. Journal of Hydroinformatics, 2014, no. 70, pp. 123-131.

[15] Giustolisi O., Berardi L., Laucelli D. Journal of Hydroinformatics, 2012, no. 14, pp. 562-573.

[16] Wood D.J. Journal of Environmental Engineering, 2005, vol. 131, no. 8, pp. 1123-1131.

Volkov V.Yu. graduated from Bauman Moscow State Technical University, Cand. Sc. (Eng.), first rank design-engineer of JSC OKB GIDROPRESS. Author of over 25 scientific publications. e-mail: vasya-volkov@yandex.ru

Golibrodo L.A. graduated from Bauman Moscow State Technical University, first rank design-engineer of JSC OKB GIDROPRESS. Author of over 20 scientific publications. e-mail: luka.golibrodo@mail.com

Zorina I.G. graduated from Bauman Moscow State Technical University, Cand. Sc. (Eng.), Assoc. Professor of Higher Mathematics Department, BMSTU. Author of over 20 scientific publications. e-mail: zorina.ig@gmail.ru

Kudryavtsev O.V. graduated from Bauman Moscow State Technical University, first rank design-engineer of JSC OKB GIDROPRESS. Author of over 20 scientific publications. e-mail: ovkudr@yandex.ru

Krutikov A.A. graduated from Bauman Moscow State Technical University, Cand. Sc. (Eng.), first rank design-engineer of JSC OKB GIDROPRESS. Author of over 35 scientific publications. e-mail: krutikov_aa@mail.ru

Skibin A.P. graduated from Bauman Moscow State Technical University, Cand. Sc. (Eng.), JSC OKB GIDROPRESS head of the group. Author of over 90 scientific publications. e-mail: askibin@yandex.ru