modeling_2008

для системы уравнений (3.32) одномерного плоского движения однородного совершенного газа или смеси совершенных газов постоянного состава (γ = cp/cv ):

δ

x

δ , [A]

2

γ−1 2

{

− }

F

= [A]

U

=

S−1

[Λ] [S] , [Λ] = diag u, u + c, u

c ,

[S] =

cu +

γ

2

u2

c + (1

γ)u γ 1

,

−

c

+

2

u

(1

−

γ)u

γ

−

1

−c + (1 − γ)u

γ − 1

cu + γ 2

u2

−

−1

−

−

1

1

1

c2

u 2c

2

u 2c

2

S−1

− u

1

1

−

c2

2c2 + 2c

−

2c

=

2c2

.

−u2

u2

+

u

1

u2

u

−

1

+

+

2c2

4c2

2c

2(γ−1)

4c2

2c

2(γ−1)

Для расчета потоков массы компонентов смеси на границе необходимо интерполировать на нее значения парциальных плотностей ρk. Это легко сделать, дополнительно применив процедуру вышеописанного вида для массовых долей Yk; для них, как характеристических переменных системы (4.4), матричное преобразование в процедуре реконструкции не требуется.

Второй или третий порядок аппроксимации по t в методе повышенной точности достигается соотвественно при двухэтапной (5.9) – (5.10) и трехэтапной (5.11) схеме обновления решения на шаге по t.

Можно записать эти схемы через разностный пространственный оператор Lh ({Uni }); применительно к (4.5) он имеет вид

Lh ({Uin}) =	x · Fi h(Fx · F )in−		21	− (Fx · F )in+ 21 i + Fi Sin.
	1			1

В описываемом методе для расчета потоков по данным, интерполированным на границы, могут применяться как процедура точного расчета РПР, так и приближенные процедуры. Экономичная «характеристическая» процедура может быть построена на системе скалярных уравне-

ний d±I± = 0, d0I0 = 0 для инвариантов I± = u + α±p, Io = p + αoρ

системы уравнений, полученной линеаризацией системы в форме (4.7) для плоского одномерного движения газа (τw ≡ 0, qw ≡ 0, F = const)

−

1

−

0 2 0

+

1

+

(5.33)

d

u −

d

p = 0, d p − c d p = 0, d

u +

d

p = 0,

ρc

ρc

где α− = −1/ρc, αo = −c2, α+ = 1/ρc — постоянные.

171

	Un+1
	i−1		n+1		t = tn+1
			U
			i		Un+1
					i+1
xi−23	(1)		n+1	(1)
	(Fx)i−21		Ui	(Fx)i+ 21		xi+ 23
			(1)	шаг корректор
			Si
		xi−21	Sn	xi+ 21
	(Fx)n		i	(Fx)n
		1	Uin		1
	i−	2		i+	2
	Ui(1)−1		(1)	t	= t(1)	= tn+1
			Ui		Ui(1)+1
xi−23		xi−21	Ui(1)	xi+ 21		xi+ 23
				шаг предиктор

Sni

(Fx)n	1	Uin	(Fx)n	1
i−	2		i+	2

Uin	1	t = tn
−		Un
		i
		Un
		i+1

Рис. 5.7. Схема обновления решения в ячейках по методу Годунова повышенной точности (двухэтапный вариант; сравн. с рис. 5.6 на с. 168)

172

Поясним способ вычисления параметров в некоторой точке D (например, на границе ячеек) по значениям инвариантов на прибывающих в точку характеристиках системы (5.33). Для простоты возьмем случай (u − c) < 0, u > 0, (u + c) > 0, показанный на рис. 5.8:

I−D = uB + α−pB , I0D = pC + α0ρC , I+D = uA + α+pA.

Простейший способ расчета коэффициентов в инвариантах — α− = αB , α0 = αC , α+ = αA; он и предпочтителен для дозвуковых течений. По величинам инвариантов и коэффициентов вычисляются параметры в точке D (рис. 5.8):

p =	I+ − I−	, u = I	+ −	α p, ρ =	I0 − p	,
	α+ − α−			+	α0

после чего по УС определяется температура T = p/(Rρ), где постоянная R или массовые доли Y1, . . . , YK−1 для ее расчета берутся для той смеси, которая, согласно знаку скорости потока, течет через границу.

t

D

I+		I−
α+	I0	α−
	α0
A	C	B

x

Рис. 5.8. Шаблон для решения линеаризованной характеристической задачи при (u − c) 6 0, u > 0, (u + c) > 0

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

173

5.5.3. Экономичный одноэтапный метод. Инварианты линеаризованной системы уравнений одномерной газодинамики дают возможность построить весьма экономичную монотонную схему для численного решения одномерных нестационарных задач. Идея и описание такого класса схем содержится в работе [14]. Опишем применяемую нами схему такого вида.

Используется то обстоятельство, что при скалярном виде уравнений для этих инвариантов (и для массовых долей компонент Yk) процедура «монотонизированной» кусочно-параболической реконструкции решения в ячейке не требует матричного преобразования:

e

2Δxe

−

e (Δx)e

−

In(x) = In +

Iin + Iin

(x

x

) + ϕ

Iin − Iin

(x

x

)2,

i

i

i

2

i

Iin = minmod(ΔIin, b Iin),

Iin = Iin+1 − Iin,

e

Iin),

Iin = Iin

Iin−1.

e Iin = minmod( Iin, b

−

Здесь в набор зависимых переменных I включены и массовые доли K − 1 компонентов газовой смеси:

I = [I−, I0, I+, Y1, . . . YK−1]T .

Параметры смеси и газодинамические потоки на границах (рис. 5.9) определяются в данном методе с помощью характеристической процедуры, подобной описанной в подразд. 5.5.2 . В данном же случае для расчета потоков на границах i-й ячейки нужно определить параметры потока в узлах сетки (i ± 12 , n + 12 ), как в точке D на рис. 5.8, т. е. по значениям инвариантов и коэффициентов α в точках A, B, и C, определяя эти ве-

личины по кусочно-параболическим распределениям In (x) на «старом»
временном´ слое (рис. 5.8).						i
В результате получается одноэтапная схема второго порядка ап-
проксимации по			времени и		повышенного	порядка	аппроксимации
по пространственной переменной x следующего вида:
n+1	n		t		n+ 21	n+ 21	t n+ 21
Ui	= Ui +			(Fx · F )i− 21 − (Fx · F )i+ 21 +				Si	, (5.34)
			x · Fi				Fi

в которой потоки «на полушаге» (Fx)n+1/2 вычисляются описанным вы-

i+1/2

ше способом, а источниковый член может быть задан (для сохранения второго порядка по времени) как Sni +1/2 = Sni + Sni +1 /2.

174

	Un+1				t	= tn+1
	i−1		Un+1
			i			Un+1
							i+1
			Un+1			n + 1
			i
	n+	1	Sn+1		n+		1
		2		i			2
	(Fx)i−21				(Fx)i+ 21
			Sn
			i
						n
xi−23		xi−21	Uin		xi+ 21		xi+ 23
	I+				t	= tn
	I−
	I0
Рис. 5.9. Схема обновления решения в ячейках экономичным одноэтапным
методом (сравн. c. рис. 5.7). I+, I− и I0				— инварианты линеаризованных
	соотношений вдоль характеристик (5.33)

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

Изложим численный метод, использованный в расчетах пространственных турбулентных течений по технологии МКВ (с. 104) для решения уравнений, описывающих трехмерное движение газа на «надсеточном» масштабе. Основное требование к методу для расчетов в рамках МКВ — малая численная диссипация.

Для удобства описания такого метода решения уравнений вида УНС удобна условная «векторная» форма их записи

	∂U	+	∂Fx	+	∂Fy	+	∂Fz	= 0.	(5.35)
	∂t		∂x		∂y		∂z

175

Пространственные производные уравнений (5.35) аппроксимируются разностным оператором

L

n

(Fx)in−21 , j, k − (Fx)in+ 21 , j, k

h

U

i, j, k

=

+

x

(F )n

(F )n

(F )n

(F )n

1

, k −

1

1

−

1

+

y

i, j−

2

y i, j+

2

, k

+

z i, j, k−

2

z i, j, k+

2

y

z

на равномерной по пространству прямоугольной сетке (рис. 5.10) с сохранением свойства консервативности.

					zk−21			y	x
(Fz)n		(Fy)		n					z
	1			i, j	+ 21 , k
i, j, k−	2					(Fx)n		1
						i+			, j, k
								2	zk+ 21
n
Ui, j, k

(Fx)n	1	, j, k		n	n
i−	2		(Fy) i, j−21 , k		(Fz) i, j, k+ 21	xi+ 21

yj−12 xi−12

Рис. 5.10. Конечный объем — ячейка пространственной «декартовой» сетки

«Невязкая» составляющая (3.27) газодинамических потоков на границах ячеек вычисляется из решения задачи о РПР на границе ячейки с применением приближенной характеристической процедуры (с. 171);

n(−)	n(+)
для плотностей потоков Fx — по значениям Ui+ 21 , j, k	, Ui+	21 , j, k	и т. п.
Улучшенная аппроксимация по пространственным		переменным

для членов (3.27) достигается при использовании кусочно-параболичес- ких распределений параметров решения в ячейках по x, y и z. Процедура интерполяции консервативных переменных (вектора неизвестных системы уравнений однокомпонентной смеси U = [ρ, ρvx, ρvy , ρvz , ρE]T )

176

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

Так, для x-границ применяются соотношения с использовани-

ем матрицы [Sx]	и	обратной	ей		матрицы	Sx−1 преобразования
						а [Λx] — диагональная
δxFx = [Ax] δxU,	где	[Ax] = Sx−1			[Λx] [Sx],

матрица собственных значений [Ax], а также ограничительной функции (5.32) с указанными на с. 170 параметрами:

−)

n

1 + ϕ

e

1

ϕ

e

Un(

1

= Un

+ S

−1

i, j, k

Wn

+

−

x

Wn

,

4

i+

2

, j, k

i, j, k

x

x

i, j, k

4

i, j, k

(+)

n

1 + ϕ

1

−

ϕ

Un

1

= Un

S

−1

e

Wn

+

e

Wn

,

где

, j, k

i, j, k

x

i, j, k

4

i, j, k

4

x

i, j, k

i− 2

−

x

xWi,n j, k = minmod(ΔxWi,n j, k, b xWi,n j, k),

e xWi,n j, k = minmod(

xWi,n j, k, b xWi,n j, k),

xWn

k

= [Sx]n

xUn

,

xUn

= Un

−

Un

,

i, j, k

i, j, k

i, j, k

i+1, j, k

i, j, k

i, j,e

xWi,n j, k = [Sx]i,n j, k xUi,n j, k, xUi,n j, k = Ui,n j, k − Uin−1, j, k .

Диагональная матрица собственных значений:

[Λx] = diag {u, u, u, u + c, u − c} ,

матрица преобразования (χ = γ − 1):

−0

−1

−0

0

2

(u

+ v −v

−

χ

2

2 + w2)

c2

χu

χv

χw χ

x

χ2 (u2 + v2−+ w2)

−

cu χu + c

χv χw χ

χ

(u2

−χu c

−χv

−χw χ

[S ] =

+ v2 + w2) + cu

,

2

−

−

−

−

обратная матрица:

1

0

0

1

1

x

c2

u

2c2

1

u

2c2

1

− w

w

w

− u

−

v

0 1

v

v−

2c

c2

2c2

2c2

c2

0 0

v

2

2c2 + 2c

v 2

2c2

−1

−v

2

1 0

.

S

=

c2

2 +

2c2

2

2c2

2

v

w

+

+

| |

| |

u

1

| |

u

1

−

−

2c

4c

2c

2(γ−1)

4c

2c

2(γ−1)

177

Парциальные плотности компонентов ρk = ρYk на границах ячеек (для многокомпонентной смеси) вычисляются затем по интерполированным на границы величинам плотности ρ с применением также кусочно-параболической интерполяции на границы величин массовых долей компонентов Y = [Y1, . . . , YK−1]T , не требующей матричного преобразования:

−)

1 + ϕ

1

−

ϕ

Yn(

1

= Yn

+

Yn

+

Yn

,

1 + ϕ e

1

ϕ e

n(+)

i, j, k

x

i, j, k

4

i, j, k

i+

2

, j, k

4

e

x

где

= Yn

−

Yn

+

−

e

Yn

,

Y

1

i−

2

, j, k

i, j, k

4 x

i, j, k

4

x

i, j, k

e

xYi,n j, k = minmod(ΔxYi,n j, k , b xYi,n j, k ),

xYi,n j, k ).

e xYi,n j, k = minmod( xYi,n j, k , b

Плотности потоков масс компонентов, тензора потока импульса и вектора потока энергии, соответствующие эффективным потокам диффузии, вязкости и теплопроводности в системе (5.35), представляются в методе обычными центральными аппроксимациями второго порядка по пространственным координатам, и складываются с «невязкими» потоками. Так, например, плотность эффективного потока тепла через x-границу с индексами i + 12 , j, k на n-м временном´ слое задается соотношением

ý n

1

= − (

κý

n

1

Tin+1, j, k − Ti,nj, k

,

è ò

(q

)i+

x

.

)i+

2

, j, k

2

, j, k

.

По временной´ координате уравнения (5.35) интегрируются двухэтапной схемой второго порядка аппроксимации с одним «промежуточным» слоем и уточняющим пересчетом при переходе на новый слой по времени (схема «предиктор-корректор»):

				n
Ui, j, k	Ui,n+1j, k = Ui,n j, k + t · Lh						Ui,n j, k	,
	= 2 hUi, j, k + Ui, j, k + t · Lh nUi, j, ko i .
n+1	1				n+1				n+1

Описанный метод лег в основу программ для ЭВМ, использованных в расчетах пространственных течений, результаты которых представлены в разд. 7.3 и 7.4. Данный метод представляет собой явный монотонный консервативный метод — обобщение на случай 3D метода типа Годунова (подразд. 5.5.2 ) повышенной точности для равномерной сетки.

178

5.7.О методологии и методах решения уравнений гидродинамики в CFD-пакетах

Изложенный в разд. 5.6 метод для расчета пространственного течения сжимаемого рабочего тела пригоден для решения как УНС, так и уравнений, замкнутых в предположении о расчете «в режиме» МКВ. Однако данный метод, использующий простейшую «декартову» равномерную сетку ячеек (конечных объемов), например, затруднительно приспособить для корректного решения реальных задач детального моделирования РП ТД, в которых очертания расчетной области произвольные и зачастую переменные во времени.

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

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

методологии:

•рассмотрение исходной системы уравнений модели течения среды как набора обобщенных интегральных уравнений переноса вида (3.28) — что обеспечивает необходимую гибкость независимого выбора частных моделей и набора совместно решаемых уравнений переноса при использовании:

–предположений о сжимаемом, слабо- и несжимаемом течении;

–уравнений состояния — от ρ = const до задаваемых пользователем;

– моделей	одно-	и	многокомпонентной,	реагирующей
или инертной смеси,			в т. ч. модели многофазного течения
и излучения;
– моделей	влияния	эффектов турбулентности		на перенос

(в т. ч. — в пристенной зоне), химические реакции, межфазный обмен и излучение;

179

•дискретизация уравнений по пространству (не выше 2-го порядка аппроксимации) — в контексте МКО, рассчитанная на применение сеток с ячейками в форме многогранников (рис. 5.11) — для обеспечения корректности и удобства:

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

–реализации произвольного измельчения ячеек сетки для выявления детальной местной структуры решения, в первую очередь — в приcтенных зонах с повышенными градиентами характеристик потока;

•дискретизация уравнений по времени — через применение неявных методов, позволяющих:

–существенно повысить допустимую величину временного´ шага t (невысокую для явного аналога того же метода — из-за ограничения на t по условию CF L для мелких ячеек, наличия ячеек сетки «плохой» формы и др. факторов, снижающих устойчивость разностной схемы), а также в некоторой степени ослабить вредное влияние «жесткости» конкретных частных моделей и сочетания конкретных УЧП в решаемой системе;

–в расчетах, нацеленных на получение стационарной

(stationary) картины течения — использовать шаги по времени, характеризуемые CF L 1 и достигать сходимости к стационарному решению за весьма умеренное число

итераций (iterations);

–в расчетах нестационарных (transient) течений — использовать желаемые (адекватные для ячеек в «ядре» потока) значения CF L ≈ 1 без потери устойчивости от вышеупомянутых факторов, и с неплохими шансами (см. ниже) на получение физически достоверной картины развития решения по времени (time-resolved solution).

•решение систем линейных уравнений, вытекающих из соотношений неявных методов для уравнений переноса — численное, итера-

180