На правах рукописи
япунов Сергей Владимирович
Разработка и Исследование ЧИСЛЕННЫХ СХЕМ ВЫСОКОГО ПОРЯДКА ТОЧНОСТИ ДЛЯ РЕШЕНИЯ УРАвнений газовой динамики НА неструктурированных сетках
05.13.18. - Математическое моделирование,
численные методы и комплексы программ
Автореферат
Диссертации на соискание ученой степени доктора
физико-математических наук
Москва - 2008
Работа выполнена в Федеральном государственном унитарном предприятии Центральном Аэрогидродинамическом институте им. проф. Н.Е. Жуковского (ФГУП ЦАГИ)
Официальные оппоненты: доктор физико-математических наук,
профессор Крайко Александр Николаевич,
доктор физико-математических наук,
профессор Хлопков Юрий Иванович
доктор технических наук,
доцент Босняков Сергей Михайлович
Ведущая организация - Вычислительный центр РАН (ВЦ РАН, Москва)
Защита диссертации состоится 26 июня в 15.00 на заседании диссертационного совета Д 215.001.01 в Военно-воздушной инженерной академии имени профессора Н.Е.Жуковского по адресу 125190, г.Москва, ул.Планетная, д.3.
С диссертацией можно ознакомиться в библиотеке Военно-воздушной инженерной академии имени профессора Н.Е.Жуковского
Автореферат разослан 2008 г.
Ученый секретарь
Диссертационного совета Д 215.001.01
Кандидат физико-математических наук
А.С.Ненашев
ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ
Актуальность проблемы. В настоящее время все больший объем информации о течениях жидкостей и газов и силах, действующих на движущиеся в них тела и ограничивающие поверхности, получается с использованием методов вычислительной аэродинамики - численных методов решения систем уравнений Эйлера и Навье-Стокса. Применение этих методов позволяет сократить объем промышленных экспериментальных исследований, лучше понять физические особенности течений и, в ряде случаев, позволяет получить информацию, которую крайне сложно, а порой и невозможно, получить в эксперименте. Главенствующее место среди численных методов занимают сеточные методы, основанные на дискретной аппроксимации уравнений газовой динамики.
Основным требованием, предъявляемым к таким сеточным методам, является, прежде всего, обеспечение высокой точности (малой численной ошибки) получаемых результатов при минимально необходимых ресурсах ЭВМ (времени и объеме памяти). Кроме того, желательно максимально автоматизировать процесс генерации вычислительной сетки, обеспечить возможность генерации сетки вокруг объектов сложной геометрии, обеспечить точное описание (лразрешить) особенности течений (скачки уплотнения, пограничные слои, отрывные зоны и т.п.), устойчивую сходимость к решению для максимально возможного числа случаев обтекания (робастность).
Перспективными подходами к построению численных методов, удовлетворяющих перечисленным требованиям, являются применение неструктурированных вычислительных сеток, численных схем высокого порядка точности, адаптация сетки и численной схемы к решению. В настоящей диссертационной работе предложены и исследованы методы адаптации неструктурированных сеток и локальной адаптации порядка точности численного метода к решению на базе метода конечного элемента Галеркина с разрывными функциями. Впервые этот метод был предложен в работе (Reed, Hill 1973) для решения уравнения переноса нейтронов. Дальнейшие многочисленные исследования были посвящены анализу математических аспектов метода, таких как скорость сходимости и пр. (LeSaint, Raviart, 1974, Johnson, Pitkaranta, 1986), а также развитию метода (Cockburn, Shu, 1989, Cockburn, Lin, Shu, 1989, Cockburn, Hou, Shu, 1990, Cockburn, Shu, 1998), (Bassi, Rebay 1997), (Warburton, Lomtev, Kirby, Karniadakis 1998, Lomtev, Karniadakis 1997).
Преимущества метода Галеркина с разрывными функциями заключаются в следующем.
- Данный метод легко адаптируется к неструктурированным сеткам и, следовательно, удобен для работы с областями сложной геометрии.
- Порядок точности метода зависит от максимальной степени полиномов базисных функций, использующихся для аппроксимации численного решения. Метод может быть сформулирован формальным образом для произвольного порядка точности на гладких решениях путем расширения подпространства базисных функций и увеличения максимального порядка полиномов.
- Метод обладает большой гибкостью, поскольку порядок базисных функций может меняться от элемента к элементу, что важно с точки зрения адаптации метода к решению.
Таким образом, актуальность работы определяется потребностью создания высокоэффективных численных методов решения уравнений газовой динамики, позволяющих получать решение задач обтекания конфигураций сложной геометрии с высокой точностью при минимальных затратах памяти и времени работы ЭВМ.
Практическая значимость работы состоит в разработке принципов адаптации неструктурированных сеток и порядка точности численной схемы к решению. Проведены методические исследования, включающие решения модельных задач, уравнений Эйлера, Навье-Стокса и Рейнольдса, демонстрирующие эффективность данных подходов с точки зрения повышения точности численных решений. Разработана научно-методическая основа для реализации предложенных подходов в промышленных программах решения уравнений газовой динамики.
Цель диссертационной работы состоит в теоретической и методической разработке методов ускорения расчета и повышения точности результатов численного решения уравнений газовой динамики - уравнений Эйлера (невязкий случай) и уравнений Навье-Стокса и Рейнольдса (вязкий случай) путем адаптации не только неструктурированной расчетной сетки, но и локального порядка точности численной схемы. Данные подходы основаны на модификации метода Галеркина с разрывными функциями (DG - Discontinuous Galerkin в англоязычной литературе), который является комбинацией метода конечного элемента и метода конечного объема типа метода Годунова. Особое внимание уделяется анализу порядка точности получаемых схем, выявлению их преимуществ по сравнению со стандартными схемами типа Годунова. Подробно рассмотрены возможности, которые обеспечивает схема DG с точки зрения адаптации к решению на неструктурированных сетках. Приведены результаты исследований на примерах одномерных и двумерных задач. Выбор этих задач в значительной степени обусловлен интересом к анализу порядка ошибки численного решения, что требует знания точного решения.
Общая методика выполнения исследований состоит в
- аналитической формулировке численной схемы решения различных законов сохранения, включая уравнения Эйлера, Навье-Стокса и Рейнольдса на базе варианта метода конечного элемента - метода Галеркина с разрывными функциями, обеспечивающего произвольный порядок точности численной схемы,
- разработке вычислительных программ, реализующих этот численный метод для различных задач с использованием неструктурированных разностных сеток и процедур адаптации сеток к решению,
- анализе порядков точности метода на примерах тестовых задач и решений уравнений газовой динамики,
- выявлении положительных особенностей метода с точки зрения адаптации схемы к решению путем адаптации вычислительной сетки (h-refinement) или порядка точности схемы (p-refinement).
Научная новизна работы заключается в разработке и исследовании новых подходов к созданию эффективных численных схем решения уравнений газовой динамики, в том числе на неструктурированных сетках, обеспечивающих высокую точность решения при умеренных затратах ресурсов ЭВМ и широкие возможности адаптации к решению, не только путем адаптации сетки, но и путем локального изменения порядка точности численной схемы.
Автор защищает следующие результаты:
- Принципы построения численных схем решения законов сохранения, в частности, уравнений газовой динамики (уравнений Эйлера, Навье-Стокса и Рейнольдса), обеспечивающих уменьшение времени расчета и повышение точности результатов, по сравнению с классическими методами типа метода Годунова, путем различных способов адаптации неструктурированной сетки и локального порядка точности разностной схемы к решению на базе численной схемы конечного элемента - метода Галеркина с разрывными функциями.
- Результаты исследований предлагаемых подходов к адаптации, представляющих научно-методический базис для реализации разработанных подходов в промышленных программах решения уравнений газовой динамики.
Практическая ценность и реализация работы. Предлагаемые в диссертации способы повышения эффективности численных методов решения уравнений газовой динамики апробированы на модельных задачах и доведены до стадии начала промышленной реализации. Ряд подходов реализован в виде вычислительных программ, которые используются для проведения расчетных исследований обтекания различных элементов самолетов (крыловые профили, взлетно-посадочная механизация) при выполнении НИР ЦАГИ по контрактам с Роспромом и при проведении инициативных исследований.
Апробация работы. Методы прошли тщательное тестирование путем сравнения результатов расчетов с имеющимися точными решениями и результатами других авторов. Результаты проведенных исследований позволили сделать ряд выводов относительно возможности реализации высокого порядка точности на гладких решениях, особенностей применения рассматриваемого подхода для тел с криволинейной границей, возможностей адаптации схемы к решению путем модификации расчетной сетки, локального изменения порядка схемы или модификации системы базисных функций. Автор полагает, что рассмотренные подходы могут послужить основой для разработки эффективных численных методов и промышленных программ расчета вязких трехмерных течений около тел сложной геометрии.
Основные результаты проведенного автором исследования содержатся в 31 публикациях, опубликованных в российских научных изданиях и за рубежом, а также докладывались на международных и российских научно-технических конференциях, в том числе, на 6-ой международной конференции по генерации сеток для вычислительной аэродинамики (1998), 16-м Конгрессе Международной ассоциации математического и компьютерного моделирования IMACS (2000), 3-ей Европейской конференции по механике жидкости EUROMECH (1997), 5-ой Российско-Китайской конференции по аэродинамике и механике полета (1997), Международных Конгрессах по авиационным наукам ICAS (1990, 1992, 1996), школах-семинарах Аэродинамика летательных аппаратов (1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007), Международной конференции молодых ученых и специалистов Современные проблемы аэрокосмической науки и техники (2000), франко-российском семинаре ONERA-ЦАГИ (2002) и др.
В данную диссертацию включены исследования, поддержанные РФФИ (Проекты № 98-01-00032-а, 1998, №00-01-00070-а, 2000, №02-01-00124-а, 2002, №03-01-00236-а, 2003, №06-01-00283-а, 2006).
Объем работы. Диссертация состоит из введения, пяти глав, заключения и списка литературы, включающего 93 наименования. Общий объем - 127 страниц.
Краткое содержание работы
Во введении изложены общие подходы, используемые в вычислительной аэродинамике, требования к вычислительным сеткам и принципам построения схем дискретизации законов сохранения. Приведен обзор преимуществ и недостатков неструктурированных сеток, схем высокого порядка точности, отмечена необходимость адаптации сетки и схемы к решению. Дан анализ подходов к построению численных схем на основе метода конечного объема и метода конечного элемента. Кратко описано место метода Галеркина с разрывными функциями среди других численных схем, как способа получения схем произвольного порядка точности, в том числе на неструктурированных сетках. Сформулирована цель диссертационной работы и приведена краткая аннотация ее глав.
В Главе I рассмотрена общая теория метода Галеркина с разрывными функциями (DG) для численного решения одномерных законов сохранения. Метод представляет собой симбиоз метода конечного элемента и конечного объема. Метод применен к решению модельных задач для одномерных уравнений конвекции, теплопроводности и уравнения Бюргерса с вязкостью. Данные примеры являются простейшими моделями конвективных, вязких членов уравнений газовой динамики, разрывных решений. Основное внимание уделено анализу порядка ошибки численного решения.
В з 1.1 рассмотрена общая теория метода Галеркина с разрывными функциями применительно к решению одномерных законов сохранения. Конкретно рассматривается следующая задача определения функции u с начальным условием и периодическими граничными условиями для одномерного уравнения конвекции-диффузии с произвольной функцией невязкого потока f(u)
в области Q= (0,1) × (0,T),
Здесь μ - постоянный коэффициент диффузии.
Это уравнение второго порядка сводится к системе уравнений первого порядка. Данная система уравнений приводится к полудискретной системе обыкновенных дифференциальных уравнений путем умножения на базисные функции в каждом элементе (отрезке сетки) и интегрирования по элементу по частям (метод Галеркина). Базисными функциями являются полиномы Лежандра внутри каждого элемента. Полученная система обыкновенных дифференциальных уравнений по времени решалась с помощью метода Рунге-Кутта. Рассмотрены различные способы аппроксимации невязкого и вязкого потоков на границах элементов, где численное решение является разрывным. Рассмотренные способы приводят к схемам, как с центральными разностями, так и со смещенными разностями. В последующих параграфах рассмотрены результаты решения конкретных модельных задач и оценка порядка ошибки численного решения. Критерием выбора модельных задач являлось моделирование определенных свойств уравнений газовой динамики (конвекция, вязкость, разрывы решения) и наличие точного решения, что было важным условием для анализа порядка точности решения.
В з1.2 рассмотрено численное решение одномерного уравнения конвекции с начальным условием в виде П-образного импульса. Точное решение соответствует перемещению этого импульса с единичной скоростью без изменения его формы. Результаты расчетов показывают, что диссипативность схемы уменьшается с увеличением порядка полиномов, а точность решения повышается, что демонстрирует преимущества применения численных схем высокого порядка точности.
На примере задачи для того же одномерного уравнения конвекции с гладким начальным условием были проведены исследования порядка ошибки численного решения в С-норме (модуль максимальной разницы между численным и точным решениями). Порядком ошибки называется величина N в асимптотической оценке ошибки Err=O(hN), где h-шаг сетки. Эти исследования показали, что порядок ошибки близок к k+1, т.е. при кусочно-постоянной аппроксимации численного решения (k = 0) получаем схему первого порядка точности, при кусочно-линейной аппроксимации (k=1) - второго порядка и т.д. Эти же результаты свидетельствуют о возможности получения схемы произвольного порядка точности в рамках единообразного подхода.
В з1.3 приведены результаты аналогичных исследований для одномерного уравнения теплопроводности вида с тем же выводом - порядок численной ошибки близок к k+1. Здесь же приведены особенности аппроксимации вязких членов в рамках схемы Галеркина с разрывными функциями.
В з1.4 приведен анализ решения одномерного уравнения Бюргерса с вязкостью . Это уравнение при нулевом коэффициенте вязкости μ дает разрывное решение при гладких начальных условиях и может моделировать решения со скачками уплотнения в газовой динамике. И в этом случае порядок ошибки равен k+1 вне области разрыва решения.
В Главе II рассматриваются вопросы применения схемы DG к решению модельной линейной двумерной задачи конвекции-диффузии, моделирующей криволинейный слой смешения в вязкой жидкости. Проведен анализ порядка ошибки численного решения, показаны преимущества схем высокого порядка точности, продемонстрированы дополнительные возможности, которые обеспечивает схема DG в результате локального выбора базисных функций.
В з2.1 приведена постановка модельной задачи решения линейного уравнения конвекции-диффузии вида в единичном квадрате с разрывным условием вблизи левой границы квадрата. На границах квадрата, где течение втекало в область решения (вектор (1, направлен внутрь области), использовались граничные условия для решения, соответствующие аналитическому решению. Аналитическое решение приведено на рис.1, где красному цвету соответствует значение решения +1, а синему - значение -1. Решение содержит криволинейный слой смешения, форма которого определяется выбором функции f(x), с переходом решения от -1 к +1. Ширина слоя определяется величиной коэффициента вязкости μ. С точки зрения газовой динамики эта задача моделирует развитие вязкого следа за профилем. Численная схема Галеркина с разрывными функциями применялась на неструктурированной сетке из треугольников, которая, на начальном этапе представляла собой триангуляцию структурированной сетки, как показано на рис.1.
Рис.1. Точное решение модельной задачи для двумерного уравнения конвекции-диффузии и расчетная сетка.
Общая теория построения дискретной схемы аппроксимации данного двумерного уравнения с помощью метода Галеркина с разрывными функциями приведена в з2.2. Результирующие уравнения представляют собой систему линейных уравнений с блочной разреженной, вообще говоря, неструктурированной матрицей. Эта система уравнений решалась с помощью обобщенного метода минимальных невязок GMRES с предобуславливателем в виде неполного LU разложения, описанного в з2.3.
В з2.4 приведены результаты численных исследований. На рис.2 представлены в логарифмическом масштабе зависимости ошибки численного решения в норме L2 (средняя интегральная ошибка) и в норме С (максимальная ошибка) уравнения конвекции-диффузии в зависимости от числа степеней свободы (неизвестных переменных) при расчетах на разных сетках. Приведены кривые для различных порядков k полиномиальных базисных функций (отмеченные на рис. как DG(k)). Эти зависимости иллюстрируют одно из преимуществ схем высокого порядка точности, а именно, начиная с некоторого числа переменных, схема высокого порядка точности будет давать меньшую ошибку, чем схема низшего порядка. Анализ порядка ошибки, который равен удвоенному тангенсу угла наклона приведенных зависимостей, показал, что он (порядок) близок к k+1.
Рис. 2. Зависимости ошибок численного решения двумерного уравнения конвекции-диффузии от числа степеней свободы при различных порядках k полиномов базисных функций.
Рассматриваемый метод Галеркина с разрывными функциями обладает большой гибкостью с точки зрения приспосабливания схемы к особенностям решения или расчетной сетки. В частности, в данном параграфе была рассмотрена возможность использования полиномов различных порядков в различных ячейках. Изменение порядка схемы в различных областях может быть использовано, как средство адаптации схемы к решению. Такой способ адаптации в англоязычной литературе называется p-refinement, в отличие от h-refinement - адаптации с помощью изменения размеров и других геометрических характеристик ячеек.
Приведенные ниже результаты расчетов демонстрируют возможность улучшения разрешения особенностей течения (тонкого сдвигового слоя) и повышения точности при наличии плохих сильно вытянутых треугольников путем варьирования порядков полиномов в различных ячейках.
Пример сетки с плохими ячейками приведен на рис. 3. Здесь узлы сетки в одном сечении по оси x (x = 0.55) были смещены по определенному закону. В результате получается искаженная сетка с двумя столбцами плохих ячеек. Три других графика на рис. 3 представляют поля ошибки решений, полученных различными способами. Численное решение, полученное с помощью кусочно-линейной аппроксимации (DG(1), верхний правый график), дает максимальную ошибку 0.52. Отметим, что максимальная ошибка, полученная с DG(1) на неискаженной сетке, равнялась 0.32. Это иллюстрирует влияние плохих ячеек, вытянутых в направлении градиента решения, на его точность. Решение, полученное с помощью кусочно-квадратичной аппроксимации (DG(2), нижний правый график), уменьшает максимальную ошибку до 0.31. Применение DG(1) во всей вычислительной области, за исключением двух столбцов плохих ячеек, где использовался метод DG(2) (квадратичные полиномы), как показано на нижнем левом графике, дало максимальную ошибку, очень близкую к случаю применения метода DG(2) во всей области течения, а именно, 0.32. Решение перед столбцами плохих ячеек очень близко к решению с методом DG(1) во всей области течения, в то время как за этими столбцами - близко к решению методом DG(2) во всей области течения. Таким образом, применение методов более высокого порядка в областях с плохими элементами и более низкого - в остальной области течения позволяет повысить порядок точности решения.
Рис. 3. Исследование влияния локального изменения порядка полиномов в случае сетки с плохими ячейками.
Увеличение порядка полиномов может также применяться локально не только в плохих ячейках, но и в областях с большими градиентами решения, где желательна более высокая точность. Результаты такого рода исследований приведены на рис. 4. Здесь представлены поля величин ошибок для метода DG(1) (верхний левый рис.) и DG(2) (верхний правый рис.). Левый нижний рис. представляет поле ошибки решения, когда метод DG(1) применялся всюду, за исключением ячеек, которые лежат в полосе ширины 0.05 (приблизительно два ряда ячеек) около линии скачка решения. В последнем случае ошибка в норме L2 и максимальная ошибка стали гораздо меньше, чем в случае применения метода DG(1) во всей области (см. таблицу на рис. 4) при небольшом увеличении числа неизвестных (10%). Результаты аналогичного подхода при ширине полосы с квадратичной аппроксимацией решения (DG(2)), равной 0.1, очень близки к результатам при применении метода DG(2) во всей области (которые помечены в таблице "ширина полосы - Q"), в то время как число неизвестных гораздо меньше, чем в случае применения метода DG (2) во всей области.
Рис.4. Влияние локального увеличения порядка полиномов в области с большими градиентами решения.
Еще один пример возможности адаптации метода DG к решению приведен ниже. Если решение медленно меняется в каком-либо направлении, возможно использование полиномов низкого порядка для аппроксимации решения в этом направлении и полиномов более высокого порядка для аппроксимации решения в ортогональном направлении. Это направление может меняться от элемента к элементу. В приведенном ниже примере в каждом элементе вводилось направление градиента решения (N-направление) и нормальное к нему (S-направление). Численное решение аппроксимировалось полиномами по координатам S и N, причем аппроксимация по координате S была кусочно-постоянной, а по координате N - кусочно-квадратичной. Три базисные функции были равны v1=1, v2=N, v3=N2. Таким образом, число неизвестных в таком подходе равно числу неизвестных в стандартном методе DG(1). На рис. 5 приведено сравнение норм ошибок, полученных на равномерных сетках с помощью метода DG(1) (кусочно-линейная аппроксимация), DG(2) (кусочно-квадратичная аппроксимация) и DGm - аппроксимация с помощью модифицированного базиса, описанного выше. Точность, полученная с помощью метода DGm, значительно выше, чем по методу DG(1) при одинаковом числе неизвестных на данной сетке и близка к точности, получаемой по методу DG(2). Такой подход к выбору набора базисных функций можно рассматривать, как способ адаптации схемы к решению.
Рис.5. Сравнение норм ошибок методов DG(1), DG(2), и DGm с адаптированным к решению набором базисных функций.
В Главе III рассмотрены вопросы применения метода Галеркина с разрывными базисными функциями (DG) к решению уравнений Эйлера - уравнений динамики идеального газа.
з3.1 содержит изложение особенностей применения метода DG к решению уравнений Эйлера для плоских течений. Уравнения Эйлера формулируются относительно консервативных переменных в следующем виде:
Здесь консервативные переменные u и декартовы составляющие f(u) и g(u) функции потока F(u) задаются в виде
(1)
где ρ - плотность, u и v - компоненты скорости, p - давление, e - полная энергия единицы массы, h - полная энтальпия единицы массы, γ - отношение удельных теплоемкостей.
Специфичной чертой данного метода является то, что потоки на границах элементов, где численное решение является разрывным, вычисляются с помощью методики Роу приближенного решения задачи о распаде разрыва. Этот прием внесения численной диссипации является для численных схем типа схем Годунова (конечного объема). В остальном технология получения дискретных уравнений близка к изложенной выше. Метод решения дискретных уравнений - неявный метод Эйлера с линеаризацией на верхнем временном слое. Результирующая система линейных уравнений с неструктурированной матрицей решается методом GMRES в безматричной форме с предобуславливателем в виде неполного LU разложения.
В отличие от традиционных методов численного решения задач газовой динамики с помощью схем конечного объема, метод DG оказался очень чувствительным к геометрическому описанию криволинейных границ обтекаемого тела. В з3.2 представлены результаты расчета обтекания кругового цилиндра потоком идеального газа при числе Маха 0.1. Как обычно принято в методе конечного объема, криволинейная граница обтекаемого тела представлена полигоном (многоугольником).
Типичные результаты расчета распределения давления на поверхности цилиндра на сетке 128*64 изображены на рис.6а. Здесь штриховой линией представлено точное решение для несжимаемой жидкости. Сплошная линия соответствует решению, полученному с помощью стандартной, классической схемы конечного объема второго порядка точности, обозначенному как схема КО 2-го порядка. Сплошная линия с кружочками соответствует расчету по исследуемой DG(1) схеме с аппроксимацией решения с помощью кусочно-линейных функций. Решение методом DG(1) носит осциллирующий характер. Отметим, что это решение соответствует невязке, уменьшенной до машинного нуля, т.е. полностью сошедшееся решение. В классической схеме решение плавное и без осцилляций и близко к точному решению, в то время как решение с использованием DG схемы сильно осциллирующее. Отметим, что мы имеем большие ошибки в градиентах решений в каждом элементе, в то время как, в среднем, решение довольно близко к точному.
Решение осциллирует только вблизи поверхности цилиндра. На рис.6а приведено решение по DG(1) схеме на расстоянии R=1.165 от центра цилиндра (радиус самого цилиндра равен 1). Видно, что в этом случае осцилляций нет. Таким образом, можно предположить, что причина плохого решения по DG схеме лежит в неадекватном представлении границы обтекаемого тела.
а) | б) |
Рис.6. Распределение давления на поверхности цилиндра при М=0.1. Сетка 128*64.
Возможным способом решения этой проблемы является использование вблизи криволинейной границы т.н. параметрических элементов с одним криволинейным ребром. Теоретические аспекты применения таких элементов содержатся в з3.3, а результаты расчетов - в з3.4. Результат расчета распределения давления по поверхности цилиндра на той же сетке, что и на рис.6а, приведен на рис.6б. Видно, что применение параметрических элементов обеспечивает высокое качество решения. Схема DG(1) с кусочно-линейной аппроксимацией решения обеспечивает также меньшую диссипацию, чем схема конечного объема второго порядка на той же расчетной сетке. Использование схемы КО приводит к несимметричному обтеканию, являющемуся результатом возникшего отрыва потока в задней части цилиндра, в то время как картина течения, полученная по DG схеме, выглядит более симметричной, в соответствии с точным решением.
Были также проведены исследования порядка точности различных схем путем анализа поведения коэффициента сопротивления на последовательности сеток с удваивающимся количеством узлов в каждом направлении. Результаты представлены на рис.7, где показана зависимость коэффициента сопротивления от количества неизвестных N в логарифмическом масштабе для a) схемы конечного объема первого порядка (кусочно-постоянная реконструкция, схема Годунова), b) схемы конечного объема второго порядка (кусочно-линейная реконструкция), c) DG(0) схемы (кусочно-постоянные элементы), d) DG(1) схемы (кусочно-линейные элементы), e) DG(2) схемы (кусочно-квадратичные элементы). Для всех DG схем расчеты были выполнены с квадратическими параметрическими элементами на границе. Как и ожидалось, обе схемы первого порядка демонстрируют первый порядок ошибки в сопротивлении. Обе схемы второго порядка дают очень близкие значения сопротивлений при одном и том же количестве неизвестных, в то время как порядок величины ошибки асимптотически стремится к 3. Схема DG(2) дает порядок ошибки, асимптотически стремящийся к 5 (на рис.7 приведены треугольники, углы наклона гипотенузы которых соответствуют указанному порядку точности схемы). Такие высокие порядки величин могут быть следствием высокой регулярности используемых сеток. Следует отметить высокую точность метода DG(2) - величина коэффициента сопротивления составляет 10-5, что весьма трудно получить при использовании схем меньшего порядка точности при разумном числе неизвестных задачи.
Рис.7. Коэффициент сопротивления в зависимости от количества неизвестных N в случае обтекания цилиндра при М=0.1
Как указывалось выше, схема DG на неструктурированных сетках обеспечивает широкие возможности в плане адаптации к решению, в частности, адаптации расчетной сетки. При этом автоматически генерируется расчетная сетка с минимальным количеством узлов при обеспечении требуемой точности решения. В з3.5 на ряде примеров иллюстрируется применение метода DG на адаптивных сетках для решения уравнений Эйлера. Методика адаптации сеток, состоящих из треугольных ячеек, разработана А.М.Сорокиным. Критерием адаптации (дробления ребра сетки), является превышение изменения числа Маха вдоль данного ребра некоторого порогового значения. Эта операция не применяется к ребрам, ориентированным вдоль градиента решения. Такой подход способствует, в целом, вытягиванию ячеек вдоль линий уровня параметра адаптации, что является положительным свойством сетки. В процессе адаптации возможно формирование избыточно мелкой сетки. Для борьбы с этим применяется процедура укрупнения ячеек. Такая процедура адаптации применяется последовательно с расчетом течения, т.е. осуществляется расчет на крупной сетке, затем применяется процедура адаптации, которая приводит к увеличению числа ячеек. Это составляет одну большую итерацию адаптации.
Ниже приведены два примера решения уравнений Эйлера на адаптивных сетках. На рис.8 представлено решение задачи о сверхзвуковом обтекании биплана Буземана при числе Маха невозмущенного потока M = 2.766. Начальная сетка около данной конфигурации состояла из 1741 ячеек. На рис. представлена сетка с 9224 ячейками и решение после 7 итераций адаптации сетки. Данный пример иллюстрирует возможность получения решения по методу DG на весьма нерегулярных сетках, а также существенное улучшение качества решения (разрешения скачков и волн разрежения) в результате адаптации сетки к решению.
Рис. 8. Сетка и решение (поле числа Маха) после 7 итераций адаптации сетки к решению для биплана Буземана. Число узлов сетки 4695, 9224 ячеек.
Еще один пример расчета на адаптивных сетках представляет собой расчет обтекания профиля NACA0012 при M=0.85, α=1.25. Использовался метод DG(1) с линейной реконструкцией решения и было выполнено 17 итерации адаптации сетки. Изменение коэффициента подъемной силы Cy на последних итерациях не превышает 0.7%, коэффициента сопротивления Cx - 2%. На рис. 9 приведены сетка, распределение давления и поле чисел Маха при расчете на сетке, полученной в результате адаптации с числом элементов 22111. Наблюдается хорошее качество решения, высокое разрешение скачков уплотнения, в частности, слабого скачка на нижней поверхности профиля.
Рис. 9. Сетка после 17-ти итераций адаптации и решение (коэффициент давления и поле числа Маха) для профиля NACA0012, M=0.85, α=1.25.
В Главе IV рассматривается решение уравнений Навье-Стокса для плоских ламинарных течений с использованием конечно-элементного метода Галеркина (DG) с разрывными базисными функциями.
В з4.1 изложены особенности применения метода DG к решению уравнений Навье-Стокса, которые записываются для вектора консервативных переменных u в следующей форме:
Консервативные переменные u и декартовы составляющие fi(u) и gi(u) определяются согласно (1). Декартовы компоненты fv(u) и gv(u) функции вязкого потока Fv(u) задаются в виде
где μ - коэффициент динамической вязкости, Pr - число Прандтля, λ=-2/3.
Наличие вязких напряжений в уравнениях, которые выражаются через градиенты консервативных переменных, приводит к необходимости введения дополнительного вектора переменных, равного градиенту консервативных переменных и соответствующему расширению системы уравнений. Вязкие потоки на границах элементов, где численное решение разрывное, определяются как полусумма значений в элементах, примыкающих к данному ребру. В остальном технология получения дискретных уравнений близка к изложенной выше. Метод решения дискретных уравнений - неявный метод Эйлера с линеаризацией на верхнем временном слое. Результирующая система линейных уравнений с неструктурированной матрицей решается методом GMRES в безматричной форме с предобуславливателем в виде неполного LU разложения.
В з4.2 приведены примеры решения уравнений Навье-Стокса с помощью метода DG с различным порядком точности схемы. Рассмотрены задачи обтекания плоской пластины при числе Рейнольдса Re=106 при малых скоростях M=0.1, кругового цилиндра при Re=40 при М=0.1, а также два случая ламинарного обтекания профиля NACA0012 для режимов М=0.5, α=0, Re=5000 и М=0.8, α=10, Re=73. В качестве примера на рис.10 приведены результаты расчета кругового цилиндра при Re=40, M=0.1 на достаточно редкой сетке 20*10 узлов при различных порядках базисных полиномов k в виде картины линий тока. Данный режим обтекания по экспериментальным данным и результатам расчетов других авторов является стационарным и характеризуется наличием развитой отрывной зоны за цилиндром. Представленные на рис.10 данные демонстрируют симметризацию картины течения и улучшение описания отрывной зоны по мере увеличения порядка базисных функций от 0 до 4.
В таблице 1 приведено общее количество переменных (степеней свободы) при различных сочетаниях сеток и величин порядка базисных полиномов k. В таблице 2 приведено угловое значение положения точки отрыва для различных сеток и различных значений k. Здесь же дано экспериментальное значение этого угла 53.5, приведенное в работе (Coutanceau, Bouard, 1977), а также расчетное значение из работы (Hauke, Hughes, 1998). Анализ результатов, представленных на рис. 10, а также в таблицах 1, 2 показывает, что решения близкого качества и степени согласования с экспериментом могут быть получены на мелкой сетке 80*40 при k=1 и на крупной сетке 20*10 при k=3 или 4. В последних случаях общее число переменных в 3-5 раз меньше, чем в первом, что свидетельствует о преимуществах применения схем высокого порядка точности.
В з4.3. приведен пример решения уравнений Навье-Стокса с использованием метода DG с адаптацией сетки к решению. Методика адаптации неструктурированной сетки, описанная в главе 3, была применена к решению задачи вязкого обтекания профиля NACA0012 при двух режимах: дозвуковое обтекание при числе Маха М = 0.5, α=0, Re=5000 и сверхзвуковое обтекание при числе Маха M=2.0, α=10, Re=106. Расчеты данных случаев обтекания приведены также в работе (Bassi, Rebay, 1997). Использовался метод DG(1) с линейной реконструкцией решения. Сетка адаптировалась к числу Маха. На рис 11 приведены сетка и поле чисел Маха для второго из указанных случаев после 14 итераций адаптации сетки к решению. Наблюдается, в частности, хорошее разрешение головной ударной волны.
Рис.10 Круговой цилиндр. Сетка 20*10. Re=40, M=0.1. Линии тока.
Таблица 1. Общее число переменных.
сетка | k=0 | k=1 | k=2 | k=3 | k=4 |
20*10 | 400 | 1200 | 2400 | 4000 | 6000 |
40*20 | 1600 | 4800 | 9600 | 16000 | |
80*40 | 6400 | 19200 | 38400 |
Таблица 2. Круговой цилиндр, Re=40. Угловое положение точки отрыва.
сетка | k=0 | k=1 | k=2 | k=3 | k=4 |
20*10 | 36.0 | 48.1 | 54.0 | 53.2 | 53.4 |
40*20 | 45.0 | 55.6 | 54.0 | 53.7 | |
80*40 | 49.2 | 54.4 | 53.8 | ||
Эксперимент (Coutanceau, Bouard, 1977) | 53.5 | ||||
Расчет (Hauke, Hughes, 1998) | 54.0 |
Рис. 11. Сетка после 14 итераций адаптации и решение (поле числа Маха) для профиля NACA0012, M=2.0, =10, Re=106.
В Главе V рассматривается применение метода Галеркина с разрывными базисными функциями к расчету вязких турбулентных плоских течений, описываемых уравнениями Навье-Стокса, осредненными по Рейнольдсу (уравнения Рейнольдса). Для замыкания системы уравнений Рейнольдса используется модель турбулентности Спаларта-Алмараса (Spalart, Allmaras, 1992).
В з5.1 описаны рассматриваемые уравнения Рейнольдса, кратко описана модель турбулентности Спаларта-Альмараса и особенности применения метода DG в данном случае. Уравнения Рейнольдса и уравнение Спаларта-Алмараса могут быть записаны совместно в консервативной форме. Эта система уравнений имеет следующий вид:
(2)
Консервативные переменные u и декартовы составляющие fi(u) и gi(u) невязкого потока Fi(u) задаются в виде
Здесь, как и ранее, ρ - плотность, u и v - компоненты скорости, p - давление, e и h - полная энергия и энтальпия единицы массы, - рабочая переменная в уравнении Спаларта - Алмараса, определяющая величину турбулентной вязкости.
Декартовы компоненты fv(u) и gv(u) вязкого потока Fv(u) задаются в виде
Источниковый член S в уравнении (2) содержит дополнительные слагаемые только в уравнении для турбулентной вязкости, связанные с производством, диффузией и диссипацией турбулентной вязкости, в соответствии с работой (Spalart, Allmaras, 1992).
Здесь компоненты тензора вязких напряжений и вектора теплового потока Qx, Qy определяются обычным образом в соответствии с гипотезой Буссинеска:
, , ,
,
Составляющие градиента удельной внутренней энергии Е равны
,
Параметры модели турбулентности определяются в соответствии с рекомендациями авторов модели, которые в деталях воспроизведены в диссертации.
Наличие дополнительного уравнения переноса турбулентной вязкости потребовало некоторой модификации стандартной схемы Роу решения задачи о распаде разрыва. В остальном технология получения дискретных уравнений совпадает с изложенной методикой решения ламинарных уравнений Навье-Стокса.
В з5.2 приведены некоторые примеры решения уравнений Рейнольдса помощью метода DG. Рассмотрены задачи турбулентного обтекания плоской пластины при числах Рейнольдса Re=106, 107 и малых скоростях M=0.1, кругового цилиндра при Re=3.6⋅106 и малых скоростях М=0.1, а также турбулентного обтекания профиля RAE2822 при M=0.725, =2.92, Re=6.5⋅106 и M=0.725, =3.19, Re=6.5⋅106, соответствующих безотрывному и отрывному режимам обтекания.
Как указывалось выше, одним из потенциальных преимуществ схемы DG является возможность варьирования порядка точности схемы в различных областях течения (изменение порядка полиномов, аппроксимирующих решение), а также адаптации схемы к решению с помощью изменения указанного порядка (p-refinement). Повышение порядка точности схемы будет наиболее эффективным в областях, где решение является гладким, но быстро меняющимся, например, в пограничных слоях. В качестве первого шага в этом направлении рассмотрено численное решение задачи о вязком турбулентном обтекании профиля RAE2822 при условиях M=0.725, =2.92, Re=6.5⋅106 (CASE-6).
На рис.12 представлены профили скорости в пограничном слое в трех сечениях на верхней поверхности профиля. Представлено 3 решения. Решение №1 получено на грубой сетке из 3660 ячеек по схеме DG(1) 2-го порядка точности (10980 неизвестных задачи). Данное решение не обладает достаточной точностью. Решение №2 получено на измельченной сетке из 6448 ячеек по той же схеме (19334 неизвестных). Результаты проведенных расчетов на различных сетках позволяют утверждать, что это решение является достаточно точным. Решение №3 получено на той же сетке, что и решение №1, однако в области пограничного слоя использовалась схема повышенного 3-го порядка. Число неизвестных в этом случае равнялось 13869, т.е. меньше, чем в расчете№2, а формы профилей скорости было получены с высокой точностью и практически совпадают с профилями скорости, полученными в расчете №2. Этот пример является иллюстрацией эффективности использования схемы повышенного порядка в областях гладкого быстро меняющегося решения, например в пограничном слое.
Рис 12. Профили скорости в пограничном слое.
В качестве примера расчета обтекания тела сложной геометрии ниже приведены результаты расчета обтекания трехэлементного профиля корпорации Mcdonnell Douglas (MCD). Данный профиль представляет собой сечение крыла с отклоненными на 30 предкрылком и закрылком в посадочной конфигурации. Этот профиль был подробно экспериментально исследован в малотурбулентной аэродинамической трубе НИЦ им. Лэнгли НАСА, и экспериментальные результаты, использованные ниже, взяты из работ [Chin, и др., 1993, Spaid, Lynch, 1996]. Общий вид профиля представлен на рис.16. Эксперимент был проведен при числе М=0.2 и числе Re=9⋅106. Расчеты были проведены по схеме DG(1) с линейной реконструкцией решения в ячейках, с моделью турбулентности Спаларта-Альмараса и применением методики адаптации сетки к решению, описанной выше.
На рис.13 приведены исходная (на первой итерации) и финальная (на седьмой итерации) расчетные сетки. Видно, что сгущение и обогащение сетки происходят в областях пограничных слоев и вязких следов. Результаты расчета поля числа Маха и его фрагментов приведены на рис.14. Обтекание характеризуется наличием заметных отрывов потока на нижней поверхности предкрылка и в нише основного профиля. Наблюдается хорошее разрешение вязких следов, сходящих с различных элементов профиля, что свидетельствует о достаточной эффективности процедуры адаптации. Данная процедура позволяет также разрешать и другие особенности течения, такие, например, как скачки уплотнения. Это иллюстрируется данными, представленными на рис.15, где приведены фрагменты сетки и поля чисел Маха вблизи передней кромки предкрылка при α=24, что соответствует критическому углу атаки (максимальной подъемной силе). Данный режим характерен наличием малой местной сверхзвуковой зоны на передней кромке предкрылка, с максимальным местным числом Маха М=1.3 и замыкающим скачком уплотнения. Эта тонкая особенность течения также хорошо разрешена в результате адаптации сетки к решению.
Рис.13. Расчетные сетки около профиля MCD на первой и последней итерации адаптации сетки к решению.
Рис.14. Поле числа Маха при обтекании профиля MCD при M=0.2, α=18, Re=9⋅106.
Рис.15. Фрагменты адаптивной сетки и поля чисел Маха вблизи передней кромки предкрылка профиля MCD при M=0.2, α=24, Re=9⋅106.
Экспериментальные исследования обтекания профиля MCD включали измерения профиля полного давления в сечении x/c=0.975 при трех значениях угла атаки (α=8, 16, 21). На рис.16 результаты измерения коэффициента полного давления сравниваются с результатами расчета в сечении, обозначенном буквой А. Наблюдается хорошее согласование как по положению областей дефекта полного давления, обусловленных вязкими следами, сходящими с элементов профиля, расположенных выше по потоку, так и по величине этого дефекта. Следует отметить, что данная характеристика является довольно тонкой, и согласование результатов свидетельствует о хорошем качестве решения.
С практической точки зрения представляющей интерес характеристикой профиля с взлетно-посадочной механизацией является зависимость подъемной силы от угла атаки и величина максимальной подъемной силы. На рис.17 приведено сравнение этих зависимостей, полученных экспериментально [van Dam, 2002] и в расчете. Наклоны этих зависимостей на линейном участке и величины критического угла атаки согласуются хорошо. Разница в величине максимальной подъемной силы составляет около 2%.
Рис.16. Сравнение экспериментальных и расчетных профилей коэффициента полного давления на закрылке профиля MCD при x/c=0.975, α=8, 16, 21.
Рис.17. Сравнение расчетной и экспериментальной зависимостей коэффициента подъемной силы от угла атаки трехэлементного профиля MCD.
Для иллюстрации возможности применимости рассматриваемого метода к расчету пространственных течений, в работе приведен пример расчета трехмерного обтекания крыла LANN при М = 0.82, Re = 7.3⋅106, α = 0.6. Течение считалось полностью турбулентным. Расчетная сетка была построена с помощью промышленного пакета прикладных программ и она состояла из 190213 гексаэдральных элементов. Расчет проводился по схеме DG(1) с кусочно-линейной реконструкцией решения в ячейках, что дает общее количество неизвестных задачи равное 760852. На рис.18 представлено распределение коэффициента давления по верхней поверхности крыла и сравнение с экспериментальными данными для сечения z = 0.475. Наблюдается хорошее согласование результатов расчета с экспериментом.
Рис.18. Сравнение распределения коэффициента давления для крыла LANN в сечении z = 0.475 с экспериментальными данными.
Основные результаты работы и выводы
- На основе перспективного метода дискретизации законов сохранения - метода Галеркина с разрывными функциями - предложен новый подход к решению уравнений газовой динамикиЦ уравнений Эйлера, Навье-Стокса и Рейнольдса. Отличительными чертами данного подхода являются:
- возможность построения схем произвольного порядка точности унифицированным образом;
- удобство работы с неструктурированными сетками;
- широкие возможности адаптации схемы к сетке и/или решению путем локального повышения порядка точности в проблемных областях (лплохие ячейки, области больших градиентов и т.п.);
- адаптация путем локальной модификации набора базисных функций;
- автоматическая адаптация сетки к решению.
Центральным отличием предлагаемого подхода от других способов численного решения уравнений газовой динамики является возможность адаптации схемы к решению путем локального изменения порядка точности численной схемы.
- На ряде примеров решений уравнений Эйлера, Навье-Стокса и Рейнольдса показаны преимущества адаптивных схем высокого порядка, позволяющие уменьшить время расчета при заданной точности в несколько раз по сравнению с классическими методами типа метода Годунова. Продемонстрировано успешное применение процедуры адаптации неструктурированной сетки к решению, обеспечившее высокое разрешение различных особенностей течения (область диффузорного отрыва, головная ударная волна, локальные особенности обтекания многоэлементных профилей)
- Метод описан, протестирован на многочисленных примерах и подготовлен к практическому применению.
Основные результаты опубликованы в работах
- А.В.Волков, С.В.Ляпунов. Исследование эффективности использования численных схем высокого порядка точности для решения уравнений Навье-Стокса и Рейнольдса на неструктурированных адаптивных сетках. ЖВМиМФ, т.46, №10, 1894-1907, 2006.
- А.В.Волков, С.В.Ляпунов. Применение конечноэлементного метода Галеркина с разрывными базисными функциями к решению уравнений Рейнольдса на неструктурированных сетках. Ученые записки ЦАГИ, , т.XXXVIII, №3-4, 2007.
- S.V.Lyapunov, A.V.Wolkov. Application of Discontinuous Galerkin finite element method to the solution of partial differential equations. Part I. 2D scalar conservation laws. 16th IMACS World Congress. Lausanne-Switzerland. August 21-25, 2000.
- A.V.Wolkov, S.V.Lyapunov. Application of Discontinuous Galerkin finite element method to the solution of partial differential equations. Part II. System of nonlinear equations: Euler equations. 16th IMACS World Congress. Lausanne-Switzerland. August 21-25, 2000.
- N.B.Petrovskaya, A.V.Wolkov, S.V.Lyapunov. Modification of basis functions in high order discontinuous Galerkin schemes for advection equations. University of Birmingham, School of mathematics, Preprint 2006/26, 2006.
- N.B.Petrovskaya, A.V.Wolkov, S.V.Lyapunov. Modification of basis functions in high order discontinuous Galerkin schemes for advection equations. Applied Mathemetical Modelling, Elsevier, 2006.
- V.S.Sakovich, A.M.Sorokin, A.V.Wolkov, S.V.Lyapunov. Anisotropic unstructured grid generation for 3D flow simulation problems. 6th Int. Conf. on numerical grid generation in computational fluid simulation, 1998.
- А.В.Волков, С.В.Ляпунов, А.Н.Храбров. Влияние установившегося вращения на аэродинамические характеристики профиля при наличии отрыва потока. Ученые записки ЦАГИ, т.XXXIV, №3-4,2003
- S.V.Lyapunov, A.V.Wolkov. A separated flow calculation about airfoils and high-lift systems on basis of viscous-inviscid interaction methods. EUROMECH, 3rd European Fluid Mechanics Conf., 1997
- S.V.Lyapunov, A.V.Wolkov. Calculation of separated flows about airfoils and high-lift systems using viscous-inviscid interaction approach. Proc. Of the 5th russian-chinise simposium on aerodynamics and flight dynamics, 1997
- S.V.Lyapunov, A.V.Wolkov Application of Viscous-Inviscid Interaction Methods for a Separated Flow Calculation About Airfoils and High-Lift Systems. ICAS Proceedings, ICAS-96-1.10.2,1996
- S.V.Lyapunov, A.V.Wolkov. Application of the Viscous-Invisid Interaction Model to Calculation of Two-Dimensional Separated Flows. TsAGI Journal, vol.2, 1,1996
- S.V.Lyapunov, A.V.Wolkov. Numerical Prediction of Transonic Viscous Separated Flow Past an Airfoil. Theoretical and Computational Fluid Dynamics, vol.6, №1, 1994
- А.В.Волков, С.В.Ляпунов. Метод расчета трансзвукового вязкого обтекания профиля с учетом изменения энтропии на скачках уплотнения. Ученые записки ЦАГИ, т.XXIV, №1, 1993
- С.В.Ляпунов. Неплоские крылья минимального индуктивного сопротивления. Изв.РАН, МЖГ, №2,1993
- S.V.Lyapunov. Accelerated Method of the Euler Еquation Solution in Transonic Airfoil Flow Problem. ICAS Proceedings, ICAS-92-4.2.3,1992
- С.В.Ляпунов. Ускоренный метод решения уравнений Эйлера в задаче о трансзвуковом обтекании профиля. Мат.моделирование, N4, 1991
- S.V.Lyapunov. Convergence acceleration and wave drag determination in transonic airfoil calculations. ICAS Proc. ж ICAS-90-6.9.2, 1990
- N.A.Vladimirova, V.V.Vyshinsky, S.V.Lyapunov, Ya.M.Serebriisky Convergence acceleration of computational methods for two- and three-dimensional flows about bodies. Fluid Mechanics, v.2, 1986
- Н.А.Владимирова, В.В.Вышинский, С.В.Ляпунов, Я.М.Серебрийский Об ускорении сходимости методов расчета плоского и пространственного трансзвукового обтекания тел в неограниченном потоке. Ученые записки ЦАГИ, т.16, №4, 1985
- М.А.Брутян, С.В.Ляпунов. О второй вариации функционала Рябушинского. ДАН СССР, т.258, №4, 1981
- В.В.Вышинский, С.В.Ляпунов. Расчет околозвукового осесимметричного обтекания мотогондолы с учетом вязкости. Ученые записки ЦАГИ, т.19, №3, 1988
- А.В.Волков, В.С.Сакович, А.М.Сорокин, Ю.Б.Лифшиц, С.В.Ляпунов. Применение неструктурированных сеток в вычислительной аэродинамике. Материалы IХ школы-семинара УАэродинамика летательных аппаратовФ, ЦАГИ, 1998.
- А.В.Волков, С.В.Ляпунов. Применение ориентированных схем второго порядка для решения двумерных уравнений Эйлера на неструктурированных сетках. Материалы Х школы-семинара УАэродинамика летательных аппаратовФ, ЦАГИ, 1999.
- А.В.Волков, С.В.Ляпунов. Исследование схемы аппроксимации высокого порядка уравнений газовой динамики. Материалы ХI школы-семинара УАэродинамика летательных аппаратовФ, ЦАГИ, 2000.
- А.В.Волков, С.В.Ляпунов. Исследование схемы аппроксимации высокого порядка уравнений газовой динамики. Современные проблемы аэрокосмической науки и техники. Международная научно-техническая конференция молодых ученых и специалистов, 2000.
- А.В.Волков, С.В.Ляпунов. Применение схемы высокого порядка точности для решения двумерных уравнений Эйлера и Навье-Стокса на неструктурированных сетках. Материалы ХII школы-семинара УАэродинамика летательных аппаратовФ, ЦАГИ, 2001.
- А.В.Волков, С.В.Ляпунов. Численный метод высокого порядка точности для решения системы уравнений Навье-Стокса и Рейнольдса. Материалы ХIII школы-семинара УАэродинамика летательных аппаратовФ, ЦАГИ, 2002.
- А.В.Волков, С.В.Ляпунов. Разработка новой монотонной схемы высокого порядка точности на компактном шаблоне для решения уравнений Эйлера и Навье-Стокса. Материалы ХVI школы-семинара УАэродинамика летательных аппаратовФ, ЦАГИ, 2003.
- А.В.Волков, С.В.Ляпунов. Применение метода конечных элементов к расчету вязких турбулентных течений на неструктурированных адаптивных сетках. Материалы ХV школы-семинара УАэродинамика летательных аппаратовФ, ЦАГИ, 2004.
- А.В.Волков, С.В.Ляпунов. Разработка новой монотонной схемы высокого порядка точности на компактном шаблоне для решения уравнений Эйлера и Навье-Стокса. Материалы ХVI школы-семинара УАэродинамика летательных аппаратовФ, ЦАГИ, 2005.
- А.В.Волков, С.В.Ляпунов. О решении акустической задачи методом Галеркина с разрывными базисными функциями. Материалы ХVII школы-семинара УАэродинамика летательных аппаратовФ, ЦАГИ, 2006.