CC BY
52
________УЧЕНЫЕ ЗАПИСКИ ЦАГИ_________________
Том XXV 1994 №3-4
УДК 629. 734. 7. 015. 3
РАСЧЕТ ТРАНСЗВУКОВОГО ОБТЕКАНИЯ СИСТЕМЫ ДВУХ ТЕЛ НА СЕТКАХ РАЗНЫХ ТОПОЛОГИЧЕСКИХ
СТРУКТУР
В. В. Вышинский, С. А. Кравченко, А. М. Сорокин
На примере достаточно сложной аэродинамической задачи расчета трансзвукового обтекания системы двух тел исследуется влияние выбора сеток при реализации одного и того же алгоритма решения краевой задачи для нестационарных уравнений Эйлера на точность результатов, физические особенности и поведение получаемых решений.
Накопление и обобщение опыта создания расчетных методов выдвигают самостоятельную задачу выработки критериев сравнения и оценок используемых сеток, неких априорных рекомендаций, позволяющих заранее, еще на этапах постановки задачи и разработки метода, в соответствии с физическими особенностями явления сформулировать требования, предъявляемые к расчетным сеткам. При этом важно удержаться в разумных пределах точности и не пытаться исправить недостатки используемых моделей и приближенных подходов (например, модели изоэнтропического скачка уплотнения при решении потенциальных задач трансзвукового обтекания или модели турбулентности при решении вязких задач) повышением точности сеточного решения, тем более в рамках идеального газа получить «точные» решения для явлений, где вязкость существенна.
Многолетний опыт создания расчетных методов в различных областях знаний и та роль, какую играют в этих методах расчетные сетки, выделяют методы построения сеток в самостоятельную специфическую область. При этом следует помнить, что решение задачи построения сетки не является единственным: один алгоритм может использовать множество сеток и одна сетка может быть использована для множества алгоритмов. С другой стороны, оценка сетки не может быть выполнена вне рамок использующего ее метода и решаемой с ее помощью краевой задачи.
По типу решаемых в аэродинамике краевых задач имеют свою специфику сетки для расчета потенциальных течений (уравнение для по-
тенциала скорости с условиями непротекания на теле); для вихревых течений (уравнения Эйлера с условиями непротекания на теле), при этом сетки должны обеспечивать достаточное разрешение вне тела для описания крупномасштабных отрывных явлений и вихревых структур; для вязких течений (уравнения Навье—Стокса с условиями прилипания на теле). В последнем случае сетки должны иметь надлежащее, в зависимости от числа Рейнольдса Ие, сгущение вблизи поверхности тела. Последние сетки малопригодны для первых двух случаев, но одни и те же сетки порой используются для расчета как потенциальных, так и вихревых задач.
В рамках одной краевой задачи различие значений параметров течения также накладывает отпечаток на выбор сеток: так входит физическая сущность явления в решение вопроса выбора сеток. Например, чем выше число Яе, тем сильнее должно быть сгущение сетки у поверхности тела в тонкой «вязкой» области. Течения при больших числах Маха набегающего потока Мм > 1, когда перед телом образуется головная волна и область возмущенного течения является полубеско-нечной («гиперболическая» область зависимости решения), требуют расчетных сеток иного вида (например, типа С), чем в случае дои околозвуковых скоростей М„ < 1, когда область зависимости решения носит в целом эллиптический характер.
Существенные и достаточно ясные требования налагают формы обтекаемых тел и областей течения: степень связанности расчетной области, размерность задачи, наличие и отсутствие границ потока (внутренние течения в каналах и внешние задачи обтекания тела неограниченным или ограниченным потоком).
Свои коррективы вносит степень гладкости обтекаемых поверхностей и получаемых решений: наличие или отсутствие в решении разрезов, областей больших градиентов параметров, связанных с физической стороной явления (скачки уплотнения, отрывы потока, выдув струй и т. п.) или с используемой схематизацией (например, наличие разреза скачка потенциала в потенциальных задачах, моделирующих циркуляционное обтекание). Необходимость обеспечить высокую точность задания граничных условий на поверхности тел требует, чтобы сетка была ортогональной вблизи границ, в частности, должна быть связана с поверхностями тел, и гладкой: шаг сетки по каждой из криволинейных координатных линий должен меняться достаточно плавно, в противном случае возникает так называемая аппроксимационная диффузия [1].
При создании сеток должны быть учтены потребные ресурсы ЭВМ (расход памяти, время вычисления). Удачный выбор сетки может повысить скорость сходимости данного алгоритма и/или сократить потребный объем памяти при сохранении точности за счет более удачного расположения расчетных узлов. В работе [2] приведен пример ускорения сходимости обычного релаксационного метода за счет применения многосеточного алгоритма в трехмерном случае. Самостоятельное значение имеет эффективность алгоритма построения сеток. Особенно это важно в методах адаптивных сеток, настраивающихся на поток
в процессе решения, или сеток, настраиваемых на форму поверхности тела при мнопяфатном изменении ее, например, в процессе вариации поверхности тела при решении задач оптимизации формы или задач учета вязкости в приближении пограничного слоя путем наращивания толщины вытеснения.
Отдельное замечание следует сделать по поводу времени в задачах и временных сеток при решении нестационарных задач. Стационарные задачи используют сетки, вообще говоря, отличные от сеток нестационарных задач. В последнем случае, в частности, неприменимы адаптивные сетки. Ситуация может осложниться тем, что в зависимости от параметров течения решение может принимать как стационарный, так и нестационарный (например, квазипериодический) характер. При решении нестационарных задач шаг по времени должен быть постоянен по пространству и обеспечивать единую метрику времени. При решении стационарных задач методами установления шаг по времени из соображений ускорения сходимости (при выполнении условия Куранта—Фрвдрихса—Леви [3]) может быть переменным по пространству при отсутствии единой и постоянной метрики времени.
В данной работе решается внешняя задача трансзвукового (0,8 5 М* <. 1,7) обтекания системы двух тел на основе численного интегрирования уравнений Эйлера для идеального газа. Исследуется влияние топологии расчетной сетки и связанной с ней постановки граничных условий на точность и физическую значимость получаемых решений.
1. Постановка задачи и обсуждение используемой модели. Рассматривается существенно нестационарная, пространственная задача торможения в диапазоне чисел Маха набегающего потока от М^ = 1,7 до 0,8 системы двух тел: спускаемый аппарат (СА) — парашют. При построении математической модели явления делается ряд упрощающих допущений.
Прежде всего, осуществляется переход к стационарной постановке в инерциальной системе отсчета: задача рассматривается в обращенном движении при обтекании системы двух тел стационарным равномерным потоком газа. Не учитываются проницаемость ткани парашюта (этот вопрос исследуется в работе [4]) и динамика его раскрытия, наполнения и последующего изменения формы купола в зависимости от величины перепада давления на нем (см., например, [5]). Пренебрега-ется наличием строп, продольными, поперечными колебаниями и вращением парашюта из-за большого числа степеней свободы, но неста-ционарность, связанная с формированием отрывных структур, пульсацией решения при М„ = const и неподвижном жестком парашюте, при этом учитывается. Такая постановка применяется при моделировании явления в аэродинамической трубе (АДТ), что позволяет использовать результаты опыта для тестирования получаемых решений. В рассматриваемом случае при численном моделировании влияние вязкости газа и стенок трубы не учитывается.
Задача решается в осесимметричной постановке, что исключает неосесимметричносгь наполнения купола и наличие в парашюте перфорации неосесимметричной формы.
Таким образом рассчитывается отрывное обтекание системы двух жестких тел (рис. 1) в рамках нестационарных уравнений Эйлера:
5 + у{рк} = 0, М1 + у{рКг} = -У/>,
ді
■^ + у{(р^ + ^)к} = о, 01
(і)
где р — плотность, р — давление, V — скорость, Е = е + КК/2, в случае идеального газа е = р/р/ (ае -1), ае — показатель адиабаты, с условиями непротекания на поверхности дС1 СА и обеих поверхностях парашюта
(2)
г)
Рис. 1. Внутренние фрагменты расчетных сеток
Форма поверхности парашюта дП1 является внешним параметром задачи и задается заранее в соответствии с формой жесткого парашюта, исследуемого в АДТ или, в случае мягкого парашюта, по фотографиям из опыта. На «бесконечности слева» в соответствии с условиями обращенного движения задается невозмущенный набегающий поток:
где / — орт оси х. На «бесконечности справа» задаются нулевые производные по нормали N к границе расчетной области:
Уравнения движения переписываются в безразмерном виде, где в качестве характерных величин выбраны параметры набегающего потока рх, р.*,, и диаметр купола парашюта 2). Посредством интегри-
ностях парашюта определяются величины сил, которые при определении коэффициентов сопротивления СХСА И схп Делятся на скоростной
расчета определяется также величина донного сопротивления СА сх д и
парашюта схд. При этом интегрирование давления в первом случае
производится по донному срезу СА, во втором — по внешней поверхности купола парашюта.
В работе используется схема первого порядка точности. Выбор ее обусловлен сложностью решаемой задачи. Из-за сильного взаимодействия ударных волн и отрывов потока при сделанных допущениях ее решение является весьма приближенным, хотя и позволяет выявить основные свойства течения. Повышение порядка аппроксимации схемы не приводит к существенному уточнению получаемых результатов, но значительно усложняет метод и увеличивает время вычисления. Кроме того, разностная аппроксимация на шаблонах с меньшим числом узлов уменьшает вредное влияние, обусловленное резким изменением шага сетки (см. [1]). Использование схемы первого порядка позволяет наряду с сетками типа О, для которых поверхность обтекаемых тел является одним из координатных уровней, использовать планшетные сетки типа Н, топологическая структура которых аналогична координатным линиям декартовой системы координат. Причина обращения к подобным сеткам заключается не только в простоте их построения, но и в том, что задание граничных условий на ступенчатой линии, ап-
Р\х-*~во Ао>
Р|х-»-оо —*'Р«в> |х->-00 ~* * Ко >
(3)
(4)
рования распределений давления на поверхности СА и обеих поверх-
напор дп = /2 и характерную площадь, равную для СА тиР/л, где
й — диаметр донного среза СА, для парашюта — я2)2/4. В результате
проксимирукяцей поверхность тела, приводит к возникновению дополнительной завихренности в потоке, которая может оказать существенное влияние на физический характер получаемого решения (см., например, [7]).
2. Алгоритм построения расчетных сеток. Координаты узлов двумерных регулярных сеток типа О находятся как решение задачи о минимизации дискретного функционала. Для этого осуществляется отображение некоторой параметрической области (? простой формы с координатами т| на исходную область П с координатами х, г. При этом отображении границы областей соответствуют друг другу, а искомая сетка получается как образ равномерной декартовой сетки, заданной естественным образом в области (?. Если такое отображение конформно, то оно реализует минимум функционала:
2 2 2 2 *Х7 +К + ДГ + Г
г т -Г и, т I
/-1—1Ц-2—1(5)
1 в
где I = х^гц- хцг^ — якобиан преобразования. Поэтому координаты
узлов сетки определяются исходя из минимума конечно-разностной аппроксимации функционала (5), предложенной в работе [6]:
4 + *1 + г!
V V («)
»7 к=1 •'
При аппроксимации параметрическая область (г разбивается на треугольники, и для каждого треугольного элемента подынтегральное выражение в формуле (5) аппроксимируется отношением суммы квадратов длин сторон треугольника, выходящих из данного узла: х| +
и *2 + /^ к площади этого треугольника /Л. Величины х|,
являются фактически конечно-разностными аппроксимациями соответствующих производных. Каждому узлу (/, у) соответствуют четыре
треугольных элемента, у которых из этого узла по координатным линиям выходят две стороны.
В данной работе при численной реализации функционал (6) трансформируется к виду
к ’
где,функции р = р{%, тО и 9 = д(§, т]) служат для управления сгущением сетки, а вместо площади треугольного элемента используется величина /Л - . В качестве е,у выбирается наименьшее значение площади
треугольных элементов, имеющих данный узел своей вершиной. Если
все площади положительны, то еу = 0. Такая модификация функционала позволяет работать с самопересекающимися сетками в качестве начального приближения.
Минимизация функционала осуществляется методом градиентного спуска при последовательном обходе узлов. При этом в сумме (7) каждый раз присутствует лишь 12 слагаемых по числу треугольных ячеек, имеющих данный узел в качестве вершины. Положение узлов на внешней границе области и на поверхности парашюта определяется из минимизации этого функционала.
На рис. 1 в разреженном виде приведены внутренние фрагменты трех сеток типа О, используемых в расчетах. Все они имеют 100 лучей по 50 узлов в каждом и одинаковую топологию у границы расчетной области. На внешних границах расчетной области линии уровня £ = const и т| = const ортогональны между собой.
Отличительной особенностью первой сетки, используемой для расчета обтекания изолированного парашюта (L -> оо, где L — расстояние между донным срезом СА и входным сечением парашюта), является ортогональность линий уровня у оси симметрии г = 0.
Единая сетка для расчета обтекания системы двух тел: СА — парашют при расстоянии между ними L = L/D = 1 в отличие от предыдущей у оси симметрии имеет неортогональные линии уровня (скошенные ячейки). Они неортогональны также у поверхности С А (конус-цилиндрического тела), особенно в донной его Части, и, как и в первом случае, у поверхности парашюта, в наибольшей степени вблизи входного сечения.
Третья сетка, используемая для расчета обтекания системы двух тел при L =1,75, лишена двух из указанных выше недостатков: она орто-гонализована у поверхности донного среза СА и у оси симметрии.
Обе последние сетки используются для расчета обтекания изолированного СА. В этом случае разрез с условиями непротекания, моделирующий в расчете парашют, в поле течения игнорируется. Размеры расчетной области по оси х составляют около 20, а по оси г — 10 длин СА
Наряду с регулярными сетками типа О в расчетах использовались нерегулярные планшетные (цекартовые ортогональные) сетки типа Н в цилиндрической системе координат с осью х, связанной с осью симметрии задачи (см. рис. 1, г). Все поле течения при этом делится на две расчетные области: область СА и область парашюта. В каждой области используется 6 блоков, в пределах которых сетка имеет постоянный шаг по пространству (различный по направлениям х и г), но в месте сшивания блоков величина шага сетки по направлению вдоль линии сшивания одинакова для обоих смежных блоков. Общее число узлов в сетке для системы двух тел равно 2x50x50, для изолированных СА — 50x50 и парашюта — 50x60. Таким достаточно простым путем удается построить единую сетку, имеющую при удалении от тел достаточно большие шаги по пространству, что позволяет реализовать краевые условия (3) и (4) на конечном (порядка 10D) расстоянии в рамках единой задачи.
3. Метод решения краевой задачи. Решение задачи на обеих сетках ведется в осесимметричном приближении в цилиндрической системе координат х, г, 0 с осью х, связанной с осью симметрии течения. При реализации алгоритма расчет выполняется для одного продольного меридионального сечения 0 = const области течения. Используется метод [7], в основу которого положена схема конечных объемов первого порядка точности по пространству и времени [8]. В расчетной схеме отсутствует искусственная вязкость, так как используемый метод с разностями против потока является устойчивым без введения дополнительных диссипативных членов.
Наличие схемной вязкости неявным образом модифицирует исходные уравнения Эйлера. Эти изменения можно смоделировать добавлением в правые части последних двух уравнений (1) некоторых членов, которые по аналогии с физической вязкостью можно записать в виде тензора (см., например, [9]), так что в численной схеме решается уже иная задача. Схемная вязкость порождает численную диффузию, которая подавляет короткие волны в решении и приводит к уширению любого разрыва (скачка уплотнения, слоя смешения) на величину 66ль-шую, чем шаг сетки. Длинные волны — а сеточные методы изучают в основном длинноволновое приближение — минимально подвержены воздействию численной диффузии, что позволяет рассматривать решения с ударными волнами, внутри которых (теперь уже за счет схемной вязкости) происходит переход кинетической энергии коротких волн в тепловую энергию. Точно так же происходят диссипация энергии коротких волн и порождение завихренности на слое смешения, разделяющем отрывную область циркуляционного движения от внешнего потока.
Как показано в работах [7], [8] и [10], отрывное обтекание с зонами возвратно-циркуляционного течения может быть успешно смоделировано в рамках идеального газа при решении нестационарных уравнений Эйлера, содержащих в себе механизм переноса завихренности. При этом наряду с физическим механизмом порождения завихренности на искривленных скачках уплотнения с переменной по фронту энтропией в численной схеме имеется нефизический механизм порождения завихренности за счет схемной вязкости. Последняя обусловлена аппроксимацией уравнений движения, а также граничных условий, особенно при использовании сетки типа Н, где условие непротекания (2) накладывается на аппроксимирующей поверхности тела ступенчатой ломаной, совпадающей с ближайшими к нему линиями уровня сетки X,- = const, Гу = const. При этом в каждой угловой точке ступеньки 0> у) v \r=rj= 0» и\х=х, = 0 в малом моделируется условие прилипания
(ujj = Vjj = о), в окрестности нее образуется зона смешения и происходит интенсивное вихреобразование.
Расчет выполнен на сравнительно редких сетках, особенно если учесть необходимость обеспечить высокую степень разрешения решения в окрестности СА и парашюта и, в частности, между ними. Последнее обстоятельство становится определяющим в выборе сеток при моделировании изменения расстояния L между СА и парашютом, что
используется на практике для изменения сопротивления системы двух тел при так называемой «аэродинамической рифовке» парашюта.
В данном методе акцент был сделан на выбор сеток, обеспечивающих учет влияния отрывного обтекания первого тела (СА) на отрывное обтекание второго тела — жесткого парашюта, чего в работах [4] и [5] не делалось.
Краевые условия на бесконечности (3) и (4) для всех используемых сеток ставятся на границе расчетной области, что, в отличие от работы [7], ограничивает снизу число Маха набегающего потока М*,. Такой подход является традиционным (см., например, [8]). Вопрос о достаточности размеров расчетной области при решении трансзвуковой задачи в рамках уравнений Эйлера исследуется в работе [11]. В частности, в двумерном циркуляционном случае достаточным признано расстояние порядка 50 хорд профиля. В бесциркуляционном пространственном случае это расстояние может быть существенно уменьшено из-за более быстрого затухания возмущений при удалении от тела. В работе [7] вопросы, связанные с асимптотикой решения на бесконечности, решаются во внешней по отношению к внутренней вихревой области потенциальной задаче с помощью хорошо отработанного конечно-разностного метода для уравнения относительно потенциала скорости на сетке типа О с отображением внешности нулевой линии тока, ограничивающей обтекаемое тело и замкнутую отрывную область за ним, на внутренность круга, а бесконечности — в его центр. Использование этого подхода позволяет корректно выбрать границу расчетной области (порядка 10 характерных размеров задачи) при постановке на ней граничных условий.
Метод реализован на ЭВМ с быстродействием 0,25 мегафлопс. На сетке типа О используется переменный по пространству безразмерный шаг по времени т = А/Кте/2):
где т0 = 0,0013, а а{х,г) — площадь текущей ячейки. Для выполнения N = 1000 шагов по времени требуется около одного часа машинного времени.
В случае применения сетки типа Н с числом узлов 2x50x50 с минимальным шагом сетки А = й/2) = 0,02 используется постоянный шаг по времени т0 = 0,01, так что число Куранта Си = (|К| + тоА > гДе
а — скорость звука меньше единицы. При этом время выполнения N = 1000 итераций около 40 мин.
Таким образом, из-за сильного сгущения сетки в первом случае условие устойчивости [3] Си < 1 на порядок по сравнению с планшетной сеткой уменьшает величину т0 и, следовательно, существенно увеличивает потребное число шагов по времени N для достижения сходимости решения. Использование переменного по пространству времен-
(8)
нбго шага (8) позволяет несколько повысить эффективность метода — при этом ощутимое увеличение шага т происходит лишь на больших расстояниях от тел,— но это ограничивает применимость метода случаем построения стационарных решений. Анализу сходимости метода и исследованию точности получаемых результатов посвящены работы [7] и [12].
4. Результаты расчетов. В качестве первого примера на рис. 2 приведены результаты расчета обтекания изолированного СА — конус-цилиндрического тела с углом раствора конуса ср = 150°, удлинением носовой части Хи = Ьн/й = 0,12 и удлинением цилиндрического участка Х.ц = Ьп/ё = 2,5, где 1^ — длина конического, а Хц — цилиндрического участков СА Цилиндрический участок заканчивается донным уступом, на изломе которого моделируется условие прилипания и естественным образом фиксируется положение отрыва. При этом от-
Рис. 2. Результаты расчета обтекания изолированного спускаемого аппарата на сетках типа О — фрагменты а — в и типа
Н— фрагмент г
рыв возникает из-за инерции жидкой частицы и наличия схемной вязкости, не позволяющей сразу, «вовремя» реализовать условие непроте-кания (2) на поверхности тела. Образовавшаяся каверна в силу сплошности среды заполняется газом, образующим зону циркуляционного течения. Физическим источником завихренности ю = rot F в соответствии с теоремой Крокко
V х 5 = -TgradS'
являются скачки уплотнения с искривленным фронтом (переменной вдоль фронта энтропией S). Таковыми являются головные ударные волны при Мте > 1 и внутренние скачки уплотнения при наличии в течении местных сверхзвуковых областей.
На рис. 2 представлены результаты расчета на сетке типа О с косыми ячейками у оси симметрии (см. рис. 1, б) при числах Мм =0,8; 1,0; 1,7 и на сетке типа Н при Мм = 1. Наряду с полуконтуром тела представлены линии тока и изомахи М = const, приведены значения донного cX(j и полного схСА сопротивлений. Величина последнего монотонно растет с ростом Мк, при этом область влияния перестраивается из эллиптической в гиперболическую, и при Мте > 1 могла бы быть использована более экономная сетка, но с ограниченной сверхзвуковыми скоростями набегающего потока областью применения.
При Мю =0,8 (рис. 2, а) источником завихренности является только схемная вязкость. Наряду с отрывной зоной в кормовой области наблюдается длинная отрывная зона на боковой цилиндрической поверхности СА. С ростом числа Мда и появлением в потоке физической завихренности последняя исчезает, а отрывная зона в кормовой области теряет свою выпуклость, т. е. наличие отрицательной (по часовой стрелке) завихренности в потоке облегчает обтекание тупого и прямого углов, уменьшая инерцию жидкой частицы и способствуя реализации условия непротекания.
Результаты расчета на сетке типа Н (при этом контур полутела СА аппроксимируется 18 ячейками) представлены на рис. 2, г. В поле течения едва различима отрывная зона на цилиндрической поверхности, что объясняется дополнительным внесением завихренности в поток за счет схемной вязкости при ступенчатой аппроксимации конической поверхности головной части.
Расчет обтекания изолированного СА при Мда = 1 проведен также на сетке типа О с ортогонализованными ячейками (рис. 1, в). При этом получено cx<i = 0,21, схса = 1,33. Последние результаты ближе к расчету на сетке типа Н как по значениям интегральных характеристик сХ(], схСА> так и по форме линий уровня М = const и нулевой линии тока
в кормовой области, чем те же результаты для сетки типа О с косыми ячейками. Отличительной особенностью последнего расчета является смыкание через узкий перешеек области отрыва на боковой поверхности с кормовым отрывом. В дальнейшем эта сетка используется для расчета обтекания системы двух тел.
На рис. 3 приведены фотографии, полученные методом Теплера в аэродинамической трубе при испытаниях с числами Мда = 0,86 (рис. 3, а) и 1,1 (рис. 3, б). Число Рейнольдса по диаметру тела в обоих случаях приблизительно одинаково Яе^ = 0,5 • 106. Как видно, и в опыте наблюдается ослабление отрыва на боковой поверхности с ростом Мж.
Результаты расчета хорошо согласуются с экспериментальными данными по донному сопротивлению cxj (см., например, [12]) и сопротивлению носовой части схн = cxqa ~ cxd работы [13], где при числах М*, = 1, Rerf = 1,8 106 получено схн = 1,0.
В качестве второго примера на рис. 4 приведены результата расчета обтекания изолированного парашюта при числе Маха набегающего потока Мж =1,05. Форма купола задана по фотографиям из опыта в
АДТ при Rej3=3-106 и Мда = 1,5 в предположении, что с уменьшением Мда она меняется мало. Расчет выполнен на двух сетках типа О (рис. 4, а) и типа Н (рис. 4, б). Наряду с формой купола, линиями тока, изомахами М = const приведены результаты сходимости решения по времени t - т0/ для значений полного сопротивления схп
(схп -> 1,83 на первой сетке и схп -> 1,50 — на второй) и в случае расчета на первой сетке — донного сопротивления cxj) ->0,57, а также значений давления p/pq для четырех указанных на рис. 4,а точек — двух внутри купола (77, ТЗ) и двух снаружи (Т2, Т4). Как видно, расчеты на планшетной сетке (при этом полуконтур парашюта аппроксимируется 35 ближайшими ячейками) сходятся быстрее: требуется всего
Рис. 3. Результаты испытаний в АДТ-112 ЦАГИ при числе Рейнольдса Rerf = 0,5 х Ю6
Рис 4. Результаты расчета трансзвукового обтекания изолированного парашюта на сетках типа О — фрагмент а и типа Н — фрагмент б
1 — 2 тысячи итераций, в то время как на сетке типа О с близким числом узлов за счет меньшего т0 требуется уже 5—10 тысяч итераций. При этом средняя по всему полю невязка для плотности
уменьшается на три порядка.
Изолированный парашют вносит очень большие возмущения в поток, что на порядок по сравнению с СА замедляет сходимость решения. В случае изолированного СА при расчете на сетке типа О было достаточно 100 — 200 итераций при М„ > 1 и до 1000 при Мда =0,8, причем сопротивление носовой части схн сходилось почти в два раза быстрее, чем Схс[.
Сравнение результатов расчета обтекания парашюта на обеих сетках показывает, что, как и в первом примере на сетке типа Н, решение дает меньшие величины сопротивлений и размер отрывной зоны, что связано с дополнительной завихренностью, вносимой в поток за счет ступенчатой аппроксимации поверхности купола на планшетной сетке.
Полученные результаты удовлетворительно согласуются с экспериментами и расчетами по определению коэффициента сопротивления изолированного непроницаемого жесткого парашюта работы [5], где в диапазоне чисел Маха 1,15 Мв ^ 3,0 найдена величина
ь м
схп = 1*5 -1,6, а с увеличением числа Мв от 0,6 до 0,9 опыты выявляют монотонный рост схп от 1,3 до 1,45. Опыты с изолированным
жестким парашютом той же формы, выполненные с целью дальнейшей проверки результатов расчета в АДГ-112 ЦАГИ при числах Мм = 1,05 + 1,50 и 1^0 = 1,1 • 106, дали значения схп = 1,33 т-1,43, причем максимум сопротивления приходится на Мм =1,15. Эти результаты хорошо согласуются с расчетом на селю типа Н: расчет завышает величину схп на 2 — 996. Расчеты на сетке типа О согласуются хуже, они на
35 — 40% превышают величину сопротивления, найденную из опыта. При этом следует помнить, что эксперименты проводились с кормовой державкой, наличие которой не моделировалось в расчетах на обеих сетках.
На рис. 5 приведены результаты расчета сверхзвукового обтекания с числом Мда =1,5 системы двух тел СА — парашют при расстоянии
Рис. 5. Результаты расчета сверхзвукового обтекания системы двух тел: спускаемый аппарат — парашют на сетках типа О — фрагмент а и типа
Н— фрагмент б
между ними L = L/D = 1,75. Соотношение диаметров парашюта и СА D/d = 6. Наличие переднего тела существенно изменяет картину обтекания. На рис. 5, а приведены результаты расчета на сетке типа О, на рис. 5, б — на сетке типа Н. В данном случае система двух тел обтекается по так называемой схеме с открытым следом, когда СА и парашют образуют единую отрывную зону, начинающуюся на первом теле и заканчивающуюся за вторым, так что сопротивление парашюта (при расчете на сетке типа О) схп = 0,58 в данном случае существенно ниже,
чем в случае L -> оо, и определяется в основном донным сопротивлением cxj) =0,43. Для сравнения: при расчетах на той же сетке изолированный парашют имеет схп =1,76, схд = 0,39.
Донное сопротивление СА в этом случае отрицательное cxd = ~ 0,21, в этой части реализуется избыток давления
ср =(р- Роо)/я<х, > 0- Сопротивление носовой части схн = 1,33 от присутствия парашюта меняется мало. Для изолированного СА расчеты на той же сетке при Мм =1,5 дают cX(j =0,25, среднее значение донного
давления Cpj =-0,25, схн =1,31, т. е. величина схн СА практически
не меняется при выпуске парашюта, а величина донного сопротивления сильно изменяется, вплоть до смены знака при приближении парашюта к СА, что связано со сменой типа обтекания от схемы с закрытым следом, когда за СА образуется самостоятельная отрывная зона, как в случае изолированного СА, и сХ(/ > 0, к схеме с открытым следом, как в данном случае, когда в общей с парашютом внутренней отрывной зоне реализуется < 0.
На рис. 5, а приведены линии тока, помеченные стрелками, и изомахи М = const, а также полуконтуры продольных сечений СА и парашюта и распределения давления вдоль них ср(»у), где 0 £ S < 1 —
длина дуги. Приведены также результаты сходимости решения по числу итераций J для полного сопротивления СА схса, донного cxd
и полного схп сопротивлений парашюта. Как видно, установление решения происходит при N = 7*8 тысяч итераций (7 = 9* 10).
Аналогичный расчет на сетке типа Н (см. рис. 5, б) при той же схеме течения с открытым следом дает меньшую величину отрывной области за парашютом и несколько меньшую величину его полного сопротивления (чем и объясняется лучшее согласование с опытом, где парашют установлен на кормовой державке, наличие которой уменьшает величину донного сопротивления) схп =0,50. Значения полного
Сд-са = 1)07 и донного cxj =-0,25 сопротивлений СА также достаточно хорошо согласуются с результатами, полученными на сетке типа О, но установление решения по числу итераций в последнем случае происходит значительно быстрее — при N = 1000 (/=10). Распределение давления по донному срезу СА ср(г) (см. рис. 5,6) соответствует участку S = 0,85 * 1,0 эпюры с^са^ ) на рис. 5, а.
Результаты расчета на обеих сетках хорошо согласуются между собой по схса и схп> а также с экспериментальными результатами для
мягкого парашюта, полученными в АДТ-109 ЦАГИ при Яед = 3 -106, Мю =1,5. Хорошее согласование результатов расчета между собой наблюдается для схСА при М* = 1,3 и для схп при М„ = 1. При остальных значениях Мда расчет на сетке типа О дает несколько ббльшие величины сопротивлений. Опыты в АДТ-109 с жестким парашютом при меньшем, чем в расчете, 1)Д/ = 4,6 дают при Мж = 1,5 существенно ббльшую величину сх п =0,92. При меньших значениях Мда величина схп, найденная из опыта, лежит между расчетными результатами, полученными на сетках типа О и типа Н. При Мда = 1 с точностью до 10% они совпадают.
Расчеты обтекания системы двух тел при меньших числах Маха набегающего потока, выполненные на сетке типа О, демонстрируют нестационарный квазипериодический характер, где в решении наряду с_ быстропериодическими колебаниями с безразмерным периодом Т = Т - = 0,2 имеются медленнопериодические колебания с од-
ной или двумя частотами. Например, при Моо=1,0 величина полного сопротивления носит периодический характер с периодом Т = АЖ-то = 3,4 и амплитудой колебаний от гшпсх:п=0,50 до шахсхп = 1,30 при почти постоянной реличине донного сопротивления схи, что связано с пульсациями внутренней отрывной зоны. При этом полное сопротивление СА сх са носит характер модуляции короткопериодического изменения с малой амплитудой длиннопериодическим — с большей амплитудой. Соответствующий характер носят и зависимости от времени значений давления р/р<) в донной области СА и особенно внутри купола парашюта, причем их колебания, как коротко-, так и длиннопериодические, происходят в противофазе. В носовой части СА и в донной области парашюта давление почти не меняется.
Как и в работе [7], где получено нестационарное решение для обтекания сферы звуковым потоком, колебания исследуемых параметров происходят примерно в одни и те же моменты времени, не зависящие от измельчения сетки и изменения шага т0. При сравнении с опытными данными и стационарными решениями, полученными на сетке типа Н, используются осредненные по периоду медленных колебаний величины сопротивлений.
Длиннопериодическое движение связано с размыканием и смыканием в единую отрывную область зон отрыва с тупого угла носовой части и прямого угла донного среза и «дыханием» внутренней отрывной зоны, причем минимум сха наблюдается, когда внутренняя отрывная область вспучена, максимум,— когда она сужена. При этом купол парашюта на фиксированном расстоянии играет роль резонатора, и период колебаний изменяется в зависимости от расстояния между телами и скорости набегающего потока, так, при Ь = 1 он ниже и при Мд, = 1 равен Т =3,3.
Последние результаты интересны тем, что демонстрируют возможность получения нестационарного квазипериодического решения при стационарной постановке краевой задачи, но при нестационарном алгоритме получения решения. Нестационарность в решении может появиться и при уменьшении числа Маха набегающего потока М*,, как при L =1,75, и при уменьшении расстояния между телами: например, решения при L = 1 для всех скоростей в диапазоне 0,8 £ М„ 21,7 носят нестационарный характер. Очевидно, что нестационарный характер решения исключает возможность применения адаптивных сеток, хотя, как показывают результаты работы [14], точность расчета на адаптивных сетках сложных течений, содержащих поверхности разрыва, заметно выше, чем на обычных, фиксированных в процессе расчета сетках.
Выбор сетки обусловлен не только используемым алгоритмом, который является критерием эффективности применяемой сетки, но и решаемой краевой задачей, при этом математическая модель ставит меру и определяет 1раницу возможного уточнения решения, а также физической стороной явления, ибо результаты эксперимента используются в качестве окончательного критерия точности метода.
В данной работе неединственность выбора сеток при решении аэродинамических задач сохранена до самого конца, что позволяет сделать вывод, основываясь на сравнении полученных решений с опытом.
Так, например, при Всех своих недостатках планшетная сетка имеет и ряд преимуществ. Наряду с простотой построения она позволяет легко и с контролируемой точностью изменять расстояние L между телами путем добавления нужного числа уровней х = const с тем же шагом, а также элементарным образом учесть наличие границ потока (стенок АДТ) при моделировании решения поставленной задачи в условиях трубного эксперимента. В используемой сетке в пределах блоков шаг постоянен по обоим направлениям, что исключает появление аппроксимационной диффузии (см. [1]), однако на стыках блоков имеется скачок шага сетки в направлении, перпендикулярном границе раздела блоков, что неминуемо приводит к ошибкам аппроксимации.
К достоинствам используемых сеток типа О следует отнести лучшую аппроксимацию уравнений и граничных условий: при том же числе узлов большую охватываемую сеткой расчетную область и меньшую схемную вязкость вблизи тел, что приводит к повышению разрешающей способности метода во внутренней вихревой области, хотя и за счет замедления сходимости решения.
Имеются две причины различия решений для одного алгоритма на двух разных сетках с одинаковым числом узлов. Первая состоит в топологическом различии сеток типа О и Ни связанного с этим различием уровня аппроксимации решения. Так, на сетке типа О из-за меньшей схемной вязкости реализуются нестационарные решения. Нестационарность проявляется в вице взаимодействия отрывов на боковой поверхности с кормовым отрывом. Вторая причина состоит во вне-
схемном различии, привнесенном извне в рассматриваемую задачу. На сетке типа О гладкие участки поверхностей тел аппроксимируются гладким образом — эта сетка лучше соответствует решаемой краевой задаче. Ступенчатая аппроксимация криволинейных участков поверхностей и участков, не совпадающих с сеточными линиями, на сетке типа Н приводит из-за реализации в малом условия прилипания к дополнительному порождению завихренности на изломах ступенек, что уменьшает размеры отрывной области за парашютом и величину донного сопротивления. В ту же сторону действует наличие кормовой державки в опыте, чем объясняется лучшее согласование с ним расчетов на сетке типа Н.
Сравнение с экспериментальными результатами показывает, что использование сетки типа О позволяет хорошо воспроизвести физические особенности явления: в опыте визуально наблюдается отрыв с тупого угла, получаемый в расчетах на сетке типа О и едва различимый на сетке типа Н с тем же числом узлов. Внесение первым телом в поток дополнительной завихренности (при расчете на сетке типа Н еще и нефизической — за счет ступенчатой аппроксимации конической поверхности головной части СА) уменьшает сопротивление парашюта. Это же наблюдается и в опыте, где в случае остроносого СА(ср = 22°) схп больше, чем у тупоносого (ср = 150°) при всех значениях М„ и L, за исключением малых L = (2- 3)d, при сверхзвуковых скоростях (М. > 1,1), когда остроносый СА за счет Х-скачков на его цилиндрической поверхности и отрывов потока из-под них вносит в поток большую завихренность, расширяя след за СА.
Авторы выражают благодарность А. Г. Ерезе и В. С. Маркину за любезно предоставленные ими результаты экспериментальных исследований.
ЛИТЕРАТУРА
1. Т о м п с о н Дх. Ф. Методы расчета сеток в вычислительной аэродинамике // Аэрокосмическая техника.— 1985, № 8.
2. Ш м и л о в и ч А., К о ф и Д. А. Многосеточный метод расчета трансзвуковых течений около комбинации крыло — фюзеляж — хвостовое оперение // Аэрокосмическая техника.— 1986, № 2.
3. Courant R., Friedrichs К. О., Lewy Н. Ober die paitiellen Differenzengleichungen der mathematischen Physik // Math. Annalen.—
1928, Bd. 100.
4. Гувернюк С. В., Лоханский Я. К., Савинов К. Г. О численном исследовании аэродинамики тонкостенных проницаемых э!фанов // Парашюты и проницаемые тела,— М.: Изд-во МГУ, 1987.
5. Б а к ш е е в С. В. Нестационарная аэроупругая модель наполнения проницаемого осесимметричного парашюта в потоке идеального сжимаемого газа // Киев: КНИГА, 1988.
6. Иваненко С. А., Чарахчьян А. А. Алгоритм построения криволинейных сеток из выпуклых четырехугольников // ДАН СССР.— 1987.
Т. 295, N° 2.
7. Вышинский В. В., Кравченко С. А. Расчет отрывного осесимметричного обтекания тел на основе решения уравнений Эйлера во внутренней вихревой области // Труды ЦАГИ.— 1990. Вып. 2494.
8. Белоцерковский О. М., Давыдов Ю. М. Метод крупных частиц в газовой динамике // М.: Наука, 1982.
9. П о т т е р Д. Вычислительные методы в физике // М.: Мир, 1975.
10. Белоцерковский О. М. Численное моделирование в механике сплошных сред // М.: Наука, 1984.
11. Т h о m a s J. L., S а 1 a s М. D. Far-field boundary conditions for transonic lifting solutions to the Euler equations // AIAA Paper 85-0020.— 1985.
ИВышинский В. В., Кравченко С. А. Расчет пространственных отрывных течений в рамках нестационарных уравнений Эйлера // Ученые записки ЦАГИ,— 1991. Т. 22, № 2.
13. Петров К.П. Аэродинамика элементов летательных аппаратов // М.: Машиностроение, 1985.
14. Н а к а х а с и К., Д е й у э р т Д ж. С. Новый метод построения трехмерных адаптирующихся сеток // Аэрокосмическая техника.— 1987, № 1.
Рукопись поступила 11/V1993 г.
