CC BY
87
УЧЕНЫЕ ЗАПИСКИ ЦА Г И Том XVII 1986
№ 2
УДК 517 1-318 : 519 : 533.6
ОБОБЩЕНИЕ,МОДИФИЦИРОВАННОЙ СХЕМЫ С. К. ГОДУНОВА НА ПРОИЗВОЛЬНЫЕ НЕРЕГУЛЯРНЫЕ СЕТКИ
Н. И. Тилляева
Описывается модификация монотонной разностной схемы С. К. Годунова [1], сохраняющая аппроксимацию при численном решении задач газовой динамики на произвольных нерегулярных сетках, адаптированных к рассматриваемой задаче. Важным элементом предлагаемой схемы, названной схемой Годунова — Колгана (СГК), является использование введенного В. П. Колганом [2—4] принципа минимальных значений производной (точнее, приращений — ПМП). В задачах с двумя независимыми переменными, а при большем числе измерений для регулярных сеток из однотипных (в двумерных нестационарных задачах — четырехугольных) ячеек произвольных размеров при отсутствии конечных (не стремящихся к нулю с измельчением сетки) изломов сеточных линий СГК практически совпадает со схемой [2—4]. Приводятся примеры, подтверждающие эффективность СГК.
1. Характерной особенностью большинства задач, встречающихся в приложениях, является нерегулярность границ рассматриваемых областей течения. К тому же на линиях излома обтекаемых идеальным газом поверхностей возникают пучки волн разрежения, аккуратный расчет которых невозможен без использования разностных сеток, адаптированных к каждому излому. Последнее вместе со сложностью рассчитываемых областей приводит к необходимости вести счет на существенно нерегулярных сетках. С другой стороны, широко распространенные разностные схемы, как правило, обеспечивают аппроксимацию уравнений течения только на достаточно хороших сетках, не имеющих изломов сеточных линий и линий скачкообразного изменения размеров разностных ячеек. Отмеченные обстоятельства делают актуальной задачу построения разностных схем, обеспечивающих хотя бы с первым порядком аппроксимацию исходных дифференциальных или интегральных уравнений течения на произвольных сетках. Именно такой схемой является описываемая ниже СГК, которая, как уже отмечалось, представляет собой модификацию схемы Годунова (СГ). Одним из важных элементов предлагаемой модификации является ПМП, сохраняющий монотонность схемы. Монотонность разностной схемы важна прежде всего потому, что при «сквозном» счете разрывов по монотонным схе-
мам отсутствуют нефизические осцилляции, затрудняющие восприятие и снижающие точность результатов (погрешности из-за размазывания разрывов, снижающие при сквозном счете порядок аппроксимации любой схемы в общем случае до первого [5], пропорциональны ширине всей структуры разрыва, а отнюдь не только участка максимального градиента). Мало того, указанные осцилляции подчас делают невозможным и само получение результатов, например, из-за забросов давления р или плотности р в область отрицательных значений. Благодаря монотонности СГ оказалась, пожалуй, самой работоспособной (безаварийной) схемой, что подтверждается обилием и многообразием задач, решаемых с ее помощью. Вместе с тем СГ не лишена недостатков. Это прежде всего низкий (первый) порядок аппроксимации дифференциальных уравнений даже на хороших сетках и сильное размазывание слабых скачков и тангенциальных (контактных) разрывов.
В СГК для задач, решаемых установлением, порядок аппроксимации дифференциальных уравнений установившегося решения на хороших сетках повышается до второго с сохранением монотонности в процессе установления. Последнее достигается благодаря тому, что в СГК распаду разрыва при определении так называемых «больших величин» предшествует интер- или экстраполяция параметров в середину каждой границы или грани ячейки (с обеих ее сторон). По этой причине в СГК в любой ячейке кроме средних параметров используются их пространственные производные. Интер- или экстраполяция применяются в зависимости от величин приращения, что и составляет ПМП.
Хотя в [2—4] были продемонстрированы монотонность и работоспособность предложенной там модификации СГ и показано, что она значительно меньше, чем СГ, размазывает слабые скачки и контактные разрывы (см. также работы [6—8]), описанная в [2] модификация СГ не получила должного распространения и развития. Из причин, помешавших этому, главными представляются следующие. Во-первых, возможность повышения точности, связанная со вторым порядком аппроксимации уравнений установившегося течения, не реализовывалась в полной мере из-за, размазывания всех скачков и недостаточно точной аппроксимации граничных условий. Во-вторых, на тех же сетках модификация СГ из [2] из-за меньшего максимально допустимого шага интегрирования по времени требует увеличения времени счета. В-третьих, практически все расчеты в рамках указанной модификации велись на почти равномерных сетках, не адаптированных к данной задаче. В то же время на основе соображений [2] и прежде всего ПМП легко предложить почти очевидное обобщение [2], которое, сохраняя монотонность, обеспечивает первый порядок аппроксимации полной системы уравнений нестационарного течения на произвольных сетках. Далее под СГК понимается именно такое обобщение, существо которого разъясняется ниже на примере разностной схемы для нестационарных течений.
Сначала, не заботясь о монотонности и консервативности предлагаемой схемы, покажем, как на любой сетке можно обеспечить разностную аппроксимацию уравнений. Для этого рассмотрим произвольную ячейку, не ограничивая числа ее сторон в двумерном случае или граней— в пространственном. Наряду со значениями параметров в некоторой ее точке О на уже известном п-и временном слое способом, описанным ниже, найдем с погрешностями О (/г) все их пространственные производные. По ним с помощью отрезков рядов Тейлора найдем на том же слое с погрешностью О (/г2) отличия от параметров в точке О их значений в «центрах тяжести» (ЦТ) граней (сторон) ячейки. Най-
денные величины используем затем, взяв за О ЦТ ячейки, при записи для нее на временном интервале т интегральных законов сохранения массы, импульса и энергии. Несложный анализ показывает, что при этом погрешности их разностной аппроксимации есть 0[т/гу (/г+т)] с у=2 и 3 соответственно в двух- и трехмерном случае, а погрешности в имеющих порядок х приращениях параметров при переходе с я-го на (я+1)-й слой оказываются 0[т(/г + т)]. При установлении интегральные законы сохранения потоков, каждый из которых на отдельной грани есть О (Л"-1), записываются с погрешностью О (Л''+1). Данные оценки показывают, что и в нестационарном случае, и после установления для любой сетки имеет место аппроксимация уравнений с первым порядком. Если сетка равномерна, то (Л + т) из-за частичной компенсации ошибок заменится на (/г2+т), что при установлении повышает порядок аппроксимации до второго.
Различие параметров и производных в соседних ячейках при описанном способе разностной записи уравнений вело бы в подобластях непрерывности течения к отличию в 0[т/г (Л + т)] потоков с разных сторон каждой грани, нарушая, без влияния на аппроксимацию, консервативность схемы. Если в областях непрерывности течения эти нарушения невелики, то на размазанных разрывах отличие указанных потоков было бы того же порядка, как сами потоки, что недопустимо. Поэтому для обеспечения консервативности схемы будем определять потоки по большим величинам, которые, как и в СГ, будем находить из задачи о распаде произвольного разрыва. При этом в отличие от СГ в упомянутой задаче за исходные возьмем не параметры в ЦТ примыкающих к грани ячеек, а параметры в ЦТ грани с разных ее сторон.
Пространственные производные от параметров в точке О и их приращения бф между значениями в точке О и в ЦТ всех граней данной ячейки, вопрос определения которых пока оставался открытым, вычисляются так. На я-м слое берется К ячеек, примыкающих по граням к рассматриваемой. Их ЦТ пронумеруем числами &=1,...,/С. С погрешностью О (А2) имеем: Дсрк = сру. 0 Дл*. Здесь — координаты, по всем \ от 1 до V проводится суммирование, <?} = ду1дх]г Дфк = срк — <р0 с нижними индексами, приписываемыми величинам в соответствующих точках.
В выписанных равенствах А<рк и Ах) известны, что при /С>2у позволяет найти не менее двух «односторонних» значений каждой производной фзо, различающихся в областях непрерывности параметров на
О (/г). По каждому набору ф, 0 вычислим приращение бф от точки О до ЦТ соответствующей грани. Если значения бф, отвечающие разным <р;0, имеют один знак, то для такой грани из двух бф выберем минимальное по модулю. В этом и состоит ПМП. В случае бф разных знаков заменим бф нулем. Подобная замена, как показал С. К- Щипин, не ухудшая аппроксимации и монотонности схемы, увеличивает ее шаг хт до 2/3 от хт для СГ. При v = 2 для четырехугольных ячеек естественно брать К=4, а затем находить ф^о и бф, решая попарно уравнения, которые связывают параметры в точке О и в верхних правых и нижних левых точках &.
Многочисленные расчеты продемонстрировали монотонность и работоспособность описанной выше версии СГК в широком круге задач газовой динамики. В дополнение к этому и к незначительному размазыванию ударных волн практически любой интенсивности ширина контактных (тангенциальных) разрывов в СГК растет только на начальном участке [8, 9]. К достоинствам СГК относится и простота перехода к
ней от СГ. При этом в программах, использующих установление, блок перехода от СГ к СГК включается после того, как от нескольких десятков до сотни шагов просчитано по СГ. Кроме того, в таких программах нет необходимости вычислять cpjo и 6ф на каждом слое, что ведет к значительному ускорению вычислений. Подобное замораживание на несколько слоев производных (или приращений) обеспечивает монотонный характер установления решения и дает вбзможность вести устойчивый счет с тем же временным шагом, что и в СГ. Это позволяет на одинаковых сетках получать по СГК существенно более точные результаты практически с теми же затратами времени, что в СГ. Заметим, наконец, что описанный выше элементарный прием вычисления производных, распространенный на любой временной слой, может оказаться полезным при аналогичном обобщении на произвольные сетКи и других разностных схем.
Реакцию СГ и СГК на неравномерность сетки демонстрирует рис. 1, где представлены результаты расчета обтекания затупленной по торцу пластины потоком с М.х = 3. Здесь и далее газ совершенный с показателем адиабаты х— 1,4. Использовались две сетки: сравнительно хорошая и преднамеренно деформированная (рис. 1 , а). Деформация состояла в образовании елочки, пересекающей расчетную область. Изобары р/р00 = const, ‘ полученные на деформированной сетке по СГ (рис. 1,6) и по СГК (рис. 1,е), даны сплошными линиями, а на исходной— штрихами. Видно, что СГК лишь слегка реагирует на весьма грубую деформацию сетки. Это и понятно, поскольку вблизи изломов в СГК снижается лишь порядок аппроксимации. Напротив, изломы изобар на рис. 1,6 указывают на отсутствие в тех же случаях аппроксимации в СГ. Заметим, кстати, что наблюдаемые погрешности СГ на нерегулярной сетке имеют место, несмотря на то, что для уравнений, записанных в декартовых координатах СГ, как и СГК, является согласованной схемой в смысле [10]. Напомним, что согласованные схемы на любых сетках не искажают равномерного поступательного потока.
Для исходной недеформированной сетки СГК за счет более высокого порядка аппроксимации во всей рассчитываемой области также обеспечивает лучшую точность, чем СГ. Так, в данном примере ошибки по'р/рх и по полной энтальпии / в точке торможения, достигающие в СГ 8,6 и 5%, в СГК на той же сетке снизились до 0,1%. С удалением
по торцу от точки торможения погрешности обеих схем растут, достигая наибольших значений перед изломом (в СГ — 12 и 8%, в СГК— 1,3 и 1,5%, т. е. примерно на порядок меньше). Это и понятно, так как в данной задаче производные pv = dpjdy,... при приближении к излому стремятся к бесконечности [И]. В промежуточных точках торца СГК, как и в точке торможения, уменьшает погрешности в определении /?/р* и / в десятки раз. Достигнутое уменьшение ошибок связано с тем, что в общем случае в окрестностях точек торможения компоненты вектора скорости и, v, w близки к нулю. Поэтому, если т — порядок аппроксимации схемы, то, как нетрудно показать, там, где и, ... = 0(h), погрешности определения плотности р и температуры Т, а следовательно, и их функций — piрх и I есть О (hm~l). Отсюда ясно, что схемы первого порядка определяют перечисленные величины в точке торможения с погрешностями порядка самих этих величин.
Еще одно подтверждение малой чувствительности СГК к качеству сетки дает рис. 2. На нем для обтекания с числом М<х, = 2 заостренного плоского тела (схема на рисунке) с полууглом при вершине 45° дано распределение по телу ошибок величины I. На рассчитанном режиме скачок отошедший, и рх, ру,... в точке торможения обращаются в бесконечность, что, естественно, затрудняет получение аккуратных результатов. Сплошной ломаной на рис. 2 показаны ошибки расчета по СГК на простейшей сетке, образованной лучами w = const с равномерным разбиением отрезков между стенкой и скачком (ю—угол полярных координат с центром внутри тела, число ячеек jV= 10X10, головной скачок везде, как и ранее, выделялся). Ошибки расчетов, выполненных в [12] на подобной сетке с N=19x19 и на специальным образом улучшенной (адаптированной) сетке с N= 11x11, даны на рис. 2 пунктиром и штрихами. На хороших сетках схема [12] дает весьма точные результаты. Например, в задаче обтекания кругового цилиндра она слегка превосходит использованную в настоящей работе реализацию СГК. Тем не менее в более сложном случае (рис. 2) ни почти четырехкратное увеличение N, ни специальная адаптация сетки не позволили с ее помощью достичь той точности, которую и без подобных мер дает СГК.
Следующие два примера иллюстрируют возможности СГК при решении более сложных задач. Первый из них (рис. 3) относится к обтеканию конфигурации, содержащей цилиндр с осью, перпендикулярной ПЛОСКОСТИ течения, сверхзвуковым ПОТОКОМ С Мех. = 6. Помимо цилиндра исследованная конфигурация содержит клин, поворачивающий поток на 15°. Клин помещен ниже и впереди цилиндра. Часть косого скачка, идущего от вершины клина, дана на рис. 3, а и б отрезком прямой, который приходит в точку 1. Из-за присутствия клина течение перед криволинейным отошедшим от цилиндра скачком оказывается кусочно-
Рис. 3
равномерным (за косым скачком М~4), что, как известно [13—16], вызывает качественную перестройку обтекания цилиндра (по сравнению со случаем цилиндра в равномерном потоке). В частности, теперь на цилиндре при некоторых положениях клина появляются зоны аномально высоких давлений и тепловых потоков.
В рассматриваемом примере по СГК рассчитывался ударный слой, ограниченный отошедшим скачком с двумя тройными точками (точки
1 и 2 на рис. 3,а). Внутри слоя располагаются тангенциальные разрывы Т1 и Т2, начинающиеся соответственно в точках 1 и 2, ударная волна выходящая из точки 2, и пучок волн разрежения с центром в I — точке падения № на В / скорость потока за больше или равна скорости звука. Прилегающая к Т± часть сверхзвукового потока после прохождения через № и разворота в пучке волн разрежения течет вдоль верхней поверхности цилиндра, разгоняясь до больших сверхзвуковых скоростей. Хотя в СГК перечисленные выше внутренние особенности не выделялись, они неплохо видны на рис. 3, где изображены линии тока (рис. 3, а, цифры — значения функции тока) и изомахи (рис. 3,6). Те и другие нарисованы через равные интервалы. Сгущения .изомах, наиболее сильные у тела, указывают на размазанные тангенциальные разрывы Т1 и Т2 и скачок V? и на область сильно неравномерного (из-за центробежного эффекта у стенки) сверхзвукового потока.
В точку торможения 5 приходит разделительная линия тока (рис. 3, а), прошедшая через три косых скачка. Число Маха сразу за третьим скачком V? не намного больше единицы. Именно, М=1,7 в точке 1 при 1<М<1,7 в точке /. Эти значения М найдены из анализа с помощью ударных поляр задач о пересечении и о расщеплении скачков в точках 1 и 2. Значения давления торможения р8, отвечающие М=1,7 и 1 за скачком № в предположении отсутствия дополнительных потерь полного давления между V? и точкой 5, равны 8,9 и 6,8. Если же принять, что в точке пересечения с разделительной линией скачок V? прямой, то рв = 5,7. Здесь и на рис. 3, в давление отнесено к роо УоЛ где V — модуль скорости. Для Мос = 6 обезразмеренное таким способом давление за прямым скачком равно 0,93.
Учитывая пикообразный характер зависимости р от полярного угла и, следует признать, что максимальная величина р = 5,9 на сплошной кривой рис. 3, в, дающей найденное по СГК распределение р по телу, находится в удовлетворительном согласии с приведенными выше цифрами, а также с экспериментальным значением р8^6,7 из [15]. То же самое можно сказать и о всей экспериментальной кривой из [15] (пунктир на рис. 3,в), тем более что показанные на рис. 3 результаты получены для идеального газа без поправки на толщину вытеснения пограничного слоя 3*. Отметим, что на рассмотренном режиме 5* много меньше толщины ударного слоя. Поэтому здесь нет необходимости проводить решение в рамках полных уравнений Навье—Стокса, как это делалось в [14—16], где, кстати, найденное из расчета при более чем втрое большем (по сравнению с данной работой) числе расчетных точек значение = 5,7. Дополнительное подтверждение эффективности СГК дает рассмотрение штриховой кривой на рис. 3, в, полученной по СГ на той же сетке, что и в СГК.
Погрешности СГ по другим параметрам в этом примере еще больше, особенно над цилиндром, где, помимо прочего, сказывается более сильное, чем в СГК, размазывание тангенциального разрыва. Так, при в> = 15° полное давление на стенке по СГ составляет лишь 0,47 р8, а по СГК — 0,81^ (числа М при этом оказываются равными
1,9 и 3,1).
Результаты, представленные на рис. 3, получены на равномерной по о>, а при каждом о>—и по г сетке с Ы=72X20 (72 — число ячеек по (о) при угле раствора расчетной области по о) — 160°. Расчет велся установлением по времени с использованием алгоритма выделения скачка, разработанного в [17, 18]. В данном примере способность СГК к аппроксимации на произвольных сетках существенна на лучах <о = = сойв!, приходящих в точки 1 и 2: На них сеточные линии противопо-
ложного семейства претерпевают излом. В разностной схеме, примененной в [14—16], такие изломы ведут к потере аппроксимации.
Последний пример относится к обтеканию сверхзвуковой нерасчетной струей осесимметричного тела с закрытым протоком. Здесь область, которая рассчитывалась по СГК, ограничена слева криволинейным скачком, а сверху — границей струи. Набегающая нерасчетная струя истекала из сопла с Мо=1,01 при ро/ре = 4, где ро и ре — давление соответственно на срезе сопла и в окружающем пространстве. Струя рассчитывалась по маршевой схеме [19] и запоминалась для получения интерполяцией параметров перед искривленным скачком. Результаты расчета предствлены на рис. 4, где в области за скачком показаны линии тока (рис. А, а) и изобары (рис. 4,6). Изобары построены через постоянный интервал, р отнесено к ре. Линии тока при 0,01 <ф<0,1, где ф—функция тока, нарисованы через 0,01, а при 0,1<4'<1 — через 0,1. Наряду с этим на рис. 4 даны граница всей струи и висячий скачок. Струя считалась без выделения висячего скачка, который в данном примере сравнительно слаб. В силу этого линия скачка строилась по точкам максимума \др/ду\. Несмотря на отмеченное обстоятельство и малое число расчетных точек (270) в расчетной области, рис. 4 дает вполне правдоподобное представление о рассмотренном течении.
В заключение автор благодарит А. Н. Крайко за полезные обсуждения, а также Э. М. Дорофееву, Л. И. Дербеневу, В. Е. Макарова и В. А. Вострецову за помощь в работе.
ЛИТЕРАТУРА
1. Годунов С. К-, 3 а б р о д и н А. В., Иванов М. Я., К Р а й-к о А. Н., Прокопов Г. П. Численное решение многомерных задач газовой динамики. — М.: Наука, 1976.
2. К о л г а н В. П. Применение принципа минимальных значений производной к построению конечно-разностных схем для расчета разрывных решений газовой динамики. — Ученые записки ЦАГИ, 1972, т. III, № 6.
3. К о л г а н В. П. Конечно-разностная схема для расчета двумерных разрывных решений нестационарной газовой динамики. — Ученые записки ЦАГИ, 1975, т. VI, № 1.
4. К о л г а н В. П. Численный метод решения пространственных задач газодинамики и расчет обтекания тела при наличии угла атаки. —■ Ученые записки ЦАГИ, т. VI, № 2, 1975.
5. И в а н о в М. Я-, Крайко А. Н. Об аппроксимации разрывных решений при использовании разностных схем сквозного счета. — Ж. вычислит, матем. и мат. физ., т. 18, № 3, 1978.
6. Косых А. П., М и н а й л о с А. Н. Исследование методов сквозного счета для задач сверхзвуковой аэродинамики. — Ученые записки ЦАГИ, 1976, т. VII, № 1.
7. Минайлос А. Н. О значении монотонности конечно-разностных схем в методах сквозного счета.—-Ж. вычислит, мат. и мат. физ., 1977, т. 17, № 4.
8. Лобановский Ю. И. О монотонизации конечно-разностных решений в методах сквозного счета. — Ж. вычислит, мат. и мат. физ., 1979, т. 19, № 4.
9. Копчен о в В. И., Крайко А. Н. Монотонная разностная схема второго порядка для гиперболических систем с двумя независимыми переменными.— Ж. вычислит, мат. и мат. физ., 1983, т. 23, № 4.
10. Flores J., Holst Т. L., Kwak D., Batiste D. M. A new-consistent spatial differencing scheme for transonic full-potential equation.—AIAA J., vol. 8, 1984 (Русский пер.: Аэрокосмическая техника, № 2, 1985).
11. V a g 1 i o-L a u r i n R. Transonic rotational flow over a convex corner. — J. Fluid Mech. 1960, vol. 9, N 1.
12. Rai М. М., Anderson D. A. Application of adaptive grids to fluid-flow with asymptotic solutions. AIAA J., vol. 20, N 4, 1982.
13. Edney В. E. Anomalous heat transfer and pressure distribution on blunt bodies at hypersonic speeds in the presense of an impinging shock. — The Aeronautical Research Institute of Sweden, Stockholm, FFA Rep. 115, 1968.
14. T a n n e h i 11 J. C., Holst T. L., R a к i с h J. V. Numerical computation of two-dimensional viscous blunt body flows with an impinging shock. — AIAA J., 1976, vol. 14, N 2.
15. Tannehill J. C., Holst T. L., Rakich J. V., Keyes J. W. Comparison of a two-dimensional shock impingement computation with experiment.— AIAA J., 1976, vol. 14, N 4.
16. T a n n e h i 11 J. C., VigneronY. C., Rakich J. V. Numerical solution of two-dimensional turbulent blunt body flows with an impinging shock. — AIAA J., 1979, vol. 17, N 12.
17. Тилляева H. И. Численный метод расчета обтекания плоского воздухозаборника сверхзвуковым потоком на режимах с выбитой ударной волной. — Ученые записки ЦАГИ, т. X, № 2, 1979.
18. Крайко А. Н., Макаров В. Е., Тилляева Н. И. К численному построению фронтов ударных волн. — Ж. вычисл. мат. и мат. физ., 1980, т. 20, № 3.
19. Иванов М. Я., Крайко А. Н., М и х а й л о в Н. В. Метод сквозного счета для двумерных и пространственных сверхзвуковых течений, ч. I. — Ж. вычислит, мат. и мат. физ., 1972, т. 12, № 2.
Рукопись поступила 15/IV 1985
