книги / Метод крупных частиц в газовой динамике

с твердыми стенками представляет большой практический интерес. В зависи­ мости от конфигурации канала, интенсивности взрыва и первоначального местонахождения его центра в различные моменты времени на твердой стенке могут реализоваться случаи как регулярного, так и нерегулярного отражения ударных волн.

В работе Ю. М. Давыдова и С. П. Шевырева [198] расчеты такого рода взрывных задач проводились с помощью метода крупных частиц. Применение этого подхода граничит здесь с проведением численного эксперимента.

Кратко о постановке задачи. Пусть в начальный момент времени в ка­ нале ABCD (рис. 10.8, а) шириной 20Ау и длиной 58Дх находится покоящийся

т

3

F

А \УУ/ЛУ7/>7/7////////у£V7/7/77777/7//777/7/,3

а)

б)

вышенное давление с энергоперепадом E JE \. В последующие моменты времени от указанного центра будет распространяться взрывная волна. При взаимо­ действии ударной волны со стенкой в зависимости от ее конфигурации и ин­ тенсивности могут реализоваться случаи как регулярного, так и нерегуляр­ ного отражения.

Поскольку задача имеет плоскость симметрии KL и ось (в осесимметрич­ ном случае v = l) или плоскость (в плоском случае, v=0) симметрии FE , ее до­ статочно решать в одной четверти расчетного поля, например, в области FBKO, заштрихованной на рис. 10.8, а. Ряд расчетов для контроля правильности постановки граничных условий проводился в половине расчетного поля F B C E , при этом на верхней границе ВС ставились условия для твердой стенки, на

оси (плос ости) симметрии FE — условия симметрии потока, на открытых границах ts t и Lh — условие свободного вытекания — поток экстраполиро­ вался вовне. Результаты расчетов показали полную симметрию картины те-

чения относительно плоскости /CL, т. е. подтвердили правильность постановки граничных условии.

На следующих рис. 10.8—10.10 приводятся результаты решения задачи, где изображены изобары p=const (в невозмущенном газе р = 0,4). На рис. 10.8, 10.9 даны поля давления для осесимметричной задачи распространения сфе­ рической (в начальные моменты времени) взрывной волны в цилиндрической трубе. На рис. 10.8, б—г приводятся результаты расчетов в последователь­ ные моменты времени *=0,50 (рис. 10.8,6), *=0,75 (рис. 10.8, в) и *=1,00

Рис. 10.9. Распространение сферической взрывной еолны в цилиндрической трубе [198] при £ 2/^ I = Ю000. а), б) линии p=const при ^=0,25 и 0,50 соответственно.

(рис. 10.8, г) при первоначальном энергоперепаде jE2/£i=500; на рис. 10.9, а ,б — в моменты *=0,25 (рис. 10.9, а) и *=0,50 (рис. 10.9,6) при E J E X= 10 000. В обоих случаях шаг по времени Д*=0,001. На рис. 10.10 изображены линии уровня давления в момент *=0,50 для плоской задачи при распространении

ЛЯЦИРасчеты проводились на ЭВМ БЭСМ-6, обработка велась на ЭВМ М-220. 5 Г С Романов и А. В. Тетерев [199] исследовали с помощью метода круп­ ных частиц задачу о нестационарном осесимметричном движении газа при сосредоточенном энерговыделении у его поверхности. Изучались особенности процесса передачи энергии в полупространство, заполненное холодным плот-

НЫМ Такой' тип течений имеет для времени (U - характерное время за­ дачи, определяемое начальными условиями) предельный режим - автомодель-

ное движение второго рода [200]. Выход на этот режим и само автомодельное движение изучены, как отмечается в [199], в одномерном плоском случае; для двумерных движений имеются лишь частные решения. В работе [199]^методом крупных частиц получено достаточно общее численное решение такой задачи в двумерном случае — рассматривается динамика развития осесимметричного

t=o,i

Ь,09-70~2

6,09-70~г °>7°3/ 0,703/

0 ,2 5 8 ,

г> 10

\ г

Рис. I0 .ll. Пространственное распределение давления при сосредо­ точенном энерговыделении у поверхности (вариант «20» [199]).

движения при сосредоточенном энерговыделении у поверхности газа при раз­ личных начальных условиях.

Несколько слов о постановке задачи (по [199]). В момент /=*0 с границей полупространства z*^0, заполненного холодным газом с плотностью р0, при­ водится в контакт цилиндр (ударник) с высотой h0, радиусом г0 и энергией Е0.

вижный цилиндр полностью погружен в холодный газ) при £=0 вся энергия Е0 сосредоточена в тепловой ее компоненте.

Расчеты проводились для трех соотношений высоты и радиуса цилиндра:

Во всех вариантах использовались безразмерные величины: начальные плот­ ности цилиндра рю= 1» плотность среды р2о=0,001, энергия £ 0=7861.

При численных расчетах помимо полей газодинамических переменных оп­ ределялись интегральные величины: изменения массы M x(t) при 2 ^ 0 и M 2(t) при 2>*0 относительно начального распределения; импульсы h (£) — в отри­ цательном и / а (0 — в положительном направлениях оси z; радиальная со­ ставляющая импульса / з(£) части среды, лежащей по одну сторону плоскости симметрии Ог; энергии £i(£) при z<0 и E 2(t) при z>0 (а также тепловая и ки­ нетическая ее компоненты по обе стороны границы раздела).

Приведем результаты и анализ расчетов Г. С. Романова и А. В. Тетерева,

Рис. 10.12. Временные зависимости интегральных вели­ чин при сосредоточенном энерговыделении у поверхности [ 199]. Е2— энергия выброшенного газа, / 2 — положитель­ ная составляющая осевого импульса, М2 — масса выбро­ шенного газа.

х= 5/3 касается поверхности z=0. Авторы отмечают, что поля течения раз­ личных вариантов становятся качественно подобными к моментам £~1, од­ нако различия ~10% в численных значениях параметров и интегральных величин, а также в характерных размерах областей, охваченных движением, сохраняются вплоть до очень больших времен £~100.

На рис. 10.12 приведены временные зависимости интегральных величин: энергии выброшенного газа Е 2 (рис. 10.12, а), положительной составляющей осевого импульса / а (рис. 10.12, б), массы выброшенного газа М 2 (рис. 10.12,в). Видно, что существует асимптотическая область £^10-М00, в которой M 2~ t a, I ^ t b (здесь также / 3~£с). Как известно [200], при выходе процесса

затель степени в законе распространения фронта ударной волны в среде x<j>~ а 0< 0,4. В полученных расчетах величина а 0, определенная по М (£),

I(t) при £^10, составила

0,36 < а 0< 0,42,

причем большие а 0 получены из Af2(£), минимальные — из / 3(£). Сопоставле­ ние результатов варианта 1а с другими данными *) для Л(£) указывает на их совпадение с графической точностью. При этом максимальное значение £р=1„сек, найденное в указанной работе, отвечает рассчитанному значению при £=1,6. По-видимому, можно утверждать, что в рассматриваемом интер­

При реализации метода расчетная область разбивалась на ячейки с раз­ мерами Дг=Дг=0,1 с максимальным числом их 66x34. На оси симметрии за­ давалось условие непротекания. При подходе возмущений от границ авто­ рами работы [199] производилось удвоение размера расчетной области — уве­ личение линейного масштаба ячейки в два раза. Полное число удвоений в не­ которых вариантах достигало 10. Изучение влияния укрупнения сетки пока­ зало, что и указанное разбиение расчетной области обеспечивает приемлемое пространственное разрешение для рассматриваемого круга задач. Например, при энерговыделении в однородной среде решение в области закономерностей точечного взрыва [190, 191] не отличалось более чем на 3% от автомодельного решения.

Приведенные данные вполне адекватно описывают особенности процесса передачи энергии в полупространство, заполненное холодным плотным газом. Указанный подход может широко использоваться и для численного модели­ рования двумерных нестационарных (быстропротекающих) газодинамических процессов [199 и др.]

Следует заметить, что метод крупных частиц сравнительно просто обоб­ щается и на случай магнитной гидродинамики. Так была исследована динамика распространения ударных волн в ближнем космическом пространстве от хромосферных вспышек на Солнце [203]. Сравнение результатов расчетов с экспе­ риментальными спутниковыми данными говорит о высокой точности вычисле­ ний и в случае нестационарной МГД-задачи.

6. Таким образом, с помощью метода крупных частиц удалось провести исследование весьма тонких нестационарных структур.

Следует заметить, что в этом случае, видимо, более предпочтительными являются полностью консервативные разностные схемы, введенные А. А. Са­ марским [63, 201, 202], в которых соблюдается детальный баланс как полной, так и внутренней энергии. Однако для рассматриваемого в настоящей работе класса задач (однокомпонентный химически нейтральный идеальный газ, те­ чения без учета сложной физики: электропроводности, магнитных полей и т. д.)

расчетные данные, полученные по консервативной (где соблюдается баланс для полной энергии) и по полностью консервативной схемам, весьма близки **). Зна­ чительно больший эффект оказывает здесь процедура уменьшения шага по времени А/. Кроме того, построение полностью консервативных алгоритмов вызывает много дополнительных трудностей (особенно для многомерных про-

*) См. сборник «Высокоскоростные ударные явления».— М.: Наука, 1973 (гл. 3, рис. 3.9).

**) Так, дисбаланс внутренней энергии при расчете ярко выраженной нестационарной задачи о дифракции ударной волны на угле не превышал нескольких процентов [185]. В за­ дачах «на установление» это расхождение еще меньше.

ГЛ. X] РАСЧЕТ НЕСТАЦИОНАРНЫХ ЗАДАЧ 267

странственных задач): они сложнее обычных консервативных схем, их введе­ ние изменяет характер устойчивости вычислительного процесса, увеличивает объем программы и требует большего времени счета. Поэтому в большинстве случаев экономически оправданным является проведение расчетов сложных многомерных задач газовой динамики по традиционным консервативным схе­ мам, описанным в гл. I I I .

Сделаем еще одно замечание. Г. С. Романов, М. В. Сузденков и А. В. Те­ терев при изучении данного класса задач провели полезные методические ис­ следования *). Они изучали влияние точности аппроксимации и различных видов искусственной вязкости в рамках метода крупных частиц на результаты численного решения двумерных газодинамических задач. Рассматривались схемы первого и второго порядков точности в рамках метода крупных частиц, а также схемы второго порядка точности метода «потоков» [55].

Схемы второго порядка точности метода крупных частиц были построены, аналогично [55], путем использования для расчетов на лагранжевом этапе перетоков массы через границы эйлеровых ячеек несимметричных формул. При этом для получения устойчивого решения необходимо на первом этапе вводить, как уже отмечалось в гл. III, явные члены с искусственной вяз­ костью q.

Проводя расчет одномерной задачи об отражающейся ударной волне и двумерных задач о сильном взрыве, авторы пришли к выводу, что схемы пер­ вого и второго порядка точности не дают принципиальных различий в резуль­ татах. Более того, для данного класса задач схемы первого порядка точности метода крупных частиц следует признать более предпочтительными, так как они позволяют более точно определить такие характеристики движения, как скорость распространения ударной волны, зависимость положения фронта разрыва от времени и т. д. Эти схемы обладают большей вычислительной ус­ тойчивостью, а также более экономичны при реализации.

*) Р о м а н о в Г. С., С у з д е н к о в М. В., Т е т е р е в А. В. О точности метода крупных частиц для задач газовой динамики.— Ж- вычисл. матем. и матем. физ., 1981, т. 21, № 3, с. 798—803.

Г Л А В А XI

ТЕЧЕНИЯ С ИЗЛУЧЕНИЕМ

Аэродинамика тел, движущихся в верхних слоях атмосферы, имеет ряд особенностей [370, 371, 389, 395, 408 и др.1. При полетах космических аппара­ тов в атмосфере Земли и других планет Солнечной системы заметно возрастает влияние эффектов радиационной газовой динамики, поэтому стали активно раз­ рабатываться различные методы решения задач с излучением. Это связано с тем, что при скоростях входа в атмосферу порядка второй космической тем­ пература за ударной волной может достигать 1000К и более. В этом случае перенос тепла излучением вносит существенный вклад в суммарный теплооб­ мен. При позышении скорости полета и, следовательно, температуры роль излучения в переносе тепла возрастает и в дальнейшем становится превали­ рующей.

В первой части настоящей главы (§§ 1—5) метод крупных частиц разви­ вается на случай расчета обтекания затупленных тел потоком излучающего газа [236, 290, 291, 332, 351]. Приводятся явные аппроксимационные формулы для давления и температуры воздуха в зависимости от плотности и удельной внутренней энергии, которые использовались для исследования течений дис­ социированного и ионизированного воздуха методом крупных частиц. Пред­ ставлены результаты систематических расчетов обтекания тел различной формы (цилиндра с торцевой передней частью, пластины конечной толщины, сферы) в широком диапазоне скоростей и высот.

Проблемы, связанные с исследованием взаимодействия излучения с вещест­ вом, являются в настоящее время весьма актуальными также и для задач физики плазмы. В то же время очевидна большая сложность указанного типа явлений — построение адекватных физических и математических моделей, численных алгоритмов расчета, определение механизма взаимодействия и т. д. Следует иметь в виду, что по своей природе данные явления носят ярко выра­ женный нестационарный характер и требуют рассмотрения, вообще говоря,

пространственной постановки.

Фундаментальные исследования для задач физики плазмы проводятся А. А. Самарским и его учениками [309—313 и др.]. Были детально изучены проблемы, связанные с различными механизмами взаимодействия излучения с веществом, гидродинамикой сжатия, рассматривались мишени и оболочки различной формы, исследовалась физика сжатия лайнеров, многослойные мишени с различным подбором оболочек и т. д. Проведены также аналитиче­ ские исследования ряда автомодельных задач, что позволило проследить за поведением важных физических параметров и получить рекомендации для повышения эффективности численного расчета.

Во второй части гл. XI (§§ 6—7) будет описан алгоритм метода крупных частиц для расчета некоторых задач взаимодействия лазерного излучения с веществом и приведены результаты расчетов.

Исследование задач радиационной газовой динамики проводилось сов­ местно с В. Н. Фоминым, расчеты производились В. П. Скотниковым (§ 5) и С. А. Кутасовым (§ 7).

Как известно [238, 239], система уравнений радиационной газовой ди­ намики является интегро-дифференциальной. В связи с этим для расчета об­ текания затупленных тел с учетом излучения было предложено и использо­ вано приближение объемного высвечивания [4, 240—242 и др.], которое приво­ дит к дифференциальным уравнениям, что существенно упрощает решение. В дальнейшем для учета излучения применялись /\-е и следующие приближе­ ния метода сферических гармоник и приближение плоского слоя. В частности, было показано, что нормальная составляющая вектора лучистого теплового потока намного больше касательной [243], а максимальное отличие Р3-го и Р6-го приближений от /Vго не превышает 2,5% для температур и 5% для лу­ чистых тепловых потоков [244]. Эти результаты подтверждают возможность учета излучения в приближении плоского слоя для расчета обтекания затуп­ ленных тел.

Для решения подобных задач использовались методы интегральных со­ отношений [26, 245], прямых [246], сеточно-характеристический [247, 248] и т. п., характерной особенностью которых является явное выделение ударной волны.

В настоящее время существует ряд эффективных подходов для расчета течения излучающего воздуха около затупленного тела [245—254 и др.].

При больших величинах лучистого теплового потока, поступающего к по­ верхности тела, происходит абляция теплозащитного покрытия. Исследова­ нию этих режимов и расчету течений с распределенным вдувом и испарением теплозащитного покрытия посвящены работы [246, 255—258 и др.].

Было показано, что при слабых вдувах собственное излучение паров теплозащитного покрытия превышает их экранирующее действие. Это приво­ дит к увеличению лучистого теплового потока, поступающего к поверхности обтекаемого тела. В режимах сильного вдува до поверхности тела доходит лишь длинноволновая составляющая лучистого теплового потока, поступаю­ щего к контактной поверхности (коротковолновая составляющая поглощается парами теплозащитного покрытия).

Более полный обзор современного состояния проблемы обтекания затуп­ ленных тел с учетом излучения можно найти в работах [259—263 и др.], где приводится подробная библиография по этому вопросу.

Следует отметить, что в настоящей главе рассматривается лишь одна из многих областей современной радиационной газовой динамики — задача гиперзвукового обтекания затупленного тела с учетом излучения. В настоя­ щее время ведутся работы и по другим направлениям и проблемам. Так, на­ пример, в ВЦ АН СССР Ю. Д. Шмыглевским и его сотрудниками проводятся исследования внутренних течений излучающего газа и рассматриваются ра­ циональные способы численного решения уравнения переноса лучистой энер­ гии [264—265 и др.].

§ 2. Аппроксимации термодинамических функций равновесного воздуха

1. При обтекании затупленных тел гиперзвуковым потоком воздуха урав­ нение состояния газа за фронтом ударной волны существенно отличается от уравнения состояния идеального газа с постоянным показателем адиабаты

K —CJ Cv—1,4. Так, например, при температурах выше 1000К начинает сказы­ ваться зависимость теплоемкости от температуры, а при Т>2000К происходит процессы диссоциации молекул и ионизации атомов, сильно влияющие на все термодинамические'функции [266, 337, 338 и др.].

данные.

Первые приближенные аналитические выражения термодинамических функций равновесного воздуха были предложены Хансеном [268] и В. В. Ми­ хайловым [269, 270] (согласно [271] их погрешность составляет 4,5 и 10% со­ ответственно). В работе С. С. Квашниной, В. П. Коробейникова [266] пред­ ставлены две формулы зависимости энтальпии h от скорости звука С и энт­ ропии 5 для температур от 1000 до 10 000К. Максимальная погрешность равна 5% для первой и 7% для второй формул. И. Н. Наумова [271] записывает уравнение состояния в виде В (р, Т)р=рТ. Представлены аналитические ап­ проксимации энтальпии h(p, Т) и В(р, Т) в диапазоне давлений / ? = 1 0 " 103 атм и температур Т<16 800К. При Т<2000К относительная ошибка аппрок­ симации не превышает 0,70%, а при 2000К<7,<;16800К для В и h — 0,47% и 0,84% соответственно.

Аналитическое представление энтальпии как функции давления и плот­ ности в интервале температур от 13 000 до 20 000К и давлений от 0,04 до 10 атм предложено В. К- Вертушкиным [272]. Относительная ошибка аппрок­ симации в основном не превышает ± 3% . В работе А. Н. Крайко [273] получены аналитические выражения для удельной энтальпии h и плотности р через дав­ ление и температуру в интервалах 10~3-4-103 атм и 400^-20 000К. Относи­ тельные отклонения от табличных данных не превышают ±3% для h и 1,5% для р.

Ю. Н. Дьяконов [274] аппроксимирует функцию х(р, q)=pC2lp (представ­ ляющую эффективный показатель изэнтропы) по q полиномами второй степени (по участкам), а по р — линейно по значениям lg /?. Параметр q является функ­ цией давления и плотности. Точность аппроксимации — порядка 1%. Даль­ нейшее развитие этих аппроксимаций содержится в работе Ю. А. Полежаева,

X.А. Рахматулина [275].

Вработе С. П. Синченко [276] проведена аппроксимация плотности и ско­ рости звука равновесного диссоциированного воздуха в зависимости от дав­

ления и энтальпии для диапазона давлений 10“ 3-н10 атм и температур от 300 до 11 000К (при давлении 10~3 атм) и до 20 000К (при давлении 10 атм), до­ пускающая экстраполирование с приемлемой точностью на меньшие давления (~10~4 атм) до температур ~5000К. Для давлений 10“ 3-М0 атм и для значе­ ний Н от —3,6 до 2,8 относительная ошибка аппроксимации функции Z = =Aphlp (где А =3,493) не превосходит 0,4%, а относительная ошибка аппрок­ симации скорости звука— 1,8% . Через Н здесь обозначена величина In h—

— In hN, где hN представляет собой энтальпию, соответствующую началу диссо­ циации азота.

В работе И. Н. Мурзинова [277] указывается, что уравнение состояния равновесного диссоциирующего воздуха можно приближенно представить в виде

Pl(ph) = fi (h) + ft{h)\gp,

где fu /2 — таблично задаваемые функции энтальпии h. При объеме каждой из этих таблиц в 30 чисел и использовании линейной интерполяции для вы­ числения функций в промежуточных точках максимальная погрешность апапроксимации уравнения состояния указанной зависимости не превосходит 2% в диапазоне 0,001 атм ^ р ^ 50 атм и h ^ 3 5 000 ккал/кг. В. И. Комаровский, Ю. П. Марасанов [278] описывают алгоритм, позволяющий рассчитывать для воздуха зависимость h=h(p, р) для диапазона изменения давления р= 10~3-^ •|* 102 атм и диапазона изменения плотности, соответствующего температурному