CC BY
52
Здесь пробелы в матрице С означают, что соответствующий ее элемент с¡}- =<х>, а в матрице Я - это Г] = ]■
Начальные операции У0 = {1, 4, 5}, конечные Ук = {3, 9} ■ Ранние сроки начала выполнения операций:
Т1рн = Т4рн = Т5рн = 0, Т2рн = 2, Т3рн = 14;
Т6рн = 1, Т7рн = 5, Т8рн = 1, Т9рн = 9, ?1р0н = 5,
соответственно ранние сроки окончания всех операций:
Т1ро = 2; Т2ро = 4; Т3ро = 17;
Т4ро = 4; Т5ро =1; Т6ро = 5; Т7ро = 9;
Т8ро = 4; Т9ро = 11; Т|ро = 14.
Критическое время tk = шах{С513, С519} = = С = 17
5,13 11 ■
Критические операции определяются с помощью матрицы Я
Г15,3 _ б; 1б,3 _ 10; Г20,3 _ 3
Таким образом, полученные в результате осуществления тернарных операций (1) и (2) матрицы С и Я позволяют определить все параметры сетевого графика.
В Ы В О Д
Разработан новый эффективный алгоритм отыскания критических путей в графе, не использующий графического представления последнего, а основанный на введенном понятии осуществления тернарных операций над матрицей дуговых весов.
Полученный алгоритм может использоваться для решения широкого круга задач, допускающих формулировку в терминах теории графов. Рассмотрены его приложения для решения задач упорядочения элементов конечного множества и задач сетевого планирования.
Л И Т Е Р А Т У Р А
1. Форд, Л. Р. Потоки в сетях / Л. Р. Форд, Д. Р. Фал-керсон. - М.: Мир, 1966. - 276 с.
2. Конвей, Р. В. Теория расписаний / Р. В. Конвей,
B. Л. Максвелл, Л. В. Миллер. - М.: Наука, 1975. - 360 с.
3. Floyd, R. W. Algorithm 97: Shortest Path / R. W. Floyd // Communication of ACM. - 1962. - № 5(6). - 345 p.
4. Корзников, А. Д. Тернарные операции в задачах оптимального преобразования сети / А. Д. Корзников // Современные прикладные задачи и технологии обучения в математике и информатике: сб. науч. ст. - Минск, 2004. -
C. 92-98.
5. Корзников, А. Д. Оптимизация системы маршрутов в транспортных сетях / А. Д. Корзников // Проблема прогнозирования и государственного регулирования социально-экономического развития. - Минск, 2004. -С. 46-49.
Поступила 6.11.2007
УДК 532.135:532.5.013.4
О РЕАЛИЗАЦИИ МЕТОДА КОРРЕКЦИИ ДАВЛЕНИЯ ПРИ РЕШЕНИИ УРАВНЕНИЙ ГИДРОДИНАМИКИ НЕСЖИМАЕМОЙ ЖИДКОСТИ
МАКАРОВ И. А.
ЗАО «БелХард»
Во многих исследований в области вычислительной гидродинамики значительную часть составляют работы, рассматривающие модель несжимаемой жидкости. Причина этого - мно-
гообразие технологических процессов, в которых применима данная модель, в частности при добыче нефти в Республике Беларусь, в том числе нефти с аномальной вязкостью
[1], где уместна модель неньютоновской несжимаемой среды [2-5]. При этом приходится сталкиваться с проблемой наличия давления в уравнениях гидродинамики, привносящего в задачу эллиптические свойства, требующие специальных мер при решении.
Приведем наиболее известные способы учета давления со сложностями и ограничениями их применения.
Метод искусственной сжимаемости. Уравнение неразрывности заменяется соотношением:
+ 8 — = 0
либо
dt
divv + гр = 0.
Центральным и сложным вопросом этого метода является правильное задание параметра искусственной сжимаемости е. Он должен быть достаточно мал, чтобы обеспечить близость модели к несжимаемой жидкости, но не настолько, чтобы нарушилось известное условие Куранта - Фридрихса - Леви, гласящего, что скорость распространения информации вдоль вычислительной области не должна превышать аналогичной физической или модельной скорости процесса [8, 9].
Метод коррекции давления (МКД) в виде наиболее известного [10-12] SIMPLE-алго-ритма (semi-implicit method for pressure-linked equations, полунеявный метод для решения уравнений, связанных с давлением) предполагает решение уравнения Пуассона для поправки к давлению. Тот факт, что это уравнение сравнительно трудоемко для численного решения и его нужно решать многократно на каждом временном слое, создает проблемы с производительностью вычислений [12]. Ниже мы приведем иную реализацию МКД, докажем ее сходимость.
Нахождение давления с помощью уравнения Пуассона. Метод, в отличие от двух предыдущих, не добавляет новых модельных свойств в уравнения гидродинамики, а по найденному на текущем шаге полю скоростей позволяет найти давление из соотношения [13]:
д
Ар =--.
дхх
] '
Как отмечается в [8], решение для давления в этом случае очень чувствительно к ошибке аппроксимации скоростей, так что реально это соотношение необходимо усложнять, вводя члены, пропорциональные дивергенции скорости, не равной фактически нулю ввиду ошибок аппроксимации.
Изложение процедуры МКД. Для конкретности рассмотрим плоский случай нестационарного течения несжимаемой жидкости, не делая каких-либо специальных предположений о связи напряжений ссо скоростями деформаций:
du du dv
--+ u--+ v— = -
dt dx dy dx dx
dP + !*n i dCT12 .
dy
dv dv dv
--+ u--+ v— = -
dt dx dy dy
dp dou da22 ,
dx
dy
du dv — +— = 0. dx dy
(1)
(2)
(3)
Применяя явную временную аппроксимацию по времени для уравнений движения, перепишем их в конечно-разностной форме, содержащей разностные операторы:
1!П+1 ИП I IT 1ИП t," ,-T"
uij = uij + (Lu (uij , vij , a1
11ij , a12ij , a22ij
) -
-KPj )At;
n+1 n /Т A-.n n n n n \ v„. = Vij + (Lv (uj , Vij , СТЩ , a12j , a22ij ) -
V V ■
-A yPj )At.
(4)
(5)
В предположении, что поле скоростей на временном слое п удовлетворяет условию неразрывности, полученное новое поле скоростей, вообще говоря, ему не удовлетворяет, так что необходимо уточнение. Для этого зададим давление на новом временном слое с помощью соотношения
и+1'
= Pin+1 -»divv;+1,
»>0,
(6)
и вновь определим скорости с помощью (4), (5) и скорректируем давление с помощью (6). Таким образом, приходим к следующей итерационной процедуре (вместо индексов т/, номеров узлов по горизонтали и вертикали, используем
далее индекс т, номер итерации):
= р ; (7)
Т
-А (АхРт-АхРТ-1); (8)
-А (АуРт-АуРт-1). (9)
При сходимости этой процедуры, очевидно, дивергенция стремится к нулю, и получаемые поле скоростей и давление на следующем временном слое близки к таковым для несжимаемой жидкости.
Докажем сходимость итерационной процедуры (7)-(9). Применяя операции Дх, Ду к соотношениям (8), (9) соответственно и складывая их, получим, имея также ввиду (7):
p
T+1
„ n+1 „ n+1
uT+1 = uT
и+1 ,,n+1
vT+1 = vT ■
divvZl - divvn+1 At
= ^(AxcdivvT"+1 + AyydivvT+l) (10)
- соотношение, являющееся конечно-разностной аппроксимацией уравнения теплопроводности для дивергенции:
d(divvn+1) dt
= ^A(divvTn+1),
(11)
где А - оператор Лапласа. Если на входной и выходной границах давления задаются постоянными, очевидно, мы должны решать (11) с граничными условиями, обращающими искомую дивергенцию в 0 на входной и выходной границах. Можно легко убедиться, что решение этой задачи асимптотически обращается в 0 во всей области при большом времени. Что касается конечно-разностной аппроксимации (10), в случае применения центральных разностей в А** и Ayy она хорошо изучена [8] и устойчива при выполнении определенного условия на соотношение шагов по времени, пространству, а также параметра итерационного процесса ц. Устойчивость итерационного процесса будет обсуждена ниже.
Численный эксперимент и обсуждение.
Для исследования реальной эффективности алгоритма (7)-(9) было выбрано нестационарное плоское осесимметричное течение ньютоновской жидкости по прямому каналу, при котором получаемые результаты легко анализировались с помощью простых физических соображений и сравнения со стационарной асимптотикой. Компоненты тензора напряжений конкретизировались следующим образом:
1
f
а_ =
2Re
ÖX:
dv __
dx,
Л
i, _ = 1,2. (12)
Задача (1)-(3), (12) решалась в замкнутой области 0 < х < 1, 0 < у < 1 (у = 0 - ось симметрии). Здесь х, у - безразмерные координаты, нормированные на полуширину канала; I - безразмерная длина канала в том же масштабе, скорости нормируются на масштаб скорости
задачи, равный
input
р
где pt
input'
р - размер-
ные входное давление и плотность жидкости. Далее имеем нулевые начальные условия для скоростей и давления и следующие граничные условия:
du dv dp о dy dy dy
при y = 0
(осевая симметрия),
u=v = 0 при y = 1
(условия прилипания),
p = 0 при x = 0,
P =
at 1 + at
при x = l,
(13)
(14)
(15)
(16)
т. е. давление на выходной границе выдерживается постоянным, а на входной границе оно нарастает, асимптотически стремясь к Р!„рШ ■
При этом течение должно асимптотически стремиться к известному стационарному решению с квадратичным профилем
(t, x, y) = R^ (1 - y2) p
input
(17)
(здесь коэффициенты согласованы с (12)). В численных экспериментах рЫри1 и I задавались равными 100 и 10 соответственно, чтобы порядок получаемых значений скорости был наиболее удобен для представления и анализа. Последнее граничное условие - условие на стенке для давления (у = 1). В случае ньютоновской жидкости обычно рекомендуется [8] использовать условие
дР п 1
— = 0 при у =1,
ду
(18)
следующее из уравнений Навье - Стокса, неразрывности для несжимаемой жидкости и условия прилипания, что и было применено в данном исследовании. В перспективе реализации рассматриваемого алгоритма применительно к неньютоновским жидкостям также использовался наиболее общий вид условия на стенке, следующий из (1)-(3), условий несжимаемости и прилипания:
др = дст12 ду дх
да-.
ду
при у = 1. (19)
Во всех описываемых ниже экспериментах использование (18) и (19) давало результаты, совпадающие в пределах аппроксимации.
Задача (1)-(3), (12)-(18) решалась конечно-разностной схемой, явной, первого порядка точности по времени. При этом градиент давления и дивергенция скорости описывались центральными разностями (второго порядка). Конвективные члены, равные нулю в данной постановке, но тем не менее способные давать возмущения, с экспериментальной целью удерживались и аппроксимировались разностями, направленными против потока [8, 10, 11]. Результаты, представленные на рис. 1, получались на сетке с 53 узлами вдоль оси х и 23 вдоль оси у, шаг по времени & = 0,0001. На такой сетке использованная реализация процедуры МКД сходилась за 20-30 итераций при значении параметра ц около 15. Отметим, что еще лучшая сходимость достигалась на «квадратной» сетке с равными шагами по обеим ко-
ординатам: вплоть до 15 итераций при ц» 30. Значение ц, обеспечивающее наилучшую сходимость, в то же время являлось и пороговым: при дальнейшем его увеличении происходила потеря сходимости.
1,0
у
0,8 0,7 0,6 0,5 0,4 0,3 0,2 0,1
Асимптотический профиль
0
1,0 2,0 3,0 4,0 5,0 и
Рис. 1. Профили течения по прямому каналу при Яе = 0,5 и а = 1 в моменты времени 1, 2, 3 (в порядке нумерации кривых) в сопоставлении с асимптотическим стационарным течением
При числах Рейнольдса, больших единицы (особенно если также а > 1), спустя некоторое время (несколько единиц) после начала течения решение «тонуло» в ошибках осцилляционного характера. К числу возможных причин этого следует отнести достаточно большую крутизну профиля скорости при Яе > 1 (рис. 2).
1.2 1,1
1,0 0,9 0,8 0,7 0,6 0,5 0,4 0,3 0,2 0,1
0 1,0 2,0 3,0 4,0 5,0 6,0 7,0 8,0 9,0 10
Рис. 2. Профили течения по прямому каналу при Яе = 1 и а = 1 в моменты времени 1, 2, 3 (в порядке нумерации кривых) в сопоставлении с асимптотическим стационарным течением
Решение, резко меняющееся во времени и пространстве, как известно, «обогащается» гармониками высокого порядка, не разрешаемыми конечномерной аппроксимацией. Путями борьбы с этим являются специфический подбор схем, а также использование сетки, более мелкой в области резкого изменения решения. Характерно, что в таких ситуациях итерационный процесс (7)-(9) часто продолжал сходиться до тех пор, пока возмущения не принимали катастрофический характер.
Наконец, была проанализирована стабилизация течения со временем по мере стабилизации входного давления. В качестве меры стабилизации N рассматривалось осредненное по всем расчетным узлам абсолютное значение аппроксимации временной производной компонент скоростей (иЩ+1 - иЩ )/А^ (ущ+1 - у^ )/At■
Рис. 3 показывает, что в диапазоне небольших чисел Рейнольдса временная зависимость стабилизации жидкостей практически одинакова и определяется условиями течения, а не природой жидкости, что согласуется с результатами [14]. 1.0 0,9 0,8 0,7 0,6 0,5 0,4 0,3 0,2 0,1
0 ...................
1 3 5 7 9 11 13 15 17 19 21 22
Рис. 3. Зависимость от времени степени отличия течения от стационарной асимптотики для трех значений числа Рейнольдса и а = 0,5
В Ы В О Д
Представленный сходящийся алгоритм коррекции давления в уравнении Навье - Стокса, как показано, работоспособен и эффективен в большом диапазоне параметров решенной гидродинамической задачи. Он может быть при-
меним и к жидкостям со сложной вязко-упругой реологией, являющимся предметом продолжения представленного исследования. Полученные же результаты применения алгоритма к простому и ясному примеру течения демонстрируют точность и экономичность предложенной реализации метода коррекции давления.
Л И Т Е Р А Т У Р А
1. Журавель, Н. Г. Об особенностях гидродинамического моделирования транспорта высоковязкой нефти / Н. Г. Журавель, И. А. Макаров, А. А. Маханек // Поиски и освоение нефтяных ресурсов Республики Беларусь. -Гомель: БелНИПИнефть, 1999. - Т. 3. - C. 270-278.
2. Макаров, И. А. Моделирование фазы установления реометрического течения вязко-упругой жидкости / И. А. Макаров // ИФЖ. - 2001. - Т. 74, № 3. - C. 173-176.
3. Макаров, И. А. Об устойчивости реометрического течения вязкоупругой жидкости к сдвиговым возмущениям / И. А. Макаров // Известия Российской Академии наук. Серия Механика жидкости и газа. - 2003. - № 2. -C. 6-12.
4. Time-Weissenberg number superposition in 4:1 planar contraction flow of a viscoelastic fluid / M. K. Ju [et al.] // Journal of the Society of Rheology. - 2005. - Vol. 33, № 4. -P. 191-207.
5. Joo, Y. L. Observations of purely elastic instabilities in the Taylor-Dean flow of a Boger fluid / Y. L. Joo // Journal of Fluid Mechanics. - 1994. - Vol. 262. - P. 27-42.
6. Флетчер, К. Численные методы на основе метода Галеркина / К. Флетчер. - М.: Мир, 1988. - 352 c.
7. Shen, Jie. Pseudo-compressibility methods for the unsteady incompressible Navier-Stokes equations / Jie Shen // Proceedings of the 1994 Beijing Symposium on Nonlinear Evolution Equations and Infinite Dynamical Systems. - Beijing: Zhong Shan University Press, 1997. - P. 803-814.
8. Роуч, П. Вычислительная гидродинамика / П. Роуч. -М.: Мир, 1976. - 612 с.
9. Поттер, Д. Вычислительные методы в физике / Д. Поттер. - М.: Мир, 1975. - 392 с.
10. Chung, T. J. Computational fluid dynamics / T. J. Chung. -Cambridge University Press, 2002. - 1022 p.
11. Ferziger, J. Computational methods for fluid dynamics / J. Ferziger, M. Peric. - Springer, 2002. - 431 p.
12. Anderson, J. D. Computational fluid dynamics / J. D. Anderson. - McGraw-Hill, 1998. - 563 p.
13. Ландау, Л. Д. Теоретическая физика в 10 т. / Л. Д. Ландау, Е. М. Лифшиц // Гидродинамика. - М.: Физматлит, 2002. - Т. 6. - 736 с.
14. Макаров, И. А. Моделирование реометрического течения ньютоновской жидкости методом конечных элементов / И. А. Макаров // ИФЖ. - 2000. - Т. 73, № 5. -C. 927-931.
Поступила 10.01.2008
