«И ТЕПЛООБМЕНА В ЗАДАЧАХ С КОНВЕКТИВНОЙ НЕУСТОЙЧИВОСТЬЮ И НЕЕДИНСТВЕННЫМ РЕШЕНИЕМ ...»

Министерство Образования Украины

ДНЕПРОПЕТРОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

На правах рукописи

Кудинов Павел Иванович

УДК 532.529

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ГИДРОДИНАМИКИ И ТЕПЛООБМЕНА

В ЗАДАЧАХ С КОНВЕКТИВНОЙ НЕУСТОЙЧИВОСТЬЮ

И НЕЕДИНСТВЕННЫМ РЕШЕНИЕМ

01.02.05 – механика жидкости, газа и плазмы Диссертация на соискание научной степени кандидата физико-математических наук

Научный руководитель Приходько Александр Анатольевич доктор физ.-мат. наук, проф.

Днепропетровск –

ОГЛАВЛЕНИЕ

ОГЛАВЛЕНИЕ

ВВЕДЕНИЕ

РАЗДЕЛ 1 ОБЗОР ЛИТЕРАТУРЫ И ВЫБОР НАПРАВЛЕНИЙ

ИССЛЕДОВАНИЯ

1.1. СОВРЕМЕННОЕ СОСТОЯНИЕ ЧИСЛЕННЫХ МЕТОДОВ МОДЕЛИРОВАНИЯ

ПРОЦЕССОВ ГИДРОДИНАМИКИ И ТЕПЛОМАССООБМЕНА В НЕСЖИМАЕМОЙ

ЖИДКОСТИ

1.2. МЕТОДЫ ОЦЕНКИ РАБОТОСПОСОБНОСТИ И ЭФФЕКТИВНОСТИ РАЗНОСТНЫХ

СХЕМ И АЛГОРИТМОВ

1.3. СВОБОДНОКОНВЕКТИВНЫЙ ТЕПЛООБМЕН В ГОРИЗОНТАЛЬНЫХ

КОЛЬЦЕОБРАЗНЫХ КАНАЛАХ С ГОРЯЧИМ ВНУТРЕННИМ ЭЛЕМЕНТОМ.................. 1.4. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ НЕСТАЦИОНАРНЫХ ОТРЫВНЫХ ТЕЧЕНИЙ....

1.5. ОСОБЕННОСТИ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ЗАДАЧ С НЕЕДИНСТВЕННЫМ

РЕШЕНИЕМ

1.6. ВЫВОДЫ

РАЗДЕЛ 2 ПРИМЕНЕНИЕ МЕТОДА КОНТРОЛЬНЫХ ОБЪЕМОВ К

РАСЧЕТУ ПРОЦЕССОВ ГИДРОДИНАМИКИ И ТЕПЛООБМЕНА В

НЕОРТОГОНАЛЬНЫХ КРИВОЛИНЕЙНЫХ КООРДИНАТАХ...........

2.1 ОБОБЩЕННЫЙ ЗАКОН СОХРАНЕНИЯ ФИЗИЧЕСКОЙ ВЕЛИЧИНЫ В

КОНТРОЛЬНОМ ОБЪЕМЕ

2.2. ДИСКРЕТНЫЙ АНАЛОГ ОБОБЩЕННОГО ЗАКОНА СОХРАНЕНИЯ В

КРИВОЛИНЕЙНОЙ НЕОРТОГОНАЛЬНОЙ СИСТЕМЕ КООРДИНАТ

2.2.1. Метод контрольного объема в случае осевой симметрии и при наличии закрутки потока

2.2.2. Особенности дискретизации уравнений сохранения количества движения в ортогональных криволинейных системах координат...........

2.3 SIMPLE МЕТОДЫ В КРИВОЛИНЕЙНОЙ НЕОРТОГОНАЛЬНОЙ СИСТЕМЕ

КООРДИНАТ

2.3.1. Вывод конечно-разностных соотношений на разнесенной сетке... 2.3.2. Процедуры SIMPLE, SIMPLEC, PISO и SIMPLER в криволинейной неортогональной системе координат

2.3.3. Вывод конечно-разностных соотношений на частично совмещенной сетке

2.3.4. Аппроксимация граничных условий

2.4. ВЛИЯНИЕ ТОЧНОСТИ ВЫЧИСЛЕНИЯ КОНВЕКТИВНОГО ПОТОКА НА РЕШЕНИЕ

ТЕСТОВОЙ ЗАДАЧИ О СДВИГОВОМ ТЕЧЕНИЕ В КРИВОЛИНЕЙНОЙ СИСТЕМЕ

КООРДИНАТ

2.5. РАЗНОСТНЫЕ СХЕМЫ ДЛЯ АППРОКСИМАЦИИ КОНВЕКТИВНОГО И

ДИФФУЗИОННОГО ПОТОКА

2.5.1. Схема аппроксимации диффузионных потоков

2.5.2. Схема QUIC

2.5.3. Схемы TVD для дозвуковых течений

2.6. ОСНОВНЫЕ РЕЗУЛЬТАТЫ

РАЗДЕЛ 3 ВЕРИФИКАЦИЯ И АНАЛИЗ ПРИМЕНИМОСТИ

РАЗРАБОТАННЫХ МЕТОДОВ И АЛГОРИТМОВ

3.1. СРАВНЕНИЕ ЭФФЕКТИВНОСТИ РАЗНОСТНЫХ СХЕМ И SIMPLE- ПОДОБНЫХ

АЛГОРИТМОВ

3.1.1. Развитие циркуляционного течения в каверне

3.1.2. Применение метода контрольного объема для описания свободноконвективного движения жидкости в приближении Буссинеска

3.1.3. Развитие свободной конвекции в квадратной области при подогреве сбоку

3.1.4. Оценка влияния кривизны сеточных линий на точность расчета локальных и интегральных характеристик течения и теплообмена.......

3.2. ТЕПЛООБМЕН В РЕЖИМЕ СВОБОДНОЙ КОНВЕКЦИИ В ГОРИЗОНТАЛЬНЫХ

КОЛЬЦЕВЫХ КАНАЛАХ

3.2.1. Свободная температурная конвекция в полости между горизонтальными цилиндрами

3.2.2. Задача о свободной конвекции в горизонтальной цилиндрической полости с внутренним, концентрично расположенным стержнем....... 3.3. ВЫВОДЫ.

РАЗДЕЛ 4. ОСОБЕННОСТИ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

НЕСТАЦИОНАРНЫХ ОТРЫВНЫХ ТЕЧЕНИЙ

4.1. МЕТОДИКИ РАСЧЕТА СТАЦИОНАРНЫХ И НЕСТАЦИОНАРНЫХ ЗАДАЧ.......... 4.1.1. Аппроксимация по времени и счет на установление

4.1.2. Развитие нестационарного отрыва потока при обтекании цилиндра.

4.1.3. Периодическое течение за цилиндром.

4.2. ВЫВОДЫ.

РАЗДЕЛ 5. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ НЕУСТОЙЧИВЫХ

ТЕЧЕНИЙ С НЕЕДИНСТВЕННЫМ РЕШЕНИЕМ

5.1. ТЕЧЕНИЕ ЗА СИММЕТРИЧНЫМ ВНЕЗАПНЫМ РАСШИРЕНИЕМ

5.2. СВОБОДНАЯ КОНВЕКЦИЯ В V-ОБРАЗНЫХ КАНАЛАХ.

5.3. ВЫВОДЫ.

ВЫВОДЫ

СПИСОК ИСПОЛЬЗОВАННОЙ ЛИТЕРАТУРЫ

ВВЕДЕНИЕ

Всестороннее исследование дозвуковых течений несжимаемой жидкости и газа является актуальной задачей гидромеханики, поскольку такие течения часто встречаются в научных и технических приложениях. Например, в различных теплообменниках [24], [8], [56], градирнях [94], [64], [63], камерах сгорания дисперсного топлива [118], [96], [119], [40], [89], вихревых устройствах [114], сепараторах и распылителях [128], [97], [119], [129], [130], [114], [6], [41].

энергетических и других ресурсов в масштабах всей Украины.

При проектировании технологического оборудования часто возникают задачи, в которых течение теряет устойчивость, или существует возможность реализации нескольких режимов течения при одном и том же наборе критериальных параметров. Среди известных задач такого типа можно назвать развитие различного рода возмущений течения в неустойчивом пограничном слое [126], [177], [69], [71], свободноконвективные течения с неустойчивой стратификацией [24], [163], [126], [72], [186], [182], [126], [70], [153], [36], [68], [37], [15], течения в плоских симметричных диффузорах и внезапных расширениях цилиндрических и сферических течениях [20], [172], [69], [2], [28], [16], [43], [4], [3].

исследований и фондов развития производства, экспериментальные работы по созданию новых технологий практически неосуществимы. В связи с этим, задача развития и совершенствования численных методов исследования течений вязкой несжимаемой жидкости становится особенно актуальной.

Наличие неединственности и неустойчивости решений задачи о течении вязкой несжимаемой жидкости требует повышения точности и эффективности произвольной формы. Особенно важно обеспечить свойство консервативности и монотонности численного метода в произвольной системе координат, чтобы свести к минимуму возмущения, вносимые методом в решение за счет особенностей разностной аппроксимации и кривизны сеточных линий. Для этой цели наиболее подходит семейство SIMPLE алгоритмов на основе метода контрольного объема [76], [161], [12], [59], [47] и схемы типа TVD, модифицированные для расчета несжимаемых течений [187].

Для оценки эффективности неявных методов, записанных в интегральной векторной форме для произвольной неортогональной системы координат, общепринятой методикой является решение тестовых задач. Существует необходимость в едином критерии эффективности, который позволял бы проводить сравнительную оценку, одновременно учитывающую точность, устойчивость и скорость сходимости метода.

неустойчивостью и неединственным решением, был выбран ряд задач, имеющих практическое и теоретическое значение.

Свободная конвекция в кольцеобразных каналах [56], [121] вызывает интерес, поскольку в этой задаче одновременно существуют зоны с устойчивой и неустойчивой стратификацией, что может привести к достаточно резким изменениям структуры течения и теплообмена [15], [53].

Нестационарное обтекание цилиндра является ярким примером развития неустойчивого отрывного течения [177]. В настоящей работе численно исследуются причины и механизмы развития вихревой неустойчивости в данной задаче.

Экспериментальные [171], [20], [150], [177], [181], [165] и расчетные данные [180], [142], [138], [115], [11], [149], [10] позволяют провести верификацию численного метода и подтвердить достоверность полученных результатов.

Экспериментальные данные [147], [148], [175], [146], [141] свидетельствуют неединственное решение. Благодаря обширному экспериментальному материалу, а также результатам расчетов [147], [175] данная задача может быть использована как тестовая для проверки применимости численного метода решения задач с неустойчивостью и неединственным решением.

Свободная конвекция при неустойчивой стратификации в горизонтальных каналах с V-образным поперечным сечением возникает при анализе течений в полости канала квадратного сечения с внутренним стержнем и продольными ребрами. Данная задача изучена мало, поэтому большой интерес представляет всестороннее исследование и поиск всех ее решений.

Для анализа ветвящихся решений задачи о течении вязкой несжимаемой жидкости в области произвольной формы представляет интерес вопрос о применимости различных экстремальных принципов типа минимума диссипации Гельмгольца [61], [174] и Пригожина [80] или теоремы Лагранжа об устойчивости положения равновесия механической системы [120].

Связь работы с научными программами, планами Диссертационная работа выполнена в соответствии с планом научноисследовательских работ Днепропетровского государственного университета.

Цели и задачи исследования Цель работы - совершенствование численных методов решения задач гидродинамики и теплообмена вязкой несжимаемой жидкости, разработка методик численного моделирования и исследование новых задач с неустойчивостью и неединственным решением, установление зависимостей интегральных и локальных характеристик от критериальных параметров.

Для достижения поставленной цели необходимо решить следующие задачи.

1. Выбрать универсальную форму записи исходных уравнений сохранения, независящую от используемой системы координат.

2. Обобщить SIMPLE- подобные алгоритмы на случай криволинейных неортогональных сеток.

3. Разработать методики уменьшения влияния кривизны сеточных линий на численное решение.

4. Разработать разностные алгоритмы с TVD свойствами для несжимаемых течений.

5. Создать и тестировать комплекс программ для расчета гидродинамики и тепломассообмена несжимаемых течений на основе уравнений Навъе-Стокса.

6. Разработать методики верификации и сравнения эффективности разностных алгоритмов и программ. Провести сравнение и добиться повышения точности и эффективности разностных схем и SIMPLE- подобных алгоритмов.

7. Определить закономерности развития свободноконвективных течений в горизонтальных каналах при устойчивой и неустойчивой стратификации.

8. Выяснить причины и механизмы потери устойчивости нестационарных отрывных течений. Определить пределы применимости для данных задач уравнений пограничного слоя при больших числах Рейнольдса.

9. Разработать методики получения решений задач с неединственным решением.

Исследовать зависимости интегральных характеристик течений от критериальных чисел в задачах с неединственным решением.

Научная новизна полученных результатов Научная новизна полученных результатов состоит в следующем:

1) предложено обобщение класса методов типа SIMPLE на разнесенных сетках для случая криволинейной неортогональной системы координат, на основе метода контрольного объема и TVD схем высокого порядка точности для расчета течений вязкой несжимаемой жидкости;

2) предложена поправка для вычисления конвективного потока через грань контрольного объема, благодаря которой линейное сдвиговое поле скорости тождественно удовлетворяет дискретным аналогам исходных уравнений;

3) численно исследована структура свободноконвективного течения в полости цилиндра с внутренним квадратным стержнем при числах Рэлея больше 105;

4) для задачи о нестационарном обтекании цилиндра при больших числах Рейнольдса исследованы пределы применимости уравнений пограничного слоя и механизмы развития неустойчивого течения в отрывной зоне;

5) исследовано влияние различных возмущений на процесс потери симметрии в задаче о течении за внезапным симметричным расширением, показано, что применение немонотонных схем может приводить к потере симметрии численного решения;

6) численно исследована задача о свободной конвекции с неустойчивой стратификацией в V-образных каналах, получено симметричное и два зеркальных несимметричных решения данной задачи, показано, что симметричное решение неустойчиво при числах Грасгофа больших некоторого критического значения, при числах Грасгофа меньше критического было получено единственное симметричное решение, устойчивое к возмущениям с масштабами сравнимыми с масштабами основного течения;

7) исследован вопрос о применимости экстремальных принципов типа теоремы Лагранжа об устойчивости положения механической системы и принципов минимума диссипации Гельмгольца и Пригожина для анализа решений нелинейных задач гидродинамики с неединственным решением.

Практическое значение полученных результатов Результаты настоящей работы могут найти широкое применение в области численного моделирования важных прикладных задач гидродинамики и теплообмена вязкой несжимаемой жидкости в областях произвольной формы, при наличии неустойчивости и неединственности решения.

Апробация результатов диссертации Материалы диссертационной работы докладывались, обсуждались и получили положительную оценку на:

гидродинамике” (г. Донецк, 1994 г.);

- международных конференциях ученых России, Белоруссии и Украины “Прикладные проблемы механики жидкости и газа” (г. Севастополь, 1995 – 1998 г.);

- международной конференции “Сучаснi проблеми водопостачання i знешкодження стiчних вод” (г. Львов, 1996 г.);

- международном симпозиуме “Advances in Computational Heat Transfer” (г.

Чешме, Турция, 1997 г.);

- 2-й республиканской научно-технической конференции “Гидроаэромеханика в инженерной практике” (г. Киев-Черкассы, 1997 г.);

- 7-м международном симпозиуме “Методы дискретных особенностей в задачах математической физики” (г. Феодоссия, 1997 г.);

- международной конференции “Наука i освiта - 98” (г. Днепропетровск, 1998 г.), - итоговых конференциях Днепропетровского государственного университета (г. Днепропетровск, 1995 - 1999 г.);

- объединенном семинаре Института технической механики НАН Украины (г. Днепропетровск, 1998 - 1999 г.);

- объединенном семинаре кафедр аэрогидромеханики и прикладной газовой динамики Днепропетровского государственного университета (г. Днепропетровск, 1999 г.);

- семинаре кафедры технической механики Днепропетровского государственного университета (г. Днепропетровск, 1998 - 1999 г.);

- всеукраинской молодежной научно-практической конференции “Людина і космос” (г. Днепропетровск, 1999 г.);

- республиканском семинаре Института гидромеханики НАН Украины под руководством академика НАН Украины В.Т. Гринченко (г. Киев, 1999 г.).

Публикации. Основное содержание диссертации опубликовано в работах.

Структура и объем диссертации. Диссертация состоит из введения, пяти глав, заключения и списка использованной литературы. Она изложена на страницах, из них 129 страниц текста, 131 рисунок и 6 таблиц, библиография содержит 188 наименований.

Обзор литературы и выбор направлений Бурное развитие численных методов и постоянное увеличение мощности ЭВМ предоставляют все более широкие возможности для численного исследования задач механики жидкости и газа. В настоящем разделе дан краткий обзор численных методов расчета несжимаемых течений, методик оценки эффективности и точности численных алгоритмов, а также работ, посвященных актуальным вопросам гидромеханики, связанным с конвективной неустойчивостью и неединственностью решения.

1.1. Современное состояние численных методов моделирования процессов гидродинамики и тепломассообмена в несжимаемой Течения жидкости и газа с дозвуковыми скоростями широко используются в науке и технике. На базе таких течений построены многочисленные теплообменники и технологические процессы [24], [8], [7], [79], [56], [123], [48], [90], [160], [55], [77], [87], [102], градирни и брызгальные бассейны [94], [64], [63].

Внутренние дозвуковые течения газовых смесей используются в химических лазерах [57], [58], при сжигании твердого и жидкого дисперсного топлива в энергетической промышленности [124], [118], [96], [119], [40], [89], [113] и т.д.

Вихревые сепараторы, распылители и инерционные классификаторы полидисперсных сред также работают при дозвуковых скоростях несущей жидкости [128], [97], [119], [129], [130], [114], [107], [92], [66], [5], [112], [6], [41].

Актуальной задачей является расчет силы сопротивления в судостроении [59], а также аэродинамических характеристик автомобилей [117], [159].

коэффициентов аэродинамического сопротивления летательных аппаратов и автомобилей, повышения эффективности сепарации и сжигания дисперсных сред определяющее значение имеют локальные характеристики течения, которые чаще всего невозможно получить интегральными методами расчета. Поэтому актуальной задачей является создание надежного численного метода расчета течений вязкой несжимаемой жидкости, который позволял бы учитывать реальную геометрию промышленных установок [49], [47], [46], [48], [51], [86], [83], [53], [82], [81], [84], [85], [169].

Наиболее общей математической моделью для описания течений сплошной вязкой жидкости являются уравнения Навье-Стокса [61], [17], [79]. Несмотря на то, что реальные процессы чаще всего протекают с большими числами Рейнольдса, использование приближения пограничного слоя [126], [61], [117], [180], [142], [138], [184], [115], хотя и позволяет существенно снизить вычислительные затраты, но не всегда физически оправдано [62], [69]. Четкие границы применимости уравнений пограничного слоя для внутренних течений установить достаточно сложно. При решении конкретных задач часто проводится сравнение результатов с экспериментом или расчетом по уравнениям Навье-Стокса [62].

Модель несжимаемой жидкости обычно используется для описания течений с числом Маха много меньше единицы, поскольку прямое применение методов, основанных на модели сжимаемого газа, становится не эффективным, из-за жесткости исходной системы уравнений [116], [98], [57], [74].

Исторически так сложилось, что разработка методов расчета течений вязкой несжимаемой жидкости шла по двум направлениям. Методы, использующие преобразованные переменные завихренность-функция тока (-) [30], [90] конкурировали с методами, основанными на записи исходных уравнений относительно естественных переменных скорость-давление (u, v, p) [116], [98], [57], [133], [76], [161], [12], [185], [157], [158], [154], [156], [23], [10]. Методы решения уравнений Навье-Стокса, записанные в преобразованных переменных (-), позволяют уменьшить количество неизвестных в двумерных задачах, однако для трехмерных течений количество неизвестных значительно увеличивается [133] по сравнению с методами использующими переменные скорость - давление. В методах (-) существуют определенные трудности с постановкой физически обоснованных граничных условий для завихренности на твердых стенках, а также для функции тока в многосвязных областях.

Методики расчета дозвуковых течений в переменных скорость - давление, основанные на модели несжимаемой жидкости и идее расщепления по физическим процессам, в свою очередь делятся на SIMPLE подобные алгоритмы [76], [161], [185], [12], [157], [158], метод искусственной сжимаемости [116], [98], [57], [23], [10], и смешанные Лагранжево - Ейлеровы методы типа метода частиц в ячейках [154], [156], [159]. Подробный обзор работ и анализ данных методов можно найти, например, в [74].

Одной из главных проблем, как для методов (-), так и методов (u, v, p), остается задача расчета поля давления. Источником трудностей является система уравнений несжимаемой жидкости, в которой отсутствует уравнение для давления.

При переходе к переменным завихренность-функция тока давление исключается из расчетов, поэтому определение силового воздействия на профиль в рамках таких методов представляется проблематичным. В методах, использующих естественные переменные, для расчета давления обычно используется комбинация уравнения неразрывности и количества движения. В таких алгоритмах как SIMPLE [76], SIMPLEC [161], [185], PISO [12] и в методе искусственной сжимаемости [116], [98], [57], фактически ставится эволюционная задача для давления, при этом отбрасываются слагаемые, которые обеспечивают прямую связь между массовыми силами и полем давления. Из анализа системы уравнений несжимаемой жидкости следует, что поле давления мгновенно реагирует на изменение полей вектора скорости и массовой силы и не зависит явно от распределения давления в предыдущий момент времени. В алгоритме SIMPLER реализована именно такая схема взаимодействия полей. Это приводит к заметному повышению эффективности алгоритма SIMPLER по сравнению с другими методами [76], [161], [185], [12], [157], [158], [49], [50], особенно для задач, где определяющее значение имеют массовые силы, в частности, для свободной температурной и концентрационной конвекции, а также течений с закруткой потока. Недостатком алгоритма SIMPLER можно считать то, что он реализуется только на разнесенной шахматной сетке, в то время как методы искусственной сжимаемости [116], [98], [57], [23], [10] и алгоритмы SIMPLE, SIMPLEC, PISO могут быть реализованы на совмещенных и частично совмещенных сетках [76], [161], [185], [12], [157], [158], [164], [140], [91], [144], [59], [31], [48]. На сегодняшний день наиболее популярным среди специалистов по вычислительной гидродинамике является семейство SIMPLE алгоритмов.

Реальная геометрия промышленных установок и устройств, требует, вопервых, возможности применения адаптивных сеток, а во-вторых, независимости свойств расчетного метода от выбранной системы координат [86], [83], [53], [82], [46], [47], [48], [49], [51]. С этой точки зрения наибольшей гибкостью обладает метод конечных элементов [167], [44], [159], [102], [21], [111]. Не смотря на последние достижения в области повышения точности аппроксимации в методе конечных элементов [44], [159], конечно-разностные алгоритмы, основанные на методе контрольного объема [76], превосходят метод конечных элементов по эффективности и точности [50]. В конечно-разностных алгоритмах проще реализуются возможности повышения точности результатов расчетов за счет применения TVD схем высокого порядка [42], [187], [155], [188], [178], [179].

Первая проблема, с которой столкнулись попытки применения SIMPLE процедур в неортогональных системах координат, было рассогласование полей скорости и давления с образованием так называемых шахматных полей давления и скорости [76], [12], [164]. Основная идея SIMPLE- подобных алгоритмов заключается в том, что для расчета давления используется разностное уравнение, полученное из дискретных аналогов уравнений количества движения и неразрывности. При этом большое значение имеет способ аппроксимации градиента давления в уравнении количества движения. В исходных алгоритмах [76], [161], [185], [12], [157], [158], записанных для декартовой системы координат на шахматной сетке, градиент давления вычисляется с помощью значений давления в двух соседних узлах. В этом случае не возникает рассогласование полей скорости и давления [164], а сеточный шаблон для давления включает 2D+1 узел, где D размерность задачи. В результате на регулярной сетке получается матрица системы алгебраических уравнений с 2D+1 диагональю, а для решения таких систем разработаны эффективные итерационные методы [65], [110], [27], [106].

Если уравнение количества движения записывается в проекции на базисные векторы декартовой системы координат, то при переходе к неортогональной системе координат на сетке близкой к ортогональной система дискретных уравнений для скорости и давления рассыпается на две независимые системы.

Большое число работ посвящено устранению этого эффекта. В [156] и [140] использована частично совмещенная сетка, где составляющие скорости хранятся в углах контрольного объема для давления, при этом авторами отмечались трудности с получением сходящегося решения из-за осцилляций поля давления. В [91] предложен метод коррекции давления на полностью совмещенной сетке, использующий предположение о близости сетки к ортогональной. В работе [144] предложен метод конечных объемов, использующий в качестве зависимых переменных контравариантные составляющие скорости, однако уравнения движения записаны не в строго консервативной форме. В работе [31] дана формулировка метода SIMPLER для ортогональных криволинейных координат. В работах М.П. Лобачева [59] уравнения движения записаны в консервативной форме на частично совмещенной сетке относительно контравариантных составляющих вектора скорости. Однако, в разностном уравнении для давления часть членов перенесена в правую часть, что привело к необходимости введения внутреннего итерационного цикла для решения уравнения Пуассона на каждой глобальной итерации. В любом случае использование совмещенных сеток делает невозможным применение алгоритма SIMPLER, который, как будет показано ниже, по эффективности превосходит остальные SIMPLE методы.

Вторая проблема, которая возникает при переходе к криволинейным координатам это влияние кривизны сеточных линий на точность численного решения, и в частности выполнение преобразования Галилея (тест на невозмущение однородного потока) для дискретных уравнений движения. Можно показать, что при записи дискретных аналогов для дифференциальной формы исходных уравнений и использовании в качестве искомых переменных контра- или ковариантных компонент вектора скорости [131] однородный поток, в общем случае, не удовлетворяет дискретным аналогам. Кроме того, можно показать, что выполнение на дискретном уровне преобразования Галилея еще не гарантирует, что такое простейшее течение, как линейное сдвиговое течение Куэтта может быть точно рассчитано на произвольной неортогональной сетке.

В работе [162] Леонард предложил схему второго порядка QUIC, которая на протяжении многих лет была очень популярна среди исследователей вязких несжимаемых течений [162], [8], [7], [59], [147], [9].

Работы В.П. Колгана [42] положили начало развития схем, которые впоследствии получили название TVD (Total Variation Diminishing) схем. На сегодняшний день монотонные схемы очень популярны среди специалистов по вычислительной аэродинамике сжимаемого газа [155], [188], [178]. Однако практически нет работ, в которых были бы отчетливо показаны преимущества схем типа TVD по сравнению с другими схемами для задач несжимаемой жидкости. В этом ключе также возникает задача разработки методологии сравнения различных разностных алгоритмов.

1.2. Методы оценки работоспособности и эффективности разностных На сегодняшний день создано огромное количество разностных алгоритмов. Каждый алгоритм обладает своими уникальными свойствами, которые играют решающую роль при численном исследовании прикладных задач. При создании нового численного метода или модификации существующего возникает вопрос об определении свойств метода.

В численном моделировании течений вязкой несжимаемой жидкости можно условно выделить два класса алгоритмов – итерационные методы решения нелинейных разностных уравнений, основанные на идеях расщепления по физическим процессам [76], [74] и разностные схемы аппроксимации исходных интегро-дифференциальных уравнений [93], [65], [74].

В качестве основной характеристики итерационных методов чаще всего используется скорость сходимости итерационного процесса решения нелинейных уравнений [161], [65], [166], [18], [48]. Для разностных схем наиболее важной оценкой является качество (точность) получаемого численного решения, при этом свойства схемы могут существенно влиять и на скорость сходимости (количество итераций) итерационного решения.

Качество численного решения определяется, в основном, не порядком аппроксимации разностной схемы [93], [65], а свойствами дифференциального приближения [127]. К сожалению, дифференциальное приближение трудно получить для нелинейных уравнений в криволинейной неортогональной системе координат [127]. Поэтому общепринятой практикой выяснения свойств разностных схем стало исследование на тестовых задачах, решение которых хорошо изучено. В частности, для внутренних двумерных течений вязкой несжимаемой жидкости такими тестовыми задачами стали: рециркуляционное течение в прямоугольной каверне с подвижной крышкой [7], [152], [18], [47], [137], [168], [135], [136], [145], [139], [173], [151], [122] и свободная конвекция в прямоугольной области при подогреве сбоку [99], [44], [79], [26]. Для задач внешнего обтекания – стационарное и нестационарное отрывное течение за круговым цилиндром [8], [7], [180], [142], [138], [184], [115], [20], [11], [149], [150], [159], [10], [151], [122], [177], [88], [33], [176], [181], [165], [171].

Обычным приемом при исследовании свойств разностных схем на тестовых задачах является сравнение по локальным и интегральным характеристикам решений, полученных на последовательно измельченных сетках [18].

Для сравнения итерационных алгоритмов и оценки влияния разностных схем на скорость сходимости строят зависимость от критериальных чисел количества итераций, необходимого для получения решения [161], [49], [47].

Отношение производительности к затратам является наиболее адекватной характеристикой эффективности любого вычислительного алгоритма. Количество итераций и оценка точности численного решения по отдельности не могут охарактеризовать эффективность алгоритма в целом, в терминах затраты – производительность. Поэтому нужны обобщенные критерии эффективности методов, основанные на характеристиках численного решения.

Важной с практической точки зрения задачей является определение степени влияния кривизны сеточных линий на расчет локальных и интегральных характеристик решения [101], [46], [47].

1.3. Свободноконвективный теплообмен в горизонтальных кольцеобразных каналах с горячим внутренним элементом Свободная конвекция в горизонтальных кольцеобразных каналах с горячим внутренним стержнем имеет большое практическое значение. Подобные задачи представляют особый интерес при проектировании ядерных реакторов, электронного оборудования, теплоизоляции передающих кабелей, систем аккумулирования тепла с неоднородным нагревом трубопроводов по угловой координате, а также приемников солнечных коллекторов тепла.

теоретических исследований. Подробный обзор литературы, анализ и детальные экспериментальные и численные исследования естественноконвективного переноса в зазоре между горизонтальными цилиндрами можно найти в работах Кьюэна и Голдстейна [160], [56], [55]. С помощью интерферометра Маха-Цандера ими были установлены распределения температур и найдены местные коэффициенты теплопередачи для воздуха и воды. Картины конвективных течений для воздуха, а также зависимости возникающих режимов от числа Грасгофа и соотношения диаметров полости рассматривались в работе [77]. Численное исследование данной задачи, как в ламинарном, так и в турбулентном режиме можно найти в работах [78], [123], [108], [87]. Имеется значительное число теоретических и экспериментальных работ по проблемам переноса в эксцентрических цилиндрических кольцеобразных областях. Так, были проведены измерения полных коэффициентов теплопередачи в воздухе для различных значений вертикальной и горизонтальной эксцентричностей внутреннего цилиндра [35], экспериментальные исследования влияния эксцентриситета и числа Рэлея на коэффициенты средней и местной теплоотдачи [56]. В работе [132] показано, что аналогичная задача возникает при анализе процесса плавления среды вокруг нагретого горизонтального цилиндра. В [95] проведено исследование теплоотдачи в кольцевом канале при наличии инверсии плотности.

Задача о свободной конвекции между наружным круговым цилиндром и внутренним стержнем играет важную роль для понимания процесса развития свободноконвективного пограничного слоя за острыми кромками, а также влияния этих кромок на характеристики теплоотдачи. В работе [121] численно и экспериментально анализировалась структура течения около квадратного стержня, расположенного концентрично внутри горизонтального кругового цилиндра при числах Рэлея меньше 105. Было получено, что в исследованном диапазоне чисел Рэлея отсутствует отрыв пограничного слоя при обтекании острых кромок.

Основной особенностью свободноконвективных течений в кольцеобразных каналах является одновременное присутствие зон устойчивой и неустойчивой температурной стратификации. Как известно [99], [51], [24], [163], [26], [34], [72], [186], [182], [183], в задачах с неустойчивой стратификацией можно ожидать резких изменений структуры течения, или даже ветвления решений при относительно малых приращениях числа Грасгофа. Другой особенностью течений с неустойчивой стратификацией является трехмерность течения, которая проявляется все интенсивнее с ростом числа Грасгофа и может существенно влиять на гидродинамику и теплообмен.

Наличие обширных экспериментальных данных по свободноконвективным течениям в горизонтальных кольцеобразных каналах [56], [123], [160], [55], [77], [78], [108], [35], [132], [95], [121], [87], [20] позволяет проверить работоспособность численного метода для различных критериальных чисел и геометрических размеров, и обоснованно делать численные прогнозы в неисследованных диапазонах изменения параметров [51].

1.4. Численное моделирование нестационарных отрывных течений нестационарные и периодические течения [122], [34], [126], [60] в которых малые возмущения могут привести к конечным изменениям структуры течения. К нестационарным течениям с неустойчивостью можно отнести турбулентные течения [104], [103], [105], [67], [20], нестационарный срыв потока [20], [122], слои смешения [104], [103], [105], [20], [122], [126], свободноконвективные течения с неустойчивой стратификацией [24], [34], [100], [25], [72], [186], [182], [183], [70], [153] и многие другие.

Описание механизма нестационарного срыва потока и зарождения турбулентного пограничного слоя является актуальной задачей для многих практически важных приложений [122], [34], [126], [60], [70], [153], [69], [172], [125]. В частности, резкие изменения аэродинамических характеристик профилей, крыльев, летательных аппаратов и автомобилей при малых изменениях угла атаки [122], [10], [159], [117], [88], [69], [125], а также динамические нагрузки на различные конструкции под действием постоянного и резко меняющегося набегающего потока [109], являются следствием нестационарного обтекания и срыва потока.

В качестве модельной задачи для данного класса течений часто используется развитие нестационарного пограничного слоя на поверхности Экспериментальные и численные исследования различных аспектов этих задач можно найти во многих работах [8], [7], [180], [142], [138], [184], [115], [20], [11], [149], [150], [159], [10], [151], [122], [177], [88], [33], [176], [181], [165], [171].

Накопленный экспериментальный материал и в частности подробные исследования нестационарных течений с помощью визуализации потока [20], а также измерения нестационарных характеристик потока в периодических течениях [171], [181], [165], позволяют подтвердить работоспособность численного метода.

Для решения задачи о нестационарном обтекании цилиндра, мгновенно приведенного в движение из состояния покоя, при больших числах Рейнольдса часто используется приближение пограничного слоя [180], [142], [138], [184], [115], [33]. В связи с этим представляет интерес задача определения пределов применимости уравнений пограничного слоя при больших числах Рейнольдса.

Большой интерес представляют причины, характер и детали зарождения и развития вихревой неустойчивости при нестационарном обтекании затупленных тел. Если начальная стадия зарождения неустойчивости достаточно хорошо описывается в рамках теории устойчивости пограничного слоя и метода малых колебаний [126], [61], [34], [25], [47], [70], [153], [69], [172], то дальнейшее развитие неустойчивости и переход к турбулентности представляет собой достаточно сложную задачу [104], [103], [105]. Эта задача может быть рассмотрена с помощью численных методов, использующих уравнения Навье-Стокса [47] в рамках ламинарной модели течения, до тех пор, пока вихри подсеточного масштаба [74] не начнут оказывать заметное влияние на решение.

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

1.5. Особенности численного моделирования задач с неединственным решением В прикладных задачах механики жидкости и газа самым распространенным видом движения сплошной среды является турбулентное течение, которое одновременно служит ярким примером неустойчивого взаимодействия вихревых структур различного масштаба [67], [126], [61], [104], [103], [105], [70], [153].

Практически любое ламинарное течение при увеличении числа Рейнольдса становится неустойчивым и малые возмущения, всегда присутствующие в реальных течениях, могут вызывать конечные изменения структуры потока и даже переход к турбулентному режиму течения [126], [122], [24], [34], [25], [99], [69], [172], [38]. К таким течениям можно отнести нестационарное обтекание цилиндра [171], [11], [177], [176], [181], [165], течение за симметричным внезапным расширением [147], [148], [170], [175], [146], [62], [141], [45] и в симметричном диффузоре [126], [22], [20], свободно-конвективные течения при неустойчивой стратификации [20], [24], [163], [186], [182], [183], [25], [47], [34], [100].

Отличительной особенностью перечисленных задач является возможность реализации нескольких решений при одном и том же наборе критериальных чисел.

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

Значительные результаты в этом направлении были достигнуты с помощью теории устойчивости ламинарного пограничного слоя и метода малых колебаний [126], [61], [34], [25], [47], [70], [153], [69], [172]. Применение теории устойчивости и метода малых колебаний сводится к следующему [126], [172], [25], [47]: задается аналитическое представление течения, которое исследуется на устойчивость, на основное течение накладывается возмущение в виде колебания с определенной длинной волны, которое удовлетворяет уравнению типа Орра-Зоммерфельда [126], [172], [25], [47], ставится задача на отыскание собственных значений. В зависимости от величины найденных собственных значений делается вывод о затухании или возрастании данного возмущения с течением времени. Основным ограничением данного подхода является необходимость в аналитическом представлении для решения, исследуемого на устойчивость [25]. Кроме того, существуют задачи, которые имеют несколько устойчивых к малым возмущениям решений при одном и том же наборе критериальных параметров (например, течение за внезапным симметричным расширением [147], [148], [170], [175], [146], [62], [141], [45] и свободная конвекция в квадратной области при подогреве снизу [47], [186]). Из анализа устойчивости одного из устойчивых решений таких задач не следует, что существуют или не существуют другие решения.

В механике твердого деформируемого тела методы решения задач с неединственным решением получили широкое развитие благодаря применению теории ветвлений решений нелинейных уравнений [19]. В частности, теория устойчивости оболочек при неосесимметричной деформации [1] позволяет строить решение во всем диапазоне нагрузок, определять и классифицировать особые точки решения, устанавливать их связь с критическими нагрузками и несущей способностью рассматриваемых конструкций.

Для нахождения точки ветвления решения используются идеи, заложенные в задаче о неявных функциях [19], [1]. Пусть F (x, y ) – нелинейный оператор, где x – искомое решение, y – параметр, и пусть (x0, y 0 ) является решением уравнения F ( x0, y 0 ) = 0. Возникает вопрос о решениях уравнения где x, y – приращение (вариация) решения.

Запишем разложение Тейлора для оператора F в окрестности (x0, y 0 ) :

где A = A(x0, y 0 ) = – матрица Якоби, N – высшие члены разложения по x и Если матрица Якоби невырожденная, то решение продолжается по параметру y единственным образом. В особой точке неединственно и, согласно теореме о неявных функциях [19], матрица Якоби A – вырождена:

Если F (x, y ) сеточный оператор, а (1.2) представляет собой линеаризацию задачи (1.1) [1], явную относительно приращения сеточного решения, то матрица линейных алгебраических уравнений для приращения сеточного решения является матрицей Якоби A. Таким образом, в особой точке (x0, y 0 ) определитель матрицы A обращается в ноль [1].

Когда найдена особая точка, можно построить решения выходящие из этой точки. В работе [1] приведены примеры построения решений, выходящих из точки ветвления для задач устойчивости оболочек.

Вернемся теперь к задачам с неединственным решением в механике жидкости. Наиболее эффективные численные методы вычислительной гидродинамики являются многошаговыми неявными или полунеявными [76], [161], [185], [12], [157], [158]. Поэтому возникает проблема построения матрицы Якоби для этих методов. Далее встает алгебраическая проблема вычисления определителей и собственных значений матриц большой размерности [110], [27], [106], которая тесно связана с общими ограничениями возможностей современных ЭВМ. В частности, обычная сетка для двумерной задачи насчитывает порядка узлов. В зависимости от конкретной задачи, на каждом итерационном шаге в каждом узле должно хранится от 5 до 10-ти переменных. Учитывая, что матрица Якоби в случае неявного метода может быть полностью заполненной, то она будет сдержать от 2.5·109 до 1010 вещественных чисел. Даже если использовать одинарную машинную точность (что совершенно недопустимо при вычислении определителя матрицы такого размера [110], [27], [106]), то для хранения в памяти такой матрицы понадобится от 1010 до 4·1010 байт, или от 10 до 40 гигабайт. Это примерно в тысячу раз больше, чем доступный объем оперативной памяти и в 2– раз больше чем объем всего дискового пространства современного персонального компьютера.

Учитывая все вышесказанное, представляется целесообразным разбить решение проблемы нахождения всех решений задачи о течении вязкой несжимаемой жидкости в произвольной области на два направления. В рамках первого направления следует развивать и тестировать высокоэффективные численные методы, которые вносили бы минимальные возмущения в численное решение при расчете течений вязкой несжимаемой жидкости в областях произвольной формы. В рамках второго направления следует разрабатывать методики применения теории ветвления решений нелинейных уравнений [19] для поиска особых точек и построения ветвящихся решений на базе, созданных в рамках первого направления, численных методов.

Важной с практической точки зрения задачей является сравнительный анализ устойчивости различных решений. А именно анализ возможности перехода с одной ветви решения на другую вследствие внешнего воздействия определенной формы и величины. В связи с этим возникает вопрос о применимости различного рода экстремальных принципов типа теоремы Лагранжа [120] об устойчивости положения равновесия механической системы, или ее аналога для механики твердого деформируемого тела [1], где устойчивость положения равновесия связана со знаком второй вариации потенциальной энергии деформации. Также возникает вопрос о применимости широко известных принципов минимума диссипации Гельмгольца [61], [174] и Пригожина [120].

Проведенный анализ современного состояния численных методов расчета несжимаемых течений с конвективной неустойчивостью и неединственностью решения позволяет сделать следующие выводы.

1. Течения жидкости и газа с дозвуковыми скоростями часто встречаются в науке и технике, что делает актуальной задачу создании новых и совершенствование существующих численных методов расчета, позволяющих учитывать реальную геометрию промышленных установок.

2. Расчетные методы должны основываться на наиболее общей математической модели сплошной вязкой жидкости - уравнениях Навье-Стокса.

3. Наиболее перспективным направлением на пути создания универсальной вычислительной методики является обобщение семейства SIMPLE подобных алгоритмов на случай произвольной криволинейной системы координат, а также повышение эффективности и точности алгоритмов, за счет применения TVD схем высокого порядка аппроксимации.

4. Актуальной задачей является исследование влияния кривизны сеточных линий на точность численного решения и разработка методик уменьшения этого влияния.

5. Необходимы методики адекватной и всесторонней оценки эффективности разностных алгоритмов, учитывающие одновременно точность, скорость сходимости по сетке и затраты машинного времени и памяти в виде единого универсального параметра эффективности.

6. Свободная конвекция в горизонтальных кольцеобразных каналах с горячим внутренним элементом имеет большую практическую важность. Основной особенностью таких течений является одновременное присутствие зон устойчивой и неустойчивой температурной стратификации, что может привести экспериментальных данных позволяет проверить работоспособность неисследованных диапазонах изменения параметров.

турбулентного пограничного слоя является актуальной задачей для многих практически важных приложений. В качестве модельной задачи для данного класса течений целесообразно использовать развитие нестационарного пограничного слоя на поверхности цилиндра. Большой интерес представляют собой причины, характер и детали развития неустойчивости в нестационарном течении, а также пределы применимости уравнений пограничного слоя для нестационарного обтекания цилиндра при больших числах Рейнольдса.

8. Прикладные задачи механики жидкости и газа часто связаны с течениями, в которых малые возмущения вызывают конечные изменения структуры потока.

Отличительной особенностью таких задач является возможность реализации нескольких решений при одном и том же наборе критериальных чисел.

Учитывая огромную практическую ценность такого класса течений, актуальной задачей является создание вычислительного аппарата, позволяющего находить все возможные решения задачи и отмечать наиболее вероятные (устойчивые) варианты. Наиболее значительные результаты в этом направлении были достигнуты с помощью теории устойчивости ламинарного пограничного слоя и метода малых колебаний. Основным ограничением данного подхода является необходимость в аналитическом представлении для решения, исследуемого на устойчивость.

9. На пути применения теории ветвления решений нелинейных уравнений существуют многочисленные проблемы, связанные в частности с недостаточной мощностью современных ЭВМ. Поэтому решение проблемы нахождения всех решений задачи о течении вязкой несжимаемой жидкости в произвольной области предлагается разбить на два направления. Первое направление является одной из целей настоящей работы, а именно, создание высокоэффективных численных методов для расчета течении вязкой несжимаемой жидкости в областях произвольной формы. В рамках второго направления следует разработать методики применения теории ветвления решений нелинейных уравнений для поиска особых точек и построения решений на базе созданных численных методов.

10. Для анализа задач с неединственным решением интересным является вопрос о применимости и ценности различного рода экстремальных принципов типа теоремы Лагранжа и принципов минимума диссипации Гельмгольца и Пригожина.

Применение метода контрольных объемов к расчету процессов гидродинамики и теплообмена в неортогональных криволинейных координатах Большое число технологических процессов основано на использовании тепломассообмена в дозвуковом течении жидкости. Сложность реальной геометрии промышленных установок требует от метода расчета независимости от выбора системы координат.

В ряде работ [7], [164], [154], [156], [140], [91], [144], [59], [31] обсуждаются методы расчета дозвуковых течений областях сложной формы, и в частности с подвижными границами [154], [156]. В работах [7], [164], [140], [91], [144], [59], [31] обсуждается несколько вариантов обобщения методов типа SIMPLE на случай криволинейной системы координат. Основными недостатками предложенных подходов являются ограничения связанные с требованием ортогональности сетки и невозможность реализовать метод SIMPLER [76] на частично или полностью совмещенных сетках. Метод, используемый в работе [59], позволяет вести расчеты на неортогональной сетке, однако применение частично совмещенных сеток позволяет использовать только метод SIMPLE или SIMPLEC. В работе [31] дана формулировка метода SIMPLER для случая ортогональных криволинейных координат.

В настоящей работе предложено обобщение метода SIMPLER на случай неортогональной криволинейной системы координат. Дискретизация исходной системы уравнений сохранения проведена на основе метода контрольного объема таким образом, что дискретные аналоги удовлетворяют преобразованию Галилея.

Исследовано влияние кривизны сеточных линий на точность численного решения.

Показано, что линейное сдвиговое течение может быть точно рассчитано на произвольной неортогональной сетке, только при повышении порядка точности квадратурной формулы для вычисления конвективного потока через грани контрольного объема. Описаны различные схемы разностной аппроксимации конвективных и диффузионных потоков, в частности схема Леонарда QUIC и схемы TVD второго и третьего порядка.

2.1 Обобщенный закон сохранения физической величины в Для получения универсального алгоритма численного исследования процессов гидродинамики и тепломассообмена в криволинейных неортогональных координатах воспользуемся законами сохранения физических величин в контрольном объеме (КО).

Рассмотрим течение вязкой несжимаемой жидкости в области произвольной формы. В качестве исходных уравнений выберем законы сохранения, записанные в интегральной векторной форме.

где W контрольный объем, S поверхность W, n внешняя нормаль к S, V вектор скорости, плотность (считается постоянной в рамках модели несжимаемой жидкости), Р тензор напряжений, G вектор массовых сил.

Для реализации численных методов часто используется обобщенное уравнение переноса в дифференциальной форме [76], [7], [164], [140], [91], [144], [31]. При использовании криволинейных систем координат такой подход менее удобен [47], [59], [23], чем применение интегральной векторной формы записи исходных уравнений в виде обобщенного закона сохранения физической величины в контрольном объеме W где коэффициент диффузии, S мощность источников. Если в (2.3) положить =1, =0, S=0, то получим уравнение неразрывности (2.1), если внутренняя энергия, = коэффициент теплопроводности, S источники тепла, то (2.3) будет представлять собой закон сохранения внутренней энергии несжимаемой жидкости. Покажем, что уравнение сохранения количества движения (2.2) также может быть записано в форме (2.3). Согласно обобщенной гипотезе Ньютона [61] где µ коэффициент вязкости, I единичный тензор, p давление, div V = 0 в несжимаемой жидкости, контравариантные компоненты тензора скоростей деформаций равны где g ik метрический тензор, Vi ковариантные компоненты вектора скорости.

Используя формулу Остроградского-Гаусса, выражения (2.4) и (2.6), преобразуем в (2.2) слагаемое с тензором напряжений Учитывая, что ковариантная производная от контравариантного метрического дифференцирования [14], [39], получим Если вязкость µ постоянна и divV=0, то Подставим (2.7)-(2.9) в (2.2) получим Очевидно, если в (2.10) положить =V и =µ, то (2.10) по форме совпадает с (2.3).

Если коэффициент вязкости и плотность не постоянны, то в правую часть (2.10) следует добавить слагаемое Уравнение (2.3) это универсальная форма записи обобщенного закона сохранения физической величины (скалярной или векторной) в контрольном объеме. На основе уравнения (2.3) удобно строить унифицированные численные алгоритмы, которые существенно облегчают разработку многоцелевых пакетов программ.

2.2. Дискретный аналог обобщенного закона сохранения в криволинейной неортогональной системе координат Будем требовать, чтобы дискретный аналог (ДА) уравнения (2.3) обладал свойством консервативности независимо от выбора системы координат. Тогда решение =const должно тождественно удовлетворять дискретному аналогу (2.3) в произвольной системе координат. Нетрудно видеть, что для =V=const это условие эквивалентно условию выполнения принципа Галилея. Кроме того, будем стремиться избегать появления в дискретном аналоге уравнения количества движения символов Кристоффеля, ибо погрешность в аппроксимации производных от компонент метрического тензора может привести к нарушению принципа Галилея для разностного решения.

Пусть в расчетной области введена криволинейная неортогональная система координат xi, i=1,,D, где D размерность задачи. Разобьем расчетную область на множество непересекающихся контрольных объемов. Введем локальную нумерацию узлов и граней контрольного объема следующим образом: узлу в центре присвоим номер 0, узлу, смещенному относительно нулевого на шаг вперед вдоль i-ой координаты, присвоим номер k=2i, а узлу, смещенному на шаг назад номер k=2i-1, грань, лежащая между нулевым и k-ым узлом, будет иметь номер k.

Таким образом, если в трехмерной области введена регулярная сетка, то узлу (l, m, n) соответствует локальный узел 0, аналогично для соседних узлов - (l-1, m, n)=1, (l+1, m, n)=2, (l, m-1, n)=3, (l, m+1, n)=4, (l, m, n-1)=5, (l, m+1, n+1)=6. В двумерной области узлу (n, m) соответствует локальный узел 0, для соседних узлов - (n-1, m)=1, (n+1, m)=2, (n, m-1)=3, (n, m+1)=4 см. рис.2.1.

Если область не односвязная, то можно применять блочные сетки, О сетки, или ставить условия периодичности.

Рассмотрим особенности построения дискретного аналога обобщенного закона сохранения (2.3). Запишем аппроксимацию слагаемого, отвечающего за конвективный перенос где k=1,,2D, Sk площадь грани с номером k, величины в скобках должны вычисляться в центре k-ой грани.

Следует заметить, что от способа определения значения ( ) k на грани по значениям k в соседних узах во многом зависит точность и устойчивость метода.

Для начала, примем линейную зависимость ( ) k от k Рис.2.1. Нумерация узлов и граней контрольного объема Рис.2.2. Контрольные объемы, схема расположения и нумерации узлов и граней.

а) - контрольный объем для давления P; б, в) - контрольные объемы для контравариантных компонент вектора скорости- V 1, V 2.

Черными точками обозначены узлы для давления P, крестиками - для V 1, кружками для V здесь ck коэффициент схемы, который может зависеть от локального сеточного числа Пекле Например, для центрально-разностной схемы для противопоточной схемы имеет место зависимость Центрально-разностная схема (CR) становится неустойчивой при сеточном числе Пекле больше 2. Противопоточная схема (UDS) вносит в решение чрезмерную схемную диффузию, но устойчива при любых числах Пекле и обеспечивает диагональное преобладание матрицы системы алгебраических уравнений. Поэтому неявную часть (2.12) запишем по противопоточной схеме (2.15), а повышения порядка аппроксимации [65], [93] можно добиться, корректируя источниковые члены [8], [7]. Например, центрально-разностная схема может быть записана следующим образом:

где UDS, n значения на текущей и предыдущей итерации, рассчитанные по соответственно.

Представим вектор скорости в виде где e j контравариантные векторы локального базиса криволинейной системы координат, V j контравариантные компоненты вектора скорости. Вектор (n)k внешней единичной нормали к поверхности Sk равен где k=2i-1, 2i; ei i-ый ковариантный вектор локального базиса, суммирования по i нет. Проекция скорости на нормаль к поверхности Sk будет равна Окончательно выражение для конвективного слагаемого будет иметь вид Рассмотрим слагаемое, отвечающее за диффузионный перенос в (2.3) здесь k=2i-1, 2i; в последнем слагаемом суммирование по j, ji, d k схемный коэффициент, зависящий от сеточного числа Пекле (2.14); g ij = ei e j компоненты метрического тензора; j вдоль k-ой грани контрольного объема.

Используя выражения (2.20) и (2.21), запишем дискретный аналог (2.3) где n1 - значение переменной на предыдущем временном слое, неразрывности. В ходе итерационного решения нелинейной системы уравнений неразрывности, т.е. второе слагаемое в (2.24) не равно нулю Для того чтобы матрица системы алгебраических уравнений (2.22) обладала свойством диагонального преобладания - необходимым условием сходимости итерационных методов решения систем линейных алгебраических уравнений, второе слагаемое в (2.24) перенесем в правую часть (2.22) [76]. Тогда (2.24) и (2.25) примут вид Пусть источник в обобщенном уравнении сохранения (2.3) S=0.

Подставим n 1 = = const в (2.22). Если поле скорости удовлетворяет уравнению неразрывности, то (2.22) обращается в тождество, независимо от выбора системы координат. Следовательно (2.22) обладает свойством консервативности. Отсюда также следует, что решение = V = const тождественно удовлетворяет дискретным аналогам уравнений количества движения и неразрывности на произвольной криволинейной неортогональной сетке, т.е. выполняется принцип Галилея для дискретного решения системы (2.22). Чтобы получить систему алгебраических уравнений относительно компонент вектора скорости в некотором базисе, следует спроектировать (2.22) на вектора этого базиса (см. раздел 2.3). При этом в уравнение (2.22) не входят символы Кристоффеля [131] (производные от компонент метрического тензора), которые трудно аппроксимировать с достаточной точностью на произвольной криволинейной сетке.

2.2.1. Метод контрольного объема в случае осевой симметрии и при наличии закрутки потока Пусть имеет место осевая симметрия течения и угловая скорость вращения жидкости относительно оси симметрии не равна нулю. Рассмотрим уравнение сохранения количества движения (2.10). Из формулы Остроградского - Гаусса следует, что [61]:

где a - может быть скаляром или тензором.

Рассмотрим второе слагаемое в правой части (2.10).

Первые два слагаемые в (2.31) по форме такие же, как и в декартовой прямоугольной системе координат, за исключением формулы для элементарной площадки:

Следовательно, и дискретный аналог этих слагаемых будет такой же как и для плоского случая, за исключением выражений для вычисления объемов и площадей. В случае осевой симметрии течения площадь грани, и объем вычисляются по формулам где l - ширина грани.

последнее слагаемое в (2.31) более подробно исключением производных от базисных векторов) и выражения (2.34) и (2.35) можно упростить.

Таким образом, уравнение сохранения количества движения в случае осесимметричных течений с закруткой потока можно представить в виде где Введем в расчетной области цилиндрическую ортогональную систему координат (x, r, ) с базисом (ex, er, e) и цилиндрическую неортогональную систему координат (x1, x2, ) с базисом (e1, e2, e). Используем выражения для компонент скорости в этих системах координат здесь суммирование по i, а eix = e i e x, eir = e i e r - проекции контравариантных векторов базиса криволинейной системы координат на базисные векторы цилиндрической системы координат.

Учитывая выражения (2.39) перепишем (2.38) в виде:

Рассмотрим грань контрольного объема (в системе (x1, x2, ) грани КО лежат на поверхностях xi =const, угол между касательной к грани и направлением радиальной оси обозначим ) где величины с индексом k определяют положение центра тяжести грани, l ширина грани. Площадь грани равна а объем вычисляется по формуле (2.33).

Таким образом, в случае осевой симметрии течения с закруткой потока для вывода дискретного аналога можно применять уравнения сохранения количества движения в форме (2.37) с источниковым членом в форме (2.38) для ортогональной и (2.40) для неортогональной системы координат.

При выводе дискретного аналога обобщенного уравнения сохранения (2.22) в цилиндрической системе координат формулы (2.12)-(2.21) остаются в силе за исключением выражения для объема (2.33) и площади грани КО (2.42).

Если вязкость и плотность жидкости – переменные величины, то выражение (2.40) принимает вид 2.2.2. Особенности дискретизации уравнений сохранения количества движения в ортогональных криволинейных системах координат Рассмотрим особенности дискретизации уравнений сохранения количества движения в ортогональных криволинейных системах координат. Покажем, что в общем случае, следуя традиционному подходу к построению дискретных аналогов на основе дифференциальной формы записи исходных уравнений, невозможно получить точное решение в виде невозмущенного течения V=const.

Пусть в расчетной области введена ортогональная криволинейная система уравнение количества движения примет вид или в тензорной форме коэффициенты Ляме, g ij = e i e j – компоненты метрического тензора. Даже если V=const, компоненты вектора скорости в криволинейной системе координат V i, в общем случае являются нелинейными функциями координат xi. Заменяя частные производные в уравнении (2.45) центральными разностями мы получим дискретный аналог уравнения (2.45), решение которого, в общем случае, не будет совпадать с решением уравнения (2.45), поскольку центральные разности точно аппроксимируют производную для полинома второй степени. Т.е. V=const, в общем случае, не будет удовлетворять дискретному аналогу уравнения сохранения количества движения.

Для примера рассмотрим случай, когда в расчетной области удобно ввести полярную систему координат (r, ) (например, при обтекании кругового цилиндра).

Запишем уравнения сохранения количества движения (2.10) в дифференциальной форме [61] для случая идеальной несжимаемой жидкости Пусть имеет место стационарное, одномерное течение жидкости V x = 1, V y = 0.

Тогда компоненты вектора скорости в полярной системе координат будут нелинейными функциями угла а уравнения (2.46) и (2.47) примут вид Введем в расчетной области сетку (Ri = R0 + (i–1) R; k = 0+ (k–1)·) и заменим в (2.50), (2.51) производные конечными разностями С помощью прямой подстановки нетрудно убедиться, что (2.48), (2.49) являются решением дифференциальных уравнений (2.50), (2.51). Проверим, являются ли выражения (2.48), (2.49) решением дискретных аналогов (2.52), (2.53).

Подставим (2.48), (2.49) в уравнения (2.52), (2.53), получим Из (2.54), (2.55) следует, что одномерный невозмущенный поток не является решением дискретных аналогов (2.52), (2.53). Как следствие нарушается свойство консервативности метода. Таким образом, применение стандартного подхода к выводу дискретного аналога может внести существенную погрешность в решение, особенно на грубых сетках.

Ранее было показано, что в случае применения, предлагаемого в настоящей работе подхода (пп.2.2), для получения дискретных аналогов исходных уравнений, решение в виде невозмущенного потока тождественно удовлетворяет разностным уравнениям (2.22), а метод обладает свойством консервативности на любой сетке. В пп.2.4 будет показано, что при минимальных поправках предлагаемая методика дает точное численное решение задачи о линейном сдвиговом течении на произвольной сетке.

2.3 SIMPLE методы в криволинейной неортогональной системе Основная идея SIMPLE- подобных алгоритмов заключается в том, что для расчета давления используется разностное уравнение, полученное из дискретных аналогов уравнений количества движения и неразрывности. Большое значение имеет способ аппроксимации градиента давления в уравнении количества движения. В исходных алгоритмах [76], [161], [185], [12], [157], [158], записанных для декартовой системы координат, на шахматной сетке градиент давления вычисляется с помощью значений давления в двух соседних узлах. В этом случае не возникает рассогласование полей скорости и давления [164], а сеточный шаблон для давления включает 2D+1 узел, где D размерность задачи. В результате на регулярной сетке получается матрица системы алгебраических уравнений с 2D+ диагональю, а для решения таких систем разработаны эффективные итерационные методы [65].

Если уравнение количества движения записывается в проекции на базисные векторы декартовой системы координат, то при переходе к неортогональной системе координат для аппроксимации градиента давления необходимо использовать более двух узловых точек. Увеличение количества узлов, участвующих в аппроксимации градиента давления, ведет к расширению сеточного шаблона уравнения для давления. В результате разностное уравнение Пуассона будет иметь матрицу с 3D диагоналями. Прямое решение такой системы алгебраических уравнений может оказаться достаточно обременительным, поэтому часть членов обычно переносят в правую часть, чтобы получить матрицу с 2D+ диагональю. Однако в этом случае требуется введение внутреннего итерационного цикла для решения уравнения Пуассона, что также снижает эффективность алгоритма. Кроме того, без введения специальных поправок или процедур, поле давления оказывается слабо связанным с полем скорости, возникает эффект шахматного поля давления. Причина возникновения этого эффекта в том, что задача расчета сеточного поля давления по заданной величине градиента давления имеет бесконечное множество решений, если для аппроксимации градиента используется более двух сеточных узлов.

относительно контравариантных компонент вектора скорости на шахматной разнесенной сетке. В этом случае градиент давления определяется по значениям давления в двух соседних узлах, что позволяет избежать рассогласования полей скорости и давления и без изменений перенести исходные алгоритмы [76], [161], [185], [12], [157], [158], разработанные и тестированные на декартовой сетке, на случай областей сложной геометрической формы. Кроме того, появляется возможность применять метод SIMPLER [76] на криволинейных неортогональных сетках, который показал свою высокую надежность, устойчивость и эффективность [76], [161], [49], [47], [166].

2.3.1. Вывод конечно-разностных соотношений на разнесенной сетке Рассмотрим метод, позволяющий получить разностное уравнение для давления на шаблоне с 2D+1 узлом, в произвольной криволинейной системе координат. В выражении для дискретного аналога обобщенного закона сохранения (2.22) последовательно положим =1, =0, и =V, =µ, получим дискретный аналог исходной системы уравнений (2.1), (2.2) где В случае осевой симметрии в (2.58) входит слагаемое (2.40), площади граней КО в (2.56), (2.57) вычисляются по формуле (2.42) а объем КО вычисляется по формуле (2.33).

Слагаемое B p представляет собой дискретный аналог слагаемого с градиентом давления Запишем проекции (2.59) на e i контравариантные векторы локального базиса криволинейной системы координат. Пусть область разбита сеткой с шахматной схемой расположения узлов [76], [154], показанной на рис.2.2.

Выберем в (2.59) в качестве W контрольные объемы для V i (рис.2.2.б, рис.2.2.в) и умножим (2.59) на (e i ) где P1, P2 значение давления в узлах на границах контрольных объемов рис.2.2.б, xi расстояние между узлами P1, P2. Итак, если использовать рис.2.2.в, разнесенную шахматную сетку и проектировать уравнения количества движения на контравариантные векторы локального базиса криволинейной системы координат, то градиент давления определяется по значениям давления в двух точках - P1, P2.

Учитывая выражение (2.60) запишем проекции (2.57) на (e i ) слагаемые с произведением e j e i, ji перенесем в источниковый член и введем обозначения здесь ji. Перепишем (2.61) с учетом обозначений (2.62) в виде Нетрудно видеть, что (2.56), (2.63) по форме совпадают с уравнениями, которые были использованы в [76], [161], [185], [12], [157], [158], для вывода SIMPLE- подобных алгоритмов в декартовой системе координат, что позволяет оставить эти алгоритмы без изменений для произвольной системы координат.

Уравнения (2.63) обладают также и другими важными свойствами. Во-первых, V=const очевидно является решением (2.63), во вторых, в уравнениях не участвуют символы Кристоффеля, несмотря на то, что искомыми переменными являются контравариантные компоненты вектора скорости [131].

Введем обозначение тогда (2.63) примет вид Подставим (2.65) в уравнение неразрывности (2.56) где где a0 k коэффициент a0 в (2.63), определенный на k- ой грани контрольного объема для давления.

Итак, задав начальное приближение для скорости V и поле массовых сил, мы можем непосредственно рассчитать давление P по (2.66). Возьмем в качестве P и рассчитаем промежуточное поле скорости V, которое P в (2.63) удовлетворяет разностному уравнению количества движения (2.63) и, в общем случае, может не удовлетворять уравнению неразрывности (2.56).

Пусть V, P удовлетворяют уравнениям (2.56), (2.63). Обозначим поправку давления Последовательно подставим V, P и V, P в уравнение (2.63), а затем вычтем из первого получившегося уравнения второе Пренебрегая в (2.69) подчеркнутым членом получим Уравнение для поправки давления получим из уравнения (2.70) и уравнения неразрывности (2.56), где A0, Ak определяются по формулам (2.67).

2.3.2. Процедуры SIMPLE, SIMPLEC, PISO и SIMPLER в криволинейной неортогональной системе координат Предложенный в предыдущем пункте подход позволяет оставить без изменения SIMPLE- подобные процедуры (детально описанные в [76], [185], [161], [157], [158], [12]) при переходе к криволинейным неортогональным координатам.

Поэтому мы лишь кратко опишем данные процедуры, и остановимся более детально на их вычислительных свойствах.

Процедура SIMPLE состоит из следующих этапов [76]:

1. Задается начальное приближение для поля давления Р*.

2. Рассчитывается V по формуле (2.63).

3. Вычисляется поправка давления P по формуле (2.71).

4. Корректируется поле скорости по формуле (2.70) и поле давления по формуле (2.68).

5. Используя скорректированное давление в качестве P*, возвращаемся ко 2-му этапу. Цикл продолжается до достижения сходимости.

Процедура SIMPLEC отличается от SIMPLE только способом вычисления коэффициентов уравнения для поправки давления (2.71) [185], [161], [7]. В частности в (2.67) рассчитанные на k- ой грани КО для давления.

В процедуре PISO [157], [158], [12] первая часть совпадает с пунктами 1- процедуры SIMPLE. После выполнения 4-го пункта процедуры SIMPLE поле скорости может не удовлетворять уравнению количества движения. Поэтому в [157], [158], [12] предложено проводить еще одну коррекцию, так, чтобы скорректированные поля скорости и давления удовлетворяли линеаризованному уравнению количества движения (2.63) и уравнению неразрывности (2.56).

Процедура SIMPLER [76]:

1. Задается начальное приближение для поля скорости V.

2. Рассчитывается P по формуле (2.66).

3. Рассчитывается V по формуле (2.63).

4. Вычисляется поправка давления P по формуле (2.71).

5. Корректируется поле скорости по формуле (2.70) (поле давления не корректируется).

6. Если сходимость не достигнута, возвращаемся ко 2-му этапу. Цикл продолжается до достижения сходимости.

Исследуя свойства изложенных процедур можно разделить их на две группы. К первой группе можно отнести процедуру SIMPLER, в которой по заданному полю скорости и массовым силам непосредственно вычисляется давление. Во второй группе методов, к которой относятся SIMPLE, SIMPLEC, PISO, давление рассчитывается в результате последовательности поправок относительно начального приближения, которые вычисляются из условия выполнения уравнения неразрывности.

Из уравнений движения несжимаемой жидкости следует, что поле давления не имеет предыстории и мгновенно реагирует на изменение полей вектора скорости и массовой силы. В алгоритме SIMPLER реализована именно такая модель взаимодействия полей. Во второй группе методов фактически ставится эволюционная задача для давления, попутно отбрасываются члены, которые обеспечивают прямую связь между массовыми силами и давлением. В результате методы второй группы менее эффективны для нестационарных задач, задач типа проточного течения в канале или трубе переменного сечения с мягкими граничными условиями на выходе из области, течений, где определяющее значение имеют массовые силы, или источниковые члены отвечающие за переносное ускорение в цилиндрической системе координат. В частности, для задач о свободной температурной и концентрационной конвекции, а также для течений с закруткой потока и в задачах с мягкими граничными условиями для скорости и давления процедура SIMPLER обеспечивает большую устойчивость и скорость сходимости итерационного решения, по сравнению с остальными методами.

Следует отметить, что для стационарных задач решения, полученные по каждому из рассмотренных методов, между собой ничем не отличаются, поскольку решаются одни и те же разностные уравнения. Различия проявляются только в скорости сходимости (количество итераций необходимых для достижения заданной погрешности решения по отношению к решению, полученному на предыдущей итерации) и устойчивости итерационного процесса (способность получать сходящееся решение при минимальных значениях параметров релаксации).

SIMPLEC и PISO, как улучшенные варианты процедуры SIMPLE, обладают более высокой скоростью сходимости. На простых задачах PISO приближается по скорости сходимости к SIMPLER. Выигрыш в суммарном времени счета по методу SIMPLER достигается благодаря существенному уменьшению количества итераций по сравнению остальными методами, в то время как одна итерация SIMPLER требует не намного больше времени, чем SIMPLEC или PISO.

Недостатком алгоритма SIMPLER можно считать то, что он реализуется только на разнесенной шахматной сетке, в то время как методы второй группы могут быть реализованы на совмещенных и частично совмещенных сетках [164], [140], [91], [144], [59], что позволяет сократить требования к объему оперативной памяти ЭВМ. Однако, учитывая, что темпы удешевления микросхем памяти значительно опережают рост производительности процессоров, данный недостаток не следует считать существенным.

2.3.3. Вывод конечно-разностных соотношений на частично совмещенной сетке При переходе к криволинейной системе координат объем вычислений существенно возрастает. Использование разнесенной сетки предполагает составление и решение системы алгебраических уравнений для каждой компоненты скорости отдельно, а также использование интерполяции для определения недостающих компонент скорости в узлах разнесенных сеток. Если все компоненты скорости будут храниться в одном узле, то разностные уравнения сохранения для компонент скорости будут отличаться только правыми частями, также отпадет необходимость в интерполяции вектора скорости с одной сетки на другую. К сожалению, в этом случае становится невозможным применение процедуры SIMPLER. Однако, для относительно простых задач, можно добиться некоторой экономии машиной памяти и времени счета, применяя SIMPLEC или PISO, записанные для частично совмещенной криволинейной сетки.

Приведем основные выражения, необходимые для реализации SIMPLEподобных методов на частично совмещенной сетке, при условии, что матрица системы алгебраических уравнений для давления будет иметь 2D+1 диагональ на регулярной сетке. Узлы для компонент скорости расположим в углах КО, а узлы для давления, температуры и других переменных в центре КО (рис.2.3).

Уравнения сохранения массы (2.56) и количества движения (2.57) не меняются. Выражение (2.60) для аппроксимации слагаемого с давлением запишем в виде где Pk = ( Pk 1 + Pk 2 ) 2 - среднее значение давления на k- ой грани КО для скорости (см. рис.2.3).

При выводе уравнения для поправки давления будем действовать аналогично подразделу 2.3.1. Вместо выражения (2.69), записанного в проекциях на базисные векторы запишем то же выражение, но в векторном виде (см. 2.3.2) где P конечно-разностная аппроксимация градиента поправки давления.

Рис.2.3. Контрольные объемы, схема расположения и нумерации узлов и граней.

а) - контрольный объем для давления P; б) - контрольный объем для вектора скорости- V.

Черными точками обозначены узлы для давления P, крестиками для вектора скорости V Рис.2.4. Схема нумерации узлов вблизи границы S расчетной области для двух типов Пренебрегая подчеркнутым слагаемым, интерполируем (2.73) в центр k- ой грани КО для давления (см. рис.2.3) Умножим (2.74) на ( 1) k e i e i, i=1,2; k=2i-1, 2i, и проинтегрируем по поверхности КО для давления здесь ij. Заметим, что на сетке близкой к ортогональной произведение e j e i 0.

Если итерационный процесс сходится, то V 0 и вторым слагаемым в левой части (2.75) можно пренебречь [91], учитывая, что V удовлетворяет уравнению неразрывности, получим уравнение для поправки давления на частично совмещенной сетке:

которое по форме полностью совпадает с (2.71) полученным для разнесенной сетки, и имеет 2D+1 диагональ на регулярной сетке. Коэффициенты в (2.76) вычисляются по формулам При реализации SIMPLE процедур на частично совмещенной сетке последовательность операций, изложенная в пункте 2.3.2, остается неизменной.

Уравнения также не изменяются, за исключением двух формул: вместо выражения (2.60) следует использовать (2.72), а вместо уравнения для поправки давления (2.71), используется уравнение (2.76), где коэффициенты рассчитываются по формулам (2.77).

2.3.4. Аппроксимация граничных условий Точность численного решения во многом зависит от метода постановки граничных условий. Применительно к задачам конвективного тепломассообмена можно выделить условия в невозмущенном потоке, условия на твердой поверхности, условия симметрии течения и условия периодичности, которые применяются для расчета в областях с регулярно повторяющейся по пространству геометрией, например теплообменники с пакетами труб или пластин [8], [7], [9].

Подробное обсуждение особенностей постановки граничных условий при решении задач внешнего обтекания и течений в каналах с помощью SIMPLE методов можно найти в [76], [8], [7], [9].

Рассмотрим постановку граничных условий на твердой стенке. Основная проблема заключается в уменьшении толщины пограничного слоя по мере роста числа Рейнольдса, что приводит к необходимости использовать более подробные разностные сетки в пристеночной области.

На рис.2.4 изображены два возможных варианта взаимного расположения границы области и узлов расчетной сетки. В первом случае (см. риc.2.4,а) для задания значения скорости или температуры на твердой поверхности S, необходимо ввести аппроксимацию для профиля соответствующей искомой переменной, вблизи границы. Наиболее распространенным является метод, основанный на линейной аппроксимации, т.е.

где i - значение функции в граничном узле, i 1 - значение функции в ближайшем к границе узле (см. рис.2.4). Соотношение (2.78) имеет первый порядок точности и приводит к известному соотношению Тома для завихренности на стенке [90]. Повысить порядок точности можно, используя квадратичный или кубический полином, аппроксимирующий профиль вблизи стенки [8], [7], [9]. В частности для квадратичной аппроксимации получаем где где xi, xi 1, xi 2 - координаты соответствующих узлов, а x S координата границы (см. рис.2.4,а). В работах [8], [7], [9] приводится обзор более сложных вариантов постановки граничных условий, основанных на анализе уравнений пограничного слоя.

составляется таким образом, что значение давления в граничном узле Pi не участвует в расчетной процедуре. Само значение Pi можно найти из условия равенства нулю нормального к стенке градиента давления [76], [8], [7], [9]. Отсюда следует, что методе в SIMPLER поправка давления в граничных узлах должна быть равна нулю, а в SIMPLE- подобных алгоритмах нужно задавать условие равенства нулю нормальной производной к границе от поправки давления.

Граничные условия в набегающем потоке ставятся в зависимости от типа задачи. Если расход известен заранее, то задается величина и направление вектора скорости на участке границы, где скалярное произведение скорости и внешней нормали к границе области меньше нуля. На той части границы, где жидкость вытекает из расчетной области, для вектора скорости и остальных неизвестных обычно ставятся так называемые "мягкие" граничные условия типа Однако, в этом случае на выходной границе области поле скорости может не удовлетворять дискретному уравнению неразрывности (2.56), что в свою очередь, вносит в возмущение в расчет поля давления и течения во всей области.

Поэтому для аппроксимации граничных условий следует использовать комбинацию уравнений (2.81) и (2.56). Те компоненты вектора скорости, узлы которых находятся непосредственно на границе области рис.2.4.б, должны рассчитываться из дискретного уравнения неразрывности (2.56), а компоненты скорости, находящихся в фиктивных узлах рис.2.4.а – из уравнения (2.81).

Если известен перепад давления и требуется определить расход, то на внешней границе ставятся мягкие условия для скорости. Для давления можно применять итерационную процедуру типа где где PS - известное среднее значение давления на некотором участке границы, суммирование в (2.83) ведется по всем узлам i-1 (см. рис.2.4), принадлежащим этому участку границы. Очевидно, в результате итерационного процесса (2.82) Подчеркнем, что граничные условия ставятся для значений давления в узлах i- (см. рис.2.4), поскольку вычислительная процедура построена таким образом, что поле давления не зависит от значения давления в граничных узлах i [76], [8], [7], [9].

2.4. Влияние точности вычисления конвективного потока на решение тестовой задачи о сдвиговом течение в криволинейной системе вычисления интегрального конвективного потока через грани КО на численное решение тестовой задачи о сдвиговом течение в криволинейной системе координат.

Пусть имеет место чисто сдвиговое стационарное течение при отсутствии градиента давления и массовых сил. Направление вектора скорости постоянно, а модуль является линейной функцией пространственных переменных. Потребуем, чтобы сдвиговое течение точно удовлетворяло дискретному аналогу уравнения количества движения (2.3) при = V, S=0 в криволинейной неортогональной системе координат. Нетрудно видеть, что для линейного сдвигового течения диффузионное слагаемое тождественно равно нулю, как в дифференциальной, так и в конечно-разностной форме. Тогда уравнение сохранения количества движения примет вид При переходе к дискретному аналогу (2.84) необходимо аппроксимировать поверхностные интегралы по граням КО. Обычно для этого используется квадратурная формула прямоугольников, точная для полиномов первой степени, если значение подынтегральной функции вычислено в центре отрезка интегрирования.

В данном случае применение формулы прямоугольников будет вносить погрешность, поскольку подынтегральная функция Vn( V ) представляет собой произведение линейно меняющихся функций (суммирование ведется по граням КО).

В общем случае, величина ( V ) k на k- ой грани вычисляется по значениям V в центрах соседних с гранью ячеек. Для противопоточной схемы ( V ) равно значению V в центре ячейки, находящейся вверх по потоку. В работе [101] в соответствии с принципом минимальных значений производной [42], для криволинейной системы координат, предложено брать в качестве минимальное по модулю из вычисленных значений V справа и слева в центре тяжести грани.

Найдем такую поправку к квадратурной формуле прямоугольников, чтобы формула была точной для произведения двух линейно меняющихся функций. Пусть a и b линейно зависят от переменной s, тогда где первое слагаемое соответствует квадратурной формуле прямоугольников, а второе слагаемое представляет собой погрешность.

С учетом (2.86) формула (2.85) в двумерном случае примет вид (VnV) k определяются в центре тяжести грани, ( Vn), ( V ) приращения где вдоль грани, которые рассчитываются по формуле (2.86).

На рис.2.5 приведено решение тестовой задачи о течении Куэтта между параллельными пластинами. В центре области сетка была преднамеренно искривлена (рис.2.5.а). На рис.2.5.б отчетливо видны возмущения продольной составляющей скорости Vx при использовании обычной формулы прямоугольников.

На рис.2.5.в видно, что изолинии Vx, полученные с использованием формулы (2.87), не зависят от кривизны сеточных линий.

Рис.2.5 Тест на не возмущение течения Куэтта на криволинейной сетке а) – криволинейная сетка; б) изолинии Vx (формула прямоугольников);

в) изолинии Vx (скорректированная формула прямоугольников) прямоугольников, независимо от порядка точности, с которым конечно-разностная схема аппроксимирует значение порядок аппроксимации конвективного потока будет не выше первого.

Погрешность порядка ( Vn) ( V ) S 12 может стать существенной в областях с большими градиентами скорости. И даже при линейном изменении поля вектора скорости (как при чисто сдвиговом течении) на криволинейной сетке, или при не совпадении направления течения с направлением сеточных линий в декартовой системе координат, в решение будут вноситься возмущения, зависящие от кривизны сеточных линий. Применение поправки (2.87) позволяет точно рассчитывать конвективные потоки в случае сдвигового течения на криволинейной сетке.

2.5. Разностные схемы для аппроксимации конвективного и 2.5.1. Схема аппроксимации диффузионных потоков Повышения точности численного решения можно добиться двумя путями: с помощью увеличения количества сеточных узлов, и повышением порядка точности разностной аппроксимации [93], [65] исходных уравнений.

Для аппроксимации конвективных и диффузионных потоков нужно определить зависимость схемных коэффициентов d k и c k в (2.23) от сеточного числа Пекле. Для этого можно воспользоваться одномерной задачей конвективнодиффузионного переноса [76] с решением Из (2.89), (2.90) следует, что для линейной одномерной задачи конвективнодиффузионного переноса типа (2.88) на равномерной сетке схемные коэффициенты в (2.23) можно вычислять по формулам где Pek - сеточное число Пекле (2.14). При этом численное решение будет совпадать с аналитическим решением (2.89).

Выражения (2.91), (2.92) определяют “экспоненциальную” схему [76]. При Pek 0 эта схема переходит в центрально-разностную схему для конвективного слагаемого ck = 0.5, а для диффузионного При Pek > 10 экспоненциальная схема переходит в противопоточную схему (2.15) для конвективного слагаемого, а вкладом диффузионного потока в перенос можно пренебречь поскольку Вычисление экспонент занимает много времени, и к тому же, для двумерной нелинейной задачи, решение (2.89), в общем случае, не является точным и вносит сильную численную диффузию в решение [8], [7]. Поэтому в [76] было предложено аппроксимировать решение (2.89), (2.90) полиномом 5-го порядка.

В настоящей работе предлагается вычислять приближенные значения (2.91), (2.92) с помощью полиномов первого и второго порядка Выражения (2.95), (2.96), (2.97) определяют схему, которая обладает всеми свойствами экспоненциальной схемы. Эффект схемной вязкости особенно сильно проявляется в случае отрывных, циркуляционных течений, когда не удается адаптировать сетку к структуре потока, так чтобы одно из семейств сеточных линий повторяло форму линий тока.

Результаты серии тестовых расчетов циркуляционных течений подтвердили, что для вычисления диффузионного потока на грани КО выражение (2.96) дает приемлемые результаты, а для конвективного слагаемого в двумерных задачах необходимо применять другие подходы, нежели (2.91) или (2.95). Обзор литературы по разностным схемам для несжимаемых течений можно найти в работах [8], [7].

В работе [162] Леонард предложил схему второго порядка QUIC. В схеме QUIC для аппроксимации значения одномерного случая имеет вид а в двумерном случае [8], [7], [147], [9] Рис.2.6. Нумерация узлов при вычислении ограничителя потоков на k-ой грани где, - координаты локальной системы координат, начало которой находится в узле 0 при (Vn )k 0, и в узле k при (Vn )k < 0 (см. рис.2.6). С помощью такого приема повышается устойчивость разностной схемы, поскольку коэффициенты интерполяционного многочлена (2.99) вычисляются с учетом того факта, что область влияния численного решения расположена вверх по потоку. Согласно схеме QUIC для двумерного случая (2.99) величина ( )k вычисляется по формуле [162], [8], [7], [147], [9] где нижние индексы соответствуют номерам узлов на рис.2.6.

антидиффузионным эффектом. В [188] показано, что при расчете по схеме QUIC нестационарного одномерного уравнения конвективного переноса скаляра с разрывами в начальных условиях возникают осцилляции решения в окрестности разрывов. Уравнения движения вязкой несжимаемой жидкости не допускают существования разрывных решений [61], поэтому антидиффузионные свойства схемы QUIC проявляются в виде завышения значений локальных максимумов, особенно на грубых сетках и в областях с большими градиентами искомых переменных. Причина этого явления заключается в завышенной или заниженной оценке конвективного потока в окрестности локального максимума или минимума, что приводит к неоправданному росту локального экстремума. Частично устранить данный эффект можно вводя ограничение на значения локального экстремума. В частности можно потребовать, чтобы промежутку значений справа и слева от k-ой грани: [ k, 0 ], если в узле 0 или k находится локальный экстремум функции Ф. Также можно вычислять ( )k в окрестности локального экстремума по противопоточной схеме (2.15) или схеме (2.95). Однако даже после введения такой поправки схема QUIC дает заметную погрешность при расчете значений Ф в окрестности локальных экстремумов. Эта погрешность становится особенно заметной в области вторичных течений с линейными размерами, сравнимыми с шагом сетки (см. рис.7.а). Кроме того, скорость сходимости численного метода при использовании схемы QUIC значительно (в зависимости от задачи в 2-4 раза) меньше, чем при использовании противопоточной схемы. На основании перечисленных недостатков схемы QUIC, можно сделать вывод, что эффективные численные алгоритмы следует строить, основываясь на монотонных схемах повышенного порядка точности.

2.5.3. Схемы TVD для дозвуковых течений После публикации работы В.П. Колгана [42] в 1972 начало развиваться новое поколение схем аппроксимации конвективного потока, которые впоследствии получили название TVD (Total Variation Diminishing) схем. Отличительным свойством этих схем является монотонность получаемого решения. Сочетание повышенного порядка точности, высокой скорости сходимости и отсутствие осцилляций в окрестности разрывов сделало схемы TVD очень популярными среди специалистов по вычислительной газовой динамике [155], [188], [178]. Однако до недавнего времени для расчета несжимаемых течений чаще всего использовалась схема QUIC [162], [8], [7], [59], [147], [9], которая не обладает свойством монотонности.

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

Можно показать [187], что практически любая схема TVD может быть записана в виде где (UDS )k - значение функции на грани КО рассчитанное по противопоточной схеме (2.15), k нелинейный ограничитель потоков (flux limiter [187], [188], [178]), который является функцией значений в окрестности k-ой грани контрольного объема (см. рис.2.6).

Запишем выражения для вычисления ( )k по схеме Колгана MinMod [42] эта задача имеет три решения: симметричное и два зеркальных не симметричных.

Симметричное решение неустойчиво к малым начальным возмущениям при Gr>Grk, в то время как не симметричные решения устойчивы к малым возмущениям, по крайней мере, при числах Грасгофа Gr Grk) необходимо, чтобы метод решения системы линейных алгебраических уравнений (СЛАУ) давал симметричное решение на каждой итерации. Поэтому при прогонках по вертикальным сеточным линиям применялась следующая последовательность перебора: от оси симметрии к краям области и обратно.

В результате расчетов с разными (симметричными и несимметричными) последовательностями прогонок было получено, что асимметрия прогонок является источником возмущений, которые могут стать причиной перехода от неустойчивого симметричного решения к несимметричному, но само стационарное асимметричное решение не зависит от реализации метода решения СЛАУ.