Лекции - Математические модели и численные методы в динамике жидкости и газа

Подождите немного. Документ загружается.

4.3. КГД уравнения в цилиндрической системе координат 71

ϕϕ

= r

+ 2

∂u

∂ϕ

−

− ζ

div ~u

+ τ ρ

∂u

∂r

+ u

∂u

∂ϕ

+ u

+ ru

∂u

∂z

∂p

∂ϕ

+ τ

∂p

∂r

∂p

∂ϕ

+ u

∂p

∂z

+ γp div ~u

ϕz

= rΠ

= η

∂u

∂z

∂u

∂ϕ

+ τρu

∂u

∂r

∂u

∂ϕ

+ u

∂u

∂z

∂p

∂z

= Π

= η

∂u

∂z

∂u

∂r

+ τρu

∂u

∂r

∂u

∂ϕ

−

+ u

∂u

∂z

∂p

∂r

zϕ

= rΠ

= η

∂u

∂z

∂u

∂ϕ

+ τρ

∂u

∂r

+ u

∂u

∂ϕ

+ u

+ ru

∂u

∂z

∂p

∂ϕ

= Π

= η

∂u

∂z

−

div ~u

+ τ ρu

∂u

∂r

∂u

∂ϕ

+ u

∂u

∂z

∂p

∂z

+ τ

∂p

∂r

∂p

∂ϕ

+ u

∂p

∂z

+ γp div ~u

В выражениях для компонент тензора вязких напряжений слагаемые

с множителем η соответствуют тензору Навье–Стокса, а слагаемые с мно-

жителем τ являются КГД добавками.

72 Глава 4. КГД уравнения и системы координат

И, наконец, выпишем компоненты вектора теплого потока ~q.

= q

= −η

γ − 1

P r

∂

∂r

−

− τ ρu

γ − 1

∂

∂r

∂

∂ϕ

+ u

∂

∂z

¶¸

−

− τ ρu

∂

∂r

∂

∂ϕ

+ u

∂

∂z

¶¸

= rq

= −η

γ − 1

P r

∂

∂ϕ

−

− τ ρu

−

∂

∂r

∂

∂ϕ

+ u

∂

∂z

¶¸

−

− τ ρu

∂

∂r

∂

∂ϕ

+ u

∂

∂z

¶¸

= q

= −η

γ − 1

P r

∂

∂z

−

− τ ρu

γ − 1

∂

∂r

∂

∂ϕ

+ u

∂

∂z

¶¸

−

− τ ρu

∂

∂r

∂

∂ϕ

+ u

∂

∂z

¶¸

Отметим, что по аналогии с выписанными выше компонентами

и Π

в выражениях для ~q слагаемые с множителем η соответствуют вектору

теплового потока Навье–Стокса, слагаемые с множителем τ являются

КГД добавками.

Глава 5

Численные алгоритмы решения

нестационарных задач газовой динамики

Данная, самая объемная глава, посвящена построению эффективных

конечно-разностных алгоритмов решения КГД уравнений для численно-

го моделирования нестационарных газодинамических течений с исполь-

зованием ортогональных сеток.

Разностные аппроксимации строятся в потоковой форме непосред-

ственно для векторов плотности потока массы

, вектора теплового

потока ~q и тензора вязких напряжений Π. Это соответствует записи урав-

нений газовой динамики в виде законов сохранения и делает алгоритм

достаточно компактным и экономичным. Устойчивость разностного ал-

горитма обеспечивается введением искусственной диссипации, вид кото-

рой определяется КГД добавками.

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

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

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

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

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

мерных задач об осесимметричном течении в окрестности цилиндра и

плоском течении в канале с прямой и обратной ступенькой.

В основу этого раздела положены материалы работ [7], [29], [42], [44],

[45].

5.1 Система уравнений для плоских двумерных те-

чений

Запишем систему КГД уравнений в декартовых координатах для дву-

мерного случая в виде, удобном для численной реализации:

∂ρ

∂t

∂j

∂x

∂j

∂y

= 0, (5.1.1)

74Глава 5. Численные алгоритмы решения нестационарных задач газовой динамики

∂(ρu

)

∂t

∂(j

)

∂x

∂(j

)

∂y

∂p

∂x

∂Π

∂x

∂Π

∂y

, (5.1.2)

∂(ρu

)

∂t

∂(j

)

∂x

∂(j

)

∂y

∂p

∂y

∂Π

∂x

∂Π

∂y

, (5.1.3)

∂E

∂t

∂(j

∂x

∂(j

∂y

∂q

∂x

∂q

∂y

∂

∂x

(Π

+ Π

) +

∂

∂y

(Π

+ Π

) , (5.1.4)

где E — полная энергия единицы объема и H — полная удельная энталь-

пия:

E = ρ

+ u

γ − 1

, H =

(E + p)

, p = ρRT. (5.1.5)

Компоненты вектора плотности потока массы

вычисляются по

формулам:

= ρ(u

− w

), j

= ρ(u

− w

), (5.1.6)

где

∂(ρu

)

∂x

∂(ρu

)

∂y

∂p

∂x

∂(ρu

)

∂x

∂(ρu

)

∂y

∂p

∂y

. (5.1.7)

Компоненты тензора вязких напряжений Π определяются с помощью

выражений:

= Π

+ u

∗

+ R

∗

, Π

= 2η

∂u

∂x

−

η div ~u,

= Π

+ u

∗

, Π

= Π

= η

∂u

∂x

∂u

∂y

= Π

+ u

∗

, (5.1.8)

= Π

+ u

∗

+ R

∗

, Π

= 2η

∂u

∂y

−

η div ~u,

5.1. Система уравнений для плоских двумерных течений 75

здесь Π

, Π

— компоненты навье-стоксовского тензора

вязких напряжений.

Величины w

∗

, w

∗

, R

∗

вычисляются по формулам:

∗

= τ

ρu

∂u

∂x

+ ρu

∂u

∂y

∂p

∂x

∗

= τ

ρu

∂u

∂x

+ ρu

∂u

∂y

∂p

∂y

, (5.1.9)

∗

= τ

∂p

∂x

+ u

∂p

∂y

+ γpdiv~u

а дивергенция вектора скорости div ~u равна:

div ~u =

∂u

∂x

∂u

∂y

. (5.1.10)

Компоненты вектора теплового потока ~q находятся по формулам:

= q

− u

, q

= q

− u

, (5.1.11)

= τ ρ

γ − 1

∂

∂x

γ − 1

∂

∂y

+ pu

∂

∂x

+ pu

∂

∂y

¶¸

где навье-стоксовские слагаемые q

и q

вычисляются как:

= −æ

∂T

∂x

, q

= −æ

∂T

∂y

. (5.1.12)

Зависимость коэффициента динамической вязкости η от температуры

выберем в следующем виде

η = η

∞

, (5.1.13)

где η

∞

= η(T

∞

) — известное значение η при температуре T

∞

Коэффициент теплопроводности æ и релаксационный параметр τ свя-

заны с коэффициентом динамической вязкости η соотношениями:

æ =

γR

(γ − 1)P r

η, τ =

p Sc

η, (5.1.14)

где Pr — число Прандтля, Sc — число Шмидта.

76Глава 5. Численные алгоритмы решения нестационарных задач газовой динамики

5.2 Система уравнений в цилиндрической геометрии

В цилиндрических координатах в двумерном случае система КГД урав-

нений имеет вид:

∂ρ

∂t

∂(rj

)

∂r

∂j

∂z

= 0, (5.2.1)

∂(ρu

)

∂t

∂(rj

)

∂r

∂(j

)

∂z

∂p

∂r

∂(rΠ

)

∂r

∂Π

∂z

−

φφ

, (5.2.2)

∂(ρu

)

∂t

∂(rj

)

∂r

∂(j

)

∂z

∂p

∂z

∂(rΠ

)

∂r

∂Π

∂z

, (5.2.3)

∂E

∂t

∂(rj

∂r

∂(j

∂z

∂(rq

)

∂r

∂q

∂z

∂

∂r

[r (Π

+ Π

)] +

∂

∂z

(Π

+ Π

) , (5.2.4)

где E — полная энергия единицы объема и H — полная удельная энталь-

пия:

E = ρ

+ u

γ − 1

, H =

(E + p)

, p = ρRT. (5.2.5)

Здесь u

и u

— проекции вектора скорости ~u на оси r и z. Влияние

внешних сил не учитывается.

Компоненты вектора плотности потока массы

вычисляются по

формулам:

= ρ(u

− w

), j

= ρ(u

− w

), (5.2.6)

где

∂

∂r

(rρu

) +

∂

∂z

(ρu

) +

∂p

∂r

∂

∂r

(rρu

) +

∂

∂z

(ρu

) +

∂p

∂z

Компоненты тензора вязких напряжений Π выпишем в виде, удобном

5.3. Граничные условия 77

для программной реализации:

= Π

+ u

∗

+ R

∗

, Π

= 2η

∂u

∂r

−

η div ~u,

= Π

+ u

∗

, Π

= Π

= η

∂u

∂z

∂u

∂r

= Π

+ u

∗

, (5.2.7)

φφ

= Π

φφ

+ R

∗

, Π

φφ

= 2η

−

η div ~u,

= Π

+ u

∗

+ R

∗

, Π

= 2η

∂u

∂z

−

η div ~u,

здесь Π

, Π

φφ

, Π

— компоненты навье-стоксовского тен-

зора вязких напряжений.

Величины w

∗

, w

∗

, R

∗

вычисляются по формулам:

∗

= τ

ρu

∂u

∂r

+ ρu

∂u

∂z

∂p

∂r

∗

= τ

ρu

∂u

∂r

+ ρu

∂u

∂z

∂p

∂z

∗

= τ

∂p

∂r

+ u

∂p

∂z

+ γ p div ~u

а дивергенция вектора скорости div ~u равна:

div ~u =

∂(ru

)

∂r

∂u

∂z

Компоненты вектора теплового потока ~q находятся по формулам:

= q

− u

, q

= q

− u

(5.2.8)

= τ ρ

γ − 1

∂

∂r

γ − 1

∂

∂z

+ pu

∂

∂r

+ pu

∂

∂z

¶¸

где навье-стоксовские слагаемые q

и q

задаются формулами:

= −æ

∂T

∂r

, q

= −æ

∂T

∂z

. (5.2.9)

5.3 Граничные условия

Приведенные системы уравнений следует дополнить начальными и гра-

ничными условиями. Постановка этих условий определяется конкретной

78Глава 5. Численные алгоритмы решения нестационарных задач газовой динамики

решаемой задачей. Приведем пример граничных условий для задачи о

течении газа в окрестности цилиндра в осесимметричной геометрии (см.

например рис. 5.5).

Профиль течения на входной границе, расположенной в плоскости

z = z

(вертикальной стенке), можно задать как

ρ = ρ

∞

, u

= u

z∞

, u

= 0, p = p

∞

. (5.3.1)

На оси симметрии, совпадающей с осью z, задаются условия симмет-

рии:

∂ρ

∂r

= 0,

∂u

∂r

= 0, u

= 0,

∂p

∂r

= 0. (5.3.2)

На свободных границах, где предполагается, что газ вытекает из рас-

сматриваемой области, ставятся так называемые "мягкие" граничные

условия, или условия "сноса". Здесь предполагается равенство нулю нор-

мальных производных плотности, давления и компонент скорости. На-

пример, если граница находится в плоскости z = z

, то такие условия

запишутся как:

∂ρ

∂z

= 0,

∂u

∂z

= 0,

∂u

∂z

= 0,

∂p

∂z

= 0. (5.3.3)

Пусть одной из границ является непроницаемая твердая стенка. Раз-

ложим скорость ~u вблизи границы на две компоненты: нормальную u

тангенциальную u

. На твердой стенке можно поставить условия непро-

текания u

= 0 и прилипания u

= 0 (или скольжения ∂u

/∂n = 0) для

скорости, и некоторое условие для температуры.

В КГД системе уравнений вектор плотности потока массы вычисля-

ется как

QGD

= ρ~u − ρ ~w = ρ~u − τ

div(ρ~u ⊗ ~u) +

∇p

. (5.3.4)

Для того чтобы поток массы, вычисляемый по формуле (5.3.4), был ра-

вен нулю на границе, условия непротекания для скорости необходимо

дополнить условием для давления вида ∂p/∂n = 0.

Пусть у нас имеется цилиндр, расположенный вдоль оси z. На поверх-

ности r = r

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

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

будем полагать, что поверхность непроницаема. Тогда граничные усло-

вия запишутся следующим образом

= 0, u

= 0, T = T

∂p

∂r

= 0. (5.3.5)

5.4. Безразмерный вид уравнений 79

Выпишем выражения для тепловой потока и силы трения на границе,

для которой выполняется условия непротекания для скорости u

= 0.

Пусть граница расположена в плоскости z = z

(торец цилиндра). Си-

ла трения на этой границе определяется компонентой тензора вязких

напряжений Π

, а тепловой поток равен q

= q

− u

= τ ρ

γ − 1

∂

∂r

γ − 1

∂

∂z

+ pu

∂

∂r

+ pu

∂

∂z

¶¸

= Π

+ τu

ρu

∂u

∂z

+ ρu

∂u

∂r

∂p

∂r

. (5.3.6)

С учетом условия непротекания на стенке (u

= 0), сразу получаем

выражение для теплового потока и коэффициента трения на стенке в

виде:

= q

= −æ

∂T

∂z

= Π

= η

∂u

∂z

. (5.3.7)

То есть для КГД уравнений выражения для теплового потока и силы

трения на твердой стенке (5.3.7) совпадают с соответствующими величи-

нами для уравнений Навье–Стокса.

5.4 Безразмерный вид уравнений

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

представить их в безразмерном виде. Это позволяет, во-первых, опериро-

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

выделить численные коэффициенты в уравнениях, от которых зависит

решение задачи — выделить так называемые параметры подобия.

Выберем в качестве основных размерных параметров характерный

линейные размер R

(например, радиус цилиндра), плотность набегаю-

щего потока ρ

∞

и скорость звука в набегающем невозмущенном потоке

∞

Запишем соотношения между размерными и безразмерными величи-

80Глава 5. Численные алгоритмы решения нестационарных задач газовой динамики

нами (знак "∼ " над переменной относится к безразмерным величинам):

ρ = ˜ρ ρ

∞

, u = ˜u c

∞

, c = ˜c c

∞

, p = ˜p ρ

∞

x = ˜x R

, t =

∞

, (5.4.1)

T =

Rρ

˜p ρ

∞

R ˜ρ ρ

∞

˜pγ

˜ρ

Rγ

∞

γR

Определим числа Маха и Рейнольдса:

∞

, Re

∞

. (5.4.2)

Уравнение полной энергии примет вид:

E =

ρu

(γ − 1)

, ⇒

E =

˜ρ˜u

˜p

(γ − 1)

. (5.4.3)

При обезразмеривании скорость звука преобразуется как

c =

γRT , ⇒ ˜c =

T . (5.4.4)

Уравнение состояния:

p = ρRT, ˜p ρ

∞

= ˜ρ ρ

∞

γR

⇒ ˜p =

˜ρ

. (5.4.5)

Таким образом уравнения связи (5.4.4) и (5.4.5) после обезразмерива-

ния изменили свой вид.

Подставляя соотношения (5.4.1) в уравнения КГД системы, можно

убедиться, что приведение к безразмерному виду не меняет вида исход-

ных уравнений. При этом безразмерные коэффициенты вязкости, тепло-

проводности и релаксационный параметр вычисляются как

˜µ =

, (5.4.6)

˜æ =

(γ − 1)P r

. (5.4.7)

˜τ =

ReSc

˜p

. (5.4.8)

Далее, используя безразмерные КГД уравнения и уравнения связи,

знак "∼" ("тильда") будем опускать.