О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии

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

Full description

Saved in:

Bibliographic Details
Date:	2014
Main Author:	Сирик, С.В.
Format:	Article
Language:	Russian
Published:	Інститут проблем математичних машин і систем НАН України 2014
Series:	Математичні машини і системи
Subjects:	Моделювання і управління
Online Access:	http://dspace.nbuv.gov.ua/handle/123456789/84458
Tags:	Add Tag No Tags, Be the first to tag this record!
Journal Title:	Digital Library of Periodicals of National Academy of Sciences of Ukraine
Cite this:	О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии / С.В. Сирик // Математичні машини і системи. — 2014. — № 4. — 139-155. — Бібліогр.: 30 назв. — рос.

Institution

Digital Library of Periodicals of National Academy of Sciences of Ukraine

id	irk-123456789-84458
record_format	dspace
spelling	irk-123456789-844582015-07-09T03:02:10Z О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии Сирик, С.В. Моделювання і управління Рассмотрены вопросы выбора и динамической (адаптивной) настройки стабилизирующих параметров весовых функций в конечноэлементном методе Петрова-Галеркина при интегрировании уравнений конвективно-диффузионного типа. На ряде вычислительных примеров продемонстрирована эффективность соответствующего метода Петрова-Галеркина. Розглянуто питання вибору та динамічного (адаптивного) налаштування стабілізуючих параметрів вагових функцій у скінченноелементному методі Петрова-Гальоркіна при інтегруванні рівнянь конвективно-дифузійного типу. На ряді обчислювальних прикладів продемонстровано ефективність відповідного методу Петрова-Гальоркіна. The choice and dynamical tuning of stabilization parameters of the weighting functions in finite-element Petrov-Galerkin method for solving convection-diffusion problems are considered. The efficiency of the corresponding Petrov-Galerkin method is illustrated and confirmed by a number of computational examples. 2014 Article О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии / С.В. Сирик // Математичні машини і системи. — 2014. — № 4. — 139-155. — Бібліогр.: 30 назв. — рос. 1028-9763 http://dspace.nbuv.gov.ua/handle/123456789/84458 519.63; 536.252; 004.75 ru Математичні машини і системи Інститут проблем математичних машин і систем НАН України
institution	Digital Library of Periodicals of National Academy of Sciences of Ukraine
collection	DSpace DC
language	Russian
topic	Моделювання і управління Моделювання і управління
spellingShingle	Моделювання і управління Моделювання і управління Сирик, С.В. О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии Математичні машини і системи
description	Рассмотрены вопросы выбора и динамической (адаптивной) настройки стабилизирующих параметров весовых функций в конечноэлементном методе Петрова-Галеркина при интегрировании уравнений конвективно-диффузионного типа. На ряде вычислительных примеров продемонстрирована эффективность соответствующего метода Петрова-Галеркина.
format	Article
author	Сирик, С.В.
author_facet	Сирик, С.В.
author_sort	Сирик, С.В.
title	О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии
title_short	О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии
title_full	О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии
title_fullStr	О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии
title_full_unstemmed	О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии
title_sort	о выборе стабилизирующих параметров в конечноэлементном методе петрова-галеркина при решении задач конвекции-диффузии
publisher	Інститут проблем математичних машин і систем НАН України
publishDate	2014
topic_facet	Моделювання і управління
url	http://dspace.nbuv.gov.ua/handle/123456789/84458
citation_txt	О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии / С.В. Сирик // Математичні машини і системи. — 2014. — № 4. — 139-155. — Бібліогр.: 30 назв. — рос.
series	Математичні машини і системи
work_keys_str_mv	AT siriksv ovyborestabiliziruûŝihparametrovvkonečnoélementnommetodepetrovagalerkinaprirešeniizadačkonvekciidiffuzii
first_indexed	2025-07-06T11:27:05Z
last_indexed	2025-07-06T11:27:05Z
_version_	1836896722514608128
fulltext	© Сирик С.В., 2014 139 ISSN 1028-9763. Математичні машини і системи, 2014, № 4 УДК 519.63; 536.252; 004.75 С.В. СИРИК* О ВЫБОРЕ СТАБИЛИЗИРУЮЩИХ ПАРАМЕТРОВ В КОНЕЧНОЭЛЕМЕНТНОМ МЕТОДЕ ПЕТРОВА-ГАЛЕРКИНА ПРИ РЕШЕНИИ ЗАДАЧ КОНВЕКЦИИ-ДИФФУЗИИ * Национальный технический университет Украины "КПИ", Киев, Украина Анотація. Розглянуто питання вибору та динамічного (адаптивного) налаштування стабілізую- чих параметрів вагових функцій у скінченноелементному методі Петрова-Гальоркіна при інтег- руванні рівнянь конвективно-дифузійного типу. На ряді обчислювальних прикладів продемонстро- вано ефективність відповідного методу Петрова-Гальоркіна. Ключові слова: метод Петрова-Гальоркіна, метод скінченних елементів, задачі конвекції-дифузії, вагові функції, стабілізуючі параметри, стабілізовані методи. Аннотация. Рассмотрены вопросы выбора и динамической (адаптивной) настройки стабилизи- рующих параметров весовых функций в конечноэлементном методе Петрова-Галеркина при ин- тегрировании уравнений конвективно-диффузионного типа. На ряде вычислительных примеров продемонстрирована эффективность соответствующего метода Петрова-Галеркина. Ключевые слова: метод Петрова-Галеркина, метод конечных элементов, задачи конвекции- диффузии, весовые функции, стабилизирующие параметры, стабилизированные методы. Abstract. The choice and dynamical tuning of stabilization parameters of the weighting functions in finite- element Petrov-Galerkin method for solving convection-diffusion problems are considered. The efficiency of the corresponding Petrov-Galerkin method is illustrated and confirmed by a number of computational examples. Keywords: Petrov-Galerkin method, finite element method, convection-diffusion problems, weighting (test) functions, stabilization parameters, stabilized methods. 1. Введение В настоящее время метод Петрова-Галеркина (МПГ) в форме метода конечных элементов (МКЭ) является одним из наиболее универсальных методов построения численных схем для решения задач математической физики [1]. При использовании МПГ для расчета задач с доминированием конвективных процессов ключевым вопросом является корректный вы- бор весовых функций (ВФ), который предотвращал бы появление в численном решении нефизических осцилляций и неустойчивостей (обеспечил стабилизацию [1] численного решения) при сохранении приемлемой его точности. В работе [2] предложены структуры кусочно-полиномиальных ВФ МПГ для решения двухмерных задач конвекции-диффузии (КД), а в работе [3] – для интегрирования задач с тремя пространственными переменными. Для настройки формы функций используются стабилизирующие параметры (СП), связан- ные с ребрами сетки разбиения. Анализ аппроксимации, устойчивости и сходимости для данных ВФ [2, 3] (а также их кусочно-степенных обобщений) для нестационарных одно- мерных случаев проведен в работе [4]. Данные ВФ впоследствии применялись для числен- ного решения различных нестационарных задач КД (в том числе и в тех случаях, когда скорость в конвективном слагаемом с течением времени резко изменяется как по величи- не, так и по направлению), а также нелинейных уравнений Бюргерса [5, 6] и уравнений магнитной гидродинамики [7]; в цитированных работах также были предложены методы (способы) выбора СП, реализующие идею адаптивности (возможности подстраивать их под эволюционирующее во времени решение для обеспечения устойчивости счета). В настоящей статье результаты работ [5–7], касающиеся выбора СП, получили дальнейшее развитие, и на примерах (считающихся "тяжелыми" для численного решения 140 ISSN 1028-9763. Математичні машини і системи, 2014, № 4 по причине наличия в них доминирующих конвективных процессов, пограничных или внутренних слоев, решений солитонного типа и прочих особенностей, доставляющих ос- ложнения при численном счете) показано эффективность ("робастность" [1]) предложен- ных вариантов выбора СП в плане точности и устойчивости численного решения. 2. Постановка задачи Будем считать, что в качестве базисных функций )}({ xNi � (каждая функция )(xNi � ассо- циирована с соответствующим узлом ix � ) используются кусочно-линейные финитные функции [1–7]. В качестве ВФ )(xWi � используем функции работ [2–4]. Обозначим через )(iΩ многоугольник, образованный объединением тех элементов kΩ , для которых узел i является одной из вершин. Множество )(iΩ является носителем для функции iN , а также для функции iW , которая определяется как [2–4]: )()()( ),(, xWxNxW iKk kikiii �� ∑ ∈ α+= , (1) где iK – множество номеров узлов k , каждый из которых соединен одним ребром с узлом i , ki,α – СП, соответствующий ребру ),( ki и позволяющий регулировать изгиб ВФ на нем. Подробная конструкция функций ),( kiW , а также всей ВФ iW , приведена в [2, 3]. Целью данной работы является выбор совокупности СП }{ ,kiα и их динамическая настройка, при котором бы обеспечивалась стабилизация (устранение нефизических осцилляций) числен- ного решения при сохранении приемлемой точности получаемого решения. 3. Динамическая настройка стабилизирующих параметров При интегрировании полудискретных аппроксимаций [1–7] в виде системы обыкновенных дифференциальных уравнений (СОДУ), построенной МПГ для аппроксимации уравнений КД, на 1+n -м шаге интегрирования предлагается выбирать СП }{ ,kiα ВФ )(xWi � в виде ( ))( )( , )1( , )1( , )1( n ki n ki n ki nF γαθ=α +β ++ , (2) где )1( +β n – некоторое фиксированное неотрицательное число (зависящее от номера шага), функция    >β ≤ββ ≡β 1\|\|,sign 1\|\|, )( zz zz zF , функция ( ) coth( ) 1/α γ ≡ γ − γ (при 0=γ полагаем, что 0 (0) lim(coth( ) 1/ ) 0 γ→ α ≡ γ − γ = ), число )~2/()( )1( , )()( , +κ⋅=γ n iki nn ki hu � � , ,i k k ih x x≡ − � � � , ( )nu � – вектор скорости (векторная величина, которая умножается скалярно на градиент неизвестной ис- комой величины в рассматриваемом уравнении типа КД) в i -м узле на предыдущем n -м шаге интегрирования СОДУ (в качестве )(nu � также можно взять взвешенное среднее дан- ных векторов в узлах i и k [7]), )1(~ +κ n i – коэффициент, регулирующий величину вводимой в расчетную схему искусственной вязкости относительно i -го узла на 1+n -м шаге интег- рирования (соответственно, в силу определения ( )α γ , чем меньше значение )1(~ +κ n i , тем большей будет искусственная вязкость). При изотропном (диагональном) диффузионном тензоре, в большинстве случаев, как показывают расчеты для устранения нефизических осцилляций, выбор данного коэффициента вязкости в виде κ (где κ – диффузионный ко- ISSN 1028-9763. Математичні машини і системи, 2014, № 4 141 эффициент при диффузионных членах решаемого уравнения КД) оказывается вполне дос- таточным. Параметр )1( , +θ n ki – некоторая константа, увеличение или уменьшение которой приводит соответственно к увеличению или уменьшению вводимой в численную схему искусственной диссипации. В большинстве случаев для нужд практических расчетов вполне достаточно выбирать ее из промежутка ]2;0[ . Увеличение этой константы может понадобиться при возникновении в решении больших градиентов и нефизических осцил- ляций в соответствующих тонких слоях (поскольку увеличение искусственной диссипации приводит к "монотонизации" численных схем и сглаживанию неоднородностей, разрывов и скачков). Заметим, что при 0)1( , =θ +n ki по определению получаем классический МКЭ Га- леркина. Положив в (2) +∞=β + )1(n и 1)1( =β +n , получим соответственно ( ) ( )( 1) ( 1) ( ) ( 1) ( ) , , , , ,sign ( ) signn n n n n i k i k i k i k i k + + +α = θ α γ = θ γ , (3) )( )( , )1( , )1( , n ki n ki n ki γαθ=α ++ . (4) Отметим, что выражения (2)-(4) для СП пригодны не только при интегрировании полудис- кретных аппроксимаций, но без изменений могут быть использованы при нахождении ре- шений разностных схем (при переходе от предыдущего n -го временнóго слоя схемы на текущий 1+n -й). Естественно, при интегрировании систем уравнений настройка СП в ВФ для каждого уравнения системы осуществляется индивидуально. Для определения СП при решении стационарных задач используются те же формулы (2)-(4) (однако, теперь они ап- риори являются независимыми от n ). При 1)1( , =θ +n ki формула (4) дает выражение, которое для стационарного линейного одномерного уравнения КД с постоянными коэффициентами обеспечивает совпадение численного и точного решений в узлах равномерной сетки (то есть, получаем схему Ильина-Аллена-Саусвелла [1]). Дополнительные замечания, касаю- щиеся выбора СП, упрощения их вычисления, а также стабилизации решений в окрестно- стях фронтов ударных волн приведены в [5] (см. замечания 1-5 в разделе 2 работы [5]). 4. Численное моделирование Пример 1. Рассмотрим граничную задачу для стационарного двухмерного уравнения КД uuv ∆κ=∇⋅� (тут v � – вектор скорости, κ – диффузионный коэффициент) в единичном квадрате Ω (центр квадрата находится в начале координат). Для упрощения записи будем обозначать 1xx ≡ , 2xy ≡ . Поле скоростей имеет вид Txyyxv );(),( −=� (вращение против часовой стрелки), 810−=κ . На границе квадрата Ω решение u равно нулю, а внутри квад- рата, на прямой 0=x , 02/1 ≤≤− y имеется разрез, на котором )(),0( ywyuu == , где 2/)1)4(cos()( +π+π= yyw . Постановка граничных условий и линия разреза изображены на рис. 1 [8]. На рис. 2 показан график функции )(yw при 02/1 ≤≤− y . На рис. 3 показано решение данной задачи [8] – "воронка", которая образуется при вращении графика функ- ции )(yw против часовой стрелки относительно начала координат. Данная задача исполь- зовалась в [8] для тестирования точности и работоспособности стабилизированных мето- дов SUPG, GLS, USFEM и RFB (однако, в [8] было взято 610−=κ ). Поскольку процесс конвекции в задаче является доминирующим над процессом диффузии ( κ – малое), то при вращении графика )(yw его высота должна оставаться практически такой же, какой она была на разрезе (то есть график не "проседает" под действием диффузии, распределенной в данном случае изотропно) [8]. Это, в свою очередь, для оценки точности и качества чис- 142 ISSN 1028-9763. Математичні машини і системи, 2014, № 4 ленного метода позволяет использовать максимальную разницу между высотой )(yw на линии разреза и высотой графика вычисленного решения на вертикальной линии узлов, находящихся слева от линии разреза на ближайшем расстоянии (для расчетов будем ис- пользовать сетку с прямоугольными треугольниками, ориентированными вдоль коорди- натных осей). Обозначим данную разницу в текущей задаче через maxerr . Отметим, что в силу особенностей данной задачи ее можно рассматривать как "стресс-тест" для проверки того, вводит ли численный метод избыточную диффузию поперек потока [1, 6, 8]. В случае утвердительного ответа при полном повороте вокруг начала координат график под воздей- ствием диффузии будет "разглажен", а величина maxerr , соответственно, будет большой [8]. Таким образом, данная задача предъявляет высокие требования к численным методам (в особенности, к их диссипативности) и считается весьма тяжелой для численного реше- ния [8]. Рис. 1. Постановка граничных Рис. 2. График )(yww = . Рис. 3. Пространственный условий (пример 1) профиль решения (пример 1) Для метода SUPG используем поэлементную стабилизацию [1, 3] с СП }{ τα вида ( ) / (2 \|\| \|\|)h Pe vτ τ τα = α � , (5) где τh – диаметр элемента τT , )2/(\|\|\|\| κ= ττ hvPe � , \|\|\|\|v � – евклидова норма вектора v � , вы- численного в центре элемента τT [8] по поводу данного выбора параметров. Данный выбор параметров в SUPG является оптимальным для одномерных стационарных задач КД с по- стоянными коэффициентами (дает решение, совпадающее с точным решением в узлах) и считается также вполне подходящим для большинства задач (в том числе, многомерных) [1, 8, 9] (особенно таких, в которых процесс конвекции не сильно доминирует над осталь- ными процессами). Используем сначала равномерную сетку с треугольными элементами и числом узлов 2121× . Для МПГ с ВФ (1) (используем здесь, для определенности, кусочно- квадратичные ВФ) и при выборе параметров в виде (4) для всех 10/1, =θ ki получаем 097,0max ≈err , а график решения практически совпадает с истинным решением (рис. 3); при возрастании или убывании ki,θ от данного значения погрешность растет. Так, для 0, =θ ki (МКЭ Галеркина) получаем 145,0max ≈err (а в численном решении имеются ос- цилляции, что противоречит физическому смыслу задачи), при 2/1, =θ ki будет 109,0max ≈err , а при 1, =θ ki получим 119,0max ≈err . Метод SUPG с выбором }{ τα в виде (5) дает 129,0max ≈err . В работе [8] для решения данной задачи также применялся метод SUPG с выбором }{ τα в виде )( )sincos31( sincos \|\|\|\|2 τ τ τ α φφ+ φ+φ=α Pe v h � , где φ – угол (локально- ISSN 1028-9763. Математичні машини і системи, 2014, № 4 143 го) направления потока (относительно декартовой системы координат, см. подробно в [8]). В работе [8] показано, что для стационарных двухмерных задач такой выбор стабилизи- рующих параметров в методе SUPG обеспечивает наилучшую точность в случае домини- рования конвекции (по сравнению не только с другими вариантами выбора параметров, но, как показали расчеты [8], также и с другими стабилизированными методами, в особенно- сти, методами GLS, USFEM и RFB); метод SUPG с таким выбором }{ τα дает 117,0max ≈err . Рассмотрим теперь результаты расчетов для равномерной сетки 3131× узлов. Сно- ва же, при 10/1, =θ ki получаем минимальную погрешность 043,0max ≈err . При 0, =θ ki получим 072,0max ≈err , при 2/1, =θ ki — 046,0max ≈err , при 1, =θ ki – 048,0max ≈err . Для метода SUPG с (5) получим 05,0max ≈err , а для вышеописанного оптимального выбора па- раметров }{ τα – 046,0max ≈err . Для МПГ с ВФ (1), параметрами (4) и 10/1, =θ ki , а также для метода SUPG c вышеописанным оптимальным выбором параметров при использова- нии разбиения 2121× узлов число обусловленности (в спектральной норме [10]) 126≈ (для разбиения 3131× узлов приближенно равно 206 ); для метода SUPG c параметрами (5) и использовании разбиения 2121× число обусловленности 210≈ (для разбиения 3131× узлов – приближенно равно 352 ). Таким образом, расчеты свидетельствуют, что предлагаемая в данной работе версия МПГ способна в данной задаче обеспечить бóльшую точность решения по сравнению с другими стабилизированными методами и, кроме того, по сравнению с методом SUPG, формирует систему уравнений с меньшим числом обу- словленности и не вводит избыточную диффузию поперек потока (что следует из опреде- ления метода SUPG и рассчитанных значений maxerr ). Пример 2. Рассмотрим задачу для нестационарного двухмерного уравнения КД:       ∂ ∂κ ∂ ∂+      ∂ ∂κ ∂ ∂= ∂ ∂+ ∂ ∂+ ∂ ∂ 22112 2 1 1 x u xx u xx u v x u v t u (6) в прямоугольной области, ii Lx ≤≤0 ( 2,1=i ) при постоянном коэффициенте κ и векторе скорости const);( 21 == Tvvv � . Начальное условие задается в виде 1),,0( 21 xexxu = . Гранич- ные условия являются продолжением по непрерывности начальных условий на границу области: 1),0,( 2 =xtu , 1),,( 21 LexLtu = , 1),,()0,,( 211 xeLxtuxtu == . Аналитическое решение данной задачи для (6) выражается в виде [2] 121 1 1 2 2 1 1 21 sinsin),,( x i j tat ij ijcxbxat e L jx L ix ee a H exxtu ij +ππ      − −Ω = ∑∑ ∞ = ∞ = Ω−−++ , где ( ) 2222 2 1 2222 1 1)1( 222 )1(1 )1( )1(1 )1(4 21 jcL e ibL e acbijH jcLibL ij π+ −+⋅ π+− −++κ+−κπ= +−+− , ))/()/(( 2 2 2 1 LjLiij π+πκ=Ω . Определим следующие значения параметров задачи: 121 == LL , 5021 == vv , 110−=κ . Поскольку 707/\|\|\|\| ≈κv � , то в данной задаче конвективные процессы значительно преобладают над процессом диффузии. 144 ISSN 1028-9763. Математичні машини і системи, 2014, № 4 Для интегрирования возникающих СОДУ во всех приведенных примерах нестацио- нарных задач использовались стандартный явный метод Рунге-Кутты 4-го порядка [11] и явный адаптивный метод 3-го порядка из работы [12]. Шаг по времени в текущем примере положим 310−=τ (начальный для метода из [12]. При этом в настройках выставлялись значения допустимой абсолютной и относительной ошибок [12] интегрирования, также равные 310− ). Отметим, что в приведенных примерах оба метода интегрирования приводят к практически идентичным результатам. Отметим также, что рассмотренные и проанали- зированные в работе (стабилизированные) методы конечных элементов, а также предла- гаемый вариант МПГ, осуществляют дискретизацию задач по пространственным перемен- ным, в результате чего получается СОДУ. Для интегрирования получившихся СОДУ в большинстве случаев оказывается вполне достаточным использования методов низкого порядка [1, 8, 9, 13]. Однако, поскольку погрешности и неустойчивости при применении конечноэлементных методов к решению задач КД обусловлены преимущественно дискре- тизациями пространственных членов и, следовательно, имеют скорее "пространственное" происхождение, чем "временнóе" [1–9, 13], то в данном разделе при рассмотрении приме- ров нестационарных задач для интегрирования СОДУ были использованы методы высоких порядков (3-го и 4-го). Это было сделано для того, чтобы минимизировать влияние дис- кретизации по временнóй переменной на окончательный результат и в явном виде выде- лить преимущество МПГ с ВФ (1) над другими стабилизированными методами простран- ственной дискретизации, рассмотренными в примерах. Рассмотрим сначала результаты расчетов на равномерной сетке 1515× узлов с тре- угольными элементами. Гипотенузы треугольников ориентированы поперек потока (из- вестно [8], что уменьшение угла между направлением потока и гипотенузами/диагоналями элементов приводит, как правило, к уменьшению погрешности, поэтому мы будем тести- ровать методы на наиболее "неудобном" с данной точки зрения расположении сетки). Гра- фик численного решения МКЭ Галеркина в зависимости от переменной 1xx ≡ при 4,0=t и фиксированном 5,02 =x приведен на рис. 4 (здесь и в дальнейшем на рисунках данного примера тонкой линией будет обозначаться истинное решение, а жирной – вычисленное). Видно, что численное решение является колебательным и не имеет ничего общего с ис- тинным решением. Обозначим через locerrmax величину максимального уклонения на графи- ке истинного решения от вычисленного, а через globerrmax – величину максимального (пото- чечного) уклонения на всей области. Тогда для МКЭ Галеркина имеем 47,1max ≈locerr , а 34,3max ≈globerr . Для метода SUPG с выбором параметров (5), а также для методов GLS и USFEM получаем 72,0max ≈locerr , а 55,1max ≈globerr (рис. 5). При этом для интегрирования СОДУ адаптивным методом [12] понадобился 341=M шаг, а показатель λ жесткости системы (отношение максимального модуля действительных частей собственных чисел к минимальному) равен 76,3 . Видно, что данные методы сильно сглаживают решение (стремятся его "выровнять" в окрестности пограничного слоя), что свидетельствует об их чрезмерной диссипативности в данной (нестационарной) задаче. Отметим, что использо- вание метода SUPG с оптимальным для двухмерных стационарных задач выбором пара- метров [8] (см. пример 1) не приводит к существенному улучшению решения. Для него 64,0max ≈locerr , 54,1max ≈globerr , 88,1≈λ , 232=M , однако решение имеет осцилляции в окре- стности пограничного слоя (рис. 6). ISSN 1028-9763. Математичні машини і системи, 2014, № 4 145 Рис. 4. Решение МКЭ Рис. 5. Решение SUPG Рис. 6. Решение SUPG c Галеркина (пример 2) c выбором СП в виде (5) оптимальным выбором СП Рис. 7. Решение МПГ с СП (4) Рис. 8. Решение МПГ с СП (4) Рис. 9. Решение МПГ с СП (4) и 5,1, =θ ki ( 1515× узлов) и 75,1, =θ ki ( 1515× узлов) и 75,1, =θ ki ( 2525× узлов) Для МПГ с ВФ (1) (для определенности используем кусочно-квадратичные функ- ции), параметрами (4) и 1, =θ ki получим 61,0max ≈locerr , 53,1max ≈globerr , 64,1≈λ , 233=M . Однако решение будет иметь небольшие осцилляции возле приграничного слоя (график решения практически совпадает с графиком на рис. 6, поэтому здесь не приводится). При 5,1, =θ ki получим 66,0max ≈locerr (рис. 7), 252=M и 45,2≈λ . Видно, что в окрестности по- граничного слоя имеется тенденция к осцилляции (образуется небольшая "ступенька"), хо- тя и значительно меньшая, чем для метода SUPG (рис. 6). Увеличение ki,θ позволяет уда- лить эту "ступеньку". Так, график решения для 75,1, =θ ki приведен на рис. 8 (при этом 67,0max ≈locerr , 52,1max ≈globerr , 98,2≈λ и 290=M ). Отметим, что использование параметров (3) в данном примере дает слегка худшие результаты. Так, для 75,1, =θ ki будет 66,0max ≈locerr , 363=M , однако в окрестности пограничного слоя наблюдается "ступенька" (аналогично ситуации на рис. 7), убрать которую можно еще увеличив ki,θ (однако при этом будет расти locerrmax). Измельчение сетки ведет к уменьшению погрешности (хотя со- отношения между погрешностями (в плане "больше" или "меньше") протестированных ме- тодов с различными вариантами СП, а также характер поведения решения в целом сохра- няются). Так, для разбиения 2525× узлов для МКЭ Галеркина получаем 64,0max ≈locerr , 47,3max ≈globerr , 02,2≈λ и 551=M . Для метода SUPG с (5) будет 5,0max ≈locerr , 49,1max ≈globerr , 69,3≈λ , 678=M . Для рассматриваемого МПГ с ВФ (1), параметрами (4) и 1, =θ ki получим 34,0max ≈locerr , 47,1max ≈globerr , 65,1≈λ , 502=M ; при 5,1, =θ ki будет 146 ISSN 1028-9763. Математичні машини і системи, 2014, № 4 42,0max ≈locerr , 48,1max ≈globerr , 36,2≈λ и 583=M ; при 75,1, =θ ki будет 45,0max ≈locerr , 48,1max ≈globerr , 96,2≈λ и 630=M (рис. 9). Пример 3. Рассмотрим начально-краевую задачу при ];[ bax ∈ для одномерного уравнения Бюргерса [5, 6]: f x u xx u u t u +      ∂ ∂κ ∂ ∂= ∂ ∂λ+ ∂ ∂ (7) при 1=λ и 0=f с известным аналитическим решением )))1(4/(exp()8/1exp(/)1(1 )1/( ),( 2 +κκ++ += txt tx xtu . Данная задача представляет собой мо- дель ударной волны [14] и применялась в [14–17] для тестирования предложенных там численных методов, основанных на методе RKPM и являющихся его дальнейшим усовер- шенствованием [14], а также методах, основанных на применении В-сплайнов, эрмитовых сплайнов и методов коллокаций [14–17]. Начальное (при 0=t ) и граничные условия полу- чаются из выписанного аналитического решения путем его непрерывного продолжения на гиперплоскости с 0=t и ax = , bx = соответственно. Для оценки уклонения численного решения u~ от аналитического u здесь и в дальнейших примерах используем величину \|),(),(~\|maxmax txutxuerr iii −≡ . Положим, 410−=κ , 0=a , 1=b . При использовании (3) положим [5, 6] ( 1) , n i k m+θ = θ ⋅ , где \|)(\|max )( , , n ki ki m γα≡ , константа 1=θ , а при использова- нии (4) положим 1)1( , =θ +n ki . В табл. 1 приведено значение величины maxerr (при 2=t ) в зависимости от выбора СП )1( , +α n ki и количества N равноотстоящих узлов сетки. Таблица 1. Значение maxerr (пример 3) )1( , +α n ki Количество узлов пространственной сетки, N 200 250 300 350 (3) 0,33203 0,32471 0,31836 0,31259 (3)* 0,12784 0,11491 0,08867 0,06558 (4) 0,19181 0,18664 0,18053 0,15949 (4)* 0,28647 0,28646 0,28644 0,28643 * при расчете использовалось сосредоточение [5, 18]. Здесь и в дальнейших примерах на рисунках жирной линией обозначен график вы- численного решения в зависимости от x , тонкой – график известного аналитического. На рис. 10 представлен результат расчета задачи МПГ (здесь и в дальнейших примерах в МГП используем кусочно-квадратичные ВФ (1)) при использовании формулы (3) и сосредото- чения [5, 18], 350=N . На рис. 11 область укручения решения показана крупным планом. На рис. 12 показана область укручения решения при использовании формулы (4), 350=N , а на рис. 13 – при использовании МКЭ Галеркина с сосредоточением. Видно, что в по- следнем случае решение является осциллирующим (в силу преобладания конвекционных процессов и нелинейности); отметим, что в случае, когда сосредоточение не используется, осцилляции будут гораздо бóльшими, поскольку сосредоточение, как известно [5, 18], уве- личивает устойчивость схем и осуществляет их "монотонизацию". При использовании (3) без сосредоточения, а также при использовании (4) с сосредоточением погрешность боль- ISSN 1028-9763. Математичні машини і системи, 2014, № 4 147 ше (по сравнению с другими вариантами, перечисленными в табл. 1) и с возрастанием N она медленно убывает. При использовании МПГ с постоянными СП (не зависящими от индексов узлов сетки и момента времени t ) решение получается неустойчивым и колеба- тельным в окрестностях тонких слоев, а в случае, когда применяется сосредоточение — чрезмерно диссипативным и сглаженным. Так, на рис. 14 изображен результат решения при прежних условиях, с использованием сосредоточения и СП, тождественно равными единице. Видно, что решение чрезмерно сглажено в окрестности слоя. На рис. 15 показан соответствующий результат, когда сосредоточение не используется. Как видим, в таком случае появляются осцилляции. Увеличение СП не приводит к улучшению решения и ис- чезновению осцилляций (а при использовании сосредоточения решение становится еще более сглаженным, чем на рис. 14). Уменьшение СП также не приводит к улучшению ре- шения (а при стремлении их к нулю получим МКЭ Галеркина, который порождает осцил- ляции независимо от того, использовалось сосредоточение или нет подобно ситуации на рис. 13). При интегрировании СОДУ здесь и в дальнейших примерах использовались ме- тоды, описанные в примере 2 (в данной задаче с шагом по времени 410−=τ ). Рис. 10. Решение МПГ с СП (3) и Рис. 11. Решение МПГ с СП (3) и использованием сосредоточения сосредоточением (область укручения) Рис. 12. Решение МПГ с СП (4) Рис. 13.Решение МКЭ Галеркина с (область укручения) использованием сосредоточения Рассмотрим теперь данную задачу при условиях 005,0=κ , 0=a , 2,1=b с шагом по времени 01,0=τ и по пространству 0012,0=h (то есть 1001=N ). При таких значениях параметров задача рассматривалась в [14] (однако там был взят еще более мелкий шаг h , 001,0=h ). Тогда при использовании (3) ( 1=θ ) с сосредоточением для моментов времени 7,1=t , 5,2=t , 3=t , 5,3=t имеем 5 max 1067 −⋅≈err , 5 max 1062 −⋅≈err , 5 max 1089 −⋅≈err , 5 max 1036 −⋅≈err соответственно. Для метода работы [14] с использованием пространств кусочно-линейных функций соответствующие значения maxerr равны 4 max 1047,13 −⋅≈err , 4 max 1055,15 −⋅≈err , 4 max 1053,15 −⋅≈err , 4 max 1022,15 −⋅≈err [14]. 148 ISSN 1028-9763. Математичні машини і системи, 2014, № 4 Рассмотрим еще задачу при условиях 005,0=κ , 0=a , 1=b с шагом по времени 01,0=τ и по пространству 005,0=h . При таких значениях параметров задача рассматри- валась в [14, 15]. Тогда, аналогично, при использовании (3) с сосредоточением для момен- тов времени 7,1=t , 5,2=t , 25,3=t имеем 5 max 1048 −⋅≈err , 5 max 105,18 −⋅≈err , 5 max 1063 −⋅≈err соответственно. Для метода QRKM (Quadratic Reproducing Kernel function Method), предложенного в работе [14], соответственно имеем 3 max 10092,0 −⋅≈err , 3 max 10115,0 −⋅≈err , 3 max 108 −⋅≈err (табл. 11 на стр. 433 в [14]). Для методов QBCM (Qua- dratic B-spline Collocation Method) и CBCM (Cubic B-spline Collocation Method), предло- женных в работе [15] и основанных на использовании В-сплайнов, имеем соответственно значения [14, 15] 3 max 10312,0 −⋅≈err , 3 max 10189,0 −⋅≈err , 3 max 1098,8 −⋅≈err и 3 max 10577,27 −⋅≈err , 3 max 10152,25 −⋅≈err , 3 max 10049,21 −⋅≈err . Таким образом, приве- денные расчеты показывают, что предлагаемая версия МПГ с адаптивным выбором СП в виде (3) и использованием сосредоточения дает лучшие результаты (в том числе даже и на более грубых сетках), чем методы, предложенные в работах [14, 15]. Рис. 14. Решение МПГ с СП, равными Рис. 15. Решение МПГ с СП, единице, и сосредоточением равными единице Для сравнения c методами, предложенными в работах [16, 17], рассмотрим данную задачу при условиях 005,0=κ , 0=a , 1=b , 001,0=τ (при таких условиях она рассматри- валась в указанных работах). Тогда при использовании (3) ( 2/1=θ ) с сосредоточением при 61=N , 81=N , 101=N и 121=N для фиксированного момента 6,3=t имеем 3 max 10,2380 −⋅≈err , 3 max 10,1340 −⋅≈err , 3 max 10,0880 −⋅≈err , 3 max 10,0590 −⋅≈err соответ- ственно. При использовании (4) (при прежних условиях) соответственно имеем 3 max 10,3540 −⋅≈err , 3 max 10,1460 −⋅≈err , 3 max 10,1560 −⋅≈err , 3 max 10,0460 −⋅≈err . Для ме- тода работы [16] будем иметь 3 max 10,470 −⋅≈err , 3 max 10,270 −⋅≈err , 3 max 10,170 −⋅≈err , 3 max 10,120 −⋅≈err соответственно (см. [16]), а для метода работы [17] — соответственно получаем 3 max 10,450 −⋅≈err , 3 max 10,310 −⋅≈err , 3 max 10,230 −⋅≈err , 3 max 10,160 −⋅≈err [16, 17]. Таким образом, данные расчеты показывают, что предлагаемая версия МПГ с адап- тивным выбором СП в виде (3) и использованием сосредоточения, а также с выбором па- раметров в виде (4) (без использования сосредоточения) дает в данном примере более точ- ные результаты, чем методы работ [16, 17]. Отметим, что для интегрирования полудис- кретных аппроксимаций в [16] использовались методы Рунге-Кутты SSP-RK43 и SSP- RK54 (3-го, 4-го и 5-го порядков) (детально в [16]). ISSN 1028-9763. Математичні машини і системи, 2014, № 4 149 Пример 4. Рассмотрим начально-краевую задачу при ];[ bax ∈ для одномерного уравнения Бюргерса (7) при 1=λ и 0=f с известным аналитическим решением )cos()exp(2 )sin()exp(2 ),( 2 2 xt xt xtu π⋅κπ−+ π⋅κπ−⋅πκ= . Данная задача применялась в [16, 19] для проверки и тес- тирования предложенных там численных методов решения уравнения Бюргерса. Положим, 0=a , 1=b , 40/1=h , 410−=τ . Тогда при использовании (3) с сосредоточением (подроб- ный выбор СП описан в примере 3) в фиксированный момент времени 310−=t для 10/5=κ , 10/2=κ и 10/1=κ имеем 5 max 101,2741 −⋅≈err , 6 max 102,2121 −⋅≈err и 6 max 100,5987 −⋅≈err соответственно. При использовании формул (4) соответственно полу- чим 5 max 10,07524 −⋅≈err , 6 max 10,66536 −⋅≈err , 6 max 101,6787 −⋅≈err . Соответствующие значения величины maxerr для метода работы [16] равны 5 max 10,447 −⋅≈err , 5 max 101,22 −⋅≈err и 6 max 10,083 −⋅≈err , а для метода работы [19] – 5 max 10,002 −⋅≈err , 6 max 10,003 −⋅≈err и 6 max 10,002 −⋅≈err (табл. 3 работы [19] на стр. 575). Таким образом, приведенные расчеты показывают, что предлагаемая версия МПГ с адаптивным выбором СП в виде (3) и использованием сосредоточения дает в данном примере более точные ре- зультаты, чем методы, предложенные в работах [16, 19]. Пример 5. Рассмотрим начально-краевую задачу при ];[ bax ∈ для одномерного уравнения Бюргерса (7) при 1=λ и 0=f с известным аналитическим решением η η + α−µ+µ+α= e e xtu 1 )( ),( , где κ αβ−µ−=η )( tx . Данная задача применялась в [16, 20, 21] для проверки и тестирования предложенных там численных методов решения уравнения Бюргерса. Ее решение представляет собой движущуюся "ступеньку" (решение типа "кинк"). Положим, 4,0=α , 125,0=β , 6,0=µ , 210−=κ , 210−=τ , 36/1=h [16, 20]. Тогда при использовании (3) с сосредоточением (подробный выбор СП описан в примере 3) в фиксированный момент времени 2/1=t при различных 0=θ , 4/1=θ , 2/1=θ , 4/3=θ и 1=θ имеем соответственно 0,03693max ≈err , 0,03045max ≈err , 0,03239max ≈err , 0,04max ≈err , 0,04670max ≈err . Аналогично, при использовании (4) и )1( , +θ n ki , приобретаю- щем значения 0 , 4/1 , 2/1 , 4/3 и 1 соответственно, имеем 0,005max ≈err , 0,00409max ≈err , 0,00604max ≈err , 0,00928max ≈err , 0,01246max ≈err . Поскольку в данной задаче нет сильного доминирования конвекции (число Пекле имеет порядок 210 ) и резких неоднородностей в решении, то, соответственно, нет нужды во введении сильной сглажи- вающей диссипации, поэтому наиболее эффективными оказались "мягкие" варианты вы- бора СП (в особенности, (4)) и при небольших значениях параметров θ и )1( , +θ n ki . Для ме- тода работы [16] при тех же условиях расчета имеем 0,00597max ≈err , а для конечноэле- ментного метода работы [20] – 0,045max ≈err [16, 20]. Положим, теперь 18/1=h и 310−=τ (такие значения параметров использовались в [21]). Тогда при использовании (4) и )1( , +θ n ki , приобретающем значения 0 , 4/1 , 2/1 , 4/3 и 1 соответственно, имеем 0,02494max ≈err , 0,01347max ≈err , 0,02165max ≈err , 0,03432max ≈err , 0,0481max ≈err . При указанных значениях параметров для конечноэле- 150 ISSN 1028-9763. Математичні машини і системи, 2014, № 4 ментного метода SGA (Standard Galerkin Approach) имеем 0,096max ≈err , для конечноэле- ментного метода PAG (Product Approximation Galerkin) 0,082max ≈err , для метода CD (Compact finite Difference) 0,151max ≈err (описание данных методов и указанные значения maxerr в [20, 21]). Таким образом, расчеты показывают, что и в данном примере предла- гаемая версия МПГ с адаптивным выбором СП (3)–(4) способна дать более точные резуль- таты, чем методы, предложенные в работах [16, 20, 21]. Пример 6. Рассмотрим начально-краевую задачу для одномерного уравнения Бюр- герса (7) при 1=λ и 0=f на отрезке ]1;0[ с известным аналитическим решением µ+κ− κ−κ− µ+κ− −µκ +κ= )ch( )sh( )ch( 1 ),( 2 tx tx tx xtu , где µ – некоторая константа. Данное решение имеет солитонный тип [22] и применялось в работах [16, 22] для проверки и тестирования предложенных там численных методов. Положим, 51075 −⋅=κ , 3−=µ , 210−=h , 210−=τ (при таких условиях задача использовалась в [16, 22]). Тогда при использовании (3) (c 1=θ , пример 3) с сосредоточением в фиксированный момент времени 20/1=t имеем 13 max 10−≈err . При использовании (4) ( 1)1( , =θ +n ki ) имеем 9 max 102 −⋅≈err . Для метода рабо- ты [16] получаем 8 max 102,1 −⋅≈err (табл. 10.1 в [16]), а для 12-точечной схемы работы [22] – 4 max 101,058 −⋅≈err (табл. 6 в [22]). Пример 7. Рассмотрим начально-краевую задачу для уравнения (7) при 1=λ и 0=f на отрезке ]1;0[ с нулевыми граничными условиями первого рода и начальным ус- ловием )sin()(0 xxu π= . Отметим, что задача в такой постановке дает описание нелинейных стоячих медленных волн в корональных магнитных петлях [23] в солнечной плазме. Ана- литическое решение начально-краевой найдено в [20, 23]. Возьмем сначала значение 410−=κ . В силу необходимости выполнения нулевого граничного условия вблизи точки 1=x возникает тонкий приграничный слой. Сравним погрешность предлагаемой версии МПГ с погрешностями методов, использованных в работе [20] при 2/1=t на наборе точек 5,0=x , 556,0=x , 611,0=x , 667,0=x , 722,0=x , 778,0=x , 833,0=x , 889,0=x , 944,0=x (именно такой тестовый набор точек был взят в [20]). Взяв 216/1=h и 8/1=τ [20] для ко- нечноэлементного метода работы [20], получим 097,0max ≈err , для конечноэлементного метода работы [24] с кубическими В-сплайнами получим 078,0max ≈err (табл. 3 в [20]), а для предлагаемой версии МПГ с (3) (при 1=θ ) и использования сосредоточения получим 0,022max ≈err (при 3,1=θ получим 0,02max ≈err ). Положим, теперь 18/1=h , 310−=τ . То- гда для метода SGA получим 0,124max ≈err , для метода PAG – 0,105max ≈err , а для мето- да CD – 0,087max ≈err (по поводу данных методов см. табл. 3 в [20]). Для МПГ с (3) (при 1=θ ) и сосредоточением получим 0,0617max ≈err . Для МПГ с сосредоточением и посто- янными СП, равными единице, получим 0,063max ≈err (а решение получается чрезмерно сглаженным возле приграничного слоя), а без использования сосредоточения получим 0,119max ≈err . Для МКЭ Галеркина решение получается колебательным вблизи пригра- ничного слоя и имеет большую погрешность (с использованием сосредоточения будем иметь 0,839max ≈err , а без его использования – 0,975max ≈err ). Пример 8. Рассмотрим начально-краевую задачу для уравнения Бюргерса (7) при 1−=λ и 0=f на отрезке ]1;0[ с известным аналитическим решением ISSN 1028-9763. Математичні машини і системи, 2014, № 4 151 ))4/(exp())2/(ch( ))2/(sh( ),( κ−+κ κ= tx x xtu . Данная задача использовалась в [25] для тестирования предложенных там неявных численных схем. Положим, 51075 −⋅=κ , 20/1=h , 210−=τ (при таких условиях проводились расчеты в [25]). Тогда при использовании (3) ( 1=θ ) с сосредоточением в фиксированный момент времени 20/1=t имеем 2 max 105 −⋅≈err . При использовании (4) ( 1)1( , =θ +n ki ) имеем 2 max 105,6 −⋅≈err . Для 8-точечной неявной схемы ра- боты [25] имеем 2 max 107 −⋅≈err , а для симметричной схемы Кранка-Николсона – 1 max 105 −⋅≈err (табл. 4 в [25]). Пример 9. Рассмотрим начально-краевую задачу для двухмерного уравнения Бюр- герса [6, 26–28]: f x u xx u xx u u x u u t u +      ∂ ∂κ ∂ ∂+      ∂ ∂κ ∂ ∂= ∂ ∂λ+ ∂ ∂λ+ ∂ ∂ 22112 2 1 1 (где ),( xtii �λ=λ , ),( xt �κ=κ , ),( xtff �= ) при 121 =λ=λ , 20/1=κ , 0=f , ]1;0[]1;0[),( 21 ×∈xx с известным аналитическим решением 1 2121 )))2/()exp((1(),,( −κ−++= txxxxtu [26–28]. Данное решение служит моделью ударной волны, волновой фронт которой движется вдоль линии txx −+ 21 (интерпретацию реше- ния и его гидродинамическую суть см. в [26, 27]). Используем равномерное разбиение 2121× узлов и положим 310−=τ . Тогда, при использовании (4) (с 2,0)1( , =θ +n ki ) для момен- тов времени 2/1=t , 4/3=t , 1=t , 4/5=t , имеем соответственно 0,00032max ≈err , 0,00097max ≈err , 0,0007max ≈err , 0,00056max ≈err . Отметим, что увеличение или умень- шение )1( , +θ n ki , как показывают расчеты, ведет к (медленному) росту ошибки (по сравнению с ошибкой для 2,0)1( , =θ +n ki ); это, как уже отмечалось ранее, связано с тем, что коэффици- ент вязкости κ не является малой величиной, следовательно, конвекция в данной задаче не является доминирующей, поэтому небольшие значения коэффициентов )1( , +θ n ki приводят к более точным решениям (однако при нулевых значениях )1( , +θ n ki погрешность увеличива- ется). Для meshless-метода ELMFS (Eulerian-Lagrangian Method of Fundamental Solutions) работы [26] при тех же условиях интегрирования задачи имеем соответственно 0,00061max ≈err , 0,00105max ≈err , 0,00299max ≈err , 0,00303max ≈err (табл. 1 в работе [26] на стр. 399). Рассмотрим теперь задачу при разбиении 2020× узлов и положим 210−=τ (при таких условиях задача рассматривалась в [27]). Тогда на промежутке 4/50 ≤≤ t для МПГ с (4) имеем 007,0max ≈err . Для методов MFS-DRM (Method of Fundamental Solution - Dual Reciprocity Method) и Kansa работы [27] при тех же условиях решения имеем соответ- ственно 0703,0max ≈err и 0719,0max ≈err (стр. 218 в [27]). Пусть теперь 1=κ , 410−=τ и используется разбиение 1111× узлов. Тогда при 2/1=t для МПГ с (4) ( 1)1( , =θ +n ki ) имеем 7 max 1072,8 −⋅≈err , а для метода LMAPS (Local Method of Approximate Particular Solutions) работы [28] – 6 max 109,6 −⋅≈err (табл. 1 в [28]; отметим, однако, что авторы [28] не указали в своей работе ни метода интегрирования своих полудискретных аппроксимаций, ни кон- 152 ISSN 1028-9763. Математичні машини і системи, 2014, № 4 кретного значения шага по времени, с которым они решали данный пример). Таким обра- зом, предлагаемая версия МПГ с выбором СП в виде (4) в данной задаче дает более высо- кую точность, чем методы работ [26–28]. Отметим, что для многомерных задач использо- вание сосредоточения является "рискованным" (в силу того, что может ввести диффузию поперек потока [1], которая приведет к возникновению больших погрешностей), а вариант (3) выбора СП без сосредоточения показывает плохие результаты даже в одномерном слу- чае (пример 3), поэтому при решении многомерных задач вариант (4) выбора СП является предпочтительным. Пример 10. Рассмотрим задачу для системы Бюргерса [26, 28]:              ∂ ∂κ ∂ ∂+      ∂ ∂κ ∂ ∂= ∂ ∂λ+ ∂ ∂λ+ ∂ ∂       ∂ ∂κ ∂ ∂+      ∂ ∂κ ∂ ∂= ∂ ∂λ+ ∂ ∂λ+ ∂ ∂ 2 2 21 2 12 2 22 1 2 11 2 2 1 21 1 12 1 22 1 1 11 1 , x u xx u xx u u x u u t u x u xx u xx u u x u u t u при 121 =λ=λ , 100/1=κ , ]1;0[]1;0[),( 21 ×∈xx с аналитическим решением [26] 1 1 1 2 1 2( , , ) 3 / 4 (1 exp(( 4 4 ) / (32 ))) / 4u t x x x x t −= − + − + − κ , 1 2 1 2 1 2( , , ) 3 / 4 (1 exp(( 4 4 ) / (32 ))) / 4u t x x x x t −= + + − + − κ . Данное решение служит моделью ударной волны (гидродинамическую интерпрета- цию решения см. в [26]) и использовалось в работе [26] для проверки работоспособности предложенного там метода ELMFS (пример 9). Используем равномерное разбиение 1111× узлов и положим 005,0=τ . Тогда, при использовании (4) (с 2,0)1( , =θ +n ki , см. комментарии к примеру 9) для моментов времени 210−=t , 2/1=t , 2=t для 1u имеем соответственно 5 max 1007,3 −⋅≈err , 3 max 1032,7 −⋅≈err , 3 max 1054,7 −⋅≈err , для величины 2u погрешности такие же (с точностью до малых более высокого порядка). Для метода ELMFS работы [26] при тех же условиях интегрирования задачи соответствующие погрешности для 1u имеют вид 4 max 1039,5 −⋅≈err , 3 max 1066,7 −⋅≈err , 2 max 1003,1 −⋅≈err (для величины 2u они будут такими же, см. табл. 3 и 4 в работе [26] на стр. 404-405). Для метода LMAPS работы [28] (см. комментарии к примеру 9) при тех же условиях и 2/1=t имеем 3 max 104,7 −⋅≈err (табл. 2 в [28]). Таким образом, предлагаемая версия МПГ с выбором стабилизирующих параметров в виде (4) в данной задаче дает более высокую точность, чем методы работ [26, 28]. Пример 11. Рассмотрим важный случай системы уравнений магнитной гидродина- мики (в предположении адиабатичности и бесконечной проводимости среды [7, 29]), в ко- торой все величины зависят от одной пространственной координаты, например, координа- ты z [29]. Тогда имеем вектор скорости ( ( , ), ( , ), ( , ))x y zu u t z u t z u t z=� , магнитную ин- дукцию ( ( , ), ( , ), ( , ))x y zB B t z B t z B t z= � , плотность ( , )t zρ = ρ и давление ( , )p p t z= . К этой ситуации сводится математическое описание распространения волн при возмущении какой-либо из МГД величин на границе рассматриваемой области, а также движение плаз- мы в конфигурациях, которые могут быть описаны (по крайней мере, в некотором при- ближении) единственной пространственной координатой (например, явление спикул [29]). Введем сетку iz , ni ..1= на отрезке ];[ 21 LLz ∈ , 11 Lz = , 2Lzn = . Начальные усло- вия, соответствующие стационарному равновесному состоянию плазмы, имеют вид ISSN 1028-9763. Математичні машини і системи, 2014, № 4 153 00 == uu �� , 00 == BB �� , ρ=ρ== CCpp 00 (где C — некоторая константа), Cgzez / 00 )0()( −ρ=ρ=ρ , где сила тяжести (направлена противоположно оси z ) const=g . Соответственно, 00 ρ= Cp . По компоненте скорости zu на левом конце области дадим следующее синусоидальное возмущение: sin(2 / (3 / 20)), 0 3 / 20, ( ) 0, 3 / 20.z t t u t t π ≤ ≤ δ =  > Без ограничения общности, для демонстрации работоспособности обсуждаемой версии МПГ положим при расчетах 2/1=C , 1)0(0 =ρ , 1=g и коэффициент вязкости 210−=υ [7]. Тогда при использовании МПГ с выбором СП в виде (3) с сосредоточением, 1=θ (пример 3), 210−=h (используем плавающую сетку работы [30]) и 210−=τ графики плотности ρ и компоненты zu скорости в момент времени 15,0=t в зависимости от коор- динаты z изображены на рис. 16 и 17 соответственно. На рис. 18 и 19 показаны графики соответствующих величин в момент времени 25,0=t . Рис. 16. График плотности ρ при 15,0=t Рис. 17. График скорости zu при 15,0=t Рис. 18. График плотности ρ при 25,0=t Рис. 19. График скорости zu при 25,0=t Видно, что возмущение )(tuzδ сохраняет свою форму, но с течением времени зату- хает по амплитуде, что объясняется действием диссипации (вязкость с коэффициентом υ ) и нелинейности конвективного члена, причем, чем бóльшим будет значение коэффициента вязкости υ , тем сильнее будет затухать амплитуда скорости с течением времени. Возму- щение скорости приводит к возмущению плотности (и, в свою очередь, давления) на фоне равновесного установившегося распределения плотности Cgzez / 00 )0()( −ρ=ρ . Как видно из рис. 16 и 18, происходит "скатывание" возмущения со стороны левого конца области. С течением времени данное возмущение затухает (поскольку затухает возмущение перемен- ной скорости) и плазма, при отсутствии других внешних возмущений, постепенно перехо- дит в равновесное состояние. Таким образом, полученные результаты расчетов полностью согласуются с физической картиной поведения рассматриваемого процесса. 154 ISSN 1028-9763. Математичні машини і системи, 2014, № 4 5. Заключение В настоящей статье получили дальнейшее развитие результаты работ [5–7], касающиеся выбора и динамической (адаптивной) настройки стабилизирующих параметров в конечно- элементном методе Петрова-Галеркина с предложенными в работах [2–4] весовыми (пове- рочными) функциями. На многочисленных примерах (разд. 4) продемонстрирована высо- кая эффективность метода Петрова-Галеркина с весовыми функциями (1) и предложенны- ми вариантами выбора стабилизирующих параметров (разд. 3) при решении различных за- дач конвекции-диффузии (стационарных и нестационарных линейных задач, а также нели- нейных уравнений и систем типа Бюргерса). Во всех приведенных примерах рассматри- ваемый метод Петрова-Галеркина оказался способным дать более точное решение, чем ме- тоды работ [8, 13–17, 19–22, 24–28], в которых данные примеры (те или иные) также при- менялись для проверки точности и работоспособности предложенных там методов. Расче- ты показали, что в целом, в случае одномерных задач с наличием тонких слоев и ударных волн (особенно задач для уравнений Бюргерса), среди вариантов (3) и (4) выбора стабили- зирующих параметров более эффективным (в плане точности и подавления осцилляций), как правило, оказывается "жесткий" вариант выбора (3) с одновременным использованием сосредоточения. В случае же, когда конвекция не является доминирующей над процессом диффузии, а также для многомерных задач, "мягкий" вариант (4) выбора стабилизирую- щих параметров в приведенных примерах оказывается более эффективным по сравнению с вариантом (3). Благодарности. Автор выражает глубокую благодарность к.т.н., с.н.с. ИКИ НАНУ Сальникову Н.Н. и д.т.н., проф. Молчанову А.А. за плодотворное обсуждение работы и ценные замечания к ней. СПИСОК ЛИТЕРАТУРЫ 1. Roos H.-G. Robust numerical methods for singularly perturbed differential equations / H.-G. Roos, M. Stynes, L. Tobiska. – Berlin, Heidelberg: Springer-Verlag, 2008. – 604 p. 2. Сальников Н.Н. О построении конечномерной математической модели процесса конвекции- диффузии с использованием метода Петрова-Галеркина / Н.Н. Сальников, С.В. Сирик, И.А. Тере- щенко // Проблемы управления и информатики. – 2010. – № 3. – С. 94 – 109. 3. Сальников Н.Н. Построение весовых функций метода Петрова-Галеркина для уравнений кон- векции-диффузии-реакции в трехмерном случае / Н.Н. Сальников, С.В. Сирик // Кибернетика и системный анализ. – 2014. – № 5. – С. 173 – 183. 4. Сирик С.В. Выбор весовых функций в методе Петрова-Галеркина для интегрирования линейных одномерных уравнений конвекции-диффузии / С.В. Сирик, Н.Н. Сальников, В.К. Белошапкин // Управляющие системы и машины. – 2014. – № 1. – С. 38 – 47. 5. Сирик С.В. Численное интегрирование уравнения Бюргерса методом Петрова-Галеркина с адап- тивными весовыми функциями / С.В. Сирик, Н.Н. Сальников // Проблемы управления и информа- тики. – 2012. – № 1. – C. 94 – 110. 6. Молчанов А.А. Выбор весовых функций в методе Петрова-Галеркина для интегрирования двух- мерных нелинейных уравнений типа Бюргерса / А.А. Молчанов, С.В. Сирик, Н.Н. Сальников // Математичні машини і системи. – 2012. – № 2. – С. 136 – 144. 7. Сальников Н.Н. О построении конечномерных математических моделей для двухмерных про- цессов магнитной гидродинамики с использованием метода Петрова-Галеркина / Н.Н. Сальников, С.В. Сирик, В.К. Белошапкин // Управляющие системы и машины. – 2014. – № 4. – С. 23 – 32. 8. Harari I. Streamline design of stability parameters for advection-diffusion problems / I. Harari, L.P. Franca, S.P. Oliveira // Journal of Computational Physics. – 2001. – Vol. 171 (1). – P. 115 – 131. 9. John V. Finite element methods for time-dependent convection–diffusion–reaction equations with small diffusion / V. John, E. Schmeyer // Comput. Methods Appl. Mech. Engrg. – 2008. – Vol. 198. – P. 475 – 494. 10. Хорн Р. Матричный анализ / Р. Хорн, Ч. Джонсон; пер. с англ. – М.: Мир, 1989. – 655 с. 11. Калиткин Н.Н. Численные методы / Калиткин Н.Н. – М.: Наука, 1978. – 512 с. ISSN 1028-9763. Математичні машини і системи, 2014, № 4 155 12. Скворцов Л.М. Простые явные методы численного решения жестких обыкновенных дифферен- циальных уравнений / Л.М. Скворцов // Вычислительные методы и программирование. – 2008. – Т. 9. – С. 154 – 162. 13. Brooks A.N. Streamline 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. – Vol. 32 (1 – 3). – P. 199 – 259. 14. Xie S.-S. Numerical solution of one-dimensional Burgers' equation using reproducing kernel function / S.-S. Xie, S. Heo, S. Kim [et al.] // J. Comput. Appl. Math. – 2008. – Vol. 214. – P. 417 – 434. 15. Dag I. B-spline collocation methods for numerical solutions of the Burgers equation / I. Dag, D. Irk, A. Sahin // Math. Problems in Eng. – 2005. – Vol. 5. – P. 521 – 538. 16. Mittal R.C. Numerical solutions of nonlinear Burgers equation with modified cubic B-splines colloca- tion method / R.C. Mittal, R.K. Jain // Applied Mathematics and Computation. – 2012. – Vol. 218. – P. 7839 – 7855. 17. Korkmaz A. Polynomial based differential quadrature method for numerical solution of nonlinear Burgers equation / A. Korkmaz, I. Dag // Journal of the Franklin Institute. – 2011. – Vol. 348 (10). – P. 2863 – 2875. 18. Сирик С.В. О применении сосредоточения в методе конечных элементов Петрова-Галеркина при решении задач конвекции-диффузии / С.В. Сирик, Н.Н. Сальников // Доповіді НАН України. – 2014. – № 5. – С. 39 – 44. 19. Ganaie I.A. Numerical solution of Burgers equation by cubic Hermite collocation method / I.A. Ga- naie, V.K. Kukreja // Applied Mathematics and Computation. – 2014. – Vol. 237. – P. 571 – 581. 20. Dogan A. A Galerkin finite element approach to Burgers' equation / A. Dogan // Applied Mathematics and Computation. – 2004. – Vol. 157. – P. 331 – 346. 21. Christie I. Product approximation for nonlinear problems in the finite element method / I. Christie, D.F. Griffiths, A.R. Mittchell [et al.] // IMA. J. Numer. Anal. – 1981. – Vol. 1. – P. 253 – 266. 22. Hesameddini E. Soliton and numerical solutions of the Burgers Equation and comparing them / E. He- sameddini, R. Gholampour // Int. J. Math. Anal. – 2010. – Vol. 4 (52). – P. 2547 – 2564. 23. Ruderman M.S. Nonlinear damped standing slow waves in hot coronal magnetic loops / M.S. Ruder- man // Astronomy & Astrophysics. – 2013. – Vol. 553 (A23). – P. 1 – 6. 24. Ali A.H.A. A collocation solution for Burgers equation using cubic B-spline finite elements / A.H.A. Ali, G.A. Gardner, L.R.T. Gardner // Comput. Meth. Appl. Mech. Eng. – 1992. – Vol. 100. – P. 325 – 337. 25. Tabatabaei A.H. Some implicit methods for the numerical solution of Burgers equation / A.H. Tabata- baei, E. Shakour, M. Dehghan // Applied Mathematics and Computation. – 2007. – Vol. 191. – P. 560 – 570. 26. Young D.L. The Eulerian-Lagrangian method of fundamental solutions for two-dimensional unsteady Burgers' equations / D.L. Young, C.M. Fan, S.P. Hu [et al.] // Engineering Analysis with Boundary Ele- ments. – 2008. – Vol. 32. – P. 395 – 412. 27. Li J. Numerical comparisons of two meshless methods using radial basis functions / J. Li, Y.C. Hon, C.S. Chen // Engineering Analysis with Boundary Elements. – 2002. – Vol. 26. – P. 205 – 225. 28. Zhang X. Local method of approximate particular solutions for two-dimensional unsteady Burgers' equations / X. Zhang, H. Tian, W. Chen // Computers and Mathematics with Applications. – 2014. – Vol. 66. – P. 2425 – 2432. 29. Ладиков-Роев Ю.П. Математические модели сплошных сред / Ю.П. Ладиков-Роев, О.К. Черем- ных. – К.: Наукова думка, 2010. – 551 с. 30. Молчанов О.А. Проекційні методи інтегрування рівнянь магнітної гідродинаміки / О.А. Молча- нов, М.М. Сальніков, С.В. Сірик // ІІІ наукова конференція магістрантів та аспірантів «Прикладна математика та комп'ютинг»: тези доповідей, (Київ, 13-15 квітня 2011). – К.: Просвіта, 2011. – С. 272 – 278. Стаття надійшла до редакції 21.10.2014

О выборе стабилизирующих параметров в конечноэлементном методе Петрова-Галеркина при решении задач конвекции-диффузии

Institution

Similar Items