CC BY
43
УДК 519.63
Д. Б. Жамбалова 1, С. Г. Черный 2
Новосибирский государственный университет ул. Пирогова, 2, Новосибирск, 630090, Россия
Институт вычислительных технологий СО РАН пр. Акад. Лаврентьева, 6, Новосибирск, 630090, Россия
E-mail: 1 zhambalova@gmail.com; 2 cher@ict.nsc.ru
МЕТОД ИНТЕРПОЛЯЦИОННОГО ПРОФИЛЯ РЕШЕНИЯ УРАВНЕНИЙ ПЕРЕНОСА *
Представлен метод интерполяционного профиля, одновременно сочетающий высокий порядок аппроксимации и учет структуры решений гиперболических уравнений. Предложены модификации метода для решения уравнений в дивергентной форме и для многомерного случая. Для повышения разрешающей способности метода применяется подход, заключающийся в решении уравнения переноса образа тангенциального преобразования искомой функции. Приводятся результаты расчетов одномерных и двумерных задач, в том числе имеющих разрывные решения.
Ключевые слова: метод интерполяционного профиля, гиперболические уравнения, дифференциальное приближение, метод характеристик.
Введение
В ряде задач вычислительной гидродинамики, таких как, например, моделирование течений с подвижными границами, требуется с большой точностью находить решения уравнений конвективного переноса. В обычных конечно-разностных методах при определении производных по значениям сеточных функции в узлах теряется информация о градиенте функции внутри ячейки. Это приводит к нежелательному сглаживанию решения. Поэтому в последнее время широкое распространение получил так называемый метод интерполяционного профиля (МИП) - Constrained Interpolation Profile (CIP), свободный от указанного недостатка. В нем функция аппроксимируется не только по своему значению, но и по значению производной, согласованной с решаемым уравнением переноса. Метод заключается в построении на каждом интервале сетки кубического полинома, удовлетворяющего условию гладкого сопряжения в узлах. Тогда на последующий временной слой решение переносится по характеристике дифференциального уравнения. Метод был разработан в середине 1980-х гг. группой японских ученых [1]. Приведем краткий обзор основных работ, посвященных МИП.
Работа [2] посвящена общему обзору различных семейств МИП. После описания основной идеи метода на примере одномерной задачи представлены его вариации, такие как: формулировка в многомерном случае; интерполяция, уменьшающая осцилляции; метод сохранения первоначального профиля. Представлена формулировка МИП для общих задач гидродинамики. На нескольких примерах показаны приложения метода для моделирования течений сжимаемой и несжимаемой жидкости. Для многофазной среды предложена числен-
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 11-01-00475-а).
ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2012. Том 10, выпуск 1 © Д. Б. Жамбалова, С. Г. Черный, 2012
ная формулировка задачи сохранения свойств поверхностного натяжения, упругопластиче-ских эффектов.
В работе [3] описываются основные подходы к решению задач вычислительной гидродинамики. Формулируются уравнения движения для сжимаемой и несжимаемой жидкости в трехмерной постановке, базовые определения и понятия. Для решения приведенных уравнений рассмотрены конечно-разностные методы, метод конечного объема, отдельно приведено несколько модификаций противопоточной схемы, описаны общие принципы МИП. В двумерном и трехмерном случаях применен метод дробных шагов, позволяющий на каждом шаге решать одномерные уравнения. Авторами также предложен способ расщепления, предупреждающий возникновение больших градиентов в решении.
В работе [4] приведена одномерная формулировка метода в консервативном виде. Для аппроксимации интегрального закона сохранения предложено использовать полиномы второго или четвертого порядка. Так, применение полинома четвертой степени улучшает аппроксимацию, но усложняет вычисления при решении многомерных задач. Приложение метода, использующего полиномы второго порядка, к многомерной модели не требует больших вычислительных затрат. В работе также представлены расчеты для нелинейных уравнений переноса при большом числе Куранта.
Недостаток описанных в [4] консервативных модификаций МИП - отсутствие в них механизма взаимодействия характеристик на разрыве решения и, как следствие, несоответствие численного решения точному решению уравнения переноса. В работе [5] предложена еще одна модификация метода, аппроксимирующая уравнение в консервативном виде - характеристический МИП, который рассмотрен на примере системы уравнений мелкой воды. Метод характеристик заключается в решении уравнений переноса для инвариантов Римана. Далее значения исходных функций находятся по найденным значениям инвариантов. В предложенном методе использована аппроксимация МИП полиномами второго (для нахождения инвариантов) и третьего (для вычисления исходной функции) порядков.
Работа [6] посвящена приложению консервативной модификации МИП к решению многомерных гиперболических уравнений. Рассмотрена одномерная модификация метода, основанная на построении полинома второго порядка. Ее продолжение на многомерный случай получено с использованием метода дробных шагов. Представлены результаты расчетов двух-и трехмерных задач.
В работе [7] описывается модификация метода в двумерном случае на неоднородной сетке. Вводится понятие сетки «Соробан», которая позволяет сохранить точность метода. Соро-бан представляет собой совокупность прямых линий в направлении координаты у при каждом заданном х. На линиях в начальный момент времени фиксируются точки, далее отслеживается их движение. Число точек на различных линиях может варьироваться. В результате вид расчетной области напоминает счеты - по-японски соробан. Для нахождения на п-м временном слое точки, из которой выходит характеристика, проводится линейная интерполяция по каждому направлению с помощью одномерной формулировки МИП.
В работе [8] МИП сравнивается с кубической интерполяцией Лагранжа и сплайн-интерполяцией на примере одномерного гиперболического уравнения. В отличие от этих методов МИП показывает лучшие численные результаты при различных начальных профилях.
Для интерполяции он требует значений функции только в двух соседних точках и хорошо передает градиенты решений. Также в работе представлено приложение МИП для двумерного уравнения Власова, приведены результаты расчетов.
Работа [9] посвящена приложениям метода для задач гидродинамики. Рассматриваются следующие постановки: численное моделирование обрушения столба воды в случае ламинарного и турбулентного течений, возникновение и распространение боры. Представлена модификация метода, являющаяся объединенной процедурой для вычисления движения сжимаемой среды, состоящей из двух или более фаз. В этом случае основные уравнения движения среды решаются в два этапа: сначала решаются однородные уравнения переноса для плотности, скорости и давления, потом для полученных решений вводятся неоднородные члены. В книге также представлен код программы на языке С.
Хотя метод интерполяционного профиля в настоящее время получил широкое применение, остаются нерассмотренными многие вопросы, связанные с его свойствами. Авторы ука-
занных работ не проводили общий анализ метода, не были доказаны аппроксимация и устойчивость. Также для случая переменной скорости переноса не был точно сформулирован метод нахождения точки п-го временного слоя, из которой выходит характеристика, приходящая в данный узел (п + 1)-го слоя. В настоящей работе проведен оригинальный анализ МИП, основывающийся на аппарате дифференциального приближения [10], позволивший доказать устойчивость метода, определить порядок аппроксимации. Было найдено изменение фазы волны для МИП и проведено сравнение с аналогичным параметром противопоточной схемы и точного решения, которое показало низкую дисперсию рассматриваемого метода. Также предложена модификация метода, основанная на более аккуратном методе решения уравнения для характеристик гиперболического уравнения. Она заключается в решении дифференциального уравнения для характеристики с повышенной точностью. Рассмотрены и проанализированы аппроксимации первого и второго порядков.
В работе предложен собственный вариант модификации алгоритма, осуществляющий взаимодействие характеристик на разрыве. Приводятся результаты расчетов одномерных задач, имеющих разрывные решения. Также предложен способ более точного вычисления потоков для МИП в консервативном виде. Предлагаемая модификация обобщается на двумерные задачи согласно методу дробных шагов. Для двумерной консервативной задачи был применен метод сохранения первоначального профиля, помогающий избежать влияния численной диссипации на решение. Указанный метод заключается в замене исходной функции на ее тангенциальное преобразование в уравнении переноса.
При выполнении данной работы ее авторы создали собственный вариант программы реализации численного алгоритма [11; 12]. В дальнейшем метод предполагается использовать в численной модели движения гомогенной трехфазной среды.
Метод интерполяционного профиля для одномерных уравнений
Начально-краевые задачи для одномерных уравнений переноса
и их точные решения
При построении МИП существенно используется свойство решения уравнения переноса с однородной правой частью сохраняться вдоль его характеристики. Приведем формулы для точных решений начально-краевой задачи для линейного, квазилинейного уравнений переноса и уравнения с переменным коэффициентом. Эти формулы будут положены в основу МИП.
Линейное уравнение переноса с однородной правой частью
Рассматривается линейное уравнение переноса с положительным постоянным коэффициентом. В этом случае дополнительные условия задаются на левой и на нижней (начальные условия) границах области, и задача формулируется как
ut + aux = 0, a = const > 0, 0 < x < l, 0 < t < T, t x (1) u (0, t) = ^(t), u (x, 0) = ф( x), ц(0) = ф(0).
Характеристики уравнения переноса в задаче (1) задаются уравнением
dx
~Г = a, dt
а точное решение задачи есть
Гф(х - М), если t < х / а; и( х, t) = •{
[ц(, - х / а), если t > х / а.
Квазилинейное уравнение переноса с однородной правой частью
При предположении, что и > 0, начально-краевая задача для квазилинейного уравнения переноса ставится следующим образом:
и, + ии = 0, 0 < х < I, 0 < t < Т, t х (2) и (0, t) = ), и (х, 0) = ф( х), ц(0) = ф(0).
Точное решение задачи (2) имеет вид
Гф(х - ш,), если t < х / и; и( х, t) = ] (4)
[ц(, - х / и), если t > х / и.
Отличительной особенностью решения (3) задачи (2) является то, что этот способ становится непригодным в случае, если с линии задания начальных данных в точку (х, ^ приходит сразу несколько характеристик либо не приходит ни одной. Этот случай характеризуется наличием разрыва решения х = Х(0, поведение которого описывается уравнением
ёХ ик + иь
Л 2
где ик и иь - решения задачи (2) справа и слева от разрыва.
Уравнение переноса с переменным коэффициентом и однородной правой частью
Следует отдельно выделить случай начально-краевой задачи для уравнения переноса с переменным коэффициентом и однородной правой частью. Пусть для определенности коэффициент уравнения а(х,t) > 0 и задача будет выглядеть как
и, + а(х,,)ш = 0, 0 < х < I, 0 < t < Т,
^ 4 ' 7 ^ ' ' ' (4)
и (0, ^ = ц(0, и (х, 0) = ф( х), ц(0) = ф(0). Решение задачи (4), так же как и задач (1) и (2), будет переноситься вдоль характеристик
^ = а( х, 0, (5)
Л,
которые в данном случае не являются прямолинейными. Это осложняет реализацию МИП ввиду необходимости численного решения уравнения (5).
Уравнение переноса с неоднородной правой частью
Наконец, наиболее общий случай заключается в рассмотрении начально-краевой задачи для уравнения переноса с переменным коэффициентом и неоднородной правой частью:
и + а(х, ,)ш = /(х, 0, а(х, ^ > 0, 0 < х < I, 0 < t < Т,
х (6) и (0, ^ = ц(0, и (х, 0) = ф( х), ц(0) = ф(0).
Точное решение задачи (6) выписать затруднительно, но МИП может быть обобщен и на этот случай.
Неконсервативный метод интерполяционного профиля Линейное уравнение переноса с однородной правой частью Формулы МИП
Предположим, что на п-м временном слое заданы значения функции и пространственной производной в узлах равномерной сетки (х,. = xi_1 + И; 7 = 1,..., N; И = I / N}. Тогда решение на интервале [ х{ —1, х{ ] (для удобства назовем его 7-м интервалом) представляется в виде 7-го кубического полинома
и(х,) = Ъгп(х) = а7х3 + Ъ х2 + с х + , х е [х1 —1,х1 ], удовлетворяющего условиям согласования
Ц (х7_1) = и7—1 , ^ (х7) = и , ---= (их )1_1, ---= (их Л • (7)
ах ах
Из (7) следует система уравнений на коэффициенты полинома.
Значения функции и ее пространственной производной на следующем слое по времени переносятся с п-го слоя по характеристике уравнения переноса. Предположим, х - точка п-го слоя, из которой выходит характеристика, приходящая в узел (х7, 7п+1). Для случая положительной постоянной скорости имеем
х = х7 — аъ
Интервал [хк—1, хк ], содержащий х, определяется величиной числа Куранта г =
ах .
к = 7 — I, где I - целая часть г.
Продифференцировав исходное уравнение переноса (1) по переменной х, получим уравнение переноса для их:
(их X + а(их )х =
имеющее структуру исходного уравнения. Окончательно получаем формулы для функции и ее производной в узле (х1, 7п+1):
ип+1 = Ркп (х, — ах), (их Г = аРк" (Г ах) • (8)
ах
Для случая числа Куранта 0 < г < 1, положив в (8) к = 7, £ = —ах, получим
иГ1 = Ъп (х + £) = а7 (х7 + £)3 +Ъ (х + £)2 +с (х7 + £) + а • Сгруппировав члены по степеням перепишем правую часть последнего равенства в виде
и п+1 = аД3 + (3а7х7 + Ъ )£2 + (3а7х2 + 2Ъixi + с 7 )£, + (а7х3 + Ъ7х2 + с7х7 + )• Заметим, что выражения для свободного члена и коэффициента при £ в первой степени совпадают соответственно с определенными в (7) значениями функции и ее производной в узле сетки (,7п). Приняв во внимание это замечание и обозначив коэффициент при за Ъ*
ЪГ = 3ах + Ъ,
получим расчетную формулу метода
и"++ = аД3 + Ъ*£2 + (их )п £ + < • (9)
Для производной получим аналогично
(их) п+1 = 3а. £2 + 2Ъ,*^ + (их )п Коэффициенты а7, Ъ* найдены из системы (7) и имеют вид
а = (их )п + (их )п—1 — 2 < — и;—1 2(их )п + (их )п—1 — 3 и п — и;—1 (10)
а = И2 2 И3 , Ъ = И И2 • ( )
Рассмотрим случай числа Куранта г > 1. Положим для начала 1 < г < 2 . Тогда характеристика, приходящая в узел (х, ^п+1), выходит из 7 - 1 интервала п-го слоя по времени. Расчетная формула для функции принимает вид
и Г1 = Ъ—1( х — ах)
Обозначив £ = И — ах, получим
""+1 = КЛX - ах) = % "1(х,_1 + к - ах) = = а,-1 (х,-1 + + Ъ;_1 (х,-1 + + с,-1 (х,-1 + £) + ^-1. Раскроем в последнем выражении скобки и сгруппируем члены по степеням 2,:
"Г1 = а,-1^3 + (3а,-1 х,-1 + Ь,.^2 + (3а,-1 х,2- + 2Ь,- х,- + с,,^ +
+(а,-1 х3-1 + Ь,-1 х2-1 + с,- х,-1 + ^,-1). Аналогично случаю 0 < г < 1 выражения для свободного члена и коэффициента при 2, в первой степени совпадают соответственно с определенными в (7) значениями функции и ее производной в узле сетки (х,._15Г). Итоговая расчетная формула принимает вид
""+1 = + Ь,*_1^2 + ("х )"-1 ^ + «"-1,
где
Ь,- =3а,-1х,-1 + Ь-1 •
Для производной получим соответствующее выражение
("х )"+1 = 3а,-¿2 + 2Ь,^ + ("х )"-1. Для к -1 < г < к, где к - целое положительное число, используя те же выкладки, выводим расчетные формулы
""+1 = а,-к _1^3 + Ь*_к-£2 + ("х )"-к-1 ^ + ""-к-1, ("х )"+1 = 3а,-к-¿2 + 2Ь,-к-1^ + ("х )"-к-1,
где ^ =(к - 1)к - ах, ЬГ-к-1 =3а,-к-1х,-к-1 + Ь,-к-1, к =12 -
Анализ МИП
Порядок аппроксимации и устойчивость
Построим для МИП дифференциальное приближение (ДП) [7]. Для этого в расчетной формуле (9) заменим коэффициенты а{, Ь* их выражениями из (10):
"+1 ". =
/ / \" . / \" " " \ ("х ), + ("х )"-1 - 2 " " - ",-1
(-ах)3
2("х )" + ("х )"-1 - 3 ""-К
(-ах)2 + ("х)"^ + "" ,
и подставим сюда вместо сеточных функций достаточно гладкие функции. Далее, применим формулу Тейлора к "(xi,¿"+1), "(х,-1,^) и "х(х,-1,^) для представления их в узле (xi,^). В результате получим исходное уравнение с погрешностью аппроксимации МИП
2 2 3 2
х а х х а х
- аи = — "„
-".... --
2хх ttt
6 6
а2 хк2
--"
2 01-11-х
" -
а2хк3 ( 3ах
-11 "V +.
24 24 ^ к ) 120 60 ^ 2к
В погрешности аппроксимации производные по времени преобразуются в производные по пространству. В результате получается дифференциальное приближение МИП
Ч2|'1
" + а"х = —
а2 хк2
(1 - г)2К +
а2 хк3
(1 - г)2
- г \"х + .
(11)
24 30 ^2
Главный член в погрешности аппроксимации в ДП содержит производную четвертого порядка, что говорит об очень низком уровне диссипации построенного метода. Из (11) следует, что МИП имеет порядок аппроксимации О(х2 + к4) и для него выполняется необходимое условие устойчивости.
= 1 С ^ , * 2, (13)
Изменение фазы волны
Решение и( х, 7) дифференциального уравнения (1) можно разложить в ряд Фурье по гармоникам
и( х, 7) = ехр [— 7(ю7 — кх)]^ (12)
Эти гармоники являются решениями уравнения при частоте волны ю = ак, где к - волновое число.
Для дифференциального уравнения вида
^ д'и
+ аих =
7 х £2 ' дх'
спектральный образ оператора шага р представляет собой такой множитель, при умножении на который гармоники (12) получается гармоника на временном слое 7 + х, являющаяся решением уравнения (13):
и( х, 7 + х) = ри (х,7 )•
Легко найти вид р:
р
р = ехр —7гф + х^ (7к)' с1 ,
_ '=2 _
где ф = кИ .
За время т гармоника изменяется на величину р, а фаза волны (аргумент гармоники) изменяется на величину
ф = arg р = -гф + xlm
É (ik)'c •
1=2
Величину ф называют изменением фазы волны (ИФВ) в решении уравнения (14) за время т.
Для уравнения переноса спектральный образ оператора шага и изменение фазы волны соответственно принимают вид
р = exp [—Гф], ф = arg р = —гф^ (14)
Аналогично вводятся понятия гармоники, множителя перехода и изменения фазы волны для эволюционной разностной схемы. Например, для противопоточной схемы, аппроксимирующей уравнение (1) с первым порядком,
ип+х — и" и" — и",
-L-L + -= 0, (15)
х h
имеем
и" = р" exp [7ф/ ], иП+1 = pU"", р(ф) = 1 — r + r cos ф — 7r sin ф^
Тогда ИФВ противопоточной схемы:
ф r sin ф
ф = а^р = —--• (16)
1 — r + r cos ф
Получим ИФВ схемы (15), исходя из ее дифференциального приближения. ДП противопоточной схемы имеет вид
ah ,, ah2 ,, J1
и + аих =y(1—Ф» ——(1—r)12—r Iи- +■
тогда ИФВ ДП противопоточной схемы
Фh =—гф + -Т(1 — r)I1 — r |ф3• (17)
Для выведенного ранее ДП МИП (11) ИФВ запишется как
Ф" =—гф + ¿(1 — r)2Í1 — r |ф5• (18)
Рис. 1. Зависимости ИФВ от величины ф при числах: г = 1 (1), 0,7 (2), 0,49 (3), 0,2 (4); (-) - уравнение (14); (■) - схема (16); (X ) -ДП схемы (17); (о ) - ДП МИП (18)
ll(X,t) 1
u(x,t) 1
0 0.2 0.4 0.6 0.Е
0.2 0.4 0.6 0.Е
а
т = 0,05; h = 0,25
б
т = 0,125; h = 0,067
Рис. 2. Точное (-) и численное, полученное МИП (■, □), решения задачи (1), (19) при числах Куранта 0,2 (а) и 1,875 (б)
На рис. 1 представлены зависимости ИФВ уравнения переноса, ДП МИП, противопоточ-ной схемы первого порядка и ДП противопоточной схемы от числа ф при различных числах Куранта. Видно, что МИП хорошо сохраняет фазу решения дифференциального уравнения, а следовательно, обладает малой дисперсией.
Результаты расчетов Гладкое решение
Входные данные задачи (1) конкретизировались следующим образом:
a = 1, l = 1, T = 1, ) = - sin(2^), ф( x) = sin(2nx). (19)
Точным решением этой задачи является функция
u (x, t) = sin (2л( x -1)).
Для МИП также необходимо задать начальное и краевое распределения производной, которые находились из точного решения задачи.
На рис. 2 приведены сравнения точного решения задачи с решениями, полученными МИП при различных числах Куранта в момент времени ^ = 0,5. Черными квадратами отмечены значения численного решения в узлах сетки, по величинам в которых функции и ее производной строились полиномы. Каждый интервал между двумя соседними узлами сетки разбивался на 20 подынтервалов, значения полинома в точках этого дополнительного разбиения с шагом 0,05к отмечены светлыми квадратами. Видно, что метод обладает низкими диссипа-тивными свойствами и остается устойчивым при г > 1.
Разрывное решение
Входные данные задачи (1) в этом случае задавались как
П, если х = 0;
а = 1; I = 1; Т = 1; = 1; ф(х) = \ (20)
[0, если 0 < х < 1.
Точное решение этой задачи представляет собой движущуюся вправо с единичной скоростью ступеньку.
Начальное и краевое распределения производной искомой функции полагались нулевыми. Альтернативные варианты их задания не приводили к улучшению решения.
и(х,1)
и(х,1)
0.6 0.8
0.6 0.8
т = 0,025; к = 0,05
т = 0,1; к = 0,05
Рис. 3. Точное (—) и численные, полученные МИП (■), а также по схеме (21) (X ), решения задачи (1), (20) при числах Куранта 0,5 (а) и 2 (б)
а
На рис. 3 приведены сравнения точного решения задачи с решениями, полученными МИП и по разностной схеме, при различных числах Куранта в момент времени ^ = 0,5. Квадратами отмечены значения численного решения МИП в основных узлах сетки. Штрих-пунктирной линией нанесено численное решение по МИП, построенное в узлах дополнительного разбиения с шагом 0,05к. Крестиками отмечено решение, полученное по противопоточной схеме второго порядка аппроксимации как по времени, так и по пространству:
и"~1 — 4и" + 3ип+1 и— 4м""1",1 + 3и"+1
-1-1-1— + а-^-^1-]— = 0. (21)
2т 2к
Как видно, МИП обладает некоторыми дисперсионными свойствами, которые выражаются в малых осцилляциях у разрыва.
Уравнение переноса с переменным коэффициентом
и однородной правой частью
В случае уравнения переноса с переменным коэффициентом необходима модификация построенного для линейного уравнения МИП в двух его частях.
Во-первых, характеристики в этом случае становятся криволинейными и для отыскания точки х на "-м слое, из которой выходит такая характеристика, приходящая в узел (, , необходимо численно решать уравнение (5).
Во-вторых, усложняется уравнение переноса для "х и требуется модификация МИП в части определения ее на (" + 1)-м слое.
Численное решение уравнения на характеристики
Для нахождения точки х аппроксимируем дифференциальное уравнение (5), описывающее характеристики, методом Эйлера первого
х1 = х + ха С х,^ ) (22)
или второго
х1 = х + ха | х + 2 а(х, t"), Г + (23)
порядков. Каждое из этих равенств может быть представлено в виде обобщенного уравнения для координаты х точки на "-м слое, из которой выходит характеристика, приходящая в узел (х,, Г+1),
ф, (х) = 0, (24)
где
ф,, (х) = х1 - х - ха( х, t")
в случае метода Эйлера (22) и
( х х
ф,,(х) = х1 -х-хаI х + 2а(х,Г),Г + —
в случае обобщенного метода Эйлера (23). При решении уравнения (24) в зависимости от соотношения знаков значений ф,, (х,-1), ф,, (х,) и а(х, t) выделяются три случая, для каждого из которых отыскивается интервал, содержащий корень уравнения. Методом дихотомии корень определяется с заданной точностью.
Метод дробных шагов решения уравнения для "х
Уравнение для "х в случае уравнения переноса с переменным коэффициентом (4) принимает вид
("х X + а(х>t)("х )х + ах (х>t)"х = 0. (25)
Для сохранения идеологии МИП уравнение (25) решается методом дробных шагов. На первом дробном шаге из уравнения
( " х X + а (х, t)( й х ) х = 0
посредством МИП находится предварительное значение производной "х . На втором шаге решается уравнение
("х X =-ах (х't)"х (26)
при начальном условии
"х[= а х.
Уравнение (26) решается конечно-разностным методом
^ = -ах (х(, t")(" х), •
т
Уравнение переноса с переменным коэффициентом и неоднородной правой частью
В данном случае метод дробных шагов применяется для решения как уравнения для самой функции (6), так и для ее производной
(их X + а( х, t )(их) х + ах (х, t )их = /х (х, t).
Метод дробных шагов решения уравнения переноса с неоднородной правой частью
На первом дробном шаге из уравнения
и t + а( х, t )й х = 0
посредством МИП находится предварительное значение функции и . На втором шаге решается уравнение
Щ = I (х, t) (27)
при начальном условии
и\ „ = и.
Уравнение (27) решается по схеме
я+1
= I (х, ).
х
Метод дробных шагов решения уравнения переноса производной
На первом дробном шаге из уравнения
(и x)t + а( х, t )(и х) х = 0
посредством МИП находится предварительное значение производной их . На втором шаге решается уравнение
(их X =-ах (х, t )их + /х (28)
при начальном условии
их[= и х.
Уравнение (28) решается конечно-разностным методом
(ихГ - (их), а (х .)(_ ) + I(х+1,) -1(х-1,)
-= -ах (х, ^ )(и х), +-—-.
х 2п
Результаты расчетов
Уравнение переноса с переменным коэффициентом и однородной правой частью
Входные данные задачи (4) конкретизировались следующим образом: а(х) = 1 + 0,5зш I — I; 1 = 200; Т = 10; = 0;
Г1, если 40 < х < 60; (29)
ф( х) = \
[0, если 0 < х < 40,60 < х < 200. Начальное и краевое распределения производной их полагались нулевыми. На рис. 4 приведены решения задачи (4), (29), полученные МИП с применением метода Эйлера (22) (о) и обобщенного метода Эйлера (23) (▲), а также по схеме (15) (◊), в момент времени t = 100.
u(x,t)
u(x,t)
т = 0,1; h = 2
б
т = 1; h = 1
Рис. 4. Численные, полученные МИП (о , ▲), а также по схеме (15) (◊) решения задачи (4), (29) в момент времени t = 100; уравнение на характеристики решается простым (о ) и обобщенным (▲) методами Эйлера
u(x,t)
u(x,t)
0.2 0.4 0.6 0.8
0.2 0.4 0.6 0.8
а
t = 0,85; т = 0,025; h = 0,033; A = 0,05
б
t = 0,86; т = 0,0045; h = 0,005; A = 0,01
Рис. 5. Точное (-) и численные, полученные МИП (о , ▲), решения задачи (6), (30). Уравнение для характеристики решается простым (о ) и модифицированным (▲) методами Эйлера
Уравнение переноса с переменным коэффициентом и неоднородной правой частью
Входные данные задачи (6) конкретизировались следующим образом: а( х) = х2, I = 1, Т = 1, / (х, t) = Г 0,5 + х2 ^ 1
2 A cosh2 )
Точное решение задачи (6), (30) имеет вид
f
u (x, t) = 0,5 1 + tan h
x - xA + 0,5t
лл
A
(30)
(31)
Функции ^), ф(х), а также начальное и краевое распределения производной "х определяются из (31). При решении уравнения (28) производная /х вычисляется не по конечно-разностной формуле, а явно. При стремлении параметра А к нулю решение (31) стремится
а
к ступеньке, движущейся влево. Начальное положение ступеньки определяется вторым параметром х0. Далее х0 полагается равным значению 0,95.
На рис. 5 представлены сравнения точного решения задачи (6), (30) с решениями, полученными МИП при различных значениях параметра А и шагов сетки.
Решение на рис. 5, а получено в момент времени 7 = 0,85. Значение параметров: х0= 0,95 и А = 0,05. Решение на рис. 5, б получено в момент времени 7 = 0,86. Значение параметров: х0= 0,95 и А = 0,01. Видно, что решения, найденные при помощи аппроксимации уравнения (5) с первым и вторым порядками, мало отличаются друг от друга.
Консервативный метод интерполяционного профиля
Предложенный выше метод обладает несколькими существенными недостатками, затрудняющими его прямое применение для решения задач гидродинамики.
Отсутствие механизма взаимодействия характеристик на разрыве
в неконсервативном МИП
В неконсервативном МИП при определении интервала сетки, из которого выходит характеристика, приходящая в нужный узел нового слоя, анализируется скорость переноса. Если скорость в узле (, ^) положительная, то полагается, что искомый интервал расположен слева от этого узла. И наоборот, если скорость в узле (х1, Г) отрицательная, то искомый интервал расположен справа, т. е. при вычислении решения учитывается информация, находящаяся только по одну сторону от узла. Этот подход неверен в случае разрывных решений, когда возникают ситуации, аналогичные представленной на рис. 6. В этом случае в (7 - 1)-й узел на любом слое по времени будут переноситься только положительные значения решения, а в 7-й узел - отрицательные. Данный алгоритм не описывает механизм взаимодействия характеристик на разрыве решения. Это приводит к тому, что при разности значений модулей скорости в этих узлах разрывное решение не будет передвигаться. Значит, численное решение не будет соответствовать точному.
Проиллюстрируем этот недостаток на примере. Рассмотрим начально-краевую задачу для квазилинейного уравнения Бюргерса с однородной правой частью (2). Ее входные данные конкретизировались следующим образом:
Г-4 х + 3, если 0 < х < 1; I = 2; Т = 2; и (0,7) = 3; и (2,7) = -1; ф( х) = <| (32)
[-1, если 1 < х < 2.
Рис. 7. Точное (-) и численное (• — •) решения задачи (2), (32) в моменты времени Г = 0,0875 (1), Г = 0,25 (2), Г = 1,125 (3); т = 0,0125, к = 0,05
: 2 3' 3
\ V 1
0.5
1.5
Рис. 7. Внесение механизма взаимодействия характеристик на разрыве: т = 0,0125, к = 0,05. Точное (-) и численные (• — •) решения задачи (2), (32) в моменты времени Г = 0,0875 (1), Г = 0,25 (2), Г = 1,125 (3, 3')
2 X
Задача (2), (32) описывает возникновение в момент времени Г = 0,25 в точке х = 0,75 разрыва и его движение с единичной скоростью при Г > 0,25. Было положено нулевое распределение производных в начальный момент времени во всей области. Уравнение на характеристики
ёх . .
— = и( х, Г) &
решалось методом Эйлера
х1 = X + хи (X).
На рис. 7 представлены точное и численное решения в моменты времени, соответствующие гладкому профилю, возникновению разрыва и его движению в дифференциальной задаче.
Число Куранта не превышает единицу во всей расчетной области. Как видно, после появления разрыва в решении, полученном МИП, профиль прекращает движение (см. рис. 7, 2), тогда как точное решение продолжает «бежать» вправо (см. рис. 7, 3). Чтобы обойти эту проблему, авторы предложили модификацию алгоритма, позволяющую учитывать направление скорости переноса в соседних узлах. Рассмотрим данную модификацию на примере задачи (2), (32).
Внесение механизма взаимодействия характеристик на разрыве на примере квазилинейного уравнения
Модификация МИП в случае решения уравнения (2) заключается во введении дополнительного шага алгоритма перед определением точки х выхода характеристики, переносящей искомую функцию в узел следующего (п + 1)-го слоя по времени. На этом дополнительном шаге проводится пересчет функции и в целых узлах на п-м слое по формуле
и( хк, 7п) =
^к (Хк-1/2 ) + ^к+1 (Хк+1/2 ) 2 '
Отсутствие свойства консервативности
С использованием предложенной модификации были произведены расчеты задачи (2), (32). На рис. 8 представлены те же решения и те же моменты времени, что и на рис. 7.
Заметим, что полученное с использованием предложенной модификации численное решение движется после возникновения разрыва (см. рис. 8, 3'). Однако оно существенно отстает от точного решения, и со временем это отставание увеличивается. Данное свойство полученного решения является недостатком МИП и связано с тем, что рассматриваемый метод не обладает свойством консервативности. Следовательно, следующим шагом в модификации МИП станет построение консервативного метода.
Построение консервативного МИП
на основе дивергентного уравнения переноса
Возьмем одномерное уравнение переноса в дивергентной форме с однородной правой частью:
и{ + фх = 0, ф = ф(и). (33)
Предположим, что ф(и) - дифференцируемая функция. Обозначив
а(и) = Сф,
Си
получим уравнение переноса с переменным коэффициентом и однородной правой частью
и{ + а (и (х, 7)) их = 0. (34)
Поставим для (34) задачу аналогично (4) и рассмотрим два варианта метода интерполяционного профиля, обладающих консервативностью [1].
Полином четвертого порядка
Предположим, как и в неконсервативном случае, что на п-м временном слое заданы значения функции и пространственной производной в узлах равномерной сетки. Решение на интервале [ х{-1, х{ ] представим в виде /-го полинома четвертой степени
и (х, 7п) = ^п(х) = д.х4 + Дх3 + С,х2 + 0,х + Е,, х е[х1-1,х1], (35)
удовлетворяющего четырем условиям гладкого сопряжения на концах интервала
Р," (х-1) = ип-1, (х) = и., ^Щг-1 = (их)п-1, ^^ = (их)п-
ах ах
Пятое условие обеспечит выполнение на сетке дискретного аналога закона сохранения
х, с
| [ип -ип-1 ]ах = | (-ф, +ф{-1)ж, (36)
х,-1 С-1
соответствующего дивергентному уравнению (33). Тем самым будет достигнуто свойство консервативности метода.
Полагая
х х;
| (х)ёх = | ипёх,
х—1 х—1
получим из (36) пятое условие на коэффициенты полинома (35)
х, х, Гп Гп
| (х) ёх = | Гп—1'(х) ёх + | ф1—1 ёГ — | ф,
где интегралы по времени вычисляются по формуле трапеции
2
фП + ФП—1
П1
интеграл
| ф, ёГ = ■■ т,
| Гп—1(х) ёх
известен по предыдущему шагу.
Для обеспечения механизма взаимодействия характеристик на разрыве уравнение на характеристики
ёх = а (и (х, Г))
аппроксимируется методом Эйлера первого порядка
- аП 1 + аП
х. = х + —1-- т,
2
где х е [ хк—1, хА ].
Характеристическая скорость а в целочисленных узлах сетки рассчитывается по значениям функций ф(и) в серединах соответствующих интервалов:
ап = фк+1/2 —Фк—!/2 , (37)
ик+1/2 — ик —1/2
где фк±1/2 = ф(ик±1/2 )> ик±1/2 = Ек (хк±1/2 ) .
х
, —1
Рис. 9. Точное (-) и численное (• — •) решения задачи (2), (32) посредством консервативного МИП в моменты времени Г = 0,0875 (1), Г = 0,25 (2), Г = 1,125 (3, 3'); т = 0,0125, к = 0,05
В случае квазилинейного уравнения (2) модификация (37) позволяет точно описывать движение разрывного решения, так как полученное таким образом выражение для скорости
а" +1/2 )2 / 2 -1/2 )2 / 2 _ Ч+1/2 + ик_у2
к ~ - ~ 2
ик+1/2 ик-1/2 2
согласуется с соотношением на линии разрыва для точного решения (3) уравнения (2)
ёХ и1 - и,2 ив + ит
& 2(иК - ит) 2 где иК и ит - значения функции соответственно справа и слева от разрыва.
Результаты расчетов
На рис. 9 представлены решения задачи (2), (32) с помощью консервативного МИП.
Видно, что консервативный МИП точно передает скорость движения разрыва, при этом в полученном численном решении отсутствуют диссипация и осцилляции (см. рис. 9, 3).
Как было сказано во введении, применение полинома четвертой степени улучшает аппроксимацию, но усложняет вычисления при решении многомерных задач. Поэтому для решения двумерных задач предложен консервативный МИП, основанный на построении полиномов второго порядка. Описание этого варианта представлено в следующем пункте.
Полином второго порядка
Решение на интервале [ х1 х1 ] п-го временного слоя представляется в виде /-го полинома второй степени
и(х,) _ (х) _ д.х2 + Вгх + Сг, х е [х1 х1 ], (38)
удовлетворяющего двум условиям согласования в узлах
V"(х,-1) _ и"-1, (х,) _ и".
Третье условие обеспечивает выполнение на сетке дискретного аналога закона сохранения
| (х) ёх _ | Ъ"-1(х) ёх + | (-ф, +ф,-1) Л.
Сравнение результатов применения полиномов второго и четвертого порядков
Линейное уравнение переноса с однородной правой частью
Сравнение было проведено на примере задачи (1). Входные данные конкретизировались следующим образом:
Г1, если 0,1 < х < 0,3;
а(х) _ 1; I _ 1; Т _ 1; *) _ 0; ф(х) _<| (39)
[0, если 0 < х < 0,1; 0,3 < х < 1.
Начальное и краевое распределения производной их полагались нулевыми.
На рис. 10 приведено сравнение точного решения задачи с решениями, полученными двумя вариантами консервативного МИП, в момент времени * = 0,6. Решения получены при шагах сетки И = 0,0167, т = 0,0067 и числе Куранта г = 0,4. Решение МИП, основанное на построении полиномов четвертого порядка, имеет большие осцилляции.
Метод интерполяционного профиля для многомерных уравнений
Метод интерполяционного профиля для двумерного уравнения переноса
В данном разделе консервативный МИП обобщается на двумерное уравнение переноса. За его основу принимаются полиномы второго порядка.
u(x,t)
Рис. 10. Точное (—) и численные решения задачи (1), (39), полученные консервативным МИП с применением полиномов (35) (X) и (37) ( о ): 7 = 0,6, т = 0,0067, к = 0,0167
Двумерное уравнение переноса в дивергентной форме с однородной правой частью можно представить в виде
u +Фх + Vy =0 Ф= Ф(й X V = V(uX (40)
0 < x<l, 0 < y <l, 0 < t <T. ( )
x' y'
Пусть на n-м временном слое заданы значения функции и пространственных производных в узлах равномерной сетки
{(x,,У,-)l x, = ^_ + hx;y- = y-_ + hy;i = 1,...,Жх; j = 1,...,Жу;hx = lx /Жх;hy = ly /Жу}.
Обобщение МИП проводится методом дробных шагов. На каждом дробном шаге решаются одномерные уравнения для функции и ее производной.
Метод дробных шагов решения двумерного уравнения переноса
На первом дробном шаге решаются одномерные уравнения вдоль линий x = const
й t + V y = 0. (41)
Посредством консервативного МИП на промежуточном слое по времени t вычисляются
предварительные значения функции й в узлах сетки и значения интеграла
y-
| й dy
y-_i
на каждом j-м интервале разбиения линий x = const. На промежуточный слой по времени требуется также перенести значения интеграла
x,
J й dx. (42)
xi_i
Для этого проинтегрируем уравнение (41) по переменной x на интервале x е [xi_i, x, ]
xi
J [ut + Vy] dx = 0, xi _1
после перестановки знаков интеграла и производной получим уравнение
d J й dx d J V dx
■ = 0, (43)
dt dy
являющееся уравнением переноса интеграла (42).
Решая уравнение (43) посредством консервативного МИП, получим на промежуточном слое по времени значения интегралов (42) на каждом интервале сетки [x._1?x.], а также значения двойного интеграла
Уj xi
II и dxdy.
yj_i xi_i
На втором дробном шаге решаются уравнения вдоль линий y = const
U +Ф x = 0
при начальном условии
и\ . = и.
It=t
Определяются значения функции и"+1 в узлах сетки и значения интегралов
xi
I и"+1 dx
xi_1
на каждом i-м интервале разбиения линий y = const.
Далее аналогично первому дробному шагу необходимо для каждого интервала [y _1? y ]
вычислить значения интеграла
yj
I и"+1 dy. (44)
yj_1
Для этого решается уравнение переноса
y, y,
д | udy д | ф dy
■ + -= 0
д* дх
посредством консервативного МИП. Получим на (" + 1) слое по времени значения интеграла (44), а также значения двойного интеграла
У] х1
| | и"+1 ёхёу.
У]-1 х,-1
Результаты расчетов
Входные данные задачи (40) конкретизировались следующим образом:
ф_и, ф_и, 1х _1, 1у _1, Т _1. (45)
Начальное распределение для этой задачи представлено на рис. 11. Функция и во всей расчетной области равна нулю, кроме квадрата, расположенного в левом нижнем углу расчетной области, где функция равна 1. Ширина стороны квадрата - 0,25.
^ г-
0.8 -
0.6 -0.4 -
0.2 -
Рис. 11. Начальное распределение для задачи (40), (45) Г...................
0 0.2 0.4 0.6 0.8 1
X
Задача (40), (45) описывает движение квадрата вверх по диагонали расчетной области со скоростью 1,41.
На рис. 12 представлены результаты расчета задачи консервативным МИП для двумерного уравнения переноса.
Рис. 12. Численное решение задачи о переносе квадрата (40), (45), полученное МИП
Решению на рис. 12 соответствуют разбиения кх = 0,0167, ку = 0,0167, т = 0,0067, момент времени г = 0,6. Видно, что применение МИП приводит к некоторому размазыванию численного решения. Поэтому для повышения разрешающей способности метода применяется подход, заключающийся в решении уравнения переноса образа тангенциального преобразования искомой функции вместо исходного. Полученное таким образом численное решение сохраняет резкость начального профиля.
Тангенциальное преобразование
Положим, есть функция, принимающая значения 0 и 1 в некоторой области й:
1, если (х, у, г) е О;
и ( х, г) =
0, если (х, у, г) £ О.
Назовем тангенциальным преобразованием / (х, у, г) следующую функцию:
Е(/) = [(1 -в)ж(/ - 0,5)], (46)
где 8 - малое положительное число.
Множитель 1 - 8 позволяет получать значения Е в окрестности -да при/ = 0 и в окрестности при / = 1. Обратное преобразование есть
/ = ^Е (/) ^ + 0,5. (47)
(1 -в)я
Допустим, функция / (х, у, г) удовлетворяет уравнению переноса
/ + аV/ = 0,
тогда для Е справедливо аналогичное уравнение переноса
+ а ^ (/) = 0.
ОТ
Следовательно, для нахождения Е (/) можно применить метод интерполяционного профиля. После вычисления Е (/) значения функции / определяются по формуле (47).
Результаты применения тангенциального преобразования для одномерного уравнения
Для начала данный подход был применен для решения одномерного линейного уравнения переноса (1) с входными данными (39). Решение задачи (1), (39) проводилось консервативным МИП, основанным на построении полиномов второго порядка. Однако вместо уравнения для искомой функции решалось уравнение переноса ее образа тангенциального преобразования.
0.6 0.7 0.8 0.9
8 = 0,7
0.8 0.9
8 = 0,2
>оооо<£
0.8 0.9
8 = 0,05
Рис. 13. Численные решения задачи (1), (39): * = 0,6, т = 0,0067, И = 0,0167. Результат применения тангенциального преобразования для одномерного консервативного МИП
На рис. 13 представлены сравнения точного решения задачи (1), (39) с численными решениями при различных значениях параметра 8. Отметим, что при малом значении параметра 8 численное решение хорошо аппроксимирует точное.
Результаты применения тангенциального преобразования для двумерного уравнения
Вернемся к задаче (40), (45) с начальным распределением, представленным на рис. 11
Рис. 14. Численное решение, полученного консервативным МИП с применением тангенциального преобразования 8 = 0,05
1
На рис. 14 приведены численные решения задачи (40), (45) консервативным МИП с тангенциальным преобразованием (46) при значении параметра 8 = 0,05 и решения, полученного по ТУЭ-схеме. Заметим, что разрешающая способность предложенного варианта МИП близка к ТУЭ-схеме, однако МИП обладает преимуществом, заключающимся в простоте его реализации. Напомним, что исследуемый метод является явным и абсолютно устойчивым.
Заключение
В данной работе представлен метод интерполяционного профиля, предложена его консервативная модификация и осуществлено обобщение на многомерный случай. Для одномерно-
го случая предложен метод решения уравнения для характеристик гиперболического уравнения. Проведен априорный сравнительный анализ свойств решения дифференциального уравнения и численных решений, полученных МИП и с помощью разностных схем. Установлено, что метод обладает низкими диссипацией и дисперсией, высоким порядком аппроксимации. Разработан механизм описания взаимодействия характеристик на разрыве решения. Проведены численные расчеты линейных, квазилинейных, нелинейных задач. Решения, полученные МИП и по разностной схеме, сравнены с точными решениями дифференциальных задач. С использованием метода дробных шагов проведено обобщение МИП на многомерные гиперболические задачи. Предложен подход, повышающий разрешающую способность метода. Проведены численные расчеты задачи двумерного переноса. Полученные результаты дают возможность рекомендовать метод интерполяционного профиля для решения многомерных задач течения многофазных сред.
Список литературы
1. Takewaki H., Nishiguchi A., Yabe T. The Cubic-Interpolated Pseudo-Particle (CIP) Method for Solving Hyperbolic-Type Equations // J. Comput. Phys. 1985. Vol. 61. P. 261.
2. Yabe T., Xiao F., Utsumi T. Constrained Interpolation Profile Method for Multiphase Analysis // J. Comput. Phys. 2001. Vol. 169. P. 556-593.
3. Hu C. An Introduction to CFD. Lecture Note for M5351 Dynamics of Ocean Systems I / Research Institute for Applied Mechanics. Kyushu University, Japan, 2011. 44 p.
4. Tanaka R., Nakamura T., Yabe T. Exactly Conservative Semi-Lagrangian Scheme (CIP-CSL) in One-Dimension // NIFS-685. 2001. P. 1-12.
5. Toda K., Ogata Y., Yabe T. Multi-Dimensional Conservative Semi-Lagrangian Method of Characteristics CIP for the Shallow Water Equations // J. Comput. Phys. 2009. Vol. 228. P. 49174944.
6. Nakamura T., Tanaka R., Yabe T., Takizawa K. Exactly Conservative Semi-Lagrangian Scheme for Multi-Dimensional Hyperbolic Equations with Directional Splitting Technique // J. Comput. Phys. 2001. Vol. 174. P. 171-207.
7. Yabe T., Takizawa K., Chino M., Imai M. and Chu C. C. Challenge of CIP as a Universal Solver for Solid, Liquid and Gas // Int. J. Numer. Meth. Fluids. 2005. Vol. 47. P. 655-676.
8. Lesur M., Idomura Y., Tokuda S. Kinetic Simulations of Electrostatic Plasma Waves Using Cubic-Interpolated-Propagation Scheme // JAEA-Research-089. 2006.
9. Furuyama S., Chanson H. A Numerical Study of Open Channel Flow Hydrodynamics and Turbulence of the Tidal Bore and Dam-Break Flows / Department of Civil Engineering. The University of Queensland, 2008.
10. Шокин Ю. И., Яненко Н. Н. Метод дифференциального приближения. Применение к газовой динамике. Новосибирск: Наука, 1985. 364 с.
11. Жамбалова Д. Б. Анализ метода ограниченного интерполяционного профиля на основе его дифференциального приближения // Материалы XLVII Междунар. науч. студ. конф. «Студент и научно-технический прогресс». Новосибирск, 2009. С. 259-260.
12. Жамбалова Д. Б. Обобщение метода ограниченного интерполяционного профиля на двумерные задачи с уравнениями в дивергентном виде // Материалы XLVIII Междунар. науч. студ. конф. «Студент и научно-технический прогресс». Новосибирск, 2010. С. 248.
Материал поступил в редколлегию 19.07.2011
D. B. Zhambalova, S. G. Cherny METHOD OF INTERPOLATION PROFILE FOR SOLVING TRANSPORT EQUATIONS
Method of interpolation profile that combines high order of approximation and consideration of structure of hyperbolic equations' solutions is presented. To increase the resolution of the method the special approach is applied. The main idea of this approach is to solve a transport equation for tangential transformation of initial function. The results of computations of one and two dimensional problems including having discontinuous solutions are given.
Keywords: constrained interpolation profile method, hyperbolic equations, differential approximation, method of characteristics.
