УДК 532.501.1:533.6.011 Л. И. Скурин
Вестник СПбГУ. Сер. 1, 2003, вып. 4 (№25)
ПАРАЛЛЕЛЬНАЯ СХЕМА ИТЕРАЦИОННО-МАРШЕВОГО МЕТОДА ИНТЕГРИРОВАНИЯ УРАВНЕНИЙ НАВЬЕ—СТОКСА*
В работах [1,2] и в других работах автора развивается итерационно-маршевый метод (ИММ) численного интегрирования уравнений Навье—Стокса для газа и жидкости. Метод применим для решения многомерных стационарных и нестационарных, внутренних и внешних задач. Отличительная особенность метода в том, что решение всех типов задач основано на единой процедуре — маршевой по пространственной координате для двумерных задач и по двум координатам для трехмерных задач. Другими словами, для всех задач используется единственная алгебраическая процедура — решение системы алгебраических уравнений для вектора искомых сеточных функций на лучах, поперечных по отношению к маршевым осям. Аналитически показано, что метод имеет безусловные устойчивость и сходимость для задач гидродинамики. Применимость метода к расчету течений сжимаемой жидкости продемонстрирована на многочисленных примерах [2,3].
В настоящей работе предлагаются в рамках ИММ вычислительные схемы для стационарных и нестационарных задач, пригодные для численного решения с помощью параллельного алгоритма. Показано, что эти схемы также обладают безусловной сходимостью. При их использовании вычисления, связанные с упомянутой алгебраической процедурой, могут осуществляться для всех лучей одновременно на разных процессорах в пределах одного временного шага или глобальной итерации. Таким образом, практически весь объем вычислений (в пределах каждого временного шага либо глобальной итерации) представляется в виде суммы большого числа независимых между собой и близких по объему частей.
Рассмотрим вначале схему для стационарных задач. В рамках ИММ для них оказывается эффективной [1, 2] модифицированная схема, представляющая собой схему основанную на принципе установления по времени. В соответствии с ней каждый временной шаг реализуется применением одной маршевой процедуры. Чтобы построить схему, реализуемую параллельным алгоритмом, изменим несколько аппроксимацию производных от составляющих скорости по маршевой координате. Рассмотрим следующий конечно-разностный аналог двумерной нестационарной системы уравнений Навье—Стокса для несжимаемой жидкости, записанной в декартовой системе координат и с "замороженными коэффициентами":
~ изт-1 т ~ _ д
Дж 2Ду ~ '
— и^т—1 . | \ и<^т — т—1 — г+-1 и^+1т — и^—1т ,
-М-+ Ы-Ах-+ -2Ау-+
па +1 —па иа —2иа + иа и а+1 —2исг + и а—1
1 зт _ аз + 1го Т щ а]т-1 Т "дго-!
Ах ~ ЯеАу2 ЯеАх2
* Работа выполнена при финансовой поддержке Минобразования РФ (грант Е02-4.0-133). © Л. И. Скурин, 2003
,, ,Vjm ~ SlVjm-l ~ S2Vjm-l ^ Vj+lm-Vj-lm + M-T-- +«0-гГл~.-+
Ai 10 Ax 2Ay
na — na va — 2va + va Va+1 — 2Va + va-1
Pj+lm Pj-lm _ uj+lm ZUjm > uj-lm Ujm-1 jm T ujm-1 n
2Ay ~ ReAy2 ReAx2 ' U
Si = (1 + u0/\u0\)/2, S2 = (1 - «o/|«o|)/2,
где индексы a, j — номера узловых точек соответственно вдоль маршевой оси x и ортогональной ей y, m — номер временного слоя, u,v —проекции вектора скорости на оси x,y соответственно, p —давление, Re — число Рейнольдса.
Схема (1) отличается от модифицированной схемы, рассмотренной в [1], тем, что значения сеточных функциий в точке с индексом a — 1 берутся со слоя m — 1 (а не с m, как в [1]).
Модифицированная схема на каждом временном слое реализуется с использованием маршевой процедуры, то есть путем последовательного отыскания решения на лучах a = 2, 3,...,N — 1. Отыскание решений на этих лучах для системы (1), как нетрудно видеть, можно осуществлять в произвольном порядке, а значит и с использованием параллельного алгоритма. Однако важно выяснить имеет ли место для этой системы устойчивость. Рассмотрим этот вопрос, используя метод разделения переменных. Подставим в эту систему вектор неизвестных в виде
dajm = (u*,v*,p*)TЛт exp[%a + aj)].
Матрица системы уравнений для u* ,v*, p* будет иметь вид
Л — e-%1 Xiry sin a 0 \
XR1 — R2 0 eÍY — Л I , (2)
0 ЛR1 — R2 ЛГ sin a J
где
R1 = rt + \uo\ + 2/(ReAx) + 2ry(1 — cos a)/(ReAy) + voryi sin a, R2 = rt + |uo|(S2eÍY + Sie-ÍY) + (eÍY + e-ÍY )/(ReAx), rt = Ax/At, ry = Ay/At. Этой матрице соответствует характеристическое уравнение
(Л — e-ÍY )(Л — eÍY )^R1 — R2) + (ЛR1 — R2)X2A2 =0, A2 = r2y sin2a.
Оно имеет три корня:
Д2 COS7 / cos2 7 1
Al~ñ? '3 ~ ТТа2 у (i + А2)2 ~ ТТА2' ( )
Рассмотрим модуль первого корня, обозначив 2ry(1 — cosa)/(ReAy) =
а
а
2 _ Irt + |ц01(5*2 + Sí) COS7 + 2/(ReAx) COS7 + г|цо|(52 - Sí) sin7|2
\rt +\и0\ +2/(ReAx)+£ + v0ryi sin а\2
(rt + 2/(ReAx) COS7 + |uo| cos7)2 + \u0\2 sin2 7 = (■n + 2/(ReAx) + £ + bol)2 + v^A2 =
_ [rt + 2/(ReAx)cos7]2 + 2[rt + 2/(ReAx) cos7)]|w0| C0S7 + |w0|2 = [rt + 2/(ReAx) + £]2 + 2[rt + 2/(ReAx) + £]K| + Ы2 + v2A2 '
Поскольку £ > 0, rt > 0, то
[rt + 2/(ReAx) +£] — [rt + 2/(ReAx) cos 7] > 0
и, следовательно, всегда |Ai| < 1.
Рассмотрим оставшиеся два корня.
| cos 7 ± г Vх 1 + А2 - cos2 7I 1
|Л2'э| --1 + А2-- VTTW - L
Таким образом, все корни не превышают единицы и, значит, схема (1) безусловно устойчива.
Очевидно, что использование в (1) сдвинутых по времени аппроксимаций пространственных производных приводит к потере точности аппроксимации. Поэтому следует ожидать, что при расчетах на одном процессоре эта схема потребует большего числа шагов до достижения сходимости, чем модифицированная. В соответствии с первыми результатами ее апробации она требует затрат машинного времени на 5-20% больших по сравнению с модифицированной схемой (рассчитывалась обтекание обращенной назад ступеньки, использовалась консервативная форма записи уравнений).
Решения нестационарных задач базируются на основной схеме ИММ [1, 2]. В этом случае решение на каждом временном слое находится с помощью глобальных итераций (ГИ), каждая из которых представляет собой маршевую процедуру. Введем в соответствующую схему изменения, аналогичные изменениям для модифицированной схемы. Получим:
Ujm ~ (Ujm )S Vj+lm ~ Vj-lm _ q
Ax 2Ay ~ '
uljm — u<jm-1 , 1 ,UJm — S1(ujm1)S 1 — S2(uim+1)S 1 uj+1m — uj-1m ,
-At-+ Kl-Ax-+ -2Ay-+
(Pjm )S 1 ~ Pjm _ uj+1 то ~ + U j-1 то (Ujm ~ + (Ujm
Ax ~ ReAy2 ReAx2
Vlm-Vjm-I , , tfm-SlivJ-1)"-1- Stiv^1)'-1 , V?+lm - V«_lm _
-At-+ Kl-Ai-+ -2Ay-+
Pj+lm ~Pj-lm _ Vj+lm ~ ^Vjm + vj-lm (Vjm 1 ~ ^Vjm + (vjm У$ ,дч
2 Ay ~ ReAy2 ReAx2 ' 1 '
где s —номер ГИ и предполагается, что все искомые функции относятся к текущей s-й ГИ за исключением тех, что имеют индекс s — 1.
Схема (4) отличается от основной схемы, рассмотренной в [1], тем, что значения сеточных функций в точках с индексом а — 1 берутся с (s — 1)-й ГИ (а не с s-й, как в [1]). Благодаря этому, отыскание решений на поперечных лучах по системе (4), как нетрудно видеть, можно осуществлять с использованием параллельного алгоритма.
Исследуем устойчивость ГИ. В соответствии с методом разделения переменных, представим вектор неизвестных на данном временном слое в виде
(j)s = (u*,v*,p*)T Xs exp[*(7a + aj)].
Этот вектор должен удовлетворять системе, которая следует из (4) с учетом того обстоятельства, что величины с индексом m — 1 фиксированы.
Матрица системы для u*,v*,p* будет иметь вид (2), где меняется лишь выражение для К2. А именно, в этом выражении исчезает положительное слагаемое rt. Таким образом, в рассматриваемом случае мы имеем те же три корня (3) и каждый из них не превышает единицу. То есть, схема (4) безусловно устойчива.
Некоторым усложнением, связанным с использованием схемы (4), является необходимость удерживания в памяти машины результатов расчета (s — 1)-й итерации при расчете s-й итерации.
Аналогично (1), (4) в рамках ИММ могут быть построены параллельные схемы для случая аппроксимации маршевых производных со вторым порядком [2], для задач о течениях сжимаемой жидкости [3], для трехмерных задач [4].
Summary
Skurin L. I. The parallel scheme for Navier—Stokes equations on the base of an iterative space-marching method.
Computational schemes that suitable for solving steady and unsteady problems with parallel algorithm a presented. It is proved that they have unconditional stability. The proofs are given for incompressible fluid schemes that are special cases of compressible ones for a zero Mach number value. When these schemes are used the computations along each line that is transverse to marching axes can be performed simultaneously on different processors on any time layer or global iteration.
Литература
1. Скурин Л. И. Итерационно-маршевый (по пространству) метод решения задач механики жидкости и газа // Математическое моделирование. 2000. Т. 12. № 6. С. 88-94.
2. Skurin L.I. Iterative Space-Marching Method for Incombressible and Compressible Full Navier-Stokes Equations / N. Satofuka (Ed.) // Proceedings of the First International Conference on Computational Fluid Dynamics, ICCFD, Kyoto, Japan. 10-14 July 2000. P. 319-324.
3. Полянский А. Ф., Скурин Л. И. Численное моделирование вязких до- и сверхзвуковых течений газа с использованием итерационно-маршевого метода // Вестн. С.-Петербург. ун-та. Сер. 1. 1998. Вып. 3 (№15). С. 111-115.
4. Полянский А. Ф., Скурин Л. И. Решение пространственных задач гидродинамики итерационно-маршевым методом // Вестн. С.-Петербург. ун-та. Сер. 1. 1999. Вып. 2 (№8). С. 90-96.
Статья поступила в редакцию 3 марта 2003 г.