Главная Журналы Поэтому сначала выполняются 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 |