Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения

Oписан алгоритм численного решения полной системы нестационарных уравнений Навье-Стокса для несжимаемой жидкости, который базируется на методе конечных объемов. Уравнения решаются в неподвижной декартовой системе координат. Используются разностные схемы, имеющие второй порядок точности по пространст...

Повний опис

Збережено в:

Бібліографічні деталі
Опубліковано в: :	Прикладна гідромеханіка
Дата:	2010
Автор:	Малюга, В.С.
Формат:	Стаття
Мова:	Russian
Опубліковано:	Інститут гідромеханіки НАН України 2010
Онлайн доступ:	https://nasplib.isofts.kiev.ua/handle/123456789/87748
Теги:	Додати тег Немає тегів, Будьте першим, хто поставить тег для цього запису!
Назва журналу:	Digital Library of Periodicals of National Academy of Sciences of Ukraine
Цитувати:	Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения / В.С. Малюга // Прикладна гідромеханіка. — 2010. — Т. 12, № 4. — С. 45-62. — Бібліогр.: 62 назв. — рос.

Репозитарії

Digital Library of Periodicals of National Academy of Sciences of Ukraine

id	nasplib_isofts_kiev_ua-123456789-87748
record_format	dspace
spelling	Малюга, В.С. 2015-10-24T16:17:27Z 2015-10-24T16:17:27Z 2010 Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения / В.С. Малюга // Прикладна гідромеханіка. — 2010. — Т. 12, № 4. — С. 45-62. — Бібліогр.: 62 назв. — рос. 1561-9087 https://nasplib.isofts.kiev.ua/handle/123456789/87748 532.516 Oписан алгоритм численного решения полной системы нестационарных уравнений Навье-Стокса для несжимаемой жидкости, который базируется на методе конечных объемов. Уравнения решаются в неподвижной декартовой системе координат. Используются разностные схемы, имеющие второй порядок точности по пространству и времени. Для дискретизации конвективных членов используется TVD схема Chakravarthy-Osher. Численный алгоритм разработан для трехмерной неструктурированной сетки. Полученная система нелинейных алгебраических уравнений линеаризуется и решается по алгоритму PISO. Для решения соответствующих систем линейных алгебраических уравнений применяются итерационные солверы, использующие методы сопряженных/бисопряженных градиентов с предобусловливанием, которые можно найти в свободно распространяемых библиотеках, доступных в репозитарии NetLib [62]. Описанный численный алгоритм использован для прямого численного моделирования течения жидкости в канале с двумя последовательно расположенными стенозами при различной ширине канала в области между стенозами. Установлено, что при достаточно малой ширине межстенозной части канала в межстенозной области генерируются симметрично расположенные вихри. А при достаточно большой ширине в определенном диапазоне чисел Рейнольдса два ряда вихрей, генерируемых течением в межстенозной области, организуются в шахматном порядке. Oписано алгоритм чисельного розв'язання повної системи нестаціонарних рівнянь Навьє-Стокса для нестисливої рідини, що грунтується на методі скінчених об'ємів. Рівняння розв'язуються в нерухомій декартовій системі координат. Використовуються різничні схеми, що мають другий порядок точності у просторі і за часом. Для дискретизації конвективних членів використовується схема Chakravarthy-Osher. Чисельний алгоритм розроблено для тривимірної неструктурованої сітки. Отримана система лінійних алгебраїчних рівнянь лінеарізується і розв'язується за алгоритмом PISO. Для розв'язання відповідних систем лінійних алгебраїчних рівнянь застосовуються ітераційні солвери, що базуються на методах спряжених/біспряжених градієнтів з передобумовленням, які можна знайти у відкритих бібліотеках, доступних в репозитарії NetLib [62]. Описаний чисельний алгоритм використано для прямого чисельного моделювання течії рідини в каналі з двома послідовно розташованими стенозами за різної ширини канала в області між стенозами. Встановлено, що за досить малої ширини міжстенозної частини канала в міжстенозній області генеруються симетрично розташовані вихори. А за досить великої ширини у певному діапазоні чисел Рейнольдса два ряди вихорів, що генеруються течією в міжстенозній області, організуються в шаховому порядку. A numerical solution algorithm for the unsteady incompressible Navier-Stokes equations is described in the paper. The discretization procedure is based on the finite volume method. The equations are solved in the fixed Cartesian coordinate system. Stabilized and bounded differencing schemes, second-order accurate in both space and time, are used. The TVD differencing scheme of Chakravarthy-Osher is employed for the discretization of the convective terms. The numerical algorithm is designed for three-dimensional arbitrarily unstructured meshes. The system of non-linear algebraic equations obtained from the discretization procedure is linearized and solved with the PISO algorithm. The corresponding systems of linear algebraic equations are solved using the iterative solvers based on the preconditioned conjugate/biconjugate gradient methods. These solvers are available from the repository NetLib [62]. The numerical algorithm described in the paper is employed for the direct numerical simulation of the flow in a duct containing two sequential stenoses. The calculations are performed for various widths of the inter-stenoses segment of the duct. It is established that if the duct is sufficiently narrow between the stenoses, in some range of Reynolds number, the vortices generated in the inter-stenoses domain are situated symmetrically. But if the duct is sufficiently wide between the stenoses, the two rows of the vortices generated in the inter-stenoses domain are chequer-wise situated. ru Інститут гідромеханіки НАН України Прикладна гідромеханіка Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения Numerical investigation of the flow in a duct with two serial stenoses. Algorithm of the solution Article published earlier
institution	Digital Library of Periodicals of National Academy of Sciences of Ukraine
collection	DSpace DC
title	Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения
spellingShingle	Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения Малюга, В.С.
title_short	Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения
title_full	Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения
title_fullStr	Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения
title_full_unstemmed	Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения
title_sort	численное исследование течения в канале с двумя последовательно расположенными стенозами. алгоритм решения
author	Малюга, В.С.
author_facet	Малюга, В.С.
publishDate	2010
language	Russian
container_title	Прикладна гідромеханіка
publisher	Інститут гідромеханіки НАН України
format	Article
title_alt	Numerical investigation of the flow in a duct with two serial stenoses. Algorithm of the solution
description	Oписан алгоритм численного решения полной системы нестационарных уравнений Навье-Стокса для несжимаемой жидкости, который базируется на методе конечных объемов. Уравнения решаются в неподвижной декартовой системе координат. Используются разностные схемы, имеющие второй порядок точности по пространству и времени. Для дискретизации конвективных членов используется TVD схема Chakravarthy-Osher. Численный алгоритм разработан для трехмерной неструктурированной сетки. Полученная система нелинейных алгебраических уравнений линеаризуется и решается по алгоритму PISO. Для решения соответствующих систем линейных алгебраических уравнений применяются итерационные солверы, использующие методы сопряженных/бисопряженных градиентов с предобусловливанием, которые можно найти в свободно распространяемых библиотеках, доступных в репозитарии NetLib [62]. Описанный численный алгоритм использован для прямого численного моделирования течения жидкости в канале с двумя последовательно расположенными стенозами при различной ширине канала в области между стенозами. Установлено, что при достаточно малой ширине межстенозной части канала в межстенозной области генерируются симметрично расположенные вихри. А при достаточно большой ширине в определенном диапазоне чисел Рейнольдса два ряда вихрей, генерируемых течением в межстенозной области, организуются в шахматном порядке. Oписано алгоритм чисельного розв'язання повної системи нестаціонарних рівнянь Навьє-Стокса для нестисливої рідини, що грунтується на методі скінчених об'ємів. Рівняння розв'язуються в нерухомій декартовій системі координат. Використовуються різничні схеми, що мають другий порядок точності у просторі і за часом. Для дискретизації конвективних членів використовується схема Chakravarthy-Osher. Чисельний алгоритм розроблено для тривимірної неструктурованої сітки. Отримана система лінійних алгебраїчних рівнянь лінеарізується і розв'язується за алгоритмом PISO. Для розв'язання відповідних систем лінійних алгебраїчних рівнянь застосовуються ітераційні солвери, що базуються на методах спряжених/біспряжених градієнтів з передобумовленням, які можна знайти у відкритих бібліотеках, доступних в репозитарії NetLib [62]. Описаний чисельний алгоритм використано для прямого чисельного моделювання течії рідини в каналі з двома послідовно розташованими стенозами за різної ширини канала в області між стенозами. Встановлено, що за досить малої ширини міжстенозної частини канала в міжстенозній області генеруються симетрично розташовані вихори. А за досить великої ширини у певному діапазоні чисел Рейнольдса два ряди вихорів, що генеруються течією в міжстенозній області, організуються в шаховому порядку. A numerical solution algorithm for the unsteady incompressible Navier-Stokes equations is described in the paper. The discretization procedure is based on the finite volume method. The equations are solved in the fixed Cartesian coordinate system. Stabilized and bounded differencing schemes, second-order accurate in both space and time, are used. The TVD differencing scheme of Chakravarthy-Osher is employed for the discretization of the convective terms. The numerical algorithm is designed for three-dimensional arbitrarily unstructured meshes. The system of non-linear algebraic equations obtained from the discretization procedure is linearized and solved with the PISO algorithm. The corresponding systems of linear algebraic equations are solved using the iterative solvers based on the preconditioned conjugate/biconjugate gradient methods. These solvers are available from the repository NetLib [62]. The numerical algorithm described in the paper is employed for the direct numerical simulation of the flow in a duct containing two sequential stenoses. The calculations are performed for various widths of the inter-stenoses segment of the duct. It is established that if the duct is sufficiently narrow between the stenoses, in some range of Reynolds number, the vortices generated in the inter-stenoses domain are situated symmetrically. But if the duct is sufficiently wide between the stenoses, the two rows of the vortices generated in the inter-stenoses domain are chequer-wise situated.
issn	1561-9087
url	https://nasplib.isofts.kiev.ua/handle/123456789/87748
citation_txt	Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения / В.С. Малюга // Прикладна гідромеханіка. — 2010. — Т. 12, № 4. — С. 45-62. — Бібліогр.: 62 назв. — рос.
work_keys_str_mv	AT malûgavs čislennoeissledovanietečeniâvkanalesdvumâposledovatelʹnoraspoložennymistenozamialgoritmrešeniâ AT malûgavs numericalinvestigationoftheflowinaductwithtwoserialstenosesalgorithmofthesolution
first_indexed	2025-11-25T06:10:05Z
last_indexed	2025-11-25T06:10:05Z
_version_	1850509225227190272
fulltext	ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 УДК 532.516 ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ТЕЧЕНИЯ В КАНАЛЕ С ДВУМЯ ПОСЛЕДОВАТЕЛЬНО РАСПОЛОЖЕННЫМИ СТЕНОЗАМИ. АЛГОРИТМ РЕШЕНИЯ В. С. М АЛ Ю Г А Институт гидромеханики НАН Украины, Киев Получено 23.12.2009 Oписан алгоритм численного решения полной системы нестационарных уравнений Навье-Стокса для несжимаемой жидкости, который базируется на методе конечных объемов. Уравнения решаются в неподвижной декартовой си- стеме координат. Используются разностные схемы, имеющие второй порядок точности по пространству и времени. Для дискретизации конвективных членов используется TVD схема Chakravarthy-Osher. Численный алгоритм ра- зработан для трехмерной неструктурированной сетки. Полученная система нелинейных алгебраических уравнений линеаризуется и решается по алгоритму PISO. Для решения соответствующих систем линейных алгебраических уравнений применяются итерационные солверы, использующие методы сопряженных/бисопряженных градиентов с предобусловливанием, которые можно найти в свободно распространяемых библиотеках, доступных в репозитарии NetLib [62]. Описанный численный алгоритм использован для прямого численного моделирования течения жидко- сти в канале с двумя последовательно расположенными стенозами при различной ширине канала в области между стенозами. Установлено, что при достаточно малой ширине межстенозной части канала в межстенозной области генерируются симметрично расположенные вихри. А при достаточно большой ширине в определенном диапазоне чисел Рейнольдса два ряда вихрей, генерируемых течением в межстенозной области, организуются в шахматном порядке. Oписано алгоритм чисельного розв’язання повної системи нестацiонарних рiвнянь Навьє-Стокса для нестисливої рiдини, що грунтується на методi скiнчених об’ємiв. Рiвняння розв’язуються в нерухомiй декартовiй системi коорди- нат. Використовуються рiзничнi схеми, що мають другий порядок точностi у просторi i за часом. Для дискретизацiї конвективних членiв використовується схема Chakravarthy-Osher. Чисельний алгоритм розроблено для тривимiрної неструктурованої сiтки. Отримана система лiнiйних алгебраїчних рiвнянь лiнеарiзується i розв’язується за алгори- тмом PISO. Для розв’язання вiдповiдних систем лiнiйних алгебраїчних рiвнянь застосовуються iтерацiйнi солвери, що базуються на методах спряжених/бiспряжених градiєнтiв з передобумовленням, якi можна знайти у вiдкритих бiблiотеках, доступних в репозитарiї NetLib [62]. Описаний чисельний алгоритм використано для прямого чисель- ного моделювання течiї рiдини в каналi з двома послiдовно розташованими стенозами за рiзної ширини канала в областi мiж стенозами. Встановлено, що за досить малої ширини мiжстенозної частини канала в мiжстенознiй обла- стi генеруються симетрично розташованi вихори. А за досить великої ширини у певному дiапазонi чисел Рейнольдса два ряди вихорiв, що генеруються течiєю в мiжстенознiй областi, органiзуються в шаховому порядку. A numerical solution algorithm for the unsteady incompressible Navier-Stokes equations is described in the paper. The discretization procedure is based on the finite volume method. The equations are solved in the fixed Cartesian coordinate system. Stabilized and bounded differencing schemes, second-order accurate in both space and time, are used. The TVD differencing scheme of Chakravarthy-Osher is employed for the discretization of the convective terms. The numerical algorithm is designed for three-dimensional arbitrarily unstructured meshes. The system of non-linear algebraic equations obtained from the discretization procedure is linearized and solved with the PISO algorithm. The corresponding systems of linear algebraic equations are solved using the iterative solvers based on the preconditioned conjugate/biconjugate gradient methods. These solvers are available from the repository NetLib [62]. The numerical algorithm described in the paper is employed for the direct numerical simulation of the flow in a duct containing two sequential stenoses. The calculations are performed for various widths of the inter-stenoses segment of the duct. It is established that if the duct is sufficiently narrow between the stenoses, in some range of Reynolds number, the vortices generated in the inter-stenoses domain are situated symmetrically. But if the duct is sufficiently wide between the stenoses, the two rows of the vortices generated in the inter-stenoses domain are chequer-wise situated. ВВЕДЕНИЕ Течения вязкой жидкости в канале с последо- вательно расположенными сужениями поперечно- го сечения (стенозами) широко встречаются как в природе, так и во многих технических устрой- ствах. Наличие таких препятствий в канале опре- деляет кинематику потока, а также характер и уровень акустических шумов, генерируемых пото- ком. Поэтому разработка методов расчета таких течений представляет практический интерес. Сло- жности расчета таких течений связаны с тем, что на различных участках такое течение обладает ра- зличной кинематической природой. Так, в опреде- ленной области его нужно рассматривать как об- текание препятствия. В частности, перед первым сужением образуется так называемая область под- пора или застойная зона. А далее, в области ме- жду первым и вторым сужениями такое течение напоминает детально изученное течение над пря- моугольным углублением. В частности, при доста- точно больших числах Рейнольдса может иметь место отрыв пограничного слоя на первом стенозе c© В. С. Малюга, 2010 45 ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 и его снос в область между стенозами. В данной статье мы детально описываем численный метод, использованный нами при расчете таких течений. Этот метод базируется на хорошо известном ме- тоде конечных объемов, однако ряд принципиаль- ных моментов должен быть обсужден отдельно. Это, в первую очередь, касается построения кон- вективной разностной схемы и выбора метода ре- шения получаемой нелинейной системы алгебраи- ческих уравнений. Многие процессы, наблюдаемые в природе и тех- нике, которые изучает механика жидкости, хара- ктеризуются высокими числами Рейнольдса, что подразумевает доминирование конвективных эф- фектов. В то время как базовые понятия конечно- объемной дискретизации хорошо известны [1 – 3], дискретизация конвективных членов продолжает оставаться предметом интенсивных дискуссий. В рамках второго порядка точности приемлемой схе- мой дискретизации могла бы быть центрально- разностная схема. Однако, как известно, в многих случаях применение центрально-разностной схе- мы для дискретизации конвективных членов при- водит к неустойчивости и появлению нефизиче- ских осцилляций решения. Дело в том, что, со- гласно теореме Годунова [4] (см. также учебную литературу, например, [2]), все линейные монотон- ные схемы для уравнения переноса имеют толь- ко первый порядок, т. е. линейные схемы второ- го и более высоких порядков не обладают свой- ством монотонности. Для достижения устойчиво- сти дискретизационной процедуры были предло- жены схемы первого порядка [5 – 7]. Однако такие схемы имеют неудовлетворительный порядок то- чности, что привело к необходимости разработки устойчивой схемы второго порядка. Такая схема была предложена Lax и Wendroff [8]. Позже появи- лось большое семейство схем Lax-Wendroff, в кото- рых устойчивость достигается за счет комбиниро- вания пространственной и временной дискретиза- ции, что привело к множеству двухшаговых [9, 10] и неявных схем [11, 12]. Однако в случае стацио- нарных задач комбинированная пространственно- временная дискретизация вводит нереалистичную зависимость решения от величины шага по вре- мени. Для разрешения этой проблемы в работах [13–15] было разработано семейство схем второ- го порядка с независимым интегрированием по времени. Хотя такой подход исключает зависи- мость пространственной точности от размера вре- менного шага, разностные схемы семейства Beam- Warming могут вызывать нефизичные осцилля- ции в решении, что приводит к потере ограни- ченности решения и существенно снижает при- менимость схемы. Предпринималась также по- пытка решить проблему ограниченности решения посредством введения члена искусственной диф- фузии [2]. Однако и такой подход не гарантирует ограниченности решения, а искусственная диффу- зия снижает точность схемы, особенно в областях высоких градиентов. Единственная конвективная схема, которая га- рантирует ограниченность и монотонность - это встречно-поточная схема [1]. Это происходит бла- годаря введению чрезмерной численной диффу- зии, что изменяет природу задачи. Задача, в кото- рой доминировали конвективные эффекты, стано- вится задачей, в которой конвективные и диффу- зионные эффекты сбалансированы. Такая потеря точности является неприемлемой и, следователь- но, было предложено несколько решений данной проблемы. Не претендуя на исчерпывающую об- щность данного обзора, выделим следующие кате- гории: • схемы, использующие точное или прибли- женное одномерное решение уравнения конвекции-диффузии для определения зна- чения неизвестных функций. К ним можно отнести схему LOADS (locally analytic di- fferencing scheme), предложенную в [16], а также схему со степенным законом [1]; • схемы, основанные на разностях против пото- ка. Например, противопоточная взвешенная схема [17], имеющая первый порядок, а так- же противопоточные схемы более высоких по- рядков: линейная противопоточная схема [18] и схема QUICK, предложенная Leonard [19]. Следует также упомянуть противопоточные схемы для сложных физических областей с сетками со скошенными ячейками или так называемые skew-upwind differencing schemes [20, 21]. В [22] представлен более детальный обзор, посвященный противопоточным схе- мам; • гибридные (или комбинированные) схемы [23]. Суть гибридной схемы состоит в том, что в области, где сеточное число Пекле меньше 2, используются центральные разности, в про- тивном случае – разности против потока. Как указано в [25], такой подход был впервые пре- дложен Allen и Southwell [24]; • смешанные схемы (blended differencing) [26]. Считая, что достаточное условие ограничен- ности решения является слишком жестким, Perić предложил подход, основанный на "сме- шении" схем, т. е. противопоточные разности 46 В. С. Малюга ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 используются в комбинации со схемами более высокого порядка (центральные или линей- ные противопоточные разности) таким обра- зом, чтобы была достигнута ограниченность решения. Все перечисленные выше схемы либо сильно жертвуют точностью, либо не гарантируют моно- тонности решения. Дальнейший поиск ограничен- ных и достаточно точных разностных схем при- вел к разработке концепции ограничения потока. Boris и Book [27] ввели понятие ограничителя по- тока в своей FCD (Flux Corrected Transport) ра- зностной схеме, которая в дальнейшем была обоб- щена в [28] на многомерные задачи. В дальнейшем эта идея была использована van Leer в серии ра- бот "Towards the ultimate conservative differencing scheme"[29–33]. В конечном итоге эти методы при- вели к появлению класса TVD схем (Total Variati- on Diminishing). TVD схемы развивали Harten [34, 35], Roe [36], Chakravarthy и Osher [37] и другие. Общая процедура построения TVD схемы описа- на Osher и Chakravarthy в [38], а Sweby [39] пред- ставил графическую интерпретацию ограничите- лей (диаграма Sweby) и исследовал точность ме- тода. Основная идея TVD метода состоит в том, что вклад схемы более высокого порядка и мо- нотонной схемы первого порядка зависят от ло- кальной формы самого решения, что делает схему нелинейной. Общие принципы построения ограни- ченных конвективных схем высокого порядка рас- смотрены в [40]. Эта работа содержит детальную классификацию и анализ большинства нелиней- ных скалярных конвективных схем, разработан- ных на данный момент. Анализ включает обзор и сравнение двух наиболее часто используемых под- ходов: ограничители потоков и нормализованные переменные (normalized variables). Согласно классификации [40], TVD схема, используемая в данной статье, соответствует обоб- щенной кусочно-линейной схеме Chakravarthy- Osher [37]. Такая схема вносит минимальную чи- сленную диффузию, что может быть существенно именно в рассматриваемой в данной статье задаче. Кроме рассмотренной выше проблемы построе- ния неосциллирующих и достаточно точных кон- вективных схем, имеют место еще две проблемы, с которыми сталкиваются исследователи при дис- кретизации уравнений Навье-Стокса. Это нели- нейность уравнения сохранения импульса и со- гласование полей скорости и давления. Нелиней- ный член уравнения имеет вид ∇ · UU и, сле- довательно, дискретная форма этого выражения будет квадратичной относительно скорости. Т. е. получаемая в результате система алгебраических уравнений будет нелинейной. Есть два возмож- ных пути решения проблемы: использовать сол- веры для нелинейных систем уравнений или лине- аризовать конвективный член. Принимая во вни- мание сложность и затратность нелинейных сол- веров, во многих случаях предпочтение отдается именно процедуре линеаризации. В этом случае конвективный член трактуется как перенос скоро- сти потоком. При этом значения потока берутся с временным лагом, т. е. уже известные, рассчитан- ные на предыдущем шаге по времени. Затем для уточнения значений скорости используется итера- ционная процедура, т. е. значения скорости, полу- ченные после решения соответствующей линейной системы, используются для расчета нового, более точного поля потока. И так число итераций повто- ряют необходимое число раз. Одна из основных проблем при численном реше- нии системы уравнений Навье-Стокса для несжи- маемой жидкости состоит в слабом согласовании полей скорости и давления. Так, при использова- нии итерационных процедур решения, линейные алгебраические уравнения, получаемые из уравне- ния сохранения импульса, используются для вы- числения дискретного набора значений скорости. Тогда естественно было бы использовать урав- нение неразрывности для вычисления давления. Однако в рамках модели несжимаемой жидкости давление или плотность не входят в уравнение неразрывности. Таким образом, если используе- тся итерационная процедура решения, необходимо найти способ как ввести давление в уравнение не- разрывности. Такие методы можно разделить на два класса. Один из них - это метод искусственной сжимаемости, в рамках которого течению припи- сывается слабая (но конечная) сжимаемость [41] и затем решается задача на установление. Другой класс методов - это методы, основанные на вычис- лении давления (pressure-based methods, PBM). PBM метод был введен Harlow и Welch [42] для расчета неустановившегося течения. Позже метод был распространен и на случай установившегося течения [1, 43]. Основная идея данного класса ме- тодов состоит в том, чтобы сформулировать урав- нение Пуассона для коррекции давления и затем вычислять новые значения полей давления и ско- рости до тех пор, пока не будет получено поле ско- рости, удовлетворяющее условию несжимаемости. В этот класс методов входят MAC (marker and cell) метод [42], алгоритм SIMPLE (semi-implicit method for pressure linked equations) [43], вклю- чая его усовершенствования – SIMPLER [1] и SI- MPLEC [44], а также алгоритм PISO (pressure В. С. Малюга 47 ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 implicit with splitting of operators), предложенный Issa в [45], который отличается двухшаговым кор- ректором. Как было показано в [45, 46], алгоритм PISO представляет собой устойчивый, хорошо схо- дящийся алгоритм, который вo многих задачах требует меньше расчетного времени, чем все упо- мянутые выше разновидности алгоритма SIMPLE. Особенно его преимущества видны при расчетах нестационарных течений. В силу вышесказанно- го именно алгоритм PISO использовался в данной статье. Следует также отметить, что PBM метод был имплементирован на двух различных типах ор- ганизации сетки: разнесенные сетки (с различ- ными контрольными объемами для скорости и давления) [47] и совмещенные сетки (с теми же контрольными объемами для всех искомых по- лей) [48]. Разнесенные сетки часто использовались на более ранних стадиях развития расчетной ги- дромеханики. Их имплементация не представляет сложности, если сетки структурированы. Однако с тех пор, как начали активно развиваться алгорит- мы для неструктурированных сеток, совмещенные сетки стали использоваться все чаще и чаще. В на- стоящее время все современные численные пакеты используют совмещенные сетки. В данной статье, поскольку мы строим численный алгоритм для не- структурированной сетки, все поля (и скорость, и давление) приписываются в центроидах контроль- ных объемов. Цель данной работы состоит в применении по- строенного численного алгоритма к решению за- дачи о движении несжимаемой жидкости в пло- ском канале с последовательно расположенными двусторонними внезапными сужениями его попе- речного сечения (стенозами), а также анализе ки- нематических особенностей такого потока. В дан- ной статье мы рассматриваем ламинарные режи- мы течения и показываем, как качественно изме- няется кинематика течения при уменьшении глу- бины полостей, образованных последовательными сужениями. 1. ПОСТАНОВКА ЗАДАЧИ В данной работе рассматривается задача о те- чении жидкости в плоском канале при наличии следующих друг за другом двух стенозов. Также предполагается, что в межстенозном пространс- тве может иметь место некоторое сужение канала. Расчетная область и принятые обозначения пред- ставлены на рис. 1. Предполагается, что стенки канала, а также стенки стенозов неподвижны и абсолютно жесткие, а поток жидкости попадает в расчетную область через левую границу x = 0; 0 ≤ y ≤ D1 и покидает расчетную область на правой границе x = L1; 0 ≤ y ≤ D1. Вверх по потоку от первого стеноза (0 ≤ x ≤ L2), а также вниз по потоку от второго стеноза (L2 + L3 + L4 + L5 ≤ x ≤ L1) канал имеет ширину D1, в то время как в межстенозном пространстве (L2 +L3 ≤ x ≤ L2 +L3 +L4) ширина канала равна D4 ≤ D1. Т.е. в межстенозном пространстве может иметь место некоторое сужение канала. Задача решается в рамках модели вязкой несжи- маемой ньютоновской жидкости с плотностью ρ и, следовательно, основным параметром задачи слу- жит число Рейнольдса. Поскольку в задаче при- сутствует не один, а несколько заданных размеров, любой из них может быть в принципе принят за масштаб длины и, следовательно, число Рейнольд- са может быть определено различным образом. Остановимся более детально на выборе масштаба длины. Объектом нашего исследования являются особенности движения среды во всей расчетной области, однако наибольший интерес представля- ет течение в межстенозной области и, прежде все- го, на границах межстенозных полостей (L2+L3 ≤ x ≤ L2 + L3 + L4, D1/2− D4/2 ≤ y ≤ D1/2− D2/2 и L2 + L3 ≤ x ≤ L2 + L3 + L4, D1/2 + D2/2 ≤ y ≤ D1/2 +D4/2) и струи (L2 +L3 ≤ x ≤ L2 +L3 +L4, D1/2 − D2/2 ≤ y ≤ D1/2 + D2/2), где (как было показано в [49]) образуются два сдвиговых слоя, которые формируются как пристенные сдви- говые слои на поверхностях верхней и нижней ча- стей первого стеноза, а затем отрываются и снося- тся течением в область между стенозами, где мо- гут уже рассматриваться как свободные сдвиговые слои. В связи с этим, наиболее естественным было бы определить число Рейнольдса по длине межсте- нозных полостей, вдоль которых происходит ра- звитие сдвиговых слоев: Re = V2L4/ν , где ν – кине- матическая вязкость среды. При этом масштабом длины является расстояние между стенозами L4, масштабом скорости – скорость V2, т.е. скорость потока в отверстии первого стеноза, осредненная по вертикальному сечению, масштабом времени – величина L4/V2, а масштабом давления – удвоен- ный скоростной напор в отверстии стеноза ρV 2 2 . В рамках принятой модели процесс описыва- ется нестационарной системой уравнений Навье- Стокса. В безразмерных физических переменных их можно представить в тензорной форме следу- ющим образом: ∂U ∂t + ∇ · UU = 1 Re ∇ · ∇U−∇p , (1) 48 В. С. Малюга ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 V 1 x y D 1 L 1 D 2 D 3 L 2 L 3 L 4 L 5 V 2 0 D 4 D 1 Рис. 1. Геометрия расчетной области ∇ · U = 0 , (2) где p – скалярное поле давления; U – векторное поле скорости; UU – тензор второго ранга, опре- деленный как внешнее произведение векторов. По- компонентное представление этих уравнений мож- но найти в учебной литературе (см., например, [53]). На границе расчетной области задавались сле- дующие граничные условия для скорости: равно- мерный поток на входе в расчетную область при x = 0, условие прилипания на твердых поверхно- стях канала и стенозов, а также «мягкие» грани- чные условия типа линейной экстраполяции, со- ответствующие равенству нулю нормального гра- диента скорости, на выходе из расчетной области при x = L1. Для давления условие равенства ну- лю нормального градиента формулировалось по всей границе расчетной области [2] за исключени- ем выхода из расчетной области. На выходе зада- валось постоянное давление. Хорошо известно, что двумерная постановка за- дачи позволяет сформулировать ее в переменных завихренность – функция тока (см., например, [54]). Такой подход позволяет исключить давле- ние из уравнений движения и таким образом избе- жать сложностей, связанных с определением поля давления и согласованием полей. Однако распро- странение подхода с использованием завихренно- сти и функции тока на трехмерные задачи осло- жнено тем, что для трехмерного течения нельзя ввести скалярную функцию тока. Использование же векторного потенциала приводит к необходи- мости решать на каждом временном слое три па- раболических и три эллиптических уравнения, что существенно усложняет задачу [25]. По этой при- чине в данной работе, так же как в работе [49], за- дача решалась именно в физических переменных скорость – давление. Значения геометрических параметров канала выбирались такие же, как в экспериментальной работе [50], а именно: D1 = 18 мм, D2 = D3 = 2 мм, L1 = 150 мм, L2 = 20 мм, L3 = L5 =1 мм, L4 = 5.8 мм. Кинематическая вязкость прини- малась равной вязкости воздуха при температуре 20◦С: ν = 1.5 · 10−5 м2/с. Определенные выше размеры области при расчетах не изменялись, а ширина канала в межстенозной области D4 варьи- ровалась в пределах 6 – 18 мм. 2. ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ЗАДАЧИ Для численного решения сформулированной за- дачи используется метод конечных объемов. Хо- тя основные положения конечно-объемной дискре- тизации хорошо известны и изложены в учебной литературе (см. [1 – 3]), проблемы, связанные с дискретизацией конвективных членов, согласова- нием полей скорости и давления, а также с приме- нением неструктурированных и неортогональных сеток, по-прежнему являются предметом интен- сивных дискуссий. В данной работе, основываясь на идеях, изложенных в [53], мы строим конечно- объемную схему второго порядка, которая бази- руется на вычислении потока среды сквозь грани конечного объема и применима к неструктуриро- ванным сеткам, которые могут быть в принципе и неортогональными. В данной работе в силу прямо- угольной геометрии рассматриваемой области мы используем ортогональную сетку с прямоугольной формой ячейки, которую, однако, все равно рас- сматриваем как неструктурированную сетку. В. С. Малюга 49 ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 m n f.. .M N Рис. 2. Контрольные объемы m и n, имеющие общую грань f Численная схема строится для общего случая трехмерного течения. А двумерные задачи рас- сматриваются как частный случай, когда в на- правлении z имеется только один контрольный объем и компонента скорости Uz = 0. Таким обра- зом, расчетная область разбивается на контроль- ные объемы, имеющие форму прямоугольного па- раллелепипеда. Два типичных соседних контроль- ных объема с номерами m и n показаны на рис. 2. Грань, разделяющую эти два объема, обозначено через f . Точки M и N являются центроидами кон- трольных объемов m и n соответственно, т. е. ∫ Vm (x − xM )dV = ∫ Vn (x − xN )dV = 0 . (3) Интегрирование в (3) происходит по соответству- ющему контрольному объему. Для того, чтобы численная схема имела второй порядок точности, неизвестная функция должна иметь линейную вариацию в контрольном объеме: U(x) = UM + (x− xM ) · (∇U)M , (4) где индекс M означает значение функции в точке M , т. е. UM = U(xM ). Тогда в силу (3) имеем ∫ Vm UdV = UMVm . (5) Основная идея используемого метода конечных объемов состоит в том, что уравнения движения (1), (2) интегрируются по контрольному объему m: ∂ ∂t ∫ Vm UdV + ∫ Vm ∇ ·UUdV = = 1 Re ∫ Vm ∇ · ∇UdV − ∫ Vm ∇p dV , (6) ∫ Vm ∇ · UdV = 0 . (7) А для вычисления интегралов в (6), (7) использу- ется обобщенная форма теоремы Остроградского- Гаусса, а именно, три следующих тождества: ∫ Vm ∇ ·A dV = ∫ Sm dS · A , (8) ∫ Vm ∇A dV = ∫ Sm dSA , (9) ∫ Vm ∇a dV = ∫ Sm dS a , (10) где Sm – поверхность контрольного объема m. Вычисляя интеграл по контрольному объему от дивергенции скорости, из (8) получим дис- кретную форму второго порядка теоремы Гаусса. Учитывая (3), имеем (∇ · U)MVm = ∑ f Sf ·Uf . (11) Здесь Sf обозначает вектор внешней нормали к грани f , модуль которого равен площади грани Sf . (На практике, при составлении компьютерной программы, для определенности принимают, что вектор нормали к грани f направлен от контроль- ного объема с меньшим номером в сторону кон- трольного объема с большим номером. Если таким образом определенный вектор является внешним к ячейке m, то при суммировании в (11) берется знак плюс, в противном случае - минус.) Сумми- рование в правой части проводится по всем гра- ням контрольного объема m. Вычисляя интеграл по контрольному объему от градиента давления, из (10) получаем следующее дискретное соотно- шение второго порядка: (∇p)MVm = ∑ f Sfpf . (12) Дискретизация уравнений движения (1), (2) осно- вывается на дискретном аналоге второго порядка теоремы Гаусса (11) и (12). Дискретизация конвективных членов приводит к соотношениям ∫ Vm ∇ · UUdV = ∑ f Sf · UfUf = ∑ f FUf , (13) где F представляет собой поток среды (flux) сквозь грань f , т.е. F = Sf · Uf . Вычисление это- го потока F будет детально обсуждаться ниже. В формулу (13) входит значение скорости на грани 50 В. С. Малюга ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 f , которое определяется из значений скорости в центроидах соседних ячеек M и N посредством конвективных разностных схем. В рамках используемого метода в качестве схе- мы дискретизации для конвективного члена мож- но было бы предложить центрально-разностную схему второго порядка: Uf = CUM + (1 − C)UN , (14) где постоянный коэффициент определен как отно- шение расстояний C = \|fN \|/\|MN \|. В [53] показа- но, что даже для неоднородной сетки эта схема имеет второй порядок точности. Однако при определенных условиях, например, если течение при больших числах Рейнольдса на- бегает на выступающий угол, центральные разно- сти могут приводить к неустойчивой схеме дис- кретизации и таким образом вызывать нефизи- чные осцилляции, что существенно снижает ка- чество численного решения или даже делает его неприемлемым вообще. Единственная конвектив- ная разностная схема, которая является безуслов- но устойчивой и гарантирует ограниченность ре- шения, – встречно-поточная разностная схема [1]: Uf = { UM , for F ≥ 0 , UN , for F < 0 . (15) Однако эта схема имеет лишь первый порядок то- чности, а устойчивость схемы и ограниченность решения достигается за счет внесения значитель- ной доли численной диффузии. И хотя измельче- ние сетки улучшает точность решения, требуемое число контрольных объемов делает эту схему не- приемлемой для практики. Было предложено несколько методов решения этой проблемы. В настоящее время активно ра- звиваются схемы с ограничителями потока, кото- рые привели к выделению класса схем TVD (Total Variation Diminishing) [40]. Основная идея состоит в том, что схема дискретизации зависит от локаль- ной формы решения. Она предоставляет достато- чную точность (более чем первого порядка) и в то же время гарантирует ограниченность реше- ния (подавляет нефизичные осцилляции, которые генерируют классические схемы второго порядка, такие как центрально разностная схема). Любая TVD схема может быть представлена в виде суммы Uf = (Uf) UD + Ψ [ (Uf) CD − (Uf) UD ] , (16) где (Uf) UD – значение скорости на грани кон- трольного объема, рассчитанное по противопото- чной схеме первого порядка (15); (Uf) CD – значе- ние скорости на грани контрольного объема, рас- считанное по центральноразностной схеме второго порядка (14); Ψ – нелинейный ограничитель по- тока, который является функцией значений ско- рости в окрестности грани контрольного объема. Следуя [30, 37, 39], предполагается, что ограни- читель потока оказывается функцией отношения градиентов на границах между данной ячейкой и ее соседями вверх и вниз по потоку. В данной рабо- те мы используем TVD схему, имплементирован- ную в [55] (limitedLinear), для которой ограничи- тель потока определен как Ψ(r) = max ( 0, min ( 2 k r, 1 )) , (17) где r = (UC − UU)/(UD − UC). Здесь U, C, D – три последовательно расположенные точки (U – вверх по потоку от C, D – вниз по потоку от C). Согласно классификации [40] такой ограничи- тель соответствует обобщенной кусочно-линейной схеме Chakravarthy-Osher [37]. Коэффициент k за- дается в диапазоне 0 < k ≤ 1. Его нулевое значе- ние k → 0 соответствует более высокой точности, а значение k = 1 – более высокой устойчивости счета. Обычно принимается значение k = 1. Диффузионный член дискретизируется анало- гичным образом: ∫ Vm ∇ · ∇UdV = ∑ f Sf · (∇U)f . (18) Учитывая линейную вариацию U и ортогональ- ность сетки, правую часть (18) можно представить следующим образом: Sf · (∇U)f = Sf UN −UM \|MN \| . (19) Для дискретизации производной по времени используется неявная трехточечная несимметри- чная схема с разностями назад (backward di- fferencing), которая имеет второй порядок точно- сти: ∫ Vm ∂U ∂t dV = ( ∂U ∂t ) M Vm = = 3 2 U k M − 2Uk−1 M + 1 2 U k−2 M ∆t Vm , (20) где U k M = UM (k ∆t) – значения скорости на k-том временном слое, а значения на предыдущих вре- менных слоях U k−1 M и U k−2 M известны. Такая схема В. С. Малюга 51 ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 применима в том случае, если сетка остается непо- движной и, следовательно, положение центроида M не изменяется. Итак, уравнение (1) в дискретной форме можно записать следующим образом: 3 2 U k M − 2Uk−1 M + 1 2 U k−2 M ∆t Vm + ∑ f FU k f − − 1 Re ∑ f Sf · (∇U k)f = −(∇pk)MVm = = − ∑ f Sfpk f , (21) а уравнение (2) – в виде: ∑ f Sf · Uk f = 0. (22) Рассмотрим дискретную имплементацию грани- чных условий. Условие Дирихле: задано значение скорости на границе. В этом случае конвективный член также дискретизируется по формуле (13). Значения скорости Uf и потока F на грани f , ле- жащей на границе области, заданы. Диффузион- ный член задается формулой (18), а значение нор- мального градиента скорости на границе рассчи- тывается из известных значений скорости в цен- троиде ячейки и на грани, лежащей на границе области: Sf · (∇U)f = Sf Uf − UM \|Mf \| . (23) Условие Неймана: задан нормальный градиент скорости на границе. В этом случае значение ско- рости на границе Uf , которое необходимо для вычисления правой части (13), рассчитывается из значения скорости в центроиде контрольного объема и заданного значения нормального гради- ента на грани: Uf = UM + \|Mf \| Sf Sf · (∇U)f . (24) Что же касается диффузионного члена, то нор- мальный градиент скорости на границе Sf · (∇U)f в правой части (18) задается непосредственно гра- ничным условием. В принципе, интерполируя в (21), (22) значения искомых полей на гранях ячеек через их значе- ния в центроидах соседних ячеек, используя опи- санные выше схемы интерполяции, можно полу- чить систему алгебраических уравнений, которая и должна быть решена численно. Однако в дан- ном случае следует уделить особое внимание двум проблемам: нелинейности задачи и согласованию полей скорости и давления. Причем, первая про- блема следует непосредственно из уравнений дви- жения и не зависит от выбранного метода реше- ния. Действительно, в силу нелинейности конве- ктивного члена (∇ · UU) уравнения (1) получае- мая система алгебраических уравнений также бу- дет нелинейной (поток F является функцией ско- рости U). Поскольку нелинейные солверы требу- ют колоссальные компьютерные ресурсы, обычно прибегают к линеаризации конвективного члена. Процедура линеаризации состоит в том, что для расчета потока F используются известные зна- чения скорости на предыдущем временном слое, а затем полученные значения уточняются путем итераций. Проблема определения поля давления, а так- же согласования полей скорости и давления сле- дует из необходимости использовать в практиче- ских задачах итерационные, а не прямые мето- ды решения системы линейных алгебраических уравнений (СЛАУ). В то время как дискретное уравнение сохранения импульса (21) используе- тся для нахождения компонент скорости, урав- нение для нахождения давления отсутствует, так как в рамках модели несжимаемой жидкости дав- ление не входит в уравнение неразрывности (22). Следовательно, при использовании итерационных солверов СЛАУ необходимо найти способ ввести давление в уравнение неразрывности. В данной статье используется процедура PISO, предложен- ная Issa [45], которая является процедурой типа предиктор-корректор с двумя корректорами. Уравнение для давления выводится следующим образом. Полностью дискретный аналог уравне- ния движения (21) можно разрешить относитель- но U k M : U k M = H o + H k − A(∇pk)M . (25) Здесь член H k содержит скорости в центроидах соседних ячеек U k N , а член H o содержит известные значения поля скорости на предыдущих времен- ных слоях U k−1 M и U k−2 M . Скорости на гранях ячейки получаются путем интерполяции на грань: U k f = H o f + H k f − Af(∇pk)f . (26) Умножая скалярно на вектор нормали Sf , полу- чаем выражение для потока через грань f : F = Sf ·Uk f = Sf · [ H o f + H k f − Af(∇pk)f ] . (27) 52 В. С. Малюга ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 Окончательно, дискретный аналог системы не- сжимаемых уравнений Навье-Стокса (21), (22) имеет вид: U k M = H o + H k − A Vm ∑ f Sfpk f , (28) ∑ f AfSf · (∇pk)f = ∑ f Sf · [ H o f + H k f ] . (29) Алгоритм PISO для нестационарных течений состоит из следующих шагов: 1. Предиктор. Для получения первого прибли- жения нового поля скорости U k решается дис- кретное уравнение сохранения импульса (28). Поле давления pk на данном этапе неизве- стно, поэтому используются известные зна- чения на предыдущем временном слое pk−1. Также выражение H(U) линеаризуется. Для этого вместо F подставляются его известные значения на предыдущем временном слое. Та- ким образом, уравнение принимает вид: U ∗ M = H o + H ∗ − A Vm ∑ f Sfpk−1 f , (30) где верхний индекс ∗ обозначает первое при- ближение. 2. Корректор 1. Используя полученные значе- ния U ∗, вычисляется оператор H ∗ и решается уравнение для давления (29), которое теперь принимает вид: ∑ f AfSf · (∇p∗)f = ∑ f Sf · [ H o f + H ∗ f ] . (31) Таким образом получаем первое приближе- ние поля давления p∗, которое используется для коррекции поля скорости. Следует отме- тить, что всякий раз как новое приближение для давления получено, следует также пере- считать поле потока F , который входит в ко- эффициенты оператора H. Для этого мож- но использовать уравнение (27), где в правую часть вместо H k, pk подставлены их после- дние полученные приближения. Коррекция поля скорости производится по явной схеме: U ∗∗ M = H o + H ∗ − A Vm ∑ f Sfp∗f . (32) 3. Корректор 2. Использование явной процеду- ры коррекции подразумевает, что поле скоро- сти корректируется за счет учета новых зна- чений F и p. Но при этом значения скоро- сти U ∗ N в центроидах соседних ячеек, кото- рые входят в H ∗, берутся неоткорректирован- ные, т. е. те, которые были получены в первом приближении. Поэтому необходимо вернуться к шагу 2 и повторить петлю PISO, т. е. прове- сти вычисление оператора H, решить уравне- ние для давления, вычислить поток F и про- вести явную коррекцию поля скорости. И та- ким образом получаем значения полей U ∗∗∗, p∗∗ после второго корректора. Часто используется именно описанная схема, со- стоящая из предиктора и двух корректоров. Одна- ко петлю PISO можно повторять и большее число раз или повторять до тех пор, пока точность не достигнет заданного значения. Рассмотрим теперь решение получаемых систем линейных алгебраических уравнений (30) и (31). Такая СЛАУ может быть записана в виде: AMxk M + ∑ N ANxk N = bM , (33) где xk M – неизвестные, а суммирование произво- дится по всем контрольным объемам n соседним с контрольным объемом m, т. е. имеющими об- щую грань с m. Система (33) может быть реше- на различными методами, которые можно разбить на две главные категории: прямые и итерацион- ные методы. Прямые методы приемлемы для не- больших СЛАУ с заполненными матрицами. Чис- ло операций, необходимых для получения реше- ния, возрастает как квадрат числа неизвестных, поэтому при работе с большими разреженными матрицами прямые методы становятся непозволи- тельно затратными. Итерационные методы в этом случае являются более экономными. В рамках ите- рационных методов решения мы стартуем с неко- торого начального приближения и затем посред- ством итерационной процедуры улучшаем реше- ние до тех пор, пока не достигаем некоторой напе- ред заданной точности. Также, в отличие от пря- мых методов, используемые итерационные методы сохраняют разреженность матрицы, что является важным их свойством, позволяющим существенно снизить требования к памяти компьютера. Следует отдельно оговорить вклад конвектив- ного члена. Известно, что диагональное домини- рование в матрице обеспечивает только встречно- поточная схема UD. Любые другие схемы созда- ют отрицательные коэффициенты, что может не- гативно сказываться на сходимости итерационно- го метода. Для того, чтобы улучшить качество матрицы для разностных схем высоких порядков, В. С. Малюга 53 ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 Khosla и Rubin [57] предложили метод отложенной коррекции (deferred correction implementation) для конвективного члена. Согласно этому методу лю- бая разностная схема рассматривается как некий апгрейд схемы UD. Та часть конвективного чле- на, которая соответствует UD, рассматривается неявно, т. е. встраивается в матрицу, а оставша- яся часть переносится в источниковый член, т. е. в правую часть системы. Таким образом, обеспе- чивается диагональное доминирование в матрице системы. Классические градиентные методы, такие как метод наискорейшего спуска и метод минималь- ных невязок, очень быстро минимизируют фун- кционалы на первых итерациях, а потом начи- нают “буксовать”, т. е. дальнейшее итерирова- ние показывает очень медленную сходимость, де- лающую применение градиентных методов неэф- фективным. Это особенно проявляется, когда соб- ственные значения матрицы А сильно различны. Этот недостаток эффективности градиентных ме- тодов устранен в методе сопряженных градиентов, первый вариант которого был предложен Hestens и Steifel [58] в 1952 году. Алгoритмы метода со- пряженных градиентов oтнoсятся к числу нaи- бoлee эффeктивных методов рeшeния СЛAУ бoль- шoй рaзмeрнoсти, вoзникaющих при числeнoм рeшeнии зaдaч мeхaники сплoшных срeд. Они ре- шают систему уравнений за конечное число ите- раций, нe прeвoсхoдящee числa нeизвeстных. При хорошем начальном приближении число итераций резко сокращается. Также к резкому сокращению числа итераций ведет предобусловливание. Для симметричных матриц в данной работе применя- ется солвер ICCG (Incomplete Cholesky preconditi- oned Conjugate Gradient), т. е. метод сопряженных градиентов с предобусловливанием типа неполной факторизации Холецкого. Метод детально описан Jacobs [59]. Для асимметричных матриц использу- ется солвер Bi-CGSTAB, представленный van der Vorst [60]. Используемые солверы реализованы в свободно распространяемых библиотеках, досту- пных через электронный сервис Netlib [62], а их подробное описание может быть найдено в [56, 61]. Представленный алгоритм был протестирован на хорошо известных задачах: задаче о течении в прямоугольной полости, вызванном движением одной из стенок (lid-driven cavity) [51], при различ- ных числах Рейнольдса, а также задаче об отрыв- ном течении за круговым цилиндром (вихревая дорожка Кармана) [52]. Тест показал хорошее сов- падение полученных результатов с результатами других авторов. 3. АНАЛИЗ РЕЗУЛЬТАТОВ РАСЧЕТОВ Согласно постановке задачи, в начальный мо- мент времени жидкость покоилась, а давление было постоянным во всей расчетной области. За- тем на входе в расчетную область задавался рав- номерный поток. Возникающая разность давлений на входе и на выходе приводила к формированию течения в расчетной области. Рассмотрим сначала этот переходной процесс в случае, когда ширина межстенозной части канала мала, а именно: за- дадим D4 = 6 мм. В [49] было показано, что на верхней и нижней поверхностях первого стеноза формируются пограничные слои, которые затем отрываются на передней кромке первого стеноза и сносятся потоком в межстенозную область. Таким образом, на границе струи и межстенозных поло- стей образуются два сдвиговых слоя, имеющие за- вихренность противоположного знака. Было по- казано, что при относительно малых значениях числа Рейнольдса (Re = 2088) течение является стационарным, а сдвиговые слои в межстенозной области устойчивы. Поэтому в данной статье пе- рейдем сразу к следующему значению числа Рей- нольдса, рассмотренному в [49]. Во всех приведенных ниже примерах скорость на входе в расчетную область задавалась V1 = 0.9 м/с, что соответствует значению числа Рей- нольдса Re = 3132. На рис. 3 представлено поле завихренности для пяти моментов времени. Цвето- вая шкала справа от рисунка показывает соответ- ствие цвета и значения завихренности1. Здесь кра- сный цвет соответствует положительной завихрен- ности, а синий - отрицательной. Как видно, в мо- мент времени t = 10−3 с (рис. 3, а) на жестких по- верхностях, включая поверхности стенозов, фор- мируются ламинарные пограничные слои. Причем завихренность имеет противоположные знаки на верхней и нижней жестких поверхностях. Затем сформировавшиеся пограничные слои отрываются от передних кромок как первого, так и второго сте- нозов, образуя две вихревые пары. Первая вихре- вая пара выносится потоком из отверстия перво- го стеноза в межстенозное пространство, а вторая пара - в канал за стенозами. В момент времени t = 2·10−3 с (рис. 3, б) первая пара попадает в меж- стенозные полости, где крупные вихри начинают формироваться у внутренней стенки второго сте- 1К сожалению, цветные рисунки представлены только в электронной версии статьи, размещенной в интернете на сайте Института гидромеханики НАН Украины. В печа- тном же варианте поле завихренности представлено в гра- дации серого, что не позволяет передать информацию о на- правлении завихренности. Автор приносит свои извинения и рекомендует обратитьсяк электронномуварианту статьи. 54 В. С. Малюга ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 0 0.5 -0.5 -1 1 X10 4 0 0.5 -0.5 -1 1 X10 4 а б X10 4 0 0.5 -0.5 -1 1 0 0.5 -0.5 -1 1 X10 4 в г 0 0.5 -0.5 -1 1 X10 4 0 0.5 -0.5 -1 1 X10 4 д е Рис. 3. Поле завихренности. Глубина каждой межстенозной полости равна высоте отверстий в стенозах, т. е. D2 = D3 = 2 мм, D4 = 6 мм: а – момент времени t = 10 −3 с, б – t = 2 · 10 −3 с, в – t = 4 · 10 −3 с, г – t = 5 · 10 −3 с, д – t = 9 · 10 −3 с, е – развитое течение при достаточно больших временах. Переходной процесс формирования течения завершен не только в межстенозной области, но и в канале за стенозами В. С. Малюга 55 ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 0 0.5 -0.5 -1 1 X 0 0.5 -0.5 -1 1 X10 а б 0 0.5 -0.5 -1 1 X10 0 0.5 -0.5 -1 1 X10 в г 0 0.5 -0.5 -1 1 X10 0 0.5 -0.5 -1 1 X10 д е Рис. 4. Поле завихренности. Глубина каждой межстенозной полости равна четырем высотам отверстий в стенозах, т. е. D2 = D3 = 2 мм, D4 = 18 мм: а – момент времени t = 3, 5 · 10 −3 с, б – t = 4, 2 · 10 −3 с, в – t = 7, 5 · 10 −3 с, г – t = 1, 3 · 10 −2 с, д – t = 2, 3 · 10 −2 с, е – t = 3 · 10 −2 с. Переходной процесс формирования течения завершен в межстенозной области 56 В. С. Малюга ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 ноза. В это же время за вторым стенозом начинает формироваться вторая вихревая пара. На следую- щем рисунке (рис. 3, в) представлена картина за- вихренности для момента времени t = 4 · 10−3 с. Видно, что в левой части межстенозных полостей, т. е. у поверхности первого стеноза, формируется еще одна пара крупных вихрей, причем завихрен- ность их будет иметь противоположный знак по отношению к завихренности крупных вихрей, ра- сположенных у поверхности второго стеноза. Так, если в верхней полости крупный вихрь, располо- женный вблизи второго стеноза, имеет положи- тельную завихренность, то вихрь, расположенный вблизи первого стеноза, будет иметь отрицатель- ную. На следующем рисунке (рис. 3, г) видно, что в момент t = 5 · 10−3 с процесс формирования вто- рого крупного вихря в каждой из межстенозных полостей уже завершен. Также на границе струи и межстенозных полостей формируются верхний и нижний сдвиговые слои, имеющие завихренность противоположных знаков. Причем эти сдвиговые слои имеют волнистую форму и возмущаются при- мерно посередине между первым и вторым стено- зом, т. е. между первыми и вторыми крупными вихрями, расположенными в межстенозных поло- стях. В этой области в верхнем и нижнем сдви- говых слоях формируется симметричная вихревая пара, которая затем сносится потоком в отверстие второго стеноза. Это отчетливо видно на рис. 3, д, на котором изображена картина завихренности в момент времени t = 9 · 10−3 с. Можно говорить, что к этому моменту времени переходные процес- сы формирования течения в межстенозной обла- сти завершены, а само течение является перио- дическим во времени. На рис. 3, е представлено поле завихренности для достаточно большого вре- мени, когда течение сформировалось во всей ра- счетной области, включая канал за вторым сте- нозом. В области за стенозами поле скорости те- чения не будет симметричным относительно оси канала y = D1/2, так как струя, истекающая из отверстия в плоский канал, совершает автоколе- бания, медленные по сравнению с автоколебани- ями в межстенозной области. Несмотря на это, в межстенозной области по прежнему доминирует симметричная часть поля скорости. Теперь рассмотрим такой же переходной про- цесс для случая широкой межстенозной части ка- нала, а именно: D4 = 18 мм, т. е. для случая, ко- гда ширина каждой межстенозной полости в че- тыре раза больше ширины отверстия в стенозе. Естественно, на входе задано то же значение ско- рости V1 = 0.9 м/с, что соответствует значению числа Рейнольдса Re = 3132. В начальные мо- менты времени течение подобно тому, что имело место в случае узкой межстенозной части кана- ла, описанному выше. На жестких поверхностях, включая поверхности стенозов, формируются ла- минарные пограничные слои. Причем завихрен- ность имеет противоположные знаки на верхней и нижней жестких поверхностях. Затем сформиро- вавшиеся пограничные слои отрываются от пере- дних кромок как первого, так и второго стенозов, образуя две вихревые пары. Первая вихревая па- ра выносится потоком из отверстия первого сте- ноза в межстенозное пространство, а вторая па- ра – в канал за стенозами. Поле завихренности в этом случае подобно тому, что представлено на рис. 3, а для области с узкой межстенозной частью канала. В дальнейшем вихревая пара, сформиро- ванная в межстенозной области, так же, как и в предыдущем примере, попадает в межстенозные полости. С этого момента течение в области с ши- рокой межстенозной частью канала существенно отличается от течения в области с узкой межсте- нозной частью канала, рассмотренного выше. Так вихри, попавшие в межстенозные полости, движу- тся вдоль поверхности второго стеноза, в результа- те чего в каждой из полостей образуется своя ви- хревая пара. На рис. 4, а видно, как эти вихревые пары, описав некоторые траектории в полостях, в момент времени t = 3, 5 · 10−3 с опять возвра- щаются к струе примерно посередине между пер- вым и вторым стенозом. Таким образом, две ви- хревые пары сверху и снизу ударяются в струю и возмущают ее. На рис. 4, б видно, что возмущения в верхнем и нижнем сдвиговых слоях, вызванные ударом по струе двух вихревых пар, будут симме- тричны относительно оси канала y = D1/2. Затем процесс повторяется. Однако на этот раз видно (рис. 4, в), что вихревые пары ударяют не в струю, а в поверхность первого стеноза. Это существен- но уменьшает их кинетическую энергию и таким образом на этот раз возмущения струи будут суще- ственно слабее, чем при первом ударе. На третий раз (рис. 4, г) вихревая пара также ударяет в по- верхность первого стеноза, причем вихри, распо- ложенные дальше от оси канала y = D1/2, уходят вглубь межстенозных полостей, формируя втори- чный вихрь в углу полости вблизи первого стено- за. Таким образом, кинетическая энергия вихре- вой пары расходуется на формирование течения во всей области межстенозных полостей. Удар по струе будет настолько слабым, что он практически не вносит никаких возмущений. На рис. 4, г видно, что течение в межстенозном участке струи можно считать стационарным, а сдвиговые слои невозму- щенными. К моменту t = 2, 1 · 10−2 с течение уже В. С. Малюга 57 ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 0 0.5 -0.5 -1 1 X10 4 Рис. 5. Поле завихренности. Глубина каждой межстенозной полости равна двум высотам отверстий в стенозах, т. е. D2 = D3 = 2 мм, D4 = 10 мм сформировано практически по всей области меж- стенозных полостей. На рис. 4, д видно как вихри, движущиеся в межстенозных полостях, в четвер- тый раз приближаются к струе. На этот раз видно, что в сдвиговых слоях на границе струи и полостей начинают формироваться возмущения, причем на этот раз верхний и нижний сдвиговые слои уже не будут симметричными относительно оси канала. Далее, с течением времени эти возмущения усили- ваются, и в момент времени t = 3 ·10−2 с (рис. 4, е) уже видно, как по мере приближения ко второму стенозу в сдвиговых слоях сворачиваются два ря- да вихрей, которые располагаются в шахматном порядке по отношению друг к другу. Рассмотрим промежуточный случай D4 = 10 мм. Опустим подробное описание переходного процесса формирования течения. На рис. 5 пред- ставлено поле завихренности уже сформировавше- гося течения. Видно, что вихри, периодически сво- рачивающиеся в верхнем и нижнем сдвиговых сло- ях, также будут располагаться в шахматном по- рядке относительно друг друга. Перейдем к сравнительному анализу течений, рассмотренных в [49] и в данной статье. На рис. 6, а,б представлена картина линий тока и по- ля завихренности, соответственно, для случая уз- ких межстенозных полостей (D2 = D3 = 2 мм, D4 = 6 мм). Как хорошо видно, в данном случае течение существенно отличается от того, которое было описано и проанализировано в [49], т. е. от течения, в области с широкими межстенозными полостями (D2 = D3 = 2 мм, D4 = 18 мм). Пре- жде всего следует отметить, что в [49] в каждой межстенозной полости присутствовал только один крупный вихрь, а более мелкие вихри располага- лись вблизи углов полости, наиболее удаленных от струи. В случае же узких межстенозных полостей (D2 = D3 = 2 мм, D4 = 6 мм), рассмотренных в данной статье, в каждой полости отчетливо видны два горизонтально расположенных относительно друг друга крупных вихря. Поскольку большие слоистые вихри, образующиеся в межстенозных нишах, являются по сути каналом обратной связи, то можно ожидать, что изменение размера, фор- мы и количества больших вихрей в межстенозных полостях может существенным образом изменить возмущение самой струи. Действительно, в [49] было продемонстрирова- но, что как только энергии в канале обратной свя- зи оказывается достаточно, чтобы влиять на фор- му струи (см. картину линий тока) и на сдви- говый слой (см. картину поля завихренности), в сдвиговых слоях начинает генерироваться хара- ктерная последовательность вихрей. При этом ря- ды вихрей в верхнем и нижнем сдвиговых сло- ях располагаются относительно друг друга в ша- хматном порядке. Набегая на отверстие второго стеноза, эти ряды вихрей вызывают несимметри- чные (относительно оси струи) колебания верти- кального профиля скорости в выходном сечении S2. В случае же узких межстенозных полостей (D2 = D3 = 2 мм, D4 = 6 мм), рассмотренных в данной статье, видно, что возмущения струи не приводят к потере симметричности поля течения относительно оси. Последовательность вихрей, ге- нерируемая в сдвиговых слоях и набегающая вме- сте с потоком на отверстие второго стеноза, остае- тся расположенной симметрично относительно оси потока и вызывает симметричные колебания вер- тикального профиля скорости в выходном сечении S2, в то время как профиль скорости во входном сечении S1 остается неподвижным (рис. 6, в, г). Естественно возникает вопрос: являются ли ко- лебания профиля скорости в сечении S2 периоди- ческими? Для ответа на него мы исследовали зави- симость изменения профиля скорости от времени, и оказалось, что при рассматриваемом значении числа Рейнольдса эти колебания имеют периоди- ческий характер и симметричны относительно оси струи. Таким образом, в рассматриваемой гидро- динамической системе действительно могут возни- кать автоколебательные явления. С целью иллю- страции этих явлений на рис. 7, приведены дан- ные, которые показывают как в течение одного пе- риода T процесса автоколебаний в системе изме- няются картины поля завихренности и осцилли- рующие части профилей скорости Vx в выходном сечении S2. Заметим, что при разложении вели- чины Vx в ряд Фурье по времени стационарная (не зависящая от времени) часть профиля скоро- 58 В. С. Малюга ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 0 0.5 -0.5 -1 1 X10 4 а б 0 2 4 6 8 10 12 8 8.5 9 9.5 10x 10 -3 Ux y (м) (м/с) 0 2 4 6 8 10 12 8 8.5 9 9.5 10 x 10 -3 Ux y (м) (м/с) в г Рис. 6. Вид течения в момент времени t = 9.295 · 10 −3 с: а – линии тока, б – поле завихренности, в, г – профили скорости в сечениях S1 и S2, соответственно сти описывается нулевым членом ряда, а осцил- лирующая часть, которая и представлена на ри- сунке, – суммой остальных членов ряда. Следует также отметить, что в данной задаче период сим- метричных автоколебаний будет почти вдвое боль- ше, чем период колебаний в соответствующей за- даче, рассмотренной в [49]. Так, при V1 = 0.9 м/с период антисимметричных колебаний был равен T = 2.77 · 10−4. В случае же симметричного те- чения, рассматриваемого в данной статье, период T = 5.2 · 10−4. ЗАКЛЮЧЕНИЕ 1. Представлен алгоритм численного решения системы уравнений Навье-Стокса для несжима- емой жидкости, разработанный на основе тех- ники конечно-объемной дискретизации. Исполь- зуемые разностные схемы имеют второй поря- док точности по пространству и времени. Для дискретизации конвективных членов используется TVD схема Chakravarthy-Osher. Уравнения реша- ются в неподвижной декартовой системе коорди- нат. При этом сетка, которая разбивает расчетную область на контрольные объемы, рассматривается как трехмерная неструктурированная. Это сдела- но с целью дальнейшего развития представлен- ных численных алгоритмов на случай сеток, состо- ящих из полихедральных контрольных объемов с произвольным числом соседей, что существен- но упрощает проблему генерации сетки для сло- жных геометрий. Для решения полученных СЛАУ использовались солверы, основанные на методе (би-)сопряженных градиентов с предобусловлива- нием, доступные в открытых библиотеках. 2. Представленный алгоритм протестирован на хорошо известных задачах: задаче о течении в прямоугольной полости, вызванным движением одной из стенок (lid-driven cavity), при различных числах Рейнольдса, а также задаче об отрывном течении за круговым цилиндром (вихревая дорож- ка Кармана). Тест показал хорошее совпадение по- лученных результатов с результатами других ав- торов. 3. С использованием представленного алгорит- ма проведено прямое численное моделирование те- чения вязкой несжимаемой жидкости в плоском канале при наличии двух последовательно распо- ложенных двусторонних сужений (стенозов). Рас- считаны поля скорости и давления, представлены В. С. Малюга 59 ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 -3 -2 -1 0 1 2 3 8 8.5 9 9.5 10 x 10 -3 Ux y (м) (м с)/ а б 8 8.5 9 9.5 10 x 10 -3 y (м) -3 -2 -1 0 1 2 3 Ux (м с)/ в г x 10 -3 (м) 8 8.5 9 9.5 10 y -3 -2 -1 0 1 2 3 Ux (м с)/ д е 8 8.5 9 9.5 10 x 10 -3 y (м) -3 -2 -1 0 1 2 3 Ux (м с)/ ж з 8 8.5 9 9.5 10 x 10 -3 y (м) -3 -2 -1 0 1 2 3 Ux (м с)/ и к Рис. 7. Осциллирующая часть профиля скорости Vx в сечении S2 (а, в, д, ж, и), и поле завихренности (б, г, з, е, к) для моментов времени: а, б – t = 8.775 · 10 −3 с, в, г – t = 8.905 · 10 −3 с, д, е – t = 9.035 · 10 −3 с, ж, з – t = 9.165 · 10 −3 с, и, к – t = 9.295 · 10 −3 с 60 В. С. Малюга ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 картины поля завихренности и проанализированы особенности движения среды в зависимости от ши- рины канала в межстенозной части. 4. Установлено, что при достаточно малых чи- слах Рейнольдса (в рассмотренном примере Re = 3132) и при достаточно узкой межстенозной части канала (в рассмотренном примере ширина отвер- стия D2 и ширина полости (D4 − D2)/2 соотноси- лись как 1/1) в сдвиговых слоях в межстенозной области будут периодически во времени генериро- ваться вихревые структуры симметричные отно- сительно оси канала y = D1/2. А, следовательно, и в колебаниях профиля скорости в отверстии вто- рого стеноза будет доминировать симметричная часть. Также показано, что с увеличением шири- ны межстенозной части канала (в рассмотренных примерах высота отверстия D2 и глубина полости (D4 − D2)/2 соотносились как 1/2 и 1/4) карти- на течения принципиально меняется. В сдвиговых слоях в межстенозной области образуются два ря- да вихрей. Эти вихри располагаются относитель- но друг друга в шахматном порядке. Как пока- зано в [49], такие вихри, набегая на второй сте- ноз, вызывают антисимметричные периодические колебания профиля скорости в отверстии второго стеноза. 1. Патанкар С. Численные методы решения задач те- плообмена и динамики жидкости.– М.: Энергоато- миздат, 1984.– 152 с. 2. Hirsch C. Numerical computation of internal and external flows.– Oxford: Butterworth-Heinemann, 2007.– 656 p. 3. Versteeg H. K., Malalasekera W. An introduction to computational fluid dynamics. The finite volume method.– New York: Longman, 1995.– 258 p. 4. Годунов С.К. Разностный метод численного ра- счёта разрывных решений уравнений гидродина- мики // Мат. сборник.– 1959.– 47, N 3.– С. 271– 306. 5. Courant R., Isaacson E., Rees M. On the solution of non-linear hyperbolic differential equation by fnite differences // Comm. Pure Appl. Math.– 1952.– 5.– P. 243–255. 6. Lax P.D. Weak solutions of non-linear hyperbolic equations and their numerical computation // Comm. Pure Appl. Math.– 1954.– 7.– P. 159–193. 7. Gentry R. A., Martin R. E., Daly B. J. An Eulerian differencing method for unsteady compressible flow problems // J. Comp. Physics.– 1966.– 1.– P. 87–118. 8. Lax P.D., Wendroff B. Systems of conservation laws // Comm. Pure Appl. Math.– 1960.– 13.– P. 217–237. 9. MacCormack R. W. The effect of viscosity in hypervelocity impact cratering // AIAA Paper 69- 354, 1969; доступно также в Caughey D. A. (ed), Hafez M. M. (ed) Frontiers of computational fluid dynamics.– London: World Scientific, 2002.– 524 p. 10. Lerat A., Peyret R. Noncentered schemes and shock propagation problems // Computers and Fluids.– 1974.– 2.– P. 35–52. 11. MacCormack R. W. A numerical method for solving the equations of compressible viscous flow // AIAA Paper 81-0110, 1981 12. Lerat A. Implicit methods of second order accuracy for the Euler equations // AIAA Journal.– 1985.– 23.– P. 33–40. 13. Beam R. M., Warming R. F. An implicit finite- difference algorithm for hyperbolic systems in conservation-law form // J. Comp. Physics.– 1976.– 22.– P. 87–110. 14. Beam R. M., Warming R. F. An implicit factored scheme for the compressible Navier-Stokes equati- ons // AIAA Journal.– 1978.– 16.– P. 393–402. 15. Jameson A., Schmidt W., Turkel E. Numerical si- mulation of the Euler equations by finite volume methods using Runge-Kutta time stepping schemes // AIAA Paper 81-1259, 1981 16. Wong H. H., Raithby G. D. Improved finite-difference methods based on a critical evaluation of the approxi- mating errors // Numerical Heat Transfer. Part B: Fundamentals.– 1979.– 2.– P. 139–163. 17. Raithby G. D., Torrance K. E. Upstream-weighted schemes and their application to elliptic problems involving fluid flow // Computers and Fluids.– 1974.– 2.– P. 191–206. 18. Warming R. F., Beam R. M. Upwind second order difference schemes and applications in aerodynamic flows // AIAA Journal.– 1976.– 14.– P. 1241–1249. 19. Leonard B. P. A stable and accurate convecti- ve modelling procedure based on quadratic upstream interpolation // Comp. Meth. Appl. Mech. Engineering.– 1979.– 19.– P. 59–98. 20. Raithby G. D. A critical evaluation of upstream di- fferencing applied to problems involving fluid flow // Comp. Meth. Appl. Mech. Engineering.– 1976.– 9.– P. 75–103. 21. Raithby G. D. Skew upstream differencing schemes for problems involving fluid flow // Comp. Meth. Appl. Mech. Engineering.– 1976.– 9.– P. 153–164. 22. Van Leer B. Upwind and High-Resolution Methods for Compressible Flow: From Donor Cell to Residual- Distribution Schemes // Commun. Comput. Phys..– 2006.– 1.– P. 192–206. 23. Spalding D. B. A novel finite-difference formulati- on for differential expression involving both fi- rst and second derivatives // Int. J. Num. Meth. Engineering.– 1972.– 4.– P. 551–559. 24. Allen D. N. de G., Southwell R. V. Relaxation Methods Applied to Determine the Motion, in Two Dimensions, of a Viscous Fluid Past a Fixed Cyli- nder // Quart. J. Mech. and Appl. Math..– 1955.– 8.– P. 129–145. 25. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычи- слительная гидромеханика и теплообмен, Т.2.– М: Мир, 1990.– 392 с. 26. Perić, M. A Finite Volume method for the prediction of three-dimensional fluid flow in complex ducts. PhD Thesis.– London: Imperial College, 1985.– 394 p. 27. Boris J. P., Book D. L. Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works // J. Comp. Physics.– 1973.– 11.– P. 38–69. 28. Zalesak S. T. Fully multidimensional flux-corrected transport algorithms for fluids // J. Comp. Physics.– 1979.– 31.– P. 335–362. 29. Van Leer B. Towards the ultimate conservative differencing scheme // Proceedings of the Third International Conference on Numerical Methods in Fluid Mechanics, volume 1, ed. by Cabannes, H. and Temem, R..– Springer, 1973.– P. 163–168. В. С. Малюга 61 ISSN 1561 -9087 Прикладна гiдромеханiка. 2010. Том 12, N 4. С. 45 – 62 30. Van Leer B. Towards the ultimate conservative di- fferencing scheme. II. Monotonicity and conservati- on combined in a second-order scheme // J. Comp. Physics.– 1974.– 14.– P. 361–370. 31. Van Leer B. Towards the ultimate conservative differencing scheme. III. Upstream-centered finite- difference schemes for ideal compressible flow // J. Comp. Physics.– 1977.– 23.– P. 263–275. 32. Van Leer B. Towards the ultimate conservative di- fferencing scheme. IV. A new approach to numerical convection // J. Comp. Physics.– 1977.– 23.– P. 276– 299. 33. Van Leer B. Towards the ultimate conservative differencing scheme. V. A second-order sequel to Godunov’s method // J. Comp. Physics.– 1977.– 23.– P. 101–136. 34. Harten A. High resolution schemes for hyperbolic conservation laws // J. Comp. Physics.– 1983.– 49.– P. 357–393. 35. Harten A. On a class of high resolution total- variation-stable finite-difference schemes // SIAM J. Numer. Analysis.– 1984.– 21.– P. 1–23. 36. Roe P. L. Large scale computations in fluid mechani- cs, Part 2 // Lectures in Applied Mathematics, volume 22.– Springer, 1985.– P. 163–193. 37. Chakravarthy S. R., Osher S. High resolution appli- cation of the OSHER upwind scheme for the Euler equation // AIAA Paper 83-1943, 1983 также доступно в Chakravarthy S.R., Osher S. High resolution appli- cation of the Osher upwind scheme for the Euler equation // Proc. AIAA Comp. Fluid Dynamics conference.– Danvers, MA, 1983.– P. 363–372. 38. Osher S., Chakravarthy S. High resolution schemes and the entropy condition // SIAM J. Numer. Analysis.– 1984.– 21.– P. 955–984. 39. Sweby P. K. High resolution schemes using flux li- miters for hyperbolic conservation laws // J. Numer. Anal..– 1984.– 21.– P. 995–1011. 40. Waterson N.P., Deconinck H. Design principles for bounded higher-order convection schemes – a unified approach // J. Comp. Physics.– 2007.– 224.– P. 182– 207. 41. Chorin A. J. A numerical method for solving incompressible viscous flow problems // J. Comp. Physics.– 1967.– 2.– P. 12–26. 42. Harlow F. H., Welch J. E. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface // Phys. Fluids.– 1965.– 8.– P. 2182– 2189. 43. Patankar S. V., Spalding D. B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows // Int. J. Heat Mass Transfer.– 1972.– 15.– P. 1787–1806. 44. Van Doormal J. P., Raithby G. D. Enhancements of the simple method for predicting incompressible fluid flows // Numer. Heat Transfer.– 1984.– 7.– P. 147– 163. 45. Issa R. I. Solution of implicitly discretised fluid flow equations by operator-splitting // J. Comput. Phys..– 1986.– 62.– P. 40–65. 46. Jang D. S., Jetli R., Acharya S. Comparison of the PISO, SIMPLER, and SIMPLEC algorithms for the treatment of the pressure-velocity coupling in steady flow problems // Numer. Heat Transfer, Part A: Applications.– 1986.– 10.– P. 209–228. 47. Arakawa A. Computational design for long-term numerical integration of the equations of fluid moti- on: Two-dimensional incompressible flow. Part I // J. Comput. Phys..– 1966.– 1.– P. 119–143. 48. Rhie C. M., Chow W. L. Numerical study of the turbulent flow past an airfoil with trailing edge separation // AIAA J..– 1983.– 21.– P. 1525–1532. 49. Вовк И. В., Гринченко В. Т., Малюга В. С. Особен- ности движения среды в каналах со стенозами // Прикл. гiдромех.– 2009.– 11, N 4.– С. 17–30. 50. Басовский В. Г., Вовк И. В., Вовк О. И. О воз- можности генерирования тональных звуковых ко- лебаний потоком воздуха в бронхах со стенозом // Акуст. Вiсник.– 2003.– 6, N 1.– С. 3–21. 51. Бруяцкий Е. В., Костин А. Г. Чиссленное исследо- вание течения жидкости в закрытой прямоуголь- ной полости с движущейся верхней крышкой // Прикл. гiдромех.– 2009.– 11, N 1.– С. 3–15. 52. Приходько А. А., Редчиц Д. А. Численное мо- делирование нестационарного течения в следе за цилиндром на основе уравнений Навье-Стокса // Прикл. гiдромех.– 2005.– 7, N 1.– С. 56–71. 53. Ferziger J. H., Perić M. Computational methods for fluid dynamics, 3rd rev. ed.– Berlin: Springer, 2002.– 424 p. 54. Роуч П. Вычислительная гидродинамика.– Мю: Мир, 1980.– 612 с. 55. http://www.opencfd.co.uk/openfoam/ 56. Barrett R., Berry M., Chan T. F. et al. Templates for the solution of linear systems: Building blocks for iterative methods, 2nd Edition.– Philadelphia: SIAM, 1994.– 107 p. 57. Khosla P. K., Rubin S. G. A diagonally dominant second-order accurate implicit scheme // Computers & Fluids.– 1974.– 2.– P. 207–209. 58. Hestens M. R. and Steifel E. L. Method of conjugate gradients for solving linear systems // J. Res.– 1952.– 29.– P. 409–436. 59. Jacobs D. A. H. Preconditioned Conjugate Gradient methods for solving systems of algebraic equations // Central Electricity Research Laboratories Report, RD/L/N193/80.– Leatherhead, Surrey, 1980.– P. 31. 60. Van Der Vorst H. A. Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems // SIAM J. Scientific Computing.– 1992.– 13(2).– P. 631–644. 61. Van Der Vorst, H.A. Iterative Krylov methods for large linear systems.– Cambridge: Cambridge univ. press, 2003.– 221 p. 62. http://netlib.org 62 В. С. Малюга

Численное исследование течения в канале с двумя последовательно расположенными стенозами. Алгоритм решения

Репозитарії

Схожі ресурси