Вычислительный метод оценки поля скорости ветра в пограничном слое атмосферы на основе уравнения Навье - Стокса Текст научной статьи по специальности «Математика»

_МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ_

УДК 519.6

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

© 2005 г. Н.И. Каргин, В.И. Наац

The Navie-Stoks vector equation, with the help of that may be state the value of wind velocity field consideration terms in it, which are formed of turbulent viscosity in the frontier layer of the atmosphere, is consider. The construction of calculation scheme for this equation, which is based on coordinate wise splitting method, is described.

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

В соответствии с аэродинамической моделью, компоненты Vx(P,t), Vy(P,t), Vz(P,t) поля скорости ветра

V(P, t)= VX(P, t)i+ Vy(P, t) j+ V(P, t) k удовлетворяют системе уравнений Навье - Стокса следующего вида [1, 2]:

+ V, (P,,). ^ + V, (P,,). + V (P,,). ^ =

= FPt) __L_ .№0 ,V2K ( P, ,), (1)

р( P, t) р( P, t) dx р( P, t) x

dV, ( P, t) dV, ( P, t) dV, ( P, t) dV, ( P, t)

+V, (P, t). -+V, (P, t). -Z--+Vz (P, t). -Z--=

dt dx dy dz

=Fy--L-.дЕШ(P,t),

р( P, t) р( P, t) dy P( P, t) y

д^Ш+к(р,,)+у (р,,)+у(р,,)=

дt дх ду дг

=(р, t),

р( Р, t) р( Р, t) дг р( Р, t)

где Р(х,у,1) еП, (П- некоторая область, ограниченная П ); р(Р, ^ - атмосферное давление в точке Р, кг-м-1-с-2; ¥х(Р, ¥у(Р, О, ¥/Р, ^ - компоненты силового поля Ж(Р^), кг-с-2-м-2, действующего на единичный объем в пределах области П; р(Р, кг-м-3 - плотность воздуха; /, кг-м-1-с-1 - динамический коэффициент вязкости. Помимо вязкости / может быть ис-

/ 2 -1

пользована и величина V =-, м -с , именуемая кинематическим ко-

Р( Р, t)

эффициентом вязкости. Система уравнений (1) может быть записана в векторной форме:

дУ+(у .у)у = --Ур + ЛУ V (2)

дt Р Р Р

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

Если считать неизвестными функции Ух, Уу, у в системе (1), то следует учесть наличие между ними дополнительной функциональной зависимости. Речь идет об условии непрерывности вида:

^ + ^ + = о. (3)

дх ду дг р Л

Обычно считается, что плотность воздуха р(Р, ^ остается неизменной во времени в пределах некоторого мыслимого эксперимента, и тогда (3) превращается в более простое условие луУ=0. При подобном подходе остается вопрос об источнике нестационарности в уравнениях (1), особенно если учесть, что компоненты ¥х, ¥у, ¥2 силового поля Ж в пределах

¥

пограничного слоя в естественных условиях малы. В лучшем случае —

Р

принимается равным g - ускорению земного притяжения, остальные компоненты ¥х и ¥у практически равны нулю. Давление р и плотность воздуха р связаны друг с другом уравнением состояния газа

р=рЯ-Т,

в котором Я - универсальная газовая постоянная; Т - температура. Если р считать постоянной, то Ур=р-Я-УТ, и тогда УТ является определяющим векторным полем для уравнений (1) (то же (2)).

Поскольку структура для всех трех уравнений в системе (1) аналогична, то более подробно рассмотреть можно одно из них, например, первое, переписав его предварительно в виде:

дУх „ дУх дУх дУх м о ^ 1 др —- + Ух--- + Уу--- + ---- — ■У1¥г = —----—. (4)

д/ дх ду дг р Р Р дх

Это уравнение входит в систему подобных еще двух уравнений, выписанных аналогичным образом на основе второго и третьего, и является достаточно сложным с точки зрения построения соответствующей вычислительной схемы. Естественным с первого взгляда кажется метод последовательных приближений, который позволяет, например, решать уравнение (4) относительно компоненты Ух независимо от двух остальных Уу, Уг, которое в свою очередь может служить далее источником получения информации об этих компонентах. Если обозначить через V номер соответствующей итерации, то уравнение (4) можно писать в виде дТ/ (м1) дУ (м1) дТ/ (м1)

ЕТх-+ /М -+ У<У) .ЕКх— +

д/ дх ду

( п (5)

+¥М) дУх -Е.у1уМ+1) = -1 др, г дг р х р р дх

0,1,2... Итерационная схема (8) для определения Ух выписана таким образом, чтобы решаемое уравнение на V -м шаге было линейным относи-

(у+1)

тельно искомой функции Тх , что требует предварительного знания

(V) (V) (V)

Ух , Уу и Уг . Аналогичным образом могут быть выписаны и две остальные вычислительные схемы, позволяющие определить соответствен-

^+1) ^+1) (V) (V) (■+1) (■+1) (■+1) (V)

но Уу по Ух , Уу , Уг и Уг по Ух , Уу , Уг . Сказать что-либо определенное о сходимости и устойчивости подобного итерационного алгоритма затруднительно.

В пределах настоящего исследования обратимся к методу расщепления [4, 5], где исходным является предположение, что переменная / меняется в пределах некоторого элементарного интервала [/,/4], 6=0,1,2., на

котором величина А/=/6+1-/6 может быть сделана сколь угодно малой. В

&+1) 6+1) 6) связи с этим Ух может быть заменено на Ух = Ух(Р,/6+1), Уг =Ух(Р,/).

Это же замечание касается и двух других искомых функций Уу и Уг. В

этом отношении рассматриваемые здесь методы в определенной степени

близки друг другу. Что же касается линеаризации первого уравнения, то в

методе расщепления она носит локальный характер по переменной /, т.е.

относится лишь к интервалу [/,/4], и потому ее следует признать более

приемлемой и обоснованной. Прежде чем построить соответствующую

вычислительную схему, примем следующее предположение. Выражение

п о

— V Ух в дальнейшем будем писать в виде

р

д

Кх

]+-*-[ ^к„

дх у дх ) ду у ду ) дг у дг где функции Кхх(Р,(), Кху(Р,1) и Кхг(Р,() соответствуют так называемым коэффициентам турбулентного обмена в направлениях осей выбранной координатной системы. Если эти функции не зависят от временной и пространственных переменных, то из этого следует, что необходимо приме-

п 2

нять выражение — ■У Ух. Следует подчеркнуть, что в подобном случае

Р

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

дУх дУу ЗУ, считая, что —— +—— + -

= 0.

дх ду dz Из этого выражения следуют соотношения вида

ду

дх

iy + ду

ду дz

\

5Уу_

ду

зу +дК

дх 3z

ду

dz

У +V ^

dx dy

Введем их в уравнения исходной системы (1). Тогда

У dt

дУу dt

- Vx

-y У у

3yL+y

ду dz

dVy, дх

т. дУх „ дУх + Vy—- + Vz—-

--V

y dy дУх

дх dz

Fx-dp+м'У2Ух

дх .

+ V,

dz

A = PfFy+

dz р\/ ду

(6)

dVz „ dVz „ dVz „ —-+Vx—-+Vv—- - Vz х у ду z

dt

дх

дVx dV,

л

+

дх ду

= рр {Fz-f + ^'V2Vz

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

dV,

(V+1)

dt

--К

(v+1)

dV/v) dV.

(v)

ду

dz

dV V+1) + V V) ^-+

у

ду

+V(v) ^х^1 -^.V2F(V+1) = 1 fF -дР z ^ х I х ^

dz р р\ дх

v = 0, 1, 2, ...,

&+1)

в которой прогноз значения Ух осуществляется уже по паре функций

(V) (V)

Уу и У2 . Таким образом, осуществлена более корректная линеаризация первого уравнения системы относительно искомой функции Ух. Это достигнуто за счет введения в исходную систему дополнительного условия (3).

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

& (Л о - (^ -f

р l дх

(-)=р (F-I

dVy dVz

ayz =— +—-

dy dz

dVx dVz axz = —- + —-дх dz

a (P.t) = P ( Fz-l

dVx dVy a = —— + -

"xy '

в которых система (6) примет вид

V „ „ dVx „ dVx —- - ayzVx + Vy—- + V—-dt yz x y dy z dz

dx

K-

dx

V

dx

dy

--Ч Ky V 1-i.fKxz ^ | = ßx(P,t),

dy | ^ dy ) dz | xz 5z 1 x

dVy, dt

- + Vx

dVy_ dx

-- a„V„ + V,

dVy_ _d_ dz dx

(

Ky

Л

dx

d( dVv Л д

dy

K

yy

dy

dz

K

yz

dVy Л dz

(7)

= Qy (P, t),

V „ dvz w d¥z d( d¥z —L + Vx—L + Vv—L - a_Vz--1 Kzx—z

Л, Л у XV Z -Л ZX

dt dx öy dx l dx

AI K V

dy l * dy

■4 k, dVz

dz l zz dz

= Qz ( p, t).

Напомним, что при решении первого уравнения (7) значения Уу и Уг считаются известными. Для метода расщепления это Уу и У2 соответственно. В силу идентичности структуры уравнений, входящих в систему (7), достаточно ограничиться записью вычислительной схемы для одного из них. В частности в первом уравнении системы (7), обозначая через ф искомую функцию Ух(Р,/), будем иметь

ф - «угф + Уу Ф Уг Кх Щ-Ц Ку д-Ф\-

у у ду г дг дх { хх дт \ ду \ х ду) (8)

д ( „ дф

-iz 1 ^ ^ l = Qx (P't).

Для трехмерной задачи (8) метод расщепления строится по следующей схеме: задача 1:

-( Кхх р дх У дх

Ф -аУФ --I Кх^1 | = • Qx(P,t,M), (9)

задача 2:

(Pt ) /х (P,0h если t = 0 (j = 0), P n

Ф(P, tj) = \ (Pt ) -IT P eQ

J I ф3 (P, tj), если j = 1,2,...,

¿(P,t}) = Vx(P,t}), если P eQ, t- *t j

задача 3:

¿>2 +Уy д^-д^Кхуд^ l = ®2 •Qx(P,t), (10)

i^(P, t-+1), если P eQ,

%(P,t.) H- -

j |Vx (P, t}), если P eQ,

V (P, t1), если P eQ, j = 1,2,., Уy (P, t) = y _

у V(P,t}), если P e QP,

tj j>

дФз д( дф

Фз + K-3^~dzI Kxz-3^ | = ®3 • Qx(P,t), (11)

i^2(P, tj+1), если P eQ,

Фз(Лtj) 4- -

j \VX (P, t}), если P eQ,

Vz (P, t}), если P eQ, j = 1,2,

У (Р, ') = ,

г У (Р, ), если Р еД

^ а>1 +ю2 +ю3 = 1.

В качестве решения уравнения (8) принимается

р(Р,()=р3(Р,() при < <(]+1, Ух(Р,(]+1)=р (Р,(]+) Аналогично вычисляются значения функций Уу(Р,^+1) и Уг(Р,^-+1). Обозначая, как и раньше, рР^)=Уу(Р,(), выпишем исходное уравнение, аналогичное (8), в виде

• ТЛ др т. др д ( др\

р+Ух—-а„р+уг——I кух— I-

т х ^ хгт г ^ -л ух

дх дг дх у дх)

-йкуРтйКгЩ^ <12>

Схема расщепления для (12): задача 4:

vf -I (Qy ()•

Я\( P, t,) =

\Vy (P,0), если t = 0 (j = 0),

задача 5:

j [<p3(P, t), если j = 1,2,..., ^(P,t3) = Vy (P,t3), если P eQ,,

\VX (P, t,), если P eQ, j = 1,2,

V (P, t) = \- _

x \Vx (P, t}), если P eQ,

tj *t *tj+i;

P eQ,

Ч>2 _ axzV2'

dy

Kyy = 1 - Qy (P, t X

<Рг(P,t,) = \

U( P, tj+1), если P eQ, (P, tj), если P eQ,

задача 6:

tj *t * tj+i;

Фз + Vz ^ _±f * M 1 = 13 - Qy (P, t)

n z dz dz ( yz dz 1 3 y

(13)

(14)

(15)

Гр2(Р,+1), если Р еП, \Уу (Р, tJ), если Р еП, (Р,t]), если Р еП, ; =1,2,... (Р, tJ), если Р еП,

^ <t <^+1, а>1 +ю2 +юъ ■ 1.

Решение уравнения (12) р(Р^)=р3(Р^) при ^ <t Уу(Р,^+1)=р (Р,^+1).

Для определения третьей компоненты У2 в момент времени t=tj+1 исходим из уравнения

р + V,( Уу(* р-А(^^-АГк (

' Х ^ у Л Х^ ^ -Х у, —у I

дх ду ох ^ дх) ду { ду,

-I Г г )■ с- - р, '',

где, как и выше, через (обозначено У2^,М). Имеем поэтапно задача 7:

задача <

Ф + Vx д-Фф ~дтf Kzxдp-1 = ®1 • Qz (P, t), (17)

дх дх ^ дх )

V (P,0), если t = 0 (j = 0) Ф (P, t,) = i z , P e Q,

m pP,tj), если j = 1,2,...

ф( P, t3) = V (P, t3), если P eQ,

|Vx (P, tj), если P eQ, j = 1,2,.

Vx (P, t) H- -

x \Vx (P, t}), если P eQ,

tj j>

pдф1

ду ду l ^ ду

P2 + Vy-дг-^Куу-^ l = -2 • Qz (P,t), (18)

\p (P, t..+1), если P eQ,

<P2(P,t.) - ,

V (P, tj), если P eQ,

¥у (P, t}), если P eQ, j = 1,2,

задача 9:

V (P, t) = ,

у \Vy (P, t}-), если P eQ

tj *tj+i;

Pp3 - ахуФз -dz l Kzz dp |=®3 • Qz (P, t), (19)

Гр2(Р,tj+l), если Р еД,

Рз(Р,ti) Н- -

У (Р, t]), если Р еД,

^ <t а1 + ю2 +ю3 = 1.

Решение уравнения (16) определяется соотношениями

р(Р, t) = Рз(Р, t) при ^ < <^+1, Уг (Р, ^+1) = р(Р, ^+1).

На этом завершается процесс прогноза У(Р,^+1) по значениям У(Р,^ в пределах описанного алгоритма.

Сделаем несколько замечаний к обоснованию изложенного выше алгоритма. Правомерность применения вычислительных схем (9)-(11), (13)-

(15) и (17)-(19) для определения Ух(Р,^+1) на основе уравнений (8), (12) и

(16) обосновывается следующим образом. Из (8) путем интегрирования в пределах интервала имеем

Pi(t)-Pi(tj )= J

ау2Ф1 +дх I Kxx ~Pl\ + a1 • Qx

dt. (20)

Поскольку р1(^=р(^ для всех внутренних точек Р из области П, то выражение (20), если к интегралу справа предварительно применить теорему о среднем, можно переписать в виде

a( Kxx ^Т" |+1 - Qx

\t=Cj '(tj+1 _tj) + 0()2), (21)

р(1) = р((/ь . .-ж .

дх ^ дх

где tj '^ < t. Применяя далее методы разложения дифференцируемых функций в степенные ряды, для (21) получим

T(tj+1) ~ Тз (tj) +

avT +—| Kxx | + i - Qx yz^ dx l ~ dx 1 1 x

<(tj+1 _tj) + 0 ((+1t, )2 )•

(22)

Теперь сделаем замену ( (t,) в квадратных скобках (22) на (р3 (tj)

(tj+1) ~(3(t,) +

д

дх

д(

ayz(3 + т-| Kxx-:— | + 1 'Q

дх

<tj+1 _tj) + 0 ((+1t, )2 )•

(23)

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

(2(tj +1) ~T(tj +1) +

X(tj +1 _ tj ) + 0 ( +1tj)2) (3 (tj+1 ) ~(2(tj+1) +

_Vy д( + ±[ Ky ( | + i - Qx

y dy dy ( xy dy 1 2 x

_Vz ( + ±(Kxz ( | + i3 - Qx

z 5z 5z ( xz 5z | 3 x

(24)

(25)

х(^.+1 - ^) + О ()2 ). Теперь делаем последовательную подстановку +1) из (23) в правую часть (24), а затем р2 ) в правую часть (25). В итоге

(3(tj+1) ~(3(tj)"

ayz(3 + Kxx ( | +

dx

dx

+1 - Qx _ Vy д( + ^|Kxy^ | + 1 - Qx _

y dy dy ( x dy

V (+-fKxz ( |+13 - Qx

z dz 5z ( xz 5z | 3

t=tj -(tj+ _tj) + О(+1tj)2)• Последнее выражение переписывается в виде

P3(tj+1) - Фз (tj ) (tj+1- tj)

3_( дфз дх I хх дх

ауФз +—| Kxx— | +

ду ду 1 у ду

^ •Qx - Vy P+d-l K:у p\+«2 •Qx - (26)

dP3 , д (г дФз

-V"b+dz 1 Kxz ir^ • Qx

Заменяя далее в (26) М) на ^(t, М), учитывая условие со\ + ю2 + + со3 = 1 и устремляя Д^-0, получим уравнение

Ф = аугф + Оф + дх, (27)

в котором дифференциальный оператор имеет вид

о = -Уу ф -Уг ф +АГк„^ФФА+Цк ^ФФА+ЦКхг Ф.

у ду г дг дх { хх дх ) ду { ху ду ) дг У хг дг )

Уравнение (27) полностью соответствует исходному (8). Таким образом, найденная функция ф(., М) в соответствии с изложенным алгоритмом действительно дает решение (8), а, стало быть, и приближенное значение первой компоненты искомого вектора скорости ветра Ух(Р,^+1). Аналогично доказывается состоятельность последующих вычислений для уравнений (12) и (16), ведущих к оценке функций Уу(Р,^+1) и Уг(Р,^+1).

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

К 1 по особенно с учетом дополнительных погрешностей, связанных с локальной линеаризацией задачи. Один из простейших способов коррекции прогноза: полученные значения Ух(~ 1+1-1, Уу(1+1), У/1+1) считать некими предварительными оценками указанных величин. Обозначим их

~ (1+1) _ _

через V х , Уу( 1+1), У2 (1+1). Вычисления теперь можно вновь повторить, введя в вычислительную схему вместо значений Ух 1+1), Уу^1+1), Уг(+1)

новые, такие как 2 У(1) + 1+1) ], | У(1 + Ру(1+1) ], 2 У(1 + У(1+1)

соответственно. Поскольку исходные уравнения в системе (1) представи-мы в виде ф = б\у ]ф + б, где б\к ] - дифференциальный оператор, зависящий от вектора V, а роль ф играет одна из искомых компонент этого

вектора в порядке следования уравнений в (1), с учетом указанной выше замены имеем конечно-разностное уравнение вида

р(j+1) -р(j)

-2 [ D( j) + D(j+1) ]р( j) + Q( j). (28)

Аt 2[

Подобная коррекция аналогична той, которая используется при решении дифференциального уравнения вида у' = /(х,у) первого порядка методом Эйлера. В соответствии с этой коррекцией в качестве решения Ук ±1 принимается значение, определяемое формулой

Л +А- Л = 2 [[ хк, Ук) + /(хк+1, -Ук+1)], (29)

где ук+1 = ук + Ах ■ /(хк, ук). Этот прием, как известно, позволяет достичь

второго порядка точности по Ах. Структурная аналогия формул (28) и (29) разумеется еще не доказывает того утверждения, что и (28) во всех случаях имеет порядок не ниже второго по Аt. Однако бесспорно, что эта про-

„ „ ./3+1) ./>)

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

Еще один аспект, связанный с практическим использованием рассмотренного выше алгоритма, состоит в том, что коэффициенты ауг, аХ2 и аху, перечисленные в порядке следования уравнений в системе (17), таковы, что для каждого момента ^ их сумма должна равняться нулю. Указанные коэффициенты, будучи функционально связанными, увеличивают взаимную обусловленность системы. Поскольку линеаризация в предлагаемой вычислительной схеме осуществляется, прежде всего, через эти коэффициенты, то о «качестве» получаемых приближений, скажем Ух((з+1), ¥у(}+1), Уг(з+1), можно судить по величине

a(Z+1)+ a(3+1) + a(j+1)

которая естественно не должна превышать некой достаточно малой величины е. Значение е согласуется тем или иным способом с порядком аппроксимации по переменной t, принятой в вычислительной схеме.

Литература

1. Марчук Г.И., Дымников В.П. и др. Мат. моделирование общей циркуляции атмосферы и океана. Л., 1984.

2. Ландау Л.Д., Лифшиц Е.М. Гидродинамика. М., 1988.

3. ЛайтхманД.Л. Физика пограничного слоя атмосферы. Л., 1970.

4. МарчукГ.И. Мат. моделирование в проблеме окружающей среды. М., 1982.

5. Марчук Г.И. Методы расщепления. М., 1988.

Северо-Кавказский государственный технический университет 20 апреля 2005 г.