на правах рукописи
Журин Сергей Викторович
Методика численного моделирования конвективного
теплообмена на телах сложной формы с использованием метода
эффективной длины
Специальность: 05.13.18 математическое моделирование, численные методы и комплексы программ
АВТОРЕФЕРАТ
Диссертации на соискание учёной степени кандидата физико-математических наук
Москва 2009 “Ракетно
Работа выполнена в открытом акционерном обществе космическая корпорация “Энергия” имени С.П. Королёва”.
кандидат физико-математических наук,
Научный руководитель:
доцент А.В. Белошицкий.
Официальные оппоненты:
член.-корреспондент РАН, д.ф.-м.н., профессор Егоров Иван Владимирович к.ф.м.н. Аксёнов Андрей Александрович
Ведущая организация: Федеральное государственное унитарное “Центральный предприятие Научно-Исследовательский Институт Машиностроения” (г. Королёв, Московской области).
Защита состоится 10 марта 2010 г. в 15 часов на заседании диссертационного совета Д 212.156.08 при Московском Физико-Техническом институте по адресу: 141700, Московская область, г. Долгопрудный, МФТИ, Институтский пер., 9, Главный корпус, аудитория
С диссертацией можно ознакомиться в читальном зале библиотеки МФТИ Атореферат разослан “” 2010 г.
Учёный секретарь диссертационного совета кандидат физико-математических наук, доцент Коновалов В.П.
I.
ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ
Определение внешних тепловых нагрузок, действующих на поверхность космических возвращаемых аппаратов, является важным этапом решения задачи создания их тепловой защиты и определения температурных режимов конструкции.
В настоящее время существует несколько подходов к оценке конвективного теплообмена сверхзвуковых летательных аппаратов. Одни обладают достаточно хорошей точностью, но требуют большого времени для вычислений. Другие основаны на упрощённых инженерных методиках, требуют малых затрат расчётного времени, но специфика существующих алгоритмов быстрого счёта позволяет оценивать тепловые потоки только на телах достаточно простой формы. Цель данной диссертационной работы - разработать удобный метод оценки конвективных тепловых потоков инженерными методами для тел произвольной формы.
Наиболее простой подход оценки тепловых потоков предполагает определение характеристик пограничного слоя методами локального подобия, т.е. для пограничного слоя на поверхности сложной формы создаётся аналогия с телом наиболее простой формы, например пластины или конуса. При этом для каждого участка сложной поверхности подбираются геометрические параметры простых тел, закономерности развития пограничного слоя на которых известны. Для определения тепловых потоков таким способом необходимо знать распределение газодинамических параметров на внешней границе пограничного слоя по поверхности обтекаемого тела, т.е. решить задачу невязкого внешнего обтекания.
Необходимо отметить, что такой подход правомерен только там, где толщина пограничного слоя достаточно мала, т.е., как правило, на наветренной стороне летательного аппарата. В отрывных зонах такой метод может давать качественно неверный результат, т.к. отрыв потока имеет вязкую природу.
В проектных работах для оперативных инженерных оценок теплообмена наиболее предпочтителен подход использования методов локального подобия. В данной работе для оценок теплового потока применяется метод эффективной длины, разработанный академиком В.С. Авдуевским. В этом методе в качестве геометрического параметра используется длина плоской пластины, пограничный слой на которой имеет те же характеристики, что и в интересуемом месте на поверхности исследуемого тела.
За последнее десятилетие появились на рынке и успешно развиваются программные комплексы с адаптивной прямоугольной расчётной сеткой Adaptive Mesh Refinement.
Типичными представителями этого класса программ являются AeroShape3D (производитель фирма Nika) и FlowVision (производитель фирма ТЕСИС). Основным их преимуществом является быстрое автоматическое построение и адаптация расчётной сетки. Геометрия исследуемого тела может быть практически любой. При этом у пользователя, как правило, нет возможности тонкой настройки сетки в местах предполагаемых особенностей течения – скачков уплотнения, волн разрежения, пограничных слоёв. Это обстоятельство неизбежно сказывается на точности получаемых результатов. Упомянутые выше программы наиболее эффективны при поисковых проектных работах, когда нужно оперативно получить результат для большого числа вариантов с достаточной для этого варианта точностью.
В данной работе для расчётов полей течения, в основном, используется программа AeroShape3D. Алгоритм программы основан на численном интегрировании полных уравнений Навье-Стокса, но возможности вычислительных средств в КБ не позволяют в нужной степени разрешать пограничные слои, так как для этого необходимо огромное количество расчётных ячеек. При недостатке расчётных ячеек для разрешения пограничного слоя, получаемое решение близко к невязкому полю обтекания. Таким образом, AeroShape3D можно с успехом применять для получения распределения газодинамических параметров на поверхности обтекаемого тела и считать их параметрами на внешней границе пограничного слоя.
Объектом исследования является процесс конвективного теплообмена на поверхности гиперзвуковых летательных аппаратов.
Предметом исследования является математическая модель пограничного слоя -метод эффективной длины. При этом используются модели ламинарного и турбулентного пограничного слоя. Вопросы ламинарно-турбулентного перехода в работе не рассматриваются.
Основной целью исследования является разработка метода расчёта теплообмена методом эффективной длины без существенных ограничений на геометрию исследуемого тела.
Актуальность данной диссертационной работы состоит в важной практической потребности в удобном методе быстрых оценок тепловых потоков без существенных ограничений на геометрию исследуемых тел, что необходимо для определения тепловых нагрузок на элементы летательных аппаратов, планирования трубных экспериментов и осмысления их результатов. Невязкие газодинамические параметры при этом с целью снижения трудозатрат целесообразно получать из решения в программных комплексах с прямоугольной адаптивной сеткой.
Научная новизна исследования заключается в применении автором диссертационной работы поверхностной треугольной неструктурированной сетки для оценок тепловых потоков методом эффективной длины, что позволяет исследовать аэродинамический нагрев на телах практически любой геометрии. Разработанная технология позволяет не вводить общую криволинейную систему координат для описания геометрии обтекаемого тела и газодинамических параметров на его поверхности, а получать значения тепловых параметров в каждой ячейке отдельно в своей собственной, независимой системе координат.
Научная новизна подтверждается отличием развитой автором технологии оперативной оценки тепловых потоков с использованием треугольной неструктурированной сетки от прежних, где необходимо “вручную” подстраивать систему координат под особенности геометрии и особенности течения на поверхности, что требует гораздо больших трудозатрат и времени.
Областью применимости разработанного метода оценки тепловых потоков являются задачи обтекания тел при достаточно больших числах Рейнольдса (Re 104…105), в рамках модели тонкого пограничного слоя в областях с малым градиентом давления вдоль линии тока и при отсутствии зон отрыва.
Разработанная методика расчёта тепловых потоков на треугольной неструктурированной сетке имеет практический интерес для организаций и специалистов, занимающихся определением теплового воздействия газового потока на конструкцию технических систем.
Практическая значимость разработанной технологии метода расчёта тепловых потоков в РКК “Энергия” для формирования облика гиперзвуковых аппаратов и тепловых расчётов аэродинамического нагрева космических аппаратов и их элементов.
Предложена концепция формирования нижней поверхности крылатого летательного аппарата с пониженным нагревом кромок крыльев. Ранее в мировой практике ракетостроения подобное решение не применялось. Решение оформлено в виде изобретения и получен патент РФ.
Программы для ЭВМ, созданные автором диссертационной работы, применяются для тепловых расчётов аэродинамического нагрева космических аппаратов и их элементов. В диссертации представлены результаты тепловых расчётов ряда проектируемых и существующих изделий. В результате систематических аэродинамических и тепловых расчётов по разработанной технологии сформирована наветренная поверхность гиперзвукового летательного аппарата с пониженным нагревом кромок несущих аэродинамических поверхностей, которая легла в основу перспективного крылатого космического аппарата Клипер.
Применение разработанной технологии теплового расчёта в программах с автоматически адаптируемой прямоугольной сеткой, позволяет автоматизировать процесс расчёта обтекания и оценки параметров пограничного слоя.
Кроме оценок тепловых потоков, разработанный в диссертации алгоритм можно доработать и применять для оценок коэффициентов поверхностного трения и толщины пограничного слоя.
Предложенная методика может быть дополнена методом среднемассовых величин В.В. Лунёва, для учёта на толщине пограничного слоя неоднородности параметров во внешнем невязком потоке (ударном слое).
При работе над диссертацией автором использованы результаты расчётов газодинамических параметров на поверхности обтекаемого тела в программном комплексе AeroShape3D и расчёты других авторов. В разработанном методе исходные данные на поверхности исследуемых тел могут быть получены любым известным из литературы способом, без привязки к какой-то конкретной программе.
Апробация. Основные результаты диссертации представлены в докладах на конференциях:
• Труды 6-го международного симпозиума по аэрогазодинамике 3-7 ноября года. Версаль, Франция (1 доклад).
• Материалы XIV международной конференции по вычислительной механике и современным прикладным программным системам (ВМСППС-2005). 25-31 мая, • Научная конференциях МФТИ “Современные проблемы фундаментальных и прикладных наук” (7 докладов).
• Конференции молодых специалистов РКК “Энергия” (1 доклад).
Публикации. По теме диссертации опубликовано 15 работ, из них: 4 в отечественных рецензируемых журналах, 10 в материалах Российских и международных конференций, патент РФ на изобретение.
Автор выносит на защиту:
1. Алгоритм и его программную реализацию для расчёта тепловых потоков на треугольной неструктурированной сетке для исследования теплообмена на телах произвольной геометрии.
2. Повышение точности расчёта в окрестности критической точки за счёт использования способа “подсеточного” интегрирования вдоль линии тока внутри треугольной ячейки.
3. Метод расчёта эффективного радиуса внутри треугольной ячейки для определения радиуса эквивалентного тела вращения и оценивания тепловых потоков по методу эффективной длины.
4. Способ профилирования наветренной поверхности гиперзвукового крылатого летательного аппарата позволяющий добиться существенного снижения тепловых потоков к кромкам крыльев.
Структура и объем диссертации. Диссертация изложена на 122 страницах, содержит 91 рисунок, 3 таблицы и состоит из введения, четырёх глав, заключения, одного приложения, списка литературы из 59 наименований, списка использованных сокращений.
II. КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ
Первая глава содержит исходную математическую модель для оценок тепловых потоков методом эффективной длины. Основные соотношения были получены академиком В.С. Авдуевским. В диссертации они обобщены и приведены в виде удобном для использования.
При ламинарном режиме течения в пограничном слое При турбулентном режиме В качестве зависимости вязкости от температуры взята формула Сатерленда Для воздуха она справедлива до Т 2500 K с параметрами µ0 = 1,72e-05 Па·с и С = 110,4 К.
Во второй главе отрабатывается методика расчёта теплообмена на сфере по данным из таблиц Любимова, оценивается погрешность расчёта.
В таблицах Любимова содержаться течения невязкого газа около тупых тел (сфер, цилиндров, параболоидов, гиперболоидов). В расчётах таблиц применялись модель совершенного газа и конечно-разностный метод, основанный на принципе установления по времени.
В диссертации метод расчёта отрабатыватся, оперируя размерными величинами.
Поэтому некоторые графики приводятся также в размерном виде. Используется поле обтекания сферы при М=6, =1.4. Параметры из таблиц приводятся к размерному виду, считая, что параметры набегающего потока: р = 300 Па, Т = 250 К, Mmol = 0.029 кг/моль.
Для дальнейших расчётов считается, что Rсф = 1 м (радиус сферы), Tw = 300 K (температура стенки).
Для расчёта интегралов в формулах (2), (4) используются три различных метода:
метод трапеций, оценка абсолютной погрешности на каждом шаге max f n( x) hn метод Симпсона, оценка абсолютной погрешности на каждом шаге max f n( 4) ( x) ( 2hn ) использование сплайн интерполяции подынтегральной функции, оценка абсолютной погрешности на каждом шаге Jn Jn На рис. 1-2 представлены кривые xef/x для ламинарного и турбулентного пограничного слоя. Три различные кривые соответствуют xef/x, полученным по трём методам интегрирования, x – координата вдоль линии тока.
На оси ординат точкой показаны теоретические значения, которые получаются при x 0. В идеале, при измельчении шага, расчётные кривые должны точно подходить к этим значениям справа. С хорошей точностью мы это можем сказать о кривых полученных по методу Симпсона и с использованием сплайн интерполяции. Метод трапеций на первых шагах даёт большую погрешность.
Для точек метода трапеций показаны планки абсолютной погрешности, оцененные с использованием формулы (6). На графиках они не доходят до кривых Симпсона и сплайн интерполяции, которые в данном случае можно считать точными. Заниженную оценку погрешности интегрирования можно объяснить тем, что вместо максимальных значений вторых производных в выражении (6) брались те, которые получаются в виде разделённых разностей по трём соседним точкам. Как видно из рис. 1-2, такой способ оценки качественно правильно описывает погрешность, но может занижать её в несколько раз (для турбулентного случая примерно в два раза).
У метода Симпсона и сплайн интерполяции абсолютная погрешность достаточно мала и не показана, чтобы не загромождать рис. 1-2.
Рис. 1-2 Отношение xef/x вблизи критической точки для ламинарного и турбулентного На рис. 3-4 приведены распределения абсолютных величин тепловых потоков вдоль образующей сферы. На обоих графиках расчётные кривые не доходят до критической точки растекания, т.к. в ней невозможно рассчитать эффективную длину. Метод эффективной длины здесь не работает. При приближении к точке растекания в формулах (2) и (4) возникает неопределённость 0/0. В ламинарном случае раскрытие её по правилу Лапиталя приводит к конечному значению теплового потока. На рис. 3 оно обозначено точкой. В случае турбулентного течения в пограничном слое тепловой поток в критической точке равен нулю.
На нескольких первых расчётных узлах (рис. 3-4) тепловые потоки, рассчитанные при помощи метода трапеций, заметно отличаются от более точных методов Симпсона и сплайн интерполяции. Далее вниз по потоку кривые для различных методов практически совпадают, что согласуется с оценкой погрешности определения xef на рис. 1-2.
Рис. 3-4 Ламинарные и турбулентные тепловые потоки на сфере.
Влияния неточности входных параметров на результаты расчёта исследуется путём введения возмущения в исходные данные и анализа получаемых результатов. Проводится большое количество расчётов при различных наборах возмущённых входных параметров.
Возмущения входных данных подчиняются нормальному распределению с относительным стандартным отклонением (ОСО) 10%. На рис. 5-6 приведены распределения по образующей сферы отношения ОСО результатов расчёта (тепловых потоков) к ОСО возмущённых входных параметров. Они показывают для каждой точки на сфере во сколько раз относительная погрешность рассчитанных тепловых потоков превышает относительную погрешность входных данных. Для различных методов интегрирования кривые отличаются незначительно. Из графиков на рис. 5-6 можно сделать вывод, что влияние погрешности входных данных для всех методов интегрирования, в целом, усиливается вниз по потоку. Эту зависимость нужно исследовать для каждого течения отдельно.
На рис. 7-8 показаны графики относительных погрешностей результатов расчёта вдоль образующей сферы за счёт погрешности методов интегрирования и в результате неточности задания входных данных. Относительная погрешность входных данных принимается равной 5·10-5.На рис. 7-8 кривая для погрешности метода трапеций лежит гораздо выше оценки влияния входных данных, т.к. этот метод недостаточно использует информацию о распределении газодинамических параметров по поверхности сферы. В этом смысле методы интегрирования Симпсона и с использованием сплайн интерполяции гораздо предпочтительнее. На первых шагах интегрирования их точность примерно совпадает с точностью входных данных. Ниже по потоку ошибка интегрирования понижается, а влияние погрешности входных данных увеличивается, становясь определяющей.
Основной вклад в погрешность определения тепловых потоков по данным из таблиц Любимова вносит погрешность интегрирования методом трапеций. При использовании метода Симпсона или сплайн интерполяции определяющей становится влияние неточности задания входных данных. Машинное округление не вносит заметного вклада на фоне других погрешностей т.к. количество шагов интегрирования невелико.
Рис. 5-6 Распределение отношений ОСО тепловых потоков к ОСО входных данных по образующей сферы (М=6) для ламинарного и турбулентного случаев.
Рис. 7-8 Распределение по сфере относительной погрешности расчёта ламинарных и Третья глава содержит основные положения и алгоритм теплового расчёта на треугольной неструктурированной сетке.
Описывается один из методов триангуляции на основе расчётной сетки AeroShape3D.
На рис. 9 показана треугольная неструктурированная сетка на сфере, полученная по предложенному методу. Критическая точка находится в центре “кольца разряжения” сетки. Для лучшего восприятия сфера немного повёрнута.
Основная идея алгоритма расчёта состоит в том, что из каждого узла треугольной сетки восстанавливается обратная линия тока по направлению к точке растекания.
Начинаясь в исходном узле, она распространяется вверх по потоку, перетекая с одной фасетки на другую, достигая критической точки. Расчет тепловых потоков происходит в обратном направлении по уже пройденному пути, от точки растекания вниз по потоку.
Для ускорения расчёта обратную линию тока можно прерывать в ячейке с рассчитанными параметрами в вершинах, а не доводить её до критической точки.
Начальные данные для интегрирования можно получить здесь линейной интерполяцией.
Предлагается метод вычисления эффективного радиуса на треугольной ячейке.
Параметр Ref (эффективный радиус) служит для учёта увеличения или уменьшения толщины пограничного слоя вследствие стекания или растекания линий тока. Он входит как радиус эквивалентного тела вращения в формулы (2) и (4). Расчёт эффективного радиуса проводится на каждом шаге интегрирования по углу между местными векторами скорости на двух соседних линиях тока.
На рис. 10 показана фасетка с векторами скоростей в вершинах. Отрезок dx – один из шагов интегрирования на линии тока. На одном конце отрезка эффективный радиус известен, обозначен как R, на другом конце он изменится на величину dR, которую требуется найти.
В плоскости фасетки перпендикулярно отрезку dx делается “шаг в сторону”- dz. Длина отрезка dx принимается равной длине отрезка dz. Концы отрезка dz и восстановленные в них вектора скорости задают прямые, которые являются образующими для конуса – эквивалентного тела вращения. Зная длину отрезка dz и угол между векторами d можно определить длину образующей x. Сопоставляя x с шагом dx в подобных треугольниках, получаем:
треугольная сетка на поверхности сферы. радиуса на треугольной сетке.
Описывается способ повышения точности в окрестности точки растекания с помощью подсеточного интегрирования. Участок линии тока внутри одной треугольной ячейки разбивается на несколько шагов. Газодинамические параметры в промежуточных узлах получаем линейной интерполяцией. Измельчая шаг таким образом, можно получить практически любую точность интегрирования, но погрешность результата расчёта будет определяться точностью линейной интерполяции внутри треугольной сетки.
На рис. 11 показаны графики распределения ламинарных тепловых потоков в окрестности точки растекания, полученные различными методами. Тепловые потоки нормированы на максимальное значение в точке растекания.
Кривую, полученную с помощью сплайн интерполяции (рис. 3), считаем точной, к которой нужно стремиться.
Другие две кривые получены интегрированием вдоль линии тока методами трапеций и Симпсона. В методе Симпсона дополнительная точка берётся в середине каждого шага, газодинамические параметры в ней интерполируются линейно.
«