Об устойчивости схем расщепления и приближенной факторизации для решения систем многомерных уравнений Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 16, № 6, 2011

Об устойчивости схем расщепления и приближенной факторизации для решения систем многомерных уравнений*

В. М. Ковеня

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: kovenya@ict.nsc.ru

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

Ключевые слова: уравнения Эйлера, Навье — Стокса, разностная схема, методы расщепления, метод факторизации.

Введение

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

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

* Работа выполнена при финансовой поддержке РФФИ (грант № 11-01-00294-а) и Интеграционных проектов СО РАН № 44, 103.

числа арифметических операций на узел сетки [9], так как реализация схем сводится к матричным прогонкам, требующим обращения матриц размером (mJl) х (mJ2) х (mJ3), где Jj — число узлов в направлении ж (г = 1, 2, 3), Очевидно, что такой подход является неэкономичным уже в двумерном случае. Поэтому начиная с 60-х годов XX века для решения многомерных задач были предложены различные классы экономичных алгоритмов: метод расщепления, аддитивные схемы и схемы суммарной аппроксимации, схемы приближенной факторизации, метод предиктор-корректор и т.д. (см, [1-11]), Основная идея при их построении состоит в сведении исходной многомерной задачи (на дифференциальном или на разностном уровне) к совокупности более простых задач, например, к последовательности их одномерных аналогов, решение которых может быть получено экономичными алгоритмами.

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

1. Базовые схемы для эволюционных уравнений

Иллюстрацию базовых разностных схем для численного решения многомерных уравнений, разрешенных относительно производной по времени, проведем для системы линейных уравнений гиперболического (Б3- = 0) или параболического типа

ди _^ д д2

ж + £и = 0, 1 = (1)

3 = 1 ^ 3

где А3-, Бз — квадратные матрицы раз мерности т х тс постоянными коэффициентами, которые будем полагать перестановочными, N — размерность задачи по пространству. Аппроксимируем операторы Ь3- разностными операторами Ь^ с порядком 0(к1к). Разностная схема с весами

ип+1 - ип

_

+ Ьн(аип+1 + вип) = 0, Ьн = ЬзН, в =1 - а, (2)

3 = 1

или эквивалентная ей схема в каноническом виде

иП+1— иП

(/ + таЬн)-+ Ьнип = 0, (3)

аппроксимирует исходные уравнения (1) с порядком 0(т2 + Ьк) при а = 0.5 + 0(т) и при а > 0.5 безусловно устойчива (см., например, [2-4]), Здесь т, — шаги сетки по времени и пространству соответственно. Схема (3) реализуется матричными прогонками и неэкономична. Для получения экономичных схем рассмотрим схему расщепления. Стандартная схема расщепления с весами [1] или локально одномерная схема [2]

(ип+з/_ - ип+(з-1)/_)/т + ЬзН(аип+з/_ + вип+(з-1)/_) = 0, з = 1,..., N (4)

аппроксимируют исходные уравнения с порядком 0(т2 + Ьк) при а = 0.5 + 0(т) и при а > 0.5 безусловно устойчивы [1, 2]. На каждом дробном шаге система разностных

уравнений может быть решена векторными прогонками, требующими приблизительно 0(Л ■ ■ ■ х ) операций при переходе с п-го па (п + 1)-й слой (ем, [9]), После исключения дробных шагов она может быть записана в виде

— м™

ТТ (/ + таЬ¿н)-+ Ьнип = -ф, V = т(а2 - [32)(Ь1НЬ2Н + ... + Ьи_шЬин)ип+

т

7=1

+т2 (а3 + + ... + ¿м-1Л ¿млК + 0(т3). (5)

Таким образом, схема расщепления (4) приводится к каноническому виду лишь в двумерном случае при а = 0.5. Следовательно, стационарное решение задачи (1) при N > 2, полученное по схеме (4) методом установления, зависит от значения итерацион-т

Схема приближенной факторизации (см. [1, 5, 11])

с = —

(I + та^ )С+1/М = С, (I + та^ )£П+1 = £-1)/М,

мга+1 = + т£га+1

аппроксимирует исходные уравнения (1) также с порядком 0(т2 + ) при а = 0.5 + 0(т) и реализуется па дробных шагах подобно схеме расщепления. После исключения дробных шагов она может быть представлена в каноническом виде

П (/ + таЗД---+ Ьнип = 0. (6)

.7 = 1

Схема стабилизирующей поправки [1]

(мп+1/м — мп)/т + — мп) = — ^

(мп+2/м — мп+1/м)/т + ^(«га+2/М — мп) = 0,

(мп+1/м — -1)/м )/т + ¿^(«га+1/М — мп) = о

а=1

что и схема приближенной факторизации. Подобные свойства имеет и схема предиктор-корректор [1]

(мп+1/2М — мп)/та + = 0,

/2м — -1)/2М)/ та + ¿м^мга+2/2М = 0, (мга+1 — мп)/т + ^мга+1/2 = 0,

так как она также приводится к виду (6). Таким образом, схемы приближенной факторизации, стабилизирующей поправки и предиктор-корректор эквивалентны в целых

шагах и, следовательно, имеют одинаковые свойства. Как следствие, разностные схемы, приводимые к каноническому виду, обладают полной аппроксимацией и пригодны для решения как нестационарных, так и стационарных задач. Для гиперболических уравнений они безусловно устойчивы лишь для N < 2 при а > 0.5 и в трехмерном

а

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

т

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

2. Базовые неявные схемы для решения уравнений Навье — Стокса несжимаемой жидкости

Традиционно уравнения Навье —Стокса вязкой несжимаемой жидкости решаются в переменных функция тока—вихрь, В этом случае система уравнений сводится к параболической системе уравнений для компонент вихря и к эллиптической системе для векторного потенциала (в двумерном случае для функции тока и вихря). При решении пространственных задач использование естественных физических переменных может оказаться предпочтительнее. Для численного решения уравнений Навье —Стокса в естественных переменных часто применяются различные численные методы [1, 5, 8-17], Большая часть этих методов использует предположение, что на промежуточном слое течение предполагается солепоидальпым, а на полном шаге решается уравнение для давления, т, е, в алгоритме используется идеология расщепления. Альтернативой этим подходам является метод искусственной сжимаемости [1, 13], состоящий во введении в уравнение неразрывности производной по времени с малым коэффициентом. Получающаяся система уравнений является классической — разрешенной относительно производной по времени, но для получения решения нестационарных задач необходимо обеспечить сходимость решения модифицированных уравнений к решению исходных уравнений на каждом временном слое [1, 13], Такой подход может оказаться неэкономичным, Ниже приводятся описание и анализ классических неявных разностных схем для решения уравнений Навье — Стокса в физических переменных, приводимых к канонической форме. Некоторые схемы, основанные на методе суммарной аппроксимации, и аддитивные схемы с расщеплением многомерных операторов рассмотрены в работах [9, 16].

В качестве исходной модели будем рассматривать уравнения Навье — Стокса вязкой

несжимаемой жидкости в декартовых координатах без учета внешних сил

з

МЖ = -УГ, W = (7)

3 = 1 ^

где

( Р \

VI

V ^ )

(

Wj

\

VIV + , - тд ^ + , - Т,2

V ^ + , - , У

м

0 0 0 0

0 1 0 0

0 0 1 0

0 0 0 1

т,

л

ди- дг>,

дя, дя

Здесь р, V- — давление и компоненты вектор а скорости V ^ _ коэффициент вязкости, которую будем полагать постоянной, Плотность р принята постоянной и исключена

м

формой (7) приведем уравнения Навье — Стокса в недивергентном виде

0, В

Е в-,

-=0

(8)

где

В

о —

0 д1 д2 д3 \ 0 0 0 0 \

д1 0 0 0 В, — 0 6, 0 0

д2 0 0 0 0 0 6, 0

\ дз 0 0 0 / 0 0 0 6, /

д,

А

6!

.

В

физическим процессам и пространственным направлениям, где матричные операторы В, — 1, 2, 3) учитывают конвективные и вязкие члены в уравнениях движения по направлению я,, а оператор В0 содержит члены в уравнении неразрывности и члены с давлением в уравнениях движения. Отметим эквивалентность представления уравнений в дивергентной (7) и недивергентной (8) формах, т, е, выполнение соотношения

W — В^..

Аппроксимируем операторы В, разностными операторами В,^ по формулам

(9)

В

о^

/О

Аш Лгл

V Л3Ь

0 0 0

Л2^ 0 0 0

Лз^ \ 0

0 0

В

0 0

0 0

0

Ы 0 0

0 0

Ы 0

0 0 0

, у

Здесь , — г^Л,-

Л:,+, Л.,- = Л.,_ при V] < 0, Aj

Л,

в случае аппроксимации производных с первым порядком О (Л,) и

Л, = Л7 = (Л,-+Л,+)/2

Л,+ при V, > 0

при аппроксимации производных со вторым порядком О(Л2), где Л,± — — Т±)/Л,, — /г±1 — операторы сдвига, а вторые производные в В,^ аппроксимированы операторами Л,, со вторым порядком О (Л2).

Рассмотрим разностную схему с весами

сп+1 _ сп 3

М + та^2 в3н ) ---+ та ^ ВзН(аГ+1 + /3) Г = О

3=0 ) т 3=0

и перепишем ее в каноническом виде (с учетом соотношений (9))

сп+1 _ сп 3 3

с---= -Л¥™ = -С = М + та^В3н. (Ю)

3=0 3=0

Данная схема линейна относительно (п + 1)-го слоя и аппроксимирует уравнения (7) с порядком 0(т2 + Ьк) при а = 0.5 + 0(т), к =1, 2 в зависимости от порядка аппроксимации операторов 3 Реализация схемы без введения итераций (для тождественного выполнения разностного уравнения неразрывности) невозможна, В силу вырожденности стабилизирующего оператора С его классическая факторизация, т. е, представление в виде произведения, например одномерных операторов, невозможна, поэтому в работах [5, 14, 15] были предложены явно-неявные схемы расщепления вида

т

п+1 ~п

-—— + v/lpra+1 = О, <Иу^п+1 = О, (И)

т

аппроксимирующие исходные уравнения (8) с порядком 0(т + Ьк) и условно устойчивые в силу явной аппроксимации уравнений движения на промежуточном шаге. На полном шаге давление рп+1 находилось из условия

vn+1 _ V = _ту^(рп+1 _ рп)

которое с учетом уравнения vn+1 = 0 приводится к уравнению Пуассона

&УнЧн(рп+1-рп) = - с11улу. (12)

т

Решение (12) итерационным методом дает новое значение давления, после чего из (11) явно определяются значения vn+1, Легко видеть, что схема (11) не приводится к каноническому виду. Для повышения устойчивости схем в работах [5, 10] предложена неявная схема расщепления (например с порядком 0(т))

V _ V"- х3-«, _

т 3=1

_ ^

+ 22В^ + Ъ1грп+1 = 0, сИу^га+1 = 0, (13)

' 3 = 1

реализация которой при вычислении V сводилась к матричным прогонкам. После исключения V (13) можно представить в каноническом виде

сп+1 _ сп 3

/ + т (М + тВон)-= - ВзНГ = (14)

3=1 / т 3=0

Рассмотрим устойчивость схем. Анализ устойчивости будем проводить спектральным методом для линеаризованных уравнений — уравнений с "замороженными" коэффициентами, полагая значения коэффициентов постоянными [3]. Пусть

т

IJ1J2J3

fo exp + i = An fo exp i ^ kjih

i=i

i=i

вариация функции f0, Здесь

пт, яг — ^ Л — ешга, где Л — коэффициент роста = тах | £„

справедлива оценка ||Р"|| < |Л| операторы Aj, Aj аппроксимируют первые производные с порядком О (К) или О (к2). Тогда при аппроксимации с первым порядком

гармоник. Введем сеточную норму ||fn|| = max |j2j31 Для всех 0 < j < Ji; для которой

|fn-1| < ... < |A|n • |f0|. Будем полагать, что разностные

т Л, fn

где, например, при v > 0

dj fn,

тА3Г

d fn

Л^- fn = _r2,

^ = г(1 - /кз, ф = т(е*к^ - 1)^ = vjdj + г2, * = 1 + ¿0,

¿о = к + ¿2 + ¿з, ¿¿Д- = -г2, г2 = А(т2/Ь2) зт2(к^/2). (15)

При аппроксимации со вторым порядком Л.,- = Aj = (Aj_ + / 2 и tj = = й^-, — (тг,/Л,) 8т(к,Л,), ¿0 — ¿1 + ¿2 + ¿3, ¿0 — ¿¿0, Разностной схеме (14) удовлетворяет характеристическое уравнение

det

0 d1A d2A d3A

di[(l+i0)A-i0] (1 + t0)X — 1 0 0

d2[(l+io)A-i0] 0 (1 + to)X — 1 0

d3[(l+io)A-i0] 0 0 (l-Ko)A-l

корни которого равны Л1 — 0, Л2 — ¿0/(1 + ¿0), Л3,4 — 1/(1 + ¿0) и, очевидно, |Л,| < 1, т, е, схема (14) устойчива. Заметим, что реализация этой схемы при вычислении компонент вектора V возможна матричными прогонками, что делает данный подход неэкономичным.

Еще одним подходом при вычислении компонент скорости могут быть схемы расщепления или стабилизирующей поправки [10, 17]

(^/6 - т + ВшК+1/6 - V") — - К] В,-Л + У^

\,=1

(^+2/6 - ^+1/6)/т + В2^"+2/6 - V") — 0, (^+3/6 - ^+2/6)/т + В2^га+3/6 - V") — 0, (^ - ^+3/6)/т + у^(рга+1 - р") — 0, — 0.

Их реализация на дробных шагах сводится к скалярным прогонкам для вычисления компонент вектора скорости и решению уравнения Пуассона вида (12) для нахождения давления. После исключения дробных шагов схема приводится к канонической форме

П(1 + TBjh)(M + TBoh)

fn+1 _ fn

j=1

т

5>jfn

(16)

j=0

Ей соответствует характеристическое уравнение 0 di Л d2A

det

da [T(Л - 1) + 1] 0

0

da Л

di [T(A — 1) + 1] T(A — 1) + i0 0 0

d2[T(Л - 1) + 1] 0 T(A — 1) + t0 0

T(Л - 1)+ to

корни которого Л1 = 0, Л2 = (T — 1)/T, Л3,4 = (T — to)/T, а

Т = Д (1 + ^) = 1 + ¿0 + Я, ¿0 = ¿1 + ¿2 + ¿3, Я = ¿1*2 + ¿1*3 + ¿2*3 + ¿1*2*3 (СМ. (15)). ¿=1

При отсутствии конвективных членов в уравнениях движения (при tj = г2) получим

А2 = (Т - 1)/Т = (¿0 + Я)/(1 + ¿0 + Я), А3,4 = (Т - ¿о)/Т = (1 + Я)/(1 + ¿о + Я),

и, очевидно, | Аj| < 1 для ] = 1, 2, 3, 4, т. е. схема (16) безусловно устойчива.

Рассмотрим второй предельный случай при г2 = 0, соответствующий уравнениям Эйлера, Тогда, например, для схемы второго порядка получим

^ = г tj, Т =1 - Я1 + г(*0 - Я2), Я1 = *1*2 + ¿1*3 + ¿2*3, Я2 = ¿1*2*3

и

Л2

—Ri — i(to — R2)

|Л2

Ri + (to — R2)2

1 — Ri — i(to — R2) (1 — Ri)2 + (to — R2)2

Отсюда следует оценка т < уД/3 • h/|v|, где h = min |hj| , |v оставшихся корней характеристического уравнения получим

1 — Ri — ¿R2

< 1 при |Ri| < 0.5.

max |vj Для двух

Л

3,4

1 — Ri — i(to — R2)

и |А3,4| < 1 в двумерном случае, поскольку, как для всякой схемы приближенной факторизации (см, [6, 11]), Я2 = 0 и |А3,4| < 1 лишь при т < 4^/|3^|, Подобные оценки справедливы и для схемы первого порядка аппроксимации. Таким образом, схема приближенной факторизации вида (16) для численного решения уравнений Навье —Стокса несжимаемой жидкости (7) или (8) теряет свойство безусловной устойчивости уже в двумерном случае. Отметим, что для систем эволюционных уравнений схемы приближенной факторизации для двумерных задач безусловно устойчивы. Можно ожидать, что полученные оценки устойчивости для линейных уравнений будут справедливы и для нелинейной системы уравнений Навье —Стокса (7), Отметим, что для модели Стокса, в которой отсутствуют конвективные члены в уравнениях движения, эти схемы безусловно устойчивы для задач любой размерности,

В работе [18] была предложена схема приближенной факторизации

= -Wra,

I + та £ Bjhj C+i/2 = С (M + таВ^ )£n+i = £n+i/2,

fn+i = fn + т£ n+i

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

М + та ^ В.л « (/ + та ^ В.л ) (М + таВдл) = М + та ^ В.л + т2а2 (^ В,-л ) Вол.

¿=0 V ¿=1 / ,=0 \,=1 /

После исключения дробных шагов она приводится к канонической форме (14), т, е, безусловно устойчива. Решение разностных уравнений на первом дробном шаге п +1/2 предлагалось находить скалярными прогонками по одному из пространственных направлений для каждой компоненты еп+1/2 с введением внутренних итераций. Проиллюстрируем реализацию схемы (17) на двумерной задаче. Система разностных уравнений, например, для первой компоненты скорости может быть записана на пятиточечном шаблоне в виде

веп+1/2 = еп+1/2 + , еп+1/2 + с еп+1/2 + , еп+1/2 + е еп+1/2 = . В?,,г = ?,-1,г + Ч? + Ч? ?,+1,г + ?,,г-1 + ?,,г+1 = Л,г.

Для решения этих разностных уравнений рассмотрим итерационную схему

А(е "+1 - ев ) = вГ,

где еУ = а оператор А выберем из условий наиболее быстрой сходимости.

Подобным образом могут быть реализованы разностные уравнения и для остальных компонент скорости, В [18] для вычислений еп+1/2 были использованы методы простой итерации и Зейделя, При расчете течений в каверне для достижения сходимости требовалось несколько внутренних итераций. После определения невязок скоростей по

итерационной схеме численно решалось уравнение Пуассона

3 1

<цул = = ; ^ V

.=1 т

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

3 е *+1_е ^ 3 / 1 \

П С - = Е М - .

,=1 ¿=1 4 7

реализуемой скалярными прогонками по каждому направлению. После вычисления ерп+1/2 явно вычислялись невязки скоростей, затем из последнего шага схемы (17) явно находились новые значения функций на новом временном шаге.

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

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

[1] Яненко H.H. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.

[2] Самарский A.A. Теория разностных схем. М.: Наука, 1989.

[3] Годунов С.К., Рябенький B.C. Введение в теорию разностных схем. М.: Наука, 1962.

[4] Марчук Г.И. Методы расщепления. М.: Наука, 1989.

[5] Белоцерковский О.М. Численное моделирование в механике сплошных сред. М.: Наука, 1984.

[6] Ковеня В.М., Яненко H.H. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, 1979.

[7] Вабищевич П.Н., Самарский A.A. Аддитивные схемы для задач математической физики. М.: Наука, 1999.

[8] Роуч П. Вычислительная гидродинамика. М.: Мир, 1980.

[9] Самарский A.A., Николаев Е.С. Методы решения сеточных уравнений. М.: Наука, 1978.

[10] Толстых А.И. Компактные разностные схемы и их применение в задачах гидродинамики. М.: Наука, 1996.

[11] Ковеня В.М. Разностные методы решения многомерных задач: Курс лекций. Новосибирск: Новосибирский гос. ун-т, 2004.

[12] Ладыженская O.A. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970.

[13] Chorin A. Numerical solution of Navier — Stokes equations // Math. Comput. 1968. Vol. 22, No. 7. P. 745-762.

[14] Amsden A.A., Harlow F.N. The SMAC Method. Los Alamos Sei. Lab. Rep. N LA-4370. Los Alamos, 1970.

[15] Белоцерковский О.M., Гущин В.А., ГЦенников В.В. Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости // Журн. вычисл. математики и матем. физики. 1975. Т. 15, № 1. С. 197-207.

[16] Вабищевич П.Н. Аддитивные схемы для некоторых дифференциально-операторных уравнений // Там же. 2010. Т. 50, № 12. С. 2144-2154.

[17] Червов В.В. Моделирование трехмерной конвекции в мантии Земли с применением неявного метода расщепления по физическим процессам // Вычисл. технологии. 2006. Т. 11, № 4. С. 73-86.

[18] Ковеня В.М. Об одном алгоритме решения уравнений Навье — Стокса вязкой несжимаемой жидкости // Там же. 2006. Т. 11, № 2. С. 39-51.

Поступила в редакцию 25 апреля 2011 г.