Замечания о применении явных схем для численного решения уравнений Навье-Стокса Текст научной статьи по специальности «Физика»

УДК 532.5:519.63

С. И. Мартыненко

ЗАМЕЧАНИЯ О ПРИМЕНЕНИИ ЯВНЫХ СХЕМ ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ УРАВНЕНИЙ НАВЬЕ-СТОКСА

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

E-mail: Martynenko@ciam.ru

Ключевые слова: уравнения Навье-Стокса, численные методы, явные схемы.

В настоящее время явные разностные схемы часто применяются для численного решения уравнений Навье-Стокса на многопроцессорных компьютерах. Безытерационный характер вычислений подразумевает раздельное вычисление компонент скорости и давления. Однако даже при использовании явных схем возможен частичный учет взаимодействия скорости и давления, что позволяет повысить точность и уменьшить объем вычислений.

В [1] для совершенствования вычислительных алгоритмов решения уравнений Навье-Стокса была предложена следующая декомпозиция давления:

p(t, x, y, z) = px(t, x) + py (t, y) + pz(t, z) + pxyz(t, x, y, z), (1)

т.е. давление представляется в виде суммы N +1 слагаемых (N «одномерных» px(t,x), py(t,y) и pz(t,z) и одного «многомерного» pxyz(t,x,y,z)). Далее верхние индексы x, y, z и xyz будут показывать зависимость того или иного слагаемого от соответствующих пространственных координат. Для определения «одномерных» слагаемых px(t,x), py(t,y) и pz(t,z) используются уравнения постоянства массового расхода (интегральные формы уравнения неразрывности), которые рассматриваются как априорная информация физического характера об искомом решении.

Отметим следующие особенности декомпозиции давления (1): 1. Каждое из слагаемых в правой части (1) лишено физического смысла, лишь их сумма имеет физический смысл. В самом деле, в отдельных алгоритмах давление отыскивают с закреплением в одной точке, например в начале координат:

p(t, 0, 0, 0) = p0 = px(t, 0) + py(t, 0) + pz(t, 0) + pxyz(t, 0, 0,0).

Значение р0 может быть перераспределено между слагаемыми в правой части (1) достаточно произвольным образом, откуда следует отсутствие физического смысла каждого отдельного слагаемого. Чаще всего удобнее полагать

рху*(г, о, о, 0) = ро, рх(г, 0) = 0, ру(г, о) = о, р*(г, 0) = 0.

2. Поскольку уравнения постоянства массового расхода являются следствием уравнения неразрывности, то «одномерные» слагаемые рх(г,х), ру (г, у) и р*(г, г) определяются перед вычислением «многомерного» слагаемого рху*(г, х, у, г).

3. Несмотря на представление давления в виде суммы N + 1 слагаемых, уравнения движения всегда будут содержать градиенты только двух из них (соответствующего «одномерного» слагаемого и «многомерного»). Например, градиент давления в уравнении движения по X принимает вид

дх = дх (рх(<,х> + ру с,у) ■+ Р1 +Рху* с,х,у,г)) = % + ^

4. Аппроксимация уравнений постоянства массового расхода должна быть алгебраическим следствием аппроксимации уравнения неразрывности.

5. Применение декомпозиции давления (1) будет наиболее эффективным при моделировании течений с выделенным направлением движения среды. В этом случае градиент одного из "одномерных" слагаемых будет доминирующим.

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

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

V(п+1/2) _ V(п)

Этап I : ---= (п) + Яе-^(п);

VV (п+1/2)

Этап II: Ар =---;

Рис. 1. Каверна с движущейся крышкой и контрольные объемы Vi и V2

V(n+1) _ V(n+1/2)

Этап III : ---= -Vp,

ht

где ht есть полушаг по времени, а V(n+1/2) — промежуточное поле скорости. Нетрудно видеть, что при вычислении промежуточного поля скорости V(n+1/2) не обеспечивается бесконечная скорость распространения малых возмущений, характерная для несжимаемых сред. Кроме того, давление явно не зависит от числа Рейнольдса (Re) и определенную трудность представляет постановка граничных условий для давления.

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

Интегрирование уравнения неразрывности

du dv dx dy

по контрольным объемам V1 и V2 (см. рис. 1) позволяет получить следующие уравнения постоянства массового расхода:

1

J u(t,x,y) dy = 0, (2)

о

1

iv(t,x'y) dx=°- (3)

о

Соответственно разностные аналоги уравнений неразрывности и постоянства массового расхода на разнесенной сетке, показанной на

Рис. 2. Разнесенная сетка

рис. 2, имеют вид

(n+1) „,(n+1) „,(«+!) „,("+!)

>i+1j

j

hx

+

^j+i

j

hy

= 0,

У / ^ Uij ' = 0,

j=1

Nx hx ^

i=1

(n+1)

(n+1)

4j

= 0.

(4)

(5)

(6)

Модифицируем явную схему расщепления по физическим факторам. Для этого запишем уравнение движения по X в виде

ujn+1/2) - j + f diu!)+ d (vu)

h

■t

+

dx 1

dy

(n)

j

dpx

dx

(n+1/2)

+

'd 2u d2 u' Re \ dx2 ду2 /

(n)

j

т.е. на промежуточном временном слое компонента скорости и и одномерное" слагаемое рх будут отыскиваться совместно с привлечением разностного аналога уравнения постоянства массового расхода (5). Для удобства перепишем последнее уравнение в виде

(n+1/2) (n)

uu ht

dpx

dx

(n+1/2)

где

aij =

d (u2) d (vu) 1

dx

dy Re \dx2 dy

d2u d2u +

(7)

(n)

j

Умножая на hy и суммируя по j, получаем

1

ht

N„

N„

N„

j=1

(n+1/2)

j=1

(n)

j=1

dpx

dx

(n+1/2)

N„

а

«j-

j=1

В силу (5) левая часть данного уравнения равна нулю. Первый член в правой части преобразуем следующим образом:

(др! Уп+1/2)_ (др! \(п+1/2) _(др!ч(п+1/2)

^^Ч дх)ь V дх)ь V дх

\ / Ь \ / ь ^ \

поскольку градиент р! не зависит от ], а сумма шагов Ну равна безразмерной высоте каверны, т.е. единице. Тогда выражение для градиента "одномерного" слагаемого р! принимает вид

др! \ (п+1/2) ^

= hyУJ , (8)

откуда

dxL

j=i

ч (п+1/2) / ч (п+1/2) " / ч (п+1/2)

р!)г = (/Л^ + ТЛ^ аы, =0-

ь ь 1 ^—1 1

Вместе с тем, подставляя (8) в (7), получаем выражение для вычисления промежуточного значения компоненты скорости и

и(п+1/2) _ и(п) N

иЫ_^ = _7 Vа-- + а--

Ху

Таким образом, благодаря поправке —7У ^^ а^ промежуточное значение компоненты скорости и удовлетворяет разностному аналогу уравнения постоянства массового расхода (5).

Промежуточные значения компоненты скорости V и "одномерного" слагаемого ру вычисляются аналогичным образом.

Второй этап с учетом декомпозиции (1) записывается в виде

1 (ди ди \(п+1/2) Др!У = - — + дV

7 \дх ду

т.е. сумма "одномерных" слагаемых р! +ру и "многомерное" слагаемое р!у вычисляются на разных временных слоях.

Третий этап модифицированной схемы с учетом декомпозиции (1) принимает вид

V(п+1) _ V(п+1/2) --- = -Ур!у.

Для вычислительного эксперимента положим, что скорость крышки изменяется по закону

п

-50'

UWn) = min (—; i) . w V50 J

Вычисления проводили на равномерной разнесенной сетке 101 х 101 (Н = Нх = Ну = 1/100) при Яе = 100, величина шага по времени выбрана Н = Н/10. Для оценки точности вычислений выбраны следующие параметры:

— максимальная погрешность разностного аналога уравнения неразрывности (4)

RUV = max

ij

(n) (n) (n) (n)

ui+1j — uij_ + vij+1 — vij

hx

h

y

(9)

— максимальная погрешность разностного аналога уравнения постоянства массового расхода (5)

R(n) =

max

Ny j=1

— максимальная погрешность разностного аналога уравнения постоянства массового расхода (6)

R(n) =

max

Nx

(n)

i=1

vij

На рис. 3 показано изменение погрешностей дЦ^ и на первых ста временных слоях (п ^ 100). Погрешность « , поэтому на рис. 3 она не приводится.

Рис.3. Изменение погрешностей RUV и ВЦ^ при моделировании течения в каверне

j

У1 'u(t,x, 1)= Uw(t) v(t,x, 1)=0

1 о д о

о" 'S г-Г ■(о4

i \ S>

о II т ) о II

а сГ Ä r-f

0 u{t,x,0)=0 v(t,x,0)=0 1 X

Рис. 4. Геометрия задачи о течении в каверне (а) и между параллельными пластинами (б)

Отметим еще одно полезное свойство декомпозиции давления (1). В методе расщепления по физическим факторам давление отыскивается из решения второй краевой задачи для уравнения Пуассона, которое в двумерном случае принимает вид

д 2р д2 р 1 / ди дю\ дх2 ду2 Ы \ дх ду) Рассмотрим возможные случаи постановки граничных условий для давления на примере моделирования нестационарного течения в каверне и между параллельными пластинами (рис. 4). Используя первую формулу Грина, получаем следующее соотношение:

i i

Г dp ^ Г dp

I dx x=1 J дх

) о

x=0

dy + /f

о

У=1

л [dp dx — — J dy

о

dx = (10)

y=0

1 1

= ^ i u(t(n+1/2), 1, y)dy — 1 i u(t(n+1/2), 0, y)dy, ht J ht J

которое является условием разрешимости данной краевой задачи. При моделировании течения в каверне правая часть уравнения (10) обращается в нуль в силу граничных условий прилипания (и(Ь, 0, у) = 0 и и(Ь, 1,у) = 0). Поэтому для задачи о каверне естественной является постановка для давления однородных граничных условий второго рода

др дх

x=1

= 0, ^ дх

x=0

= 0, дг

dy

У=1

=о, дт

dy

= 0,

y=0

1

которая обращает уравнение (10) в тождество, т.е. приводит к разрешимости второй краевой задачи для уравнения Пуассона.

Несколько сложнее выглядит постановка граничных условий при моделировании течения между пластинами. В традиционных вычислительных алгоритмах компоненты промежуточного поля скорости не удовлетворяют уравнениям постоянства массового расхода, поэтому в данном случае правая часть уравнения (10) отлична от нуля. Зададим по аналогии с задачей о каверне нулевые градиенты давления на омываемой поверхности пластин. Тогда уравнение (10) принимает вид

dp дх

Ж=1

л f дР

dy - дХХ

dy

ж=0

1 Г , , 1 ЛЛ , , 1

0

1 1

= — J «(¿(п+1/2), 1, у) ¿у - - у «(¿(п+1/2), 0, у) ¿у.

0 0 Очевидна трудность задания градиентов давления во входном (х = 0) и выходном (х = 1) сечениях для обеспечения условия разрешимости второй краевой задачи для уравнения Пуассона.

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

Метод искусственной сжимаемости. В [6] предложен метод расчета стационарных течений несжимаемых сред, основанный на использовании модифицированного уравнения неразрывности

I + ^ = 0,

дЬ

т.е. вычисления проводятся по схеме

V(п+1/2) _ V(п)

ht

= -(Vp)(n) - (V(n)V)V(n) + Re-1AV(n),

р(п+1/2) _ р(п)

Р-^ = - (п)

Ъ

Численное решение уравнений Навье-Стокса проводится посредством счета на установление, физический смысл имеет только ста-

1

Рис. 5. Параллельные пластины и контрольные объемы V и У2

ционарное решение. Рассмотрим возможность ускорения сходимости метода искусственной сжимаемости на примере моделирования развивающегося симметричного течения между параллельными пластинами (рис. 5).

Интегрирование уравнения неразрывности по контрольным объемам VI и У2 (см. рис. 5) позволяет получить следующие уравнения постоянства массового расхода:

1 1

/ и(Ь,х,у) ¿у = и(Ь, 0, у) ¿у;

J v(t,x,y) dx = - J (u(t,L,0 - u(t, 0,£)) dy, 0 0 которые на равномерной разнесенной сетке аппроксимируются выражениями

Ny

hvYl

(n+1)

'У / у "13

3=1

Ny 3 = 1

Nx

h.£ v!;,+1) = -h

i=1

3 (u^1) _ u(" + 1)

У / v \^UNx + 1k U1k k=1

(11)

(12)

Поскольку уравнение (12) связывает значения двух компонент скорости, то имеет значение порядок отыскания компонент скорости из уравнений движения. В данном случае разностный аналог уравнения движения по X принимает вид

(n + 1) _ (n)

uij uij

ht

дрХ

дх

(n+1)

+ aij,

где

dpxy d(u2) д (vu) 1 f d2u d2u

(n)

а-- = —1----—1---—1 +----+ „

1 \ дх дх ду Яе \ дх2 ду2

Умножая на ку и суммируя по ], получаем

1 / \ / дрх \ (п+1) ^

к ^^и(7п+1)-) = ку[) + а1.

* V 1=1 1=1 / 1=1 ^ ' г 1=1

Отсюда с учетом (11) левая часть выражается через известные значения на входной границе х = 0. Тогда

Я„х\ (n+l) 1 / ^у ^у

др \ _ 1 h Vu(n+1) h Vu(n) д^ ) _ -h \nv/ -uij - hy / >uij I 1 hy

hy S u1n - hy u1j) +

j=i 3=1/ j=i

Соответственно компонента скорости вычисляется так:

u(n+1) _ u(n) 1

hy / J u1j - hy Z^ u13 - hy ai3 + a3 •

к* к* ,

1=1 1=1 / 1=1 Промежуточные значения компоненты скорости V и "одномерного" слагаемого ру вычисляются аналогично.

На данном временном слое вычисления завершаются расчетом "многомерного" слагаемого рху:

(рху)((1П+1) - (Рху)(п) _ ( ди + д^N (п) к* Vдх 1

Вычисления проводили на равномерной разнесенной сетке 201 х х 201 (к = кх = ку = 1/200) при Яе = 100 и Ь = 10, величина шага по времени выбрана к* = к/10. Для наглядности величина скорости на входе задавалась в виде

(п) . ( п Л („)

На рис. 6 показано изменение скорости на оси симметрии. Метод искусственной сжимаемости строился по аналогии с явными схемами, предназначенными для моделирования сжимаемых течений, поэтому скорость распространения возмущений с входной границы ограничена из-за наличия производной р* в уравнении неразрывности. Использование в модифицированном методе искусственной сжимаемости уравнений постоянства массового расхода (11) и (12) обеспечивает бесконечно большую скорость распространения малых возмущений в несжимаемой среде, несмотря на наличие производной р* в уравнении неразрывности. Указанное обстоятельство позволяет быстрее получить стационарное решение, если оно существует.

Рис. 6. Изменение скорости на оси симметрии

На рис. 7 показаны изменения погрешности дЦ^ (9) и максимальной погрешности разностного аналога уравнения постоянства массового расхода (11), определяемой как

R(n) =

max

i

Ny Ny

hy u(3 - hy ^^ u1j

3=1

3=1

В модифицированном методе величина Лц ) остается на уровне погрешностей округления.

Cжимаемые среды. Численное решение уравнений Навье-Стокса

д U д E д F _ dt дх дy

U=

/р\

pu pv

w

E=

/

pu

pu2 + p - Txx puv - Txy

+ p)u

uTxx vTxy + (Ix)

pv

puv - Txy 2

pv2 + P - Tyy

\(e + p)v - uTxyy - vTyy + QyJ

Рис. 7. Изменение погрешностей Ли? и ) при моделировании течения между

пластинами

часто проводят в переменных скорость-плотность. Одна из первых явных схем (схема Мак-Кормака [7]) состоит из двух этапов: Предиктор

TT _ TT(n) _ h (E(n) - E(пЛ - — ( F(n) - F U ij _ T ij h \Ei+1j Eij ) h \F ij+1 F ij hX hV

,(n)

(n)

Корректор

T

(n+l) _

T

(n) ij

ht

+ U — —(E(n) - E

+ ij h ij Ei-ij

(n)

_ (Fw _ F(n)

hy V ij ij-1

Данная явная схема имеет второй порядок аппроксимации как по пространству, так и по времени [8].

Для наглядности моделирование течений сжимаемых сред рассмотрено на примере течения в каверне с движущейся крышкой. Интегрирование уравнения неразрывности

др д (ри) д(ру)

м+

+

_ 0

dx dy

по контрольным объемам Vi и V (см. рис. 1) позволяет получить

уравнения постоянства массового расхода в виде

х 1 1

др (*

— АуАх + р(г,х,у)п(г,х,у) Ау = 0; (13)

0 0 0 У 1 1

д р С

— АхАу + р(Ь,х,у)ю(Ь,х,у) Ах = 0. 0 0 0

Соответственно разностный аналог уравнения постоянства массового расхода (13) имеет вид

N

hxhy ЕЕ (|) , + hy E(Pu)lr1) = 0 . (14)

k=i j=i ^ ' kj j=i

По аналогии с (7) запишем уравнение движения в виде

МГ" - (puj (dp' V"+1)

,j h^" = -l + ßj- (15)

Умножая на hy и суммируя по j, получаем

1 / Му \ / дгХ\ (п+1)

1 Ну Е(рп)!Г1} - Ьу Е(рп);;М = — дх) + Ьу Е ,

* V 3=1 3=1 ) ^ / » 3 = 1

откуда с учетом (14) следует выражение для градиента "одномерного" слагаемого рх

др^ \ ("+1) дх , г

НхНу^Е Н + Ну £ (рп)(? + Н^Е в,.

к=1 3=1 * 3=1 / 3=1

Аналогично вычисляются градиенты остальных одномерных слагаемых в (1). Далее полученные градиенты рассматриваются как источ-никовые члены.

Следует отметить, что градиенты «одномерных» слагаемых вычисляются с некоторой погрешностью, которая вызвана использованием приближения р,г+1) в уравнении постоянства массового расхода. Данная погрешность компенсируется при вычислении "многомерного" слагаемого рху следующим образом:

(р, е) ^ уравнение состояния ^ р ^ рху = р — рх — ру . Таким образом, при решении уравнений Навье-Стокса для сжимаемых

сред часть вычислений проводится в переменных скорость-давление, а часть — в переменных скорость-плотность. Полученные соотношения остаются верными и для уравнений Эйлера.

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

Работа выполнена при поддержке РФФИ (проект 09-01-00151).

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

1. Мартыненко С. И. Совершенствование вычислительных алгоритмов для решения уравнений Навье-Стокса на структурированных сетках // Вестник МГТУ. Сер. Естественные науки. - 2008. - № 2. - C. 78-94.

2. Marty nenko S. I. A physical approach to development of numerical methods for solving Navier-Stokes equations in primitive variables formulation // Int. J. of Comp. Science and Math. - 2009. - V. 2, № 4. - P. 291-307.

3. Мартыненко С. И. Замечания о вычислении давления при численном решении уравнений Навье-Стокса // Математическое моделирование. - 2010. -Т. 22, №3.-С. 105-119.

4. Белоцерковский О. М., Гущин В. А., Щенников В. В. Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости // ЖВМ и МФ. - 1975. - T. 15, №1. - C. 197-207.

5. Гущин В. А. Пространственное обтекание трехмерных тел потоком вязкой жидкости // ЖВМ и МФ. - 1976. - T. 16, № 2. - C. 529-534.

6.Chorin A. J. A numerical method for solving incompressible viscous flow problems // J. Comp. Phys. - 1967. - V. 2. - P. 12-26.

7. MacCormackR. W. The effect of viscosity in hypervelocity impact cratering // AIAA Paper 69-354, Cincinnati, Ohie. - 1969.

8. Андерсон Д., Таннехил Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен: В 2 т. / Пер. с англ. - М.: Мир, 1990.

Статья поступила в редакцию 4.12.2010

Сергей Иванович Мартыненко окончил в 1988 г. МВТУ им. Н.Э. Баумана. Канд. физ.-мат. наук, научный сотрудник отдела химмотологии и спецдвигателей ФГУП "Центральный институт авиационного моторостроения им. П.И. Баранова".

S.I. Martynenko graduated from the Bauman Moscow Higher Technical School in 1988. Ph. D. (Phys.-Math.), researcher of department for chemistry of motor combustive-lubricating materials and special motors of the Federal State Unitary Enterprise "Central Institute of Aviation Motors n.a. P.I. Baranov".