Главная  Журналы 

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 [ 81 ] 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99

Поэтому сначала выполняются 15 итераций для решения уравнения для скорости, а затем 3 итерации - для температуры. Таким образом, значение LAST = 18.

Для турбулентной вязкости используется коэффициент релаксации REGAM.

OUTPUT. Практически все детали этой процедуры вам должны быть уже знакомы. Для турбулентного течения вместо величины /Re рассчитывается коэффициент трения /, который в программе обозначен через FRIC.

Решение находится для заданного значения числа Рейнольдса, равного 10. Так как средняя скорость w меняется на каждой итерации, постоянно пересчитывается молекулярная вязкость AMU, чтобы получить заданное число Рейнольдса (с физической точки зрения это равносильно подбору такой жидкости, которая при заданном градиенте давления даст желаемое число Рейнольдса). Измененное значение АМи выводится на печать после каждой итерации. Ожидается, что при приближении к сошедшемуся решению AMU станет постоянным.

PHI. Эта подпрограмма в общих чертах похожа на PHI в примере 11. Для получения турбулентной вязкости по (11.8) нужны градиенты скорости. Для этого первая итерация выполняется для ламинарного течения с GAM(I, J) равным молекулярной вязкости AMU. Впоследствии GAM(I, J) становится равным ц + где ц, задается

по (11.8). Величины FLUXIl (J, 1) и FLUXJl (1,1) берутся в качестве касательных напряжений на стенке канала.

Для вычисления производных от скорости в (11.8) значения w на гранях контрольных объемов находятся с помощью интерполяции. По каждому направлению эти значения обозначаются как WP и WM. Так как применяется неравномерная сетка, то грань контрольного объема лежит не посередине между соседними расчетными точками. Используется линейная интерполяция, основанная на конкретных размерах каждого контрольного объема и соответствующим образом модифицированная для приграничных контрольных объемов. Более точной процедурой для получения значений w на гранях контрольного объема является использование выражений вида (2.81), однако в данном случае нет необходимости в подобных усложнениях.

Чтобы регулировать изменения д, от итерации к итерации, вводится релаксация согласно (5.65) с коэффициентом релаксации REGAM.

После первой итерации рассчитывается среднее значение GAM (I, J) и градиент давления пересчитывается таким же образом, как и в процедуре, описанной в примере 11.

В уравнении для температуры GAM (I, J) берется равным к + к,, где турбулентная теплопроводность задается по (11.15). 256



11.3.3. Дополнительные имена на ФОРТРАНе

А1, А2 АМи

AMUT(I,J)

ANU AR ARG ASUM

CAPL(Al,A2)

COND CP DEN

DFACT(ARG) DH

DPDZ DTDZ DWDX DWDY FM, FP FRIC GAMAV

GAMAVP GAMT HP PR PRT QW RE

REGAM

REL RHOCP RL

RLX, RLY

T(I, J) TAUW

аргументы функции CAPL; вязкость ja;

турбулентная вязкость ja,; среднее число Нуссельта; площадь контрольного объема dA; аргумент функции DFACT; площадь поперечного сечения; длина пути смещения по Никурадзе Л [см. (11.12)]; теплопроводность к;

TSUM TWAV W(I, J)

- теплоемкость с;

- плотность р;

- функция Ван Дриста D [см. (11.11)];

- гидравлический диаметр D;

- продольный градиент давления dpIdz;

- то же температуры дТ/dz;

- dw/dx;

- dw/dy;

- интерполяционные коэффициенты;

- коэффициент трения /;

- среднее значение эффективной вязкости е

- значение GAMAV на предыдущей итерации;

- промежуточное значение турбулентной вязкости д,;

- обогреваемый участок периметра;

- число Прандтля;

- турбулентное число Прандтля;

- плотность теплового потока на стенке qy,

- число Рейнольдса;

- коэффициент релаксации для турбулентной вязкости;

- временная переменная 1-REGAM;

- объемная теплоемкость рс;

- итоговая длина пути смещения L;

- длины пути смещения и L;

- температура Т;

- касательное напряжение на стенке т,;

- среднемассовая температура Tj,;

- jwTdA;

- средняя температура греющей стенки Т.;

- продольная скорость w;



WBAR - средняя скорость в сечении w;

WM - скорость на грани;

WP - смоченный периметр или скорость на грани;

WSUM - dA ;

XPLUS(I,J) - безразмерное расстояние д;" [см. (11.13)];

YPLUS(I,J) - то же [см. (11.14)].

11.3.4. Листинг подпрограммы ADAPT

с СССССССССССССССССССССССССССССССССССССССССССССССССССССССС SUBROUTINE ADAPT

с--EXAMPLE 13 - TURBULENT FLOW IN A SQUARE DUCT

$ INCLUDE: COMMON

DIMENSION W(NI,NJ),T(NI,NJ),AMUT(NI,NJ),XPLUS(NI,NJ), 1 YPLUS(NI,NJ)

EQUIVALENCE (F(1,1,1),W(1,1)),(F(1,1,2),T(1,1)),

1 (F(l, 1, 3) ,AMUT(1,1) ), (F(l, 1, 4),XPLUS(1,1) ) ,

2 (F(l,l,5),YPLUS(1,1) )

DFACT(ARG)=l.-EXP(-ARG/26.)

CAPL(A1,A2)=A1*(0.14-0.08*(l.-A2/Al)**2-0.06* 1 (1.-A2/A1)**4)

ENTRY GRID

HEADER=TURBULENT FLOW IN A SQUARE DUCT

PRINTF=PRINT13

PL0TF=PL0T13

CALL INTA2(NCVLX,10,NCVLY,10)

CALL DATA4(XL,0.5,YL,0.5,POWERX,3.,POWERY,3.)

CALL EZGRID

RETURN

ENTRY BEGIN

TITLE(1)= W/WBAR •

TITLE(2)=•(T-TWAV)/(TB-TWAV)

TITLE(3)= TURB. VISCOSITY

CALL INTAS(KSOLVE(1),1,KPRINT(1),1,KPLOT(1),1,KPRINT(2), 1 1, KPLOT (2) , 1, KPRINT (3) , 1, KPLOT (3) , 1, LAST, 18)

CALL DATA5(AMU,1.,DEN,1.,DPDZ,-1.,RE,1.E5,REGAM,0.8) CALL DATA4(PR, 0.7 , CP, 1.,PRT,0.9,QW,1.) COND=CP*AMU/PR RHOCP=DEN*CP GAMAVP=AMU

С- SINCE THE ZERO DEFAULT VALUES OF W(I,J) AND T(I,J) ARE C- SATISFACTORY,

С THESE ARRAYS ARE NOT FILLED HERE. 258





0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 [ 81 ] 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99