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

Физико-математические науки

113

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

© Обухов А.Г.* *, Сорокина Е.М.*

Тюменский государственный нефтегазовый университет, г. Тюмень

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

Ключевые слова нестационарные течения газа, полная система уравнений Навье-Стокса, краевые условия.

Для моделирования сложных течений жидкости и газа, обладающих диссипативными свойствами вязкости и теплопроводности, обычно используется полная система уравнений Навье-Стокса, которая в безразмерных переменных имеет следующий вид [1]:

V.

pt + V -Vр + р div V = 0, —

ур У

\+(V-V)V +—Vp+ -VT = g —2QxV + ^° -v(div V) + 3aV ' v ’ ур у p L 4 V ’ 4 _

— + V-VT + (у — l)TdivV =^°AT + ^°у(у-1)|^(^ -^) 2+

+ (ux — wz) 2+(vy — W)2]+ 3\_(uy + vx) 2+(uz + wx) 2+(vz + wy)2]j,

(1)

Система (1) имеет смешанный тип: первое уравнение - уравнение неразрывности (дифференциальная форма закона сохранения массы) - образует гиперболическую часть системы, так как определяет в течениях сжимаемого теплопроводного вязкого газа наличие слабого разрыва на контактной поверхности [1, 2]; второе и третье уравнения - уравнения движения и энергии (дифференциальные формы законов сохранения импульса и энергии соответственно) - составляют параболическую часть системы, так как

* Профессор кафедры «Высшая математика», доктор физико-математических наук, доцент.

* Старший преподаватель кафедры Естественнонаучных дисциплин Филиала военного учебнонаучного центра сухопутных войск «Общевойсковая академия Вооруженных Сил Российской Федерации».

114

ФУНДАМЕНТАЛЬНЫЕ И ПРИКЛАДНЫЕ ИССЛЕДОВАНИЯ

содержат вторые производные скорости и температуры по пространственной переменной. В приведенной системе (1) также учитывается влияние сил тяжести и Кориолиса [3, 4].

В системе (1): t — время; x, у, z — декартовы координаты; р— плотность газа; V = (и, v, w) — вектор скорости газа с проекциями на соответствующие декартовы оси; T — температура газа; g = (0,0, - g) — вектор ускорения силы тяжести; у = 1.4 — показатель политропы для воздуха; -2QxV = (av - bw, - аи, bu) — вектор ускорения Кориолиса, где а = 2Q sin/, b = 2Q cos/, Q = |Q|; Q —

вектор угловой скорости вращения Земли; /— широта точки O — начала декартовой системы координат xyzO, вращающейся вместе с Землей; V и div — операторы градиента и дивергенции по декартовым пространственным переменным, точкой обозначено скалярное произведение векторов.; безразмерный коэффициент вязкости р = 0.001, а безразмерный коэффициент теплопроводности к и 1.458333ро.

Расчетная область представляется в виде прямоугольного параллелепипеда с длинами сторон x0, у0 и z0 вдоль осей Ox, Oy и Oz соответственно.

Для плотности р на всех шести гранях параллелепипеда: x = 0, x = x0, у = 0, у = у0, z = 0, z = z0 — предлагается ставить «условие непрерывности» потока, то есть значения искомой функции на границу области сносятся линейной интерполяцией по нормали к данной граничной поверхности из внутренней части расчетной области.

Краевые условия для компонент вектора скорости газа предлагается брать соответствующими «условиям непротекания» для нормальной составляющей вектора скорости и «условиям симметрии» для двух других компонент вектора скорости течения. А именно:

где f— нормальная составляющая вектора скорости газа к поверхностям % = 0, %= %, а g — две другие составляющие вектора скорости газа, то есть тангенциальные по отношению к поверхностям %= 0, % = %.

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

Расчетная область заполняется трехмерной сеткой узлов пересечения трех семейств плоскостей x = xi, у = yj, z = zk, где xi = i-Ax, yj = j-Ay, zk = k-Az,

0 < i < L, 0 < j < M, 0 < k < N.

(2)

%=0,%=%

(3)

Физико-математические науки

115

Разностные шаги по трем пространственным переменным равны соответственно Ax = x0 / L, Ау = у0 / M, Az = z0 / Ж.

Пусть в начальный момент времени t = 0 во всех точках прямоугольного параллелепипеда все искомые функции заданы

U = U 0. (4)

1/=0

Затем с помощью явной разностной схемы

Un+l = Un + At ■ F (Un); U =

u

v

w

T V T

(5)

вычисляются значения всех искомых функций во всех внутренних точках прямоугольного параллелепипеда. После этого значения искомых функций определяются во всех внутренних точках каждой из шести граней: x = 0, x = x0,

У = 0, У = у°, z = 0, z = z0.

Для плотности эти значения определяются «по непрерывности», то есть находятся с помощью линейной интерполяции по значениям плотности в двух точках, ближних к рассматриваемой по нормали к данной грани (рис. 1 для случаев граней z = 0, z = z0):

Р( A) = 2p(B) -p(C).

Рис. 1

Компонента скорости вдоль оси Ox во внутренних точках граней x = 0, x = x0 полагается нулем. Во внутренних точках остальных четырех граней у = 0, у = у0, z = 0, z = z0 значения функции и определяются «по симметрии», то есть определяются из условия, что в этих точках производная по нормали к данной грани равна нулю (рис. 1 для граней z = 0, z = z0):

du

дП

3u( A) - 4u(B) + u(C)

= 0,

A

2h

116

ФУНДАМЕНТАЛЬНЫЕ И ПРИКЛАДНЫЕ ИССЛЕДОВАНИЯ

и, следовательно, получается такая расчетная формула

u( Л) = I [4u(B) - u(C)].

Компонента скорости вдоль оси Оу во внутренних точках граней у = 0, у = у0полагается нулем. Во внутренних точках остальных четырех граней х = 0, х = х0, z = 0, z = z0 значения функции v определяются «по симметрии», то есть определяются из условия, что в этих точках производная по нормали к данной грани равна нулю и, следовательно, получается такая расчетная формула (рис. 1 для граней z = 0, z = z0):

v( Л) = 1 [4v(B) - v(C)].

Компонента скорости вдоль оси Oz во внутренних точках граней z = 0, z = z0полагается нулем. Во внутренних точках остальных четырех граней х = 0, х = х0, у = 0, у = у0 значения функции w определяются «по симметрии», то есть определяются из условия, что в этих точках производная по нормали к данной грани равна нулю и, следовательно, получается расчетная формула (рис. 2 для случаев граней х = 0, х = х0):

w(A) = 1 [4w(B) - w(C)].

Для температуры T во внутренних точках всех шести граней х = 0, х = х0, у = 0, у = у0, z = 0, z = z0 значения определяются «по симметрии», (для температуры - это «условие теплоизоляции»), то есть определяются из условия, что в этих точках производная по нормали к данной грани равна нулю:

дГ

дп

3T (Л) - 4T (B) + T (C)

= 0,

Л

2h

Физико-математические науки

117

и получается такая расчетная формула (рис. 1 для граней z = 0, z = z0):

T (A) = 1 [4T (B) - T (C)].

В случае расчетов при заданных на плоскостях z = 0, z = z0 значениях температуры краевые условия считаются заданными, например:

TL = To(f > ^ )> Tz=z° = T(f ’ ^ )’

где n - номер временного слоя и тогда t = At-n.

Во внутренних точках остальных четырех граней x = 0, x = x0, у = 0, y = у0значения температуры определяются «по симметрии» (для температуры - это условие «теплоизоляции»), то есть определяются из условия, что в этих точках производная по нормали к данной грани равна нулю и, следовательно, получается такая расчетная формула (рис. 2 для случаев граней x = 0, x = x0):

T (A) = I [4T (B) - T (C)].

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

На двух гранях, образующих ребро, берутся по две точки (B1, C1), (B2, C2), расположенные на нормалях к ребру, проведенных через рассматриваемую точку A. Затем вдоль каждой из этих нормалей линейной интерполяцией определяются два своих промежуточных значения искомой вектор-функции в точке A :

U(A) = 2U(B)-U(C1); U2(A) = 2U(B2)-U(C2).

А затем значение на ребре в точке A берется как среднее арифметическое двух полученных промежуточных значений:

118

ФУНДАМЕНТАЛЬНЫЕ И ПРИКЛАДНЫЕ ИССЛЕДОВАНИЯ

щА) U( A) + U2( A)

В вершинах параллелепипеда значения вдоль каждого из трех ребер, приходящих в вершину, по двум ближайшим к вершине точкам линейной интерполяцией определяется свое промежуточное значение вдоль каждого ребра (рис. 4):

йг(Л) = 2U(B)-U(Q); и2(Л) = 2U(B2)-U(C2); U3(A) = iU(B3)-U(C3),

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

U (Л) =

U\( Л) + Ui( Л) + U3( Л) 2 '

Далее в этой работе проведен расчет одномерного нестационарного течения сжимаемого вязкого теплопроводного газа: в начальный момент времени плотность газа постоянна и равна единице, скорость газа равна нулю, то есть и = v = w = 0, а температура задана следующим соотношением T(t, х)| = 1 + 0. 1cos(^ х), и на границах отрезка 0 < x < 1 заданы условия

прилипания и теплоизоляции u(t, x)|

1=0,

dT

дх

= 0.

Результаты расчетов газодинамических параметров в некоторые фиксированные моменты времени для z0 = 0.05 представлены на последующих рисунках. На рисунках 5 и 6 приведены графики функций р(х, у) и T(x, у) на 1000 временном шаге. При увеличении времени счета значения плотности и температуры постепенно приближаются к единице, сохраняя при этом косинусоидальную зависимость от х.

На рис. 7 приведены результаты расчета первой компоненты и скорости для двух моментов времени. Видно, что за представленный промежуток времени параболическая зависимость этой компоненты скорости постепенно переходит в волнообразную зависимость от x.

Физико-математические науки

119

Рис. 5

Рис. 7

На рис. 8 изображены результаты расчетов второй компоненты v скорости для двух моментов времени. Как следует из приведенных рисунков, эта компонента скорости равна нулю во всех точках за исключением небольших пиков по углам. Этот расчетный эффект связан с неизбежным незначительным скачком второй компоненты скорости в вертикальных гранях расчетной области.

V|f=300Ar U|t=500Af

Рис. 8

120

ФУНДАМЕНТАЛЬНЫЕ И ПРИКЛАДНЫЕ ИССЛЕДОВАНИЯ

Следует отметить, что величина этих пиков имеет порядок 10-8 и никак не сказывается на вычислении плотности, температуры и первой компоненты скорости.

На рис. 9 приводятся результаты расчетов третей компоненты w скорости для трех моментов времени. При расчете этой компоненты скорости в начальные моменты времени вблизи плоскостей x = 0 и x = 1 возникают «иглы», величина которых 10"17М0"16 Приблизительно к 400 расчетному шагу по времени эти «иглы» исчезают и графические зависимости вертикальной составляющей скорости становятся гладкими поверхностями, а отличия значений скорости от нуля имеют величину порядка 10"12

Таким образом, в приведенном трехмерном расчете фактическая зависимость р, u, T только от одной пространственной переменной и фактическое равенство нулю функций v и w полностью соответствуют моделированию одномерного нестационарного течения [5, 6].

Список литературы:

1. Баутин С.П. Характеристическая задача Коши и ее приложения в газовой динамике. - Новосибирск: Наука, 2009. - 368 с.

2. Баутин С.П. Представление решений системы уравнений Навье-Сток-са в окрестности контактной характеристики // Прикладная математика и механика. - 1987. - Т 51, вып. 4. - С. 574-584.

3. Баутин С.П. Торнадо и сила Кориолиса. - Новосибирск: Наука, 2008. -96 с.

4. Баутин С.П., Обухов А.Г. Математическое моделирование разрушительных атмосферных вихрей. - Новосибирск: Наука, 2012. - 152 с.

5. Баутин С.П., Замыслов В.Е. Представление приближенных решений полной системы уравнений Навье-Стокса в одномерном случае // Вычислительные технологии. - 2012. - Т. 17, № 3. - С. 3-12.

6. Замыслов В.Е., Скачков П.П. Сравнение двух приближенных методов решения одной начально-краевой задачи газовой динамики с учетом вязкости и теплопроводности // Вестник УрГУПС. - 2012. - № 4 (16). - С. 29-38.

II

Рис. 9