Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена (1185899), страница 31
Текст из файла (страница 31)
В подпрограмме производятся следующие процедуры. Сначала обнуляются массивы 6 и Г, соответствующие глобальным матрице и столбцу. Далее организуется цикл по всем элементам, в котором: 151 1 2 з 4 5 е 7 3 в 19 11 !2 !3 14 15 !е 17 !3 13 29 21 22 23 24 25 26 27 23 23 29 3! 32 Зз 34 35 зе З7 33 33 49 41 42 43 44 45 -46 47 4В 43 59 51 52 53 Бпхно(>т)нк мкв(х,т,)и>м,н,ААА.ВВБ,МБ,О,Р> описАник ПАРЫ(вттов: и - число узлов н - число конкчных зллх(ентов х(м),т(м> — мАссиВН кООРдинАт узлоВ (нв(4«н) - ин>пкснАЛ ИАтвипА ААА,ВВВ - ИИКНА В>НННИХ ПОДПГОРРМВ( ! из - иивинА лкн)м лкнточнои МАТРИНН о — ип(од>п>я лвточння МАТРИПА Р - ВЫХОДНОИ ВЕКТОР-СТПЛБКП ПРАВЫХ ЧАСТЕИ ; В)П(ОДН>(К ( ПАРАМЕТРЫ о(млмз!Он х(и,т(и,!Нв(и,с(и,к(и с с вычислиник иивины лкнтн МБ с Мз.в 11 1 во ! ( >,н 1-)КО(Ы .>=1Ио(И.«!> К- НВ(Ы.«2> МБ МАХО(МЗ,1АВБ(1-)), .>АВБ(1-К> !АВБ(У-К» 1 И-Ы «4 с С ОЧИСПи МАССИВОВ ОУ с с с КАНАДО п>ло>А по ТРкугольникям с (! 1 ОО 12 Ь 1,И с д М«МБ«(2«М-МБ-И/2 ОО2 1 !У 2 С(И=9 ВО З 1=(,И 3 Р(И-9 НОМЕРА УЗЛОВ 1=1НВП).) у-!Ио(ы и К 1КО(И «2> ничисление НООРдинАт пльтРА ТРВХГОлъникА хс-(х(и«х(у >«х(к»!3 хс«(х(и«т(у)«т(к»73 Рис.
4.15 по)вичит(АеА ЛОРМНРОВАюя ИАТР>ив( иктодА конкчных ЭЛЕМЕНТОВ ДЛЯ ДВУМЕРНОИ СТАЦИОНАРНОИ ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ 54 55 56 57 58 59 66 81 62 83 84 85 68 87 88 68 76 71 72 73 74 75 76 77 78 79 89 З! 82 83 84 85 88 87 88 89 96 а! 92 93 94 95 96 97 98 98 166 !Е! 162 163 164 165 166 НРМВенение НООР>вин(У узлОВ н центру х)-хп >-хс ю х(з>-хс хз-х(н>-хс т(-т(П-тс Т2 Т(,1)-ТС тз-тлк>-тс ВЫЧИСЯЕНИЕ ИОЗЕЕИЦИИНШВ ВУННЦИИ ЕОУИЫ ВЫЧИСЛЕНИЕ ПЯОИЛ)В( ТРЕУГОЛЬНИКА А (Х2»ТЗ-ХЗ«12+Х1»Т2- «Х1» ТЗ+ХЗ»Т 1-Х2 «Т 1 >/2 вычнсяение зяпантов яокляьной илтуицы и вентоул пглв)в( члстяй ВЫЧИСЯЕНИЕ ЯОЦ)ц(и)ЕЛЬ Н(х СЯАГЫМК яян гглничных зя)мунтов Ь(Ь) 3 !Г((КР(ы.кц.е)со то 5 17(1КР(1Л.).ХС.()СО то 4 вкяля стоуоны л( 81 Т2-ТЗ С(-Хз-Хз 82 Тз-Т! СЗ-Х(-ХЗ ВЗ-Т)-Тг СЗ Хз-Х! СЛЫ.
ААЛ(АЬХ,ЛЬТ,Л(УУ,СТ,ХС,ТС> Н=4»Л Р-ЛЬУЧ Л/В 91! (Л(Л«51«В!»АЬТ«с!»С()/Н»Р С12 (Л)Л«81»82»ЛЬТ»С1«сз)/Н»Р/2 С!3 .(ЛЬХ«81«ВЗ+АЬТ«С1»СЗ)/Н+Р/2 92»2=(Л)Л«82»82»ЛЬТ«сз»сз)/Н»Р С23 ()АЛ»82»ВЗ»АЬТ«С2«СЗ)/Н»Р/2 СЗР (ЦЛ«83»83+А(Л»СЗ»СЗ)/Н+Р У! )Т»л/3 72-71 ТЪУ) СЛЬЬ Ввв(л(УЬ,ОЬ,(х(«)»Х(К))/2.(Т(«)»Т(К))/2) Н ЗСКТ((ХЗ-ХЗ)»»2«(Т2-ТЗ)«»2) Р АЬУЬ»Н/3 С22 922+В СЗЗ СЗЗ»Р Рис. 4.15.
Продолжение 197 198 :98 '19 И< И2 из И4 115 Ив 117 ив Ив 129 <г< 122 123 124 125 128 127 128 128 139 131 132 )ЗЗ 134 135 136 137 )ЗЗ !36 149 141 142 14З 144 145 148 147 148 146 159 151 152 153 154 155 !ВВ 157 156 153 С23 С23+В/2 О Сь»н/2 Р! Р<»В РЗ РЗ+В ВКЛАД СТОРОНЫ !Л 4 сньь ввв<АьРь,оь,<х<1)+х<л))/г,(х(!)+Т<л))/г! Н БОЛТ((Х1-Х2)»»2»(Х(-Х2)»»2) В-АЬЛЬ»Н/3 си~и»л сгг=сгг»п С12 С<2»9/2 О ЯЬ»Н/2 Р< РИО 72 Р2+Л ДОВАВЛПС<Е ЗЛЕМЕНТОВ ЛОКАЛЬНОМ МАТРНПЫ Н ЭЛЕМЕНТАМ ГЛОВАЛ5НОЙ 5 И1=1 И2 ! О СИ АББ)СИ 8 ТО ЬВЬ СО ТО 13 8 И2 Л Ь*С12 АЯБ!СИ 7 ТО ЬВЬ СО ТО !3 7 И2=К О С13 АЯБ1СИ 8 ТО ЬВЬ СО ТО 13 Вн! У Н2=3 О С22 АББ!СИ 3 ТО ЬВЬ СО ТО 13 3 И2=К О=С23 АББ!СИ 19 ТО ЬВЬ СО ТО 13 19 И! К Л СЗЗ АЯБ1СИ И ТО 1.81.
СО ТО 13 ТО ЕЕ ДЛН ВЕКГОРА ПРАВЫХ ЧАСТЕЙ И П 1) Р(1)»Р! Р(Л) Р<Л)»72 Р(н) Р(К)+РЗ Рис, 4.15. Продолжение 166 с НОнец Фклд пс тгютолъникьм !6! С 162 12 (А Ы+1 163 КБТОКИ !64 С Ий С ИОИВИ)МНИЕ МАПЭИВ В ВИКЕ, ПРИГОЛНСМ 186 С ДЛЯ ОБРАИЕНИЯ К ПОДПРОГРНВЕ МСНВ 187 С 168 !3 М! М1НЭ(И1,И2) 189 )и МАХЭ(И1,И2) 176 Ми 1»М2-М! 171 МБ! МБ»1 172 й(Б М-НБ 173 1ПИ1.се.ммБ)сс то 14 174 НИ 4В1» (М1-1)»)В 175 СО ТО 15 176 14 ЮВ М-М1»1 177 ИН НБ)»йй»МБ»МБ! /2-Мии» (МММ»1 ) /2МВ !78 15 С(ИИ) С(ИИ)»5 !79 СО ТО !БЬ,(8,7,8,9,18,11) 188 ВБ) Рис.
4.15. Продолжение по индексной матрице определяются глобальные номера узлов !', 1, й в данном элементе и принадлежность сторон границе; вычисляются коэффициенты функций формы элемента по (4.14) и рассчитываются коэффициенты локальных матрицы и столбца по (4.29), (4.30), значения всех коэффициентов локальных матрицы и вектор- столбца прибавляются к текущим значениям коэффициентов глобальных матрицы и вектор-столбца, индексы которых определяются глобальными номерами узлов 1, 1, /г.
В результате обхода всех У треугольников оказываются полностью сформированными глобальная матрица С и столбец правых частей Ф системы уравнений. Отметим, что при формировании матрицы С необходимо учитывать способ записи матрицы в машинной памяти для используемой стандартной подпрограммы решения системы линейных уравнений. В данном случае предполагается использование подпрограммы МСНВ из математического обеспечения ЕС ЭВМ (1б), реализующей метод квадратного корня для симметричных ленточных матриц. При зтом коэффициенты матрицы должны быть записаны в одномерный массив путем поп!едовательного обхода верхней части ленты над главной диагональю по строкам.
Такой пересчет индексов элемента матрицы в индекс одномерного массива реализован операторами 1бй †1. 155 ГЛАВА КОНЕЧНО-РАЗНОСТНЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ КОНВЕКТИВНОГО ТЕПЛООВМЕНА $ $.1. КОнечнО.РАЗнОстные схемы ДЯЯ УРАВНЕНИЯ ЭНЕРГИИ Класс задач конвективного теплообмена весьма широк (19, 23). В связи с этим в данной книге нет возможности остановиться на сколько-нибудь подробном разборе особенностей даже основных постановок задач, и поэтому мы рассмотрим методы численного расчета только для некоторых наиболее простых, но довольно часто встречающихся на практике задач решения уравнения энергии при заданном поле скоростей.
Для однофазного потока уравнение энергии можно записать в виде (5.1) где в левой части записана полная производная от энтальпии, включающая локальную и конвективную производные, а в правой части — члены, соответствующие тепловым потокам в единице объема, вызванным различными причинами. Первый член учитывает перенос теплоты за счет теплопроводности, второй — изменение энтальпии за счет работы сил давления, третий — выделение энергии за счет вязкого трения. При моделировании процессов конвективного теплообмена уравнение энергии должно рассматриваться совместно с уравнениями неразрывности, движения и состояния. При анализе многих процессов, например в случае свободной конвекции или при необходимости учета зависимости вязкости от температуры, необходимо все эти уравнения решать совместно.
Численные схемы для уравнений гидродинамики гораздо сложнее, чем рассмотренные в главе 3 схемы для уравнения теплопроводности. С ними можно познакомиться по книгам (19 — 21, 23). Мы будем считать, что поле скоростей известно, т. е. либо была решена гидродинамическая задача при условии отсутствия влияния на нее температуры, либо было проведено соответствующее экспериментальное исследование.
Поясним основные особенности численных схем для уравнения энергии, введя еще ряд упрощений в уравнение (5.1). В случае малых скоростей течения (числа Маха М (( 1) п следние два члена ср ( — + МАЧТ) = Ч ().Ч Т) . (5. 2) На первый взгляд схемы, рассмотренные в главе 3, легко перенести на уравнение (5.2). Действительно, оно отличается от уравнения теплопроводности только членом от7Т, содержащим первые производные от температуры по координатам, которые можно аппроксимировать конечными разностями.
Однако некоторые варианты такого «естественного» подхода приводят к неудачным численным схемам. Поэтому новый конвективный член вносит ряд существенных особенностей в процедуру выбора вида разностной схемы. Рассмотрим нх на примере простейшего одномерного стационарного уравнения энергии дТ д»Т срп — =)» —, (п=д„), (5.3) дх дх' Введем равномерную пространственную сетку х„= (и — 1) й, и = 1, ..., Л/. Конечно-разностное уравнение для внутренней точки будем строить методом баланса, выбрав элементарные объемы вида (х„— й/2, х„+ /»/2). Сеточную функцию численного решения обозначим, как обычно, через и„, п = 1, ..., /Ч.
Уравнения баланса п-го элементарного объема (рнс. 5.1) для единичного промежутка времени записывается так: 1 Энтальпия Я« (х„ + й/2) потока жидкости, выходящего из объе-— ма через сечение х„+/»/2 л з Энтальпия Я"( х„— — ) потока 2 ~ жидкости, входящего в объем через сечение х„ — Л/2 Количество теплоты 9'(х„ + й/2), + выносимое путем теплопроводности из элементарного объема через сечение х„+Ь/2 Количество теплоты Я'(х„— /»/2), вносимое путем теплопроводности в элементарный объем через сече- ние х„— й/2 или 1~" (х„.+ й/2) Як (х» — /г/2) = Я' (х» — й/2) Ят (х„+ й/2). (5.4) в правой части (5.1) можно не учитывать и рассматривать уравнение энергии вида Тепловые потоки ()' (х„+ й/2) и Я'(х„— И/2) аппроксимируются так же, как и в случае одномерного уравнения теплопроводностн, см. (3.44): Я'(х„+й/2) = Л "" "" Я'(х — й/2) Л "" ' "" (5.5) Дтя аппроксимации 9" (х„+ И~2) и /З«(х„— й/2) можно предложить несколько вариантов: «/«(х«+й/2) прои„»„Я'(х„— й/2) ж срои„ Я" (х„+ й/2) ж сро(и„«,+и„)/2, Я" (х„— й/2) сро (и„+ и„,)/2, Я" (х, ( й/2) срои„, 9" (х„— й/2) ж срои„..
(5.6) (5.7) (5.8) Проанализируем соотношения (5.6) — (5.8), учитывая, что жидкость течет в направлении оси х, проходя через элементарные объемы в порядке возрастания их номеров. В аппроксимации (5.6) принято, что температура втекающей в объем У„жидкости равна температуре центра этого элементарного объема, а температура вытекающей — температуре центра следующего за ним по течению объема У„„. Это приближение является весьма неудачным, поскольку получается, что на тепловой режим элементарной ячейки )', не оказывает никакого «конвективного» влияния температура жидкости в предшествующем ей объеме У„п но воздействуеттемпература следующего объема У„,о Однако физика процесса конвективного переноса такова, что температура жидкости в объеме У„должна определяться только условиями ее протекания и теплообмена в предшествующих объемах, т.
е, жидкость при течении собирает и несет дальше информацию о своем прошлом, но не знает о своем «будущемм Формулы же (5.6) противоречат этому факту, и поэтому их использование приводит во многих случаях к абсурдным результатам расчетов. При использовании формул (5.7) принимается, что в элементарный объем У„втекает жидкость с температурой, равной полусумме температур данного и предыдущего объемов, а вытекает жидкость с температурой, равной полусумме температур данного У„ и следующего за ним по течению объема У„«о Принципиальных возражений против аппроксимации вида (5.7) нет.