Численное моделирование некоторых нестационарных сверхзвуковых течений в каналах и струях тема автореферата и диссертации по механике, 01.02.05 ВАК РФ

Карпов, Алексей Валерьевич АВТОР
кандидата физико-математических наук УЧЕНАЯ СТЕПЕНЬ
Волгоград МЕСТО ЗАЩИТЫ
2005 ГОД ЗАЩИТЫ
   
01.02.05 КОД ВАК РФ
Диссертация по механике на тему «Численное моделирование некоторых нестационарных сверхзвуковых течений в каналах и струях»
 
Автореферат диссертации на тему "Численное моделирование некоторых нестационарных сверхзвуковых течений в каналах и струях"

ВОЛГОГРАДСКИЙ ГОСУДАРСТВЕНИНЫЙ УНИВЕРСИТЕТ

На правах рукописи

КАРПОВ АЛЕКСЕЙ ВАЛЕРЬЕВИЧ

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ НЕКОТОРЫХ НЕСТАЦИОНАРНЫХ СВЕРХЗВУКОВЫХ ТЕЧЕНИЙ В КАНАЛАХ И СТРУЯХ

(01.02.05. — Механика жидкости, газа и плазмы)

Автореферат диссертации на соискание учёной степени кандидата физико-математических наук

Волгоград 2005

Работа выполнена на кафедре вычислительных методов и программирования Волгоградского государственного университета.

Научный руководитель:

доктор физико-математических наук, доцент Васильев Е.И.

Официальные оппоненты:

доктор физико-математических наук,

старший научный сотрудник, доцент Мышенков Е.В.

доктор физико-математических наук, доцент Коваленко И.Г.

Ведущая организация:

Московский Государственный Университет им. М.В. Ломоносова, механико-математический факультет.

Защита состоится 22 апреля 2005 г. в 15 час. 20 мин. на заседании диссертационного совета К 212.029.05 при Волгоградском государственном университете по адресу: 400062, Волгоград, Университетский проспект, 100.

С диссертацией можно ознакомиться в библиотеке Волгоградского государственного университета.

Автореферат разослан «.../.» марта 2005 г.

49

«././.»I

Ученый секретарь диссертационного совета

Е. А. Мазепа

А. Общая характеристика работы

Актуальность темы. Исследование различных течений жидкости и газа с помощью моделирования на ЭВМ приобрело в последнее время широкое распространение. Это связано как с разработкой эффективных численных методов, так и появлением производительных ЭВМ. За последние три десятилетия, помимо разработки новых методов численного моделирования, были существенно улучшены старые.

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

Одним из важных применений численного моделирования является проверка полноты понимания физического явления путем сравнения результатов вычислений с результатами эксперимента. Если экспериментальные данные точны, то расхождение с результатами численного моделирования свидетельствует о неучтенных в данной модели эффектах или о существенных погрешностях численного метода, требующих устранения. Совершенствование математических моделей и исследование диапазона их применимости - актуальная задача и это направление активно может продвигаться лишь во взаимодействии с экспериментом. В ряде случаев численное моделирование способно заменить эксперимент, однако в большинстве случаев они дополняют друг друга. Для экспериментаторов численное моделирование позволяет лучше понять исследуемое явление, а для численного моделирования эксперимент позволяет выявлять диапазон применимости модели и изыскивать направления улучшения модели.

Актуальность вышеупомянутых задач совершенствования численных методов и математических моделей в последнее время возросла в связи с широким внедрением в практическое использование пакетов для инженер-но-прикладного моделирования, таких как ANSYS, LS-DYNA, STAR-CD. Одним из основных требований, предъявляемых этими программными- продуктами к численным методам и математическим моделям, является высокая надежность. Для численных методов это означает надежность вычисления, а для математических моделей - достоверно определенный диапазон их применения.

Цель работы. Целями данной работы являются:

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

- численное моделирование и исследование с помощью созданной методики пульсирующих течений перерасширенной струи из короткого осесим-метричного сопла;

- создание на основе сравнения с результатами эксперимента адекватной численной модели и моделирование течения, возникающего при распространении ударной волны в длинном канале с препятствием.

Научная новизна и практическая значимость.

Научная новизна результатов диссертации заключается в следующем:

- Разработаны: эффективный метод устранения флуктуации численного решения нефизического характера за фронтом сильной малоподвижной относительно расчетной сетки ударной волной; методика ускорения процедуры решения задачи Римана о распаде разрыва;

- в результате исследования задачи об истечении перерасширенной струи из короткого осесимметричного сопла обнаружены: новая разновидность ударно-волновой конфигурации в струе; диапазон двузначности структуры течения; эффект гистерезиса при вариации давления в пространстве;

- в диапазоне гистерезиса обнаружены пульсирующие режимы течения и впервые описан механизм неустойчивости, обусловленный зонами возвратного течения и трансзвуковым характером течения между ними.

Научная и практическая значимость. Предложенная автором методика вибрирующей сетки позволяет без больших затрат ресурсов ЭВМ и без сколь-нибудь значительного ущерба для производительности существенно повысить качество результатов при численном моделировании течений с сильными ударными волнами. Эта методика не затрудняет повышение порядка точности метода и выставления граничных условий. Единственное требование - необходимость реализации исходного метода на подвижной сетке. Предложенный подход к ускорению решения задачи Римана позволил повысить производительность методов, активно использующих решение этой вспомогательной задачи о распаде разрыва. Исследованные пульсирующие течения в перерасширенной струе могут возникать в струях за реальными короткими соплами с коническим насадком и имеют важное значение для оптимизации работы сопловых устройств реактивных двигателей на различных эксплуатационных режимах. Результаты последней главы могут быть использованы при проектировании коридоров гашения ударной волны в каналах и воздуховодах защитных сооружений.

Методика исследования. Исследования проводятся с помощью привлечения нескольких моделей для нестационарных течений сжимаемого газа, а именно: модели идеального газа Эйлера, модели Навье-Стокса вязкого теплопроводного газа и модели Навье-Стокса с учетом эффектов турбулентности. В качестве базового численного метода использовалась хорошо зарекомендовавшая себя W-модификация метода Годунова второго порядка точности, разработанная Васильевым Е.И. [ЖВМиМФ, 1996]. При анализе достоверности расчетов привлекаются результаты экспериментов, проводимых в рамках совместных исследований с группой ученых университета им. Бен-Гуриона (Израиль, Беер-Шева)

Результаты, выносимые на защиту. На защиту выносятся:

- Методика вибрирующей сетки, предназначенная для подавления численных нефизических флуктуации за фронтом сильных ударных волн;

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

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

Структура и объем диссертации. Диссертация содержит 130 страниц, включая 34 страницы с рисунками и 7 страниц списка литературы; состоит из введения, трёх глав, заключения и списка литературы. Работа содержит 54 рисунка. Список литературы состоит из 93-х наименований.

Апробация диссертации. Основные результаты диссертации докладывались на XVII международной школе по моделям механики сплошной среды (Казань, 2004), международной конференции «Геометрический анализ и его приложения» (Волгоград, 2004), международном семинаре «Супервычисления и мат. моделирование» (Саров, 2002), на V (2000) и VIII (2003) межвузовских конференциях молодых ученых Волгоградской области, конференциях профессорско-преподавательского состава ВолГУ (2001-2004), а также на научных семинарах математического факультета ВолГУ (рук. д.ф.-м.н. А.Г. Лосев и д.ф.-м.н. В.М. Миклюков, 2004) и (рук. д.ф.-м.н. Васильев Е.И., 2004).

Публикации. Основные результаты диссертации опубликованы в [1]-[8].

Личный вклад автора. В совместных работах [1-3, 5] автором проведена модификация и адаптация к решаемой задаче имеющегося программного кода, проведены расчеты и обработка результатов вычислений. В совместно опубликованных работах [6-8] автору принадлежат существенные модификации программного кода, использовавшегося соавторами для проведения вычислений. Постановка задач и анализ результатов численного решения осуществлялся совместно с научным руководителем.

В. Содержание работы

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

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

В § 1.1 приведены дифференциальные уравнения используемых математических моделей: 1) базовая для численного моделирования модель идеального совершенного газа, 2) модель вязкого теплопроводного газа и 3) вязкого теплопроводного газа с учетом эффектов турбулентности.

Уравнения последней модели (уравнения Навье-Стокса с уравнениями модели турбулентности) наиболее сложны из использованных, и по форме записи включают в себя остальные модели:

Здесь ? - тензор вязких напряжений, 5 - необратимая диссипация кинетической энергии* Т- поле температуры в газе. Для плоских двумерных течений дивергенция тензора вязких напряжений и диссипация кинетической энергии 5 выражаются через производные от компонент скорости в следующем виде:

Г6{циХ + 3(м{иу + Ух))у -2(м(их + гД М^Х + 3{м(иу + V,)) - 2{»{их + уу))

5 = | (2 (и, - ууГ + 2(и1 + ф + 3(и, +

Входящие в V- ? и уравнения модели ц = ця + цт, Я = Ят + Ят лекулярных и турбулентных вязкости и теплопроводности:

сумма мо-

И, =с„

Ркг

лт=-

« ' ' (г-1)Ргт

Главным элементом к-е модели турбулентности являются два дифференци-

альных уравнения для описания изменения величин к и е.

Выше приведенная полуэмпирическим модель замыкается заданием констант которые подбираться с учетом особенностей моделируемого течения.

В работе в основном использовались следующие значения: 0^=0.09,

При отбрасывании Цт и получим уравнения Навье-Стокса динамики вязкого теплопроводного газа. Бели же отбросить еще и правые, вязкостные, члены, то получим модель идеального газа (уравнения Эйлера) с постоянным отношением теплоемкостей (у= const).

В § 1.2 сформулирована разностная схема С.К. Годунова, использующая для определения числовых потоков решение задачи Римана о распаде разрыва на границах ячеек. Далее излагается используемая в работе W-модификация Васильева Б.И., идея которой заключается в добавлении двумерных поправок к аргументам задачи Римана. Используемый для вычисления повышающих порядок точности поправок переменный шаблон обеспечивает монотонность W-модификации. Далее приведены формулы поправок для 2 порядка точности в случае криволинейной подвижной сетки, используемой методикой улучшения результатов вычислений из § 1.4.

В § 1.3 предложена ускоряющая процедура приближенного решения задачи Римана, сокращающая общие затраты процессорного времени на ее решение. Суть предложенной процедуры заключается в начальном приближении с эмпирически полученной оценкой точности.

В § 1.4 рассматривается проблема побочных флуктуации численного решения за сильными разрывами. Предложены несколько сеточных методов подавления флуктуации. Идея метода ломаной границы ячейки (рис. \,а) заключается в преобразовании границы ячейки в ломаную границу при вычислении потоков. Метод расширенного шаблона (рис. 1,6) заключается в последовательном проведении расчетов на сетках, обозначенных сплошными линиями и пунктиром, и последующем осреднении полученных результатов. Эффективность методов нивелируется рядом существенных недостатков: удвоение времени вычислений; усложнение повышения порядка точности; усложнение при выставлении граничных условий; необходимость применения методов на всей расчетной области.

Этих недостатков лишена достаточно эффективная методика вибрации сетки (рис. 1,в), предлагаемая в § 1.4. Идея методики заключается в принудительной вибрации расчетной сетки относительно своего местоположения. Шаг метода состоит из четырех этапов (отмечены цифрами). На первом

этапе расчетная сетка сдвигается в каком-либо направлении, затем она возвращается в исходное положение. На третьем этапе расчетная сетка сдвигается в противоположном направлении, и после опять возвращается в исходное положение. Каждый этап представляет собой самостоятельный шаг по времени исходного численного метода. Особенностью метода является существенная разница в шагах по времени на соседних этапах - для определения шага на текущем этапе следует пользоваться шагом по времени на предпоследнем этапе. Подробнее методы изложены в работах [4, 5]. Результаты, полученные с их применением, были опубликованы в [1,2, 6 - 8].

В главе 2 проводится численное решение задачи об истечении перерасширенной струи из короткого осесимметричного сопла с прямолинейными образующими. Особенностью течения внутри такого сопла является наличие так называемого висячего скачка, возникающего из-за переразгона потока у стенки сопла. Вблизи стенки скачок слабый, но у оси симметрии он усиливается. В случае коротких сопел возможно взаимодействие этого скачка и формирующегося в струе фронта Маха. В данной главе в рамках трех моделей из § 1.1 проводится численное исследование закономерностей этого взаимодействия в зависимости от величины давления в пространстве рю в диапазоне 0.8-4.5 атм.

В § 2.1 приведена постановка задачи. Камера высокого давления (КВД) сопла представляет собой полубесконечный цилиндр, а расширяющаяся часть (камера низкого давления - КНД) имеет выход в неограниченное пространство. КВД и КНД разделены диафрагмой. Разность давлений покоящихся газов по обе ее стороны существенна. При разрыве диафрагмы в некоторый момент 10 начинается процесс запуска сопла: влево внутрь КВД распространяется волна разрежения, вправо по КНД и затем в пространство движется ударная волна.

При использовании модели Эйлера выставляются условия непротекания. В случае вязкого газа выставляется условие прилипания и условие отсутствия теплообмена со стенками. При использовании к-е модели турбулентности начальное значение кинетической энергии турбулентных пульсаций к принималось равным £ = 2-10^ а2, где ах - скорость звука в пространстве. Турбулентность внутри сопла считалась замороженной.

Длина сопла составляет 0.4 м, длина входной цилиндрической части равна 0.2 м. Радиус области пространства 0.5 м (рис. 2). Параметры газа в КВД следующие: То = 3200°А", ро = 1 кг/и1, — 9.32 атм. Параметры газа за соплом: варьировалось в диапазоне вычис-

ляется из уравнения состояния газа. Отношение теплоемкости газа

Численное моделирование изменения течения при изменении осуществлялось в рамках единого расчета. Характер изменения рх поясняется на рисунке: #»(/) уменьшалась с Л_ = 0.03 атм непрерывным образом с 1.54 атм. Уменьшение на Ар со значения р\ до нового значения р\- Ар происходило за время после чего до давление оставалось неизменным с

целью достижения стационарного течения. На участок гладкого перехода отводилось 3 мс, а на весь промежуток от ^ ДО отводилось 30 мс. В качестве единицы измерения времени взята 1 мс, для данного сопла примерно соответствующая единице безразмерного времени. Для повышения качества результатов численного моделирования в области пульсирующих ударных волн применяется методика вибрирующей сетки.

В § 2.2 приведены результаты численного расчета течения в трех различных постановках задачи и их анализ. В части 1 § 2.2 рассматривается численное решение в постановке идеального газа. На рис. 3 приведен фрагмент течения за соплом для р„ = 1.54 атм в виде изолиний числа Маха при М> 1 с наложенными поверх линиями тока. Главной особенностью течения внутри струи является необычность ударно-волновой конфигурации. Взаимодействие слабого скачка 8 с фронтом Маха приводит к его расщеплению, нижняя часть фронта И выдвинута вперед по оси относительно верхней части Б2. За фронтом И возникает протяженная приосевая дозвуковая область. Непосредственно за фронтом П в этой области образуется вихревая зона возвратного движения 71, на рис. 4 между точками А и В. Далее исследуется трансформация зоны Ъ, и всей струи в целом при изменении На рис. 5 представлена эволюция максимального числа Маха обратного течения в цирку-

ляционной зоне 71 как функция времени на 12 актов изменения р„Флуктуации М, сопровождающие акт изменения р», начиная с 1.48 атм не затухают. К моменту уменьшения р„, до 1.42 атм пульсации приобретают довольно регулярный периодический характер.

На рис. 6 приведено распределение числа Маха вдоль оси симметрии для р„ = 1.45 атм на протяжении одного периода пульсации. Отличительной

1 _Jill

s 1 АЛЛ

— и 4uuu нупг ............IT -'риршг "" 1 кпМШШ ш

48 1 /г^-45 рт~\лг ! ! р^ 1 39 д,-1.33 i Ршш 1 18

A BCD

Рис. 6. Изменение распределения числа Л/ вдоль оси симметрии.

X [мс]

Рис. 5. Зависимость от времени максимума М в вихревой зоне в процессе поэтапного уменьшения рас 1.54 атм до 1.18 атм.

чертой от стационарного режима является появление второй зоны возвратного течения Границы этой зоны отмечены буквами С и Б. Объяснить ее появление можно следующим образом: при прохождении через фронт ^ давление в газе увеличивается, однако полное давление из-

за необратимых потерь уменьшается (полное давление - давление торможения, то есть давление, получаемое в случае изэнтропического торможения частиц газа до нулевой скорости). Ударные волны, расположенные выше скачка 8, существенно слабее. В результате складывается ситуация, при которой давление газа в дозвуковой области за верхней частью фронта Маха выше, чем полное давление газа под этой областью, на оси симметрии струи, то есть даже полное торможение газа на оси струи не может привес -ти к выравниванию давления.

Такая ситуация в стационарном течении существовать не может. Возникает пульсирующий режим течения, при котором стадия торможения, накопления газа и возникновения области возвратного течения сменяется стадией выноса этой массы газа вниз по потоку. При уменьшении рх до 1.36 ат.м амплитуда колебаний увеличивается, равно как и срок жизни /2, достигающий половины периода колебания. Выяснено, что максимальный размер Ъ\ соответствует отсутствию /2, а максимальный размер Ъ2 достигается при близком к минимальному размеру /[ (см. рис. 7). Это свидетельствует о перекачке массы газа из в Ъ2 на этапе уменьшения 2\\. На этом этапе газ накапливается в /2. На этапе увеличения 2\\ происходит процесс сноса накопленной в /2 массы газа вниз по потоку. В итоге перенос массы в ядре пульсирующей струи осуществляется по своеобразной эстафете с образованием зон возвратного течения.

Рис. 7. Изолинии плотности я пульсирующей струе в моменты максимального (а) и минимального (б) размеров 2.\ при р» = 1.36 атм.

1 УшШ Ушйш

/ р„ = 1^24 ащ. /' ■" • 'чЛ 1.24 атм. \ г

1 р„ = 0.88 атм. ^

1

Рис. 8. Двузначность - стационарных численных решений в режиме изменения />„ = 1.33 -» 0.8 1.33 атм для невязкого течения

Уменьшение амплитуды колебаний в диапазоне рт 1.33+-1.27 атм может быть связано с прогрессирующим влиянием сверхзвуковой перемычки между вихревыми областями, которая периодически изолирует Ъ\ от /2, а также с периодическим исчезновением фронта Б2 с соответствующей перестройкой ударно-волновой конфигурации струи во время колебаний. Для ра = 1.15 и 1.12 атм происходит увеличение длительности существования сверхзвуковой области на оси симметрии, изолирующей Ъ\ от /2, что приводит к резкому уменьшению амплитуды колебаний, а при /Л«, = 1.09 атм формируется практически стационарное течение. При дальнейшем уменьшении пульсирующий режим отсутствует.

Было проведено сравнение стационарных конфигураций для различных значений рОбнаружено, что для вариантов с одним и тем же ро, но раз

Рис. 9. Сравнение картин эволюции течения в циркуляционной зоне с применением и без применения процедуры вибрирующей сетки для ,£>„,,= 1.18 атм.

личной предысторией получения, имеет место двузначность решения. Полученные конфигурации приведены на рис. 8.

Левая колонка соответствует режиму понижения рт. Критическим моментом является переход от р„ = 0.88 к 0.85 атм, при котором Zi исчезает и происходит переход к классической конфигурации ударных волн в перерасширенной струе с отсутствием взаимодействия между Fi и скачком S. Правая колонка кадров соответствует серии расчетов с повышением от 0.7 до 1.5 атм. В этом режиме фронт Fi приближается к точке отражения скачка S от оси, однако взаимодействие между ними достаточно слабое и качественных изменений в волновой конфигурации не происходит. С повышением рх с 1.3 атм до 1.33 атм фронт достигает точки пересечения скачка ,'» с осью симметрии, и возбуждаются колебания, близкие к периодическим.

Таким образом, для ударно-волновой конфигурации в струе при вариации рю имеет место эффект гистерезиса в диапазоне />,„6(0.85, 1.33) атм, причем на одной из ветвей петли гистерезиса в диапазоне атм особенности структуры ударно-волновой конфигурации порождают самоподдерживающиеся пульсации в струе.

Следует отметить, что, как правило, в пульсирующих течениях колебательные перемещения ударных волн около положения равновесия происходят со скоростями существенно меньшими, чем скорость набегающего потока. Это приводит к побочным численным флуктуациям. Для их устранения применялась методика вибрирующей сетки. Результаты применения приведены на рис. 9. Верхние кадры, полученные с применением методики, сравниваются с нижней последовательностью кадров, полученных на неподвижной сетке (вибрация сетки была отключена за несколько десятков шагов по времени до момента На нижнем правом кадре ^ фронт Fi перестал быть плоским, на нем появились заметные численные неоднородности из-за истончения фронта после отключения процедуры вибрирующей сетки. В совокупности с его небольшим наклоном это вызывает неоднородности решения в местах перехода через сеточные линии, вызванные флуктуация-ми скорости. Накопившиеся искажения фронта приводят к возникновению

а 1 . \ б

Рис. 10. Сравнение картин течения в стационарном струе, полученных в модели Навье-Стоксапри Не — 500 ООО (о) ивмоделигЭйлера(б)для р,ц = 1.5 атм.

вихрей вблизи оси симметрии (кадр

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

На кадре & приведены изолинии через длительный промежуток времени. Изменение формы косых скачков под воздействием вихрей спровоцировало переход регулярного взаимодействия в маховское. Вихревые возмущения на фронте Маха порождают неустойчивость Кельвина-Гельмгольца на тангенциальном разрыве. На оси наблюдается нефизический эффект струи, являющийся характерным признаком неразвитого «карбункул»-эффекта.

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

В части 2 § 2.2 рассматриваются результаты численного моделирования течения в модели вязкого теплопроводного газа. Анализ результатов предыдущей части позволяет предположить, что следует ожидать заметного влияния турбулентности в первой и второй циркуляционных областях, а также в пограничном слое на образующей в сверхзвуковой части сопла. Из-за ограниченной производительности доступных ЭВМ исследование влияния этих эффектов проводилось в двух постановках: модели вязкого газа с относительно низким числом Яе с граничными условиями прилипания на стенке сопла; модели вязкого газа с к-е моделью турбулентности без учета вязкости и турбулентности внутри сопла.

На рис. 10 приводится сравнение картин течения в первой постановке

Рис. 11. Эффект гистерезиса для стационарных конфигураций ударных волн в струе при изменении/>„ а диапазоне 1.5 -*■ 0.9 -* 1.5 ата. Число Ле = 500000.

при Яв - 500 000 с результатами, полученными в невязкой модели Эйлера, при величине />„ = 1.5 атм. Происходит отрыв потока от образующей сопла. Этот эффект размазывает скачок 8 и приводит к меньшей, чем в случае идеального газа, ширине струи. Это вызывает изменение ударно-волновой конфигурации внутри струи. В вязкой постановке при таком рю в дозвуковых зонах за фронтом не возникает обратного течения. С уменьшением точка отрыва потока от стенки сопла перемещается к выходу из сопла, и волновая конфигурация начинает постепенно приближаться к наблюдаемой в невязкой постановке задачи. Однако в вязком случае пульсирующие режимы течения во всем рассматриваемом диапазоне давлений не наблюдались. Присутствие зоны 2\ приводит к появлению эффекта гистерезиса решения по давлению аналогичного наблюдаемому в невязкой постановке задачи.

На рис. 11 изображена последовательность стационарных конфигураций при изменении рт. Левая последовательность кадров соответствуют понижению давления. Главной чертой этой ветви цикла является зона Ъг При уменьшении давления с 1 атм до 0.9 атм она исчезает. На другой ветви цикла при повышении давления восстановление зоны наблюдается прир*, порядка 1.4 атм. Таким образом, для ударно-волновой конфигурации в струе при вариации имеет место эффект гистерезиса в диапазоне

Роое(1, 1.33) атм. Однако в отличие от невязкого случая решение на обеих ветвях цикла является стационарным. В постановке задачи с моделью турбулентности учет вязкости и турбулентности проводился только в струе. Структура ударно-волновых конфигураций на входе в струю от случая идеальной постановки отличается незначительно. Однако высокая турбулентная вязкость в циркуляционной зоне за фронтом Б] уменьшает интенсивность возвратного течения и изменяет характер расположения дозвуковых

( 8000 12-, *-3»- 399 -э верхняя стенка N3 /

^—не—* N1 N2

^диафрагма ^препятствие

датчики давления

торец канала нижняя стенка

Рис. 12. Схематическое изображение канала и экспериментальной установки.

зон в струе. Приосевые дозвуковые зоны и дозвуковые зоны за верхними частями фронтов Маха (см. рис. 10) сливаются, между ними образуется локальная область сверхзвукового течения. Интенсивная диссипация на тангенциальных разрывах ведет к выравниванию параметров скорости поперек струи. Эти факторы препятствуют возникновению пульсаций в струе. Возвратного течения во второй дозвуковой зоне нет. Размеры первой зоны из-за сильного влияния вязкости значительно меньше. Турбулентная вязкость уменьшает интенсивность движения циркуляционных зон. Как следствие, во всех рассмотренных вариантах течение в обратной струе было дозвуковым, причем по мере уменьшения размеров зоны скорость обратной струи также уменьшается. Все это приводит к снижению степени сопротивляемости вихревой зоны набегающему потоку. В результате образование и исчезновение зоны /1 при циклической вариации происходило практически при одном и том же значении рх, то есть сколь-либо заметного эффекта гистерезиса и двузначности решения в рассматриваемом диапазоне обнаружить не удалось. Но это не исключает возможность гистерезиса при иной геометрии сопла.

В главе 3 проводится численное моделирование распространения ударной волны в длинном канале с одиночным препятствием - пластиной с вырезом. Препятствие используется для гашения интенсивности ударной волны. Этот прием применяется для защиты сооружений от воздействия ударной или детонационной волны (газовые трубопроводы, угольные шахты, шахты ракетных установок, воздуховоды в бомбоубежищах и др.). Экспериментальное исследование защитных свойств перфорированных экранов с различной геометрией отверстий проводилось в университете им. Бен-Гуриона (Израиль, Беер-Шева) группой под техническим руководством А.Б. Британа. Результаты численных исследований представлены в данной главе.

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

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

В § 3.1 приведена постановка задачи. Рассматривается течение в канале с препятствием в виде пластины. Правый торец канала закрыт (см. рис. 12). Газ в левой части канала (КВД) отделен от покоящегося газа справа диафрагмой, расположенной до препятствия. Препятствие закрывает 60% сечения канала. После разрыва диафрагмы в начальный момент *о формируется ударная волна с числом Маха -Мь которая движется вправо. При достижении препятствия часть ударной волны отражается от него и распространяется против потока, другая часть проходит через отверстие и распространяется по каналу вправо. За препятствием формируется струя. Затем прошедшая через препятствие ударная волна достигает закрытого торца канала, отражается от него и начинает распространяться обратно, против потока. При достижении отраженной ударной волной препятствия интенсивность втекающей струи резко снижается и она превращается в слабую дозвуковую струю (эффект заглушки струи).

Для моделирования течения использовались две модели: модель Эйлера идеального газа, а также модель вязкого теплопроводного газа с моделью турбулентности. При использовании модели Эйлера на стенках канала и на препятствии выставляются условия непротекания. При использовании к-е модели турбулентности начальное значение кинетической энергии турбулентных пульсаций к — 10"4 ао2, где ао - скорость звука в канале в начальный момент. Начальное значение скорости диссипации кинетической энергии е задавалось таким образом, чтобы цг~ 2.25 ц„. При численном решении использовалась прямоугольная равномерная сетка различного разрешения.

В модели с к-е турбулентностью условие непротекания заменялось условием прилипания. Также изменялись и граничные условия теплообмена со стенками канала для улучшения согласия с экспериментальными измерениями. Выбирались условие отсутствия теплообмена газа со стенками или условие интенсивного теплообмена между газом и стенкой, в результате которого температура газа достигает температуры стенки 298°К. В качестве начальных условий использовались следующие значения параметров: Г=298°К,/?о= 1 атм, ро вычислялась из уравнения состояния газа. Отношение теплоемкости газа

При проведении экспериментов отклонение в числе Маха набегающей ударной волны не превышало 1%, погрешность в измерении ее скорости была меньше 0.7%, а погрешность при измерении давления была в пределах ± 7%. На основании анализа экспериментальных данных был сделан вывод: за изначально нестационарным всплеском следует квазистационарное течение. Такое поведение давления типично для течений, формируемым глав

Ыу= 32 Ыу = 320

Рис. 13. Увеличенные фрагменты численного решения в идеальной постановке у препятствия для разных размеров расчетной сетки на начальном этапе развития течения. Вертикальный размер сетки Ыу соответствует 1 или 10 ячейкам на 1 мм.

Рис. 14. Результаты численного решения в вязкой постановке с к-е моделью в начальной (а) и промежуточной (б) стадиях развития течения.

ным образом струей. Течение в этой области неоднородно, неизэнтропично, и дает максимальный прирост турбулентной кинетической энергии из-за струи и волн, генерируемых препятствием.

В §3.2 содержатся результаты численного моделирования в двух различных постановках с вариациями граничных условий, их анализ и сравнение с результатами эксперимента. В части 1 §3.2 содержатся результаты численного моделирования начального этапа течения в идеальной постановке для различного количества ячеек расчетной сетки (на рис. 13 приведены результаты для минимальной и максимальной использовавшихся расчетных сеток). С увеличением разрешения сетки увеличивается нестабильность струи и появляются осцилляции в параметрах во всей области течения. Проведенный соавторами работы [3] спектральный анализ флуктуации показал увеличение их частоты и расширение диапазона частот при уменьшении размера ячейки, что подтверждает турбулентную природу течения. Увеличение разрешения расчетной сетки заставляет семейство вихрей на границе струи подстроиться под новую расчетную сетку, что приводит к изменению характера взаимодействия вихрей друг с другом и с вторичными отраженными ударными волнами. Это свидетельствует об ограниченности области применения модели Эйлера при моделировании этой задачи и необходимость учета эффектов турбулентности.

Вторая часть §3.2 содержит результаты численного моделирования течения с применением уравнений модели. Подробно проанализирован процесс развития течения вязкого газа на протяжении 3280 мкс с момента прихода ударной волны к препятствию. Результаты численного решения на расчетной сетке 2700x96 ячеек приведены на рис. 14. Из анализа результатов следует, что основным элементом формирования течения является струя. Особенностью струи является ее отрыв от нижней стенки и после-

а 6 а б

1 2 Рис. 15. Сравнение диаграмм давления, полученных в ходе эксперимента (а) и по результатам численного моделирования (б) в постановке задачи с вязким теплопроводным газом с учетом эффектов турбулентности. Кадры 1 соответствуют показаниям датчиков на нижней стенке канала, кадры 2 - на верхней стенке. Условия прилипания.

дующее отражение от верхней. Об этом свидетельствует искажение фронта отраженной от торца канала ударной волны (см. рис. 14,

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

Сравнение расчетных эпюр давления с экспериментальными данными в трех сечениях канала не позволяет однозначно отдать предпочтение конкретному виду граничных условий. При условии прилипания с экспериментом лучше согласуются датчики на нижней стенке канала, а при условии скольжения - на верхней стенке. На взгляд автора, более правдоподобны условия прилипания, но более определенного ответа нужна экспериментальная фотосъемка картины течения.

В части 3 § 3.2 содержатся результаты численного моделирования в невязкой постановке. Разность скоростей по обе стороны тангенциального

ж

Рис. 17. Фрагмент численного решения в идеальной постановке. Изолинии плотности. а б

■1 №

--#

к.

и

— м \

-41 — Ъ М

М-

Рис. 18. Диаграммы давления, построенные по результатам численного моделирования в идеальной постановке, (а) - нижняя стенка канала, (б) - верхняя стенка.

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

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

В заключении диссертационной работы перечислены наиболее важные результаты и выводы работы.

Публикации автора

(в обратном хронологическом порядке)

1. Васильев ЕЖ., Карпов А.В. Численное моделирование нестационарных сверхзвуковых течений в струях // Труда мат. центра им. Лобачевского. Т. 27. Модели мех. сплош. сред. Материалы XVII сессии межд. школы по моделям механики сплошной среды. Казань. 4-10 июля 2004 г.: изд-во. Казан, мат. общ-ва, 2004, С. 60-64.

2. Карпов А. В., Васильев ЕЖ. Численное моделирование пульсирующих течений в перерасширенной сверхзвуковой струе // Тезисы докладов межд. школы-конференции "Геометрический анализ и его приложения". 24-30 мая 2004 г., Волгоград. С. 69-70.

3. A.Britan, A. V.Karpov, E.I. Vasiliev, O.Igra, G.Ben-Dor, E.Shapiro. Experimental and Numerical Study of Shock Wave Interaction with Perforated Plates // Journal of Fluids Engineer, 2004, v. 126, p. 399 - 409

4. Карпов А. В. Исследование эффективности применения метода вибрирующей сетки в ID-задачах газовой динамики. // Вестник ВолГУ; С. 9. 2003-2004. Ч. 2. С. 20-23

5. Карпов А.В., Васильев ЕЖ. Метод искусственного излома сеточных линий для погашения численных флуктуации за фронтом сильных ударных волн в схемах сквозного счета.//Вестник ВолГУ,С. 1.2002. В. 7, с. 40—49.

6. Зенович А.В., Карпов А.В., Васильев ЕЖ. Автоколебательный режим стационарного сверхзвукового обтекания двойного клина // Вестник ВолГУ, С. 1. 2002. В. 7, С. 34—39.

7. Карпов А.В, Зенович А.В., Васильев ЕЖ. Автоколебательный режим ударно-волновой конфигурации при обтекании двойного клина стационарным сверхзвуковым потоком // Труды I Школы-семинара «Мат. моделирование, вычислительная механика и геофизика», 14-18 октября 2002 года. - Ростов-на-Дону: Изд. РГУ, 2002 - С. 117-119.

8. Васильев ЕЖ., Зенович А.В., Карпов А.В. Автоколебательные режимы сверхзвукового обтекания двойного клина // Тезисы Межд. семинара Супервычисления и мат. моделирование. 17-21 июня 2002. Саров. С. 33

Подписано в печать 14.03.2005 г. Формат 60x84/16. Бумага офсетная. Гарнитура Тайме. Усл. печ. л. 1,2. Тираж 100 экз. Заказ 54.

Издательство Волгоградского государственного университета. 400062, Волгоград, просп. Университетский, 100.

0101 - CL 0.6

 
Содержание диссертации автор исследовательской работы: кандидата физико-математических наук, Карпов, Алексей Валерьевич

Введение

1. Математические модели и методика численного моделирования.

1.1. Математические модели движения сжимаемого газа.

1.2. Повышение порядка точности схемы Годунова.

1.3. Методика ускорения решения вспомогательной задачи 24 Римана.

1.4. Сеточные методы погашения нефизических эффектов 29 численного решения за сильными малоподвижными ударными волнами.

1.4.1. Метод искусственного излома сеточных линий.

1.4.2. Метод вибрирующей сетки.

2. Численное моделирование истечения перерасширенной струи 49 газа из короткого осесимметричного сопла

2.1. Постановка задачи и методика численного решения

2.2. Результаты расчетов.

2.2.1. Результаты расчетов в модели идеального газа.

2.2.2. Результаты расчетов в модели вязкого теплопроводного 84 газа.

3. Численное моделирование распространения ударной волны в 93 канале с препятствием.

3.1. Постановка задачи и методика решения.

3.2. Результаты расчетов.

3.2.1. Начальный этап формирования течения без учета вязких 100 эффектов.

3.2.2. Результаты расчетов в модели Навье-Стокса с учетом 103 эффектов турбулентности.

3.2.3. Результаты расчетов в модели идеального газа. 115 Заключение 122 Список литературы

 
Введение диссертация по механике, на тему "Численное моделирование некоторых нестационарных сверхзвуковых течений в каналах и струях"

Исследование различных газодинамических течений с помощью моделирования на ЭВМ приобрело в последнее время широкое распространение. Это связано с появлением как простых и надежных численных методов, так и доступных производительных ЭВМ. За последние три десятилетия разработаны новые методы численного моделирования и существенно улучшены старые. Однако, несмотря на бурный рост вычислительной мощности ЭВМ и прогресс в развитии методов компьютерного моделирования газодинамических течений, задача совершенствования методов на настоящий момент остается актуальной. Так, для большинства методов второго порядка сужение зон локализации ударных волн в численном решении порождает побочные флуктуации решения за сильными ударными волнами. Этой проблеме за последние десятилетия посвящено большое количество публикаций [36, 39, 45, 53, 56, 57, 74, 75,81,87, 88, 92].

Одним из важных применений численного моделирования является проверка полноты понимания физического явления путем сравнения результатов вычислений с результатами эксперимента. Если экспериментальные данные точны, то расхождение с результатами численного моделирования свидетельствует либо о неучтенных в данной модели эффектах, либо о погрешностях в методике численного расчета. Совершенствование математических моделей и исследование диапазона их применимости -актуальная задача и это направление активно может продвигаться лишь во взаимодействии с экспериментом. В ряде случаев численное моделирование способно заменить эксперимент, однако в большинстве случаев они дополняют друг друга. Для экспериментаторов численное моделирование позволяет лучше понять исследуемое явление, а для численного моделирования эксперимент позволяет выявлять диапазон применимости модели и изыскивать направления улучшения модели.

Актуальность вышеупомянутых задач совершенствования численных методов и математических моделей в последнее время возросла в связи с широким внедрением в практическое использование пакетов для инженерно-прикладного моделирования, таких как ANSYS, LS-DYNA, STAR-CD. Одним из основных требований, предъявляемых этими программными продуктами к численным методам и математическим моделям, является высокая надежность. Для численных методов это означает надежность вычисления, а для математических моделей - надежность диапазона их применения.

Целями данной работы являются:

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

- численное моделирование и исследование с помощью созданной методики пульсирующих течений перерасширенной струи из короткого осесимметричного сопла;

- создание на основе сравнения с результатами эксперимента адекватной численной модели и моделирование течения, возникающего при распространении ударной волны в длинном канале с препятствием.

Исследования проводятся с помощью привлечения нескольких моделей для нестационарных течений сжимаемого газа, а именно: модель идеального газа Эйлера, модель Навье-Стокса вязкого теплопроводного газа и модель Навье-Стокса с учетом эффектов турбулентности. В качестве базового численного метода использовалась хорошо зарекомендовавшая себя W-модификация метода Годунова второго порядка точности, разработанная Васильевым Е.И. [6]. При анализе достоверности расчетов привлекаются результаты экспериментов, проводимых в рамках совместных исследований с группой ученых университета им. Бен-Гуриона (Израиль, Беер-Шева)

Структура и объем работы. Диссертация состоит из введения, трех глав, заключения и списка литературы; содержит 130 стр., включая 34 стр. с рисунками и 7 стр. списка литературы. В работе 93 библиографических ссылки.

 
Заключение диссертации по теме "Механика жидкости, газа и плазмы"

Основные результаты и выводы работы:

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

2. Разработана методика, уменьшающая ресурсоемкость массовой процедуры решения задачи Римана о распаде разрыва с контролем точности и увеличивающая производительность методов, использующих результаты решения этой вспомогательной задачи;

3. Проведено численное решение задачи об истечении перерасширенной струи газа из короткого осесимметричного сопла при перепаде давления в пространстве 0.8- 1.5 атм в постановке идеального газа, вязкого теплопроводного газа и вязкого теплопроводного газа с к-е турбулентностью, в ходе которого:

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

3.2. Обнаружен диапазон двузначности течения и эффект гистерезиса для ударно-волновой конфигурации в струе при вариации внешнего давления в моделях Эйлера и Навье-Стокса;

3.3. В модели Эйлера в области гистерезиса обнаружены режимы пульсирующих течений перерасширенной струи;

3.4. Показано, что пульсирующие режимы струи в модели Эйлера связаны с образованием зон возвратного течения в струе и с трансзвуковым характером потока между зонами;

3.5. Продемонстрирована высокая эффективность метода вибрирующей сетки при моделировании пульсирующих течений с ударными волнами;

3.6. Установлено, что эффекты турбулентности сильно сужают область гистерезиса и стабилизируют течение;

4. Проведено численное решение задачи о распространении ударной волны в длинном канале с препятствием в постановках идеального газа и вязкого газа с учетом эффектов турбулентности, в ходе которого:

4.1. Установлена ограниченность применимости модели Эйлера в данной задаче на длительных временах и показана необходимость учета эффектов турбулентности в головной части нестационарной струи;

4.2. На основании серии расчетов установлен наиболее приемлемый вариант к-е модели турбулентности, который обеспечивает хорошее соответствие результатов вычислений с результатами экспериментов;

4.3. Обнаружен зигзагообразный характер фарватера струи за препятствием, обусловленный отрывами струи от стенок канала и приводящий к сильной неравномерности и неоднородности движения фронта отраженной ударной волны.

Заключение.

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

 
Список источников диссертации и автореферата по механике, кандидата физико-математических наук, Карпов, Алексей Валерьевич, Волгоград

1. Баженова Т.В., Гвоздева Л.Г. Нестационарные взаимодействия ударных волн. М.: Наука. 1977. 271с.

2. Британ А.Б., Васильев Е.И. Особенности формирования течения в профилированном сопле ударной трубы // Доклады АН СССР. 1985. т. 281. № 2. с. 295-299.

3. Британ А.Б., Васильев Е.И. Исследование запуска профилированного сопла ударной трубы большого диаметра // Известия АН СССР. МЖГ. 1986. №5. с. 88-95.

4. Васильев Е.И. Нестационарное истечение струи в затопленное пространство // Известия АН СССР. МЖГ. 1984. № 1. с. 42-46.

5. Васильев Е.И. W-модификация метода С.К. Годунова и ее применение для двумерных нестационарных течений запыленного газа // ЖВМ и МФ. 1996. т. 36. № I.e. 122-135.

6. Васильев Е.И. Численное моделирование нестационарных сопловых и струйных течений запыленного газа // Вестник ВолГУ. Серия 1, математика и физика. 1996. № I.e. 55-64.

7. Васильев Е.И., Данильчук Е.В. Численное решение о развитии течения в ударной трубе при поперечном выдвижении диафрагмы // Известия АН СССР. МЖГ. 1994. № 2. С. 147-154.

8. Васильев Е.И. W-модификация метода Годунова и ее приложения в моделировании газодинамических течений с ударными волнами. Диссертация на соискание ученой степени д.ф.-м.н. Волгоград, ВолГУ. 1999.213с.

9. Ганжело А.Н., Крайко А.Н., Макаров В.Е., Тилляева Н.И. О повышении точности решения газодинамических задач // Современные проблемы аэромеханики. М.: Наука. 1985. с. 87-102.

10. Годунов C.K. Разностный метод численного расчета разрывных решенийуравнений гидродинамики // Математический сборник. 1959. т. 47.с. 271-306.

11. Годунов С.К, Забродин A.B., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. М.: Наука, 1976.400с.

12. Гринь В.Т., Крайко А.Н., Славянов H.H. Решение задачи о запуске сопла, вмонтированного в торец ударной трубы // Изв. АН СССР. МЖГ.щ 1981. №6. С. 117-123.

13. Дедеш В.В. Об одном методе построения схем типа Годунова // Доклады АН СССР. 1991. т. 321. № 1. с. 36-39

14. Зенович A.B., Карпов A.B., Васильев Е.И. Автоколебательный режим стационарного сверхзвукового обтекания двойного клина // Вестник ВолГУ, С. 1.2002. В. 7, с. 34—39.

15. Карпов A.B., Васильев Е.И. Метод искусственного излома сеточных линий для погашения численных флуктуаций за фронтом сильных ударных волн в схемах сквозного счета. // Вестник ВолГУ, Серия 1. 2002. В. 7, с. 40-49.

16. Колган В.П. Применение операторов сглаживания в разностных схемах высокого порядка точности // ЖВМ и МФ. 1978. т. 18. № 5. с. 1340-1345.24