CC BY
15
УДК 519.6
МНОГОЧЛЕНЫ БЕРНШТЕЙНА И МЕТОД НАИМЕНЬШИХ КВАДРАТОВ В ВЫЧИСЛИТЕЛЬНОЙ МОДЕЛИ УРАВНЕНИЯ НАВЬЕ-СТОКСА
© 2007 г. Н.И. Каргин, Р.А. Рыскаленко
Navier-Stokes equation is employed to ecological monitoring problems for planetary boundary layer. Its solution may estimate pole of wind speed with atmosphere turbulence. For the equation calculation algorithm are build by coordinate-wise method, least-squares method and with using Bernstein polynomials.
Построена вычислительная модель для векторного нелинейного уравнения Навье-Стокса, предназначенная для вычисления значений компонент поля скорости ветра в условиях пограничного слоя атмосферы. Подобные алгоритмы могут быть включены в соответствующие вычислительные схемы для решения задач нестационарного диффузного переноса загрязняющих примесей в атмосфере с целью последующего прогноза её экологического состояния. Уравнение Навье-Стокса, рассматриваемое в статье, предполагает пространственно-временную распределенность полей скорости ветра, турбулентности и силового поля, действующего на воздушные массы в пределах локального объема атмосферы. Предполагается, что вычислительная модель должна усваивать данные мониторинга, полученные в эксперименте. Для ее построения применялись методы покомпонентного и покоординатного расщепления, локальной линеаризации, аппроксимации на основе многочленов Бернштейна, а также вариационный метод наименьших квадратов. Вариационные методы, как известно, позволяют реализовать вычислительную схему при условии неопределенности исходных данных и при этом повысить устойчивость алгоритмов к погрешностям в них. Таким образом, решается проблема усвоения данных.
Рассмотрим основные этапы построения вычислительной схемы. Компоненты вектора скорости ветра V](Р, ^ , I = 1,3, где Р(^1,*2,хз)е ^ , удовлетворяют системе уравнений Навье-Стокса [1]:
3 V(Р t ) + x3 v. (Р, t ).3 V(Р ^
3 t
j=i
3 x.
(i)
—(F ( р, t) -Щ^ +mv2v ( р, t )).
р 3 Xi
i = 1,3 .
В (1) p(P, t) - атмосферное давление в точке P ; Fi (P, t) - компоненты силового поля F(P, t), действующего на единичный объем в пределах области W ; p (P, t) - плотность воздуха; m - динамический коэффициент вязкости; Vi (P, t) - неизвестные функции в (1), подлежащие определению. Для системы (1) полагаются справедливыми предположения [2, 3]: выполняется условие непрерывности divV = 0; Р (P, t) = const; (F3/р) = g ; F = F2 » 0; p = p R T ; R - универсальная газовая постоянная; T - темпера-
тура; Vp = р R' V T ; iL. у 2V = X
P j = 13 xj
a
/
K„
3 x,-
{Kj (Р, t)}
i, j = 1,3 - коэффициенты турбулентного
обмена.
При сделанных предположениях в рамках методов расщепления построена вычислительная схема [2, 3]:
3 V
3 V
3 V 3
an - a23v1 + v2"-n + v3u-n -J_ 31 3 X2 3 X3 3 X1
K
11
3 V1 3 x1
3 X2
/
K
3 V1
\
12
3 x2
3 V 3 x1 œ K 3 V2
K 22 . T-
3 x 2
3 x3
/
K
13
3 V1
3 x3
Qb
3 V2 + V 3 V2 V + V 3 V2 3 V1~- «13V2 + V3~- ^-
3t
J_
3 X2 3 V3
J_
3 x3 3 V3
3 x3 3 x1
œ K 3 V2
K 23 . 7-
3 x3
K
21
3V2 3 x1
= Ô2,
(2)
+ V1+ V2 3 t 3 x1 3 x2
3 X2
/
K
32
3V3 3 x 2
_3_
3 x3
a12V3
œ
3 x1
/
K
31
3V3 3 x1
K
33
3V3
3 x3
= Ô3:
^ 3 D — 3 V2 3 V3
где Qi = - (F - T*1), i = 1,3 ; «23 = T^- + T^, р 3 xi 3 x2 3 x3
3 V1 3 V3 «13 = —1 + —3
3 x1 3 x3
3 V1 3 V2 «12 = —1 + —2
3 x1 3 x2
В силу аналогичности структуры уравнений, входящих в систему (2), достаточно ограничиться построением вычислительной схемы для одного из них. В частности, в первом уравнении (2) обозначим через ^(Р, t) искомую функцию Р, t).
„ 3j1 т. 3j1 3
j 1 - «23 j 1 + V2 + V3 -
3 x 2 3 x3 3 x1
K
11
3 x1
3 X2
/
K
12
3 x2
3 x3
/
K
13
3 x3
= Qi .
(3)
Тогда для трехмерной задачи (3) в соответствии с методом расщепления имеем алгоритм [2, 3]: подзадача I:
'1 - «23 j 1
3
j1(P, tj ) =
3 x1 V1( P, 0)
K
11
1
3 x1
W1Q1,
если t = 0 (j |j 3(P, tj ), если j = 1,2,.
j1(P, tj ) = V1(P, tj ), если P £ W подзадача II:
0)
P £ W
tj £ t £ tj +1;
' 2 + V2
3 j 3 x2
2
3 X2
/
K
12
j 3 x2
= w 2 . Q1,
(4)
1
[j1(P,tjx 1), если Pе W
j 2 (P, tj) _\_ -
J [V1(P, tj), если Pе W
P е W , tj £ t £ tj +1; подзадача III:
;+ V3
д x3 д x3
Ki
13 '
д x3
= w 3 Ö1:
[j 2(P, tj +1), если P е W
] I У1(Р, tj), если Ре А р е А , tj £ t £ tj +1, ю 1 + ю 2 + ю 3 = 1.
В качестве решения уравнения (3) принимается
у (Р, :) з(Р, t) при tj £ t £ 1, ^(Р, tj+1) =
= У (Р
, : j +1). Аналогично вычисляются значения функций ^2(Р' tj + 1) и ^3(р, tj + 1) .
Основной задачей, решаемой в данной статье, является построение соответствующих вычислительных алгоритмов для уравнений, входящих в систему (4). На первом этапе реализуется так называемая процедура построения параметризованной модели (2), включающая в себя нормирование переменных и полей, входящих в нее; введение соответствующих параметров, обеспечивающих возможность их варьирования с учетом конкретной специфики решаемых прикладных задач.
Для построения параметризованной модели задаются начальные и краевые условия. Пусть переменные , Х2, Х3 и : принимают значения на отрезках
X е [0,X], Х2 е [0,Y], х, е [0,г], :е [0,Т]. Тогда начальные условия: Vг (Р,: = 0) = ¥[°\Р), г = 1,3 ; а граничные условия:
V- (хь Х2, Х3 = 0,:) = V(4)(хь Х2,:), Vг (Х1, Х2, Х3 = г,:) = V- (5)( Х1, Х2,:), vг (х1, х2 = 0, х3, :) = ^(6)(хь х3, :),
V- (Х1, Х2 = Y, Х3,:) = ^ (7)( Х1, Х3,:),
V (Х1 = 0, Х2, Х3,:) = V (8)( Х2, Х3,:),
V (Х1 = X, Х2, Х3,:) = V (9)( Х2, Х3,:). Нормирование переменных и распределений, входящих в систему (2), выполняется по формулам:
Х = Х1/Х, х2 = х21 Y, х = х^г, = :/т;
К(0)(Р) = у(0)(Р)/V*, V{к) = V(к)/V*,
V
max
P
V(0)(P), i= 1,3, k _ 4,9;
K (P, t) = K, j (P, t )/max K, j (P, t _ 0) ■ _ j _ 13.
P
F (P, t) _ F (P, t)/max F (P, t _ 0), ■ _ 13.
■p (P, t) _ p (P, t)/max p (P, t _ 0).
P
■p( P, t) = p(P, t)/max p( P, t = 0)
P
. В результате выполнения процедуры нормирования все переменные и
поля принимают значения из интервала [0, 1] . Все нормировочные коэффициенты - безразмерные величины.
В результате приходим к параметризованной модели, соответствующей системе уравнений (2):
<< V а V д у2 а V д V + а V д V +
—— - а2 ■ У1 ■ ---а3 -У1- --+ а2 ' У 2' Т-+
а: а Х2 а Х3 <т Х2
+ a 3 ' V3 ' - ß 11 д x3
l л /
l л /
12
д Х2
\ v
K
12
д Х1
д V öö
K
д V1
ö
11
д x
_1
д x1
ll
-ß
13
д x3
\ V
K
13
3V1
д x3
ö
F x д Р ö д x1
V*'T V*'T „ K11T K12T
a 2 = -V- , «3 = ß 11 ß 12 = 12
13
7 *
Ki3'T
z
2
g1
z X2
ii
F ' t x _ p T
* *, ? 1 * *
p V1 p ' X 'V1
(5)
д V2 + a V д V2 « V д V1 a V д V3 +
—— + « 1 ' V1 ' ---« 1 ' K2 ' ---a 3 ' V2 ' T-+
д t д x1 д x1 д x3
+ a 3'V3'dV2 - ß 21'
l
д x3
ll
ß 22'
l
д x2
\ v
K
22
д x1
д V2 ööV
K
д V2
ö
21
д x
_2
д x1
ll
-ß
23
д x3
\ V
K
23
д x3
ö
•j—r , дР
g 2'F2 -X 2 '
\
V t a 1 _ —-, a 3 _
23
X *
K23 T
д X2 ^ *
V3 T
z
ß _ K 21 ' T 0 _ K 22' T
ß 21 _ 2 , ß 22 _ -2—
X 2 Y 2
z
2
_ F2 'T X _ Р T
g 2 _ * X 2 _ * *
p * V2* p * Y V2*
V3
д V3 + a V д V3 + a V д V3 a V д V1
——+ a 1 ' V1 ' --+ a 2 ' V2 ' T--a 1 ' V3 ' T--
д t д X1 д X2 д X1
a 2'V3'^ - ß 31' д X2
l
l
ß 32'
l
д X2
\ V
K
32
д X1
д V3 ö ^
K
д V3
ö
\
д X2
•j—r , дp
g 3 ' F3 - X 3'T^
д x3 , *
V2* T
-ß
31
33
_3
д X1
ll
д x3
\ v
K
33
3V3
д x3
ö
V1* T a 1 _ —-, a 2 _
ß
33
X *
K33 'T
,2
g3
Y F3* T
ß _ K31' T ß _ K32 ' T ß 31 _ -ß 32 ■
X^
Y
X 3 _
p* T
* * , ^3 * ...
^ Р V р ■ X'V?
Применение многочленов Бернштейна для построения в^гчислительного алгоритма уравнения нестационарного переноса примесей в пограничном
3
слое атмосферы подробно изложено в [4]; там же даётся обоснование подобного подхода, а также приведены результаты численных исследований эффективности разрабатываемых алгоритмов. Этот же подход в данной работе, позволяющий систему дифференциальных уравнений (5) редуцировать к системе линейных алгебраических уравнений, применяется для построения вычислительных схем для векторного нелинейного уравнения Навье-Стокса. В общем случае многочлены Бернштейна имеют вид:
Bn(x) _ X f (Xk)P„k(x) _ X f (Xk)Ckxk(1 - x)n- k.
k _ 0 n- 1
k _ 0
B'n(x) _ n X [f(xk +1) - f(xk)\Chn_ 1xk (1 - x)n- k-1 k _ 0
1
(
■-П £ дР
g 1F1 -Х1'_г~
д x1
\
_Y 01, где
Y _ max
P
11 l 'p
■T- £ дР д x1
ö
д V2 д V3 'a23 _ a 2 ' т--a 3 ' -— . Тогда оно получит след x2 д x3
дующее представление:
д a V ß l д f K — - a23 ' V1 - ß11 ' -- K
11
_v1
д x1
ö
_W 1Y 01. (7)
< Х1
Далее положим, что:
3вм (х1;: ,Х1)» вм (хь,х1) - вм ('хьtj- 1,х1) <: Д: '
А: ® 0, е ,+1 ], tJ е [0,1], j = 0^.
В уравнении (7) аппроксимируем все поля представлениями, соответствующими многочленам Берн-штейна (6) следующим образом:
X 1 (Х1, tj ) = П(Х1,tj |-Х2,Х3) ,
м ,, ,,
К1(хь7/ |'х2;х3)» I X хт)'СМ''хЛ(1 - х^" т =
m_ 0
М
_ X f1(xm > tj )' PM, m (x1) _ BM (xb tj ,X1),
K11(x1; tj i x2, x3)» bm (x1> tj, K11), дKn(xb'tj |-x2,x3)
д x1
BM (xb tj, K11^
Vi (x1, tj | ■ x2, x3)» Bm (x1, tj, Vi),
01(x1,tj |"x2,x3)» Bm (xb'tj,01), i _ 1,2,3 :
д V1(xb-tj \x2,x3)
(6)
В"п (х) = п(п - 1) I [/(Хк + 2) - 2/(Хк + 1) +
к = 0
+ /(Хк )] СП- 2Хк (1 - Х) П- к- 2.
При этом полагается, что х е [0,1] /(х) = {/(Хк)}, к = 07П , /(х) » Вп (х) , т.е. Вп (х п®¥ /(х) , / (х) » В'п (х) и /"(х)» В"п (х). В качестве базисной функции выступает Рпк(х) = сП,хк(1 - х)п-к , для которой справедливы свойства Рп,0 (х = 0) = 1, Рп,0( х = 1) = 0, Р„п (х = 0) = 0, Рп, п (х = 1) = 1,
Рп,к(х = 0) = Рп,к(х = 1) = 0 . Применительно к системе подзадач (5) многочленами Бернштейна аппроксимируем все поля, входящие в них, полагая при этом, что соответствующие распределения представлены массивами дискретных данных. В качестве примера возьмем первое уравнение в (5). В нем обозначим
д x1
д 2V1( xb' tj \ x2, x3)
д x12 Тогда
BM(xbtj ,X1),
BM (xbtj ,X1).
д bm (x1;t ,x1) . B (x, f ) -Jt--a23 ' BM (x1,tj ,x1)_
- Ь 11 ' (Вм (Х1,, 'К11) В"м (Х1,, х 1) +
+ Вм(х1; tJ, ки) ■ Вм(х1,_ tJ, х 1)) = *
= ю 1? Вм (ХЬ tj, 01). (8)
Введем обозначения: ~(Х1) = 1 - Д: а23, Ь(хь-) = Д:'Ьи'В'м(Х1,,Кп),
~(Х1;) = Д: ■ Ь11' Вм (Х1,, кп), ~ *
d (Х1,) = А :■ ю 1? Вм (Х1,, 01) + Вм (хь'_ 1,^1).
С их учетом (8) преобразуется в систему линейных алгебраических уравнений относительно неизвестных
Iх 1(хт, tj, т = 0,м . С использованием заданных начальных и граничных условий х 1(х0, tj) и х 1(хм, tj): ~(Х1) ■ Вм (Х1,tj ,Х 1) - Ь (Х1,tj) ■ В'м (Х1,tj ,Х 1) -
- ~(Х1,tj) ■ В^ (Х1,,Х 1) = d(Х1,), или (6):
м
а(Х1) ■ I X 1 (Хт,) ■ Рм,т (Х1) -
т = 0
м - 1
~ М - 1 \
b (x1; tj)' X 1( xm + 1, tj ) -X1( x m' tj))'Pm- 1,m (x1) -
m_ 0
м- 2 1 ~
c(xb'tj)' X Ii 1(xml 2,tj)■ 2i1(xml1-'tj)tf 1(xm,tjPM- 2,m(x1)_ d(x1,tj)
m_ 0
м- 2 i ~
- ~(x1-ij)' I |f 1(xml 2/j)- 21(xml 1,tj)l f 1(xm>lj)|' PM- 2,m(x1): d(xb'tj),
m 0
m _ 0
которую далее можно преобразовать, если ввести обозначения:
Я0 (Х1,tj) = ~(Х1 )Рм,0 (Х1) + м ь (Х1, tj )Рм-1,0 (Х1) -
- м(м - 1)~(Х1,)Рм - 2,0 (Х1), Я1( Х1, tj) = ~(Х1) Рм ,1( х) +
+ м ■ ~(ХЬ- tJ ) [Рм - 1,1(Х1) - Рм - 1,0 ('Х1)] -
- м(м - 1)~(Х1, tj)[Рм - 2,1(х1) - 2Рм- 2,0 (х1)] ,
(Хь tj) = ~(Х1) Рм, ^ (Х1) +
+ м Ъ(x1;tJ)[Рм-1,* ('х1) - Рм-1 , * - 1 ( Х1)]-
п-
- M(M - 1)~(xbtj ) [PM_ 2,, (xi)-
- 2PM - 2, * - 1(x1) + PM - 2, * - 2(x1)] , 5 = 2, n - 2, gM - 1( xb tj ) = a (x1) PM, M - 1( x1) +
+ M b(x1ltj )[pM- 1,M- 1(x1) - PM- 1,M- 2(x1)] -
- M(M - 1)~(xb tj) [Pm - 2, M - 3( x1) - 2 PM - 2 (x1)], gm ('xbtj) = ~('x1)PM m ('x1) - m b ('xbtj )pm - 1,M - 1('x1) -
- M(M- 1)c(xb-tj)Pm- 2,M- 2(x0 , уравнение (8) примет вид:
M -1
i l(W/)'gmW,) + g(l(xb )X 1(4'Ji)+ gMШ/)i 1(xM- ) = d(x1 )
m= 1 M-1
i ' 1V i
I1 (xm-'/}gm(x1.i/)t?l(xb^1(х(Ьt/)t?M(lb'/1(xM>'/ d('xbУ-
ir1
(9)
Для нахождения неизвестных значений хт,tj)}, т = 0,М системы линейных уравнений (9) в соответствии с вариационным методом наименьших квадратов построим невязку £ (t j):
S (t,)=
1
d(xj,t, )
M -l s= 0
M
l X 1(xm>''/ )' gm (xs>''j ) - d (xs>')
m= 0
a stab
M
l (ii( Xs tj ) - X 1 (Xs - tj - !))
M X1, tj - 1) s = 0
( M
в которой
l (i1(xs -1, ) -i1(xs -1, - 1))
s = 0
(10)
является
стабилизатором, а а - параметром регуляризации в соответствии с терминологией теории некорректных задач. Построение регуляризирующего вычислительного алгоритма на основе выражения (10) позволит повысить его устойчивость к погрешностям в исходных данных. Далее, минимизируя невязку
V (tj) на каждом временном слое, получаем значения {|5 1 (хт >'tj )}М и V (Х1t Гх2,хз) =
компонент
M
= lXi(Xm -1 )CM'X1m (1 - Xi)
\M - m
соответственно.
m = 0
Для минимизации V (tj) в структуру решающего алгоритма включена процедура, соответствующая так называемому эвристическому симплексному методу минимизации п -мерного нелинейного функционала нулевого порядка [4].
Выполняя аналогичные построения для всех подзадач (4) применительно к каждому уравнению системы задач (5), в итоге получим алгоритм, в котором последовательно вычисляются значения компонент поля скорости ветра V, ^ и V?. Продемонстрируем часть алгоритма для первой компоненты V. Для первого уравнения системы (5):
Подзадача I:
Ы хЬ tj I' х2 = х1, х3 = Ч ) =
M ,,
l i 1 (Xm-tj ) CM -X1m (1--X!)
m = 0
M - m
X 1(xm , tj ) = V1(x1 = xm ,t = 0,\' x2, x3) , если j = 0, X 1(0)(xm , tj - 1) = V1(x1 = xm , tj - 1 lx2 ,x3 ) , если j > 0 , x1(x0>'tj ) = V1('x1 = x0>"tj \'x2, x3) , x 1(xM i tj ) = V1(x1 = xM ltj \'x2 ,x3 ) , ~(x1) = 1 - A t • a23, b (x1,tj) = A t • ß 11 • B'm (x1,tj, K11),
~(' x1; tj) = A t • ß 11 • Bm (xb" tj, K11), ~ *
d(x1,tj) = A t^ а 1Y Bm(xb-tj,61) + Bm(xb'tj_ 1), ~(x1) • Bm (xb" tj 1) - b(x1, tj) • B'm (xb" tj 1) -
- ~(x1, tj )• BM (xb' tj ,X1) = d (xb' tj), g 0( x1, tj) = ~(x1) Pm ,0( x1) +
+ m • , tj )pm _ 1,0 (' x1) - m (m -1)~( , tj p _ 2,0 (' x1)
g1( x1, tj) = ~(x1) Pm ,1( x) + + M • ~(x1, tj)[ Pm -1,1 (x1) - Pm -1,0 (' x1)] -
- M(M - 1)~(xb tj)[Pm - 2,1(x1) - 2Pm - 2,0 (x1)], g* (x1, tj) = ~(x1) Pm , * (xO +
+ M ~(x1,tj)[рм-1,*(x1)- рм- 1,*- 1(x1)]-
- M(M - 1)~(xb tj) [Pm - 2, * (x1) -
- 2Pm - 2, * - 1(x1) + рм - 2, * - 2(x1)] > * = 2, n - 2, gM - 1(x1, tj ) = a (x1)PM,M - 1 (x1) +
+ M•Й(xl; tj ) [рм _ 1,M - 1(x1) - PM - 1,M - 2(x1)]-
- M(M - 1)~('x1,tj) [Pm - 2,M - 3 (x1) - 2PM - 2,M - 2 ('x1)], gM ('X1,' tj) = ~ ('X1)Pm M ('x1) - M b ('x: tj )Pm _ 1,м-1' x1) -
- M (M - 1)~(xb- tj )Pm - 2,м - 2 (xO , M - 1
I X 1(xm,tj ) • gm (x1,tj ) + g0 (x1,tj )X 1(x0, tj ) +
m = 1
+ gM (x1tj )X 1 ( xm , tj) = d (x1, tj), i (tj) =
1
d(X1,t/)
M -l s= 0
M
l X 1(xm- 'j )' gm (xs-) - d (xs- 'I )
m= 0
a stab
M
l (i1(xs -'/ ) -i1(xs -- 1))
IM X1, tj - 1)| s = 0
inS (t/ ) ^ lX 1(Xm */ )\m-_ 0M
s' /-
min
0£ xj< 1
M
V1(x1-'/ I' x2 = xi ,x3 = xk ) = l X 1(xm -)'PM ,m (x1)
m= 0
m = 0,M , l = 0,L , k = 0, £ , -1 < 7 < '/ . Подзадача II:
L -L . l n ... vL-1
L l l = 0
V1(x2,''/ I'x1 = xi,x3 = xk) = l X 2(xl;'/)'Cj ''x2 '(1 - x2)
2
2
+
2
2
+
2
+
X 2 (xi' tj) = Vi(X2 = xi' t = 0,1"xb"X3), если j = 0,
X 20)(xl ' tj ) = V1(x2 = xl ' tj \X1 = xm 'x3 = xk ) при x1 = xm;x3 = xk , если j > 0,
X 2 (x0'' tj ) = V1( x2 = x0'' tj \ xb' x3),
X 2 (xL 'tj ) = V1(x2 = xL 'tj rxb'x3), ~(x2) =
¿Cx2,'tj) = a t • (b 12 • B'l(x2,tj,- a 2 • BL(Ыt}_ i,V2)),
C( x2, tj ) = A V b 12 BL (x2, tj , Г12), ~ *
d (x2;tj-) = A tw 2Y Bl (x2; tj, Q1) + Bl (x^tj. i,X 2), ~(x2) • Bl (x2tj, X 2) - b (x2, tj) • B^ (x2' tj, X 2) -
- ~(x2'' tj )B"l (x2,' tj 2) = d (x2,' tj ),
g0(x2,' tj ) = x2)Pl,0(x2) +
+ L b~(x2'tj )Pl- 10 (x2) - L(L - 1)~(x2'tj )Pl- ^x2)' x2' tj ) = x2)Py(x2) + + L~( x2' tj )[Pl- 1'1(x2) - Pl- 1'0(x2)]-
- L(L- 1)~(x2'tj)[Pl-2'1(x2)- 2Pl-2'0(x2)]' gs (x2' tj ) = x2)^^L'S (x2) +
+ L b (x2' tj) [Pl- 1'S (x2) - Pl- 1^-1 ( x2)]-
- L(L - 1)c(x2'tj) [Pl- 2'S (x2) -
- 2PL- 2'S- 1(x2) + PL- 2,s- 2 (x2 )] ' S = 2, n - 2 , gL- 1(x2' tj ) = x2) PL,L-1( x2) +
+ Lb(x2\tj) [Pl- 1'L- 1(x2)- Pl-\,l-2(x2)]-
- L(L - 1)c(x2' tj ) [PL- 2'L- 3(x2) - 2PL- 2'L- 2 (x2)] ' gL (x2'' tj ) = a(x2)PL,L (x2) - ¿'b fe tj )Pl- 1(x2) -
- L(L- 1)c(x2'tj)Pl_2'L-2(x2).
L-1
I X 2 (xi' tj )g fe; tj) + g 0 (x2' tj )X 2 (x0' tj) +
l = 1
+ gL (x2' tj)X 2( xL' tj) = d tj).
i (tj) =
1
ii2
L ~
I X 2(xl'tj)' gl(xs' tj) d (xs' tj)
^ tj )|| a stab
s= 0
\2
l = 0
X 2(x2' tj -1)||2 st 0"ZV " J
0^x2^ 1? (tj ) ^ 2(xl'tj )ll=
V1(- x2 ' tj [ x1 = xl' x3 = xk ) =
-iL
I (X 2(xs' tj ) -X 2(xs ' tj - 1))
I X 2 (xi ' tj ) CL -x2l (1 - x2 )L-1 l= 0
m = 0' M ' l = 0' L ' k = 0' £ ' ^j -1 £-t £ tj Подзадача III:
Vl(Х'З' tj \'x1 = xm 'x1 = xl) =
£
= I Xз(xk'tj-)• c£ -x3 • (1 - x3) k = 0
£ - k
X3(xk'tj) = V1(x3 = xk't = x1'x2)' если j = 0'
X 3(0)(xk'tj ) = V~1(x3 = xk' tj \'x1 = xm 'x2 = xl) при ■x1 = xm'x2 = xl' если j > 0'
X 3 (x0'tj ) = V1(x3 = x0'tj \-x1'x2) '
X 3 (x£'tj ) = V1(x3 = x£ 'tj Гxl;x2); ~(x3) = 1'
b^j) = A t• (b 13• B£('x3'tj'£1()- a3 B£(x(' j 1'V3))'
c( x(' tj) = A t b1( B£ (x(' tj, £13)' ~ *
d (xd tj) = A t^ w 3Y B£ (xd tj, Q1) + B£ (x(' tj _ b£ ();
x3) • B£ (x3'tj X 3 ) - b(x3'tj ) • B£ (x3'tj X 3) -
- x3' tj ) B£ (x(' tj 'X 3) = d (x(' tj)' g0(x3'tj ) = x3)P£'0(x3) +
+ £ • x(' tj )P£_ 10 (x() - £ (£ - 1)c(x(' tj )P£_ 2 0 (x()' g1( x3' tj ) = x3) P£ '1( x3) +
+ £ • x( ' tj )[P£_ 11 (x( ) - P£_ 10 (x()]-
- £(£ - 1)~(x('tj)[P£- 21 (x() - 2P£- 2'0 (x()]' gs(x3'tj ) = «(x3)PL;S(x3) +
+ £~( x^ tj )[p£- 1'S (x() - P£- 1'S- 1(x()]-
- £(£ - 1)C(x('tj ) [P£_ 2'S (x() -
- 2P£- 2'S- 1(x() + P£- 2'S- 2(x()] ' S = 2'П - 2 ' g£ - 1 (x3 ' tj ) = x3) P£' £ - 1( x3) +
+ £ ~(x('tj )[p£- 1(x() - P£- !'£- 2(x()]-
- £(£ - 1)C(x('tj )[P£- 2'£- 3(x3) - 2P£- 2 (x()]' g£ (x3') = x3)P££ (x3) - £ b ('x3;)P£- 1£- 1(x3) -
- £(£ - 1)C(x('tj )P£_ 2'£- 2 (x()' £ -1
I X 3 (xk' tj) • gk(x3' tj ) + g0 (x3' tj )X 3 (x0' tj) + k = 1
+ g£ (x3' tj )X 3 (x£' tj ) = d (x3' tj ) ' i (tj ) =
1
К ( к
I
' s= 0
a stab
I X 3 (xk'tj )^gk (xs'tj ) - d (xs'tj )
k = 0
£
I (X 3 (xs' tj ) -X 3 (xs't/-1))
s= 0
3 (x3' tj_ 1J min (tj) ^ {xз(xk'tj)\k__ —'
s j -
0£ x3 £ 1
£
V1(x3' tj ^x1 = xm ' x2 = xl) = I X 3(xk' tj )•P£; k (x3)
k = 0
m = 0'M ' l = 0' L ' k = 0' £ ' ^j-1 £ _t £ tj ; w1 + w 2 + w 3 = 1'
M
V1(x1' t ^з) = I X ((xm' t )• CM -x1m (1 - x1)
m = 0
при x2 = xl' x3 = xk ' ^j- 1 £ 7 £ tj .
M - m
2
+
+
t
2
2
+
Для вычисления двух других компонент поля скорости ветра алгоритм строится аналогично.
Рассмотренный выше алгоритм был реализован программно, и получены соответствующие результаты. Для проведения численных исследований в работе выполняется постановка вычислительного эксперимента, включающая в себя две части. Создается тестовый пример, в котором задаются исходные данные и на их основе моделируются все распределения, входящие в исходное уравнение Навье-Стокса, а также генерируется точное решение. На основе тестового примера выполняется тестирование вычислительного алгоритма и проводятся численные исследования свойств алгоритма, такие как сходимость и устойчивость, и затем исследуется зависимость поля скорости ветра от турбулентного состояния пограничного слоя атмосферы, определяемого коэффициентом турбулентной диффузии.
Тестовый пример строится следующим образом. Введем функции вида: иу(ху) = соб^х + су),
и2(Х2) = СОБ(Ь2Х2 + С2), из( Х 3) = С05(Ьз'Х 3 + С3),
и4(г) = cos(b4t + с4) ,0'(Х') = бш^х' + Ь'), в 2(Х2) = БШ^х + Ь2), 0'(х3) = БШ^х3 + Ьз), в4(г) = sin(a4t + Ь4),где Ьг , сг, г = 1,4 - произвольные константы. Положим далее, что значения компонент скорости ветра:
^2 (Х1;Х2 ;Х3;t ) = ^0 • и' (Х' и (Х2 )и3 (Х3 )U4(t) , V3 (Х' ,Х2 ,Х3,t ) = ^30 • в 1 (Х' )0 2 (Х2 )0 3 (Х3 в 40 ) , У'( Х', Х2, Х3,' t) =
x1
а V2 (xj; Х2; Х3 ,Г) + Щ (xj; Х2; Х3; t)
a x2
а x3
dx1 . (11)
Последнее уравнение (11) удовлетворяет условию
а V + а V2 + а V3
непрерывности —— + —— + —— - о , из которого сле-
а Х1 а Х2 а Х3
дует, что выполняются условия
avL а Х1
/
а v2 а v3
+
а Х2 ахз
av2 а х2
а V, а V
+
а Х1 ах-
ÍVl
а х3
а Vi а v
+
а Х1 ах-
Здесь и далее по тексту имеется в виду, что значения переменных и функций нормированные. Аналогично определим функции h11(t) = cos2(euT + dn),
2 2 h 12(t ) = COS (e^t + ^12) , h 13 (t ) = COS (enT + dn) '
22 h 21 (t) = cos (e21t + d21), h 22 (t) = cos (e22t + d22) ,
22 h 23(t) = cos (e23t + d23), h 31(t) = cos (e^T + d31),
22 h 32 (t) = cos (e32t + d32b h 33 (t) = cos (6331 + d33^
d, -
^ ij (t) = sin2 (ejt + dij), i = 1,3, j = 1,3,
константы; w (xj; X2; X3, t) =
= (1 + a- cos(b Xj + c))• (1 + a• sin(b X2 + c))•
■ (1 + a - sin(b • 'X3 + c))- (1 + a - cos(b • 't + c)),
с помощью которых далее распределения турбулент ного обмена, плотности воздуха и атмосферного дав ления примут вид
Kij (Х1,'х2 'х3,t) = Kj,o ' h ij (X1 )h j (x2 )h j (x3 )£ j (t) , i, j = 1,3, p (x1;x2;x3,t) = p 0 = const,
p(xb'x2;x3,' t) = Po • w(X1'X2;X3,' t). (12)
Неизвестными считаем значения функций компонент силового поля X1,'X2,'X3,' t ), F2(X1,X2'X3,' t), F3( X1; X2' X3' t). Для их определения подставим значения исходных данных в левую часть первого уравнения (5) и вычислим правую часть обозначая ее при этом "61 (X1' X2' X3'' t). Предварительно вычислим значения соответствующих производных' входящих в первое уравнение (5) путем непосредственного дифференцирования выражений' введенных выше. Функции' соответствующие производным: а Vk (Xb-X2; X3I t)
gtk (X1,X2;Хз,t ) -g ik(x1, x2; x3; t) -
gg ik(x1, x2; хз; t)-
rkik ( x1, x2, x3; t) -
а t
а Vl (Х1; Х2; Х3; t) axk '
а 2Vi (Х1; Х2; Х3; t)
а xk2 а Kik (Х1; Х2; Х3; t)
а xk
, ч ар(х',хо,'х3,t) — Чг (Х',Х2,Х3,t) = 1 а 2 л-, г, к = 1,3.
а х,
В этих обозначениях первое уравнение (5) принимает вид gtl - (а 2gg22 + а 3ggзз)Vl + а 2Кгgl2 + + а 3^gl3 - Ь 1lgl1^k11 - Ь 12gl2гк12 - Ь 13gl3rк13 -
- Ь 11K1lggl1 - Ь 12К12ggl2 - Ь 13К13ggl3 = х1,х2,х3;I) .
В результате, вычислив значения функции
ОК-х'^-хз,-г), получим Sl('xl;x2;xз;t) = у^! 1^1 - (1 • Ч1),
откуда имеем = — (£ г Ч1 + 61 • р ). !1
Аналогично вычисляем компоненты
р2 (ху, х2,хз, г) и рз (ху, х2, хз, г):
gt2 - а 1Vlg21 - (а lggl1 + а 3g33)V2 + + а 3V3 Я 23 - Ь 21^218821 - Ь 22К 22 ЯЯ 22 - Р 23К 23 ЯЯ 23 -
- Ь 21Я21гк21 - Ь 22Я22гк22 - Р 23Я23гк23 = 62('xl;x2;'x3;г )
F2 - —2'92 + Ö2'P) . 7 2
(13)
gt3 1 a 1V1g31 1 a 2V2q32 - (a 1gg33 1 a 2gg22 )V3 -
- ß 31K31gg31 - ß 32K32gg32 - ß 33K32gg32 -
- ß 31g31rk31 - ß 32g32rk32 - ß 33g33rk33 _ 03(x1,x2,x3>t) ,
F3 _ —' (f3' q31 Q3' p ) . g3
Далее каждой точке P( x^' x2,' x3, t) ставится в соответствие точка на сетке P(xm, xl, xk > tj), m _ 0,М
, l _ 0, Z , k _ 0, K , j _ 0, N, и проводятся вычисления значений соответствующих сеточных функций по формулам (11)—(13). Таким образом, значения всех массивов исходных данных считаются заданными, а значения компонент скорости ветра далее будем считать точными решениями, которые
обозначим Vi (x1;x2;x3,t) _ VT (x1;x2;x3,t), i _ 1,3 .
Дальнейшая реализация вычислительного алгоритма проводится в предположении о том, что
компоненты Vi(xi;x2,'x3'J), ^'2('х1-'х2-'хз-О. Т з('х1;х2,'х3,7) поля скорости ветра неизвестны и должны определяться в ходе вычислений согласно ре-гуляризирующему алгоритму, описанному выше при заданных исходных данных. Точность вычислений при этом оценивается значением величины отклонения приближенного решения Т](Xi:X2,'X3'J ) ОТ ТОЧНОГО Tf (х1,-х2,-х3,7) - / = Гз по формуле:
s i _ ViT (x1;x2;x3; t) - Vi (x1;x2;x3; t)
Описанный выше тестовый пример реализован при следующих значениях исходных данных:
1,0 = 6.33 м/с, f'*2.о = 5 м/с, f' з.о = 3 м/с, Kh0 = 108,72 м2/с, К2_0 = 108,72 м7с, Кхо = 108,72 м2/с, А'= Г= Z = 1916,52 м, Т = 300 с. Л/ = Z = К = 10 , j=l20 a _ [1,8 0,9 - 0,5 5\, b _ [- 1 - 2,1 0,05 7\,
Рис 1. Графическое представление точного решения для фиксированных момента времени и высоты
vi
0.06 004 002 О
402 404 406 408
Рис. 2. Графическое представление приближенного решения 1\ (х^, х2 | х3, Г), полученного при й = 0 , для фиксированных момента времени t = 30 с и высоты х3 = 192 м
с = [0Д 0'2 3 0'8]' etj = 0'2 ' dij
0'02
у - - и
На рис. 1-3 показано пространственное распределение компоненты поля скорости ветра
ViT (xi ,х2 | x3,i) - точного решения и приближенных решений l\(xi,x2 | x3,i) для двух случаев: параметр регуляризации a slui1 = 0 (рис. 2), и в результате получаем некорректное решение, сильно отклоняющееся от точного; a, sla!-, - 0,19 . Графическое представление не имеет всплесков (рис. 3).
В последнем случае относительная погрешность
vi
1
V1 (xm' xl' xk' tj) - V1( xm' xl' xk' tj)
|M' L' £ составила 0J2.
V1 (xm' xl' xk' h)
(14)
Рис. 1 рафическое представление приолиженного решения, полученного при й = 0,19 , для фиксированных момента времени : = 30 с и высоты Х3 = 192 м
2
Для г - 30 с и хз = 192 1 в таблице приведены значения параметра регуляризации и относительной погрешности (14). Очевидно, что в ходе вычислительного эксперимента можно подобрать некоторое оптимальное значение а *6-гаь , обеспечивающее минимальное V ш;п при заданных значениях исходных данных.
Влияние параметра регуляризации а ¿¡аЬ на относительную погрешность V
a stab s
0,08 0,2451630
0,10 0,2205225
0,15 0,1644680
0,16 0,1540130
0,17 0,1452507
0,18 0,1389044
0,19 0,1259973
0,20 0,1268064
1,0 0,7195334
На рис. 4, 5 приведены результаты численных исследований влияния величины коэффициента турбулентной диффузии на поле скорости ветра.
Рис. 5. Линии уровней, соответствующие графическому представлению приближенного решения ху,х2 | хз,г) первой компоненты поля скорости ветра, полученной при
значении коэффициента турбулентной диффузии Кц(ху,х2 | хз,г) - 120 м2/с при фиксированном моменте времени г - 30 с на высоте хз - 192 1
Анализ полученных результатов показывает, что при увеличении коэффициента турбулентной диффузии появляются четко выраженные (особенно на рис. 5) вихри. Структура поля при увеличении коэффициента турбулентности существенно меняется.
Таким образом, разработанные вычислительный алгоритм и соответствующий вычислительный эксперимент открывают возможности для проведения численных исследований и моделирования процессов, протекающих в турбулентной атмосфере в пределах пограничного слоя.
Литература
\.Марчук Г.И. Математическое моделирование в проблеме окружающей среды. М., 1982.
2.Каргин Н.И., Наац В.И. II Изв. вузов. Сев,-Кавк. регион. Естеств. науки. 2005. Приложение. № 5. С. 3-13.
3..Каргин Н.И., Наац В.И. // Информационные технологии в науке, проектировании и производстве: Тез. докл. в материалах
первой компоненты поля скорости ветра, полученной при
,, ^ - л.л. XV Всерос. науч.-техн. конф. (июнь 2005 г.
значении коэффициента турбулентной диффузии
. , Н. Новгород). Нижний Новгород, 2005. С. 3-4.
К\\(х\, хт | хз, г) - 60 м2/с при фиксированном моменте л тт п тж ,, тг т
11 1 ^ -3 г 4Наац В.И. // Компьютерное моделирование 2004: Тез.
времени г - 3° с на высоте хз - 192м докл. в трудах 5-й Междунар. науч.-техн. конф. (29 июня-
3 июля 2004, г. Санкт-Петербург). СПб., 2004. С. 128-129.
Рис. 4. Линии уровней, соответствующие графическому представлению приближенного решения V'(ху,х2 | хз,г)
Северо-Кавказский государственный технический университет, г. Ставрополь
1 декабря 2006 г.
