CC BY
39
(Wap (t,s)) =
Ёф(1) (t)9(1) (s) Ёф(1) (t)9(2) (s) j=i j=i
E9(2)(t)9((1)(s) E9((2)(t)9(2)(s)
U=1 j=1
так как
9(a)(t) = (eltAZca, gj) = \ f (j)(u,t)ya (u)du,
то из Uj (x,t) = U^)(0,t) - iJf(j)(u,t)ya (u)du при
0
U^)(0,t) = 0 получаем, полагая x = I,
9(a)(t) = iU^j)(^, t).
Выводы
С помощью треугольных моделей диссипативных операторов построены представления для корреляционных матриц и векторных нестационарных случайных процессов, которые определяются только спектром (лежащим в верхней полуплоскости в случае дискретного спектра или бесконечнократного спектра в нуле).
Научная новизна. Впервые дан спектральный анализ корреляционных матриц и самих векторных нестационарных случайных процессов при помощи треугольных моделей диссипативных опер аторов и линейных систем, ассоциированных с операторными узлами или комплексами.
Перспективы исследования. Для получения спектральных представлений корреляционной матрицы ВНСП можно воспользоваться спектральным разложением адекватных векторных кривых x(t), используя спектральную теорию систем коммутирующих несамосопряженных операторов, развитую в [8]. Однако соответствующие спектральные представления носят очень громоздкий характер, поэтому мы ограничились тем классом нестационарных векторных кривых, для которых все Aj совпадают (аналог стационарных и стационарно связанных случайных процессов, для которых в соответствующем гильбертовом пространстве Aj = A = A* (j = 1,2)).
Практическая значимость работы заключается в том, что, воспользовавшись универсальными моделями диссипативных операторов [4], можно рассмотреть векторные эволюционно представимые случайные процессы произвольного конечного ранга нестационарности и получить модельные представления для корреляционных матриц и соответствующих векторных случайных процессов.
Литература: 1. Розанов Ю. А. Спектральная теория многомерных стационарных случайных процессов с дискретным временем // Успехи матем. наук. 2 (13). М., 1958. С. 92-142. 2. Розанов Ю. А. Стационарные случайные процессы. М.: Физматгиз, 1963. 284 с. 3.ХеннанЭ. Многомерные временные ряды. М.: Мир, 1974. 576 с. 4. Livshits M. S. and Yantsevitch A. A. Operator colligations in Hilbert space, Wiley, New-York, 1979. 210 p. 5. Петрова А. Ю. Корреляционная теория некоторых классов нестационарных случайных функций конечного ранга нестационарности // Радиоэлектроника и информатика. 2007. № 1. С. 29-34. 6. Когут Е. А., Черемская Н. В., Янцевич А. А. О представлении резольвент вольтерровых операторов // Крайові задачі для диференціальних рівнянь: Зб. наук. праць. Київ, Ін-т математики НАН України, 1998. Вип. 1(17). С. 99-101. 7. Кирчев К. П. Об одном классе нестационарных случайных процессов // Сб. Теория функции, функциональный анализ и их приложения. Х., 1971. Вып. 14. С. 150-159. 8. Золотарев В. А. Аналитические методы спектральных представлений несамосопряженных и неунитарных операторов. Х.: Изд-во ХНУ, 2003. 342 с.
Поступила в редколлегию 14.12.2007
Рецензент: д-р физ.-мат. наук, проф. Золотарев В. А.
Янцевич Артем Артемович, д-р физ.-мат. наук, проф., зав. каф. высшей математики и информатики механикоматематического фак-та ХНУ им. В. Н. Каразина. Научные интересы: спектральная теория несамосопряженных или неунитарных операторов; линейные системы, асссо-цированные с операторными узлами; прикладная теория случайных процессов; прогнозы фильтрации случайных процессов. Адрес: Украина, 61077, Харьков, пл. Свободы, 4, тел. 707-55-42.
Петрова Анжела Юрьевна, соискатель, преподаватель кафедры информационных технологий и математики Харьковского гуманитарного университета «Народная украинская академия». Научные интересы: моделирование случайных процессов, нестационарные случайные функции. Адрес: Украина, 61000, Харьков, ул. Лермонтовская, 27, тел. 716-44-09 (доп. 2-22).
УДК 517.95 : 519.63
ОБ ОДНОМ ПОДХОДЕ К МАТЕМАТИЧЕСКОМУ МОДЕЛИРОВАНИЮ ПЛОСКИХ СТАЦИОНАРНЫХ КОНВЕКТИВНЫХ ТЕЧЕНИЙ В КОНЕЧНЫХ ОДНОСВЯЗНЫХ ОБЛАСТЯХ
АРТЮХ А.В., ГИБКИНА Н.В., СИДОРОВМ.В.
Рассматривается применение метода R-функций в сочетании с методом последовательных приближений для решения задачи расчета плоских стационарных течений
40
вязкой несжимаемой теплопроводной жидкости в конечных односвязных областях. Доказывается сходимость построенного итерационного процесса в норме пространства W23(Q) х W1(Q) к обобщенному решению исходной задачи. Получены оценки скорости сходимости. Предложенный метод протестирован на модельных областях, полученные приближенные решения сравнены с решениями, полученными другими авторами.
Введение
Актуальность задачи. Изучение законов движения жидкости играет важную роль в развитии техники и естествознания. Исследования в этой области стимулируются потребностями авиации, кораблестроения, теплоэнергетики, геофизики, биологии и т.д. за последние десятилетия сфера исследования и применения
РИ, 2007, № 4
явлений, связанных с движением жидкости, постоянно расширяется и охватывает ведущие направления промышленности (химические технологии, нефте- и газоразработка, металлургия и т. д.) и ряд естественных наук (биология, физика атмосферы и океана и др.). Во многих практически важных случаях жидкость можно с большой достоверностью считать вязкой несжимаемой ньютоновской средой, и проходящие в ней процессы могут быть промоделированы с помощью уравнений Навье-Стокса [1, 2]. Различные задачи, возникающие при изучении динамики вязкой жидкости, могут быть исследованы теоретическим путем или с помощью физического эксперимента. Однако с развитием ЭВМ все активнее используется математическое моделирование. Существует множество численных методов, применяемых при расчете вязких течений. Литература по этому направлению обширна [3-5 и др.]. В основном эти численные методы используют метод конечных разностей и метод конечных элементов. Они просты в реализации, но не обладают необходимым свойством универсальности - при переходе к новой области (особенно неклассической геометрии) необходимо генерировать новую сетку, а часто и заменять сложные участки границы простыми, составленными, например, из отрезков прямых. Точно учесть геометрию области можно, воспользовавшись конструктивным аппаратом теории R-функций, разрабатываемой акад. В.Л. Рваче-вым и его учениками [6,7 и др.]. Задачи гидродинамики решались в работах С.В. Колосовой, К.В. Макси-менко-Шейко, И.Г. Суворовой, Т.И. Шейко и др. [811], однако в основном рассматривались задачи динамики идеальной жидкости или вязкой для случаев, когда можно построить решение за счет удачного выбора ко ор динат (осесимметр ические течения, тече -ния, обладающие винтовой симметрией, и т. п.), теплопроводность при этом не учитывалась. Поэтому разработка новых, а также совершенствование существующих методов математического моделирования динамики вязкой теплопроводной жидкости на основе метода R-функций является актуальной научной проблемой.
Цели и задачи исследования. Целью данной работы является создание современного и эффективного метода математического моделирования плоских стационарных течений вязкой теплопроводной несжимаемой жидкости в конечных областях неклассической геометрии с кусочно-гладкой границей. Основные результаты по теоретическому обоснованию корректности начально-краевых и краевых задач для уравнения Навье-Стокса получены О. А. Ладыженской [12], Ж.-Л. Лионсом [13], Р. Темамом [14]. Для достижения поставленной цели решаются следующие задачи: разработка и обоснование метода расчета плоских стационарных течений теплопроводной вязкой жидкости в односвязных областях неклассической геометрической формы с кусочно-гладкой границей в переменных « функция тока - температура» (приближение Буссинекса); применение разработанных численных методов для решения модельных задач расче-
РИ, 2007, № 4
та течения вязкой несжимаемой жидкости в односвязных областях при различных числах Рейнольдса и чисел Грасгофа и Пекле.
В этой работе методы, предложенные в [24], распространяются на систему нелинейных уравнений.
1. Постановка задачи
Рассмотрим плоскую стационарную задачу о конвективном движении вязкой несжимаемой жидкости в конечной односвязной области Q с кусочно-гладкой границей 3Q. Такое течение можно описать с помощью кр аевой задачи для системы дифференциальных уравнений относительно функции тока у (x,y) и температуры T(x,y) вида (приближение Буссинеска):
2 ЭТ
Д у = ReJ(Ay, у) + Gr---ьReF(x,y) в Q ,
3x
AT = PeJ(T, у) + PeG(x,y) в Q ,
= g(s), s edQ,
4»= f <s).
on
Tl3Q = ^(s), s ecQ..
(1)
(2)
(3)
(4)
„ T.. . ЗДу Зу ЗДу 3y
Здесь J(Ay, y) =----------------, Re - число Рей-
dx dy dy dx
нольдса; Gr - число Грасгофа; Pe - число Пекле; G(x, y) - объемная плотность тепловых источников;
df
li(s) - распределение температуры на границе; —, g
ds
- некоторые распределения нормальной и касательной составляющих скорости потока соответственно; n - внешняя нормаль к q . Предполагаем, что
F, G є L2(Q), і' є W2(3Q), g, h є L2(3Q), система координат выбрана так, что сила тяжести сонаправлена с вектором (0, -1).
2. Выбор и обоснование метода решения
Для решения задачи (1) - (4) применим процесс последовательных приближений.
Обозначим через (u0, v0) решение следующей задачи:
Д2u0 = ReF(x, y) + Gr -^V° в Q , 3x
Av0 = PeG(x,y) в Q ,
I f, \ ^0
u"l - =f <s)- Ш
= g(s), s є 3Q,
V0 |aQ = h(s), s edn.
(5)
(6)
(7)
(8)
В сделанных предположениях задача (5) - (8) однозначно разрешима и (u0, v0) є W4(Q) x W2(G). Ее решение может быть построено с помощью метода, базирующегося на методе R-функций [15 - 17].
В задаче (1) - (4) сделаем замену:
41
V = Uo + U,
T = Vo + V,
(9)
(10)
+ U), U0 + U)+Gr & в n ■ (11)
v0 + v U0 + U) в Q , (12)
9u = 0 ? SQ
>= 0> 9n (13)
v| яо loQ = 0. (14)
Д2и(к+1) = ReJ(A(u0 + U(k)), u0 + U(k)) +
G
+Gr~ В Q.
ox
Av(k+1) = PeJ(v0 + v(k), u0 + U(k)) в Q , (16)
3u(k+1)
,(k+Д =
\SQ
dn
r(k+^ — 0
= 0
||Uj| < m
II IIw32(Q) =
M
IIw2( q)
< L
j = 0, 1, 2, ..., k . Из (15), (16) имеем
IM 3 <
II llw|( Q)
+ Gr
IIL^ Q)
5v(k)
dx
) <
4( Q)
^(Rejjj^д(u0 + u0 +
< CjC2Re(| ujl^^+1 Iu^II )2 + CjGrl Iv^ll <
1 2 V|1 0|™дп) У IIwI(q)' 1 II IIw2(q)
< CjC2Re(M0 + M)2 + CjGrL ; (19)
< C2C3Pe v0 + V
(k)||
,(k)||
где u, v - новые неизвестные функции. Это приводит к задаче
dv
IIw2( q) II II w к q)
< C^Pe^ + L)(M0 + M), (20)
где 1НЦQ) < M0 , Ilv0Ilw2(Q) - L0 .
Из (19), (20) следует, что итерационный процесс (15) - (18) будет давать ограниченное решение, т.е. будет выполнено
,(k+1)11
||u'“ II < M |Mk+9II < L
II llw|(q) ’ II IIw2(q) :
если
Re <
Мр
Для решения задачи (11) - (14) воспользуемся методом последовательных приближений. Пусть начальное приближение (U^, v^) задано. Тогда при известном значении (Uk^, v^) функций (u, v) на k-й
итерации приближение на (k +1) -й итерации находится как решение задачи
Gr <
C^CM + М)2 ’ М(1 -р)
Pe <-
CjL
L
(15)
(17)
(18)
Изучим вопрос сходимости итерационного процесса (15) - (18). На каждом шаге итерационного процесса нужно найти пару функций
(U(k+1), v(k+1)) є (W32(П)П W2(Q)}x W2(О). Пусть
(21)
(22) (23)
C2C3(L0 + L)(M0 + M) ’ (24)
где 0 < p < 1.
Докажем сходимость последовательности (U^, v^), k = 0, 1, 2, .... Рассмотрим разности
(5U(k+1),5v(k+1)) = (U(k+1) -U(k),v(k+1) -v(k)).
Они удо влетворят ур авнениям
ДWk+1) = Re[J(A(U0 + U(k)), U0 + uw) -
-J(A(U0 + U(k-1)), U0 + U(k~1})] + Grв q ,
ox
ASv(k+1) = Pe[J(v0 + v(k), U0 + uw) -
-J(v0 + v(kЛ U0 + Uk^9)] в Q и краевым условиям
c6Uk+1)
для всех
5Uk+1) = 0
laa ’
k = 0, 1, 2, .... Тогда
9n
= 0 sv(k+1) = 0
||5U(kH
1)
HWK n)
^ C,
3Re||J(A(u0 + Uk 9), 5u^) 5Sv(k)
+J(A§U^, U0 + mL) + C1Gr
< 2C!C2Re(M0 + M)||5U^|w^ + C1Gr||5^||w^^ , (25)
Isv^l <C3Pe|J(v0 + vk-9, sU^) +
II IIw2(q) 3 II v 0 ’ 2
3x
vk)|
q)
vk+1)ll
IIw2( q)
3Pe||j(v„ + v(k), U0 + U(k)|
+J(5v^, u0 + Uk 9
l2( q)
4( n)
42
РИ, 2007, № 4
< СЛРЄ(||v. + v(k-'ЧЦJK1,(0|
llsv^ II ||u0 + Uk ^ II ) <
II llw‘2( Q) II 0 llw|( Q)7
pe((Lo + L)| |б,,<‘>Ц o) +
+(Mo + ЧИЦ). (26)
< C2C3
Обозначим
£(k) = (|<k), |2k))T = (15uW|
A =
, ||Sv(k)|| )T
w|(Q) II llw2(Q) ’
a11 a12 2CjC2Re(M0 + M) CjGr "
_a21 a22 _ _ C2C3Pe(L0 +L) C2C3Pe( M0 + M)
< max}
= max{
Как видно из последнего неравенства, если max{au + a12, a21 + a22} = a< 1, (28)
то
1^ ^ 11 = max
Re <-
ap'
Gr <
2C'C2(M0 + M) ’ a(1 -p')
Pe <-
(29)
(30)
(31)
(W32(Q)П W2(^))X W2(Q), причем имеют место оцен-
max<||u* - uk)H
IIw|(q) I где обозначено У = max
iv*-v(k)|| 1 . \<—l, (32)
W2( ^)J 1 -a
w|(Q) II IIw
Q)l '
Докажем, что предельные функции (u*, v*) удовлетворяют дифференциальным уравнениям (11), (12). Рассмотрим уравнение (11). Удовлетворение краевым условиям (13) очевидно. Пусть u* є W2(Q). С учетом (15), (31), (32) имеем
(по смыслу задачи все входящие в вектор ^^ и матрицу A величины неотрицательны).
Тогда (25), (26) в матричном виде запишется так:
£(k+1) < A£(k). (27)
Будем рассматривать (27) как систему линейных неравенств с нормой ||£|| = max{|2}. Тогда
||^k+1)|| = max j^1k+1), |(k+1>j <
< ma^|au^14 + a12|(2k), a21^14 + a22^2^} ^
Ф11 + a12. a21 + a22}max{ ^ |(24} =
^{a11 ^ a12, a21 ^ a2^||^^| |.
Д2u* - ReJ(A(u0 + u*), u0 + u*) - Gr
dv*
dx
W2‘( Q)
<||д2(u* -u(k+1))||W_ ^ + Re||J(A(u0 + u*), u0 + u*) -
-J(A(u0 + u^), u0 + u^
W^( Q)
+ Gr
T-(v*- v(k))
ox
< u*- uk+1>IL +
wY( Q)
llW32( Q)
+ 2ReCc2 ||u0 + u(k)|| llu* -u(k)|| + Grllv* -v(k)|| <
2II 0 IIw32(q) II IIw32(q) II IIw3(q)
a
Wl(Q) , iH^lw,^6удЄт схоДИТься к нулю. Ограничения (28) на элементы матрицы A приводят к таким ограничениям для чисел Рейнольдса Re , Грасгофа Gr и Пекле Pe :
< (^ + 2ReCC2(M0 + M) + Gr)-----у—у 0 при
1 -a
k .
Рассмотрим уравнение (12). Удовлетворение краевому условию (14) очевидно. Пусть v* є W2(Q). С учетом (31), (32), (16) имеем
К- PeJ(v0 + • “0 + »■ 1,,,
SK'- '">|Ц ч +
+Pe||J(v0 + v*, u0 + u*) - J(v0 + vW, u + uW)||w£
< llv* - vk+^ + C Pe(Jv0 + v*, u* - u^)|| +
11 iiw‘2(^ 'll v 0 ;Hl2(Q)
+||j(v*-v(k), u0 + u*)|| ) <||v*-v(k+1)|| +
II V 0 /|Il2(Q) II IIw КQl
C2C3(L0 + M0 + L + M)
Приходим к выводу о сходимости последовательности (Vy v^) к некоторым функциям (u*, v*) из
+Cc2 Pe(||v0 + vil w l u*- u^l + 2 N1 0 llw‘2( Q)|| Ilw32( Q)
+||v* - v(k)|| 1 ||u0 + u*|| ) <
II IIw 2( Q)ll 0 HWK o)'
ки
C
РИ, 2007, № 4
43
< (а + Pecc2 (L0 + M0 + L + M))-у —^ 0
1 -a
при k
Выясним теперь, при каких ограничениях на числа Рейнольдса Re, ГрасгофаGr и Пекле Pe пара(u*, v*) будет единственным решением дифференциальных уравнений (11), (12). Пусть задача (11) - (14) имеет
еще одно решение (u**, v**), удовлетворяющее условиям (21). Рассмотрим разности su = u* -u** и 5v = v* - v**. Они удовлетворяют уравнениям
A25u = Re [j(A(u0 + u*), u0 + u*) -
-J(A(uo + u*), u0 + u**)]+ Gr в q ,
J ox
A5v =
= Pe[J(v0 + v*, u0 + u*)-J(v0 + v**, u0 + u**)] в
Q .
Получим:
IMIw,^ C1(Re||J(A(u0 + u*), u0 + u*)-
38v
) <
L((
) +
2( a) _ II Hw2( q)
-J(A(u0 + u**), u0 + ^1(0) + Gr
< Cj (Re ||J(A(u0 + u**), 5u) -
-J(A5u’ u0 + u*GrlNIw((q)) *
< C1C2Re||НЦa) (u0 + u*1L^) +Iu0 + u1|w((Q) +C1GrlI54Iw2(q) ^
< 2C,C2Re(M0 + M)||Hlw^Q) + C1GrNlwi(Q) , l8vlw((q) ~ ll8vllw((q) ~
< C3Pe||J(v0 + v*, u0 + u*) - J(v0 + v**, u0 + u**^(q)
= C3Pe ||J(v0 + v**, Su) + J(Sv, u0 + u**)|^м <
< C2C3Pe(||v0 + v*1|w.(q)INIw((a) +
+ INIw((Q)|Iu0 + u*1|wl(Q)) *
< C2C3Pe(L0 + L)|| 5u|
W2( Q)+ C2C3Pe(M0 + M)|| S^lwi
W2( Q)
Таким образом,
IMIw,Q) ^ 2C1C2Re(M0 + M)||8u||w2(Q) + ^гНЦ(
|МЦ(Q)^ C2C3Pe(L0 + L)||S^lwKQ) +
+C2C3Pe(M0 + M^|НЦQ). (34)
Перепишем систему неравенств (33), (34) в виде
(1 - 2C!C2Re(M0 + m))||8HIw^ * C!Gr| |8v||w2(Q),
(1 - C2C3Pe(M0 + M))||8v||w((Q) ^ C2C3Pe(L0 + L^18u
llw|(Q)-
Из (28), (30) следует, что 1 - 2CjC2Re(M0 + M) > 0 и 1 -C2C3Pe(M0 + M) > 0. Тогда
ll8ul
IN
< c,Gr ,,
w32<q) “ 1 -2CjC2R^M0 + M v W‘2(n), (35)
1 - C2C3Pe( M0 + M).
Iw((a) >
-|NIw2(q) . (36)
C2C3Pe( L0 + L)
Система (35), (36) будет иметь лишь нулевое решение
НИЦ ^=0, НИЦ ^=0, если
CjGr 1 - C2C3Pe(M0 + M)
1 - 2CjC(Re(M0 + M) C2C3Pe( L0 + L) (37)
Объединив (22) -(24) и (28) - (30), получим следующие ограничения для чисел Рейнольдса Re , Грасго-фа Gr и Пекле Pe :
Re < min
Mp
apj
C1C2(M0 + M)2’ 2CjC( (M0 + M) (38)
Gr <minJMM Ml
CjL Cj
(39)
Pe < min
a
C(C3^0 + L)(M0 + M) C(C3(L0 + M0 + L + M)
(40)
44
где p и p1 - произвольные числа из (0, 1), а константы C1, C2 и C3 зависят только от области Q .
Полученные результаты можно сформулировать следующим образом.
Теорема. При достаточно малых числах Рейнольдса Re, Грасгофа Gr и Пекле Pe, удовлетворяющих условию (3 7), последовательные приближения, формируемые по схеме (15) - (18), сходятся к единственному обобщенному решению задачи (11) -
(13), принадлежащему (W32(Q)ПW((Q))х W((Q)•
РИ, 2007, № 4
k
Условие малости для чисел Re, Gr и Pe формулируется в виде неравенства (38) - (40).
При реализации вычислительного процесса по схеме (15) - (18) функции Uk+1) и Vk+1) также могут быть найдены с помощью метода Ритца.
Как видно из (15) - (18), эта задача при известных функциях u0, v0,uW иv(k) распадается на две независимые задачи (15), (17) и (16), (18).
С задачей (15), (17) свяжем оператор А, этой краевой задачи, действующий по правилу
AjU = Д2u ,
с областью определения
D ai)
C4(Q) n C2(Q),u| an
9u
9n
= 0
. (41)
Пополнив множество (41) в норме, порожденной скалярным произведением
[uj, u^ = J AUjДи2 dx dy
Q ’
получим энергетическое пространство HI.
С задачей (16), (18) свяжем оператор Ап этой краевой задачи, действующий по правилу
Anv = -Av,
с областью определения
D(An) = {v|v є C2(Q) ПC1(Q),v|aQ = 0} . (42)
Пополнив множество (42) в норме, порожденной скалярным произведением
К vJn =jVv1Vv2dxdy,
получим энергетическое пространство Hn. Тогда при известных функциях u0 , v0, u^ и v^ задача (15),
(17) эквивалентна задаче нахождения в HI минимума функционала
Iju] = j[(Au)2 -
Q
^(k)
-2u(ReJ(A(u0 + u(k)), u0 + u(k)) + Gr——)]dxdy, (43)
dx
а задача (16), (18) эквивалентна задаче нахождения минимума в HII функционала [19]
I2 [v] = j[(Vv)2 + 2vPeJ(v0 + v(k), u0 + u(k))] dxdy (44)
Q
Выберем систему координатных элементов {ф;} такую, что
1) ф; є HI для всех ;;
2) система {ф;} полна в HI;
3) для любого N система ф1, ф2, ...,фн линейно независима;
и систему координатных элементов {тJ, такую, что
1) т; є HII для всех ;;
2) система {xj полна в HII;
3) для любого N система т1, т2, xN линейнонезависима.
Функцию uk, реализующую минимум функционала (43), будем искать в виде
а функцию vk, реализующую минимум функционала (44), будем искать в виде
£+1)
N
(46)
Согласно методу Ритца [19] задача нахождения коэффициентов ер+^, ; = 1, 2, N1, из (45) и коэффици-
ентов dp+1), ; = 1, 2, ..., n2 , из (46) приводит к решению двух систем линейных алгебраических уравнений:
bk+1 WI ? (47)
bk+1) UII ’ (48)
где
N,.N,,bIk-[ Ы'-] N„1. C ‘ =[Є;“ ] N.„. AJI =[a" ] N. ■ N. ' bIk =[ b"(k ] N..-
d(k+1) = [dp+1) ] N 1,
a' = J Дф;Дф] dxdy , (49)
bj(k+1) = Re j J(A(u0 + u(k)), u0 + u(k)^; dxdy +
,dv(k)
+Gr I------ф; dxdy
aj' = j ^T;^T j dxdy
(50)
(51)
bj1 (k+1) = -Pej J(v0 + v(k), u0 + u(k))x; dxdy (52)
Q
Как видно, матрицы систем (47), (48) являются симметричными, не изменяются от итерации к итерации, вычисляются один раз при решении задачи (5) - (8), а на каждом шаге итерационного процесса лишь пересчитываются векторы правых частей, что приво-
РИ, 2007, № 4
45
дит к вычислению Nj + N2 интегралов вида (50), (52). Используя свойства якобианов, интегралы (50) и (52) можно преобразовать к виду
bj(k+J) = ReJJ(u0 + u(k), Фі)Д(и0 + u(k))dxdy -
Q
Gr f^-yMdxdy
Idx =
b“ (k+1) = -Pej J(u0 + u(k), ті)(у0 + v(k))dxdy
Q
который уже не содержит под интегралом производных третьего порядка.
Из теорем о сходимости метода Ритца [19] и теоремы в данной статье следует, что при числах Рейнольдса Re, Грасгофа Gr и Пекле Pe, удовлетворяющих соотношениям (37) - (40), при N1 ^^,N2 ^да и k ^ да последовательность (u0 + uN^J, v0 + vNk2) сходится к функциям (у*, T*) = (u0 + u*, v0 + у*), являющимся решением (вообще говоря, обобщенным) задачи (1) - (4).
Рассмотрим задачу (5) - (8). Ее можно решить следующим образом. Пусть известны функции f є W 2(Q), g є W2(Q), h є W 2(Q), такие, что f |aQ = f , g|an = g, h|aQ = li. Тогда в соответствии со структурным методом В. Л. Рвачева [6] структуру решения задачи (5) -(8) можно записать в виде
ц, = f-ra(Djf + g) + ю2®0, (53)
У0 = h + raY0, (54)
где Ф0, Y0 - неопределенные компоненты,
D _д(й д дю д
Dj ~~дх ~дх +~dy 5у , а функция ra(x,y) удовлетворяет таким условиям:
1) ra(x,y) є W4j(Q);
2) ro(x,y) = 0 при всех (x,y) є дО;
3) ro(x,y) >0 при всех (x,y) еО;
4) |V ю| = 1 при всех (x,y) є дО.
Для аппроксимации неопределенных компонент в (53), (54) можно воспользоваться методом Ритца [19],
представив Ф0 и Y0 в виде линейной комбинации элементов какой-либо координатной системы (например, полиномов, сплайнов и пр.).
Результаты вычислительного эксперимента
Вычислительный эксперимент был проведен для прямоугольной области. Предполагалось, что массовые силы f потенциальны, т.е. rotf = 0 и внутренние источники тепла отсутствуют, а значит, F(x, y) = 0, G(x, y) = 0. В качестве базисных функций выбирались степенные полиномы, тригонометрические полиномы, полиномы Лежандра, кубические сплайны
Шенберга. Сплайны показали большую вычислительную устойчивость и им было отдано предпочтение. Шаг сетки сплайнов выбран hx = hy = 0,1. При вычислении интегралов в системах Ритца использовалась формула Гаусса с 16 узлами по каждой переменной на каждом частичном квадрате 0,1 х 0,1. Предполагалось, что жидкость ограничена твердыми неподвижными стенками. Итерационный процесс был построен для чисел Рейнольдса Re = 1, Пекле Pe = 1 и Грасгофа Gr = 50, 150 . Условие окончания итерационного процесса выбрано в виде
max<||у(k+1) -у(k)H
||Tk+1 — Tk) її
- є, є = 10 5
М q) II IM q) j
Рассмотрим прямоугольную область О = {(x,y)|0 < x <a, 0 < y <b}. Математическая модель свободноконвективного течения имеет вид
д V =
( я ( я,..л я Ґ
= Re
dx
д ( Зш Ду
dy) dy
д ( Зш Ду
3x
„ сТ
+Gr аГ в п, (55)
(56)
ДТ = Pe
т|_Af т
3x ^ 3y J dy ^ 3x
в О-
I _
^ 3n
= 0
т =f1H2x -1. у =0;
sn [ 0, (x,y) є ЗО \{y = 0}.
(57)
(58)
Итерационный процесс последовательных приближений для задачи (55), (58) строим следующим образом. Задаем начальное приближение - пару функций ^0 (x,y), (x,y). Его можно задать произвольно
или, например, взять решение задачи
II о x з дт(0) = 0 в Q , (59)
V1 = 0 loQ ’ 3n = 0 (60)
Ї1-12x -1, 0, (x,y) Є aQ y = 0; dO \{y = 0}. (61)
Т0) I =
Для решения задачи (59) - (61) воспользуемся методом R-функций. Структура решения имеет вид
y^(x,y) = ro2(x,y)®^^(x,y),
Т0) (x, y) = Т) (x, y) + ffl(x, y) Y<0) (x, y),
(62)
(63)
где
ra(x,y) =
x(a - x)
A
y(b - y)
46
РИ, 2007, № 4
To(x,y) =
2xrot (x, y)ro3 (x, y) + (2 - 2x)k>j (x, y)ro2 (x, y) ra2 (x, y)ro3 (x, y) + Qj(x, y)ra3 (x, y) + (x, y)ro2 (x, y)
ro^x,^ = x(1- x)(1 - y),
(64)
ff>2(x,y) =
юз (x,y) =
_L_fх _1, _y-
16 l, 4 '
ЇЇТГ41 ^
Ла(-y )
Ла (-y2). (65)
Итерационный процесс имеет вид
Д 2^k+4 = Re
(я ( л и л
а V к
\dx К
dy
я ( я W ^
а V ч
dy
dx
JJ
5^(k)
+G^^— в Q.
ox
(66)
дТк+1) _
Pe 5 ( tM 9y^k) ^ d ( Tk) 9y^ ^
dx 'v 9y j dy 'v dx j
(k+1) I _
(k+1)
9n
= 0
в О, (67)
(68)
Nx +1Ny +1
= Z Z 91J(x,y).
1=-1 j=-1
Вычислительный эксперимент был проведен при a = b = 1. Результаты в виде линий уровня температуры (изотерм) и линий уровня функции тока (линий тока) приведены на рис. 1,2 (Re = 1, Pe = 1и Gr = 50).
Сходимость с точностью є = 105 была достигнута за 7 итераций.
Рис. 1. Линии уровня функции тока при Re = 1, Pe = 1 и Gr = 50 (x,y) = 0,0064)
T“>| ={i-12x -i|, y = 0;
I» [ 0, (x,y) Etliljy = 0); (69)
Структура решения задачи (66) - (69) на каждой итерации имеет вид
^k+1(x,y) = ro2(x,y)^k+4(x,y), (70)
Tk+1) (x, y) = T0 (x, y) + ra(x, y) Y<k+1) (x, y), (71)
где kk+4 - неопределенные компоненты структуры, а функции ro(x,y), T0 (x,y) определяются с помощью соотношений (64), (65).
Неопределенные компоненты фМ , ,
k = 0, 1, 2, ..., в (62), (63) и в (70), (71) аппроксимировались кубическими сплайнами Шенберга [22, 23]:
Nx +1 Nv +1
1>,k)(x.y) = s У 4ЧВз
l=-1 j=-1
Nx
--1 IB.
( N,
y
Nx +1Ny +1
= Z Z cik)ф^у^
i=-1 j=-1
Nx +1N, +1
^ k)(x,y) =Z Z 4^3
i=-1 j=-1
N x
-1 IB.
(Nyy Л
,---J
V
Рис.2. Изотермы при Re = 1, Pe = 1 и Gr = 50 . 1
0.8
0.6
0.4
0.2
0
P „ м > t 1Г * ■Т ч
* " 4
_JL і і t V і \ 1 * ',Л ч
ч і V"
і * і ч
і f і і ■4- L t ' “Т / 1 TJ і...: і ! .4- І \ ъ
+ А- V . ' * ...і j t
^t _Ч, ■V .. S f \ V * а/_
/
0 0.2 0,4 0.G 0.8 1
Рис.3. Векторное поле скоростей при Re = 1, Pe = 1 и Gr = 50
РИ, 2007, № 4
47
Рис. 4. Линии равной завихренности
Полученные результаты хорошо согласуются с результатами физических экспериментов [1, 2] и результатами, полученными методом фиктивных областей [22, 23]. Расхождения составили около 5%.
Выводы
Предлагаемый метод показал свою эффективность на модельной задаче и, по мнению авторов, имеет ряд преимуществ. Матрицы систем (47), (48) являются симметричными и не изменяются от итер ации к итер а-ции, т.е. становятся универсальными. Это даёт возможность использовать одну и ту же матрицу при различных числах Рейнольдса, Грасгофа, Пекле и на каждой итерации пересчитывать только правые части системы, что снижает вычислительные затраты по сравнению с другими методами. В отличие от сеточных методов решение получается в явном аналитическом виде, что облегчает его дальнейшее использование для определения различных характеристик течения.
Научная новизна полученных результатов заключается в том, что впервые разработан итерационный метод решения систем нелинейных уравнений для функции тока и температуры в односвязных областях неклассической геометрии с кусочно-гладкой границей, основанный на совместном применении метода R-функции и метода последовательных приближений, отличающийся от известных методов более простым вычислительным алгоритмом (решение нелинейной задачи сводится к решению последовательности задач с одним и тем же линейным оператором) и универсальностью (алгоритм не изменяется при изменении геометрии области), что позволило получить приближенное решение задачи расчета этого класса течений в областях неклассической геометрии. Получены условия и оценки скорости сходимости в норме
W23 (Q) х W2 (Ц к обобщенному решению задачи расчета стационарного течения вязкой теплопроводной несжимаемой жидкости в односвязной области с кусочно-гладкой границей.
Практическая значимость полученных результатов. Разработанные методы расчета плоских течений вязкой жидкости в односвязных областях являются
простыми в алгоритмизации и более универсальными, чем используемые в данное время, поскольку при переходе от одной области к другой требуется лишь изменить уравнение границы. Полученные результаты позволяют проводить вычислительные эксперименты во время математического моделирования различных физико-механических, биологических течений (течения в инжекторах, форсунках, соплах, обтекание подводных тел и т. д.).
Литература: 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. Ладыженская О.А. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Мир, 1972. 588 с. 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. Вабищевич П.Н, Вабищевич Т.Н. Численное решение стационарных задач вязкой несжимаемой жидкости на основе метода фиктивных областей // Вычислительная математика и математическое обеспечение ЭВМ. М.: Изд-во Моск. ун-
48
РИ, 2007, № 4
та, 1985. С. 255 - 262. 24. Тевяшев А.Д., Гибкина Н.В., Сидоров М.В. Об одном подходе к математическому моделированию плоских стационарных течений вязкой несжимаемой жидкости в конечных односвязных областях // Радиоэлектроника и информатика. 2007. № 2. С. 50 -57.
Поступила в редколлегию 20.12.2007
Рецензент: д-р физ.-мат. наук, проф. Колосов А.И.
Артюх Антон Владимирович, студент группы ПМ-04-1 факультета ПММ ХНУРЭ. Научные интересы: математическое моделирование, математическая физика, численные методы. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. 702-14-36.
УДК546.28 '
АРХИТЕКТУРА И ФУНКЦИОНАЛЬНОСТЬ ДВУХУРОВНЕВОЙ СИСТЕМЫ УПРАВЛЕНИЯ ВЫРАЩИВАНИЕМ СЛИТКОВ КРЕМНИЯ
ОКСАНИЧ А.П., ПЕТРЕНКО В.Р.,
ПРИТЧИН С.Э._______________________________
Рассматриваются вопросы управления процессом выращивания слитков кремния диаметром до 200 мм. Предлагается архитектура двухуровневой системы управления процессом выращивания слитков монокристаллического кремния. Описываются функции подсистем верхнего и нижнего уровней.
Введение
В настоящее время получение слитков кремния, выращиваемых методом Чохральского, качество и диаметр которых соответствует мировым требованиям, возможно только при высоком качестве функционирования системы управления процессом выращивания.
Современные системы управления должны обеспечивать воспроизводимость результатов технологических процессов при достаточно высоких значениях показателей качества выращенных кристаллов, которые определяются техническими условиями и требованиями заказчика. Это может быть достигнуто только путем использования в контурах управления математических моделей с настраиваемыми параметрами систем оперативной идентификации этих моделей и систем оптимизации режимов. При таком подходе до минимума сводится влияние на процесс «человеческого фактора».
Вопросу создания систем управления уделялось достаточно большое внимание с самого начала промышленного производства слитков кремния. Уровень автоматизации на том или ином этапе определялся степенью развития базы и средств вычислительной техники [1-4].
Г ибкина Надежда Валентиновна, канд. техн. наук, доцент кафедры прикладной математики ХНУРЭ. Научные интересы: экономический риск, актуарная математика, математическая физика, оптимальное управление динамическими объектами. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. 702-14-36.
Сидоров Максим Викторович, ассистент каф. прикладной математики ХНУРЭ. Научные интересы: математическое моделирование, математическая физика, теория R-функций и ее приложения. Увлечения и хобби: история культуры. Адрес: Украина, 61166, Харьков, пр. Ленина, 14, тел. 702-14-36.
Достаточно подробный анализ развития автоматизации управления процессом выращивания монокристаллов по методу Чохральского представлен в работе
[5].
Первые системы динамических моделей, описывающих процесс выращивания монокристаллов [4,5], были получены на основе классического подхода к решению задачи идентификации объекта управления. Модель процесса выращивания в пространстве состояний была предложена и развита Г. А. Сатункиным [7,8].
В работах [12-14] развит подход к моделированию передаточных функций, необходимых для регулирования диаметра растущего кристалла, на основе использования класса комбинированных моделей «передаточная функция - шум» [15]. Описание процедуры синтеза схемы регулирования с применением такой модели представлено в [16].
Современные подходы к проектированию эффективных автоматизированных систем управления процессом выращивания слитков по методу Чохральского достаточно подробно изложены в работах [9, 10], а принципы построения системы управления технологическим процессом выращивания бездефектного кремния с поддержанием заданных диаметра и температуры расплава описаны в [11, 12]. Примером реализации этих подходов и принципов может служить автоматизированная система «Кремень» [17]. Именно эксплуатация этой системы позволила выявить дальнейшие пути повышения эффективности систем управления процессом выращивания. Основные из них: замена локальной стойки управления установкой в целях обеспечения возможности реализации ею всех локальных контуров регулирования параметров процесса, оставляя за основной ПЭВМ настройку параметров регуляторов; обеспечение возможности с помощью одной ПЭВМ управлять работой нескольких установок; ведение информационной базы системы на основе использования промышленных СУБД; обеспечение информационной связи с другими технологическими участками производства; применение систем идентификации управляемых процессов и систем адаптации параметров используемых моделей.
РИ, 2007, № 4
49
