СИСТЕМЫ И
Л. \
A. I
УПРАВЛЕНИЯ —
УДК517.9+532.5
ОБ ОДНОМ ПОДХОДЕ К МАТЕМАТИЧЕСКОМУ МОДЕЛИРОВАНИЮ ПЛОСКИХ СТАЦИОНАРНЫХ ТЕЧЕНИЙ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ В КОНЕЧНЫХ ОДНОСВЯЗНЫХ ОБЛАСТЯХ
ТЕВЯШЕВ А.Д., ГИБКИНА Н.В., СИДОРОВ М.В.
Рассматривается применение метода R-функций в сочетании с методом последовательных приближений для решения задачи расчета плоских стационарных течений вязкой несжимаемой жидкости в конечных односвязных областях. Доказывается сходимость построенного итерационного процесса в норме пространства w|(^) к обобщенному решению исходной задачи. Получены оценки скорости сходимости. Предложенный метод протестирован на модельных областях, полученные приближенные решения сравнены с решениями, полученными другими авторами.
Введение
Актуальность задачи. Изучение законов движения жидкости играет важную роль в развитии техники и естествознания. Исследования в этой области стимулируются потребностями авиации, кораблестроения, теплоэнергетики, геофизики, биологии и пр. За последние десятилетия сфера исследования и применения явлений, связанных с движением жидкости, постоянно расширяется и охватывает ведущие направления промышленности (химические технологии, нефте- и газоразработка, металлургия и т.д.) и ряд естественных наук (биология, физика атмосферы и океана и др.). Во многих практически важных случаях жидкость можно с большой достоверностью считать вязкой несжимаемой ньютоновской средой, и происходящие в ней процессы могут быть промоделированы с помощью уравнений Навье-Стокса [1, 2]. Различные задачи, возникающие при изучении динамики вязкой жидкости, могут быть исследованы теоретическим путем или с помощью физического эксперимента. Однако с развитием ЭВМ все активнее используется математическое моделирование. Существует множество численных методов, применяемых при расчете вязких течений. Литература по этому направлению обширна [3-5 и др.]. В основном эти численные методы используют метод конечных разностей и метод конечных элементов. Они просты в реализации,
но не обладают необходимым свойством универсальности - при переходе к новой области (особенно неклассической геометрии) необходимо генерировать новую сетку, а часто и заменять сложные участки границы простыми, составленными, например, из отрезков прямых. Точно учесть геометрию области можно воспользовавшись конструктивным аппаратом теории R-функций, разрабатываемой акад. В. Л. Рвачевым и его учениками [6, 7 и др.]. Задачи гидродинамики решались в работах С.В. Колосовой, К.В. Максименко-Шейко, И.Г. Суворовой, Т.И. Шейко и др. [8-11], однако в основном рассматривались задачи динамики идеальной жидкости или вязкой для случаев, когда можно построить решение за счет удачного выбора координат (осесимметрические течения, течения, обладающие винтовой симметрией, и т.п.). Поэтому разработка новых, а также совершенствование существующих методов математического моделирования динамики вязкой жидкости на основе метода .^-функций является актуальной научной проблемой.
Цели и задачи исследования. Целью данной работы является создание современного и эффективного метода математического моделирования плоских стационарных течений вязкой несжимаемой жидкости в конечных областях неклассической геометрии с кусочно-гладкой границей. Основные результаты по теоретическому обоснованию корректности начально-краевых и краевых задач для уравнений Навье-Стокса получены О.А. Ладыженской [12], Ж.-Л. Ли-онсом [13], Р. Темамом [14]. В данной же работе не обсуждается степень строгости, условия применимости используемых уравнений движения жидкости. Они рассматриваются как математические модели, подлежащие численной алгоритмизации. Для достижения поставленной цели решаются следующие задачи: разработать и обосновать метод расчета течений Навье-Стокса в односвязных областях в переменной «функция тока»; в целях обоснования эффективности разр аботанного метода решения задачи р асчета течения вязкой жидкости применить его для решения модельных задач при различных числах Рейнольдса.
1. Постановка задачи
Плоское стационарное течение вязкой несжимаемой жидкости в переменных «скорость - давление» описывается известной системой уравнений Навье-Стокса [1, 2]:
Л, 5V1 1 5Р f
vAVi — Vi----v V,----і------f1
1 1 SX 2 ЗУ p 3X 1:
лл, _л, 5V2 л, 3V2 1 ЭР .
vAVo — V1---v V9----1------fo
2 1 dX 2 5Y p ЭУ 2
и уравнением неразрывности:
^Vl +^vl = 0
dX дУ .
(1)
(2)
(3)
50
РИ, 2007, № 2
Здесь V = {Vj, V2} - вектор скорости жидкости, р -давление, f = {fj, f^} - вектор массовых сил, v -кинематический коэффициент вязкости, р - плотность жидкости.
Далее будем рассматривать безразмерные перемен-
ные Vj =
V1 v V2 X Y
U0 , V2 - U0 x = —, L y = —, где L
0
характерная скорость, l - характерная длина. В двумерном плоскопараллельном течении уравнение неразрывности (3) интегрируется с помощью функции
тока y(x,y), определяемой соотношениями
Зу Эш
Vj =%, V2 =-£
Тогда, введя завихренность
(4)
„ 3v2 9vi
с=—17. <5)
систему (j) - (3) в случае, когда вектор массовых сил потенциален, можно свести к системе
2. Выбор и обоснование метода решения
Рассмотрим плоскую стационарную задачу о движении вязкой несжимаемой жидкости в конечной односвязной области q с кусочно-гладкой границей 3Q. Такое течение может быть описано краевой задачей для одного нелинейного уравнения четвертого порядка относительно функции тока y(x,y):
Д2у = Re J(Ay, у) в Q.
Мап= f (s), $
3n
= g(s), s є dQ..
(9)
(10)
5Q
Здесь Re - число Рейнольдса;
ЭДу Эу ЗДу Зу
J(Ay, у) =-
3x dy dy Эх
df
ds ,
g - некоторые распределения нормальной и
касательной составляющих скорости потока соответственно; n - внешняя нормаль к 3Q. Предполагаем,
что f є w2(3Q) , g є L2(3Q).
_дс+Re|3vas_svaC| = 0
1 3y Эх Эх cyy
-Ду = С,
Для решения задачи (9), (10) построим процесс пос-
(6) ледовательных приближений.
(7) Обозначим через u0(x,y) решение следующей задачи:
где Re = -U°L
V
число Рейнольдса.
Система (6) - (7) уже не содержит давление, однако в общем случае не удается корректно задать граничные условия для завихренности Q на твердых стенках и при решении задач эти условия приходится задавать приближенно [2]. Этого можно избежать, перейдя к эквивалентному нелинейному уравнению четвертого порядка для функции тока
д V Re (аАуЗу_ЗАуЗуї 0 ^ 3y 3x dx dy j
(8)
Граничные условия для функции тока могут быть получены из условий, накладываемых на вектор v. Если жидкость примыкает к неподвижной стенке, то в этих точках скорость жидкости обращается в нуль. Это означает, что в нуль обращается нормальная и тангенциальная составляющие скорости (условие прилипания). Если же жидкость примыкает к подвижной твердой стенке, то в таких точках скорость жидкости должна по величине и направлению совпадать со скоростью соответствующей точки стенки. Таким образом, на границе области течения можно задать значение функции тока V и ее нормальной производ-
Эу
ной — .
3n
Д2u0 = 0 в Q ,
(11)
U01Xi= f (S),
= g(s), s e9Q . (12)
5Q
Задача (11), (12) в сделанных предположениях может быть решена с помощью метода, базирующегося на
методе R-функций [15-17], при этом U0(x,y) є W2(Q).
В задаче (9), (10) сделаем замену у = U0 + u , где u -новая неизвестная функция. Это приводит к задаче
Д2u = ReJ(A(u0 + u), u0 + u) в Q. 3u
uШ 0 ,
3n
= 0 .
дП
(13)
(14)
Последовательные приближения к решению задачи (13), (14) будем строить следующим образом. Пусть
начальное приближение u(0) задано. Тогда при известном значении u(k) функции u на k-й итерации следующее (k +1) -е приближение находится как решение линейной задачи
д 2u(k+1) = ReJ(A (u0 + u
,(k+1)
5Q
= 0,
3u(k+1)
3n
дП
(k)), u0 + u(k)) в Q , (15) = 0, k = 0, 1, 2, .... (16)
РИ, 2007, № 2
51
Для дальнейшего изложения нам потребуются некоторые свойства якобианов.
2
Лемма. Пусть u, v, w є W 2(Ц), тогда
1) J(u,v) = -J(v,u);
2) J(Ui,Vi) - J(U2,V2) = J(U2,V1 - V2) + J(Ui -U2^) ;
3) если U є w2(Q) , то J J(U,V)dxdy = 0 ;
Q
о
4) если w є W 2(Ц), то
J J(U,V)wdxdy = J J(V,w)Udxdy ;
Q Q
5) если U є w2(Q), V є w2(Q) , то
(u,v)IIl2(Q) - CPllw2(Q)lrllw3(Q) .
(17)
жим, что
,(j)
(k+1)
W3(ri) -c1Re J(A(u0 + u(k))’ u0 + uW) W 2^)
,(k).
L2(Q) .
Применяя неравенство (17), получаем ||u'
(k+1)| 3 < c^Rell A(U0 + U(k))|| 1 ||u0 + u(k)|| 3 <
IIw|(Q) 1 2 II 0 llwi(Q)ll 0 Hw3(
(
< C1C2Re
IKIlw^Q) +
W2(Q) (k)
IIW2(Q)
W^(^)
Пусть ||U0IW3(Q) ^ M0, тогда
, (k+1)
w|(Q) S C1C2Re(M0 + M)
Условие ограниченности полнено, если
Re <-
,(k+1)
M
c1c2(M0 + M)
2.
Докажем сходимость последовательности U(k) , k = 0, 1, 2, .... Рассмотрим разности
SU(k+1) = U(k+1) - U(k). Они удовлетворяют уравнению A2Su(k+1) = Re[J(A(u0 + U(k)), U0 + U(k)) -
-J(A(U0 + u(k_1)), U0 + u(k_1))j в Q и краевым условиям
58u(k+1)
Su(k+1)
= 0,
5Q 9n
= 0, k = 0, 1, 2, ....
дП
Тогда
Su(k+1)
W3(ri) ^ c1Re J(A(u0 + u(k)X u0 + uW) -
W 2(““)
,(k))
Доказательство леммы основывается на теоремах вложения [18].
Изучим вопрос сходимости итерационного процесса (15), (16). На каждом шаге этого процесса нужно
о
найти функцию u(k+1) є W 2(Ц) П W 2(Ц). Предполо-
-J(A(U0 + u(k-1)), U0 + u(k-1))T (n)
l2(Q) '
Используя свойство 2 из леммы, получаем
J(A(U0 + U(k)), U0 + U(k)) - J(A(U0 + U(k_1)), U0 + u(k_1)) = = J(A(U0 + u(k_1)), 8u(k)) + J(A8u(k), u0 + u(k_1)). Тогда
3 < M для всех j = 1, 2, ..., k.
W^Q)
Выясним, при каких ограничениях на число Рейнольдса Re итерационный процесс (15), (16) будет давать ограниченное в норме W^(^) решение u(k+1). Имеем
К І
< CiC^Re
W^(O)
3_ - ^2™I ||A(U0 + u(kI8u(k)|w3
A8u(k)
< 2C1C2Re
u0 + u
W2(Q) (k-1)
u0 + u
W^(Q)
(k-1)
W2(Q)H HW^Q)
у
W^(Q)
^u(k) 3 (19)
w2(Q). (19)
В силу ограниченности U0 и u(k 1) из (19) получаем
Su(k+1)
w3 ^ 2c1C2Re(M0 + M) W 2(^“)
Su(k)
W^Q) .
Пусть 2c1C2Re(M0 + M) <a< 1, т.е.
a
Re <-
2c1c2(M0 + M) '
(20)
Тогда
3 < M будет вы-
W^Q) J
Su(2) 3 < a Su(1)
W^Q)
W^Q),
(18)
Su(3) 3 < a Su(2) 3 <a2 Su(1)
W^Q) W^Q)
Таким образом, при соответствующем выборе начального приближения u(0) и при выполнении условия (18) решение u(k+1) на каждом шаге итерационного процесса (15), (16) ограничено в норме пространства W32(Q).
Su(k+1) 3 < a Su(k)
W^Q)
w2(Q)
< ... < a
W2(Q)^ •••
Su(1) 3
W^Q) •
Далее
u(k+p) - u(k) 3 ^ u(k+1) - u(k)
W^Q)
w2(Q)
52
РИ, 2007, № 2
+ u(k+2) _ u(k+1) 3 +... + u(k+p) _ u(k+p-1)
W2( Q)
5u(k+1)
3 + Su(k+2) +... + 5u(k+p)
W2( Q) W2(Q)
k+1 +... + ak+p-1) 5u(1) „.3,™ k a -a
W^CQ)
W^(Q)
k+p
W 2(0)
1 -a
(21)
где обозначено у =
5u(1)
3 . Поскольку a < 1, то
W^Q)
u(k+p) - u(k)
3 при k , т.е. последова-
W ^(Q)
тельность {u(k)} является фундаментальной. В силу полноты пространства w2,(0) это означает, что последовательность {u(k)} сходится (с геометрической скоростью), т. е. существует функция
u* є W?(0) П W2(0) такая, что lim u^ = u* .
k^<»
Устремив в неравенстве (21) p к ж, получим оценку для k -го приближения:
ak
w2(Q) _ 1 -аУ . (22)
Проведя сходные рассуждения, можно доказать, что
предельная функция u* удовлетворяет исходному дифференциальному уравнению (13) и является его единственным ограниченным решением при
u *- u(k)
Re <-
1
2c1c2(Mo + M) •
(23)
Объединяя условия (18), (20) и (23), для числа Рейнольдса получаем оценку
Re < min
M
c1c2(M + M0)2
a
2c1c2(M + M0)
(24)
причем константы c1 , c2 зависят только от области Q .
Таким образом, доказали следующую теорему. Теорема. При достаточно малом числе Рейнольдса Re последовательные приближения, формируемые по схеме (15), (16), сходятся к единственному
обобщенному решению u* є w| (Q) П W 2(0) задачи
(13), (14). Условие малости для Re формулируется в виде неравенства (24).
При реализации этого итерационного процесса возникает вопрос, как решать задачу (15), (16)? При известном u(k) задачу (15), (16) можно решить методом Ритца, используя для построения координатной пос-
ледовательности и описания геометрии области метод R-функций. Задача (15), (16) эквивалентна задаче
о
2
нахождения в W 2(0) минимума функционала [19]
I[u] = J [(Au)2 - 2uReJ(A(u0 + u(k)), u0 + u(k))]dxdy.(25) Q
Выбирая систему координатных функций {фі} , функцию u(k+1), доставляющую минимум функционалу (25), будем искать в виде
uN+1) = Е c(k+1) Фі. (26)
І=1
Согласно методу Ритца, задача нахождения коэффициентов c(k+1) i = 1, 2, N, приводит к системе линейных алгебраических уравнений
Ac(k+1) = b(k+1), (27)
где
A = [aij]NxN , b(k+1) = [b(k+1)]Nx1 , C(k+1) = [c(k+1)]Nx1 ,
aij = J Дфі Аф jdxdy, (28)
Q
b(k+1) = Re J J(A(u0 + u(k)), u0 + u(k))9i dxdy =
Q
= Re J J(u0 + u(k), Фi)A(uo + u(k))dxdy (29)
n ■ v J
Отметим, что матрица A системы (27) симметрична, не изменяется от итерации к итерации, вычисляется лишь один раз при решении задачи (11), (12) для функции u0 , а на каждой итерации лишь пересчитывается столбец b(k+1). Таким образом, на (k + 1)-й итерации нужно вычислить лишь N интегралов вида (29).
Из теорем сходимости метода Ритца [19] и теоремы следует, что при числах Рейнольдса, удовлетворяющих соотношению (24), при N ^ жо и k ^ж последовательность u0 + u® сходится в W 2(0) к функции u0 + u =ф , являющейся решением (вообще говоря, обобщенным) задачи (9), (10).
3. Результаты вычислительного эксперимента
Вычислительный эксперимент был проведен для областей трех типов: прямоугольная, имеющая форму полуэллипса, и трапециевидная. Везде предполагалось, что массовые силы f потенциальны, т.е. rotf = 0, а значит, F(x,y) = 0 . В качестве базисных функций выбирались степенные полиномы, тригонометрические полиномы, полиномы Лежандра, кубические сплайны Шенберга. Сплайны показали большую вычислительную устойчивость и им было отдано предпочтение. Шаг сетки сплайнов выбран hx = hy = 0,1. При вычислении интегралов в системах Ритца исполь-
РИ, 2007, № 2
53
зовалась формула Гаусса с 16 узлами по каждой переменной на каждом частичном квадрате 0,1 х 0,1. Предполагается, что стенки области твердые, непроницаемые, неподвижные, кроме верхней, которая движется влево со скоростью 1. Итерационный процесс был построен для чисел Рейнольдса Re = 50 , 100, 200, 300, 400.
Область А. Прямоугольник
Q = {(x,y)|0 <x <a, 0 <y <b} . Математическая модель течения имеет вид
(
Д 2у = Re
3 | . Зш —I Дш— 3x { ду
д( Зш --1 Дш—
Эу^ ТЭх
в Q, (30)
Nx +1 Ny+1
Ф(к+1)(х,у) = £ £ с(к+1)Бз
i=-1 j=-1
Nxx
- i IB
Nyy
Nx+1Ny+1 (k+1)
= S S ci,k+1)9ii(x,y).
i=-1 ІК1 1J
Условие окончания итерационного процесса выбрано
< S
в виде
y(k+1)-y(k)
L2(Q)
Вычислительный эксперимент был проведен при
a = b = 1. Было выбрано Nx = Ny = I0 , є = 10_5. На рис. 1 показаны графики сходимости
y(k+1) -y(k)
L2(Q) :
v>k+1) - vxk)
L2( ^):
5Q
-1, y = b;
0, 3Q \{y = b}.
(31)
vyk+1) - v(yk)
L2(Q):
£ (k+1) (k)
L2(Q)
Итерационный процесс последовательных приближений для задачи (30), (31) строим следующим образом. Задаем начальное приближение y(0)(x,y). Его можно задать произвольно или взять решение, сошедшееся при меньших числах Рейнольдса. Итерационный процесс имеет вид
Д2 y(k+1) = Re
Г_д_Г
dx
Ду1
(k) 9V
(k) ^
3y
д_
dy
Ду'
(k)
Vk)
3x
JJ
в Q ,
(32)
(k+1)
= 0, ^
5Q 3n
(k+1)
-1, y =b; (33)
0, 3Q \{y = b},
k = 0, 1, 2, ....
Структура решения задачи (32), (33) на каждой итерации имеет вид
V <k+1>(x.y) =
L + ro2(x_y)0 M^y,
®1(x,y) + ®2(x,y)
(34)
к нулю при Re = 200 . В таблице приведены координаты “ вихревого центра” (точки, в которой скорости равны нулю: vx = vy = 0), а также соответствующие ему значения функции тока и завихренности в зависимости от числа Рейнольдса.
Re Число итераций xv.c. yv.c. 0,7646 V v.c. C v.c.
0 0,5000 0,1001 3,3374
50 9 14 0,4230 0,7599 0,1009 3,2379
100 0,3789 0,7417 0,1035 2,8224
200 20 0,3945 0,6620 0,1093 2,5493
300 25 0,4206 0,6265 0,1127 2,4887
400 36 0,4351 0,6163 0,1143 1,9096
Как видно, имеет место сходимость в W 2(Q), но с ростом порядка производной скорость сходимости уменьшается. С ростом числа Рейнольдса сходимость итерационного процесса замедляется. В качестве начального приближения выбиралось течение Стокса, а также приближения, полученные при меньших числах Рейнольдса. С ростом числа Рейнольдса для повышения устойчивости итерационного процесса использовалась процедура взвешивания двух приближений.
где ®(k+1)
- неопределенная компонента структуры,
ro(x,y)
^x(a - x) Ла iy(b - y)
_ a _ L b J
®1(x,y) = b -y, ®2(x,y)
^x(a - x) a
ла y ’
С ростом числа Рейнольдса увеличивается зона циркуляционного течения в центре каверны (растет скорость циркуляционного движения от значения 0,205 при Re = 0 до значения 0,349 при Re = 0), растут угловые вихри, а при Re = 300 появляется третий пристеночный вихрь, вихревой центр смещается ко дну каверны, а ymax увеличивается.
ла - знак R-конъюнкции.
Неопределенная компонента ®(k+1) в (34) на каждой итерации аппроксимировалась кубическими сплайнами Шенберга [20, 21]:
Область Б. Область, имеющая форму полуэллипса:
Q = ](x,y)0<y <b, b2fx-|j +y(y-b)2 <^j.
54
РИ, 2007, № 2
а
0.012
0.01
0.008
0.006
0.004
0.002
• •
• •
10
15
20
0.0175
0.015
0.0125
0.01
0.0075
0.005
0.0025
• •
10
15
20
5
б
5
в
0.3
0.25
0.2
0.15
0.1
0.05
• •
• л
10
15
20
Рис. 1. Графики сходимости к нулю при Re = 200 :
(k+4 _ v(k)
^ k+1) o
L2( П)
(а),
L^( П)
(б),
(k+1) (k)
РИ, 2007, № 2
L2( П)
(в),
q( k+4 _c(k)
L^ П)
(г)
5
г
Математическая модель течения имеет вид (30), (31). Итерационный процесс последовательных приближений строится, как и для области А, и имеет вид (32), (33). В случае области Б
Ю(Х,у) = ®1(Х,у) Ла ®2(X,y) , ®1(х,у) = b - у ,
Ю2(х,у) =
2.2 / \2 2 От - b2 (х - О - т'у - b)2
2U2 a b
+ b21 x-+ ^(у _b)2
Вычислительный эксперимент был проведен при a = b = 1 и для Re = 50, 100, 200, 300, 400. Рассмотрим результаты для Re = 100 . В качестве начального приближения было выбрано решение Стокса. Итерационный процесс сошелся с точностью є = 10_5 за 23 итерации. Вихревой центр находится в точке (0,388347, 0,758544), а уmax = 0,0948827 .
Область В. Область, имеющая форму трапеции с вершинами в точках (0,0), (0,b), (a1,0), (a1,b):
Q = <!(x,y)
„ , „ a — a,
0 <y <b, 0 <x <-------1y + ab
a1 <a
Математическая модель течения имеет вид (30), (31). Итерационный процесс последовательных приближений строится, как и для области А, и имеет вид (32), (33). В случае области В
Ю1(х,у) = b - у ,
®2(х,у) = у Ла X лс
b a1b
у------х +—1—
a - a1 a - a1
1 +
a - a1
ю(х,у) = ®1 (х,у) ла ®2(х,у). Вычислительный эксперимент был проведен при 2
a = b = 1, a1 = 3 и для Re = 50, 100, 200, 300, 400.
Рассмотрим результаты для Re = 50, 100. В качестве начального приближения при Re = 50 было выбрано решение Стокса, а при Re = 100 - решение Стокса и решение для Re = 50 . Итерационный процесс для Re = 50 сошелся с точностью є = 10_5 за 7 итераций, а для Re = 100 - за 19 итераций, если в качестве начального приближения выбиралось течение Стокса, и за 8 итераций, если в качестве начального приближения выбиралось решение для Re = 50 . На рис. 2 приведены линии уровня функции тока и завихренности. Графики сходимости
y(k+1) -y(k) L2(Q) ’ vXk+1) - v«
v(k+1) - vyk) £(k+1) (k)
у у L2( П)’
55
к нулю имеют характер, аналогичный показанным на рис. 1. Вихревой центр для Re = 50 находится в точке (0,370542,0,787655), а уmax = 0,0896502, а для Re = 100 - в точке (0,329752, 0,773602), а
V max = 0,090799.
Как показывает анализ результатов, для непрямоугольной области с изменением числа Рейнольдса течение изменяется аналогично: увеличивается зона циркуляционного течения в центре каверны, растут угловые вихри, вихревой центр смещается ко дну каверны, а ymax увеличивается.
Полученные результаты хорошо согласуются с результатами физических экспериментов [1, 2] и результатами, полученными методом фиктивных областей [22, 23] и методом конечных элементов [4, 5].
Выводы
Предлагаемый метод показал свою эффективность на модельной задаче и, по мнению авторов, обладает рядом преимуществ. Матрица системы (27) является симметричной (это следует из (28)) и не изменяется от итерации к итерации, т.е. становится универсальной. Это дает возможность использовать одну и ту же матрицу при различных числах Рейнольдса и на каждой итерации пересчитывать только правую часть системы, что снижает вычислительные затраты по сравнению с другими методами. Кроме того, исходные уравнения не заменяются приближенными, в отличие от метода, предложенного в [8] (линеаризация по методу Ньютона-Канторовича), что повышает надежность метода и позволяет избежать громоздких выкладок. В отличие от сеточных методов решение получается в явном аналитическом виде, что облегчает его дальнейшее использование. Например, зная функцию тока, можно по формулам (4) построить поле скоростей; по формуле (7) найти завихренность; решив соответствующую задачу Неймана для уравнения Пуассона, найти давление. Предложенный в работе подход с некоторыми изменениями может быть распространен на случай многосвязной области и на системы уравнений, например, на уравнения Навье-Стокса в приближении Буссинеска в переменных « функция то ка - темпер атур а».
Научная новизна полученных результатов заключается в том, что впервые разработан итерационный метод решения нелинейного уравнения для функции тока в односвязных областях неклассической геометрии с кусочно-гладкой границей, основанный на совместном применении метода R-функций и метода последовательных приближений, отличающийся от известных методов более простым вычислительным алгоритмом (решение нелинейной задачи сводится к решению последовательности задач с одним и тем же линейным оператором) и универсальностью (алгоритм не изменяется при изменении геометрии области), что позволило получить приближенное решение задачи расчета этого класса течений в областях неклассической геометрии. Получены условия и оценки
0.2 0.4 0.6 0.8 1
Re = 50
0.2 0.4 0.6 0.8 1
Re = 100
Рис. 2. Линии уровня функции тока у и завихренности Q = -Д\р
56
РИ, 2007, № 2
скорости сходимости в норме пространств w| (Q) к
обобщенному решению задачи расчета стационарного течения вязкой несжимаемой жидкости в односвязной области с кусочно-гладкой границей.
Практическая значимость полученных результатов. Разработанные методы расчета плоских течений вязкой жидкости в односвязных областях являются простыми в алгоритмизации и более универсальными, чем используемые в данное время, поскольку при переходе от одной области к другой требуется лишь изменить уравнение границы. Разработанный метод математического моделирования вязких течений позволяет проводить многовариантный вычислительный эксперимент при математическом моделировании различных физико-механических, биологических течений (течения в инжекторах, форсунках, соплах, кровеносных сосудах, обтекание подводных тел, течения газа и нефти в трубах и пр.).
Литература: 1. Ландау Л.Ф., Лифшиц Е.М. Теоретическая физика. В 10 т. Т. VI. Гидродинамика. М.: Физматлит, 2003. 736 с. 2. Лойцянский Л.Г. Механика жидкости и газа. М.: Дрофа, 2003. 840 с. 3. РоучП. Вычислительная гидродинамика. М.: Мир, 1980. 616 с. 4. Donea J., Huerta A. Finite Element Methods for flow Problems. London: Wiley, 2003. 350 p. 5. Zienkiewicz O.C., Taylor R.L. The Finite Element Method. Vol. 3: Fluid Dinamics. Oxford: BH, 2000. 334 p. 6. Рвачев В.Л. Теория R-функций и некоторые ее приложения. К.: Наук. думка, 1982. 552 с. 7. Колодяжный ВМ., Рвачев В.А. Структурное построение полных последовательностей координатных функций вариационного метода решения краевых задач: Препр. АН УССР. Ин-т пробл. машиностр. Харьков, 1975. 75 с. 8. СувороваИ.Г. Компьютерное моделирование осесимметричных течений жидкости в каналах сложной формы // Вестн. НТУ ХПИ. Харьков, 2004. № 31. С. 141 - 148. 9. Колосова С.В. Об обтекании невязкой жидкостью цилиндра в трубе // Прикл. мех., 1971. № 7. В. 10. С. 100 - 105. 10. Максименко-Шейко К. В. Исследование течения вязкой несжимаемой жидкости в скрученных каналах сложного профиля методом R-функций // Проблемы машиностроения, 2001. Т. 4, № 3 -4. С. 108 - 116. 11. Рвачев В.Л., Корсунский А.Л., Шейко Т.И. Метод R-функций в задаче о течении Гартмана // Магнитная гидродинамика. 1982. № 2. С. 64 - 69. 12. Ладыженская О.А. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970. 288 с. 13. Лионс Ж.-Л. Некоторые методы решения нелинейных краевых задач. М.: Мир, 1972. 588 с. 14. Темам Р. Уравне-
ния Навье-Стокса. Теория и численный анализ. М.: Мир, 1981. 408 с. 15. Сидоров М.В. О построении структур решений задачи Стокса // Радиоэлектроника и информатика, №3, 2002. С. 52 - 54. 16. Сидоров М.В. Применение метода R-функций к расчету течения Стокса в квадратной каверне при малом числе Рейнольдса // Радиоэлектроника и информатика. 2002. №4. С. 77 - 78. 17. Колосова С.В., Сидоров М.В. Применение метода R-функций к расчету плоских течений вязкой жидкости // Вісн. ХНУ. Сер. Прикл. матем. і мех. № 602, 2003. С. 61 - 67. 18. Слободецкий Л.Н. Обобщенные пространства С. Л. Соболева и их приложение к краевым задачам для дифференциальных уравнений в частных производных. //Уч. зап. Ленингр. гос. пед. ин-та им. А.И. Герцена. 1958, Т. 197. С. 54 - 112. 19. Михлин С.Г. Вариационные методы в математической физике. М.: Наука, 1970. 512 с. 20. Федотова Е.А. Атомарная и сплайн-аппроксимация решений краевых задач математической физики: Дис. ... канд. физ.-мат. наук: 01.01.07. Харьков, 1985. 170 с. 21. Федотова Е.А. Практические указания по использованию сплайн-аппроксимации в программирующих системах серии “Поле”: Препр. АН УсСр. Ин-т пробл. машиностр.; 202. Харьков, 1984. 60 с. 22. Вабищевич П.Н. Метод фиктивных областей в математической физике. М.: Изд-во Моск. ун-та, 1991. 156 с. 23. Вабищевич П.Н., Вабищевич Т.Н. Численное решение стационарных задач вязкой несжимаемой жидкости на основе метода фиктивных областей // Вычислительная математика и математическое обеспечение ЭВМ. М.: Изд-во Моск. ун-та, 1985. С. 255 - 262.
Поступила в редколлегию 20.05.2007
Рецензент: д-р физ.-мат. наук, проф. Дорошенко В. А.
Тевяшев Андрей Дмитриевич, д-р техн. наук, проф., зав. каф. прикладной математики ХНУРЭ. Научные интересы: математическое моделирование, оптимальное стохастическое управление. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел.: 702-14-36.
Гибкина Надежда Валентиновна, канд. техн. наук, доцент кафедры прикладной математики ХНУРЭ. Научные интересы: экономический риск, актуарная математика, математическая физика, оптимальное управление динамическими объектами. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел.: 702-14-36.
Сидоров Максим Викторович, ассистент каф. прикладной математики ХНУРЭ. Научные интересы: математическое моделирование, математическая физика, теория R-функций и ее приложения. Увлечения и хобби: история культуры. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел.: 702-14-36.
УДК658.012
ТЕКУЩИЙ РЕГРЕССИОННЫЙ АНАЛИЗ С УЛЬТРАКОРОТКИМ СКОЛЬЗЯЩИМ ОКНОМ
ЯКУНИН А.В.__________________________
т
Предлагается модификация UDU -факторизованной схемы рекуррентного алгоритма ТРА для задачи параметрической идентификации линейной регрессионной модели по матрице неполного ранга.
1. Введение
Пусть идентифицируемый объект представлен линейно параметризованной регрессионной моделью
T с
Уп =c xn +^n ,
где xn eRN - вектор входных сигналов; yn eR1 -г. N
выходной сигнал; c eR - вектор оцениваемых параметров; £,n eR1 - помеха измерения выходного сигнала; n = 0,1,2, ... - дискретное время. Если параметры объекта изменяются во времени, то для отслеживания их динамики применяются алгоритмы теку-
РИ, 2007, № 2
57