Моделирование течения несжимаемой невязкой тяжёлой жидкости и на регулярной сетке в трёхмерном пространстве Текст научной статьи по специальности «Математика»

Д. П. Ливеровский, С. П. Шевырёв. Моделирование течения несжимаемой невязкой жидкости

identification of thermoelastic characteristics for inhomo-geneous bodies. Inzhenerno-fizicheskiy zhurnal [Journal of Engineering Physics], 2014, vol. 87, no. 1, pp. 217-224 (in Russian).

УДК 533.6.011

10. Tikhonov A. N., Goncharsky A. V., Stepanov V. V., Yagola A. G. Chislennye metody resheniya nekorrekt-nykh zadach [Numerical methods for solving ill-posed problems]. Moscow, Nauka, 1990, 230 p. (in Russian).

МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ НЕСЖИМАЕМОЙ НЕВЯЗКОЙ ТЯЖЁЛОЙ ЖИДКОСТИ И НА РЕГУЛЯРНОЙ СЕТКЕ В ТРЁХМЕРНОМ ПРОСТРАНСТВЕ

Д. И. Ливеровский1, С. П. Шевырёв2

1 Аспирант кафедры математического и компьютерного моделирования, Саратовский государственный университет им. Н. Г. Чернышевского, LiverovskiyDI@gmail.com

2Кандидат физико-математических наук, доцент, заведующий кафедрой прикладной информатики, Саратовский государственный университет им. Н. Г. Чернышевского, ShevSP@rambler.ru

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

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

1. СУТЬ МЕТОДА ДАВЫДОВА

Суть метода Давыдова (метода крупных частиц) состоит в том, что использует расщепление по физическим факторам и по координатам. В методе Давыдова выделяют три этапа [1]. Будем рассматривать декартову прямоугольную систему координат. Решение будет проводиться на регулярной разностной сетке, а это означает, что расчётная область разбивается на ячейки с шагом Нх вдоль оси X, с шагом Ну вдоль оси У и с шагом Н2 вдоль оси Z. Расчётная область будет разбита на К частей вдоль оси X, М частей вдоль оси У, N частей вдоль оси Z. Ячейки сетки будут иметь вид прямоугольных параллелепипедов. Каждая ячейка будет определяться тройкой индексов (г,],к), где г определяет положение ячейки вдоль оси Z, ] определяет положение ячейки вдоль оси У и к определяет положение ячейки вдоль оси X. Все расчётные параметры (давление и компоненты вектора скорости) будут отнесены к центру ячейки. Значения на границах ячейки будут иметь полуцелый индекс по соответствующей координате. Рассмотрим этапы метода Давыдова подробнее. Приведённые далее формулы получаются аналогично формулам для двумерного случая [2].Течение тяжелой несжимаемой невязкой жидкости в трёхмерном пространстве моделируется путем решения краевых задач для системы уравнений Эйлера и уравнения неразрывности:

du d (u2) d(uv) d (uw) 1 dp ^

dt dx dy dz p0 dx

dv d (v2) d (vu) d(vw) 1 dp ^

dt dy dx dz po dy

dw d (w2) d (wu) d (wv) dt dz dx dy

du dv dw о

dx dy dz

+ a. di = -

+ P0 dz 9l

© Ливеровский Д. И., Шевырёв С. П., 2014

На эйлеровом этапе опускают конвективные члены, таким образом разностные аналоги первых трёх уравнений системы (1) примут вид(2), где р? ^ < к, ^ < к, -о ^ < к, ^ < к — значения давления и компонентов вектора скорости в ячейке (г,.?, к) на д-м шаге по времени, т — величина шага по времени, а

Нх , Ну , Н^ — размеры ячейки регулярной разностной сетки. Давление берётся на границах

о I о

о рг ' 7 ' к + 1 + рг ' 7 ' к

расчётных ячеек и выражается по формуле рг < ^ < к+1 = ——-^-и т. д.

, 7 , к - иг , 7 , к + +1 Рг,7,к-2 = 0

Т Р0 Нх

-г 7 , к _ -г , 7 , к + +1,к 2,к = 0 (2)

Т р0 Ну '

~о о _ ро _ ро

,7,к _ <,7,к , 1 Рг+ 2 7,к Рг- 17,к

+--^-— = _£•

Т Ро Ну

Отсюда мы можем получить явные формулы для нахождения компонент вектора скорости в конце эйлерова этапа:

м0,7,к = м0,7,к _ 2р0Нх (р?,7,к + 1 _ Ро,7,к-1

= -г,7,к _ 2р0н (рг,7,к+1 _ рг,7,к-1(3) = _ 2р0Н (р9,7,к + 1 _ р0,7,к-1) _

Далее, рассмотрим четвёртое уравнение системы (1), это уравнение неразрывности. Введём следующие обозначения:

дм д- д-1 ^

7Г + я" + я" = О (4)

дх ду дг ди д- д1 ~

7Г + + = I' (5)

дх ду дг

где I — это дивергенция вектора скорости, а О — это дивергенция вектора скорости на эйлеровом этапе, играющая ключевую роль в расчетах несжимаемой жидкости. Дивергенция должна равняться нулю в любой момент времени на любом этапе расчета, как в дифференциальной форме, так и в разностной форме, это следует из четвёртого уравнения системы (1). Запишем разностные аналоги выражений (4), (5) через скорости на границах ячеек:

и .. , 1 и .. , 1 1 . 1 • 7 1 • 1-7

—о г,7,к +1 г,7,к-2 , г7 + 2,к г,7-2 ,к , г+ 2 7,к г-2 7,к

I г,7,к = -^Т-2 +--2-Н-^--1--^-Н-^-' (6)

Нх Ну Н^

о о о о о о

и .. , 1 и .. , 1 V . . , 1 , V . 1, 1. 1 • 7 1 • 1-7

„о г,7,к + 2 г7,к-2 , г7 + 2,к г,7-2 ,к , г+ 2,7,к г-2 7,к /"7\

О7,к = - ь-^ +--2—Н-2--1--2-^-2-' (7)

Нх Ну

„у

о + о о + о о + о

о иг,7,к + 1 + иг,7,к о -г,7 + 1,к + -г,7,к о 1г+1,7,к + 1г,7,к

где и^ ,1 = ——---—, . ^ = ——-—, 1°. ^ , =--—,

м г,7,к +1 2 г,7 + 2,к 2 г+1,7,к 2

Так как на каждом шаге по времени I-^ к = 0 для любых троек (г,^, к), то, подставив в (6) и—ок, -о 7 к и 7 к, выраженные из формул (3), получим следующее выражение:

/ 111 А т т

Ог ,7 , к + рг ,7 , к —^ + + —^ —--—-(рг 7 , к+2 + рг 7 , к-2)_

\Нх Ну Н^ 2Р0 4Нх Р0

т

(Рг,7+2,к + Рг,7-2,к) _ ——2-(Рг+2,7,к + Рг-2,7,к) = 0. (8)

.. 2— (Рг,7+2,к + рг,7-2,к) _ -т—2— (рг+2,7,к + Рг-2,7,к 4Ну Р0 4Н2 Р0

В том случае, когда все переменные относятся к одному временному слою, будем опускать верхний индекс. Из уравнения (8) уже можно выразить рг,7,к и получить матричное уравнение, но можно получить матричное уравнение с более простым шаблоном, использовав разложение в ряд Тейлора.

Д. П. Ливеровский, С. П. Шевырёв. Моделирование течения несжимаемой невязкой жидкости

Записав разложения для в окрестностях рг+1,з,к, рг-1,з,к, Рг,3 + 1,к, Рг,3-1,к, Рг,3,к + 1, Рг,з,к-1 > оставив только члены до первого порядка малости включительно, выразив из разложений рг+2^,к, рг-2,з,к, рг,з+2,к, рг,з-2,к, рг,з,к+2, рг,ззк-2 и подставив их в уравнение (8), получим матричное уравнение:

рг,],к —

Ро кх 2 ку 2 Б г^,к

+

1

т (кх2ку 2 + кх 2к* 2 + ^ 2 ку 2) 2(кх2ку 2 + кх 2к* 2 + ^ 2 ку 2)

х ((рг,5,к+1 + рг,з,к-\)Ьг2ку2 + (рг,5+1,к + рг,з-1,к)к*2кх2 + (рг+1,5,к + рг-1,з,к)Ь,х2ку2). (9)

Уравнение(9) может быть решено итерационным методом, в таком случае получим выражение для

(г+1)

рг з к , где г — номер итерации:

р(г+1) — рг , з , к —

(г)

Ро кх 2 ку 2 Иг, з, к

1

т (к:2 ку2 + к:2 к2 + к2 ку2) + 2(кх2 ку2 + к:2 к2 + к2 ку2)

(г) (г) Ч 2<(Лг) <Лг) м- 2 г. 2

х ((р(Г,1к + 1 + р(%-1)к*2ку2 + (р(Г,з + 1,к + р(Г,з-1,к)к*2кх2 + (р(Г+1,з,к + р(-1,3,к)кх2ку2) .

(10)

Уравнение (10) должно решаться на каждом шаге по времени перед эйлеровым этапом.

На лагранжевом этапе вычисляются потоки массы через каждую из границ ячейки. При этом полагаем, что вся масса переносится через границу только за счёт нормальной к границе составляющей скорости. Запишем потоки массы для всех 6 граней ячейки: АМг^кк-2 — Р0йг,з,к-2кук*т,

РоЩ,з,к + 2 АМг,з-2 к — Р0^г,з - 2 к кук*Т, АМг,з + 2 к — Р0$г,3 + 2 к кук*т,

АМг,3,к + 2 —

АМг-2зк — р0йг-2 зккук*т, АМг+2 зк — р0йг+ 2 зккук*т. На заключительном этапе вычисляют- 2 - 2 2 2

ся газодинамические параметры в следующий момент времени:

4+1

иг ,з, к

й1з, к 0 + й!,з, к-Л

АМ4 2

г ,з, к-2

+ < з-1 'к Л2

АМ4 . 2 ,

г ,з-2 , к

+ Щ-1 з з к к Л3

АМ4 2

г-2, з, к

Р0к х кук

+

+ ■

йЪ , к + 1

АМг,з,к+2

+ й4,з + 1,к

х10у10*

АМг,з+2 ,к

+ й1+1,з,к

АМг+2 ,з,к

°г,з,к

й,з,к 0 + Язк-Л

АМ 2

г,з,к-2

Р0кхку к*

+ з-1,' к Л2

АМ 2

г,з-2 ,к

+ Ъ1-1,з,к Л3

АМ 2

г-2 ,з,к

Р0к х кук

+

+ -

^1г,з,к + 1

АМг,з,к+2

+ ^г,з+1,к

х10у10*

АМг,з+2 ,к

+ V.

г+1,.з,к(

АМг+ 2 ,з,к

Р0 кхку к*

4+1

з0 + ^Ьк-Л

АМ 4

+ й>1 к Л2

АМ4. 2 , г,з-2 ,к

+ ™дг-1,з,к Л3

АМ4 2

г-2 ,з,к

+

АМг,з,к+2

+ ™1з + 1,к

Р0к х ку к*

-я

г,з + 2 ,к

+

АМ4 .. !

+гьг+1,з,к

АМ4 2

г+ 2 ,з,к

Р0 к х

— Р0кх ку к* - (1 - Л1)

АМ 4

- (1 - Л2)

АМ 4

- (1 - Лз)

АМ4 2

г-2 ,з,к

-(1 - )

АМ 4

- (1 - )

АМк

- (1 - Лб)

АМг+2 ,з,к

где

—

Лз —

—

0, если V4 2 г,з,к-^ < 0,

1, если й4 2 г,з,к-^ > 0,

0, если 4 й 2-й г- 2 ,з,к < 0,

1, если 4 йй 2-й г- 2 ,з,к > 0,

1, если й)4 г,з+2 ,к < 0,

0, если й}4 г,з+2 ,к > 0,

Л2 —

ёл —

Лб —

0, если й4- 2 к г ,з-2, к < 0,

1, если й4- 2 к г ,з-2 з к > 0,

если й . . , 2 г,з,к+ 2 < 0,

0, если й . . , 2 г,з,к+ 2 > 0,

1, если 4 йй л_ 2 • и г+ 2 ,з,к < 0,

0, если 4 г + 2 ,з,к > 0,

0

2

Описанные выше вычисления выполняются на каждом шаге по времени. В рассматриваемой задаче моделирования движения тяжёлой несжимаемой невязкой жидкости существуют ячейки следующих типов: ячейки, заполненные жидкостью (ячейки среды), фиктивные ячейки, представляющие собой твёрдое тело, и вакуумные ячейки, не содержащие в себе жидкость и не принадлежащие твёрдому телу. Расчёты по приведённым выше формулам проводятся только для ячеек, в которых содержится жидкость. Значения в остальных ячейках вычисляются исходя из граничных условий на твёрдой стенке и на свободной поверхности. Положение свободной поверхности меняется с течением времени. Таким образом, количество расчётных ячеек может быть различным на каждом шаге по времени. Для определения ячеек, занятых жидкостью, в рассмотрение вводят маркеры [3]. Они представляют частицы, не имеющие массы. Изначально маркеры располагаются в ячейках расчётной области, занятых жидкостью. Скорость маркеров вычисляется интерполированием между значениями локальной скорости среды в соседних ячейках сетки, которые предварительно были найдены разностным решением.

Метод Давыдова является методом первого порядка точности. Проверка аппроксимации осуществляется путём собирания всех трёх этапов разностной схемы в одну схему и разложением в ряд Тейлора в окрестности фиксированной точки пространства-времени. Устойчивость разностной схемы для классического двумерного метода Давыдова исследовалась в работе [4], в которой показана неустойчивость метода в области малых скоростей. Обобщим критерий устойчивости из статьи на трёхмерный случай:

2

т т т

иот~ + уот~ + + 2со\1т

\\х Ну Н2 у Н

т2 т2 т2

* * - С0'1 Г"2 + Т2 + Т'.2 т2 т2 т \ Нх Ну Н

2

2

2

2

Н

2

у

Н

2

т т т

ио ~Г~ + *о Н— + ^0 Н-

Нх Ну Н^

< 1.

(11)

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

2. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ

Было рассмотрено несколько задач моделирования движения тяжёлой несжимаемой жидкости в трёхмерном пространстве с разными начальными и граничными условиями. Во всех задачах столб жидкости предполагался изначально неподвижным. Затем этот столб распадался, взаимодействуя с препятствием в центре расчетной области и со стенками «аквариума» (рис. 1).

а б

Рис. 1. Начало. Моделирование падения столба жидкости в трёхмерном пространстве: а — начальный момент времени; б — 600 шагов по времени

х

в г

Рис. 1. Окончание. в — 2000 шагов по времени; г — 2500 шагов по времени (вид сверху)

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

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

Рис. 2. Задача о падении столба жидкости. Начальный момент времени

0.0 0.5 1.0 1.5 2.0 2.5 3.0

а б

Рис. 3. Сравнение результатов моделирования: а — двумерный случай; б — трёхмерный случай

Библиографический список

1. Белоцерковский О. М., Давыдов Ю. М. Метод крупных частиц в газовой динамике. Вычислительный эксперимент. М. : Наука, 1982. 392 с.

2. Ливеровский Д. И., Шевырёв С. П. Метод Давыдова для случая несжимаемой невязкой тяжёлой жидкости на регулярной сетке // Математика. Механика : сб.

науч. тр. Саратов: Изд-во Сарат. ун-та, 2011. Вып. 13. С. 161-164.

3. Поттер Д. Вычислительные методы в физике. М. : Мир, 1975. 392 с.

4. Шевырёв С. П. Исследование двумерной схемы ме-

тода крупных частиц на устойчивость // Дифференциальные уравнения и теория функций в приложении к аэродинамике и теории вероятностей : межвуз. науч. сб. Саратов: Изд-во Сарат. ун-та, 1987. Вып. 7. С. 16-21.

Simulation of Incompressible Nonviscoussevere Fluid on a Regular Grid in Three-dimensional Space

D. I. Liverovskiy, S. P. Shevirev

Saratov State University, 83, Astrakhanskaya str., Saratov, 410012, Russia,, LiverovskiyDI@gmail.com, ShevSP@rambler.ru

This study focuses on modification of the method of Davydov (large particles) in the case of incompressible liquid. We consider the simulation of a heavy incompressible nonviscous fluid by Davydov's modified method in a three-dimensional case. Besides, comparison of the received results with a two-dimensional case is carried out. Formulas of modified Davydovs's method for the case of three spatial dimensions are lead out including difference analogue of the three-dimensional Poisson equation for the pressure. The criterion of stability is generalized. The description of the algorithm for constructing a three-dimensional free surface is provided. Implemented a three-dimensional computer graphics, showing development processes in online mode.

Key words: Davydov's method, fluid flow, method of large particles, free surface, Poisson equation for the pressure, the criterion of stability of difference schemes, symmetric flows.

References

1. Belotserkovskii O. M., Davydov Yu. M. Metod krupnykh chastits v gazovoi dinamike. Vychislitel'nyi eksperiment [Method of coarse particles in gas dynamics]. Moscow, Nauka, 1972, 392 p. (in Russian).

2. Liverovskiy D. I., Shevirev S. P. Metod Davydova dlia sluchaia neszhimaemoi neviazkoi tiazheloi zhidkosti na reguliarnoi setke [Davydov's method for a case of incompressible nonviscous heavy liquid on a regular grid]. Matematika. Mekhanika: sb. nauch. trudov, Saratov, Saratov Univ. Press, 2011, iss. 13, pp. 161-164 (in Russian).

YAK 519.6, 531

3. Potter D. Computational Physics. New York, John Wiley and Sons, 1973.

4. Shevirev S. P. Issledovanie dvumernoi skhemy metoda krupnykh chastits na ustoichivost' [Research of the two-dimensional scheme of a method of large particles on stability]. Differentsial'nye uravneniia i teoriia funktsii v prilozhenii k aerodinamike i teorii veroiatnostei: mezhvuz. nauch. sb. [Differential equations and the theory of functions in the Annex to aerodynamics and the theory of probability: a collection of papers]. Saratov, Saratov Univ. Press, 1987, iss. 7, pp. 16-21 (in Russian).

ПРИМЕНЕНИЕ МЕТОДА ГАЛЁРКИНА К РЕШЕНИЮ ЛИНЕЙНЫХ ЗАДАЧ ОПТИМАЛЬНОГО УПРАВЛЕНИЯ

И. А. Панкратов

Кандидат технических наук, доцент кафедры математического и компьютерного моделирования, Саратовский государственный университет им. Н. Г. Чернышевского, PankratovIA@info.sgu.ru

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

Ключевые слова: оптимальное управление, линейная система, метод Галёркина.

1. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим управляемую систему, описываемую линейным векторным обыкновенным дифференциальным уравнением:

^ = Лх + Ви, (1)

© Панкратов И. А., 2014