Математическое моделирование течений вязкого газа в многосвязном объёме с перфорированными стенками тема автореферата и диссертации по механике, 01.02.05 ВАК РФ

ИЫ4600248

На правах рукописи

УДК 533.6.011.3

Семакин Артём Николаевич

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ТЕЧЕНИЙ ВЯЗКОГО ГАЗА В МНОГОСВЯЗНОМ ОБЪЁМЕ С ПЕРФОРИРОВАННЫМИ СТЕНКАМИ

Специальность: 01.02.05 - Механика жидкости, газа и плазмы

1 ДПР 2010

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук

Ижевск-2010

004600248

Работа выполнена в Учреждении Российской академии наук Институте прикладной механики УрО РАН

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

Липанов Алексей Матвеевич

Официальные оппоненты: доктор физико-математических наук

Горохов Максим Михайлович

доктор физико-математических наук Якобовский Михаил Владимирович

Ведущая организация: Томский государственный университет

Защита состоится « 23 » апреля 2010 г. в 14:30 часов на заседании диссертацио ного совета ДМ 004.013.01 в Учреждении Российской академии наук Инстит прикладной механики УрО РАН по адресу: 426067, г. Ижевск, ул. Т. Барамзиной, 3

С диссертацией можно ознакомиться в библиотеке ИПМ УрО РАН

Автореферат разослан «18» марта 2010 г.

Ученый секретарь диссертационного совета д.ф.-м.н. О О ^

С.П. Копысов

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы:

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

Пористая среда с геометрической точки зрения имеет довольно сложное строение, и обычно течение жидкости или газа в такой среде исследуется с помощью перехода от истинных значений гидромеханических величин к фиктивным, которые «размазываются» по всей рассматриваемой среде. Далее, для этих фиктивных величин формулируется система уравнений гидромеханики. Получающаяся в результате этого математическая модель упрощает проблему исследования течения жидкости или газа через пористую среду и в некоторых частных случаях даже позволяет находить аналитические решения. Но все результаты получаются только для осред-нённого по объёму течения, получить какие-либо данные о поведении жидкости или газа на уровне пор невозможно.

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

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

Среди работ, посвященных изучению гидромеханики, можно выделить работы таких учёных, как Лойцянский Л.Г., Седов Л.И., Шлихтинг Г. Движению газа и жидкости в пористых средах посвящены работы Лейбензона Л.С., Полубариновой-Кочиной П.Я., Чарного И.А., Щелкачёва В.Н. Численные методы решения задач гидромеханики приводятся в работах Андерсона Д., Патанкара С., Роуча П., Флет-чера К., Самарского A.A., Липанова A.M.

Объект исследования: течение вязкого сжимаемого газа в многосвязных областях.

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

Цели диссертационной работы:

1. Разработка и реализация метода численного решения уравнений гидромеханики в рассматриваемых многосвязных областях.

2. Исследование течения вязкого газа в объёме с перфорированными стенками заполненном сферическими частицами.

Для достижения поставленных целей решаются следующие задачи:

1. Создание стандартного набора конечных объёмов, на которые можно разби вать объём, заполненный сферическими частицами.

2. Разработка криволинейной системы координат для каждого выделенного ко нечного объёма.

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

4. Организация взаимодействия соседних конечных объёмов и передачи дан ных между ними.

5. Создание и тестирование программного комплекса для расчёта течения вяз кого газа в объёме с перфорированными стенками, заполненном сфериче скими частицами.

6. Проведение расчётов для исследования течения вязкого газа в объёме с сферическими частицами.

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

Достоверность научных результатов и выводов подтверждается следующим:

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

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

На защиту выносятся:

1. Методика численного расчёта течения вязкого газа в объёме со сферичесг ми частицами.

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

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

4. Результаты параметрических расчётов течения вязкого газа в объёме со сф рическими частицами.

Научная новизна работы:

1. Разработана методика численного расчёта течения вязкого газа в многосвя ных областях.

2. Впервые проведено исследование поведения вязкого газа в объёме со сфер ческими частицами на основе численного решения уравнений гидромеха] ки.

Практическая значимость:

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

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

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

Апробация работы:

Материалы диссертации апробированы на следующих конференциях: «Динамика изследования - 2008» (София, 2008 г.), «Vöda: teorie а praxe - 2008» (Praha, 2008 г.), «Predni vödecke novinky - 2008» (Praha, 2008 г.), «Naukowy potencjal swiata -2008» (Przemysl, 2008 г.), «Nastoleni modern! vödy - 2008» (Praha, 2008 г.), «Актуальные проблемы науки в России» (Кузнецк, 2008 г.), «Краевые задачи и математическое моделирование» (Новокузнецк, 2008 г.), «Актуальные вопросы современной науки» (Таганрог, 2008 г.), «Актуальные проблемы современной науки» (Самара, 2008 г.). В целом диссертационная работа докладывалась на Учёном Совете Института прикладной механики УрО РАН (2009 г.).

Публикации:

Основные результаты диссертационной работы опубликованы в 12 статьях, из них 1 статья - в журнале, рекомендованном ВАК для публикации результатов диссертации на соискание учёной степени доктора и кандидата наук по механике и 2 статьи по перечню ВАК.

Структура и объём:

Диссертация состоит из введения, трёх глав, заключения и списка литературы. Общий объём диссертации составляет 145 е., включая 19 таблиц и 80 рисунков. Список литературы содержит 100 наименований.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

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

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

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

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

Во второй главе рассмотрен численный метод решения задачи обтекания вязким газом сферических частиц в ограниченном объёме. Согласно данному методу сначала в исходной области вводится некоторая глобальная декартовая система координат (.X, Y, Z). Далее эта область делится на N подобластей или конечных объ-

ёмов (КО). Для объёма с перфорированными стенками, заполненном сферическими частицами, все эти КО можно свести к пяти стандартным видам (рисунок 1): прямоугольный, сферический, пирамидальный, кольцевой и цилиндрический. Прямоугольный объём представляет собой обычный прямоугольный параллелепипед. У сферического КО одной из шести граней является поверхность сферы, а боковые рёбра направлены перпендикулярно ей. Пирамидальный КО соприкасается одновременно с двумя телами, он имеет клинообразную форму и одна его вершина является точкой касания тел. У кольцевого КО одна грань имеет цилиндрическую форму. Цилиндрический КО представляет собой обычный цилиндр, используемый для описания входа и выходов.

а) б)

Рисунок 1 Типы конечных объёмов а) прямоугольный; б) сферический; в) пирамидальный; г) кольцевой; д) цилиндрический

В каждом КО вводится локальная декартовая система координат (JICK) (х, у, z). Каждая такая локальная система координат определяется заданием начала координат o(al,al,al) и базисных векторов i = (а,1, а,2, а\), \ = {a\,al,a]), к = (а', а], а]) в глобальной системе координат (х, У, Z).

Локальные (x,y,z) и глобальные (X,Y,Z) координаты точки связаны соотношениями:

a) (x,y,z)->{X,Y,Z) б) {XJ,z)->(x,y,z)

X = а\ + а\х + а\у + a\z , х = а\Х + a¡Y + a\Z - a¡a¡¡ - a¡üg - a¡üg,

Y = a20 + afx + a\y + a]z , (1) y = a\X + a\Y + a\Z - a\al - a\a\ - a\a\, (2)

Z = al + a\x + a\y + a\z , z = a\X + a¡Y + a\Z - a]a¿ - a]al - a]al.

Во всех видах КО, кроме прямоугольного, помимо декартовой ЛСК необходимо также вводить криволинейную систему координат в которой данный КО

можно представить в виде параллелепипеда:

1) сферический КО

Í = tJx2+у2+Z2 , т| = —у/х, q = -z/x, (3)

де 4 - расстояние от начала ЛСК до данной точки, -1 < г| < 1, -1 < <; < 1; 2) пирамидальный КО

$=(х2+у2+г2]

¿Х2+22

2+у2+г2),с; = -ф

(5)

, У] = 2у/(х

де 0 < г\ < 1/г, -1 < $ < 1, г - радиус сферы;

3) кольцевой КО

1 = 4хг + г2 , ц = у , <; = -2/х,где -1<<;21;

4) цилиндрический КО

4 = х, г| = -¡[у2 + г2 , с, = ап^(г/у)+ яЛ;, где -я < <; < я. (6)

Далее для каждого полученного КО в его локальной системе координат форму-ируется система уравнений гидромеханики в безразмерной форме, в которой кон-ективные члены представляются в симметричном виде. Эта система включает равнение неразрывности, три уравнения импульса и уравнение энергии. В локаль-ой декартовой системе координат она имеет вид:

др 1 — + —

3/ 2

дх ду дг дх ду

ди 1 —+—

д( 2

ди ди ди ди диу дш

и——+\У—+--ь—н--

дх ду дг дх ду дг

дг ) 2 1

рЫе

(7)

2 1 ^ д2и д2и д2и"

1

1

--у)т--

ркМ2 дх 2 ЗрЫе дх

Му),

ду_ ]_ д1 + 2

ду ду ду диу ду2 Эпе

и—+у—+м>—+—+—+-

дх ду дг дх ду да

1

рЛе

дх2 ду2 дг2

дх2+ ду2 + дг2

(8)

.1

ркМ2 ду 2 4 ' ЗрЯесу 4 '

ды 1 2

д\у дю д\у ди\у дш дм

и— + у— + м/—+-+-+ —

дх ду дг дх ду дг

рЯе

(д2уц д2\у

Их2

(9)

ркМ2 дг 2 К ' ЗрЯе5г 4 '

(10)

дТ ^

а 2

дТ дТ дТ диТ дуТ дwT

и-+ У-+ М>-+-+--!-——

дх ду дг дх ду дг

/

р У 2 К ' рЯе

/

ди

дх

рЯеРг 2

+ 2

гд2Т д2Т д2Тл

дг2 2

дх2 ду2

ди ду ду дх

\2

ди дн> , — + —| +

дг дх

ду д\у — + —

дг ду

\2

~(МУ))2

ду

у

+2

/ {дг,

(П)

Здесь V = (м,у, м) - вектор скорости, р,р,Т - давление, плотность и температу

ра, сй'у(у) = —+ —+ — - дивергенция вектора скорости, Яе = р0С/0/г/ц - числ дх ду дг

Рейнольдса, Рг = с^ц/^ = 1 - число Прандтля, М = и0/с0 - число Мах

со = кро/ро, к = Ср/Су=1Л - отношение изобарной и изохорной теплоёмкостей.

Коэффициенты вязкости ц и теплопроводности X полагаются константам Все переменные безразмерные. В качестве масштабов берутся: к - диаметр входно го отверстия, 1/0 - максимальная величина продольной компоненты скорости в по перечном сечении входа, р0, р0 - характерные давление и плотность, соответст вующие ио.

Безразмерное уравнение состояния:

Р = РТ. (12)

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

д1¥ д1Г. дЖ д!Г

-=--Г|х+-(13)

дх 8Е, дц дс,

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

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

При нахождении частных производных по пространственным направлениям узлах разностной сетки используется метод неопределённых коэффициентов, позво ляющий рассчитывать эти производные с любым порядком точности. Согласно это му методу для определения производной в данной точке необходимо знать значеш искомой функции в нескольких соседних узлах. Поэтому при расчёте производных точках, расположенных около границы КО необходимо «заходить» в соседние КО Обычно разностные сетки в двух смежных КО несогласованны, т.е при «заходе» соседний КО нужные нам точки могут не совпадать с расчётными точками разност ной сетки этого КО (см. рисунок 2).

I I I

Рисунок 2 Область соприкосновения двух соседних КО расчётная точка КО, + - интерполируемая точка КО

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

В данной работе для проведения интерполяции использовался метод, основанный на представлении функции в виде отрезка ряда Фурье по ортогональным многочленам. Алгоритм метода следующий.

Неизвестное значение величины / в точке х0 = (х0,у0,20) определяется через

её известные значения в точках хт = (х„,ут,гт), т = \,Ы с помощью отрезка ряда Фурье, построенного на основе системы ортогональных многочленов:

/(хо)= 2>/Л*(хо)> (14)

V'

/+./'+*=о

где

^ =Ы4/(х)). (15)

Система ортогональных многочленов (х) строится последовательно по формулам:

о®

Фооо=1> (17)

Фу* (х) = * (х) + УУ.нк (х) + (х) + £ оси„,Фтп/ (*), (1В)

и

ат„, = + 21|/,ум(х), (19)

М - множество уже созданных многочленов.

Скалярное произведение и норма определяются следующим образом:

(/(хЫх)) = 2/(хиЫхя), ||/(х)И(/(х),/(х)). (20)

Точки хт выбираются таким образом, чтобы ||ф^(х)||> Ю-4 для всех Ф,-,4(х).

Также в работе полагалось п = 2, N = 10, что эквивалентно интерполяции полиномом второй степени.

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

и = а\и + а\V + а\м>, V = а\и + + а]м, Ж = а\и + а\у + а]ы, (21)

и далее из глобальной системы координат в локальную, но уже другого КО:

и = а\и + а\У + Й13Ж, V = а\и + а\У + аг21У, = а\и + а]У + а\Ш, (22) коэффициенты перехода а{ определены выше.

В качестве начальных условий берутся следующие значения:

и = V = w = 0, р = р = Т = 1.

Граничные условия на границах многосвязной области имеют вид:

а) поверхность тел

и — V = w = 0 - условия прилипания, дТ/дп = 0 - условие адиабатичности,

б) вход

u = f(r,b0)-cr{p-pi)-f(r,bs), v = w = 0, +

где с, = 1¡кМ, рх - среднее по входному сечению давление, 50 = 1 - заданная величина энтропийной функции в ядре потока, ô0 = 0.2 - толщина динамического пограничного слоя, 5s =0.1 - толщина теплового пограничного слоя,

, ч il, r<R-b /V,5) = « , .„ , R = 0.5 - радиус входного отверстия,

[l-((r-Ä + 8)/5)\ r>R-b

в) выход

Задаём только одно условие - величину давления на выходе из объёма в следующей форме:

L

p = \ + c1(m{L,t)-msr{tty,T№ m(x,t)= Jjp udydz, msr(t) = —^m{x,t)dx, c2=0.2.

S(r) о

Если какие-либо грани KO выходят на границу многосвязной области, то для этих граней ставится одно из указанных выше граничных условий.

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

В начале главы для подтверждения правильности работы изложенного выше метода конечных объёмов приведена тестовая задача — задача обтекания сферы. В этом случае расчётная область представляет собой шар радиуса 7.8, в который помещена сфера радиуса 0.5. Данная расчётная область делится на 6 сферических КО. В качестве критерия правильности расчётов рассматривался коэффициент сопротивления сферы Сх. Число Маха M бралось 0.1. В этом случае результаты расчётов должны совпадать с данными для несжимаемой жидкости и, следовательно, значения Сх, полученные путём расчётов, можно сравнивать со стандартной кривой коэффициента сопротивления сферы несжимаемой жидкости, аппроксимируемой зависимостью:

Сх = —(l + 0.25VRß + 0.0117Re), 1 < Re < 1000. (23)

Re

Результаты расчётов представлены на рисунке 3. Из этого рисунка видно, что полученные расчётные значения хорошо соответствуют стандартной кривой коэффициента сопротивления.

При Re < 250 были измерены длина отрывной зоны за сферой и угол отрыва потока (таблица 1). Здесь Р - значения, полученные в данной работе, далее идут результаты, приведённые в работах других авторов (1 - Гущин В.А., Матюшин П.В., 2 - Jenson V., 3 - Taneda S., 4 - Rimon Y., Cheng S.I.). Из таблицы 1 видно, что полученные результаты удовлетворительно совпадают с данными других работ.

;

V

ч

О 200 400 600 800 1000

Рисунок 3 Зависимость коэффициента сопротивления от Re 1 - расчёт, 2 - стандартная кривая (23)

Таблица 1

Длина отрывной зоны I и угол отрыва потока а

Re / а,0

Р 1 2 3 4 Р 1 3 4

40 0.28 0.29 0.20 0.32 0.34 35 36.5 - 35

50 0.40 - - 0.47 - 41 - ' - -

100 0.80 0.77 - 0.91 0.86 52 53.8 53 53

250 1.12 - - - 65 - - -

Далее рассматривалось течение вязкого газа в объёме с перфорированными стенками, в который помещалось от I до 64 сфер радиуса 0.5. Объём представляет собой прямоугольный параллелепипед. Расчёты проводились при Re = 100, 500 и М= 0.6. Для 1-3 сфер течение при Re = 100 получается стационарным, а при Re = 500 - нестационарным, периодическим. Для 4-64 сфер течение в обоих случаях стационарно. Ниже приведены рассматриваемые случаи положения сфер и некоторые результаты расчётов.

Для одной сферы рассматривались три варианта ей положения:

1) сфера находится в центре канала;

2) сфера опускается на расстояние одного радиуса от центра канала;

3) сфера касается нижней плоскости канала.

Размеры области: длина - 3, ширина - 3, высота - 3. Вход располагается в центре передней грани, два выхода находятся на вертикальной оси симметрии задней грани. Вход й выходы - цилиндры радиуса 0.5.

На рисунке 4 приведены графики Cx{t), а на рисунке 5 представлены примеры полей скоростей в горизонтальной плоскости симметрии JÎZ, Y= 1.5 для варианта 1 (сфера находится в центре объёма).

При Re= 100 течение получается стационарным (рисунок 4.а). Из рисунка 5.а видно, что газ втекает в расчётную область в виде явно очерченной струи, которая сохраняет свою форму до столкновения со сферой. В процессе её обтекания струя приобретает форму купола, края которого направлены к выходам из объёма. При этом скорость движения газа уменьшается. В рассматриваемом горизонтальном сечении XZ, Y= 1.5 около боковых стенок объёма и в кормовом пространстве сферы находятся отрывные зоны. Вихри, расположенные около левой и правой боковых стенок, по размеру вдоль оси х равны длине расчётной области и занимают про-

странство от передней до задней граней рассматриваемого объёма. Два других вих ря занимают всё пространство между сферой, задней гранью объёма и набегаюшт потоком газа, т.е их длина вдоль оси х равна 1, а вдоль оси z - 0.5. Можно указав также, что поле скоростей симметрично относительно оси х. Аналогичные вихри ^ вертикальном сечении XY, Z = 1.5 отсутствуют.

При Re = 500 течение является нестационарным и периодическим (рисунок 4.6} Если рассматривать горизонтальную плоскость симметрии области XZ, Y= 1.5 (ри сунок 5.6), то картина течения следующая. Газ входит в расчётную область в вид струи, движется по направлению к сфере, наталкивается на неё и приобретает купо' лообразную форму, продолжая двигаться к выходам. Через определённые проме^ жутки времени по обе стороны от входящего потока газа образуются небольшие вихри, которые, далее, движутся вдоль потока газа от передней стенки объёма i задней и объединяются с находящимся там стационарным вихрем. Одновременно j каждой стороны струи может находится до трёх вихрей: только что образовавшийся около передней стенки, движущийся по направлению к задней грани и стационар! ный вихрь около задней стенки. Данная система вихрей заменяет большие вихри занимавшие всё пространство между передней и задней гранями объёма при Re = 100. Между сферой и задней гранью рассматриваемой области как и пр^ Re = 100 располагается отрывная зона, но в данном случае она имеет более сложно« строение. В частности, в месте отрыва набегающей струи от сферы периодически образуется маленький вихрь. Далее, он отрывается от сферы и движется к задней стенке объёма.

На рисунке 4.6 W-образный минимум отражает отрыв вихрей сначала от перед ней стенки, потом от поверхности сферы. Одинарный минимум соответствует нача лу образования вихря у поверхности сферы в месте отрыва потока, а небольшая

Рисунок 4 График Сх = при a) Re = 100, б) Re = 500

a) Re = 100 б) Re = 500

Рисунок 5 Поле скоростей в сечении XZ, Y = 1.5 для варианта 1

На рисунке 6 приведена разностная сетка для вариантов 2 и 3. Суммарное количество точек по всем КО равно 90168 (вариант 2) и 85486 (вариант 3), в пирамидальном КО - 456, в кольцевом - 700, в сферическом - 5648.

Рисунок 6 Разностная сетка для вариантов 2 (а) и 3 (б)

Для случая двух сфер рассматривались четыре варианта их взаимного положения:

1) центры сфер расположены на оси симметрии объёма, при этом сферы касаются друг друга;

2) центры сфер расположены на оси симметрии объёма, расстояние между сферами равно один, три, шесть и девять радиусов;

3) сферы касаются друг друга и нижней грани объёма;

4) расстояние между сферами равно два радиуса, сферы расположены на нижней грани объёма.

Размеры области Q в вариантах №1, 2 брались следующими: длина - 5.5, ширина - 3, высота - 3. Исключение составляет случай, когда расстояние между сферами равно 6 и 9 радиусов, здесь длина Q равнялась 7 и 8.5. Вход располагается в центре передней грани, два выхода находятся на вертикальной оси симметрии задней грани. Вход и выходы - цилиндры радиуса 0.5. Координаты центра первой сферы -(1.5;1.5;1.5).

В таблице 2 приведены некоторые, параметры течения для вариантов №1,2 при Re = 100. В этой таблице / - расстояние между сферами (радиусы), р\ и р2- значения давления в передней и задней критических точках первой сферы, р3 и - значения давления в передней и задней критических точках второй сферы, Сх\ и Сх2 - коэффициенты сопротивления первой и второй сфер, соответственно, а„ - угол отрыва потока от первой сферы, ар - угол, отсчитываемый от задней критической точки первой сферы до точки минимума давления на поверхности сферы, рт„ - минимальное давление на сфере.

Из таблицы 2 следует, что при некотором /е(1;3) 0,2 = 0, т.е. главный вектор сил, действующих со стороны набегающего потока газа на вторую сферу, равен нулю. Но данное положение неустойчиво, т.к. при любых отклонениях от него на сферу начинает действовать сила, стремящаяся увести её от данного положения.

На рисунке 7 приведены поля скоростей в вертикальной плоскости симметрии XY, Z= 1.5 для Re = 100 (стационарное течение) и Re = 500 (нестационарное периодическое течение) при I = 1.

Таблица 2

Некоторые параметры течения при 11е = 100

ь Р1 Рг Рз Ра Сх 1 Сх 2 «и, 0 Рт'т

0 2.17 1.33 1.33 1.37 - - 67 114 1.17

1 2.17 1.35 1.34 1.37 1.03 -0.10 73 119 1.17

3 2.18 1.38 1.38 1.37 0.92 0.10 81 119 1.20

6 2.19 1.39 1.40 1.37 0.94 0.15 81 119 1.21

9 2.19 1.39 1.40 1.38 0.94 0.16 81 119 1.21

а) К.е = 100 б) Ые = 500

Рисунок 7 Поле скоростей в сечении ХУ , 2 = 1.5 для двух сфер

На рисунке 8 приведены коэффициент давления и распределение температуры по поверхности первой сферы, а на рисунке 9 - по поверхности второй сферы при Яе = 100. Коэффициент давления определяется формулой:

(24)

~кМ2 2

Для первой сферы графики ср и Т при изменении 1 различаются незначительно. Для второй сферы с ростом I точки минимума и максимума ср постепенно меняются местами, а профили Г сначала растут (/ = 0, 1, 3), потом падают (/ = 6) и снова начинают расти (/ = 9). Также на рисунке 9 можно видеть, что при переходе от / = 6 к / = 9 значения ср и Т меняются значительно в меньшей степени, чем при остальных изменениях /. Это относится и к данным таблицы 2.

Рисунок 8 Распределение давления (а) и температуры (б) по поверхности первой сферы для различных / при 11е = 100

1.53

1.31

1.38

1.60

1.46

=d-1-1-1—-1 > 1.224 -1-1-H-h

■180 -90 0 90 180 6) "180 -90 0 90 180

Рисунок 9 Распределение давления (а) и температуры (б) по поверхности второй сферы для различных I при Re = 100

£1

1.224

3

а)

Когда в объёме располагалось три сферы, рассматривались три возможных варианта их положения:

1. Центры сфер располагаются в точках с координатами (1.5;1.5;2.0), (3.0; 1.5; 1.0), (4.5;1.5;2.0), т.е. центры сфер лежат в плоскости у = 1.5 и расстояние между сферами вдоль оси х равняется одному радиусу. Размеры объёма: длина - 6.0, ширина - 3.0, высота - 3.0.

2. Центры сфер располагаются в точках с координатами (1.5;1.5;1.78868), (2.0;1.5;0.92265), (2.5; 1.5; 1.78868), т.е. они лежат в вершинах равностороннего треугольника со стороной, равной диаметру сфер, в плоскости у = 1.5. Размеры объёма: длина-4.0, ширина-3.0, высота-3.0.

3. Центры сфер располагаются в точках с координатами(1.5;1.5;1.78868), (1.5;2.0;0.92265), (1.5;2.5;1.78868), т.е. они также лежат в вершинах равностороннего треугольника со стороной, равной диаметру сфер, но в плоскости, параллельной передней грани, х = 1.5. Размеры объёма: длина - 3.0, ширина - 3.0, высота - 4.0.

Здесь полагалось, что передний левый нижний угол объёма расположен в точке с координатами (0;0;0). Вход расположен в центре передней грани, два выхода лежат на вертикальной оси симметрии задней грани.

На рисунке 10 приведены поля скоростей в горизонтальной плоскости симметрии XI, Г= 1.5 для Г1е = 100 (стационарное течение) и Re = 500 (нестационарное периодическое течение) для варианта №1.

При Re = 100 для варианта №1 коэффициенты сопротивления сфер имеют следующие значения: Сл = 0.72, Сх1 = 0.63, Сх3 = 0.13, т.е. наибольшему воздействию со стороны набегающего потока газа подвергаются первая и вторая от входа сферы. В точках лобовой поверхности этих сфер, расположенных под углами 24° и 13°, отсчитываемых от передних критических точек первой и второй сфер, соответственно, в сечении У= 1.5, давление достигает своего максимума: ртш 1 = 1.99 - для первой сферы, ртаХ2 = 1.85 - для второй сферы. Далее, в точках, расположенных под углами 96° и 114°, отсчитываемых от задних критических точек первой и второй сфер, соответственно, давление принимает минимальные значения ртш 1 = 1.12 для первой сферы и рт„а = 1.27 для второй сферы. Они располагаются перед соответствующими точками отрыва потока (67° и 63°). В остальной области, включая пространство около третьей сферы, давление находится в интервале (1.3;1.4) и резко падает до 1 в выходных отверстиях из объёма.

а) 11е = 100 б) Ле = 500 I

Рисунок 10 Поле скоростей в сечении Х2, У = 1.5 для трёх сфер

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

Рисунок 11 Расчётная область для четырёх сфер

Здесь объём имеет один вход и четыре выхода. Радиус входа - 0.5, радиус выходов - 0.35, т.е. общая площадь выходов в два раза больше площади входа. Размеры объёма следующие: длина - 4.0, высота - 1.81650, ширина - 1.86603. Размеры расчётной области подобраны таким образом, чтобы его боковые стенки касались пирамиды, составленной из сфер.

Сферы, составляющие пирамиду, обозначаются следующим образом: верхняя сфера называется сферой Я, сфера основания, касающаяся левой грани параллелепипеда, - сферой Е, сферы основания, касающиеся правой стенки параллелепипеда, -Р и й, соответственно. Начало координат ГСК располагается в центре тяжести треугольника с вершинами в центрах сфер основания, ось х направлена в сторону выходов, ось_>• - вверх. Центры сфер имеют координаты (г = 0.5 - радиус сфер):

\ г

Е\ 0;0;—¡=г

I & .

л/3

Я

л/3

Ниже приведены некоторые результаты исследования течения вязкого газа для случая четырёх сфер.

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

мере, приближаются к пирамиде и отходят к боковым граням. Около верхней и нижней граней параллелепипеда вдоль лобовой поверхности сфер располагаются два вихря с осью вращения, параллельной оси г. В пустое пространство внутри пирамиды газ попадает по криволинейному каналу, образованному сферами Я, £ и О. Это вполне объяснимо, поскольку именно этот канал расположен по ходу движения струи. По каналу, ограниченному сферами основания Е, Г, (?, газ из данной области движется вниз. Одновременно газ движется в сторону выхода из пирамиды по каналам с границами, образованными сферами Н, Е, Р и Я, .Г, б.

За пирамидой располагается область, заполненная многочисленными вихрями с осями вращения, параллельными как оси у, так и оси г. С ростом числа Рейнольдса количество вихрей и их интенсивность растут (ютах = 5.77 при Яе = 100 и югаах = В.77 при Яе = 500 в пространстве за пирамидой 3<х<3.8, где сошах - максимум модуля завихренности). Данная область смещена преимущественно к основанию пирамиды. В кормовом пространстве верхней сферы Я располагаются два вихря, вытянутые к задней грани расчётной области в виде хвоста. С ростом числа Рейнольдса данные вихри удлиняются, одновременно растёт их интенсивность.

На рисунках 12-13 представлены поля скоростей в горизонтальном сечении плоскостью, проходящей только через сферы основания (рисунок 12) и все четыре сферы одновременно (рисунок 13).

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

В центре пирамиды максимальное значение скорости равно 0.8 для Яе = 100 и 1.1 для Яе = 500, т.е. величина скорости имеет тот же порядок, что и при входе в объём. Давление изменяется в интервале (2.3;2.5) для Яе = 100 и (1.4;1.5) для Яе = 500, температура - в интервале (1.4;1.5) для Яе = 100 и (1.1;1.2) для Яе = 500.

В области перед пирамидой давление лежит в интервале (2.7;3) при Яе = 100 и (1.5;1.7) при Яе = 500, после пирамиды - в интервале (1.7;2) при Яе = 100 и (1.2;1.3) при Яе = 500. В выходах оно резко падает до 1.

Температура распределена по объёму более равномерно: при Яе = 100 Г е (1,4;1.5), при Яе = 500 Т е (1.1;1.2). Только при входе газа в объём наблюдается понижение его температуры.

а) Яе = 100 б) Яе = 500

Рисунок 12 Поле скоростей в сечении XI, У = -0.26 для четырёх сфер

а) Яе = 100

б) Яе = 500

Рисунок 13 Поле скоростей в сечении Х7., У = 0.41 для четырёх сфер

?п

а)хг , Г = 0.1

б)ху , г = 0.05

Рисунок 14 Поле скоростей в центре пирамиды при 11е = 500 в плоскостях

а) хг, у = 0.1,6) ху , г = 0.05

Далее рассматривается неплотная кубическая упаковка сфер (4-64 сферы), когда в первом слое каждая сфера касается только четырёх соседних сфер, и все последующие слои повторяют первый. Расчётная область <2 представляет собой прямоугольный параллелепипед с одним входом и четырьмя выходами цилиндрической формы. Вход размещается в центре передней грани, выходы - в центрах квадрантов задней грани, образованных её горизонтальной и вертикальной осями симметрии. Пример расчётной области приведён на рисунке 15 для упаковки (3;3;3), где сначала указывается количество поперечных слоев, а затем количество сфер вдоль осей у иг. Число Рейнольдса бралось 25,100 и 500, число Маха равно 0.15.

Средние характеристики течения приведены в таблице 3. Здесь М= (Мх\Му,М2) -упаковка сфер, Мх - число сфер в направлении оси х, Му- число сфер в направлении | оси у, М2 - число сфер в направлении оси г, N - суммарное по конечным объёмам количество точек разностной сетки. Далее идут средние по поперечному сечению

значения давления р, модуля вектора скорости |у| и продольной компоненты скорости и . Индексы 1 и 2 относятся к поперечным сечениям, расположенным перед и после упаковки сфер на расстоянии 0.03, соответственно.

Сходимость по количеству точек разностной сетки исследовалась при всех Яе на упаковках (1;2;2) и (2;2;2). При Ке = 25 на упаковке (1;2;2) рассматривалось три сетки. При переходе от сетки №1 (23716 точек) к сетке №2 (66004 точек) макси-| мальное изменение средних величин составило 3.7%, а при переходе от сетки №2 к сетке №3 (203820 точек) - 1.7% (см. таблицу 3). В остальных случаях рассматривались только сетки, соответствующие сеткам №2 и №3 данного случая. Разница меж-

ду средними величинами на наименьшей и наибольшей разностных сетках не превосходит 2.2% для Яе = 25,100 и 4.7% для Яе = 500.

Физическая адекватность расчётов проверялась с помощью закона сохранения массы, согласно которому для стационарного течения в отсутствии источников и стоков должно выполняться условие /771 = т2, где т - массовый расход в поперечном сечении перед (1) и после (2) упаковки. Для всех рассматриваемых упаковок разница между Ш1 и тг не превосходит 3.4%, причём ошибка растёт с увеличением числа сфер. Данная погрешность расчётов при выполнении закона сохранения массы образуется за счёт трёх факторов: погрешность численного интегрирования при вычислении массового расхода, погрешность интерполяции (значения гидромеханических переменных передаются из одного конечного объёма в другой с помощью полиномиальной интерполяции, которая реализовывалась без учёта законов сохранения), погрешность разностной схемы (разностная схема построена на основе системы уравнений гидромеханики в симметричной форме и, следовательно, не является консервативной, т.е. реализует интегральные законы сохранения с некоторой ошибкой).

Таблица 3

Средние характеристики потока в неплотной кубической упаковке сфер

м N Рх Рг V 1 V 2 "1 и2

Яе = 25, М= 0.15

0;2;2) 23716 1.297 1.056 0.288 0.313 0.185 0.217

66004 1.283 1.054 0.288 0.320 0.187 0.225

203820 1.295 1.053 0.293 0.321 0.188 0.227

(1;3;3) 147534 1.179 1.010 0.340 0.281 0.185 0.218

(1;4;4) 261620 1.162 1.011 0.430 0.376^ 0.184 0.219

(1;5;5) 408262 1.131 0.992 0.500 0.369 0.183 0.217

(1;б;б) 587460 1.107 0.975 0.576 0.448 0.185 0.220

(2;2;2) 95484 1.559 1.068 0.285 0.360 0.188 0.272

271276 1.592 1.069 0.290 0.365 0.188 0.277

(3;3;3) 283070 1.571 1.017 0.341 0.370 0.184 0.290

(4;4;4) 627236 1.678 1.019 0.425 0.466 0.183 0.305

Яе = 100, М= 0.15

(1;2;2) 66004 1.108 1.038 0.321 0.281 0.185 0.199

203820 1.107 1.034 0.327 0.282 0.185 0.202

(1;3;3) 147534 1.050 0.997 0.359 0.247 0.183 0.197

(3;3;3) 283070 1.146 1.001 0.364 0.287 0.182 0.214

Яе = 500, М= 0.15

(1;2;2) 125570 1.065 1.028 0.428 0.455 0.183 0.191

212676 1.059 1.024 0.413 0.450 0.184 0.200

(1;3;3) 278640 1.013 0.991 0.392 0.353 0.178 0.186

(3;3;3) 511184 1.049 0.994 0.408 0.350 0.181 0.195

Далее изложим полученные результаты. Из таблицы 3 видно, что при переходе! через упаковку сфер средняя скорость газа в продольном направлении й возрастает, причём она зависит только от числа поперечных слоёв упаковки Мх. С ростом Мх также увеличивается перепад давления р^-р2, т.к. газу необходимо преодолевать1 всё большее расстояние между рассматриваемыми поперечными сечениями. Но при возрастании Му и М2 перепад давления уменьшается, что сопровождается ростом|

|у[. С увеличением Яе значения р{ - р2 падают, т.к. требуется меньше усилий для1

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

Характер течения следующий. Газ входит в объём в виде чётко очерченной' струи. Наталкиваясь на упаковку сфер, поток разбивается на несколько меньших! струй, которые движутся по направлению к выходам через криволинейные каналы, образованные сферами. При этом даже при Яе = 25 в пространстве между сферами; наблюдаются небольшие вихри (рисунок 16). На рисунке 17 приведе ны профили] продольной компоненты скорости для упаковки (3;3;3) в поперечном сечении между вторым и третьим слоями вдоль оси у. Там же приведён данный разрез, на котором! пунктиром указана линия измерения с координатой 2= 1. Из рисунка 17 видно, что1 скорость движения центральных струй больше пристеночных. Но с ростом Яе максимальная скорость в центре падает, а около стенок растёт. При этом пики графика смещаются к стенкам.

Рисунок 15 Рассматриваемая область (вид сбоку)

Рисунок 16 Поле скоростей в продольном сечении 2 = 1.11 для упаковки сфер (3,3,3) при Яе = 25

в поперечном сечении X = 2.32 вдоль оси у для упаковки (3,3,3), 1-Яе = 25,2-Б1е = 100,3-Яе = 500

ЗАКЛЮЧЕНИЕ

Основные результаты, полученные в диссертационной работе, состоят из следующих положений:

1. Разработан метод численного решения уравнений гидромеханики в многосвязных областях на примере прямоугольного объёма с перфорированными стенками, заполненного сферическими частицами. Создан набор стандартных типов конечных объёмов (КО), на которые можно разбить исходную область со сферическими частицами. Для каждого конечного объёма разработана криволинейная система координат. Для организации взаимодействия между КО предложено три метода интерполяции гидромеханических переменных (метод, основанный на вычислении определителей; метод, основанный на разложении функции по формуле Тейлора; метод, основанный на представлении функции в виде отрезка ряда Фурье по ортогональным многочленам). При проведении расчётов разностная схема, построенная исходя из дивергентной формы записи системы уравнений гидромеханики, оказывалась неустойчивой при Яе > 100. Устойчивость достигнута лишь при переходе к системе уравнений гидромеханики, в которой конвективные слагаемые записываются в симметричной форме.

2. На основе изложенных в главе 2 теоретических положений создан программный комплекс для проведения расчётов течения вязкого газа в объёме с перфорированными стенками, заполненном сферическими частицами. Он протестирован на задаче обтекания сферы в неограниченном объёме в диапазоне чисел Рейнольдса 40-^1000. Полученные значения коэффициента сопротивления сферы, угла отрыва потока и длины отрывной зоны удовлетворительно совпадают с экспериментальными данными и результатами расчётов других авторов.

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

4. В работе рассмотрено течение вязкого газа в объёме с перфорированным стенками, в котором располагаются одна, две, три и четыре сферы в различных по ложениях, а также неплотная кубическая упаковка сфер (4-64 сферы) при Re = 100 Re = 500. Для каждого случая проанализирован характер формирования вихрей, об разующихся отрывных зон и движения вязкого газа в целом для объёма в рассмат риваемых условиях. В итоге получается, что, когда в объёме располагаются одн две или три сферы, течение газа при Re = 100 становится стационарным, а пр Re = 500 - нестационарным, периодическим. Когда в объёме располагаются четыр и более сфер, течение является стационарным как при Re = 100, так и при Re = 50 Во всех случаях при переходе от Re = 100 к Re = 500 наблюдается увеличение коли чества вихрей, они становятся более интенсивными и усложняется сама их структу ра.

Основные результаты диссертации опубликованы в работах:

1. Липанов A.M., Семакин А.Н. Обтекание вязким газом сферических частиц в ог раниченном объёме // Вестник Удмуртского университета. Математика. Механи ка. Компьютерные науки. 2009. Вып. 4. С. 79-86.

2. Липанов A.M., Семакин А.Н. Обтекание вязким газом системы двух сфер объёме с перфорированными стенками // Математическое моделирование. 2009 Т. 21, №7. С. 67-74.

3. Липанов A.M., Семакин А.Н. Обтекание трёх сфер потоком вязкого газа пр Re=100 // Вестник ИжГТУ. 2008. № 4. С. 203-205.

4. Липанов A.M., Семакин А.Н. Применение метода конечных объёмов к задач обтекания сферы // Материали за 4-а международна научна практична конферен ция «Динамика изследования - 2008». 1.21. Математика. Съвременни техноло гии на информации. Здание и архитектура. София: «Бял ГРАД-БГ» ООД, 2008 С. 31-35.

5. Липанов A.M., Семакин А.Н. Структура течения, распределение давления и тем пературы при обтекании сферы, расположенной в объёме с перфорированным стенками, потоком газа // Materiäly IV mezinärodni vödecko-praktickä konferenc «VSda: teorie a praxe - 2008». Dil 14. Matematika. Vystavba a architektura. Praha Publishing House «Education and Science» s.r.o., 2008. С. 39-42.

6. Липанов A.M., Семакин А.Н. Течение вязкого газа в трёхсвязной области (дв сферы в прямоугольном объёме) // Материалы международной научно практической конференции «Актуальные проблемы науки в России». Выпуск V Т.З. Кузнецк: КИИУТ (филиал ПТУ), 2008. С. 101-103.

7. Липанов A.M., Семакин А.Н. Некоторые параметры течения при обтекании вяз ким газом системы двух сфер // Materialy IV Mi^dzynarodowej naukowi praktycznej konferencji «Naukowy potencjal äwiata - 2008». T.9. Techniczne nauki Budownictwo I architektura. Matematyka. Fizyka. Przemysl: Nauka i studia, 2008. С 81-84.

8. Lipanov A.M., Semakin A.N. Methods for mathematical modeling of a viscous g-flow in porous media // Materiäly IV mezinärodni vödecko-praktickä konferenc «Predni vSdecki novinky - 2008». Dil 6. Matematika. Moderni informaöni technolo gie. Fyzika. Vystavba a architektura. Praha: Publishing House «Education an' Science» s.r.o., 2008. С. 28-30.

9. Липанов А.М., Семакин А.Н. Обтекание вязким газом четырёх соприкасающихся сфер // Materiäly IV mezinärodni vedecko-praktickä konference «Nastoleni moderai vedy - 2008». Dil 8. Matematika. Modemi informa&ii technologie. Fyzika. Vystavba а architektura. Praha: Publishing House «Education and Science» s.r.o., 2008. C. 17-20.

10. Липанов A.M., Семакин A.H. Метод конечных объёмов: аналитические преобразования координат для различных типов конечных объёмов и примеры расчётов при Re=100 // Краевые задачи и математическое моделирование: сб. ст. 9-й Всероссийской научной конференции. 28-29 ноября 2008 г. В 3 т. Т. 1. Новокузнецк, 2008. С. 59-63.

11. Липанов А.М., Семакин А.Н. Обтекание вязким газом сферы в объёме с перфорированными стенками // Международная Интернет-конференция «Актуальные вопросы современной науки»: Сборник научных трудов. М.: Издательство «Спутник +», 2008. С. 82-85.

12. Семакин А.Н., Липанов А.М. Обтекание вязким газом сферы в прямоугольном объёме при наличии точки касания // Актуальные проблемы современной науки: Труды 4-го Международного форума (9-й Международной конференции молодых учёных и студентов). Естественные науки. Ч. 1-3: Математика. Математическое моделирование. Механика. Самара: Изд-во СамГТУ, 2008. С. 231-234.

Подписано в печать «15» марта2010 г. Бумага офсетная. Формат 60x84/16 Объем 1,44 П.ЮД.Л. Тираж 100 экз. Отпечатано в ИПМ УрО РАН 426067, г. Ижевск, ул. Т.Барамзиной, 34

Основные обозначения и сокращения.

ВВЕДЕНИЕ.

1. ОБЗОР МЕТОДОВ ИССЛЕДОВАНИЯ ТЕЧЕНИЙ ГАЗА И ЖИДКОСТИ ЧЕРЕЗ ПОРИСТЫЕ СРЕДЫ И ОПРЕДЕЛЕНИЕ КРУГА ЗАДАЧ.

1.1. Подходы к исследованию течений в пористой среде.

1.2. Методы численного решения уравнений гидромеханики.

1.3. Определение круга задач.

Актуальность темы:

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

Пористая среда с геометрической точки зрения имеет довольно сложное строение, и обычно течение жидкости или газа в такой среде исследуется с помощью перехода от истинных значений гидромеханических величин к фиктивным, которые «размазываются» по всей рассматриваемой среде. Далее, для этих фиктивных величин формулируется система уравнений гидромеханики. Получающаяся в результате этого математическая модель упрощает проблему исследования течения жидкости или газа через пористую среду и в некоторых частных случаях даже позволяет находить аналитические решения. Но все результаты получаются только для осреднённого по объёму течения, получить какие-либо данные о поведении жидкости или газа на уровне пор невозможно.

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

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

Среди работ, посвящённых изучению гидромеханики, можно выделить работы Лойцянского Л.Г. [1], Седова Л.И. [2], ШлихтингаГ. [3] и других. Рассмотрению движений газа и жидкости в пористых средах посвящены работы таких учёных, как Лейбензон Л.С. [4], Полубаринова-Кочина П.Я. [5], Щелка-чёв В.Н. [6], Чарный И.А. [7]. Численные методы решения задач гидромеханики приводятся в работах Андерсона Д. [8], Патанкара С. [9], РоучаП. [10], Самарского А.А. [11], Липанова A.M. [12], ФлетчераК. [13].

Объект исследования:

Течение вязкого сжимаемого газа в многосвязных областях.

Предмет исследования:

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

Цели диссертационной работы:

1. Разработка и реализация метода численного решения уравнений гидромеханики в рассматриваемых многосвязных областях.

2. Исследование течения вязкого газа в объёме с перфорированными стенками, заполненном сферическими частицами.

Для достижения поставленных целей решаются следующие задачи:

1. Создание стандартного набора конечных объёмов, на которые можно разбивать объём, заполненный сферическими частицами.

2. Разработка криволинейной системы координат для каждого выделенного конечного объёма.

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

4. Организация взаимодействия соседних конечных объёмов и передачи данных между ними.

5. Создание и тестирование программного комплекса для расчёта течения вязкого газа в объёме с перфорированными стенками, заполненном сферическими частицами.

6. Проведение расчётов для исследования течения вязкого газа в объёме со сферическими частицами.

На защиту выносятся:

1. Методика численного расчёта течения вязкого газа в объёме со сферическими частицами.

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

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

4. Результаты параметрических расчётов течения вязкого газа в объёме со сферическими частицами.

Научная новизна работы:

1. Разработана методика численного расчёта течения вязкого газа в многосвязных областях.

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

Практическая значимость:

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

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

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

Апробация работы:

Материалы диссертации апробированы на следующих конференциях: «Динамика изследования — 2008» (София, 2008 г.), «Veda: teorie a praxe - 2008» (Praha, 2008 г.), «Predni уёёескё novinky - 2008» (Praha, 2008 г.), «Naukowy po-tencjal swiata - 2008» (Przemysl, 2008 г.), «Nastolem modern! vёdy - 2008» (Praha, 2008 г.), «Актуальные проблемы науки в России» (Кузнецк, 2008 г.), «Краевые задачи и математическое моделирование» (Новокузнецк, 2008 г.), «Актуальные вопросы современной науки» (Таганрог, 2008 г.), «Актуальные проблемы современной науки» (Самара, 2008 г.). В целом диссертационная работа докладывалась на Учёном Совете Института прикладной механики УрО РАН (2009 г.).

Публикации:

Основные результаты диссертационной работы опубликованы в 12 статьях [14-25], из них 1 статья [23] - в журнале, рекомендованном ВАК для публикации результатов диссертации на соискание учёной степени доктора и кандидата наук по механике и 2 статьи [24, 25] по перечню ВАК.

Краткое содержание работы по главам:

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

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

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

В третьей главе приведены результаты расчётов течения вязкого газа для случаев, когда в объёме с перфорированными стенками расположено от 1 до 64 сфер. В начале главы представленный в предыдущем разделе численный метод тестируется на задаче обтекания сферы, расположенной в неограниченном объёме сжимаемого газа, в диапазоне чисел Рейнольдса 40-К000. Далее рассматривается обтекание вязким газом одной, двух, трёх и четырёх сфер в ограниченном объёме с одним входом и несколькими выходами, а также неплотная кубическая упаковка сфер (4-64 сферы) при различных Re. Указывается распределение давления и температуры по объёму в целом и по поверхности каждой сферы в отдельности, приводятся средние характеристики течения при входе газа в объём и выходе из него, углы отрыва потока от сфер, углы минимума давления и температуры, коэффициенты сопротивления сфер. Для каждого случая положения сфер рассматривается характер формирования вихрей, образующихся отрывных зон и движения вязкого газа в целом.

Автор выражает глубокую признатёльность научному руководителю работы академику РАН Липанову A.M. за консультации и обсуждение полученных результатов, а также другим коллегам по работе.

Основные результаты, полученные в диссертационной работе, состоят из следующих положений:

1. Разработан метод численного решения уравнений гидромеханики в многосвязных областях на примере прямоугольного объёма с перфорированными стенками, заполненного сферическими частицами. Создан набор стандартных типов конечных объёмов (КО), на которые можно разбить исходную область со сферическими частицами. Для каждого конечного объёма разработана криволинейная система координат. Для организации взаимодействия между КО предложено три метода интерполяции гидромеханических переменных (метод, основанный на вычислении определителей; метод, основанный на разложении функции по формуле Тейлора; метод, основанный на представлении функции в виде отрезка ряда Фурье по ортогональным многочленам). При проведении расчётов разностная схема, построенная исходя из дивергентной формы записи системы уравнений гидромеханики, оказывалась неустойчивой при Re >100. Устойчивость достигнута лишь при переходе к системе уравнений гидромеханики, в которой конвективные слагаемые записываются в симметричной форме.

2. На основе изложенных в главе 2 теоретических положений создан программный комплекс для проведения расчётов течения вязкого газа в объёме с перфорированными стенками, заполненном сферическими частицами. Он протестирован на задаче обтекания сферы в неограниченном объёме в диапазоне чисел Рейнольдса 40-^1000. Полученные значения коэффициента сопротивления сферы, угла отрыва потока и длины отрывной зоны удовлетворительно совпадают с экспериментальными данными и результатами расчётов других авторов.

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

4. В работе рассмотрено течение вязкого газа в объёме с перфорированными стенками, в котором располагаются одна, две, три и четыре сферы в различных положениях, а также неплотная кубическая упаковка сфер (4-64 сферы) при Re = 100 и Re = 500. Для каждого случая проанализирован характер формирования вихрей, образующихся отрывных зон и движения вязкого газа в целом для объёма в рассматриваемых условиях. В итоге получается, что, когда в объёме располагаются одна, две или три сферы, течение газа при Re = 100 становится стационарным, а при Re = 500 — нестационарным, периодическим. Когда в объёме располагаются четыре и более сфер, течение является стационарным как при Re= 100, так и при Re = 500. Во всех случаях при переходе от Re = 100 к Re = 500 наблюдается увеличение количества вихрей, они становятся более интенсивными и усложняется сама их структура.

ЗАКЛЮЧЕНИЕ

1. Лойцянский Л.Г. Механика жидкости и газа. М.: Дрофа, 2003. 840 с.

2. Седов Л.И. Механика сплошной среды. Т. 1-2. СПб.: Издательство «Лань», 2004.

3. Шлихтинг Г. Теория пограничного слоя. М.: Наука, 1974. 712 с.

4. Лейбензон Л.С. Движение природных жидкостей и газов в пористой среде. М.-Л.: Гостехиздат, 1947. 244 с.

5. Полубаринова-Кочина П.Я. Теория движения грунтовых вод. М.: Гостехиздат, 1952. 676 с.

6. Щелкачёв В.Н., Лапук Б.Б. Подземная гидравлика. М.-Л.: Гостоптехиздат, 1949. 523 с.

7. Чарный И.А. Подземная гидромеханика. М.-Л.: Гостехиздат, 1948. 198 с.

8. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. Т. 1-2. М.: Мир, 1990.

9. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.

10. Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 612 с.

11. Самарский А.А., Вабищевич П.Н. Численные методы решения задач конвекции-диффузии. М.: Едиториал УРСС, 2004. 248 с.

12. Липанов A.M., Кисаров Ю.Ф., Ключников И.Г. Численный эксперимент в классической гидромеханике турбулентных потоков. Екатеринбург: УрО РАН, 2001. 161 с.

13. Флетчер К. Вычислительные методы в динамике жидкостей. Т. 1-2. М.: Мир, 1991.

14. Липанов A.M., Семакин А.Н. Обтекание вязким газом сферы в объёме с перфорированными стенками // Международная Интернет-конференция «Актуальные вопросы современной науки»: Сборник научных трудов. М.: Издательство «Спутник +», 2008. С. 82-85.

15. Липанов A.M., Семакин А.Н. Обтекание вязким газом сферических частиц в ограниченном объёме // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. 2009. Вып. 4. С. 79-86.

16. Липанов A.M., Семакин А.Н. Обтекание трёх сфер потоком вязкого газа при Re=100 // Вестник ИжГТУ. 2008. № 4. С. 203-205.

17. Липанов A.M., Семакин А.Н. Обтекание вязким газом системы двух сфер в объёме с перфорированными стенками // Математическое моделирование. 2009. Т. 21, №7. С. 67-74.

18. Липанов A.M. Метод численного решения уравнений гидромеханики в многосвязных областях (первое сообщение) // Математическое моделирование. 2006. Т. 18, № 12. С. 3-18.

19. Коровин Н.В. Общая химия. М.: Высшая школа, 2000. 558 с.

20. Глинка Н.Л. Общая химия. М.: Интеграл-Пресс, 2004. 728 с.

21. Liu G. High Resolution Modeling of Transport in Porous Media: PhD thesis. Baton Rouge, 2002. 108 p.

22. Басниев К.С., Дмитриев Н.М., Розенберг Г.Д. Нефтегазовая гидромеханика. М.-Ижевск: Институт компьютерных исследований, 2005. 544 с.

23. Пирвердян A.M. Физика и гидравлика нефтяного пласта. М.: Недра, 1982. 192 с.

24. Маскет М. Течение однородных жидкостей в пористой среде. М.-Ижевск: НИЦ «Регулярная и хаотическая динамика», 2004. 628 с.

25. Механика насыщенных пористых сред / В.Н. Николаевский, К.С. Басниев, А.Т. Горбунов, Г.А. Зотов. М.: Недра, 1970. 339 с.

26. Азиз X., Сеттари Э. Математическое моделирование пластовых систем. М.: Недра, 1982. 407 с.

27. Ромм Е.С. Структурные модели порового пространства горных пород. Д.: Недра, 1985. 240 с.

28. Маркин B.C. О капиллярном равновесии в модели пористого тела с пересекающимися порами переменного сечения // Доклады АН СССР. 1963. Т. 151, №3. С. 620-623.

29. Ентов В.М., Фельдман А .Я., Чен-Син Э. Программное моделирование процесса капиллярного вытеснения в пористой среде // Программирование. 1975. №3. С. 67-74.

30. Balhoff М. Modeling the Flow of Non-Newtonian Fluids in Packed Beds at the Pore Scale: PhD thesis. Baton Rouge, 2005. 192 p.

31. Zhang W. Experimental and Computational Analysis of Random Cylinder Packings with Applications: PhD thesis. Baton Rouge, 2006. 162 p.

32. Guo G. The Effects of Local Hydrodynamics on Mass Transfer in Disordered Porous Media: PhD thesis. Baton Rouge, 2002. 99 p.

33. Зенкевич О. Метод конечных элементов в технике. М.: Мир, 1975. 641 с.

34. Коннор Дж., Бреббиа К. Метод конечных элементов в механике жидкости. JL: Судостроение, 1979. 264 с.

35. Норри Д., де Фриз Ж. Введение в метод конечных элементов. М.: Мир, 1981.304 с.

36. Сегерлинд JI. Применение метода конечных элементов. М.: Мир, 1979. 392 с.

37. Кудинов П.И. Численное моделирование гидродинамики и теплообмена в задачах с конвективной неустойчивостью и неединственным решением: дис. канд. физ.-мат. наук. Днепропетровск, 1999. 229 с.

38. Самарский А.А., Попов Ю.П. Разностные методы решения задач газовой динамики. М.: Едиториал УРСС, 2004. 424 с.

39. Самарский А.А., Гулин А.В. Устойчивость разностных схем. М.: Едиториал, 2005.384 с.

40. Самарский А.А. Введение в теорию разностных схем. М.: Наука, 1971. 552 с.

41. Самарский А.А. Теория разностных схем. М.: Наука, 1977. 656 с.

42. Liang G. Numerical Simulation of 3D, Complex, Turbulent Flows with Unsteady Coherent Structures: From Hydraulics to Cardiovascular Fluid Mechanics: PhD thesis. Atlanta, 2004. 220 p.

43. Steger J.L., Benek J.A. On the use of composite grid schemes in computational aerodynamics // Computer Methods in Applied Mechanics and Engineering. 1987. Vol. 64. P. 301-320.

44. Khier W., Breuer M., Durst F. Flow structure around trains under side wind conditions: a numerical study // Computers&Fluids. 2000. Vol. 29, № 2. P. 179195.

45. Takanashi D., Takemoto M. An automatic grid generation procedure for complex aircraft configurations // Computers&Fluids. 1995. Vol. 24, № 4. P. 393400.

46. Aoki T. 3D simulation for falling papers // Computer Physics Communications. 2001. Vol. 142, № 1-3. P. 326-329.

47. Tang H.S., Jones S.C., Sotiropoulos F. An overset-grid method for 3D unsteady incompressible flows // Journal of Computational Physics. 2003. Vol. 191, № 2. P. 567-600.

48. Pan H., Damodaran M. Parallel computation of viscous incompressible flows using Godunov-projection method on overlapping grids // International Journal for Numerical Methods in Fluids. 2002. Vol. 39, № 5. P. 441-463.

49. Wright J., Shyy W. A pressure-based composite grid method for Navier-Stokes equations // Journal of Computational Physics. 1993. Vol. 107, № 2. P. 225-238.

50. Ильин B.A., Позняк Э.Г. Аналитическая геометрия. М.: Наука, 1999. 224 с.

51. Умнов А.Е. Аналитическая геометрия и линейная алгебра. М.О.: Издание ЗАО «Оптимизационные системы и технологии», 2004. 368 с.

52. Мусхелишвили Н.И. Курс аналитической геометрии. М.: Высшая ттткола, 1967. 655 с.

53. Постников М.М. Аналитическая геометрия. М.: Наука, 1973. 754 с.

54. Гусак А.А., Гусак Г.М., Бричикова Е.А. Справочник по высшей математике. Минск: Тетрасистемс, 1999. 640 с.

55. Беклемишев Д.В. Курс аналитической геометрии и линейной алгебры. М.: Физико-математическая литература, 2000. 376 с.

56. Александров П.С. Лекции по аналитической геометрии, пополненные необходимыми сведениями из алгебры с приложением собрания задач, снабжённых решениями, составленного А.С. Пархоменко. М.: Наука, 1968. 912 с.

57. Бугров Я.С., Никольский С.М. Элементы линейной алгебры и аналитической геометрии. М.: Наука, 1980. 176 с.

58. Федорчук В.В. Курс аналитической геометрии и линейной алгебры. М.: Изд-во МГУ, 1990. 328 с.

59. Кадомцев С.Б. Аналитическая геометрия и линейная алгебра. М.: ФИЗ-МАТЛИТ, 2003. 160 с.

60. Лаврентьев М.А., Шабат Б.В. Методы теории функций комплексного переменного. СПб.: Издательство «Лань», 2002. 688 с.

61. Эйдерман В.Я. Основы теории функций комплексного переменного и операционного исчисления. М.: ФИЗМАТЛИТ, 2002. 256 с.

62. Камке Э. Справочник по дифференциальным уравнениям в частных производных первого порядка. М.: Наука, 1966. 260 с.

63. Комеч А.И. Практическое решение уравнений математической физики: Учеб.-метод. пособие. Механико-математический факультет МГУ, 1993. 155 с.

64. Русяк И.Г., Горохов М.М., Колосов С.М. Постановка задачи пространственного течения несжимаемой жидкости в криволинейных координатах // Интеллектуальные системы в производстве. 2006. № 1. С. 68-93.

65. Колосов С.М., Русяк И.Г. Течение вязкой теплопроводной жидкости в канале с криволинейной образующей // Вестник ИжГТУ. 2006. № 4. С. 17-21.

66. Федорченко А.Т. О методике численного исследования нестационарных дозвуковых течений вязкого газа в каналах. // Журнал вычислительной математики и математической физики. 1981. Т. 21, № 5. С. 1215-1231.

67. Федорченко А.Т. О расчётных моделях взаимодействия вихрей с проницаемой границей области дозвукового потока. // Доклады Академии наук СССР. 1983. Т. 273, № 1. С. 66-70.

68. Федорченко А.Т. О проблеме вывода вихрей через проницаемую границу расчётной области нестационарного дозвукового потока. // Журнал вычислительной математики и математической физики. 1986. Т. 26, № 1. С. 114129.

69. Федорченко А. Т. Численное исследование нестационарных дозвуковых течений вязкого газа во внезапно расширяющемся плоском канале. // Механика жидкости и газа. 1988. № 4. С. 32-41.

70. Березин И.С., Жидков Н.П. Методы вычислений. Т. 1. М.: Физматлит, 1962. 464 с.

71. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 1987. 600 с.

72. Калиткин Н.Н. Численные методы. М.: Наука, 1978. 512 с.

73. Кудрявцев JI.Д. Курс математического анализа (в двух томах): Учебник для студентов университетов и втузов. Т. 1-2. М.: Высшая школа, 1981.

74. Вержбицкий В.М. Основы численных методов. М.: Высшая школа, 2002. 840 с.

75. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М.: Наука, 1968. 801 с.

76. Демидович Б.П., Марон И.А., Шувалова Э.З. Численные методы анализа. М.: Наука, 1967. 368 с.

77. Бердышев В.И., Петрак Л.В. Аппроксимация функций, сжатие численной информации, приложения. Екатеринбург: УрО РАН, 1999. 297 с.

78. Носач В.В. Решение задач аппроксимации с помощью персональных компьютеров. М.: МИКАП, 1994. 382 с.

79. Липанов A.M., Кисаров Ю.Ф., Ключников И.Г. Численное моделирование развития вихревых структур в отрывных течениях // Математическое моделирование. 1994. Т. 6, № Ю. С. 13-23.

80. Липанов A.M., Кисаров Ю.Ф., Ключников И.Г. Математическое моделирование турбулентных течений // Математическое моделирование. 1997. Т. 9, №2. С. 113-116.

81. Липанов A.M., Кисаров Ю.Ф., Ключников И.Г. Численное моделированиевязких дозвуковых потоков при числе Рейнольдса 104 //Математическое моделирование. 1997. Т. 9, № 3. С. 3—12.

82. Липанов A.M., Ключников И.Г., Глухова Е. Ю. Решение модельных задач методами высокого порядка точности //Математическое моделирование. 1997. Т. 9, №2. С. 106-110.

83. Липанов A.M., Кисаров Ю.Ф., Ключников И.Г. Прямое численное моделирование трехмерных турбулентных течений методами высокого порядка точности // Современные проблемы внутренней баллистики РДТТ / Под ред. А.В. Алиева. Ижевск: ИПМ УрО РАН, 1996. С. 9-37.

84. Ключников И.Г. Прямое численное моделирование дозвуковых турбулентных течений газа: дис. д-ра физ.-матем. наук. Ижевск, 1999. 230 с.

85. Горохов М.М. Математическое моделирование обтекания и горения гранул твёрдого топлива в турбулентных потоках: дис. д-ра физ.-мат. наук. Ижевск, 2005. 258 с.

86. Двухфазные моно- и полидисперсные течения газа с частицами / Под ред. Л.Е. Стернина. М.: Машиностроение, 1980. 184 с.

87. Гущин В. А., Матюшин П. В. Численное моделирование пространственных отрывных течений // Применение математического моделирования для решения задач в науке и технике: Сб. трудов конференции (Ижевск, 1996). Ижевск: ИПМ УрО РАН, 1996. С. 44-61.

88. Jenson V. Viscous flow round a sphere of low Reynolds number // Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences. 1959. Vol. 249. P. 346-366.

89. Taneda S. Experimental investigation of the wake behind a sphere at low Reynolds numbers // Journal of the Physical Society of Japan. 1956. Vol. 11, № 10. P. 1104-1108.

90. Rimon Y., Cheng S.I. Numerical solution of a uniform flow over a sphere at intermediate Reynolds numbers // Physics of Fluids. 1969. Vol. 12, Jv|° 5. P. 949959.

91. Волков В.И., Мухин В.А., Накоряков B.E. Исследование структуры течения в пористой среде // Журнал прикладной химии. 1981. Т. 54, Jte 4. С. 838-842.