CC BY
79
_________УЧЕНЫЕ ЗАПИСКИ ЦАГИ
Том XXVII I 199 6 1
№3-4
УДК 629.7.015.3:533.6.011.6
НЕСТАЦИОНАРНАЯ ЗАДАЧА СОПРЯЖЕННОГО ТЕПЛООБМЕНА ДЛЯ ПЛОСКИХ И ОСЕСИММЕТРИЧНЫХ ТЕЛ
В. А. Башкин, С. М. Решетько
Описана математическая постановка нестационарной задачи сопряженного теплообмена для плоских и осесимметричных тел ппи обтекании их сверхзвуковым потоком совершенного газа. Изложена процедура численного анализа нестационарного уравнения теплопроводности. Работоспособность составленного комплекса программ протестирована на ряде частных задач.
При сверх- и гиперзвуковых скоростях полета на температурный режим обтекаемых поверхностей заметное влияние могут оказывать теплофизические свойства материала стенки. Для учета этого влияния необходимо решать задачу сопряженного теплообмена; при больших числах Рейнольдса она сводится к совместному решению уравнений пограничного слоя, описывающих поле вязкого течения (внешняя задача), и уравнения теплопроводности, описывающего поле температуры в оболочке тела (внутренняя задача).
Решение стационарной задачи сопряженного теплообмена показало, что теплопроводность материала стенки может оказаться эффективным средством управления максимальной температурой обтекаемой поверхности сверх- и гиперзвукового самолета (см., например, [1—3]).
В практических приложениях очень часто приходится сталкиваться с анализом нестационарной задачи сопряженного теплообмена (проведение теплового эксперимента в аэродинамических трубах с применением термоиндикаторов, движение снарядов и ракет по траектории, переходные режимы движения ЛА и т. п.). Однако численное решение этой задачи связано с выполнением больших объемов вычислений на ЭВМ. Для сокращения временных затрат часто используется предположение об одномерности внутренней задачи, связанной с решением уравнения теплопроводности. Такой подход был использован, например, в [4—6] и позволил определить начальную стадию прогрева движущегося в атмосфере тела; при этом внешняя задача решалась в квазистационарном приближении. Однако такая приближенная поста-
новка нестационарной задачи сопряженного теплообмена имеет ограниченные рамки применения.
Настоящая работа посвящена численному анализу нестационарной задачи сопряженного теплообмена для плоских и осесимметричных тел, обтекаемых сверх- и гиперзвуковым потоком совершенного газа. Сначала рассмотрены условия, при которых внешняя задача (уравнения Прандтля) может решаться в квазистационарном приближении. Затем сформулирована математическая постановка задачи и кратко описана процедура ее численного интегрирования. В заключение на ряде частных задач протестирована работоспособность созданного пакета программ.
I. Для уравнений Прандтля нестационарность течения обусловливается двумя причинами: неустановившимсяг движением тела (внешние граничные условия) и изменением температуры поверхности во времени (внутреннее граничное условие для уравнения энергии). Уравнения Эйлера можно рассматривать в квазистационарном приближении, если выполняется условие
где *1 = £/Ко — характерное время пребывания жидкой частицы около обтекаемой поверхности, — характерное время движения тела, X — характерный линейный размер тела, К, — скорость набегающего потока, вИ — число Струхаля. Для условий сверх- и гиперзвукового полета ЛА условие (1.1) обычно выполняется.
В твердом теле изменение температуры характеризуется коэффициентом температуропроводности а^|м2/с| = Хд / {сд рд), т. е. определяется теплофизическими свойствами используемого материала: теплопроводностью Хв, ПЛОТНОСТЬЮ рв и теплоемкостью Св- Пусть Ьв есть характерный линейный размер по прогреву тела (в общем случае он отличается от линейного размера Ь). Тогда с помощью указанных характерных величин можно определить характерную скорость изменения температуры а0в~ав1^в и, следовательно, характерное время прогрева тела
Для того чтобы течение газа в пограничном слое можно было считать в главных членах квазистационарным, время пребывания газовой частицы около обтекаемой поверхности должно быть много меньше времени прогрева тела
Для газовой среды характерными параметрами подобия являются числа Рейнольдса Яе и Прандтля Рг, которые связаны между собой соотношением
БЬ = « 1,
(1.1)
(1.2)
Не = _ У«Ь
Цоо ^^«0
где а2„ = *-»/(роосроо) — коэффициент температуропроводности газовой среды. Тоща условие квазистационарности (1.2) преобразуется к следующему виду:
ifА. ‘ «!
LrJ а2оо Pr.Re
(1.3)
Таким образом, если одновременно выполняются условия (1.1) и (1.3), то течение газа в пограничном слое можно рассматривать в ква-зистационарном приближении, а эффекты нестационарносги проявляются только при решении уравнения теплопроводности.
2. При постановке задачи сопряженного теплообмена для плоскот го или осесимметричного тела конечной длины, движущегося со сверхзвуковой скоростью в плотных слоях атмосферы под нулевым углом
тельно, уравнения Прантдля можно рассматривать в квазистационар-ном приближении.
Уравнения пограничного слоя записываются в системе координат xi, у, связанной с обтекаемой поверхностью: Х2 — координата, направленная вдоль направляющей или образующей тела и отсчитываемая от критической точки; у — координата, ортогональная обтекаемой поверхности и отсчитываемая от нее вдоль внешней нормали (см. рис. 1, а).
Уравнение теплопроводности записывается в ортогональной системе координат хь х2 (см. рис. 1, а), где — координата, задаваемая как расстояние линии = const от критической точки в сечении х^ — = 0. Выбор координаты х^ зависит от конкретной конфигурации внешней и внутренней траниц тела.
Для численного интегрирования системы уравнений Прандтля вводится безразмерная функция тока путем тождественного удовлетворения уравнению неразрывности, а уравнения импульсов и энергии
атаки, предполагается, что условия (1.1) и (1.3) выполнены, и, следова-
\В
с
О х,
xt
В
С
А
В х.
Рис. 1. Системы координат и области интегрирования: а — физическая область; Ь — расчетная область
преобразуются к обобщенным переменным подобия. Уравнение теплопроводности также записывается в безразмерном виде. В результате этих преобразований приходим к следующей системе уравнений:
«щах = (2Й£0)1/2, ка = к = г1Н1Н2рвсв > Ф = г]НхН2Ф,
Ма
и,к,р,\1,ср — соответственно продольный компонент вектора скорости, энтальпия, плотность, коэффициент динамической вязкости и теплоемкость при постоянном давлении газа, Я — характерный размер (радиус затупления), г — расстояние от оси симметрии, у = 0 или 1 соответственно для плоского или осесимметричного тела, Т — температура, 5са — криволинейные ортогональные координаты, х — время, Ф — мощность источников тепловыделения, #а — коэффициенты Ламе; нижние индексы 0, у/, Е и т — соответственно обозначают значения газодинамических параметров торможения на обтекаемой поверхности, внешней границе пограничного слоя и для турбулентного пото-
Рх[лг(1 + Ит)/Г]' + Л/|" - Ф? -*)-
(2Л)
V
/
(2.2)
Здесь введены обозначения:
ка; верхним индексом 0 обозначены соответствующие масштабные величины.
Начальное условие накладывается только на решение уравнения теплопроводности (2.3) и характеризует собой поле температуры в начальный момент времени: Т (ха, 0) = Т\ (д^).
Внешние граничные условия для уравнений импульсов (2.1) и энергии (2.2) выражают собой условия сращивания с внешним невязким потоком, а внутренние граничные условия для уравнения импульсов есть условия непротекания и прилипания на обтекаемой поверхности:
Внутреннее храничное условие для уравнения энергии (2.2), а также граничные условия для уравнения теплопроводности (2.3), решение которого находим в физической области АВСВА (см. рис. 1, а), представляют собой условия локального баланса тепла. В переменных *1, Х2 область интегрирования представляет собой четырехугольник, и граничные условия на 1раницах этой области можно записать в общем виде (а, (3 = 1, 2):
Не нарушая общности, граничное условие (2.4) можно записать в квазилинеаризованном виде:
где ка> зе+а, g+a — некоторые известные функции пространственных координат, температуры и времени.
В задачах внешней аэродинамики обычно предполагается, что воздух в поле возмущенного течения является оптически прозрачным и с обтекаемой поверхности в окружающее пространство происходит излучение тепловой энергии по закону Стефана — Больцмана. Тогда уравнение теплового баланса на обтекаемой поверхности тела (граница АВ) примет вид:
/[{х2, оо) = g(x2, оо) = 1, /г(х2, 0) = /^(х2,0) = 0.
кадТ/дха =зе_аТ - g_a, х -кадТ/аха=ве+аТ-8+а,
1ж--о2{т*-т*) = 0
где
Р ЕиЕРЕг1
__ « /л
-®1 = ^В Р^(-^Р£0^тахМ-£о) >
Здесь ст — постоянная Стефана — Больцмана, б — степень черноты поверхности; индекс оо обозначает параметры набегающего потока. Величины А и являются параметрами подобия и представляют собой отношение характерных масштабов кондуктивного и радиационного потоков тепла к конвективному соответственно.
j На торцевой поверхности (граница ВС) обычно ставится условие теплоизолированности F22 — дТ/дх2 = 0, а на оси симметрии (граница DA) — условие симметрии F2\ — дТ/дхъ = 0.
На внутренней поверхности тела (граница CD) условие локального баланса тепла зависит от конкретных условий решаемой задачи. Во многих практически важных случаях можно ограничиться предположением о теплоизолированности внутренней стенки: F\2 = дТ/дп = 0.
Для указанных граничных условий коэффициенты соотношений (2.5) имеют следующий вид:
При решении уравнения теплопроводности в общем случае в расчетной области допускается конечное число поверхностей разрыва, наличие которых связано с использованием неоднородных по теплофизическим свойствам материалов. На этих поверхностях ставятся условия совместности — непрерывности температуры и тепловых потоков:
_ хдТ
д. дп дп
3. Численное интегрирование уравнений Прандтля проводится конечно-разностным методом [7] с помощью стандартной программы [8]. Динамическая вязкость ц для совершенного газа вычисляется по формуле Сазерленда.
Уравнение теплопроводности записывается в строго консервативной форме с последующей дискретной аппроксимацией разностной схемой интегроинтерполяционного типа.
Методом суммарной аппроксимации [9] начально-краевая задача (уравнение (2.3) + начальные + граничные условия) сводится к цепочке одномерных задач, каждая из которых (при а = 1 и а = 2) рассматривается на временнбм интервале тт+(а_^/2 < Т < хт+а/2> тп = 1,2,...:
эе+1 = эе+2 = g+i = g+2 - о, ae_! = 4XBr^H2T^D2/Dl,
г-i = *о2(зт' + £)].
О) — ?!> 7(і)іх>хт) ~Т^2){х9^т)> 111 — 1,2,...;
Т(2){х> %т+уі) = ^(і)(*> Хт+Ці\і И = 0,1,...
и граничными условиями
^а^Т^/й^д = аЄ-вї(а) — £_а, ха — 0 5 —кадТ(а}/дха — ае+а^(а) ~ £+а» /*а
_ ^шах . Я '
Функция источников ф должна удовлетворять условию нормировки:
Ф = Ф1 +Ф2- . ;
На сетке ой = юх х юх, .•
где
со.
= \х=\хих2)ев, - £йа>/Ь іа= 2,...,І\Г„, а = 1,2
к= 1
произвольная сетка на отрезке 0 £ ха ^ х^1ах с шагами Ль,*-ъ ^ = = 2, ... ,ЛГа, .
= Ьи = тт-1 + Лтт> И = 1,..., И1шах}^
приходим к следующей разностной начально-краевой задаче:
Аут = Л(т)(сгр + (1 - ст)у) + фа , (3.1)
*ад(^д +(1-0)^,1 ) = ае-а(оУі + (1-®)л)-?-а +0,5/^д^Т)1,
*а,ЛГа (°Уха,Ка + (1 - ст)^,АГа ) = (СТ^а + (1 - <*)УКа ) ~ «+а + °>5Аа,^, Фт,ЛГ •
Здесь у — сеточная функция температуры Т, ух = (у - у)/&хт, У * У(«]а/2» У = ’ Ати = Х«+1 " г*-> АУ - (^*а )х» Фі + Ф2 = Ф>
Ё-а - Я-а + 0,5АаДфа(0, х), ї+а = £+а(ї) + 0,5*,Л-іФ«(*Г“, *)’
х = хт+ 0,5Атт , функции £, зе, # и ка берутся в момент времени т . Приведем систему уравнений (3.1) к виду для расчета на ЭВМ:
Аг Уі-1 - Сі Уі + 4-+1 Уі+1 = і = 2,..., -1,
Уі = «і у2 + цІ5 = ає2 Уі\га-і + И2>
у4,- = стАг0()//Аа>,-_2, С{ = + А,-+1 + кЬа>]-1/&1т,
— {^о,,-1/Ат — (1 — ст) {ка {/Ьа,1-1 + ^а,1+1 /^а,» )}у/ + (1 — ст) х
х (^а,( У{-\ / ^а,1-1 + ^а,»Ч1 У/'+1/^а,/) + ^а,1-1Фа,1>
«1 = в*^,2(*)/[(*а,2(*) + Аа,1 *-«(*))<* + АЛвд/(2Ат)],
®2 = (*)/[(*^,ЛГв (т) + На>М^Х^а{х))о + **в>Лв_1/(2Ах)],
_ (1 - д)[*а,2(*)У5,,2 - ае-а(т)у1] + 0,5/|ад куг/Ьт + &-а И ^а^АхД + «-«(*)]+ АадМ2Дт) ’
(1 - ст)[-£„,лга (т)Уха,ЛГа - ае+а(^)УАГа ] + °>5^а,^а-1 Фуа/А» + 1+а
^[^^(тУЛа^.!+ае+а(т)] + Ла1Лга_1^Д2Дх) ’
^а,/' = (Аа,1-1 + ^а,/)/2-
Используемая схема с весами (0 й ст < 1) имеет первый порядок точности по времени и второй порядок точности по пространственным координатам. Нелинейная задача (излучение, переменность теплофизических свойств и т. п.) решается с привлечением итерационного процесса до выполнения заданной точности. Система разностных уравнений цепочки решается методом потоковой прогонки [9], который хорошо зарекомендовал себя при решении прикладных задач.
Для сокращения вычислений пространственно-временная сетка дается с разрядкой по формулам:
Ьх* =сьхАха~1’ Ах" =сДтАхи_1,
где Сд*, сДх — некоторые постоянные, значения которых подбираются эмпирическим путем (сд* £ 1, сДх > 1); к = 2, ... — порядковый номер шага, отсчитываемый от поверхности ха = 0; и = 2, ... — порядковый номер слоя, отсчитываемый от х = 0; Дх* и Ах1 — начальные шаги. В общем случае с учетом условий ограничения Ах* < , Дх” < Дх1®3*
пространственно-временная сетка является многопараметрической; правильный выбор ее параметров позволяет получать надежные результаты при значительной экономии ресурсов ЭВМ. Значения этих параметров для конкретной задачи, близкие к оптимальным, можно определить на основе предварительных расчетов.
Для сокращения объема вычислений часто используется подход, в основе которого находится известный результат того, что местное число Стантона относительно слабо зависит от температурного фактора. Поскольку при наличии излучения местная температура поверхности меньше местной температуры восстановления, то число Стантона можно ввести следующим образом:
£ __________в»______
н реые(^ео~К)
При таком подходе уравнения Прантдля интегрируются для одного или нескольких контрольных моментов времени и определяется распределение местного числа Стантона вдоль обтекаемой поверхности. Для решения внутренней задачи в последующие моменты времени используется распределение числа Стантона, вычисленное по данным последнего контрольного момента времени. В этом случае коэффициенты соотношений (2.5) принимают вид '
®+1 = ®+2 = 8+1 = 8+2 = °> ж-1 = + РгСн Кен реые) /А»
£-1 = 'квг^Н2 РгСн Яея рЕиЕ +П2(ЪТ* + 7^)| /Бъ
где Ле# = Ищах-КРяо/м^о-- 4. Для проверки работоспособности и эффективности описанной
численной процедуры были проведены методические расчеты, которые можно разбить на две группы^
К первой группе относятся задачи о прогреве тел относительно простой конфигурации. Они при сделанных предположениях имеют точное решение (в квадратурах или в виде рядов), и на их Примере анализировалась эффективность работы программы по решению уравнения теплопроводности. < •
Ко второй группе относятся частные задачи по нестационарному сопряженному теплообмену при обтекании тела сверхзвуковым потоком совершенного газа при ламинарном течении в пограничном слое.
Отметим, что расчеты были выполнены с использованием чисто неявной схемы (сг = 1). Результаты расчетов при варьировании параметров (влияющих на точность расчета) могут служить основой для рационального выбора пространственно-временной сетки.
5. Работоспособность программы по численному анализу уравнения теплопроводности исследовалась на трех частных задачах о прогреве тел простой конфигурации с постоянными теплофизическими свойствами материала.
Первая задача связана с прогревом полубесконечного тела при заданном значении теплового потока и имеет точное решение в явном виде [10]. Сравнение численного и аналитического решения (рис. 2) показало, что погрешность, принимающая наибольшее значение вблизи начального момента времени, быстро уменьшается при Ах -> 0, так
что при Ат < 0,1 кривые на графике Ту, (т) практически сливаются в одну.
Вторая задача имеет дело с нагревом прямоугольной пластины размерами 2Ь х 23, имеющей начальную температуру Г} и помещенную в неограниченную среду с постоянной температурой Ть.
Сформулированная выше задача зависит от двух параметров подобия — теплофизического Ы = аЬ/Хд, где а — коэффициенту теплоотдачи, и геометрического В = В/Ь (характеризует область инте1рирования). Эта задача, имеющая аналитическое решение в виде рядов [10], была численно решена для следующих значений параметров подобия: = 5/3, В = 2/3. При расчетах варьировались шаг по времени Дт и число узлов пространственной расчетной сетки, при этом число узлов в продольном и нормальном направлениях принималось одинаковым. Хотя рассматриваемая задача имеет две плоскости симметрии, численное решение ее находилось сразу для всей области интегрирования.
Сопоставление точного и численного решений задачи для точки, расположенной в центре пластины, представлено в таблице. Результаты расчетов показывают, что наибольшая погрешность счета наблюдается вблизи начального момента времени и по мере прогрева пластины она уменьшается. Естественно, на более мелкой сетке точность расчета в начальные моменты времени существенно возрастает. Переход от сетки 21 х 21 к сетке 41 х 41 приводит к сравнительно небольшому повышению точности расчета (в третьей значащей цифре). Поэтому можно считать, что этот резерв повышения точности расчета практически исчерпан и целесообразно ее повышать за счет изменения величины временного шага. Уменьшение последней приводит к повышению точности расчета, и при наименьшем временном шаге Дх * » 1,7 • 10'3 численное решение с хорошей степенью точности согласуется с точным решением даже при малых временах прогрева. Сравнение результатов решения этой задачи с аналогичными данными [11], полученными с использованием локально-одномерного метода (ЛОМ) и метода последовательного энергетического баланса (МПЭБ), показывает, что в целом метод суммарной аппроксимациии более эффективен, чем названные методы.
Рис. 2. Прогрев полубесконечного тела. Условия расчета:
.а2 • . 1 _ т-4 • *..тах
Дх =10
0,1;
„ - -А*
— Дт = 0,5;------------Дт *• 0,1; - • - • —
Дт = Ю'2; ------Дт " 10'3;-------------
точное решение; Д7\, — Ткгоч — Ткчис
X Точное Метод суммарной аппроксимации ЛОМ МПЭБ
решение Дх = Дх0 = 0,006945 Дхд/2 Дх0/4 Дх = Дхо Дх = Дхо
Щ х N2 11 х 11 21 х 21 41 х 41 41 х 41 41 х 41 И х 11 11 х 11
0,0139 0,0000 0,0006 0,0005 0,0005 0,0002 0,0001 0,0017 0,0006
0,0278 0,0010 0,0038 0,0033 0,0032 0,0021 0,0016 0,0081 0,0009
0,0417 0,0061 0,0110 0,0102 0,0099 0,0081 0,0072 0,0199 0,0051
0,0556 0,0166 0,0226 0,0215 0,0213 0,0191 0,0179 0,0367 0,0105
0,1111 0,0984 0,1027 0,1019 0,1017 0,1001 0,0993 0,1345 0,0625
0,1667 0,2057 0,2067 0,2065 0,2065 0,2061 0,2060 0,2478 0,1575
0,2222 0,3119 0,3107 0,3110 0,3110 0,3115 0,3117 0,3552 0,2666
0,2778 0,4082 0,4059 0,4064 0,4065 0,4073 0,4078 0,4505 0,3807
0,3333 0,4928 0,4901 0,4905 0,4907 0,4917 0,4923 0,5330 0,4683
0,3889 0,5660 0,5632 0,5636 0,5637 0,5649 0,5655 0,6037 0,5333
0,4444 0,6290 0,6262 0,6266 0,6267 0,6278 0,6284 0,6639 0,5911
0,5000 0,6830 0,6803 0,6806 0,6807 .0,6819 0,6824 0,7150 0,6515
0,5556 0,7292 0,7267 0,7269 0,7270 0,7281 0,7287 0,7584 0,7110
0,6111 0,7687 0,7664 0,7666 0,7666 0,7676 0,7682 0,7952 0,7404
0,6667 0,8024 0,8003 0,8004 0,8005 0,8015 0,8020 0,8264 0,7728
Третья задача [12] связана с прогревом кругового цилиндра конечного размера. В указанной работе отмечено, что вследствие особенности в точке т = О решение задачи обычным методом Фурье [10] при малых временах затруднительно, поэтому эта задача была решена методом разложения в ряд по малому параметру (при малых числах Фурье). В настоящей работе рассматриваемая задача была решена методом суммарной аппроксимации для тех же значений параметров подобия: 1/Я — 1, В1 = 1 и 10. Расчеты были проведены как на равномерной, так и на неравномерной сетке по пространственным координатам; параметры сетки не варьировались, а были подобраны на основании результатов расчета рассмотренных выше задач (рис. 3). При решении задачи Коши шаг Ат был неравномерным, при этом начальный шаг был подобран достаточно малым (Ат = 10'7). В качестве примера на рис. 3 показано изменение температуры в двух точках на боковой поверхности цилиндра: г = 1 (зависимости 1) и г = 0,5 (зависимости 2).
При В1 = 1 результаты численных расчетов на обеих сетках практически совпадают между собой и хорошо согласуются с аналитическим решением задачи. Связано это с тем, что при данном числе Био возмущение температуры поверхности становится заметным при относительно больших временах по сравнению с начальным шагом интегрирования.
При В1 = 10 имеем уже другую ситуацию: возмущение температуры поверхности происходит при относительно малых временах, сравнимых с начальным шагом интегрирования. Вследствие этого при относительно малых временах численное решение заметно отличается от аналитического; с ростом времени прогрева погрешность расчета уменьшается. Как и следовало ожидать, погрешность расчета на рав-
Ш-10
1%х -В -5 -V -з Точка
• 0,9763 0327 олт 0,5301 1
0,9831 Ц9657 орт №89 2
о 0,9981 Ц9Ш тп 0,5277 1
0,3980 0,3821 00039 07253 2
0,9375 0,9150 №31 1
№75 0,95В №72 о,т 2
Г‘ 2 1
- / \3
)
21 г
Рис. 3. Прогрев кругового цилиндра конечного размера:
--- — аналитическое решение [10); численное решение: •• — сетка
82 х 102 с разрядкой: Дт‘ = 10'7; с41 = 1,2; Дх" = 5 х 10‘5, с** = 1,2,
Л*тах = 5 х 10‘2, о о — равномерная сетка 101x201: Дт° = 10‘7; сд, = 1,2
номерной сетке заметно больше, чем на неравномерной. Это еще раз подтверждает, что использование неравномерной сетки является эффективным средством снижения погрешности счета.
6. Осесимметричные оболочки в виде сферически затупленных круговых конусов конечной толщины часто используются в прикладных задачах, поэтому полная задача по нестационарному сопряженному теплообмену была рассмотрена на примере прогрева конической оболочки в сверхзвуковом потоке совершенного газа при движении под нулевым углом атаки с постоянной скоростью, соответствующей числу М„ — 6. Теплофизические свойства материала оболочки принимались постоянными. Данные по параметрам невязкого потока, необходимые для расчета пограничного слоя, заимствовались из [13].
Уравнения Прандтля интегрировались численно с помощью стандартной программы [8], которая широко используется в вычислительной практике, и влияние параметров расчетной сетки на точность расчетов, которое исследовалось на примере обтекания тел с изотермической и теплоизолированной поверхностями, а также при наличии излучения с обтекаемой поверхности, хорошо известно. Поэтому в предпринятых методических расчетах сетка в пограничном слое не варьировалась, а ее параметры были выбраны по условиям экономичности счета такими, чтобы максимальная погрешность по тепловому потоку не превышала 2%. Варьированию подлежали только параметры сетки для внутренней задачи.
Влияние параметров пространственно-временной сетки и итерационного процесса на качество получаемых результатов исследовалось
на примере обтекания конуса с полууглом раствора сок = 5* сверхзвуковым потоком (М* = 6, высота Н = 20 км) при наличии излучения с обтекаемой поверхности согласно закону Стефана — Больцмана. Так как сетка в общем случае является многопараметрической, приведем результаты исследования влияния на точность решения задачи лишь некоторых основных параметров: временного шага Ат, шага по нормальной координате Дхі и точности итерационного процесса е.
ВлияНие временного шага интегрирования Ат при фиксированной неравномерной сетке на точность расчета изменения приращения температуры в критической точке на внешней (Аі = Ту, — Т\) и внутренней (Д2 = Г» — Т{) поверхностях оболочки показано на рис. 4. Хотя величина временного шага изменялась в пять раз, это очень мало сказалось на точности расчета величин Аі и Д2: максимальное влияние имеет место при малых временах прогрева и составляет менее 1% для Аі и около 50% для Д2. В последнем случае влияние в процентах, конечно, велико, но эта максимальная погрешность связана с точностью определения времени достижения теплового возмущения внутренней поверхности оболочки. При Последующем возрастании времени проірева погрешность быстро уменьшается и при т = 0,5 максимальное влияние шага интегрирования Ат составляет менее 0,1% для Ді и около 1% для Д2.
Изменение точности є сходимости итерационного процесса в исследованном диапазоне ее изменения также оказывает слабое влияние на точность расчета величин Ді и Д2 при фиксированном временном шаге и пространственной сетке. При этом следует иметь в виду, что алгоритм решения задачи построен таким образом, что, если итерационный процесс не сходится за заданное максимальное число итераций, временной шаг делится пополам, так что величины Ат и є не являются абсолютно независимыми.
Влияние пространственной сетки на точность расчета Ді и Д2
Рис. 4. Приращение температуры в критической точке на внешней (Д[ = Тщ - 2}) и внутренней (Дг = Т. - Тх) поверхностях оболочки. Условия расчета:
Дт1 = 0,025, Дх”” = 0,15, в = НГ4; 1 - = 1,
Дх1 = 0,15; 2 - Сд, = 10, Дх1 = 0,15; 3 -с^ = 2, Дх1 =0,075
показано на рис. 4. Как и следовало ожидать, равномерная сетка (кривые Г) дает наибольшую погрешность расчета. На неравномерных сетках в указанных пределах изменения ее параметров получаются очень близкие результаты и обеспечивается существенно более высокая точность расчета величин Д1 и Д* Роль неравномерности пространственной сетки возрастает по мере уменьшения теплопроводности оболочки (уменьшение параметра Б{), поскольку возрастают градиенты температуры в оболочке тела.
Термодинамическая система поток — тело является замкнутой, поэтому естественно ожидать установления характеристик системы при т -> да. Рассмотрение этого процесса является необходимым для проверки надежности работы программы и для подтверждения возможности решения с ее помощью стационарной задачи сопряженного теплообмена. Процесс установления был рассмотрен для случая аэродинамического нагревания конической однослойной оболочки. При указанных выше условиях полета значение параметра Бг — 0,2315. Расчеты прогрева оболочки были выполнены для двух случаев, соответствующих материалам с малой и большой теплопроводностью и характеризуемых значениями параметров подобия В\ = 3,59 • ДО"5 и 3,59 соответственно.
В качестве примера на рис. 5 приведены распределения температуры вдоль обтекаемой поверхности конуса в различные моменты времени; здесь же сплошными линиями нанесены зависимости [3], полученные для стационарного режима для различных значений параметра подобия 2>1.
Для слаботеплопроводной оболочки (Х>1 = 3,59 • 10'5) наблюдается очень резкое изменение температуры внешней поверхности: в окрестности затупления она почти мгновенно выходит на соответствующее стационарное значение; на конической поверхности процесс установления протекает более медленно, и для указанного момента времени максимальное отклонение от стационарного режима составляет примерно 10%. В последующие моменты времени происходит медленный прогрев оболочки практически при заданном распределении температуры на ее внешней поверхности.
Для оболочки с большой теплопроводностью (Дг = 3,59) наблюдается существенно меньший темп
¥
Лт
о,в
О.*
о,г
Ь
ЛИ? \ 2,вС
------------
о
2,0
Рис. 5. Распределение температуры вдоль обтекаемой поверхности конуса.
Условия расчета:
М„ = 6, В = 20 км, Мг = 0,2315,
=о,25. Решение стационарной (--------
[3]) и нестационарной (- - -х 10-5) задач
А - 3,59
изменения температуры ее внешней поверхности, поскольку за счет теплопроводности происходит большой отток тепла внутрь оболочки. Температура постепенно выходит на свое стационарное значение, причем это установление происходит практически одновременно как для затупления, так и для конической поверхности; в последнем случае максимальное отличие от стационарного решения не превышает 2%.
В [14] получено решение нестационарной задачи сопряженного теплообмена для сферически затупленного кругового конуса с полууг-лом раствора юк = 15°, обтекаемого сверхзвуковым потоком совершенного газа при Re^2- Х/Хв = 0,314. Тело представляет собой оболочку
постоянной толщины из материала с постоянными теплофизическими свойствами. Для интегрирования уравнений пограничного слоя использовался метод И. В. Петухова [7], а внутренняя задача — интегрирование уравнения теплопроводности — решалась конечно-разностным методом с расщеплением по пространственным координатам.
Для тех же исходных данных задача была решена по описанной выше методике. На рис. 6 приводится сравнение распределений температуры вдоль обтекаемой поверхности для различных моментов времени. На сферическом затуплении в целом результаты расчетов согласуются между собой, а наблюдаемое различие данных в окрестностях критической точки связано, по-видимому, с выбором расчетной сетки. На конической поверхности расхождение результатов можно объяснить неполным совпадением исходных данных: например, в работе [14] не указаны продольная длина области интегрирования и число Прандтля газа.
Таким образом, приведенные результаты расчетов показывают, что разработанная программа позволяет не только проследить эволюцию прогрева оболочки, но и получить стационарное решение задачи.
ЛИТЕРАТУРА
1. С аре у Е. С. Alleviation of leading-edge heating by conduction and radiation/ARC R & M 3540. - 1968.
2. Б a hi к и н В. A., P e ш e т ь к о С. М. Расчет максимальной температуры затупления с учетом теплопроводности материала//Ученые записки ЦАГИ. — 1989. Т. 20, № 5.
3. Башкин В. А., Р е ш е т ь к о С. М. Температурный режим затупленных клиньев и конусов в сверхзвуковом потоке с учетом теплопроводности материала стенки//Ученые записки ЦАГИ. — 1990. Т. 21, № 4.
4. Б а ш к и н В. А., Р е ш е т ь к о С. М. Нестационарная задача сопряженного теплообмена для затупленных тел в сверхзвуковом потоке//Ученые записки ЦАГИ. — 1987. Т. 18, № 3.
Рис. 6. Распределение температуры вдоль обтекаемой поверхности конуса. Численное решение:
------работа [14],--------данная работа
5. Зинченко В. И., Тр о фи мчу к Е. Г. Решение неавтомодельных задач теории ламинарного пограничного слоя с учетом сопряженного теплообмена//Изв. АН СССР, МЖГ. — 1977, № 4.
6. Зинченко В. И., Путятина Е. Н. Решение задач сопряженного теплообмена при обтекании тел различной формы//ПМТФ. — 1988, №2.
7. П е т у х о в И. В. Численный расчет двумерных течений в пограничном слое//В сб.: Численные методы решения дифференциальных и интегральных уравнений и квадратурные формулы.— М.: Наука. — 1964.
8. Колина Н. П., Солодкин Е. Е. Программа на языке ФОРТРАН для численного интегрирования уравнений пространственного пограничного слоя на линии растекания и на бесконечном скользящем ци-ливдре//Труды ЦАГИ. — 1980. Вып. 2046.
9. Самарский А. А. Теория разностных схем. — М.: Наука. —
1983.
10. Л ы к о в А. В. Теория теплопроводности. — М.: Высшая школа. —
1967.
11. Расчет температурных полей узлов энергетических установок / Под ред. И. Г. Киселева. — Л.: Машиностроение. — 1978.
12. Башкин В. А., Лапина Н. Г. Расчет прогрева цилиндра конечного размера при малых и больших числах Фурье/Друды ЦАГИ. — 1974. Вып. 1582.
13. Чушкин П. И., Шулишнина Н. П. Таблицы сверхзвукового течения около затупленных конусов. — М.: ВЦ АН СССР. — 1961.
14. Ганимедов В. Л. Расчет теплового состояния оболочки в процессе сопряженного теплообмена от сжимаемого ламинарного пограничного слоя // Численные методы решения задач теории упругости и пластичности. Материалы 8-й Всесоюзной конференции, Ужгород, 24—26 мая 1983 г. — Новосибирск. —- 1984.
Рукопись поступша 8/ХІІ1994
