CC BY
6Проектирование, сооружение и эксплуатация систем трубопроводного транспорта
УДК 519.63+533.6
МЕТОД РАСПАРАЛЛЕЛИВАНИЯ АЛГОРИТМА ЧИСЛЕННОГО РЕШЕНИЯ ПОЛНОЙ СИСТЕМЫ УРАВНЕНИЙ НАВЬЕ — СТОКСА
METHOD OF PARALLELIZATION OF NUMERICAL ALGORITHM OF NAVIER — STOKES EQUATIONS COMPLETE SYSTEM
Р. Е. Волков, А. Г. Обухов
R. E. Volkov, A. G. Obukhov
Тюменский государственный университет, Тюменский государственный нефтегазовый университет, г. Тюмень
Ключевые слова: система уравнений газовой динамики; полная система уравнений Навье — Стокса; краевые условия Key words: system of gas dynamics equations; complete system of Navier — Stokes equations;
boundary conditions
Теоретические и численные исследования, проведенные в работах [1-8], подтвердили общую схему возникновения и последующего функционирования восходящего закрученного потока. В указанных работах были изучены течения воздуха в разных частях восходящего закрученного потока с использованием модели, основанной на решении системы уравнений газовой динамики. Основная идея предложенной в [9] схемы возникновения восходящего закрученного потока заключается в том, что в результате локального прогрева поверхности суши или водной поверхности появляется восходящий поток воздуха. Замещающее его радиальное течение под действием силы инерции Кориолиса приобретает осевую закрутку.
В работах [10-16] предприняты попытки исследований сложных течений, предполагающих математическое моделирование и численные расчеты трехмерных неста-
92
Нефть и газ
% 2, 2016
ционарных течений сжимаемого вязкого теплопроводного газа в целом. Для описания сложных нестационарных трехмерных течений вязкого, сжимаемого, теплопроводного газа в указанных работах используется модель сжимаемой сплошной среды, основанная на численном решении полной системы уравнений Навье — Стокса. Полная система уравнений Навье — Стокса в дифференциальной форме передает законы сохранения массы, импульса и энергии. Поэтому эта модель наиболее адекватно описывает физические процессы конвективных течений газа в замкнутых объемах при локальном нагреве подстилающей поверхности в условиях действия силы тяжести [14-16], а также в восходящих закрученных потоках газа при холодном продуве [10-12] и локальном нагреве под действием силы тяжести и Кориолиса [13].
Необходимость использования полной системы уравнений Навье — Стокса, а не системы уравнений газовой динамики при математическом моделировании подобных течений обусловлена необходимостью корректно учитывать диссипативные свойства вязкости и теплопроводности газа, как движущейся сплошной среды.
Учет диссипативных свойств газа, с одной стороны, наиболее адекватно отражает физические процессы, происходящие в сложных течениях газа. С другой стороны, наличие диссипативной функции в полной системе уравнений Навье — Стокса позволяет сгладить нефизические осцилляции типа «пилы», неизбежно возникающие при численном решении системы дифференциальных уравнений [17].
Ясно, что наличие диссипативных членов в полной системе уравнений Навье — Стокса и учет действия силы Кориолиса значительно усложняет численное моделирование течений газа, особенно в трехмерном и нестационарном случае. Усложнение численного моделирования в подобной ситуации имеет как математический, так и вычислительный аспекты. С точки зрения математики полная система уравнений Навье — Стокса представляет собой систему пяти сложных и громоздких дифференциальных уравнений с частными производными. Наличие большого числа частных производных в полной системе уравнений Навье — Стокса приводит и к значительным вычислительным проблемам, связанным с разностной аппроксимацией производных и с выбором величины шагов расчетной сетки, как по пространству, так и по времени. Кроме того, большое число значительных по объему трехмерных промежуточных и итоговых числовых массивов и необходимость выбора малых по величине расчетных шагов по времени приводит к тому, что расчеты требуют больших объемов памяти и очень большого общего времени счета.
Трудностей, связанных с большим временем счета, можно было бы избежать, при-мененив современные разработанные сложные схемы расчета (неявные разностные схемы, схемы «предиктор — корректор» и т. д.). Однако применение такого рода вычислительных схем в данных расчетах авторами работ признано несвоевременным. Поскольку численные расчеты нестационарных трехмерных течений вязкого сжимаемого теплопроводного газа на основе решения полной системы уравнений Навье — Стокса в указанных выше работах проводятся впервые, то для этого желательно использовать простые и отработанные рабочие вычислительные схемы и алгоритмы. Исходя именно из этих соображений, за основу была выбрана явная разностная схема расчета.
Другим возможным способом сокращения времени расчета может быть изменение алгоритма численного решения полной системы уравнений Навье — Стокса. В частности, попытка распараллелить вычислительную процедуру для более эффективного использования ресурсов вычислительной системы. В данной работе описывается предпринятая попытка распараллеливания вычислительного алгоритма для численного решения полной системы уравнений Навье — Стокса. Описываются и сравниваются результаты расчета одного из видов течений газа в восходящем закрученном потоке с использованием программ расчета с последовательной и параллельной структурой.
Полная система уравнений Навье — Стокса. Начальные и граничные условия. Численная реализация.
Для описания сложных течений сжимаемой сплошной среды, обладающей диссипа-тивными свойствами вязкости и теплопроводности, используется полная система урав-
% 2, 2016
Неф ть и газ
93
нений Навье — Стокса, которая, будучи записанной в безразмерных переменных с учетом действия силы тяжести и Кориолиса в векторной форме, имеет следующий вид [18]:
(1)
где Щ = 0, 001 , и Ко 1, — постоянные значения безразмерных коэффици-
458333/0
ентов вязкости и теплопроводности; / — время; х, у, г — декартовы координаты; 7 — плотность газа; V = (и, V, г) — вектор скорости газа с проекциями на соответствующие декартовы оси; Т — температура газа; g = (0,0, g) — вектор ускорения
силы тяжести; w = 1, 4
= (av bw,
au,bu)
a = 2&sin.'.срф , b = 2& cos /.срф , = ;
показатель политропы для воздуха; вектор ускорения силы Кориолиса,
^^ — вектор угловой скорости вращения Земли;
.'.(рф — широта точки О (начала декартовой системы координат Охуг , вращающейся вместе с Землей).
В качестве начальных условий при описании конкретного течения сжимаемого вязкого теплопроводного газа, в случае постоянных значений коэффициентов вязкости и теплопроводности, взяты функции, задающие точное решение [19] системы (1):
к =
lxi
00
T0 (Z1 = , kz к =-
u = 0, v = 0, w = 0, l = 0, 0065 К/м, x = 50 м, T = 288° K ,
T
00
™0 (Z(1= )kz
=
00
w g
00
= c°0nst >
(2) (3)
к
Расчетная область представляет собой прямоугольный параллелепипед с длинами сторон х0 = 1 , у0 = 1 и г0 = 0,04 вдоль осей Ох , Оу и Ог соответственно. Для
/л 0 л 0
плотности на всех шести гранях параллелепипеда ( х = 0 , х = х , у = 0 , у = у ,
г = 0 , 2 = ) ставится «условие непрерывности» потока [20]. Краевые условия для компонент вектора скорости газа соответствуют «условиям непротекания» для нормальной составляющей вектора скорости и «условиям симметрии» для двух других компонент вектора скорости течения [20]. Для температуры на всех шести гранях задаются условия теплоизоляции [20]. Продув газа через вертикальную трубу моделируется заданием вертикальной скорости течения газа в зависимости от времени Ь в виде
w()0,=003 1 Iexp 1() t
(4)
через квадратное отверстие размером 0,1^ 0,1 в центре верхней грани расчетной об- ласти. Расчеты проводились при следующих входных параметрах: масштабные раз- мерные значения плотности, скорости, расстояния и времени равны соответственно
7Т00 = 1, 2928 кг/м3, U00 = 333 м/с, Х00 = 50 м, t00 = Х00 / U00 = 0,15c .
ТУТ
и
1
94
Неф ть и газ
% 2, 2016
Разностные шаги по трем пространственным переменным — фАх = фАу = 0,005 (раз- мерное значение 0,25 м), фАг = 0, 004 (размерное значение 0,2 м), а шаг по времени —
фА/ = 0, 001 (размерное значение 0,00015 с).
Сопоставление последовательного и параллельного вычислительного алгоритма численного решения полной системы уравнений Навье — Стокса.
Расчетная область заполняется трехмерной сеткой узлов пересечения трех семейств
плоскостей у = Уj 2 = 2к уу = У фАу
х □= хг , , , =у т ,
XI □= I
фАх = I I
гк = к фАг = к к , 0 У г У Ь , 0 У у У М , 0 У к У N .
Разностные шаги по трем пространственным переменным
I = х0 / Ь , т = у0 / М, к = г0 / N . Фрагмент трехмерной расчетной сетки узлов изображен на рисунке.
слой к+1
/
/ . слой к
/ (г,],к) (¡+1,],к)
г=кк Ы-Ш/*
1 >у=?т / / # слой к-1
/ / '(г,},к-1) (¡+1,},к-1)
х=г1 /
Рисунок. Фрагмент трехмерной расчетной сетки узлов
При последовательной организации вычислительного алгоритма по известным в начальный момент времени искомым функциям с помощью явной разностной схемы вычисляются значения искомых функций во внутренних точках прямоугольного параллелепипеда. Данная последовательная вычислительная процедура для первого уравнения системы (1), записанного в скалярном виде
7Г+ + и Т + V Т,, + Ж Т + Т I Иг + V,, +
) = »
представляется следующим образом. Воспользуемся для аппроксимации производной по времени значениями функции плотности с двух последовательных временных слоев в узле (г у к)
аТ
Т
п+1
а/
фАt
а для аппроксимации производных по пространственным переменным центральными-разностями значений функций с предыдущего временного слоя (см. рисунок)
ТГуг 7п 7тР ТП
7 Т
а Т г +1У,к г 1,у ,к , а Т г, _у+1,к г,у 1,к а Т г, _У,к+1 г, _У,к 1
ах ' ау ' аг 2фАг
2фАх 2фАу
% 2, 2016
Неф ть и газ
95
n n
u u
au i+1.j.k i 1.j.k
ax
2фХх
n n
v v
av ij+1k i j 1k
2фХу
n n
w w aw i, j.k+1_ i, j.k 1
«у «2
Получим разностное уравнение для вычисления плотности газа во внутреннем узле (I I к) расчетной области
n+1 n j n
Ti,jk = Ti,jk фх
iffj. k
тП ТС. ,
i 1.j.k
2фХх
ui,j,k
n n
+ Vn Ti,j+1,k i, j,k+1
- EJT - EJT n
n i, j,k 1
i, j 1,k_
+ W ijk 2ф
+
i, j,k
+ Ti, jk
Xy
ui+1.j.k ui 1.j.k
1ф~ Xx
vi. j+1.k vi. j 1.k n n
2фХу
2фXz
wi. j.k+1 wi. j.k 1
2фХг
IX
Расчет плотности и остальных четырех газодинамических функций n + 1 временного слоя во всех внутренних узлах расчетной области ведется в тройном вложенном цикле по трем пространственным переменным. После этого, используя граничные условия, значения искомых функций рассчитываются во всех точках граней, ребер и вершин расчетной области.
Для реализации программы расчета полной системы уравнений Навье — Стокса на каждом шаге по времени в параллельном режиме была использована библиотека Task Parallel Library (TPL) на платформе .NET Framework 4.0 с использованием языка программирования C# [21. 22]. Использование указанной библиотеки считается предпочтительным способом работы с потоками в среде .NET. поскольку она динамически масштабирует степень параллелизма для наиболее эффективного использования всех доступных процессоров. Для расчета значений газодинамических функций во внутренних точках расчетной области используется механизм распараллеливания вычислений. Механизм применяется к измерению z . Максимальное количество возможных создаваемых потоков равняется N .
Для каждой конфигурации компьютера определяется количество процессоров, и создается соответствующее количество потоков. Расчет распределяется между доступными потоками. После того как потоки закончили расчет одной части расчетной области, они переходят к следующей, и так, пока не будут обработаны все внутренние точки на текущем шаге по времени. В процессе расчета значений функций во внутренних точках n -го шага по времени используются данные в соответствующих внутренних точках с предыдущего n 1 -го шага по времени, а расчет граничных условий на конкретном шаге по времени требует данных с текущего шага по времени.
Массивы для хранения данных создаются таким образом, что для Z каждый элемент является отдельным объектом, и каждый поток расчета работает в рамках одного объекта, поэтому не возникает случаев взаимоблокировки одного объекта в процессе расчета.
После завершения расчета всех граничных условий происходит смена массивов с «текущих» на «предыдущие». Далее через указанные интервалы времени происходит сохранение массивов на жесткий диск. Процесс выполняется отдельно от расчета обособленным потоком. чтобы потоки расчета не ожидали завершения записи файлов на диск.
Сравнение результатов расчета последовательной и параллельной программами.
Для проверки корректности работы распараллеленной программы были проведены расчеты течений газа в восходящем закрученном потоке, инициированном вертикальным продувом газа вверх [13] с помощью последовательной и с помощью параллельной программ при одинаковых входных данных. Что касается визуального сопоставления результатов расчетов газодинамических параметров последовательной и параллельной программами, то различия между ними не фиксируются. Сопоставление характерных безразмерных численных значений газодинамических характеристик возникающего течения на высоте 1 м в фиксированный момент времени, соответствующий 30 000 расчетному шагу, приводится в таблице.
+
n
n
n
n
n
+
+
96
Неф ть и газ
% 2, 2016
Результаты расчета последовательным и параллельным алгоритмами
Газодинамический параметр Последовательный алгоритм Параллельный алгоритм Отличие, %
^гпт 0,999641 0,999728 0,009
Т ■ 'шт 0,999408 0,999427 0,002
ршт 0,999049 0,999183 0,013
ишах 0,003078 0,003012 2,14
мшт -0,004014 -0,003868 3,63
^шах 0,003447 0,003316 3,80
^шт -0,002849 -0,002924 2,63
^шах 0,001544 0,001487 3,69
Ж 5,4274 10 8 5,4225 10 8 0,09
Ж " ху 5,3722 10 8 5,3674 10 8 0,09
Ж 4,8005 10 8 4,7985 10 8 0,04
Из сопоставления результатов расчета последовательным и параллельным алгоритмами видно, что значения основных термодинамических параметров отличаются в среднем на 0,008 %. Значения скоростных характеристик отличаются в среднем на 3,18 %. Отличие значений энергетических характеристик составляет в среднем 0,07 %. Таким образом, можно сделать вывод о корректности расчетов на основе параллельной программы, и она может быть использована при исследовании сложных течений газа.
Основным достоинством параллельной программы является значительное увеличение производительности вычислений. Так, время расчета 10 000 шагов последовательной программой на компьютере с 4 ядрами и оперативной памятью 4 Гб составило 1 025 минут. Тот же самый расчет с помощью параллельной программы составил 34 минуты. В системе с 8 ядрами и 16 Гб оперативной памяти время работы последовательной программы составляет 900 минут, а параллельной программы — 20 минут.
Выводы
Предложенный вариант распараллеливания вычислительного алгоритма по вертикальному направлению обеспечил выигрыш в производительности в 30-45 раз при практически неизменных результатах расчета значений газодинамических параметров. Столь значительное сокращение времени счета позволяет увеличить число узлов расчетной сетки по всем трем пространственным переменным и существенно детализировать характер течений газа в наиболее интересных зонах, в частности в вертикальной части восходящего закрученного потока.
Исследования поддержаны Министерством образования и науки РФ (проект № 3023).
Список литературы
1. Баутин С. П, Обухов А. Г. Математическое моделирование разрушительных атмосферных вихрей. -Новосибирск: Наука, 2012. - 152 с.
2. Баутин С. П., Обухов А. Г. Моделирование спиральных течений в придонной части восходящего закрученного потока. - Екатеринбург: УрГУПС, 2011. - 80 с.
3. Баутин С. П., Белова Е. Д., Замыслов В. Е., Крутова И. Ю., Мезенцев А. В., Обухов А. Г., Баутин П. С. Математическое моделирование природных восходящих закрученных потоков типа торнадо: труды X Всероссийского съезда по фундаментальным проблемам теоретической и прикладной механики // Вестник Нижегородского университета им. Н. И. Лобачевского. - Н. Новгород: Изд-во ННГУ им. Н. И. Лобачевского, 2011. - № 4. - Ч. 2. - 594 с.
4. Баутин С. П., Обухов А. Г. Математическое моделирование и численный расчет течений в придонной части тропического циклона // Вестник Тюменского государственного университета. Физико-математические науки. Информатика - 2012. - № 4. - С. 175-183.
5. Обухов А. Г. Математическое моделирование и численные расчеты течений в придонной части торнадо // Вестник Тюменского государственного университета. Физико-математические науки. Информатика - 2012. - № 4. -С. 183-189.
%% 2, 2016
Неф ть и газ
97
6. Баутин С. П., Обухов А. Г. Математическое моделирование придонной части восходящего закрученного потока // Теплофизика высоких температур. - 2013. - Т. 51. - № 4. - С. 567-570.
7. Баутин С. П., Крутова И. Ю., Обухов А. Г., Баутин К. В. Разрушительные атмосферные вихри: теоремы, расчеты, эксперименты. - Новосибирск: Наука; Екатеринбург: Изд-во УрГУПС, 2013. - 215 с.
8. Bautin S. P., Obukhov A. G. Mathematical Simulation of the Near-Bottom Section of an Ascending Twisting Flow // High Temperature. - 2013. - V. 51. - No. 4. - P. 509-512.
9. Баутин С. П. Торнадо и сила Кориолиса. - Новосибирск: Наука, 2008. - 96 с.
10. Абдубакова Л. В., Обухов А. Г. Численный расчет скоростных характеристик трехмерного восходящего закрученного потока газа // Известия высших учебных заведений. Нефть и газ. - 2014. - № 3. - С. 88-94.
11. Обухов А. Г., Абдубакова Л. В. Численный расчет термодинамических характеристик трехмерного восходящего закрученного потока газа // Вестник Тюменского государственного университета. Физико-математические науки. Информатика. - 2014. - № 7. - С. 157-165.
12. Абдубакова Л. В., Обухов А. Г. Численный расчет термодинамических параметров закрученного потока газа, инициированного холодным вертикальным продувом // Известия вузов. Нефть и газ. - 2014. - № 5. - С. 57-62.
13. Обухов А. Г., Баранникова Д. Д. Особенности течения газа в начальной стадии формирования теплового восходящего закрученного потока // Известия высших учебных заведений. Нефть и газ. - 2014. - № 6. - С. 65-70.
14. Обухов А. Г., Сорокина Е. М. Математическое моделирование и численный расчет трехмерного конвективного течения газа // Известия высших учебных заведений. Нефть и газ. - 2013. - № 6. - С.57-63.
15. Сорокина Е. М., Обухов А. Г. Численное исследование температурной зависимости скоростных характеристик нестационарного конвективного течения газа // Вестник Тюменского государственного университета. Физико-математические науки. Информатика. - 2014. - № 7. - С. 147-156.
16. Сорокина Е. М., Обухов А. Г. Численный расчет скоростей конвективного течения газа при кольцеобразной схеме нагрева // Известия высших учебных заведений. Нефть и газ. - 2015. - № 3 - С. 84-90.
17. Баутин С. П. Представление решений системы уравнений Навье — Стокса в окрестности контактной характеристики // Прикладная математика и механика. - 1987. - Т. 51. - Вып. 4. - С. 574-584.
18. Баутин С. П. Характеристическая задача Коши и ее приложения в газовой динамике. - Новосибирск: Наука, 2009. - 368 с.
19. Баутин С. П., Обухов А. Г. Одно точное стационарное решение системы уравнений газовой динамики // Известия высших учебных заведений. Нефть и газ. - 2013. - № 4. - С. 81-86.
20. Баутин С. П., Обухов А. Г. Об одном виде краевых условий при расчете трехмерных нестационарных течений сжимаемого вязкого теплопроводного газа // Известия высших учебных заведений. Нефть и газ. - 2013. - № 5. - С.55-63.
21. Библиотека параллельных задач (TPL) [Электронный ресурс]. - Режим доступа: https://msdn.microsoft.com/ru-ru/library/dd460717
22. Параллельное программирование с помощью языка C#. [Электронный ресурс]. - Режим доступа: http://www.microsoftvirtualacademy.com/training-courses/parallel-programming-c-sharp-rus
Сведения об авторе
Волков Роман Евстафъевич, аспирант кафедры «Алгебра и математическая логика», Тюменский государственный университет, г. Тюмень, тел. + 79129211245, е-тай: ета11@гот anvolkov.ru
Обухов Александр Геннадъевич, д. ф.-м. н.,
профессор кафедры «Бизнес-информатика и математика», Тюменский государственный нефтегазовый университет, г. Тюмень тел. 89220014998, е-тай: аоЪ1иккоу@1зо&и.ги
Information about the author Volkov R. E., postgraduate oof the chair «Algebra and mathematic logics», Tyumen State University, phone: +79129211245, e-mail: email@rom anvol-kov.ru
Obukhov A. G., Doctor oof Physics and Mathematics, professor o;f the chair «Business-informatics and mathematics», Tyumen State University, phone: 89220014998, e-mail: [email protected]
