Comparison of the MUSCL-type schemes for a gas flow calculation in de Laval nozzles.pdf При моделировании течений продуктов сгорания в газодинамических трактах ракетных двигателей широко применяются численные методы повышенного порядка точности, основанные на схеме Годунова [1-3]. Методы такого типа развиваются уже на протяжении более сорока лет, и существует огромное множество их вариантов и реализаций [4]. Впервые схема второго порядка точности по пространству, основанная на методе Годунова была предложена Колганом в 1972 году [5, 6]. Основная идея схемы Годунова-Колгана заключалась в замене кусочно-постоянного распределения функции внутри ячейки на кусочно-линейное распределение. Для удовлетворения условиям сохранения полной вариации (Total Variation Diminishing, TVD) Колган предложил использовать принцип минимального значения производной. В 1979 году подход, предложенный Колганом, получил свое развитие в работе Ван Лира [7, 8], в которой была предложена схема MUSCL (Monotonic Upstream-Centered Scheme for Conservation Laws). Схема приобрела большую популярность, и в результате её развития родилось большое семейство схем типа MUSCL [9]. Схемы типа MUSCL хорошо протестированы на течениях с разрывами (задачи: «Сода», «Лакса», «Mach 3», «Supersonic shock tube», «1-2-3», «High Mach», «Сильная ударная волна», «Подвижный контактный разрыв») [10]. Однако для решения инженерных задач газовой динамики ракетных двигателей, помимо определения структуры течения, важнейшей задачей является определение интегральных характеристик. В связи с этим требуется проведение исследований и определение комбинаций базовых схем, функций ограничителей и методов определения потоков через грани расчетных ячеек, позволяющих с высокой точностью определять энерго-тяговые характеристики двигателей. Математическая постановка Система уравнений Эйлера, описывающая нестационарное течение сжимаемого невязкого газа в одномерном приближении имеет вид [11, 12] dQ+MQL0, (1) dt dx где Q - вектор консервативных переменных, F (Q) - вектор потоков: Р Q = (2) ( Р \ ( pu pu 2 + p pu ; F (Q) = v(e+p )u Здесь p - плотность; t - время; p - давление; u - скорость движения газа; e = pe + pu2/2 - полная энергия единицы объема; е - удельная внутренняя энергия. Система замыкается уравнением состояния идеального газа: е = p/(y-1)p, где у - показатель адиабаты. Схема «Годунова» Для системы уравнений (1) схема Годунова имеет вид [12] Qi - Qi + F+1/2 - Fi-1/2 = 0 (3) Дt Дх где Дt - шаг по времени, Дх - шаг по пространству; Q1 и Qn+l - значения вектора консервативных переменных на текущем и следующем шаге по времени соответственно, F+1/2 и F-1/2 - потоки через правую и левую грани расчетной ячейки. Для расчета потоков может быть использовано точное решение задачи о распаде произвольного разрыва (определение потоков по методу Годунова) [12] либо другие методы, основанные на приближенных решениях задачи Римана: Рое, Ошера, HLL, HLLC, HLLE и другие [10]. Таким образом реализация шага по времени в схеме Годунова проводится в два этапа: Этап 1 - «Решение задачи Римана». На данном этапе по известным значениям Qin , с использованием точного или приближенного решения задачи Римана о распаде произвольного разрыва, определяются потоки FMn и Fi-1/2 через грани расчетной ячейки: F+1/2 = RimmanSolverQn, Q+,). (4) Этап 2 - «Вычисление параметров на новом временном слое». На этом этапе по известным значениям Qi, F+1/2 и F-1/2 с использованием (3) определяются параметры газа Qin+1 на новом временном слое. Схема «MUSCL» Для повышения порядка точности по пространству часто используется техника квазимонотонной интерполяции сеточных решений (MUSCL) [8, 13, 14], удовлетворяющая условию TVD. Основная идея схем MUSCL состоит в повышения точности схемы за счет изменения порядка интерполяции в пределах ячейки. В общем случае для большинства схем типа MUSCL параметры на гранях расчетных ячеек могут быть определены с использованием следующих соотношений [10]: Й+1/2 = Qn + 4[(1 -k)(VQ) + (1 + k)(AQ)]; (5) QR-1/2 = Qn -4[1 + k)(VQ) + (1 -k)(AQ)]; (6) aq = Q+ - Qn ; VQ = Qn - Q" 1, где Ql+1/2 - параметр на левой границе ячейки; QR-\/2 - параметр на правой границе ячейки; Qn - параметры в центрах ячеек; k - коэффициент схемы. Возможные варианты реализации схем типа MUSCL представлены в табл. 1. Значения коэффициента схемы k = -1,0,1 соответствуют схемам с линейной интерполяцией величин, а k = 1/3 - схеме третьего порядка точности с параболической интерполяцией. QL+1/2 = Qn + 4 (7) (8) AQ Qi-1/2 = Qi ^ I r + r21 V1 + r: , ; Van Leer тели: MINMOD - y(r) = MINj^,3-kr\; Van Albada - y(r) = - y(r ) = ^ ; Superbee - r )= MAX [0, MIN (2r,1), MIN (r,2)] и другие [9]. Таблица 1 Типы схем MUSCL Название схемы Коэффициент схемы Тип реконструкции Порядок точности MUSCL-Upwind -1 Линейная 2-й MUSCL-Fromm 0 Линейная 2-й MUSCL-CDS 1 Линейная 2-й MUSCL-CUS 1/3 Параболическая 3-й Для выполнения условий сохранения полной вариации (TVD) схема (5), (6) должна быть дополнена функцией ограничителем наклона [15, 16]: (-kн ^Vq + (1 + k) JЩ\AQ VQ (+k h^ |vq+(1 - k um\aq У Q^ ' \aq где r) - функция-ограничитель. Для семейства схем MUSCL могут быть использованы следующие ограничи Таким образом реализация шага по времени в схеме MUSCL проводится в три этапа: Этап 1 - «Реконструкция». На этом этапе по известным значениям Qi, с использованием соотношений (7), (8) определяются параметры QL , QR для каждой грани расчетной ячейки. Этап 2 - «Решение задачи Римана». На данном этапе по известным значениям QL , QR с использованием точного или приближенного решения задачи Римана о распаде произвольного разрыва определяются потоки F+i/2 и F-i/2 через грани расчетной ячейки: F41/2 - RimmanSolver(QL+1/2,QR+1/2). (9) Этап 3 - «Вычисление параметров на новом временном слое». На этом этапе по известным значениям Qt", F+1/2 и F-1/2 с использованием (3) определяются параметры газа Qin+1 на новом временном слое. Схема «MUSCL-Hancock» Одним из простых способов повышения порядка точности по времени является использование схемы MUSCL-Hancock [17,18]. Идея схемы заключается в том, что потоки F+1/2 и F-1/2 через грани расчетной ячейки определяются на половинном шаге по времени (t + 2 ), а затем эти значения используются для определения решения на следующем временном слое (t + Д). Реализация шага по времени в схеме MUSCL-Hancock проводиться в четыре этапа: Этап 1 - «Реконструкция». На этом этапе по известным значениям Q4, с использованием соотношений (7), (8) определяются параметры QL , QR для каждой грани расчетной ячейки. Этап 2 - «Эволюция». На данном этапе по известным значениям QL , QR , определяются параметры QL , QR на гранях расчетных ячеек для момента времени (t + Д/ 2): QL+1/2 -QL+1/2 + 2Д[F(QL+1/2)-F(Qf+1/2)] ; QR-1/2 - QR+1/2 + 2Д[F((+ш)-F(Qf+!/2)], где F (Q^j^ ) и F (Q-+1/2 ) определяются с использованием (2). Этап 3 - «Решение задачи Римана». На данном этапе по известным значениям QL и QR , с использованием точного или приближенного решения задачи Римана о распаде произвольного разрыва, определяются потоки Fi+1/2 и Fi-1/2 через грани расчетной ячейки: F41/2 - RimmanSolver(QL+1/2,QRU2) (10) Этап 4 - «Вычисление параметров на новом временном слое». На этом этапе по известным значениям Qi, F+1/2 и F-1/2 с использованием (3) определяются параметры газа Qi+1 на новом временном слое. Численные исследования Проведена серия расчетов квазиодномерного течения идеального газа в сопле Лаваля [12] с использованием комбинаций схем «Годунова», «MUSCL», «MUSCL-Hancock», с функциями-ограничителями: «MINMOD», «Van Albada», «Van Leer», «Superbee» и методами нахождения потоков: «Роя», «HLLC», «Точное решение задачи Римана». Сравнительный анализ проводился по результатам определения коэффициента расхода. В качестве тестовой геометрии использовано сопло JPL (рис. 1), где R - радиус минимального сечения. Рис. 1. Геометрия сопла JPL Fig. 1. JPL nozzle geometry На входе в сопло ставились граничные условия - постоянство энтальпии и энтропии, на срезе сопла - сверхзвуковое истечение. Рабочее тело: воздух. Расчет проводился методом установления. Количество узлов расчетной сетки: 500. Коэффициент расхода определяется как G* где G* - расчетный расход газа через критическое сечение сопла, Gg - теоретический расход, определенный для адиабатического течения. Результаты численных исследований по определению коэффициентов расхода для различных комбинаций базовых схем, функций ограничителей и методов определения потоков через грани расчетных ячеек представлены в табл. 2. Сравнение по базовым схемам. Для схемы «Годунова» первого порядка точности по пространству, расход газа через критическое сечение сопла отличается от теоретического на 1.33 - 1.34 %. Для схем второго порядка точности по пространству («MUSCL-Upwind», «MUSCL-Fomm», «MUSCL-CDS») - на 0.0170.08 %. Схема «MUSCL-CUS» третьего порядка точности по пространству дает значения расхода газа, отличающиеся от теоретического на 0.017-0.07 %, а схема «MUSCL-Hancock» второго порядка точности по времени и третьего по пространству - на 0.07-0.17 %. Таблица 2 Коэффициенты расходов для различных комбинаций схем типа «MUSCL» № Базовая схема Порядок схемы Ограничитель Метод решения задачи Римана Коэффициент расхода 1 Годунов 1 - Точное решение 0.986510 2 Годунов 1 - Рое 0.986682 3 Годунов 1 - HLLC 0.986602 4 MUSCL-Upwind 2 minmod Точное решение 0.999138 5 MUSCL-Upwind 2 Superbee Точное решение 1.000170 6 MUSCL-Upwind 2 VanAlbda Точное решение 0.999558 7 MUSCL-Upwind 2 VanLear Точное решение 0.999637 8 MUSCL-Upwind 2 VanLear Рое 0.999636 9 MUSCL-Upwind 2 VanLear HLLC 0.999662 10 MUSCL-Fomm 2 minmod Точное решение 0.999138 11 MUSCL-Fomm 2 Superbee Точное решение 1.000171 12 MUSCL-Fomm 2 VanAlbda Точное решение 0.999558 13 MUSCL-Fomm 2 VanLear Точное решение 0.999637 14 MUSCL-Fomm 2 VanLear Рое 0.999636 15 MUSCL-Fomm 2 VanLear HLLC 0.999662 16 MUSCL-CDS 2 minmod Точное решение 0.999138 17 MUSCL-CDS 2 Superbee Точное решение 1.000173 18 MUSCL-CDS 2 VanAlbda Точное решение 0.999558 19 MUSCL-CDS 2 VanLear Точное решение 0.999637 20 MUSCL-CDS 2 VanLear Рое 0.999636 21 MUSCL-CDS 2 VanLear HLLC 0.999662 22 MUSCL-CUS 3 minmod Точное решение 0.999138 23 MUSCL-CUS 3 Superbee Точное решение 1.000170 24 MUSCL-CUS 3 VanAlbda Точное решение 0.999558 25 MUSCL-CUS 3 VanLear Точное решение 0.999637 26 MUSCL-CUS 3 VanLear Рое 0.999636 27 MUSCL-CUS 3 VanLear HLLC 0.999662 28 MUSCL-Hancock 3 minmod Точное решение 0.998236 29 MUSCL-Hancock 3 Superbee Точное решение 0.999213 30 MUSCL-Hancock 3 VanAlbda Точное решение 0.998637 31 MUSCL-Hancock 3 VanLear Точное решение 0.998710 32 MUSCL-Hancock 3 VanLear Рое 0.998710 33 MUSCL-Hancock 3 VanLear HLLC 0.998733 Сравнение по функциям-ограничителям. Использование ограничителя «minmod» дает значения расхода газа, отличающиеся от теоретического на 0.086 % (для схем «MUSCL-Upwind», «MUSCL-Fomm», «MUSCL-CDS», «MUSCL-CUS») и 0.176 % (для схемы «MUSCL-Hancock»). Ограничитель «Superbee» дает значения расхода газа, отличающиеся от теоретического на 0,017 % (для схем «MUSCL-Upwind», «MUSCL-Fomm», «MUSCL-CDS») и 0.078 % (для схемы «MUSCL-Hancock»), ограничитель «VanAlbda» - 0.044 и 0.136 %, а ограничитель «VanLear» - 0.036 и 0.129 % соответственно. Сравнение по методам решения задачи Римана. Коэффициент расхода, вычисленный по схеме «Годунова» с использованием приближенных методов решения задачи Римана о распаде произвольного разрыва (метод Роя, метод HLLC) отличается от коэффициента расхода, вычисленного с использованием точного метода решения задачи Римана на 0.009 - 0.071%. Коэффициент расхода, вычисленный по схемам «MUSCL» второго и третьего порядка точности по пространству, с использованием приближенных методов решения задачи Римана о распаде произвольного разрыва отличается от коэффициента расхода, вычисленного с использованием точного метода решения задачи Римана на 0 - 0.002 % Заключение По результатам расчетов локальных и интегральных характеристик в сопле Лаваля, можно сделать следующие выводы. Схема «Годунова» первого порядка точности по пространству дает достаточно высокую погрешность при определении коэффициента расхода (1.4 %). Схемы «MUSCL» второго и третьего порядка точности по пространству дают достаточно низкую погрешность при определении коэффициента расхода (0.017-0.17 %) вне зависимости от функции ограничителя. Следует отметить, что результаты, полученные с использованием схем второго порядка точности по пространству, незначительно отличаются от результатов, полученных с использованием схем третьего порядка точности. Выбор функции-ограничителя и метода решения задачи Римана (точное решение задачи Римана, метод Роя, метод HLLC) влияют на точность определения интегральных характеристик в пределах сотых долей процента.

Скачать электронную версию публикации