CC BY
368
УЧЕНЫЕ ЗАПИСКИ Ц А Г И
Том VIII 1 977 М3
УДК 533.6.011.3/5
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ СЛОЖНЫХ ЗАДАЧ АЭРОГАЗОДИНАМИКИ МЕТОДОМ „КРУПНЫХ ЧАСТИЦ“
I. МЕТОД. ИССЛЕДОВАНИЕ СХЕМ
О. М. Белоцерковский, Ю. М. Давыдов
В статье дается описание метода „крупных частиц* для решения полной системы вихревых уравнений газовой динамики. Приводится постановка граничных условий. Полученные разностные схемы анализируются с помощью аппарата дифференциальных приближений.
1.1. За последнее время практика выдвинула перед исследователями ряд важных задач, решение которых требует тщательного изучения вихревых околозвуковых течений газа, течений со „вду-вом“ струи в основной поток, срывных турбулентных зон и т. п. Интересно, естественно, иметь здесь динамику развития явления, получить полные картины течений и характерные их свойства в широком диапазоне изменения параметров среды, тела и др.
В трансзвуковых областях имеют место, например, смешанные течения газа сложного характера, могут возникнуть так называемые „вторичные“ висячие скачки уплотнения, местные сверхзвуковые зоны и т. п. Механизм образования этих явлений в достаточной мере не изучен; остаются открытыми вопросы о виде минимальной области влияния и о ее перестройке при изменении параметров обтекания, о причинах и условиях возникновения вторичных скачков уплотнения, о свойствах локальных сверхзвуковых зон. Большой практический интерес вызывают также исследования характеристик течений со „вдувом“ струи в следе за телом. Необходимо изучить возможности моделирования таких явлений с помощью ЭЦВМ.
Трудности теоретического исследования таких задач определяются весьма сложной математической постановкой (краевые задачи для нелинейных уравнений, эллиптическо-гиперболического типа, сжимаемого газа). Только численные методы с использованием быстродействующих электронно-счетных машин и тщательно проведенные эксперименты позволяют получить полное решение указанных выше задач и определить необходимые характеристики
течения. Таким образом, разработка численных схем, построение вычислительных алгоритмов, проведение расчетов различных сложных задач газовой динамики, а также изучение аналитических свойств решений, построение асимптотик и т, п. представляет в настоящее время несомненный интерес.
Ни в коей мере не умаляя значения и решающей роли физического эксперимента, хотелось бы лишь отметить, что для хорошо смоделированной задачи информация, получаемая из расчетов, значительно полнее и экономически существенно дешевле соответствующих экспериментальных исследований. В ряде же случаев (трансзвуковые задачи, турбулентные потоки, течения с физико-химическими превращениями, излучением и др.) данные эксперимента или носят весьма ограниченный характер, или вовсе отсутствуют.
В частности, для вихревых плоских и пространственных трансзвуковых задач аналитические и численные методы исследования практически не были разработаны. Классические аналитические методы, развитые для плоских потенциальных течений и основанные преимущественно на использовании плоскости годографа и различного типа упрощений и линеаризаций, здесь, как правило, не проходят. Численное решение указанных задач также вызывает много трудностей и требует построения специальных схем для интегрирования уравнений эллиптико-гиперболического типа. Это объясняется, видимо, тем, что большинство „традиционных“ численных схем строилось ранее специально для газодинамических задач чисто гиперболического или чисто эллиптического типов при использовании определенных свойств этих уравнений. Случаи же смешанных течений газа (особенно для закритических областей), а также малых сверхзвуковых (или больших дозвуковых) скоростей, изучение турбулентных течений в зонах срыва и т. п. особенно трудны для численного исследования. Это связано с необходимостью построения сложных математических моделей, решением уравнений смешанного типа в неограниченных областях и постановкой специальных краевых условий задачи. Во многих случаях в процессе расчета необходимо осуществлять переход через скорость звука, проводить построение внутренних скачков уплотнения, учитывать эффекты нестационарности и вязкости газа и т. д.
1.2. Для расчета течений в сложных задачах аэрогазодинамики естественно использовать нестационарные схемы „сквозного“ счета, где вычисления проводятся без предварительного выделения особенностей, поверхностей разрыва и т. п. В ряде случаев кажется рациональным введение в алгоритмы элементов метода Харлоу „частиц в ячейках“ [1], применение которого граничит с проведением численного эксперимента.
В последние годы в Вычислительном центре АН СССР был осуществлен ряд численных экспериментов по исследованию сложных газодинамических течений. Они были проведены с помощью нестационарного метода „крупных частиц“, толчком к разработке которого послужили работы Ф. Харлоу [1], М. Рича [2] и С. Херта [3]. Основная цель этих разработок — исследование широкого класса задач современной аэрогазодинамики (от чисто дозвуковых течений до сверхзвуковых режимов, включая трансзвуковые области, зоны отрыва и т. п.). Полученные схемы используют сравнительно небольшое число расчетных точек (обычно не более двух-трех тысяч), что позволяет применять их также для проведения серийных расчетов в практике прикладных НИИ и КБ.
Описание метода „крупных частиц“ применительно к задачам газовой динамики содержится в работах [4 — 7]. Возможно, что отдельные локальные свойства численных решений будут определяться при этом недостаточно точно, но, видимо, лишь только таким путем можно получить общие характеристики сложных явлений и картину течений в целом.
Метод Харлоу „частиц в ячейках“ сочетает в себе в опредлен-ных чертах преимущества лагранжева и эйлерова подходов. Область решения здесь разбивается неподвижной фиксированной в пространстве (эйлеровой) расчетной сеткой; однако сплошная среда трактуется как дискретная модель — рассматривается совокупность „частиц“ фиксированной массы (лагранжева сетка частиц), которые движутся через эйлерову сетку ячеек. Частицы служат для определения параметров самой жидкости (массы, энергии, скорости), в то время как эйлерова сетка используется для определения параметров поля (давления, илотности, температуры).
Метод „частиц в ячейках“ позволяет исследовать сложные явления в динамике многокомпонентных сред, частицы хорошо „следят“ за свободными поверхностями и линиями раздела сред, взаимодействием разрывов и т. п. Однако дискретный метод частиц обладает и рядом недостатков. Главный из них, лежащий в самой природе метода, состоит в том, что из-за дискретного представления сплошной среды конечным числом частиц в ячейке параметры течения также определяются дискретным образом — как только частица пересечет границу эйлеровой ячейки, то масса, импульс, энергия частицы вычитаются из соответствующих величин прежней ячейки и прибавляются к новым значениям.
Такие скачки, весьма характерные для расчетов по методу Харлоу, приводят к большим нефизическим флуктуациям рассчитываемых величин (особенно плотности), в решениях появляются автоколебания и т. п. Кроме того, сами численные схемы этого метода обладают недостаточной вычислительной устойчивостью (особенно в областях застоя, при небольшом числе частиц в ячейке и др.), почему приходится вводить в схемы явные диссипативные члены с искусственной вязкостью, использовать неявные схемы первого шага, рассматривать „частицы“ различных форм [8] и т. д. Затруднительно также получение информации для сильно разреженных областей, откуда практически уходят все частицы и т. п. Значительно же увеличить число частиц не позволяют технические возможности современных вычислительных машин. По существу, метод Харлоу благодаря введению дополнительного параметра (числа „частиц“ в данной ячейке) увеличивает на единицу размерность задачи.
Для газодинамических течений однокомпонентной среды оказалось целесообразным отойти от дискретной модели „частиц“ и исходить из концепции непрерывности, рассматривая вместо „частиц“ поток массы через границы эйлеровых ячеек. Плотность газа здесь будет уже находиться не путем деления суммарной массы всех „частиц“ в ячейке на ее объем, а из закона сохранения массы, записанного для данной ячейки („крупной частицы“). При этом естественно сохранить сильные стороны метода Харлоу — эйлерово-лаг-ранжевый подход и сам процесс организации вычислений [2, 4 — 8, 11].
Таким образом, вместо совокупности ряда частиц в ячейках здесь рассматривается масса всей ячейки в целом — „крупная час-
з
тица“ (отсюда и название метода) и на основе конечно-разностных представлений законов сохранения изучаются нестационарные (и непрерывные) потоки этих „крупных частиц“ через эйлерову сетку. По существу, здесь используются законы сохранения, записанные в форме уравнений баланса для ячейки конечных размеров (как это обычно делается в процессе вывода газодинамических уравнений, но без дальнейшего предельного перехода от ячейки к точке).
Описание в различных вариантах численных схем метода „крупных частиц“ и полученных результатов содержится в работах [4— 7]*. Здесь дается наиболее полное изложение и исследование схем метода применительно к трансзвуковым задачам аэродинамики, течениям со „вдувом“ струи в основной поток и т. п. Такой подход дал возможность получить сложные картины плоских и осесимметричных тел различной формы в широком диапазоне изменения скоростей трансзвукового режима (от чисто дозвуковых до сверхзвуковых режимов, включая „закритическое“ обтекание и переход через скорость звука). Рассматриваются также внутренние течения со сложной конфигурацией скачков уплотнения, дифракционные задачи и др. Подчеркнем, что использование метода „крупных частиц“ позволило рассмотреть весь этот класс сложных течений для плоских и осесимметричных задач с единых позиций, с помощью одного численного подхода.
Для оценки надежности полученных результатов, постановки задачи, граничных условий и т. п. вычисления здесь проводились на различных расчетных сетках; имела место проверка выполнения законов сохранения; результаты сравнивались с имеющимися отдельными расчетами по другим схемам, а также с экспериментальными и аналитическими данными. В этом отношении полезной оказалась приведенная ниже асимптотика звукового обтекания, построенная для плоского и осесимметричного случаев. Практически везде наблюдалось достаточно удовлетворительное согласие.
Разработка метода „крупных частиц“ проводилась в ВЦ АН СССР с 1965 г. Результаты исследований в этом направлении докладывались на многочисленных научных семинарах и конференциях и, в частности, на III Всесоюзном съезде по теоретической и прикладной механике (Москва, январь 1968); на I (Новосибирск, август 1969), II (Беркли, США, октябрь 1970), III (Париж, Франция, июль 1972) и IV (Булдер, США, июнь 1974) Международных конференциях по численным методам в газовой динамике [9, 10] и др.
Из других исследований в этой области, использующих зйле-рово-лагранжевый подход, следует отметить схему РЫС-метода [11], разработанную Р. Джентри, Р. Мартином и Б. Дали для расчета высокоскоростных течений газа, а также работы [12—19].
1.3. Дадим вначале формальное описание нестационарного метода „крупных частиц“. Основная его идея состоит в расщеплении исходной системы дифференциальных уравнений по физическим процессам. Моделируемая среда заменяется здесь системой из жидких частиц, совпадающих в данный момент времени с ячейками эйлеровой сетки. Стационарное решение задачи получается в результате установления, поэтому весь процесс вычислений состоит из
* См. также монографию „Численное исследование современных задач газовой динамики“ (под редакцией О. М. Белоцерковского). М., .Наука", 1974. (Глава 11. Ю. М. Давыдов. Исследование трансзвуковых и сверхзвуковых течений методом .крупных частиц“.)
многократного повторения шагов по времени. Расчет каждого временного шага (вычислительного цикла) разбивается на три этапа:
I — „эйлеров“ этап, когда пренебрегаем всеми эффектами, связанными с перемещением жидкости (потока массы через границы ячеек нет); здесь на фиксированной эйлеровой сетке определяются промежуточные значения искомых параметров потока (и, V, £);
II — „лагранжев“ этап, где вычисляется плотность потока массы при движении жидкости через границы эйлеровых ячеек;
III — заключительный этап — определяются „окончательные“ значения параметров потока ф(м> Е, р) на основе законов сохранения массы, импульса и энергии для каждой ячейки и всей системы в целом.
По существу, как это и принято в методах „частиц“, на I этапе рассматривается за время А Ь изменение импульса и энергии лаг-ранжева элементарного объема жидкости („крупной частицы“), заключенного внутри эйлеровой ячейки (при этом граница объема смещается относительно начального расположения).
На II этапе моделируется движение частиц через границы эйлеровых ячеек и происходит перераспределение частиц по пространству; на III этапе лагранжев объем возвращается в первоначальное положение и происходит перераспределение массы, импульса и энергии по пространству (здесь определяется за время Д Ь изменение параметров потока в элементарной эйлеровой ячейке, полученной возвращением лагранжева объема в исходное положение). Счет фактически ведется в лагранжевых координатах с последующим пересчетом (интерполяцией) на эйлерову расчетную сетку.
Таким образом, эволюция всей системы за время АЬ осуществляется здесь путем следующего расщепления: вначале изучается изменение внутреннего состояния подсистем, находящихся в ячейках — „крупных частицах“, в предположении их „замороженности“ или неподвижности (эйлеров этап), а затем рассматривается смещение всех частиц, пропорционально их скорости и Д Ь без изменения внутреннего состояния подсистем с последующим пересчетом расчетной сетки в начальное состояние (лагранжев и заключительный этапы).
Рассмотрим движение идеального сжимаемого газа. В качестве исходных можно взять либо дифференциальные уравнения Эйлера в дивергентном виде (уравнения неразрывности, импульса и энергии)
+ (Иу(р!Г) = 0,
О I
^+ (11у(рдГ) + |^ = 0,
д Ь их
+ Шу (р®#)+4^- = 0,
о I о у
+ («V (р Е Щ + сНу (р Щ = О,
<? I
либо уравнения движения в интегральном виде в форме законов сохранения. Заметим, что при наличии скачков уплотнения в области интегрирования указанная форма записи исходных уравнений
(1.1)
в виде законов сохранения является наилучшей. Для замыкания этой системы используется уравнение состояния
Р = Р(?у А ^Е-Ш212. (1.2)
Обозначения здесь общепринятые: 1^ — скорость (и, V— составляю щие вдоль осей х, у соответственно); р — плотность; J и Е — удель’ ные внутренняя и полная энергии соответственно.
В
и Е N к
щ ///А
Фиг. 1.1
В данной работе рассматривается случай двумерного течения» хотя возможно обобщение указанного подхода и на случай трех пространственных переменных. Наложим на область интегрирования эйлерову (фиксированную в пространстве) расчетную сетку (фиг. 1.1), состоящую из прямоугольных ячеек со сторонами Ах, Ау (А г, А г—в цилиндрической системе координат).
На I этапе („эйлеровом“) изменяются лишь величины, относящиеся к ячейке в целом, а жидкость предполагается моментально
заторможенной. Поэтому конвективные члены вида фу (<р р Ш), где <р = (1, и, V, Е), соответствующие эффектам перемещения, в уравнениях (1.1) откидываются. Тогда из уравнения неразрывности следует, что поле плотности будет „заморожено“ и исходная система (1.1) здесь примет вид
Р-тт + ^-О. р-^- + <и*(р®)-0. (1.1')
д Ь ' д у ' ' д t
Здесь, как уже отмечалось, по существу
рассматривается за время Д^ изменение импульса и энергии лагранжева элементарного объема жидкости, заключенного внутри эйлеровой ячейки.
Простейшие явные конечно-разностные аппроксимации, (симметричные разности) первого порядка точности по времени и пространству для (1.1') будут выглядеть так
Иг, / == и
I, I
Д х
р* /
VI, / = V.
Г," А ,
Рг, /-Н/2 /—1/2 Ы
(1-3)
Р? .
Величины с дробными индексами, относящиеся к границам ячеек, равны полусумме значений соответствующих параметров в соседних ячейках. Например, р%+11% } = (рпи /+Р?+1.у)/2-
В ряде случаев целесообразно использование схем и более высокого порядка точности. Так, для улучшения устойчивости вычислений на этом шаге рассматривались также схемы метода интегральных соотношений, где аппроксимации подынтегральных функций проводились по N лучам (N=3, 4, 5) [20]. Производились также расчеты по неконсервативным схемам с использованием уравнения для внутренней энергии У.
Несмотря на то что приведенные численные схемы первого этапа обладают, вообще говоря, недостаточной устойчивостью, из-за простоты вычислений возможно использование аппроксимаций (1.3), так как на последующих этапах проводится соответствующая компенсация неустойчивости этого шага и, в целом, оператор полного шага оказывается устойчивым [4 — 6].
На II этапе („лагранжевом“) моделируется движение потока массы А Мп через границы эйлеровых ячеек и происходит перераспределение массы по пространству. При этом полагаем, что вся масса переносится только за счет нормальной к границе составляющей скорости. Так, например,
А Л^+1/'2_ у — (р )(■ + 1/2, / ( и” )г+1/2, / Ау А
Знак „{ )“ обозначает значения р и и на границе ячейки. Выбор этих величин имеет важное значение, так как он влияет на устойчивость и точность вычислений (для всех видов записи АМп характерен учет направления потока на данной границе) [2 — 7].
а) Значение АМп можно определять по формулам второго порядка точности. Например, рассчитывается следующим
образом:
если и1 у + и"+1, / > 0
и
й?+Ь% , = и?.,+ (ди/дх),, ,-±£- = и*; + .а"+Ч-И"-Ь/ > 0,
2 4
то
Ьм?*т.1 =(*’.!+ + уум- (1.4)
если и” ] + И?+1. ) < 0
и
«?+1А, == и?+г, , - (дфх)1+ь J ^ = «?+1, ,--< о,
& 4
то
р;+-2'Ь' ¡. ^ АуАг. (1.4')
В противном случае ДЖ?+1/2 / = 0.
Использование в (1.4), (1.4') скоростей на предыдущем временном слое I приводит здесь к неустойчивости [6]. Поэтому в данном случае предопределен порядок вычислений: сначала на эйлеровом этапе вычисляются и и V, а потом они используются на лагранжевом этапе для расчета А М. Такая последовательность вычислений дает возможность проводить устойчивый счет в целом лишь при условии устойчивости первого этапа, для чего в (1.3) вместо термодинамического давления р необходимо ввести
Р=Р + д.
АМп
1+1/2. Г-
«¿ + 1, /----------
п
Рг+1. У
Здесь д— нелинейное искусственное вязкостное давление, например, типа, введенного Ландшоффом [2, 11], д—— ВСрк-—^- при
4^<0 и ? = 0 при -^->•0, действие которого проявляется лишь
в волнах сжатия для сглаживания образующихся разрывов (В — постоянный коэффициент, С — скорость звука, /г — размер ячейки в направлении 5 — х или у). Рассматривались также другие виды искусственной вязкости: линейная, квадратичная, без „зануления“ коэффициента вязкости на волнах разрежения и т. п.
б) Значение Д Ма можно определять по формулам первого порядка точности*, тогда
р» Л1+..а1+Ы- А у Д*, если II / + иР+1, / > 0;
2 ' {1'5) р^г.цЬ + цж./ если и1/ + и?+и<0.
Расчеты и аналитический анализ, приведенные в работах [4 — 6], показали, что формулы (1.5) позволяют проводить устойчивый счет без введения явных членов искусственной вязкости, т. е. при ^ = 0. Устойчивость вычислений обеспечивается при этом внутренней структурой разностной схемы — наличием аппроксимационной вязкости. Здесь также возможно без потери устойчивости использование значений скорости на предыдущем временном слое (применение И И V лишь несколько ускоряет сходимость), поэтому в данном случае мы можем сначала проводить II этап, а потом I.
Возможна, естественно, комбинация обоих подходов [например, при расчете областей со сложной структурой целесообразно использовать формулы (1.5), а для зон „гладких“ течений — аппроксимации второго порядка точности (1.4), (1.4') и т. п.].
в) Возможно, конечно, использовать здесь и центральные разности
Д Мп — о" ип А\> А1 “ ¿±1/2, / ".’±1, /И/±1, / У
или определять Д Мп без учета направления потока:
(р иУ{± 1/2, / = (р”, / и1±1, / + ?"± 1, / и1, /)•
Однако такие подходы приводят к заметной неустойчивости вычислений всего цикла (например, в случае использования центральных разностей неустойчивость проявляется на волнах разрежения и т. д.).
г) Если использовать дискретную модель сплошной среды и рассматривать в данной ячейке совокупность частиц фиксированной массы т, то в этом случае АМП определяется как алгебраическая сумма масс всех А-частиц, пересекших за время Д£ данную границу
д мп=2 т1-
к
В этом случае мы получаем консервативный метод „частиц в ячейках“, который позволяет изучать движение многокомпонентных сред (проводились и такие расчеты).
* Аналогичная схема используется во Р1ЛС-методе [11].
Как правило, при конкретных расчетах в зависимости от характера рассматриваемого течения использовались различные формы аппроксимаций I и II этапов.
И, наконец, на III этапе (заключительном) происходит перераспределение массы, импульса и энергии по пространству, проводится регуляризация сетки (пересчет сдвинувшейся сетки в прежнее положение) и определяются окончательные поля эйлеровых параметров потока в момент времени £л+1 = -¡- Д£ („выбираются“ все не-
вязки уравнений). Как уже отмечалось, уравнения этого этапа представляют собой законы сохранения массы М, импульса Р и полной энергии Е, записанные в разностной форме
Рп+1=Рп + ^Гп^ (1.6)
где Р = (М, Р, Е).
Эти уравнения утверждают, в частности, что внутри поля течения нет источников и стоков М, Р и Е, а их изменение за время М осуществляется только за счет взаимодействия на внешней границе области течения. При этом предполагается, что потоки массы через границы ячеек Д М”, определяемые на II этапе, несут с собой промежуточные значения скорости и удельной энергии, вычисленные на I этапе (величины Д Мп играют здесь роль весовых функций).
Исходя из этого, окончательные значения параметров потока Р, X = (м, V, Е) на следующем временном слое вычисляются по формулам
р»+.‘ = р» Х'1+1 = X- I + ^ *к А-М-, (1.7)
где суммирование проводится по всем границам ячейки.
Можно показать [4], что в данной разностной схеме внутри области интегрирования имеет место строгое выполнение законов сохранения массы, импульса и энергии. Общее изменение, например,
величин Е и Р за время М равно сумме их изменений на I и III этапах (законы сохранения для всей сеточной области являются алгебраическим следствием разностных уравнений), поэтому в целом разностная схема является дивергентно-консервативной*, т. е. уравнения (1.7) удовлетворяют (1.6). Проверка выполнения законов сохранения для схемы представляется важным, так как известно, что расчетная схема дает наиболее точные результаты, когда она строго сохраняет те же величины, что имеет место и в рассматриваемом физическом процессе.
Описание алгоритма метода „крупных частиц“ дано здесь для „целых“ ячеек — внутренних ячеек поля течения (со всех сторон окруженных жидкостью) и граничных ячеек, прилегающих к обтекаемому телу так, что его контур совпадает с их границами. В случае же обтекания тел произвольной формы задача значительно усложняется: контур тела пересекает прямоугольную эйлерову сетку, поэтому необходимо в рассмотрение вводить „дробные“ ячейки [2, 6, 7, 11]. Заметим, что использование здесь дробных ячеек дает
* Анализ разностных схем, например, с точки зрения исследования их устойчивости, показывает, что целесообразно использовать здесь энергетическое равенство для полной энергии Е,
гораздо больший эффект, чем другие способы задания криволинейной границы тела (например, расчет в координатах 5, п, где 5 — длина дуги вдоль контура тела, л —нормаль к телу и др.).
1.4. Обсудим теперь постановку краевых условий задачи. В качестве начальных данных обычно использовались параметры невозмущенного потока.
Граничные условия задачи ставились следующим образом (см. фиг. 1.1): на левой границе АВ использовались условия в набегающем потоке газа; на плоскости (или оси) симметрии АО, Ю — условия симметрии потока; на теле ОЕКЬ — обычные условия на твердой стенке (условия „непротекания“ или „прилипания“); на верхней ВС и правой СО, „открытых“ границах области, проводилась экстраполяция параметров течения за рассматриваемую область. Чтобы не нарушать единообразия вычислений и не применять особые формулы для граничных ячеек, вдоль всех границ вводились слои „фиктивных“ ячеек, куда и засылались параметры из внутренних (смежных) ячеек. Число таких слоев определялось порядком разностной схемы (для схемы первого порядка — один слой и т. д.).
Следует различать два рода границ — твердая стенка (или ось симметрии) и „открытая“ граница расчетной области.
В первом случае, например при условии „непротекания“, нормальная к стенке компонента скорости меняет знак, а остальные параметры потока берутся без изменения (при условии „прилипания“ обе компоненты скорости меняют знак)*.
Через „открытые“ границы области жидкость может втекать или вытекать и здесь должны быть обеспечены некоторые условия непрерывности движения. Пусть жидкость втекает в прямоугольную сетку с левой стороны, тогда здесь и задаются параметры набегающего потока. На остальных открытых границах области проводится экстраполяция параметров потока „изнутри“, т. е. в фиктивный слой переносятся значения параметров из ближайшего (к границе) слоя (экстраполяция нулевого порядка). Возможна и более сложная постановка граничных условий или применение более точных формул экстраполяции (линейная, квадратичная и др.).
При сверхзвуковом течении характер граничных условий (экстраполяция изнутри) не вносит каких-либо осложнений, так как возмущения, распространяющиеся со скоростью звука, „сносятся“ потоком. Сложнее при дозвуковых или трансзвуковых течениях — здесь необходимо позаботиться, чтобы влияние границ было небольшим. Естественно, что внешняя граница облисти интегрирования должна выбираться „достаточно далеко“ от источников возмущения, тогда методы экстраполяции потока „во вне“ возможны. Таким образом, основной принцип постановки краевых условий заключается в том, что через открытые границы области не должны проникать сколько-нибудь заметные возмущения внутрь расчетной зоны.
1.5. Детальное исследование полученных разностных схем метода „крупных частиц“ (вопросы аппроксимации, устойчивости, образования механизма диссипации и др.)‘, обоснование постановки задачи, краевых условий и т. п. содержится в работах [4—6, 22].
* Для обоих типов граничных условий из-за „сглаживающего* влияния ап-проксимационной вязкости на поверхности тела реализуются значения скоростей, близкие к нулю (условия „прилипания“ [21]).
Известно, что при изучении в рамках идеального газа задач с разрывами приходится вводить понятие обобщенного решения. Последнее исходит из единообразного описания гидродинамического течения, при этом соответствующие разностные схемы называются однородными. Возможны два способа такого единообразного описания течений — использование интегральных законов сохранения и рассмотрение уравнений газовой динамики с диссипативными членами.
Обобщенное решение, определяемое интегральными законами сохранения в эйлеровых или лагранжевых координатах, является единообразным, поскольку уравнения газовой динамики и условия совместности на разрывах являются следствием законов сохранения. Во втором подходе рассматриваются уравнения с диссипацией (например, в уравнения газовой динамики искусственно вводятся члены с малым параметром — »вязкостью“— при старших производных), что приводит к сглаживанию разрывов. За обобщенное решение здесь принимается предел решений полученных таким образом параболических уравнений, когда коэффициент при „вязких“ членах устремляется к нулю [23, 24].
Как следует из сказанного выше, в рассматриваемом подходе применяются однородные разностные схемы, построенные по второму способу, позволяющие в рамках идеального газа проводить по единым алгоритмам „сквозной“ счет как в областях гладких течений, так и на разрывах. При этом аппроксимирующие (1.1) разностные уравнения должны обладать определенными свойствами, чтобы допускать решения, аналогичные скачку уплотнения. Достигается это путем использования соответствующих диссипативных конечно-разностных однородных схем и аппроксимацией дифференциальных уравнений (1.1), записанных в дивергентной форме.
В качестве исходных уравнений рассматривались уравнения газовой динамики для невязкого газа (1.1), однако нашей разностной схеме присущи и „вязкостные эффекты“, которые обеспечивают единообразие расчетных формул как в области гладких решений, так и на разрывах. Вязкостные эффекты могут быть обусловлены, во-первых, введением в схему явного члена с искусственной вязкостью („вязкостного давления“) ц, и, во-вторых, наличием собственно схемной (аппроксимационной) вязкости, зависящей от внутренней структуры конечно-разностных уравнений.
Приведенные выше разностные схемы являются многослойными, а разностные уравнения — существенно нелинейными с переменными коэффициентами. Это делает практически невозможным использование здесь метода Фурье для анализа устойчивости всей схемы в целом, основанного, как известно, на исследовании линейных уравнений с постоянными коэффициентами без учета градиентов.
Исследование вопросов аппроксимации, образования механизма диссипации (для возможности проведения „сквозного“ счета) и устойчивости нелинейных разностных схем с переменными коэффициентами проводится здесь с помощью разработанной в [4—6, 22] методики последовательного рассмотрения так называемых дифференциальных приближений,полученных из разложения всех членов разностных уравнений в ряды Тейлора по сеточным параметрам Ах, А1,.,., например,
и1+/ — и{х, у, ¿ + «?+:, / == и (•* + Ах> У< *)•
Члены нулеього (низшего) порядка представляют при этом исходные дифференциальные уравнения, учет в разложении членов более высокого (£ + 1)-го порядка (к = 0, 1, 2,...) позволяет определить дополнительные (к уравнениям) члены разложения, а следовательно, и структуру дифференциального приближения к-то порядка (£-го дифференциального приближения). Вид дифференциальных приближений зависит от структуры аппроксимации.
Для одномерных квазилинейных уравнений гиперболического типа исследования устойчивости разностных схем с помощью первого дифференциального приближения были проведены в работе [25]. Для нелинейных уравнений строгого математического обоснования пока нет, но, как показал Херт [3], такой подход используется и там. Рассматривая последовательно нулевое, первое и второе дифференциальные приближения, можно изучить все основные свойства разностных схем, полученные в методе „крупных частиц“ [4—6, 22, 26].
Основной смысл введения дифференциальных приближений заключается в том, что их рассмотрение позволяет весьма наглядно проводить исследование основных свойств разностной схемы. При этом можно показать на модельных примерах, что дифференциальное приближение сглаживает начальные данные и разрывы примерно так же, как и соответствующая разностная схема [24], а оператор решения разностного уравнения асимптотически совпадает с оператором решения соответствующего дифференциального приближения [23].
Исследуя порядок остаточных членов нулевого приближения, можно показать [22], что во всей области течения (как во внутренних точках, так и на границе с твердым телом) будем иметь для приведенных выше схем первый порядок аппроксимации по времени и первый-второй по пространственным переменным (в зависимости от вида формул для вычисления АМп). Надо иметь в виду, что эти оценки носят асимптотический характер в предположении достаточно гладких решений и малых шагов сетки, так что для нелинейных уравнений и реальных сеток (особенно при наличии сильных разрывов) использование схем формально более высокого порядка кажется нецелесообразным.
Выписав первое дифференциальное приближение схемы, можно определить характер вязкостных эффектов нашей схемы. Так, например, при использовании формул второго порядка точности (1.4), (1.4') оно имеет следующий вид в одномерном случае:
др , дри_ _ п .
+ дх и + • • • дщ , д(р + риЦ дq , д г да ]
дх дх дху дх) ' ' ' ’
0(Д*. Ах2). (1.8)
При использовании формул первого порядка точности (1.5) будем иметь
аР , дРи __д(ду\ с)£ ' дх дх\ дх) ^ '
др“ д (р + рцз) _ _д£_ _д_ ! дри\ , д*и_ .
дх дх дх ч дх) * ^ дх2 * ' '
д$Е . д г . . г-ч! дди . д2 (р£) , де д$Е
где е = | и I Ах/2.
+0(А^Ах%{\Я')
Слева в (1.8) и (1.8') получили точные выражения исходных дифференциальных уравнений, а справа выписаны члены, которые являются следствием наличия „вязкостных“ эффектов в разностных уравнениях, т. е. разностная схема является диссипативной. Важно отметить, что порядок уравнений дифференциальных приближений, описывающих движение, выше, чем в исходной системе. Члены с д возникают из-за явного введения в систему искусственной вязкости, а члены с е (схемная или аппроксимационная вязкость) определяются внутренней структурой приближения и возникают в результате аппроксимации точных дифференциальных уравнений конечно-разностными („ошибки разложения“)*. Таким образом, дифференциальные приближения (1.8) и (1.8') имеют различную структуру членов с аппроксимационной вязкостью, что и предопределяет в этих двух случаях разные механизмы диссипации.
Следует отметить, что при наличии разрывов построение разностных схем с диссипацией желательно, вообще говоря, проводить без введения явных членов с искусственной вязкостью д. Дело в том, что для д ф 0 искомое обобщенное решение должно получаться здесь асимптотически при д 0. Однако такой предельный переход при расчетах сложных задач провести практически оказывается невозможно, так как заметное уменьшение величины д приводит обычно к потере вычислительной устойчивости, а точную оценку границ изменения параметров искусственной вязкости получить бывает очень трудно. Кроме того, в этом случае значительно затрудняется построение алгоритма расчета для течений с двумя или тремя пространственными переменными, и особенно для тел с криволинейной образующей (способы введения искусственной вязкости в целых и дробных ячейках различны [6,7]); при дф 0 сильно „размываются“ области больших градиентов (например, ударные волны) и т. п. Таким образом, схемы, обладающие лишь внутренней диссипацией (схемной вязкостью в) и позволяющие проводить устойчивый счет без введения явных членов с искусственной вязкостью, являются, на наш взгляд, более естественными и перспективными при построении обобщенных решений многомерных задач газовой динамики.
Легко видеть, что при измельчении сетки (А(, Ах 0) значения е -»• 0 (в пределе получаем обобщенное решение) и уравнения дифференциального приближения переходят в точную систему исходных уравнений. При конкретных вычислениях (из-за конечности Ах, At>...) в разностной схеме даже при д = 0 неявным образом всегда присутствуют члены, которые содержат а и, в свою очередь, аналогичны диссипативным членам уравнений Навье—Стокса. При этом роль коэффициента реальной вязкости умод здесь играет коэффициент схемной вязкости е, зависящей от локальной скорости потока и размера разностной сетки. Подчеркнем здесь еще раз, что для однородных схем сквозного счета наличие аппроксимационной вязкости в разностных уравнениях обязательно, так как именно ее присутствие обеспечивает единообразие расчетных формул в областях гладких решений и на разрывах.
* Более строгое математическое определение аппроксимационной вязкости содержится, например, в [23].
В двумерном случае при использовании (1.4), (1.4') из уравнений импульса следует, что схемная вязкость (при <7 = 0) имеет вид тензора:
Д =
1 иАх да дх + ъАу ди ~ду~
-ГР иАх дь дх иДу ди ~3у
1 Р#д
2
где Дг = Дхг + Ду/.
Отсюда видно, между прочим, что схемная вязкость, благодаря
присутствию векторов Дг и не обладает инвариантностью относительно преобразований Галилея, причем практически проявляется она лишь в зонах больших градиентов: на ударной волне, которая „размазывается“ на несколько счетных ячеек, у поверхности тела, где образуется достаточно широкий пограничный слой, и т. п. При этом коэффициент схемной вязкости (а следовательно, и ширина получаемой „размазанной“ ударной волны) зависит от величины локальной скорости потока и размера ячеек.В областях же гладкого течения, где градиенты параметров потока сравнительно малы, влияние схемной вязкости незначительно.
Задача указанного механизма диссипации — сгладить, „размазать“ сильные разрывы, что и позволяет проводить „сквозной“ счет в рамках идеального газа по единым численным алгоритмам. При корректном построении алгоритма диссипативные схемы должны обеспечивать, как уже отмечалось, решение, аналогичное по своей структуре скачку уплотнения. На границах „размыва“ должны выполняться с определенной точностью соответствующие условия на разрывах (например, условия Гюгонио для ударных волн), причем с увеличением интенсивности разрыва „толщина“ скачка должна уменьшаться, а при г -> 0 зона „размыва“ должна асимптотически переходить в поверхность разрыва.
Указанным свойством обладают гиперболические схемы первого порядка точности, когда е — Ах. Если бы удалось построить для расчета сложных течений газа конечно-разностную схему с аппрок-симационной вязкостью, пропорциональной величине (Ах)2, то скачки в расчетах получились бы более тонкими. Однако в ряде случаев гиперболические схемы второго порядка точности являются не диссипативными, а дисперсионными (когда различные фурье-компоненты точного решения разностного уравнения распространяются с различными скоростями, но без изменения амплитуды). По этой причине они оказались менее пригодными для практических расчетов, например при исследовании трансзвуковых течений со скачками уплотнения [21, 27] и др.
Из анализа первого и второго дифференциальных приближений следует, что необходимые условия вычислительной устойчивости полной схемы (условие а-параболичности) можно получить, оценивая знак „коэффициента диффузии“ аг у диссипативных членов, содержащих частные производные второго порядка по пространственным переменным [3—6, 22]. Можно показать, что при отрицательном значении а, уравнение дифференциального приближения допускает экспоненциально возрастающее по времени неустойчивое решение, причем в линейном случае результаты анализа устойчивости с помощью дифференциального приближения и метода Фурье полностью совпадают [3,4,26]. Таким образом, требование <*(>-0
как необходимое условие вычислительной устойчивости, и является определяющим при построении соответствующих диссипативноустойчивых схем *.
Рассмотрим, какой вклад в неустойчивость дают разные формы записи уравнения неразрывности (II этап вычислений), считая, что уравнения импульса и энергии устойчивы.
Если ДМа определяется по формулам второго порядка точности (1.4), (1.4'), то, разлагая приведенные разностные схемы в ряд Тейлора и удерживая члены, содержащие д2р[дхя, получим в одномерном случае
др I дрЦ ___ Д , 2 , , , д,
В случае вычисления ДМ" по формулам первого порядка точности (5) имеем
^ л ЙВ _ д. + (Д£. | и | + (И2 _ С2) _ ^ ^ , (1.9')
1 дх 1 \ 2 1 1 2 4 дх } дх?
где Д15 Д^ — члены первого дифференциального приближения, пропорциональные Ах и содержащие первые производные. Имеем Д*» 0,071, Д*=к 0,0071, р_оо = 1, £/_оо=1.
При практических расчетах, когда возникают ударные волны, контактные разрывы и волны разрежения,
Р]«1= 1, ^ <0,3, -% Ах<2.
Отсюда следует, что в (1.9') коэффициент при д2р/дх2 положителен, а в (1.9) — отрицателен, т. е. схема (1.4)—(1.4') второго порядка точности для АМп допускает быстрорастущие решения и является неустойчивой при счете, в то время как схема первого порядка (1.5) устойчива. Это иллюстрирует также фиг. 1.2, где
* Следует не путать физическую и вычислительную неустойчивости. Если физическая неустойчивость порождается самим явлением, то вычислительная — в результате флуктуаций, возникающих при счете из-за наличия в схеме диффузионных (второго порядка) ошибок округления, которые могут возрастать под действием конвективных членов.
показаны профили плотности на линиях г = const перед цилиндрическим торцом для этих двух случаев [сплошные линии — формулы (1.4), (1.4') и (1-9), штриховые — (1.5) и (1.9')]. Видно, например, что схема (1.9') допускает устойчивый счет, в то время как для (1.9) сразу за ударной волной возникает существенная неустойчивость. В последнем случае устойчивость вычислений достигается введением в схему явного члена вязкостного давления д. Можно показать также, что использование для вычисления формул первого порядка точности (1.5) делает всю нашу разностную схему устойчивой даже без введения члена с q [22].
Итак, была поставлена задача получить схему „сквозного“ счета, позволяющую в рамках идеального газа проводить непрерывным образом расчет и поверхностей сильного разрыва. Показано, что выполнение в конечно-разностных схемах необходимых условий устойчивости для широкого класса газодинамических задач обеспечивается [при использовании формул (1.5)], как правило, лишь схемной вязкостью s (диссипативно-устойчивые схемы). При этом схема имеет решение, аналогичное по своей структуре скачку уплотнения, что достигается наличием соответствующего механизма диссипации в разностных уравнениях и аппроксимацией дивергентной системы (дивергенто-консервативные схемы). Оказалось важным также (это следует из структуры дифференциальных приближений) использование здесь уравнения для полной (а не внутренней) энергии Е. Как показали многочисленные расчеты двумерных (по пространству) газодинамических задач, эти схемы обеспечивают устойчивый „сквозной“ счет с постоянными сеточными параметрами практически во всей области интегрирования (при этом, например, ударная волна „размазывается“ на две-три счетные ячейки и т. п.).
В работах [4—6, 21] проводились аналитические и численные оценки и обоснования постановки краевых условий задачи, а также имел место контроль процесса вычислений и полученных данных („склейка“ с асимптотикой; расчет на разных сетках аппроксимации; „продолжение“ полей течений; сравнение с расчетами по другим методикам, экспериментом и т. п.). Некоторые данные из области этих методических исследований будут приведены ниже.
1.6. Из самого характера построения расчетной схемы следует, что здесь, по существу, решается полная система нестационарных уравнений газовой динамики, причем каждый вычислительный цикл представляет собой законченный процесс расчета данного временного интервала. При этом удовлетворяются все исходные нестационарные уравнения, граничные условия задачи и определяется действительное течение жидкости в соответствующий момент времени. Таким образом, метод „крупных частиц“ позволяет получить характеристики нестационарных течений газа, а также с помощью процесса установления — их стационарные значения. Использование такого подхода кажется особенно целесообразным в задачах, когда имеет место разнохарактерное развитие по времени физического явления. Например, при изучении трансзвуковых потоков газа, расчете обтекания тел конечных размеров и др., где при довольно быстром установлении большей части поля структура течения в местных сверхзвуковых зонах, областях срыва и т, п. формируется сравнительно медленно.
В отличие от методик, применяемых зарубежными авторами [1, 2, 11, 28], метод „крупных частиц“, разработка которого проводилась независимо с 1965 г., широко используется для системати-
ческих расчетов задач обтекания тел сжимаемым газом в широком диапазоне изменения начальных данных (по единому алгоритму рассчитываются, как уже отмечалось, до-, транс- и сверхзвуковые задачи, области срыва, местные сверхзвуковые зоны и т. п.). В этом методе рассматриваются дивергентные формы исходных и разностных уравнений; используется энергетическое равенство для полной энергии Е; по-новому трактуются этапы вычислительного цикла, граничные условия задачи и др. Так, например, для определения АМп здесь применяются разностные представления различных порядков точности (в зависимости от типа течения); на заключительном этапе вводится дополнительный пересчет плотности (что способствует устранению флуктуаций и позволяет получить удовлетворительные результаты при сравнительно небольшом числе узлов расчетной сетки); имеет место строгое выполнение законов сохранения массы, импульса и энергии и т. п.
Все это дает возможность получить в методе „крупных частиц“, в отличие от упомянутых выше подходов, дивергентно-консервативные и диссипативно-устойчивые схемы, позволяющие для большого класса газодинамических задач проводить устойчивый счет без введения явных членов с искусственной вязкостью. Это кажется особенно важным при изучении обтекания тел с криволинейной образующей, так как способы введения явных членов искусственной вязкости в целых и дробных ячейках различны. Кроме того, путем изменения лишь второго этапа вычислительного цикла отсюда можно получить консервативный мотод „частиц в ячейках“, так что алгоритм расчета является достаточно общим.
В следующих статьях этого цикла будут приведены асимптотики звуковых течений, методические расчеты, а также результаты численных исследований различных задач газовой динамики.
ЛИТЕРАТУРА
1. Evans М. W., H а г 1 о w F. Н. The particle-in-cell method for hydrodynamic calculations. Los Alamos scientific Lab., Rept. N LA-2139,
1957.
2. Rich M. A method for eulerian fluid dynamics. Los Alamos scientific Lab., Rept. N LAMS-2826, 1963.
3. Hirt C. W. Heuristic stability theorie for finite-difference equations.
J. Comput. Phys., 3, N 4, 1968.
4. Белоцерковский О. М., Давыдов Ю. М. Нестационарный метод .крупных частиц“ для решения задач внешней аэродинамики . М., Препринт ВЦ АН СССР, 1970.
5. Белоцерковский О. М., Давыдов Ю. М. Нестацио« нарный метод .крупных частиц“ для газодинамических расчетов. Ж. вычисл. матем. и матем. физ., 11, М 1, 1971.
6. Давыдов Ю. М. Метод „крупных частиц" для задач газовой динамики. Канд. диссертация. М., МФТИ, 1970.
7. Давыдов Ю. М. Расчет обтекания тел произвольной формы методом „крупных частиц“. Ж. вычисл. матем. и матем. физ., II,
Ws 4, 1971.
8. Яненко H. H., Анучина H. H., Петренко В. E., Шокин Ю. И. О методах расчета задач газовой динамики с большими деформациями. Инф. бюл. „Числ. методы мех. сплошной среды*,
1, íé 1, 1970.
9. Beloíserkovskii О. М. Numerical method of some transsonic aerodynamics problems. J. Comput Phys., 5, N 3, 1970.
10. Belotserkovskii О. М., Davydov Ju.M. Numerical approach for investigating some transsonic flows. Proceeding of the Third Inter. Conference on Numerical Methods in Fluid Dynamics (July 3—7,
1972, Paris, France), Lecture Notes in Physics, Springer—Verlag, 1973.
2—Ученые записки Jè 3
17
11. Gentry R. A., Martin R. E., Daly B. J. An eulerian differencing method for unsteady compressible flow problems, J. Comput. Phys., N 1, 1, 1966.
12. Дьяченко В. Ф. Об одном новом методе численного решения нестационарных задач газовой динамики с двумя пространственными переменными. Ж, вычисл. матем. и матем. физ., 1965, 5, № 4.
13. Нох В. СЭЛ — совместный эйлерово-лагранжев метод для расчета нестационарных двумерных задач. В сб. „Вычислительные методы в гидродинамике“. М., »Мир*, 1967,
14. А н у ч и н a H. Н. О методах расчета течений сжимаемой жидкости с большими деформациями. Инф. бюл. »Численные методы механики сплошной среды“. Новосибирск, 1, № 4, 1970.
15. Херт С. Произвольный лагранжево-эйлеров численный метод. В сб. „Численные методы в механике жидкостей*. М., „Мир“, 1973.
16. Кроули У. FLAG — свободно-лагранжев метод для численного моделирования гидродинамических течений в двух измерениях. В сб. »Численные методы в механике жидкостей“. М., „Мир“, 1973.
17. H 1 r t С. W., Cook J. L,, Butler Т. D. A lagrangian method for calculating the dynamics of an incompressible fluid with free surface. Труды секции по численным методам в газовой динамике II международного коллоквиума по газодинамике взрыва и реагирующих систем (Новосибирск, 1969 г.). М., Изд. ВЦ АН СССР, т. II, 1971.
18. Франк Р. М., Лазарус Р. Б. Смешанный метод, использующий переменные Эйлера и Лагранжа. В сб. „Вычислительные методы в гидродинамике*. М., „Мир“, 1967.
19. Харлоу Ф., Амсден Э„ Херт С. Численный расчет течений жидкости при произвольном числе Маха. В сб. „Численные методы в механике жидкостей“. М., „Мир“, 1973.
20. Давыдов Ю. М. Разработка нестационарного метода „крупных частиц“ и расчет обтекания цилиндрического торца на трансзвуковых и сверхзвуковых режимах. Труды Казанского авиационного института, 1971.
21. Белоцерковский О. М., Давыдов Ю. М. Расчет методом „крупных частиц* трансзвуковых „закритических“ режимов обтекания. Ж, вычисл. матем. и матем. физ., 13, № 1, 1973.
22. Белоцерковский О. М., Давыдов Ю. М. Исследование схем метода „крупных частиц* с помощью дифференциальных приближений. В сб. „Проблемы прикладной математики и механики“. М., .Наука“, 1971.
23. Рождественский Б. Л., Яненко H. Н. Системы квазилинейных уравнений и их приложения к газовой динамике. М., „Наука“, 1968.
24. Г о д у н о в С. К., Р я б е н ь к и й В. С. Разностные схемы. М., „Наука*, 1973.
25. Яненко H. Н., Ш о к и н Ю. И. О первом дифференциальном приближении разностных схем для гиперболических систем уравнений. Сибирский матем. журнал, 10, Jé 5, 1969.
26. Давыдов Ю. М. Об одном методе исследования устойчивости нелинейных резностных уравнений. Труды МФТИ, сер, Аэрофизика, прикладная математика. М., 1971.
27. М е р м а н, К р у п п Дж. Решение трансзвукового уравнения для потенциала с помощью смешанной системы конечно-разностных схем. В сб. „Численные методы в механике жидкостей“. М., „Мир“,
1973.
28. Harlow F. H., H f r t C. W. Recent extentions to eulerian methods for numerical fluid dynamics. Ж. вычисл. матем. и матем. физ., 12, № 3, 1972.
Рукопись поступила I3jX 1975 г.
