Содержание к диссертации
Введение
1 Обзор известных исследований 7
1.1 Ламинарное установившееся движение вязкой несжимаемой жидкости в трубе эллиптического сечения 7
1.2 Линейная устойчивость течения Пуазейля в трубе эллиптического сечения 12
1.3 Результаты расчетов турбулентных течений в областях различной формы 15
1.4 Краткая сводка вычислительных подходов 21
2 Вычислительный метод 26
2.1 Постановка задачи 26
2.2 Пространственная дискретизация 29
2.3 Определение давления 35
2.4 Схема интегрирования по времени 38
3 Результаты расчета 45
3.1 Алгоритмические и геометрические параметры задачи 45
3.2 Характеристики течения в установившемся режиме 50
3.3 Вторичные течения 67
3.4 Баланс кинетической энергии пульсаций 74
Заключение 78
Список литературы
- Линейная устойчивость течения Пуазейля в трубе эллиптического сечения
- Краткая сводка вычислительных подходов
- Пространственная дискретизация
- Характеристики течения в установившемся режиме
Введение к работе
Прямой расчет, основанный на численном решении полных уравнений Навье-Стокса, превратился за последние годы в надежный метод исследования физики турбулентных течений. Многочисленные примеры свидетельствуют о хорошем согласии результатов прямых расчетов с экспериментом [1]. Данный метод не требует какой-либо априорной информации в виде эмпирических констант или подгоночных коэффициентов для расчета исследуемого течения. Такой подход зачастую является даже предпочтительным по сравнению с экспериментальными подходами, поскольку обеспечивает практически неограниченную по детальности информацию о структуре исследуемого течения.
Главная трудность на пути численного изучения турбулентности связана с величиной диапазона динамически значимых масштабов, свойственных турбулентным течениям, которые должны быть адекватно разрешены на дискретном уровне. Для расчета пристенных турбулентных течений при больших числах Рейнольдса объем вычислений пропорционален Re7'2. Последнее ограничивает численные исследования весьма невысокими числами Рейнольдса. Другое существенное ограничение связано со степенью сложности геометрии течений, доступных для прямого расчета.
В настоящее время довольно подробно изучены как численно, так и экспериментально одномерные в среднем течения. В классе пристенных течений большинство исследований проведено для плоского канала, асимптотического погранслоя над плоской пластиной, трубы круглого сечения [2]—[8].
Течения, средние характеристики которых зависят от двух координат, например, течения в некруглых трубах, интересны в прикладном и научном плане не только большей пространственной сложностью, но и наличием так называемых турбулентных вторичных течений, именуемых также вторичными течениями Прандтля 2-ого рода. Вторичные течения — это организованные движения жидкости в плоскости, перпендикулярной к направлению основного потока. В отличие от вторичных течений Прандтля 1-ого рода,
возникающих в потоках вдоль вогнутой поверхности под действием центробежных сил как в турбулентных, так и в ламинарных потоках, вторичные течения Прандтля 2-ого рода — исключительно турбулентное явление, вызываемое анизотропией компонент тензора напряжений Рейнольдса. Интенсивность турбулентных вторичных течений невелика (как правило 1-3% от средней скорости потока), однако их вклад в процессы переноса импульса, массы, примеси в поперечной к направлению потока плоскости весьма значителен. Непосредственное измерение вторичных течений в экспериментальных условиях затруднительно, поскольку их величина сравнима с точностью измерений. Отсутствие достоверных экспериментальных данных задерживает разработку приближенных методов расчета таких течений [9]. В этих условиях прямой расчет оказывается практически единственным источником надежной информации о свойствах и структуре вторичных течений в некруглых трубах. До недавнего времени численное исследование турбулентных течений с неодномерными средними характеристиками ограничивалось трубами прямоугольного сечения [10]—[12].
Эллиптическая труба является незначительной модификацией классической трубы и простейшим типом трубы некруглого сечения. Несмотря на это для турбулентных течений вязкой несжимаемой жидкости в трубах эллиптического сечения в литературе до недавнего времени отсутствовали какие-либо данные об их свойствах и структуре. Для данного типа труб были лишь подробно описаны свойства ламинарного течения вязкой несжимаемой жидкости [13] и проведено исследование линейной устойчивости ламинарного течения [14]. Целью настоящей работы являлись прямой расчет и анализ турбулентных течений в трубах эллиптического сечения, описание интегральных и пульсационных характеристик, изучение структуры вторичных течений.
Линейная устойчивость течения Пуазейля в трубе эллиптического сечения
Ламинарное течение в эллиптической трубе служит связующим звеном между течением Гагена—Пуазейля в круглой трубе и плоским течением в канале. Если отношение длин осей эллипса близко к единице (случай малой эллиптичности), течение аппроксимирует течение Гагена—Пуазейля, которое считается устойчивым к бесконечно малым возмущениям. Наоборот, когда эллиптичность велика, течение близко к течению в плоском канале, для которого критическое число Рейнольдса существует.
Хоукинг [20] рассмотрел случай больших отношений длин осей и получил асимптотический результат для критического числа Рейнольдса: Recru = 5772 + 86300/Д где А — отношение длин полуосей эллипса. В той же работе им было исследовано, как поперечная кривизна границы может стабилизировать плоское течение Пуазейля в канале. Потер, Хайнс, Рейнольде и Потер [21]—[23] описали стабилизацию плоского течения Пуазейля при помощи добавления продольного сдвигового движения граничных плоскостей. Тем самым устанавливается смесь плоского течения Пуазейля и плоского течения Куэтта. Линейная устойчивость достигается для любого числа Рейнольдса, если скорость плоскостей превышает 53% от средней скорости сквозного потока. Татсуми и Яшимура [24] изучали стабилизирующий эффект наложения боковых стенок на течение в канале с помощью исследования течения в трубе квадратного сечения. Начиная с большого отношения длин осей и переходя к пределу — квадратной трубе, авторы определили, что критическое число Рейнольдса становится бесконечным, когда отношение длин в трубе приблизительно равно 3.2.
Однако Хоукинг [20] показал лишь относительную важность для критического числа Рейнольдса кривизны поперечной границы и присутствия боковых сторон. Для прямоугольной трубы с отношением длин сторон А 1 им было найдено, что критическое число Рейнольдса возрастает лишь на величину 6844/А2 по сравнению с течением в канале. При этом в случае эллиптической трубы кривизна границы влияет намного сильнее: критическое число возрастает на 86300/А Отсюда, в частности, можно сделать вывод, что течение в эллиптической трубе линейно устойчиво при отношении длин 3.2, при котором течение в прямоугольной трубе уже теряет линейную устойчивость. Керсвел и Даффи [14] провели исследование линейной устойчивости течения Пуазейля в трубе с эллиптическим поперечным сечением с конечным соотношением длин полуосей. С помощью увеличения эллиптичности они провели плавный переход от линейно устойчивого течения в круглой трубе к течению в канале. Необходимо было установить, при каком отношении длин полуосей проявится линейная неустойчивость, то есть появится конечное критическое число Рейнольдса. В эллиптической трубе для случая конечных отношений длин полуосей и околокритических чисел Рейнольдса единственным механизмом изучения линейной устойчивости является решение двумерной задачи на собственные числа с использованием численных методов. На фиг. 1.2 сплошной линией показана полученная из расчетов нейтральная кривая устойчивости в плоскости (Re, A). Re—число Рейнольдса, базирующееся на максимальной скорости основного течения и на длине меньшей полуоси. На том же рисунке показан асимптотический результат Re = 5772 + 86300/А из работы Хоукинга [20] для большого отношения длин полуосей (линия с коротким пунктиром). Линией с длинным пунктиром обозначено постоянное значение Re = 5772, соответствующее критическому числу Рейнольса для течения в канале.
Из рисунка видно, что критическое отношение длин, при котором течение становится линейно неустойчивым, приблизительно равно 10.4. Это могло бы означать, что в случае, когда длины осей отличаются больше, чем на порядок, течение в трубе начинает больше походить на течение в плоском канале, чем в круглой трубе. Однако Хоукингом [20] показано, что для течения в трубе эллиптического сечения критическое число Рейпольдса очень чувствительно к отношению длин: при А — 50, оно все еще равно 8000.
Краткая сводка вычислительных подходов
В первых работах по численному моделированию турбулентности применялись конечно-разностные методы аппроксимации искомых функций по пространственным координатам. При этом наиболее плодотворным оказался подход, использующий перемежающиеся сетки, разработанный Харлоу и Уэл-шем [71] и обобщенный на случай цилиндрической геометрии для области между двумя концентрическими цилиндрами Вильямсом [54]. В этом подходе значения скалярных величин — давления, температуры, концентрации и т.д. отыскиваются в центрах расчетных ячеек, а значения компонент скорости отыскиваются в узлах, лежащих на гранях ячеек. Такое разнесение узлов в пространстве позволяло создавать схемы, обладающие рядом полезных аналогов законов сохранения, свойственные исходным дифференциальным уравнениям.
Альтернативным подходом являются спектральные методы. В этих методах искомые функции представляются некоторыми аналитическими выражениями со свободными коэффициентами, значения которых определяются из тех или иных критереев минимальности невязки уравнений. Наиболее популярным из спектральных методов является метод Галеркина, в котором решение краевой задачи ищется в виде разложения по некоторой полной системе базисных функций, удовлетворяющей граничным условиям, а значения коэффициентов разложения находятся из условия ортогональности невязки уравнений функциям базисной системы. В задачах гидромеханики метод Галеркина стал применяться после работы Г.И. Петрова [55], где показано, что для отыскания коэффициентов разложения можно использовать другую систему функций, являющуюся базисом в более широком функциональном пространстве. Метод Галеркина—Петрова с различными модификациями интенсивно используется при численном моделировании турбулентности. Первоначально спектральные методы использовались лишь в линейных задачах. Их неудобство в нелинейных проблемах связано с тем, что для вычисления коэффициентов разложения произведения двух функций требуется вычис лять свертки рядов. Это очень дорогая (в смысле числа необходимых операций) процедура делает алгоритмы неэффективными, если число удерживаемых членов ряда достаточно велико. Ситуация кардинально изменилась с изобретением Кули и Тыоки [56] алгоритма Быстрого Преобразования Фурье (БПФ). Значительное число идей повышения эффективности спектральных методов, использующих алгоритм БПФ, предложено Орзагом [57].
Орзагом и Келсом [58] разработан полностью спектральный алгоритм моделирования турбулентных течений в плоском канале, близкие алгоритмы применялись Рождественским и Симакиным [2], [59] и Кимом, Мойном и Мозером [3]. Спектральный метод с использованием многочленов Чебыше-ва применялся Пономаревым, Приймаком и Рождественским [34] для трубы кольцевого сечения. Полностью спектральные алгоритмы решения уравнений Навье—Стокса для трубы круглого сечения применялись Приймаком и Рождественским [42] и Приймаком [5]. Спаларт [4] для моделирования турбулентности в пограничном слое использовал многочлены Якоби. Векторные собственные функции задачи Стокса использовались Н. В. Никитиным [60],[61] при изучении развития возмущений конечной амплитуды в плоском канале. Этот базис также применялся Никитиным [62], [41], [63] для моделирования течений в трубе круглого сечения и Герценштейном и Никитиным [64] для расчета вторичных течений во вращающейся трубе.
Считается, что в случае гладких решений, конечно-разностные методы серьезно уступают в эффективности спектральным. Несколько сравнительных примеров приведено Орзагом [65]. Это имеет место, в частности, в предтурбу-лентных движениях жидкости, когда в потоке доминируют крупномасштабные структуры. Оказывается, что в развитых турбулентных течениях сходимость рядов довольно слабая, в силу чего спектральные методы теряют значительную часть своих преимуществ. Раи и Моин [66] провели моделирование турбулентности в плоском канале полностью конечно-разностным методом. Оказалось, что даже в этом, наиболее удобном для реализации спектральной идеалогии случае, конечно-разностный подход не уступает в эффективности.
Решение трехмерных уравнеий Навье—Стокса в цилиндрической области связано с дополнительными трудностями из-за геометрической особенности на оси симметрии. Конечно-разностный метод аппроксимации уравнений в радиальном направлении с удобной постановкой условий на оси симметрии был разработан Никитиным [67]. Усовершенствованный его вариант, включающий метод смягчающего оператора при интегрировании по времени, приведен в работе Никитина [7]. При сравнении спектральных и конечно-разностных методов следует учитывать также, что последние являются гораздо более простыми для реализации и более гибкими по отношению к изменению геометрии течения.
В [49], [51] Никитиным была предложена общая методика моделирования турбулентных течений, ориентрованная на расчет течений в областях, имеющих форму параллелепипеда в некоторой ортогональной системе координат. Тогда границы области течения являются координатными поверхностями, и пристенные области могут быть подробно разрешены на простой расчетной сетке. Единым принципом дискретизации уравнений Навье—Стокса по пространственным координатам является постановка условий периодичности в тех направлениях, в которых исследуемое течение может считаться однородным, и применение конечно-разностной дискретизации по пространству с сохранением ряда важных локальных и интегральных свойств исходных дифференциальных уравнений.
Пространственная дискретизация
Первое слагаемое в правой части последнего тождества равно нулю в силу условия несжимаемости, а второе дает ноль после интегрирования по Q, применения формулы Гауса—Остроградского и учета условий непротекания на границе. Таким образом, нелинейные члены и градиент периодической части давления не дают интегрального вклада в изменение энергии, а лишь перераспределяют ее между различными районами области течения. Далее используется тождество и rot ш = и) rot u + div (ш х и)
Интеграл от второго слагаемого в правой части равен нулю из-за выполнения условий прилипания (нулевые значения касательных скоростей на неподвижных твердых.стенках), и уравнение принимает вид: = -uj(wtu)2 + DpJu2, E=\j\u\2 (23)
Первое слагаемое в правой части уравнения баланса энергии описывает вязкую диссипацию, второе — мощность внешних сил (в рассматриваемом случае среднего градиента давления). Следствием условия прилипания при неподвижных границах является отрицательность вклада вязких членов. Для вывода уравнения баланса энергии в дискретном случае кинетическая энергия определяется как где суммирование проводится по всем р-узлам в области течения. Легко показать, что при таком определении вклад членов градиента периодической части давления q в (20)— (22) равен нулю. Слагаемые, отвечающие вкладу нелинейных членов, раскладываются на три пары, соответствующие трем компонентам завихренности. Каждая такая пара дает нулевой вклад. Например, вклады первого члена в (20) и второго в (21) совпадают между собой (множитель ArA5iA2 опущен):
Поэтому они нейтрализуют друг друга в уравнении баланса кинетической энергии. Для вклада вязких членов получается выражение аппроксимирующее выражение для вязкой диссипации в (23).
По известным значениям компонент скорости вычисляются завихренности, после чего можно вычислить нелинейные и вязкие члены в (20)— (22). Обозначим правые части (20)— (22) кроме членов с давлением через v = v (u) = (vr, ve, vz). Уравнения запишутся в виде:
В силу условий непротекания на границе dur/dt = 0 и, следовательно, Srq = Hf vr. Далее можно вычислить дивергенцию поля du/dt и положить ее равной нулю во всех р-узлах. Результатом будет уравнение для давления: Div Grad q = Div v, 5rq = Hf vr при r = r , (25) где Div, Grad — разностные операторы дивергенции и градиента, Div Grad — разностный оператор Лапласа. Таким образом, для давления получается задача Неймана для уравнения Пуассона. Условием разрешимости задачи для давления является условие постоянства объема жидкости в области течения (поток жидкости через границы области равен нулю), что очевидно необходимо для несжимаемости.
Для определения среднего градиента давления Dp интегрируется (в разностном смысле) последнее из уравнений (24) по всей расчетной области:
Суммы вычисляются по всем точкам определения компоненты uz. В силу несжимаемости жидкости интеграл от uz по площади поперечного сечения трубы равен расходу жидкости Q (и) и не зависит от продольной координаты, поэтому сумма под знаком производной в левой части (26) равна LzQ(u). Если расход не меняется со временем (как в данной постановке),
Характеристики течения в установившемся режиме
Пространственные дискретизации уравнений Навье—Стокса при достаточно подробном разрешении приводят к жестким системам и требуют неявных методов для интегрирования по времени. Полностью неявные методы порождают системы связанных нелинейных уравнений для зависимых переменных на новом временном слое и поэтому, как правило, неприемлемы для длительного интегрирования, в частности, в задачах прямого расчета турбулентных течений. Полунеявные методы, в которых только часть оператора Навье— Стокса выносится на новый временной слой, представляют разумный компромисс для этого класса задач. В действительности большинство успешных расчетов турбулентных течений проведено с использованием полунеявных методов. В расчетах пристенных течений чаще всего лишь линейные вязкие члены уравнений Навье—Стокса обрабатываются неявно. Получающиеся при этом системы линейных уравнений эффективно решаются с помощью метода приближенной факторизации.
Одна из наиболее популярных схем была предложена в [72] и включает метод Адамса—Бэшфорда для конвективных членов и метод Кранкла— Николсона для вязких. Известно, что метод Адамса—Бэшфорда (как и методы Рунге—Кутты первого и второго порядков) безусловно неустойчив для уравнения конвективного переноса. Однако отмеченная неустойчивость достаточно слабая, и при наличии вязких членов метод обычно работает при невысоких числах Куранта. Методы Рунге—Кутты третьего порядка точности условно устойчивы, поэтому они представляются более подходящими для конвективных членов уравнений Навье—Стокса. Несколько близких вариантов полунеявных схем, основанных на методах Рунге—Кутты третьего порядка, предложены в [бб], [74], [75], а также в [8]. Это трехшаговые схемы, требующие в три раза больше операций для продвижения на следующий временной уровень по сравнению со схемой [72]. Тем не менее, благодаря высокой точности и устойчивости они оказываются предпочтительными в сравнении со схемами, основанными на методе Адамса—Бэшфорда.
В отмеченных схемах, использующих метод Рунге—Кутты третьего порядка для конвективных членов, вязкие члены обрабатываются неявно со вторым порядком. Таким образом, суммарный порядок точности этих схем остается вторым. По всей видимости, невозможно добиться суммарно третьего порядка точности в рамках классической трехшаговой схемы, каждый шаг которой состоит из трех стадий: а) вычисления конвективных и вязких членов; б) решения уравнения Пуассона для давления (или псевдодавления); в) решения линейной системы, связаной с неявностью.
В данной работе используется полунеявная схема Рунге—Кутты для несжимаемых уравнений Навье—Стокса, обладающая третьим порядком точности [76] в отличие от других полунеявных трехшаговых схем [8], [66], [74], [75], обладающих лишь вторым порядком. Повышенный порядок точности достигается благодаря тому, что на последнем шаге выполняется одна дополнительная стадия (в). Кроме того, удовлетворение уравнению неразрывности производится в последнюю очередь, чтобы окончательное поле было бездивергентным и обладало нужным расходом. Среди стадий (а) - (в) решение уравнения Пуассона является наиболее дорогостоящей стадией по затратам машинного времени. Что касается стадии (в), то она не требует значительного времени при любом способе пространственной дискретизации при использовании метода приближенной факторизации.
Одна из проблем при построении неявных схем высокого порядка для несжимаемых уравнений Навье—Стокса это так называемая проблема давления. Она возникает вследствие того, что компоненты скорости и давление образуют связанную систему с уравнением неразрывности. Почти все численные схемы решения уравнений Навье—Стокса в простейших переменных (скорость - давление) используют методы расщепления, в которых сначала вычисляется промежуточное поле скорости без учета условия несжимаемости, а затем оно проецируется на пространство бездивергентных векторов. Формулировка граничных условий для промежуточного поля скорости и для связаных с давлением величин требует особого внимания. Эти вопросы об суждаются в литературе в течение многих лет (см., например, [77]). Представляемая схема формулируется в терминах дискретных уравнений, в которые уже включены граничные условия. При таком подходе [78], [79] промежуточные поля скорости не нуждаются в каких-либо специальных условиях, что весьма облегчает построение схемы интегрирования по времени. В применяемой схеме только окончательное поле скорости на каждом полном шаге удовлетворяет условию несжимаемости. На промежуточных этапах безди-вергентность скорости выполняется лишь с некоторой погрешностью. Такой подход позволяет относительно просто преодолеть проблему давления при сохранении желаемого порядка точности схемы.
Неявная часть является линейным оператором, включающим в себя те члены уравнений, которые производят наибольшую жесткость. Пусть L— неявный линейный оператор, 7—некоторое положительное число, вязкие члены обозначены через v = v(u). Для вычисления поля скорости u = ип+1 в момент времени t = tn+1 = tn + At по известному полю u = un в момент времени t = tn используется следующая процедура (через Q(u) обозначается операция вычисления расхода по формуле (13)):
