Версия статьи с цветными иллюстрациями размещена по адресу http://www.keldysh.ru/pages/cgraph/publications/cgd_publ.htm.

Содержание

1. Введение.......................................................................................................... 4

2. Гибридные разностные схемы....................................................................... 6

3. Анализ разностной схемы.............................................................................. 7

3.1. Пример гибридной разностной схемы.................................................... 7

3.2. Постановка задачи................................................................................... 8

3. Организация расчетов.............................................................................. 10

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

5. Обсуждение................................................................................................... 14

6. Заключение................................................................................................... 14

Список литературы........................................................................................ 166

1. Введение

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

Проблемы, задачи и пути развития научной визуализации были  подробно представлены в работах [3,4,5,11,13] на примере такого важного раздела математического моделирования как вычислительная аэрогазодинамика.

Роль методов и алгоритмов визуального представления результатов математического моделирования (и, в частности, численного моделирования в вычислительной аэрогазодинамике)  изначально рассматривалась как вспомогательная. Визуализация имела две основные вспомогательные функции: 

а) обеспечение контроля и лучшего понимания численных результатов;

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

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

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

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

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

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

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

Рассматриваемые ниже разностные схемы относятся к классу гибридных разностных схем.

2. Гибридные разностные схемы

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

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

В простом случае гибридную схему можно записать как комбинацию GS1 + (1-G)S2, где G – коэффициент гибридности, S1 и S2 – разностные схемы, обладающие различными интересующими исследователя свойствами. Например, S1 – схема первого порядка точности, а S2 – второго порядка.

Большинство применяемых для решения практических задач аэрогазодинамики являются гибридными. Согласно [6], к  гибридным схемам относятся такие широко известные алгоритмы как FCT (flux corrected transport), различные типы TVD (total variation diminishing) разностных схем, схемы типа ENO (essentially non-oscillatory) и WENO (weighted essentially non-oscillatory)  и многие другие. Подробное описание и классификация различных типов гибридных разностных схем приведена в [6]. Таким образом, использование гибридных схем чрезвычайно полезно и удобно, так как позволяет исследователю использовать наилучшие свойства разных схем. В то же время возникают и определенные трудности – надо тщательно изучать свойства и возможные ограничения коэффициентов гибридности (или весовых коэффициентов) для того, чтобы используемое свойство соответствовало физической модели рассматриваемой задачи.

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

3. Анализ разностной схемы

            3.1. Пример гибридной разностной схемы

Рассматривается WW (with weights) – разностная схема, разработанная автором в  1990-х годах [1,2,7]. Приведем ее на примере уравнения Бюргерса, являющегося одномерным аналогом уравнений Навье-Стокса:

                                               (3.1)

Введем обозначения для разностей

;    

Тогда для уравнения (3.1) WW- схема записывается следующим образом:

                            (3.2)

где  - весовые коэффициенты;  - шаг по времени и пространственной переменной.

Нетрудно заметить, что при выборе  разностная схема (3.2) является линейной комбинацией схемы Крэнка-Никольсона, обладающей вторым порядком точности по времени и пространству, и схемы Лакса, обладающей существенной искусственной вязкостью.

При выборе  добавочный член с искусственной вязкостью записывается как

     .                                             (3.3)

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

Данная разностная схема успешно использовалась для решения широкого круга модельных и практических задач вычислительной аэрогазодинамики [7,8,9,10,12,14].

Рассмотрим WW-схему на примере задачи течения в дальнем следе.

3.2. Постановка задачи

В расчетной области ABCD (рис.1) рассматривается течение вязкого сжимаемого теплопроводного газа, описываемое полной системой нестационарных уравнений Навье-Стокса. На границе AB задаются граничные условия симметрии. На границах DC и CB задаются «мягкие» граничные условия типа линейной  экстраполяции.

.       

                  Рис.1

На входной границе AD задаются распределения газодинамических параметров, полученные из расчетов обтекания осесимметричного тела и участка следа за ним. Эти начальные распределения давления P и скорости U представлены на рис.2,3.

                 Рис.2                                                 Рис.3

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

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

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

3.3. Организация расчетов.

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

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

                                                (3.4) 

 Иначе говоря, для каждого набора  нужно решить обратную задачу, варьируя   до тех пор, пока не будет выполнено условие (3.4).

Схема расчета выглядит следующим образом.  При  для каждого сеточного разбиения  решается прямая задача моделирования течения в дальнем следе с помощью WW- схемы при некотором заданном начальном значении . В счетной области определяется значение  - количество локальных экстремумов. Далее решается классическая обратная задача путем вариации  до выполнения условия (3.4). Одновременно проводится в режиме online визуальное представление  в виде поверхности .

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

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

Полученная в результате расчетов поверхность значений предельных весовых коэффициентов  представлена на рис.4. Эти же данные представлены в виде изолиний на рис.5.

                                        Рис.4

                                           Рис.5

Анализируя вид данных, представленных на рис.4,5, естественно применить преобразование данных и представить их в виде   . Вид поверхности после применения преобразования представлен на рис.6,7.   

                                             Рис.6

                                             Рис.7

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

                                          (4.1)

Или учитывая что  , ,  выражение (4.1) можно представить как

                                           (4.2)

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

5. Обсуждение

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

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

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

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

6. Заключение

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

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

Иллюстрации получены при помощи комплекса GeoGraph, разработанного под руководством В.Н.Кочина, которому автор выражает свою признательность.

Список литературы.

[1] А.Е.Бондарев «Численное решение уравнения Бюргерса в области высоких градиентов» Препринт  ИПМ им М.В.Келдыша РАН, М., № 12, 1990, 13 с.

[2] А.Е.Бондарев «Разработка метода численного исследования отрывных течений вязкого газа». Отчет ИПМ им. М.В.Келдыша РАН, N 170-90, 1990, 16 с.

[3] А.Е.Бондарев «Визуализация в вычислительной аэрогазодинамике: проблемы развития». Тр. ГРАФИКОН-99, М., МГУ, 1999, с.9-13.

[4] А.Е.Бондарев, Е.Н.Бондарев «Функции визуализации в вычислительной аэрогазодинамике» Общероссийский научно-технический журнал «Полет», М., «Машиностроение»,  N 10, 2000, с.53-60.

[5] А.К.Алексеев, А.Е.Бондарев, Ю.А.Молотилин «О визуализации сопряженных полей при идентификации и управлении трехмерным течением вязкой жидкости». Сб. «Применение методов научной визуализации в прикладных задачах», М., МГУ, 2000, с.8-18.

[6] А.Г.Куликовский,Н.В.Погорелов, А.Ю.Семенов «Математические вопросы численного решения гиперболических систем уравнений» М, Физматлит, 2001, 608 с.

[7] А.Е.Бондарев «Влияние параметров сверхзвукового потока на характерное время установления течения при обтекании уступа» Изв.АН СССР, МЖГ, N 4, 1989, с.137-140.

[8] А.K.Alexeev, A.E.Bondarev, Y.A.Molotilin  « On Inverse Problems for 3D Time-Dependent Free Convection Heat Transfer»  Proc. of  National  Heat Transfer Conference ASME,  v.10, 1995, p.113-122.

[9]  Alexeev, А.K., Bondarev, A. E., Molotilin, Y.A. «On Boundary Condition Determination Using Inverse Free Convection Problem.»  Proc. of  National Heat Transfer Conference ASME,  Oregon, USA, v.10, 1995.

[10] А.К.Алексеев, А.Е.Бондарев «Визуализация переноса погрешности при расчете поля течения» Сб. «Научная визуализация в прикладных задачах», М., МГУ, 2003, с.4-13.

[11]  А.Е. Бондарев «Конструирование и оптимизация разностных схем с применением методов визуализации» Тр. ГРАФИКОН-2003, М., МГУ, 2003, с.261-264.

[12] А.Е. Бондарев, Ю.А. Молотилин «Численное моделирование дифракционных эффектов при распространении звуковых волн» Препринт  ИПМ им М.В.Келдыша РАН, М., № 10, 1995, 13 с.

[13] А.Е.Бондарев, Е.Н.Бондарев «О трассировке вихревых структур» Сб. «Научная визуализация в прикладных задачах», М., МГУ, 2003, с.40-45.

[14] А.К.Алексеев, А.Е.Бондарев, Ю.А.Молотилин «Обратные задачи на основе уравнений Навье-Стокса для задач свободной конвекции» Препринт  ИПМ им М.В.Келдыша РАН, М., № 32, 1995, 15 с.

[15] Н.А. Черанева «Неявная разностная схема на неравномерной сетке для решения уравнений Навье-Стокса сжимаемого газа» Сб. «Численные методы эродинамике» Вып.4, НИВЦ МГУ, Изд.МГУ, 1980, сс. 3-19.