Выбор весовых функций в методе Петрова-Галеркина для интегрирования двумерных нелинейных уравнений типа Бюргерса Текст научной статьи по специальности «Математика»

УДК 519.6

А.А. МОЛЧАНОВ, С.В. СИРИК, Н.Н. САЛЬНИКОВ

ВЫБОР ВЕСОВЫХ ФУНКЦИЙ В МЕТОДЕ ПЕТРОВА-ГАЛЕРКИНА ДЛЯ ИНТЕГРИРОВАНИЯ ДВУМЕРНЫХ НЕЛИНЕЙНЫХ УРАВНЕНИЙ ТИПА БЮРГЕРСА

Анотація. Запропонована модифікація методу Петрова-Гальоркіна для чисельного розв'язання нестаціонарних двовимірних рівнянь з нелінійними конвективними членами типу Бюргерса, заснована на використанні вагових функцій, які змінюються у часі в залежності від виду еволюціонуючого розв'язку. Ефективність запропонованої модифікації для задач з домінуванням конвективних процесів над дифузійними підтверджується розрахунками.

Ключові слова: метод Петрова-Гальоркіна, адаптивні вагові функції, домінуюча конвекція, рівняння Бюргерса.

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

Ключевые слова: метод Петрова-Галеркина, адаптивные весовые функции, доминирующая конвекция, уравнение Бюргерса.

Abstract. Modification of Petrov-Galerkin method for solving transient 2D problems with nonlinear convective terms of Burgers' type based on the weight functions changing in time depending on the evolutional solution is proposed. The test calculations confirm efficiency of the method proposed for convection-dominated problems.

Keywords: Petrov-Galerkin method, adaptive weight functions, convection-dominated problems, Burgers' equation.

1. Введение

Одной из особенностей уравнений задач гидро- и газодинамики, магнитной гидродинамики [1, 2] является наличие у них негладких (разрывных) решений, соответствующих ударным волнам [3]. Большинство подобных задач в силу их нелинейности даже в простейших случаях невозможно (или же нецелесообразно) решать аналитически, поэтому их решения получают численно.

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

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

© Молчанов А.А., Сирик С.В., Сальников Н.Н., 2012 ISSN 1028-9763. Математичні машини і системи, 2012, № 2

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

Искусственную диффузию целесообразно вводить лишь в окрестностях фронтов ударных волн, однако это дополнительно порождает непростую задачу отслеживания позиции фронта волны [3]. Кроме того, данный подход плохо себя зарекомендовал при решении многомерных задач [7]. Это связано с тем, что диффузию целесообразно вводить лишь вдоль направления потока, иначе пространственный вид решения может быть искажен и, следовательно, решение будет некорректным [7, 8]. Этим недостатком обладает большинство классических численных схем, использующих введение искусственной диссипации, так как при этом вносится диффузия поперек потока (crosswind diffusion).

Другим подходом к стабилизации является использование для производных аппроксимаций против потока (upwind/upstream difference) в конвективных слагаемых [5, 6]. Однако подобные методы имеют, как правило, лишь первый порядок аппроксимации и вносят в схему чрезмерную диссипацию [7, 8]. Поэтому данный подход к стабилизации также обладает перечисленными выше недостатками.

Одним из наиболее совершенных подходов к построению и стабилизации численных схем является подход Петрова-Г алеркина в форме метода конечных элементов (МКЭ) [6, 7]. При этом приближенное решение ищется в виде элемента некоторого конечномерного пространства (являющегося линейной оболочкой совокупности так называемых базисных функций с финитным носителем), а коэффициенты разложения определяются из условия ортогональности невязки другому конечномерному пространству (пространству, образованному совокупностью так называемых весовых функций). Для одной одномерной линейной стационарной задачи конвекции-диффузии были найдены такие весовые функции, при которых численное и аналитическое решения совпали в узлах сетки [9, 7]. Однако для нелинейных задач с преобладающей конвекцией численные решения, получаемые даже с помощью метода Петрова-Галеркина, могут содержать неприемлемые осцилляции нефизического характера [6, 10]. В работе [10] была предложена версия метода Петрова-Галеркина для интегрирования одномерного уравнения Бюргерса, основанная на использовании весовых функций, вид которых изменялся в зависимости от характера эволюционирующего решения. Это позволяло избежать появления неустойчивостей и осцилляций в решении, сохраняя при этом необходимую крутизну фронта в ударных волнах. В данной статье результаты работы [10] обобщены на случай нелинейного скалярного уравнения конвективно-диффузионного типа, зависящего от нескольких пространственных переменных.

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

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

du . du . du

-----+ liu----------+12u-----= к

Эt 9xi Эх2

(~\2 л2 ^

Э u + Э u Эxl2 Эx| j

u = u (t, x), x = (xb x2 ) gO , (1)

являющегося прямым обобщением одномерного уравнения Бюргерса [1, 6, 12, 13]. Здесь О - некоторая ограниченная область в двумерном вещественном евклидовом пространст-

2

ве R , ^1, 12 и к - константы, к > 0, время t е [0; T]. Предполагается, что для уравнения (1) задано начальное условие ^0, х) = Uo( х), а на границе ЭО области О выполняется одно из стандартных граничных условий [11]. Краткое описание физических задач, приводящих к уравнениям, подобного уравнению (1) вида, дано в [2, 3, 10, 12, 13].

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

3. Построение численного решения

Предположим, что область О (черта над множеством обозначает его замыкание в ^ ) триангулирована [14], т.е. разбита на конечную совокупность треугольных элементов (треугольников) О k : О = и Оk , I іп Оі = 0 (ІП Оk - множество внутренних точек Оk ).

k k

Будем предполагать, что границей пересечения двух соседних треугольников может быть только смежное ребро этих треугольников целиком. Каждый элемент однозначно определяется своими вершинами (узлами сетки). Триангуляцию области можно осуществить, используя один из алгоритмов [15]. Считаем, что в результате триангуляции области получается N узлов, которые будем нумеровать натуральными числами от 1 до N .

Будем искать приближенное слабое решение [5, 16] уравнения (1) в виде

N

u

х) = Ъ aj 0) Nj (x),

7=1

(2)

где (Му (х)} - набор заданных непрерывных кусочно-линейных базисных функций,

(ау (^)} - подлежащие определению неизвестные коэффициенты. Каждая функция Му (х), у = 1 к N отлична от нуля только на элементах, имеющих узел с номером у своей вершиной, равна единице в узле у (т.е. Му (Ху) = 1, где Ху = (Ху1, Ху2) е О) и равна нулю в остальных узлах (рис. 1). В этом смысле функция Му и узловая точка Ху находятся во взаимно однозначном соответствии.

В соответствии с процедурой метода Галеркина [5, 6, 16] умножим равенство (1) на весовую функцию (х), соответствующую узлу сетки г, и проинтегрируем по области

О. Описание функции (х) будет приведено ниже. Подставив выражение (2) в получившееся тождество, получим для определения коэффициентов (ау (^)} следующую систему обыкновенных дифференциальных уравнений (г = 1... N):

N ( ^

ёа7

Ъ

7 =14 О

N N( „

-І2ЪЪ \\wNj —со

9x2

л

N N (гг ЭМ ^

-І1ЪЪ Э^со

7=1 і-=1^ О 7 Эх1

л

аа-к

а.аі -

к

.=1 і=14 о

N ( ЭЖг ^7

Ъ О ах^ аТ® а + к /' ■

7=14 О 2 2 у

Эх1 Эх1

а7 -

(3)

где является слагаемым, полученным после применения ко вторым производным многомерной формулы интегрирования по частям (или формул Грина для оператора Лапласа)

В формуле (4) криволинейный интеграл берется вдоль контура ЭО (границы области О), который предполагается достаточно гладким [11], п = («1,«2) - вектор единичной внешней нормали к области О .

Рассмотрим вопрос учета граничных условий в системе (3). Пусть граница ЭО

представляется в виде ЭО = ЭОи и ЭОе ( ЭОи I ЭО е =0 ), и на ЭОи задано условие типа Неймана [11]:

При получении равенства (7) предполагалось, что функция Щ (х) обращается в нуль на множестве ЭО е. Для аппроксимации граничного условия (6) можно в пробном решении (2) положить а у = g (ху), если узел ху е ЭО е (при этом соответствующая часть коэффициентов (а у (^)} становится известной и потому не входит в систему (3) в качестве неизвестных переменных [10]). Для аппроксимации начального условия и(0, х) = и0(х) можно положить а у (0) = и о (ху ) или же использовать среднеквадратическую аппроксимацию функции и0(х) по некоторому набору функций (например, по набору базисных функций Му, как осуществляется в [16]).

В качестве весовых функций Щ используем функции, предложенные в [17]. Обозначим через О(г) многоугольник, образованный объединением тех элементов Ок, для которых узел г является одной из вершин. Множество О(г) является носителем для функции N (х), а также для функции Щ (х) , которая определяется следующим образом [17]:

где Кг - множество номеров узлов к, каждый из которых соединен одним ребром с узлом г, аг к - настроечный числовой параметр (параметр изгиба), соответствующий ребру (г, к) . Этот параметр позволяет менять форму функции Щ (х) . Формулы для его вычисления приведены ниже. Носителем функции Щ(г к) (х) является объединение элементов О у и О у', для которых ребро (г, к) является общим, она непрерывна на О и квадратичная на

[5, 6]:

(4)

п • Уи = И на ЭО к,

(5)

а на ЭО е - условие типа Дирихле [11]:

и = g на ЭО е .

(6)

Тогда в формуле (4)

(7)

ЭОи

Ж; (х) = N (х) + Ъ а )(х),

(8)

іеКі

каждом из элементов Оу и О у'. В узлах г и к функция Щ к) (х) равна нулю, в середине ребра (г, к) принимает значение, равное (-0,75). Способ построения функции Щ(г к)(х)

подробно изложен в [17].

Заметим, что при выборе весовой функции Щ (х) в виде (8) интегралы от произведений весовых и базисных функций (и их производных), входящие в систему (3), линейно зависят от параметров а,к, и потому непосредственно операции интегрирования достаточно провести лишь один раз в начале расчета. Например,

\\wiNуdО = Ц NiNуdО+ ^ аг,к {{Щлк)NуdО . Аналогично и для других интегралов, О О кеКг о

входящих в систему (3).

Рассмотрим вопрос выбора настроечных параметров (аг к } в формуле (8), определяющих величину и направление изгиба весовой функции Щ (х). Заметим, что выбор параметров изгиба является, пожалуй, наиболее важным этапом при построении схемы интегрирования, так как путем надлежащего выбора параметров можно добиться отсутствия неустойчивостей и нефизических осцилляций (при сравнительно небольшом количестве узлов сетки) [18, 19]. В случае, когда все параметры а,к тождественно равны нулю, из (8)

получаем Щ = Ni, что соответствует стандартному (классическому) методу Галеркина, который в задачах с доминирующей конвекцией может давать неустойчивые осциллирующие результаты. Заметим, что в задачах с ненулевой диффузией в ряде случаев имеется условная устойчивость метода Галеркина, т.е. можно добиться приемлемой устойчивости путем значительного увеличения количества базисных функций, но при этом резко возрастает вычислительная сложность задачи. Анализ устойчивости различных вариантов метода Галеркина для линейных одномерных задач приведен в [6].

Напомним, что метод Петрова-Галеркина отличается от метода Галеркина использованием набора весовых функций (Щ(х)}, не совпадающего с набором базисных {Ni(х)} .

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

Можно показать [7-9], что при использовании набора одномерных кусочноквадратичных весовых функций (подобных функциям (8) и также зависящих от совокупности параметров изгиба (а г}) в методе Петрова-Г алеркина для интегрирования одномер, , Л du d2и ,

ного линейного стационарного уравнения конвекции-диффузии 1 — = V —— (при посто-

dх dх2■

янных 1 и V ) зависимость

аг = а(Ре) = соШ(Ре • И /2) - 2 /(Ре • И)

параметров изгиба от числа Ре /V [5, 6] обеспечивает совпадение численного решения с точным решением указанного уравнения в узлах сетки. Особенности построения различных одномерных весовых функций детально рассмотрены в [10, 17].

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

сами должны зависеть от времени (в противоположном случае решения будут получаться неустойчивыми, осциллирующими и с большими погрешностями [6]). В статье [10] подробно проанализированы недостатки стандартных методов типа Галеркина (и Петрова-Г алеркина [6, 7]) при интегрировании одномерного нестационарного нелинейного уравне-

Эи . Эи Э 2и

ния-----+ lu— = к—— и отмечена, в частности, их связь с отсутствием учета характера

dt Эх Эх2

эволюционирующего во времени решения. В связи с этим в работе [10] параметр изгиба ai, соответствующий одномерной весовой функции, на n + 1-м шаге интегрирования получившейся после применения процедуры Галеркина системы обыкновенных дифференциальных уравнений (см. выражение (3)) предложено выбирать в зависимости от решения соответствующей системы на предыдущем n -м временном шаге.

Но в многомерном случае, в отличие от описанного одномерного, весовая функция (8) характеризуется уже не одним параметром ai, а целым набором {ai k }, где k е Kt. В

связи с этим в данной работе предлагается несколько вариантов выбора параметров изгиба ai к (при к е Kt) на n + 1-м шаге интегрирования системы (3).

Вариант 1:

a(”k+1) = 0 ' a(gi,k ). (9)

Вариант 2:

a (J+1> =0-m-Fp (a(ya.)), (10)

[Pz, | bz | < 1

где функция FP (z) = < , ар - некоторое фиксированное неотрицательное

F [sign z, | pz | > 1

число. В формулах (9)-(10) функция a(g) = coth(g) -1/ g, gik = (Pe-h^ )/2, где (Pei -hik ) - евклидово скалярное произведение вектора hi к ° Хк - xi, проведенного из узла xt в xk, и вектора локального числа Пекле Pet ° (11u(n)(xt); 12u(n)(xt))/ к, характери-зирующего направление и величину потока в узле xi. В последней формуле и(n)(x) ° ~(x, tn ) - решение, определяемое выражением (2) на предыдущем n -м временном шаге интегрирования системы (3). Число m = max | a(gi k ) |. Настроечный параметр 0

i,k

предоставляет дополнительные возможности в регулировании величины ai k . Обычно в приведенных ниже расчетах его значения выбирались из промежутка 0,5 < 0 < 1. При

0 = 0 получаем метод Галеркина; увеличение параметра 0 приводит к возрастанию абсо-

лютной величины, вносимой в численную схему искусственной диссипации.

Использование скалярного произведения при определении выражений gik приводит к пропорциональной зависимости величины a(gi k) от угла между векторами Pei и

ht k . Когда угол между векторами Pei и hi k является тупым, их скалярное произведение отрицательно, и в случае их противоположной направленности максимально по абсолютной величине. При этом величина a(gi k ) (и соответственно параметр a ) также будет отрицательной, а весовая функция (8) на ребре (i, k) положительной и выпуклой вверх, что обеспечивает взвешивание выражений в (3) компонентой W(i k)(x) весовой функции в пределах своего носителя с большим весом со стороны набегающего потока (т.е. в направле-

нии, противоположном вектору Ре1 ). Аналогично, если направления векторов Ре1 и Ъ1 к

образуют острый угол, а(уг-,£) и а(?к+1) будут положительными, а Щг,к)(х) на ребре (г, к)

выпуклой вниз, и тогда при взвешивании невязки в (3) влияние выражений системы вдоль потока будет уменьшено (по сравнению с направлениями, образующими с направлением

потока тупой угол). Когда же векторы Рег и к перпендикулярны, то уг- к = 0, и тогда

а(?к+1^ = ®(уг к ) = 0 (следовательно, при взвешивании невязки коэффициент а(?к+1) не будет вносить вклад в создание величины диффузии поперек потока).

Параметр Ь в выражении (10) позволяет регулировать ширину "коридора", в котором величина а(>к+1^ линейно зависит от а(уг к). Чем больше значение параметра Ь, тем меньшей является эта ширина. В частном случае, положив формально Ь =+¥ в (10), видим, что выражение для коэффициента изгиба аг,к превращается в выражение

а(”к+1) = 0 •т • §1§п(а(У г,к )). (11)

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

4. Численное моделирование

Рассмотрим следующую краевую задачу для уравнения (1): пусть область О ° (0; 1)х(0; 1) - единичный квадрат, начальное условие ^ (хь х2 ) = 1 - хь константы 1 = 1, 12 = 2,

к = 10-6, граничное условие является условием типа Дирихле и получается из начального условия путем его непрерывного продолжения на границу области О. Поскольку диффузионный параметр к мал (в сравнении с конвективными параметрами 11, 12 и величиной начального условия ^( хь Х2) ), процесс конвекции в задаче будет доминирующим. На физическом уровне данная краевая задача описывает процесс нелинейного переноса начального условия в направлении, задаваемом локально в каждой точке х е О вектором V = (1^(х); 12и(х)) [1, 12, 13]. Поскольку 12 > 1, то в каждой точке области (локальная) скорость переноса в направлении возрастания координаты х2 будет больше скорости переноса в направлении возрастания х1 , и потому, в силу необходимости удовлетворения граничных условий, возле участка х2 = 1 границы образуется резкий перепад решения (приграничный слой) [12, 13, 10]. График решения задачи методом Галеркина в момент времени t = 0,4 приведен на рис. 2 (использовано разбиение области 30 х 30), график решения методом Петрова-Г алеркина с весовыми функциями (8), где параметры изгиба определяются с помощью выражений (11), представлен на рис. 3 (использовано разбиение области 20х20, параметр 0° 1). Для случая, когда параметры изгиба определяются выражениями (9), график численного решения практически совпадает с графиком на рис. 3 и поэтому здесь не приводится. Таким образом, в данном примере предложенные варианты выбора весовых коэффициентов показывают практически одинаковую эффективность. Заметим, однако, что в ряде случаев вариант 1 выбора аг,к (в виде выражения (9)) при ин-

тегрировании системы (3) (и при условии, что не применяются никакие дополнительные средства стабилизации численных схем [6, 10]) обеспечивает более точные и надежные результаты расчетов.

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

Рис. 2. Решение методом Галеркина в Рис. 3. Решение предложенным мето-

момент времени t = 0,4 при разбие- дом Петрова-Галеркина в момент вре-нии области на 30х30 узлов мени t = 0,4 при разбиении обдасти на

20х20 узлов

Отметим, что для реальных нелинейных задач, описываемых уравнениями, подобными уравнению (1), в большинстве случаев весьма проблематично получить аналитическое решение. Но даже в тех редких случаях, когда решение задачи все же удается получить в явном аналитическом виде, оно может выражаться через медленно сходящиеся интегралы/ряды или специальные функции, суммирование и вычисление которых само по себе является весьма сложной проблемой. Для ее преодоления часто используется замена исходных выражений их главными асимптотиками при стремлении переменной интегрирования (параметра суммирования) к бесконечности, что в конечном итоге негативно отражается на точности получаемого результата (характер возникающих проблем достаточно наглядно продемонстрирован в работе [20], а также [10]). Поэтому для подобного рода задач применение численных методов и использование физических аналогий для контроля результатов являются практически единственными доступными средствами проверки правильности получаемого решения.

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

Предложено обобщение метода численного решения краевых задач для уравнений с нелинейными конвективными членами типа уравнения (1), изложенного в работе [10] и базирующегося на проекционном подходе Петрова-Галеркина. Построена конечномерная модель в виде системы (3) обыкновенных дифференциальных уравнений. Предложены способы выбора коэффициентов, определяющих изгиб весовых функций, что позволяет избежать возникновения неустойчивостей и нефизических осцилляций при решении задач с доминированием процессов конвекции. Результаты, полученные в работе, подтверждаются расчетными данными.

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

СПИСОК ЛИТЕРАТУРЫ

1. Ладиков-Роев Ю.П. Математические модели сплошных сред / Ю.П. Ладиков-Роев, О.К. Черем-ных. - К.: Наукова думка, 2010. - 551 с.

2. Куликовский А.Г. Магнитная гидродинамика / А.Г. Куликовский, Г.А. Любимов. - [изд. 2-е,

испр. и доп.]. - М.: Логос, 2005. - 328 с.

3. Куликовский А.Г. Математические вопросы численного решения гиперболических систем уравнений / Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. - М.: Физматлит, 2001. - б08 с.

4. Roos H.-G. Robust numerical methods for singularly perturbed differential equations / Roos H.-G., Stynes M., Tobiska L. - Berlin, Heidelberg: Springer-Verlag, 2008. - б04 p.

5. Флетчер К. Численные методы на основе метода Галеркина / Флетчер К.; пер. с англ. - М.: Мир, 1988. - 352 с.

6. Finlayson B.A. Numerical methods for problems with moving fronts / Finlayson B.A. - Seattle, Washington USA: Ravenna Park Publishing, Inc., 1992. - бІЗ p.

7. Fries T.P. A Review of Petrov-Galerkin Stabilization Approaches and an Extension to Meshfree Methods / T.P. Fries, H.G. Matthies. - Germany; Brunswick: Technische Universitat Braunschweig, Informa-tikbericht-Nr., 2004. - 71 p.

8. Brooks A.N. Stremline upwind Petrov-Galerkin formulations for convection dominated flows with particular emphasis on incompressible Navier-Stokes equations / A.N. Brooks, T.J.R. Hughes // Computer Methods in Applied Mechanics and Engineering. - 1982. - N 32. - P. 199 - 259.

9. Finite element methods for second order differential equations with significant first derivatives / I. Christie, D.F. Griffiths, A.R. Mitchell [et al.] // International Journal of Numerical Methods in Engineering. - 197б. - Vol. 10. - P. 1389 - 139б.

10. Сирик С.В. Численное интегрирование уравнения Бюргерса методом Петрова-Галеркина с адаптивными весовыми функциями / С.В. Сирик, Н.Н. Сальников // Проблемы управления и информатики. - 2012. - № 1. - C. 94 - 110.

11. Тихонов А.Н. Уравнения математической физики / А.Н. Тихонов, А.А. Самарский. - М.: Наука, 1999. - 798 с.

12. Заславский Г.М. Введение в нелинейную физику: От маятника до турбулентности и хаоса / Г.М. Заславский, Р.З. Сагдеев. - М.: Наука. Гл. ред. физ.-мат. лит., 1988. - 308 с.

13. Солитоны и нелинейные волновые уравнения / Р. Додд, Дж. Эйлбек, Дж. Гиббон [и др.]; пер. с англ. - М.: Мир, 1988. - б9б с.

14. Сегерлинд Л. Применение метода конечных элементов / Сегерлинд Л. - М.: Мир, 1979. - 392 с.

15. Препарата Ф. Вычислительная геометрия: Введение / Ф. Препарата, М. Шеймос; пер. с англ. -М.: Мир, 1989. - 478 с.

їб. Дейнека В.С. Математические модели и методы расчета задач с разрывными решениями / Дейнека В.С., Сергиенко И.В., Скопецкий В.В. - К.: Наукова думка, 1995. - 2б2 с.

17. Сальников Н.Н. О построении конечномерной математической модели процесса конвекции-диффузии с использованием метода Петрова-Галеркина / Н.Н. Сальников, С.В. Сирик, И.А. Терещенко // Проблемы управления и информатики. - 2010. - № 3. - С. 94 - 109.

18. Harari I. Streamline design of stability parameters for advection-diffusion problems / I. Harari, L.P. Franca, S.P. Oliveira // Journal of Computational Physics. - 2001. - N 171 (1). - P. 115 - 131.

19. Tezduyar T.E. Stabilization and shock-capturing parameters in SUPG formulation of compressible flows / T.E. Tezduyar, M. Senga // Comput. Methods Appl. Mech. Engrg. - 200б. - N 195. - P. 1б21 -1бЗ2.

20. Hon Y.C. An efficient numerical scheme for Burgers' equation / Y.C. Hon, X.Z. Mao // Applied Mathematics and Computation. - 1998. - Vol. 95, N 1. - P. 37 - 50.

Стаття надійшла до редакції 05.04.2012