Расчет пожарных рисков / Ocenka i raschet pozharnogo riska (NIIPPB) 2012

5.3.5. Дискретизация обобщенного уравнения на криволинейных сетках

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

стве в декартовой системе координат, в которой координатные линии и по-

верхности пересекаются под прямым углом. Однако при моделировании те-

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

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

писать исходные уравнения в обобщенных криволинейных координатах, со-

гласованных с границами расчетной области. В этом случае удается более точно задать граничные условия. Для обобщенного уравнения переход к криволинейным координатам (ξ, η, ζ) дает:

∂ ( ρФ)

× J +

∂

(ρUФ) +

∂

( ρVФ)

+

∂

( ρWФ) =

∂ξ

∂η

∂ζ

∂ t

∂

G

Φ

∂ Ф

+ q12

∂ Ф

+ q13

∂ Ф

q11

∂ξ

∂η

∂ζ

J

∂ξ

(5.3.5.1)

∂

G

∂

∂

∂ Ф

+

Φ

Ф

+ q22

Ф

+ q23

q21

∂η

J

∂ξ

∂η

∂ζ

+

∂

G

Φ

∂

Ф

+ q32

∂

Ф

+ q33

∂ Ф

+ JSФ ,

q31

∂ζ

J

∂ξ

∂η

∂ζ

где J - якобиан преобразования;

U, V, W - контравариантные составляющие скорости;

qij - функции метрических коэффициентов;

J = xξ yη zζ + xη yζ zξ + xζ yξ zη − xξ yζ zη − xη yξ zζ − xζ yη zξ ,

U = ( yη zζ − yζ zη )u + (zη xζ − zζ xη )v + (xη yζ − xζ yη )w,

V= ( yζ zξ − yξ zζ )u + (zζ xξ − zξ xζ )v + (xζ yξ − xξ yζ )w,

W= ( yξ zη − yη zξ )u + (zξ xη − zη xξ )v + (xξ yη − xη yξ )w,

xξ , x η , xζ , y ξ , y η , yζ , zξ , zη , zζ

q = ( y z − y z )2				+ (z x − z x )2		+ (x y − x y )2			,
11		η ζ	ζ η	η ζ	ζ η	η ζ	ζ η
q	22	= ( y z − y z )2		+ (z x − z x )2		+ (x y − x y )2			,
		ζ ξ	ξ ζ	ζ ξ	ξ ζ	ζ ξ	ξ ζ
q = ( y z − y z )2				+ (z x − z x )2		+ (x y − x y )2		,
33		ξ η	η ξ	ξ η	η ξ	ξ η	η ξ

q12	= ( yη zζ − yζ zη )( yζ zξ − yξ zζ ) + (zη xζ − zζ xη )(zζ xξ − zξ xζ ) +
(xη yζ − xζ yη )(xζ yξ − xξ yζ ),
q13	= ( yη zζ − yζ zη )( yξ zη − yη zξ ) + (zη xζ − zζ xη )(zξ xη − zη xξ ) +
(xη yζ − xζ yη )(xξ yη − xη yξ ),			(5.3.5.2)
q23 = ( yζ zξ − yξ zζ )( yξ zη − yη zξ ) + (zζ xξ − zξ xζ )(zξ xη − zη xξ ) +
(xζ yξ − xξ yζ )(xξ yη − xη yξ ),
q21 = q12		,
q31	= q13	,
q32	= q23 .

В преобразованном математическом пространстве для упрощения за-

писи все расстояния между координатными узлами (ξ, η, ζ) выбираются еди-

ничными. Производные определяются чис-

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

Построение конечно-разностного аналога исходных дифференциаль-

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

ванной с границами расчетной области системе координат (ξ, η, ζ) обобщен-

ное уравнение (5.3.5.1) интегрируется при помощи метода контрольного объ-

ема.

164

∂ ( ρФ)

× J + ( ρUФ) - (ρUФ) + ( ρVФ) - ( ρVФ) + (ρWФ) - (ρWФ) =

∂ t

R

L

D

F

T

B

G

∂ Ф

∂ Ф

∂ Ф

G

∂ Ф

∂ Ф

∂ Ф

Φ

q11

∂ξ

+ q12

∂η

+ q13

∂ζ

-

Φ

q11

∂ξ

+ q12

∂η

+ q13

∂ζ

J

R

J

L

(5.3.5.3)

G

Φ

∂ Ф

∂ Ф

∂ Ф

G

Φ

∂ Ф

∂ Ф

∂ Ф

+

q21

∂ξ

+ q22

∂η

+ q23

∂ζ

-

q21

∂ξ

+ q22

∂η

+ q23

∂ζ

J

J

D

F

G

Φ

∂ Ф

∂ Ф

∂ Ф

G

Φ

∂ Ф

∂ Ф

∂ Ф

+

q31

∂ξ

+ q32

∂η

+ q33

-

q31

∂ξ

+ q32

∂η

+ q33

∂ζ

+ JSФ ,

J

J

∂ζ

T

B

Члены, содержащие смешанные производные вычисляются при помо-

щи значений Ф полученных на предыдущей итерации (k).

∂ ( ρФ)

× J + (ρUФ)

- (ρUФ)

+ (ρVФ)

- (ρVФ)

+ (

ρWФ)

- (ρWФ)

k +1

=

∂ t

R

L

D

F

T

B

G

Φ q

∂ Ф

-

G

Φ

∂ Ф

k +1

q

∂ξ

+

J

11

11

∂ξ R

J

L

G

Φ q

∂ Ф

-

G

Φ q

∂ Ф

k +1

22

∂η

+

J

22

D

J

∂η F

G

Φ q

∂ Ф

-

G

Φ q

∂ Ф

k +1

∂ζ

∂ζ

+

J

33

33

T

J

B

G

Φ

∂ Ф

∂

Ф

G

Φ

∂ Ф

∂ Ф

k

q12

∂η

+ q13

∂ζ

-

q12

∂η

+ q13

∂ζ

+

J

J

R

L

G

Φ

∂ Ф

∂

Ф

G

Φ

∂ Ф

∂ Ф

k

q21

∂ξ

+ q23

∂ζ

-

q21

∂ξ

+ q23

∂ζ

+

J

J

D

F

G

∂ Ф

∂

Ф

G

Φ

∂ Ф

∂ Ф

k

k

Φ

q31

∂ξ

+ q32

-

q31

∂ξ

+ q32

∂η

+ J ×(SФP )

(5.3.5.4)

J

J

∂η

T

B

Аппроксимация смешанных производных осуществляется при помощи цен-

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

		∂ Ф	= ( q12 )L ×	1	(Фi , j +1,k	+ Фi −1, j +1,k	-Фi, j −1,k - Фi−1, j −1,k )	(5.3.5.5)
	q12
				4
		∂η L

5.3.6. Дискретизация нестационарного члена

Дискретизация по времени осуществляется неявным способом одной из двух представленных схем аппроксимации временной производной.

Неявная схема Эйлера первого порядка:

(ρF)nP+1 - (ρF)nP	h	n+1		(5.3.6.1)
	+ L (F)		= 0
τ

Неявная схема Пейре второго порядка:

3×(ρF)nP+1 - 4 ×(ρF)nP + (ρF)nP−1	h	n+1		(5.3.6.2)
	+ L (F)		= 0
2τ

где τ -шаг по времени, Lh (Φ) -разностный оператор полученный после дискре-

тизации уравнения (3.1) по пространству.

5.4. Связь полей скорости и давления (SIMPLE алгоритмы)

При решении уравнений Навье-Стокса в естественных переменных для несжимаемой жидкости возникают трудности в связи со сложностью интер-

претации взаимодействия давления и составляющих скорости. Это обуслов-

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

ных уравнениях не выражается явным образом. Подход с расщеплением ис-

ходной задачи, реализованный в процедуре SIMPLE, предложенной Патанка-

ром и Сполдингом, позволяет разрешить эту проблему. Согласно этой проце-

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

сти выводится уравнение для поправки давления. Используя решение урав-

нения для поправки давления, производится коррекция поля скорости и дав-

ления. На практике используется модификация описанного выше алгоритма - SIMPLEC-процедура, применение которой приводит к значительному повы-

шению эффективности расчетов.

Шаг 1:

Запишем дискретный аналог уравнения количества движения (5.1.1.2)

A		v p	= ∑A		vsp	+ S
	pu			spu			u
						(5.4.1)

и применим операцию нижней релаксации:

166

1

A

v p* = ∑A

* + S

+

1− αu

A

v p k

pu

spu vsp

u

pu

αu

αu

,

(5.4.2)

здесь v p k , v p*

- скорость на k-ой итерации и промежуточное поле скоро-

сти. На первом шаге из уравнений (5.4.2) находится промежуточное поле скорости.

Шаг 2:

Согласно SIMPLE-C процедуре из дискретных аналогов уравнений ко-

личества движения (5.4.2) и неразрывности (5.1.1.1) выводится уравнение для поправки давления:

Ñ(τÑp¢) =	¶ρ + Ñ(ρ v )	,	(5.4.3)
	¶t

где локальный параметр τ для SIMPLE-C процедуры:

τ =			ρ J
	A
		p	- ∑Asp
	α					(5.4.4)
		u
На втором шаге находится поправка давления, после чего определяется
давление, часто с введением нижней релаксации:
pn+1 = pn + α p p¢					.	(5.4.5)

Шаг 3:

После определения поправки давления корректируется скорость, полу-

ченная на первом этапе:

	ρ vk +1 - ρ v*	= -Ñp¢
	τ
		(5.4.6)

Скорректированная скорость принимается за скорость на k+1-ой ите-

рации.

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

Условия на входе потока

На входе используется либо фиксированное значение скорости, либо массовый расход. При этом, если задан расход, то компоненты скорости u, v,

w в каждом входе рассчитываются по заданному массовому расходу, углу наклона оси входного потока к осям координат и по доли потока, заданного для данного входного «окна». Например, компонента u на входе определяет-

ся следующим образом:

u =	G × cosα
in	ρ A	(5.5.1)
где G – массовый расход жидкости на входе, ρ – суммарная плотность, А –

площадь сечения входного окна, α – угол между осью OX и осью входного потока.

Кинетическая энергия турбулентности k и скорость диссипации турбу-

лентности ε на входе оцениваются по соотношениям вида:

kin = kin0 (uin2 + vin2 + win2 ),
ε = ε0	k3/ 2 .	(5.5.2)
in	in

где константы koin, εoin подбираются эмпирически.

Обычно для развитых турбулентных потоков на входе уровень турбу-

лентных пульсации составляет порядка 1%, поэтому koin=0,01, а εoin =Cμ/L,

где L – характерный размер энергосодержащих вихрей.

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

ставу газов или жидкости.

Граничные условия для Ym задаются в виде массовых долей компонен-

тов на входе.

Условия на выходе потока

На выходной границе при решении уравнений для u, v, w, k, ε, h, Ym

ставятся условия отсутствия градиентов (так называемые «мягкие усло-

вия»):

∂ u =	∂ v =	∂ w =	∂ h =	∂ k =	∂ ε =	∂ Ym = 0
∂ n	∂ n	∂ n	∂ n	∂ n	∂ n	∂ n	,	(5.5.3)

где n – вектор внешней нормали к расчетной области.

168

Условие симметрии (скольжение)

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

ной по нормали к плоскости симметрии всех скалярных величин Φ и тан-

генциальной составляющей скорости u|| , нормальной к плоскости состав-

ляющей скорости u .

∂Φ = 0		∂u||	= 0		u = 0
∂n	, ∂n			,		(5.5.4)

где n – вектор нормали к плоскости симметрии.

Твердая стенка

Нормальную и тангенциальную компоненты скорости на стенках пола-

гаем равными нулю, что моделирует соответственно непротекание и прили-

пание.

,

(5.5.5)

Для определения турбулентных характеристик вблизи стенки исполь-

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

Для массовых долей задаются условия отсутствия диффузионных по-

токов через стенку (конвективные потоки отсутствуют в силу (5.5.5))

∂Ym = 0
∂n	(5.5.6)

где n – вектор нормали к стенки.

Теплоотдача на стенках

Граничные условия для уравнения на энтальпию можно задать двумя способами:

1) Задание постоянной температуры на стенке

Tw = Tset

(5.5.7)

2)	Задание теплового потока и теплоотдачи на стенке
	λ ∂T = α ×(T		- T	) + q	(5.5.8)
	¶n	w	∞	set
где Tw –	температура на стенке, Tset – заданная температура, T∞ –				темпера-
тура внешней среды, qset – заданный тепловой поток через стенку,					α – ко-
эффициент теплоотдачи.

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

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

∑a f φS + Sφ −φP aP = 0 f

aPi, j ,kφi, j ,k = aLi, j ,kφi−1, j ,k + aRi, j ,kφi+1, j ,k + aFi, j ,kφi, j −1,k + aDi, j,kφi, j +1,k +
+aBi, j ,kφi, j ,k −1 + aTi, j ,kφi, j ,k +1 + Sφ i, j ,k	(5.6.1)

Для решения систем линейных алгебраических уравнений в пакете реа-

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

вующей литературе:

·полинейный,

·переменных направлений,

·неполной факторизации Булеева (Булев Н.И., 1989),

·метод сопряженных невязок с факторизацией по Булееву (Булев Н.И., 1989)

·многосеточный.

170

5.7 Радиационный перенос энергии.

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

ломассопереноса каждый элементарный объем среды находится в непосред-

ственном взаимодействии со всеми другими элементарными объемами и ре-

шение интегро-дифференциальных уравнений, описывающих это явление,

очень трудоемко. Вычислительные затраты на решение уравнения радиаци-

онного теплопереноса (УРТ) (5.7.1) может существенно превосходить затра-

ты на решение всех остальных процессов.

dI (r , sˆ )

= −β (rR )I (rR, sˆ )

+ κ (rR)Ib (rR ) +

σ (r )

I (rR, sˆ′)Φ(sˆ′, sˆ )dΩ′

(5.7.1)

R

R

4π

∫

ds

4π

Граничные условия на УРТ для случая диффузионного отражения

R

R

R

ρ(r )

R ′

′

′

I (r , s) = ε (r )I

(r )+

R

(r , s )

s ∙ n

dΩ

b

π

∫

I

(5.7.2)

ˆ

sˆ ′∙nˆ <0

ˆ

ˆ

ˆ

Радиационный источниковый член в уравнении энергии

Ñ · q = k(Eb

(r )- E(r ))

(5.7.3)

R

R

E(rR) = ∫I (rR, sR)dΩ

4π

Eb (rR ) = 4 ×σ 0T 4

где I – интенсивность (Вт/(ср*м2)); Ib – интенсивность абсолютно черного тела, σ – коэффициент рассеивания, м-1; к – коэффициент поглощения, м-1; β

– коэффициент затухания (β = к + σ), м-1; W - телесный угол, рад; r – радиус-

вектор, м; s – угловое направление, м; Ф – функция рассеивания.

Коэффициенты поглощения газа к вычисляются по модели суммы се-

рых газов (WSGG модели), коэффициенты поглощения рассчитываются че-

рез оптическую плотность дыма.

Для задачи радиационного обмена в разработан ряд подходов. Это ме-

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

ним относятся различные варианты метода Монте-Карло и методы дискрет-

ных направлений.

Ряд других методов основан на представлении исходного интегро-

дифференциального уравнения переноса излучения в форме дифференциаль-

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

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

Основными достоинствами метода являются малые требования к вычисли-

тельным ресурсам и его легкая совместимость с методами расчета аэродина-

мики и теплопереноса. При расположении очага пожара в помещении не-

больших объемов целесообразно использовать диффузионное приближение,

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

объемный метод (КОМ) и дискретно-ординатный метод (ДОМ)

Конечно-объемный метод

КОМ – один методов из дискретных направлений. Данные методы не имеют ограничений по применению, их точность зависит, прежде всего, от дискретизации углового пространства, однако, требуют существенных вы-

числительных ресурсов. В КОМ переход к разностному аналогу уравнения

(5.6.1) происходит интегрированием (5.6.1) по контрольному объему и угло-

вому пространству.

	∫	∫	dI l	dVdΩ =	∫	∫ (− βI l + S l )dVdΩ
			ds
	ΔΩ l	V			ΔΩl	V
	S l (rR, s) = κ (rR)I (rR, s) +					σ (r )		I (rR, s′)Φ(s′, s)dΩ′
						R	∫
		ˆ		b	ˆ	4π		ˆ	ˆ ˆ
							4π

где l – l-ый контрольный телесный угол для КОМ.

Разностный аналог УРТ для КОМ

(5.7.4)

(5.7.5)

∑Iil Ai Dcil = (− βI l + S l ) VΔΩl	(5.7.6)
i=nb

172