CC BY
25
УДК 533.6
С. Б. Березин, В. М. Пасконов
ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ НЕСТАЦИОНАРНОГО ТРЕХМЕРНОГО ТЕЧЕНИЯ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ В КАНАЛЕ КВАДРАТНОГО СЕЧЕНИЯ НА ОСНОВЕ МОДЕЛИ НАВЬЕ-СТОКСА1
(кафедра математической физики, лаборатория моделирования процессов тепломассопереноса факультета ВМиК)
Введение. В наших работах [1, 2] был предложен метод численного исследования течений вязкой несжимаемой жидкости в рамках классической модели Навье-Стокса на основе системы уравнений, сохраняющих прямую и обратную связь динамических уравнений и уравнения энергии и при численных расчетах обеспечивающих высокую точность выполнения закона сохранения массы. Такая система уравнений получается из полной системы уравнений Навье-Стокса для вязкой сжимаемой жидкости, состоящей из трех уравнений количества движения, уравнения неразрывности, уравнения энергии, записанного относительно температуры, и уравнения состояния. Система уравнений замкнутая: в ней 6 уравнений и 6 неизвестных. Если предположить, что в этой системе плотность постоянна, то уравнение состояния дает линейную связь давления и температуры и, следовательно, существует проблема выбора: одно из пяти дифференциальных уравнений необходимо при расчете отбросить. При отбрасывании уравнения энергии теряется и уравнение состояния — в результате получаем классический подход: описание изотермического течения с давлением, которое находится из уравнения Пуассона. Более естественно, на наш взгляд, вести расчеты задач на основе системы уравнений количества движения и уравнения энергии, которые выводятся, как известно (см. [3-5]), с использованием закона сохранения массы. Вследствие этого он выполняется с достаточно высокой точностью при расчетах. Настоящая статья посвящена численному исследованию развития течения при инжекции вязкой несжимаемой жидкости в канал квадратного сечения, заполненный той же жидкостью. Рассмотрение такой задачи связано с изучением нестационарных трехмерных структур течений, которые исследовались в многочисленных теоретических работах на основе теории пограничного слоя, параболизации уравнений Навье-Стокса и использования различных моделей турбулентности [6]. Среди работ по этой проблематике следует отметить монографию [7], содержащую не только теоретические, но и экспериментальные результаты.
Постановка задачи. Рассматривается канал квадратного сечения длиной X£ (рис. 1), заполненный в начальный момент времени покоящейся вязкой несжимаемой жидкостью.
Рис. 1
Течение в канале возникает за счет инжекции такой же жидкости при х = 0. Предполагается, что за выходным сечением канала находится пространство, также затопленное жидкостью. В качестве характерной длины выбрана длина стороны квадратного сечения канала, характерной скорости — максимальная скорость инжектирующей жидкости при х = 0.
1 Работа выполнена при финансовой поддержке РФФИ, проекты № 04-01-0332, № 05-07-90378.
Ф = 2
ди\2 / ду\2 /дю дх) \ду) \ду
Для расчетов использовалась система уравнений, полученная из полной системы уравнений Навье-Стокса для вязкого, сжимаемого, теплопроводного, совершенного газа, записанной в безразмерных величинах (см. [5]), где плотность, давление и температура относятся к соответствующим величинам в инжектируемом потоке. При этом предполагается, что коэффициенты динамической вязкости и теплопроводности постоянны, массовые силы отсутствуют, плотность постоянна и равна 1, что приводит к линейной связи температуры и давления. Таким образом, полная система уравнений в безразмерных величинах, состоящая из трех уравнений количества движения, уравнений энергии, неразрывности и состояния, может быть записана в следующем виде:
ди ди ди ди 1 дТ 1 / д2и д2и д2и\
дЬ дх ду дг 7М2 дх Яе\дх2 ду2 дг2)'
ду ду ду ду 1 дТ 1 ( д2 у д2у д2у\
дх ду дг 7М2 ду Яе\дх2 ду2 дг2/1
диз диз диз диз 1 дТ 1 / д2 из д2из д2из\
дЬ дх ду дг 7М2 дх Яе \ дх2 ду2 дг2 )'
дТ дТ дТ дТ 1 (д2Т д2Т д2Т\ . М2
+ + ^ = + ^ + + 1 (4)
/ ди ду\2 / ди дги\2 / ду дги\2 \ду дх) \дг дх) \дг ду) ' ди ду дги
дх ду дг '
" = <7>
Здесь введены общепринятые обозначения: £ — время, и,у,ги — компоненты вектора скорости соответственно по координатным осям х,у,гф, р — давление, Т — температура, Ф — диссипативная функция. Параметры подобия: число Рейнольдса Яе = у ^ , число Маха М = число Прандтля
г> 2 7Рп " с'р (
гг= -у-] квадрат скорости звука а = отношение удельных теплоемкостеи 7 = -¡§- (величины со
штрихом — размерные).
В целом система (1)-(7) является переопределенной. Однако уравнение (7) не дает основания сомневаться в том, что давление (температура) должно определяться в процессе решения системы уравнений (1)-(5), а уравнение неразрывности (6) должно использоваться в качестве контрольного соотношения. Как было определено выше, в начальный момент времени £ = 0 канал заполнен покоящейся вязкой несжимаемой жидкостью — внутри канала и на его стенках и = у = ги = 0 при начальной температуре Т = 1. Течение в канале возникает за счет инжекции такой же жидкости, которая осуществляется при х = 0 по всему начальному сечению канала: и = 1, у = ги = О, Т = 1 для £ ^ 0. В качестве граничных условий на стенках канала использовались условия прилипания и = у = ги = Ои условие теплоизоляции для температуры (которые для давления также имеют естественную физическую интерпретацию):
— - О
дп '
где п — нормаль к соответствующим поверхностям стенок канала.
В плоскости выходного сечения канала х = Ьх использовались "мягкие" граничные условия для всех искомых функций
д2и д2у д2т д2Т
дх2 дх2 дх2 дх2
= 0.
Для численного решения сформулированной выше задачи использовался широко известный вычислительный метод покоординатного расщепления [8] для системы уравнений (1)-(5). При этом диссипа-
тивная функция представлялась в виде:
Ф = Фж + ФУ + Ф
х \ ^ у \ ^ z 1
(8)
(П)
(10)
(9)
Процесс расчета при переходе от предыдущего временного слоя к последующему осуществлялся в три этапа на дробных временных слоях: на первом из них решалась система четырех уравнений разностным методом для "ж-расщепления", на втором этапе — для "у-расщепления" и на третьем этапе (на следующем временном слое) — для "^-расщепления". На каждом из этих этапов проводились итерации по нелинейности с использованием уже полученных значений на предыдущем дробном шаге. В результате на следующем временном слое получалось только первое приближение решения. Окончательное решение на временном слое достигалось с помощью последующих "глобальных" итераций, содержащих все три этапа расщепления со своими внутренними итерациями. Следует отметить, что число "глобальных" итераций по нелинейности всегда превышало число итераций на дробных шагах. Точность решения проверялась путем вычисления максимума модуля интеграла от дивергенции вектора скорости по объему каждой ячейки трехмерной сетки. Численный метод, описанный выше, был реализован программно на языке Си и применялся для расчетов на многопроцессорном вычислительном комплексе "Regatta" с использованием технологии ОрепМР [9].
Обсуждение результатов расчетов. Мы будем представлять результаты расчетов в виде построенных на компьютере рисунков, на которых изображаются поверхности продольных и поперечных составляющих вектора скорости, а также поля направлений вектора скорости в поперечных сечениях канала при различных значениях х. Чтобы проследить развитие течения, приводятся значения безразмерного времени. Эти рисунки сделаны с помощью системы компьютерной визуализации результатов научных исследований [10]. Процесс развития течения при одинаковых начальных и граничных условиях существенно зависит от числа Рейнольдса.
Первоначально остановимся на описании результатов расчетов при числе Рейнольдса Re = 100, числе Маха М = 0,03 и длине канала Lx = 3,0. Расчеты были проведены на сетке 301 X 101 X 101, при этом максимум модуля интеграла от дивергенции по всем ячейкам сетки не превышал 0,00004, что показывает высокую точность выполнения закона сохранения массы. На рис. 2 при t = 1,0 приводятся поверхности продольной и и поперечной v составляющих скорости над плоскостью сечения канала для различных значений ж, по которым можно судить о том, как развивается возмущение вдоль канала. В начале канала происходит формирование "параболоида и" на квадратном основании, при этом максимальное значение продольной составляющей скорости остается в центральной части канала равным и = 1,0 (х = 0,6). После х = 0,6 с увеличением х начинается уменьшение и по всему сечению канала, меняется конфигурация поверхности "параболоида" (х = 1,0). При дальнейшем увеличении х "параболоид" становится более наполненным (х = 1,5), максимальная величина продольной составляющей скорости в центре канала с увеличением х продолжает расти и в конце канала достигает значения 0,6. Но не только изменение продольной составляющей скорости вдоль канала характеризует структуру течения. Рассматривая изменение поперечных составляющих скорости v и w, следует отметить, что их поверхности над плоскостью сечения канала повернуты на 90 градусов, а их величины меняются в одинаковых пределах. Заметим, что это относится ко всем расчетам, представленным в настоящей статье. При исследовании поведения v и w вдоль канала при i = 1,0 оказалось, что на интервале 0,7 < х < 1,5 величины v и w соизмеримы с величиной продольной компоненты, что приводит в этом интервале изменения х к образованию сложной структуры течения, о характере которой можно судить по представленным на рис. 3 полям проекции вектора скорости на плоскость поперечного сечения канала при х = 0,6 и х = 1,0.
Здесь векторы изображены в виде вытянутых треугольников, цвет которых меняется от белого (минимальное значение величины вектора) до черного (максимальное значение). Наблюдаются четкая симметрия поля относительно двух диагоналей квадрата и четыре сопряженных между собой "полу-
Рис. 2
вихря" в каждом равнобедренном треугольнике. Максимальные значения величин векторов находятся в центре квадрата, минимальные — вблизи стенок канала. При переходе от х = 0,6 к х = 1,0 происходит некоторая трансформация поля: увеличивается область малых величин векторов вблизи стенок и особенно в углах, где образуются вихревые зоны двух типов: в левом нижнем и правом верхнем — зоны "разлета", а в правом нижнем и в левом верхнем — зоны типа "фокуса". При этом отметим, что векторы полной скорости имеют общее направление вдоль канала. В последующие моменты времени такие вихревые структуры не наблюдались. Отметим также, что величина максимальной со-
ч> и и я
<7 \ч- Р"- ч'ч- ч-'- !>- Г-- Г-ч- о !•>- ¡>- ::ч>.
КР' ^ О- £> £>■ Ьь- О Сз* £>•
<5 у р ¡у
у* Л/
о- Й». &>• &>■ й*. О- О ¡.о ^ ^ ^ ^/ч
V V- Ъ ¿1
С? \> <1 V а-<3 <\, ^ Са*йа. Йа. £
в«.а»а.в»-^
I? $ У Ч ^ ^^^^^ ^ ^ 4 ~
17 ? 1 <-1 ^ Ч ^ Л <
р V ^ Я Ч ^ ЧАЛ
? V 1 1 1 ^ ^ ^ ^
У 7
У 7 V
V «л- к»'»- ^^
й IV
й Ь>
/! а
¿3 А
4 4 4 4
4 4 4
4 4 4 4 & й
4 4 4 А & к
4 4 4 4 ' '
4 ,1 4 /I
4 4 .1 1
У 7 У ? ?
1Ш1
? ? ? ?
V ? ? ? ? ?
7 ?
Р I
? '
4 4 4
444441111 4 Н I 1 4 I з 4 ккккккккк к к к
УЧ.... 7111111111 л
71111111111
£ £ иииш
р р р ? 0> > Ч V V V к
> > * ч > > V
^ ^ ^ЙЙ -«а -««-»«8 -45 ^ V
? Г ,
>н
7 ? V Г Г с V ? 7 ¡? Г ^ Р Р Я V ? Р ^ ^ ^ ^
1 ^ ? А#
р ^ ^ ''ч/ "ч: и
^ ^ -¿й -53 -«¿3 -43 ^ -«а Ч >
^ & чгз ^сз оЕЗ-^53 -«3 чй -С- -«3 ^ ^ ^
-'/ -'1 ^ V Ч? -г-: ^ -л- -.-ч <з .-г.- --ч /■
Г.'П «'(» 'Рп Ч-Л Г_ 'п !.'(! ^''ц
^ ^ I. Ь. к к &
ЪЬ Ь. I. I к к
% % I, (Ь й
Ь. Ь 6 Й
>>>>VVVЬ Ь М
^ х« > % V V Ь & ¿1
^ч; Ъ 'ч; V
¿1 4
4 С-
С'
■1. ^
Сч- чЧ- Ч/ Ч- 'V -.-ч ч/ -г/ ч/
-С/ 4
л <4, й» ^- С» о- Е»
^ ^ ч. члч ччч чч, -чч. Ч:^ч. ;.>--:чччч- ¡;;----чччс>"ч-;:'' С^ <7 <)'
<ч Ч;..^ <Ч /Чч -:'_ч_ /ч, ^ '---ч ^ч- ч-ч ^ -5" ч^' : ч/ ч/ ч/ ч/ /1
^ ^ в7 ^ 4 <1 4
<4 Ч Ч Ч
Ч Ч Ч Ч Ч Ч А ^ ч ^ ^ ^ ^
V*".*"
1 -и 1 и и 1 а ,
Ч Ч 1 1 1 1 1 1 ^ 1 Л А ^ * * * 5 1 1
^•9 1111111111 а 1 ] 1
ЧЧЧЧ 11111111 1 4 к 4 1
? 7 ? I ? 1 ? 7 ? Т М 11 ¡11 Ик
г. Ч/уф* & ^ ^
7 ?????»?! т ? т т ч * 1 Vг г г г г г г г
?????????ГГГГ 5? Г ? Г Г Г Г Г Г Г Г I, к
1? V Р ? Г Г Г Г Г
17 р г Г Г "
Р р Г Р г ? ^ Р ? Vр ? ? РРР Р Р Р Р > А ^
рр г * * >
р р рр А А
ф* ^ <? 4 с/ 4 <7
4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 й 4 4 4 4 4 4 4 4 444 44 4 4 4 44 4 4 А 4 4 4 А
а а А а А
А ¿АА & ¿&
& й Й & & & & &
% к к к к к к & к & & & Ь
Мкккккккк&й
к к к к к & ь
Ь к к к к к ь
Ъ-Ы к к к к
ь. ь
> V ь-
> V % V V V
Ъ \ V 'Ь
-л—
9 ¿> ¿> & А
/;. р р /:■. ,чЧ> -Сч -ЧЩ -:;ч -ч;_, чч-, ч:;7 4., ч ... Чч \'ч V .
'> А А /> А 4> ч:-'': .чЧ.'; ч-ч .-.'ГА чЧ-; ^ чг; чг7 -г. Ч-^ ч;> ЧЧ '-ч \ч
>-, ¿¡> ^ ¿а ж» «я -я <» <з -=)•<!--
^ ё1о 8^0 & ^о ё^з
1 "чг х? V V V V V
N5
о"
Рис. 3
ставляющей скорости и на выходе из канала со временем начинает падать (например, при £ = 2,5 и = 0,1). Затем в конце канала образуется практически почти застойная зона, левая граница которой постепенно сдвигается в сторону начала канала. В результате при £ = 15,0 происходит установление течения, которое характеризуется переходом от инжекции при х = 0 к полному покою при х > 2,0, где и = V = т = 0. В области 0 < х < 2,0 при £ = 15,0 продольная компонента скорости равномерно убывает от 1 до 0 и играет превалирующую роль по сравнению с поперечными компонентами, что приводит к образованию стационарной застойной зоны при х > 2,0. При этом температура (давление) во всем канале незначительно повышается. Таким образом, происходит процесс блокировки канала.
Рис. 4
Последующие расчеты при более высоких числах Рейнольдса показали, что блокировка канала при описанных выше начальных и граничных условиях происходит всегда на начальном участке канала и практически не зависит от его длины. На рис. 4 представлены результаты расчета течения в канале длины Ьх = 1,5 при числе Рейнольдса Яе = 1000 и числе Маха М = 0,03 на сетке 226 X 151 X 151. При этом максимум модуля интеграла от дивергенции по всем ячейкам сетки не превышал 0,000018. На верхней части рис. 4 при £ = 1,0 представлены поверхности продольной составляющей скорости и для х = 0,75; 0,9, где и > 0 во внутренних точках области, и для х = 1,2; 1,5, где вблизи стенок канала продольная составляющая скорости становится отрицательной.
Изменение решения для этого варианта происходит значительно более интенсивно как по времени, так и по пространству. В начале канала векторы полной скорости во всех внутренних точках поперечных сечений сохраняют направление вдоль увеличения продольной координаты х. В центральной области сечения продольная составляющая скорости и сохраняет свое максимальное значение. Однако при х = 0,75 определяющую роль играют уже поперечные компоненты, что приводит к формированию четырех взаимно сопряженных центральных вихрей, центры которых находятся вблизи центров стенок канала. Это показано на двух нижних рисунках, где представлены поля проекции вектора скорости на плоскость поперечного сечения канала при х = 0,75 и х = 0,9. При этом так же, как и при Re = 100, соблюдается четкая симметрия относительно диагоналей квадрата сечения. В результате исследования плавного изменения вдоль по х направления вектора полной скорости вблизи центра канала было обнаружено винтовое движение жидкости вдоль х по часовой стрелке, если смотреть со стороны сечения х = 0. На рисунке при х = 0,9 видно, что центры четырех больших вихрей сместились к центру канала, а в углах канала образовалось по два вихря, которые разрешают "конфликт" больших вихрей в угловых областях, где величины скорости существенно меньше характерной скорости. Здесь следует отметить, что вихревые структуры в углах канала согласуются с экспериментальными лабораторными данными, приведенными в [7, 11]. При этом отмечаются трудности постановки экспериментов по исследованию таких вторичных течений, связанные с тем, что величины скоростей в угловых конфигурациях на порядок меньше, чем вне приграничного слоя. Подобные трудности, возникающие и в вычислительном эксперименте, разрешаются при расчетах существенным увеличением числа узлов разностной сетки, а при анализе численных результатов с помощью графических методов — путем выделения соответствующих фрагментов данных.
В отличие от течения при Re = 100, при Re = 1000 для х > 1,1 сначала в углах канала (см. поверхность компоненты скорости и для х = 1,2), а затем вблизи стенок компоненты скорости и становятся отрицательными, и впоследствии область положительных значений и в центре канала постепенно сокращается и на выходе из канала (х = 1,5) и < —0,05 во всех внутренних точках поперечного сечения, при этом структура поля проекций полной скорости состоит из четырех сопряженных между собой вихрей и развитых двувихревых вторичных течений в углах. При этом векторы скорости направлены внутрь канала, а их модули в среднем на два порядка меньше характерной величины скорости на входе в канал.
В заключение отметим, что приведенные результаты расчетов нестационарных трехмерных течений вязкой несжимаемой жидкости в канале квадратного сечения, заполненного в начальный момент времени такой же жидкостью, дали возможность исследовать процесс развития течения. При этом было установлено возникновение в канале сложного нестационарного четырехвихревого течения с вторичными вихрями вдоль двухгранных углов канала, структура которых согласуется с экспериментальными данными. Обнаружено винтовое движение жидкости вдоль канала в процессе установления стационарного режима блокировки канала, возникающей в силу предположения о том, что за выходным сечением канала находится пространство, затопленное жидкостью. Хотя представленные здесь решения сложны и богаты вихревыми структурами, на наш взгляд, решение при Re = 1000 не моделирует турбулентное течение не только потому, что происходит установление решения по времени, но и потому, что, используя для решения уравнений Навье-Стокса разностные методы, невозможно гарантировать согласие вычислительных результатов и результатов лабораторных экспериментов по числу Рейнольдса [12, 13].
СПИСОК ЛИТЕРАТУРЫ
1. Березин С. Б., Пасконов В.М. Численное исследование вдува вязкого несжимаемого газа в плоский канал на основе уравнений Навье-Стокса // Вычислительные методы и программирование. НИВЦ МГУ. 2003. Т. 4 (1). С. 5-17.
2. Березин С. Б., П ас ко н о в В. М. Неклассические решения классической задачи о течении вязкой несжимаемой жидкости в плоском канале // Прикладная математика и информатика. М.: МАКС Пресс, 2004. № 17. С. 13-30.
3. ШлихтингГ. Теория пограничного слоя. М.: Наука, 1974.
4. ЛойцянскийЛ.Г. Механика жидкости и газа. М.: Наука, 1973.
5. Абрамович Г. Н. Прикладная газовая динамика. Ч. 1. М.: Наука, 1991.
6. Трехмерные турбулентные пограничные слои. М.: Мир, 1985.
7. Корнилов В. И. Пространственные пристенные турбулентные течения в угловых конфигурациях. Новосибирск: Наука, 2001.
8. Самарский А. А., Н и к о л ае в Е. С. Методы решения сеточных уравнений. М.: Наука, 1978.
9. Воеводин В. В., Воеводин Вл. В. Параллельные вычисления. СПб.: БХВ-Петербург, 2002.
10. Березин С. Б., Пасконов В.М. Компонентная система визуализации результатов расчетов на многопроцессорных вычислительных системах // Тр. конф. "Высокопроизводительные вычисления и их приложения". Черноголовка, 2000. С. 202-203.
11. Johnston J. On the three-dimensional boundary layer generated by secondary flow //J. Basic Engin. Trans. ASME. 1960. Ser. D. 1.
12. Пасконов В.М. Модификация уравнений Навье-Стокса для расчета вязких течений разностными методами // Обратные задачи естествознания. М.: Изд-во МГУ, 1997. С. 189-197.
13. Р a s к о п о v V. М. Modified Navier-Stockes equation for finite-difference computation of viscous flow // Сотр. Math, and Modeling. 1997. 8. N 4. P. 400-407.
Поступила в редакцию 10.10.05
УДК 519.2
0. В. Шестаков
НЕЛИНЕЙНАЯ РЕГУЛЯРИЗАЦИЯ ОБРАЩЕНИЯ ЭКСПОНЕНЦИАЛЬНОГО ПРЕОБРАЗОВАНИЯ РАДОНА
(кафедра математической статистики факультета ВМиК)
1. Введение. Задачи однофотонной эмиссионной томографии приводят к проблеме обращения так называемого экспоненциального преобразования Радона в Л2 [1], которое отличается от классического преобразования Радона наличием экспоненциального множителя у интегрируемой функции. В настоящее время существует множество методов восстановления функции по ее экспоненциальному преобразованию Радона, среди которых методы фильтрации обратных проекций и итерационные методы.
Метод разложения в вейвлет-ряды с последующей пороговой обработкой вейвлет-коэффициентов, предлагаемый в данной работе, основан на относительно недавно введенном объекте — дискретном вейвлет-преобразовании. Преимущество вейвлетов заключается в том, что они локализованы как в пространственной, так и в частотной области, что делает их хорошо подходящими для анализа и обработки нестационарных сигналов и изображений.
Предлагаемый в данной работе метод был разработан и применен для классического преобразования Радона в работе [2] и исследован численно в работе [3].
2. Экспоненциальное преобразование Радона. Пусть V = 51 X Л, где 51 — единичный круг в Л2. Экспоненциальным преобразованием Радона для ¡1 6 И. и функции /(ж, у), имеющей компактный носитель (без ограничения общности будем считать этим носителем круг II с единичным радиусом и центром в начале координат), называется преобразование вида
оо
J е»{х>е±)${х)<1х= J е^ЦзО + М^М, (0, в) 6 V.
{х,в) = 8 -ОО
Заметим, что если ¡1 = 0, то мы имеем классическое преобразование Радона. Сопряженным оператором к оператору является оператор обратного проецирования, задаваемый преобразованием вида
Е*д(х) = J е"{х'в±)д(0, {х,в))йв, х Е В2.
Б1
